The TELEMAC-MASCARET system  trunk
om0101.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE om0101
3 ! *****************
4 !
5  &(op ,dm,typdim,xm,typexm,dn,typdin,xn,typexn,d,c,
6  & ikle,nelem,nelmax,ndiag)
7 !
8 !***********************************************************************
9 ! BIEF V6P1 21/08/2010
10 !***********************************************************************
11 !
12 !brief OPERATIONS ON MATRICES WITH P1 SEGMENT.
13 !code
14 !+ D: DIAGONAL MATRIX
15 !+ C: CONSTANT
16 !+
17 !+ OP IS A STRING OF 8 CHARACTERS, WHICH INDICATES THE OPERATION TO BE
18 !+ PERFORMED ON MATRICES M AND N, D AND C.
19 !+
20 !+ THE RESULT IS MATRIX M.
21 !+
22 !+ OP = 'M=N ' : COPIES N IN M
23 !+ OP = 'M=CN ' : MULTIPLIES N BY C
24 !+ OP = 'M=M+CN ' : ADDS CN TO M
25 !+ OP = 'M=M+N ' : ADDS N TO M
26 !+ OP = 'M=TN ' : COPIES TRANSPOSE OF N IN M
27 !+ OP = 'M=M+TN ' : ADDS TN TO M
28 !+ OP = 'M=MD ' : M X D
29 !+ OP = 'M=DM ' : D X M
30 !+ OP = 'M=DMD ' : D X M X D
31 !+ OP = 'M=X(M) ' : NOT SYMMETRICAL FORM OF M
32 !+ (OLD MATSNS)
33 !+ OP = 'M=MSK(M)' : MASKS M EXTRADIAGONAL TERMS
34 !+ (OLD MASKEX)
35 !+ THE MASK IS TAKEN FROM D
36 !
37 !code
38 !+ CONVENTION FOR THE STORAGE OF EXTRA-DIAGONAL TERMS:
39 !+
40 !+ XM(IELEM,1) ----> M(1,2)
41 !+ XM(IELEM,2) ----> M(2,1)
42 !
43 !history J-M HERVOUET (LNHE) ; F LEPEINTRE (LNH)
44 !+ 05/02/91
45 !+ V5P6
46 !+
47 !
48 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
49 !+ 13/07/2010
50 !+ V6P0
51 !+ Translation of French comments within the FORTRAN sources into
52 !+ English comments
53 !
54 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
55 !+ 21/08/2010
56 !+ V6P0
57 !+ Creation of DOXYGEN tags for automated documentation and
58 !+ cross-referencing of the FORTRAN sources
59 !
60 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
61 !| C |-->| A GIVEN CONSTANT USED IN OPERATION OP
62 !| D |-->| A DIAGONAL MATRIX
63 !| DM |<->| DIAGONAL OF M
64 !| DN |-->| DIAGONAL OF N
65 !| IKLE |-->| CONNECTIVITY TABLE.
66 !| NDIAG |-->| NUMBER OF TERMS IN THE DIAGONAL
67 !| NELEM |-->| NUMBER OF ELEMENTS
68 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
69 !| OP |-->| OPERATION TO BE DONE (SEE ABOVE)
70 !| TYPDIM |<->| TYPE OF DIAGONAL OF M:
71 !| | | TYPDIM = 'Q' : ANY VALUE
72 !| | | TYPDIM = 'I' : IDENTITY
73 !| | | TYPDIM = '0' : ZERO
74 !| TYPDIN |<->| TYPE OF DIAGONAL OF N:
75 !| | | TYPDIN = 'Q' : ANY VALUE
76 !| | | TYPDIN = 'I' : IDENTITY
77 !| | | TYPDIN = '0' : ZERO
78 !| TYPEXM |-->| TYPE OF OFF-DIAGONAL TERMS OF M:
79 !| | | TYPEXM = 'Q' : ANY VALUE
80 !| | | TYPEXM = 'S' : SYMMETRIC
81 !| | | TYPEXM = '0' : ZERO
82 !| TYPEXN |-->| TYPE OF OFF-DIAGONAL TERMS OF N:
83 !| | | TYPEXN = 'Q' : ANY VALUE
84 !| | | TYPEXN = 'S' : SYMMETRIC
85 !| | | TYPEXN = '0' : ZERO
86 !| XM |-->| OFF-DIAGONAL TERMS OF M
87 !| XN |-->| OFF-DIAGONAL TERMS OF N
88 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
89 !
90  USE bief, ex_om0101 => om0101
91 !
92  USE declarations_special
93  IMPLICIT NONE
94 !
95 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
96 !
97  INTEGER, INTENT(IN) :: NELEM,NELMAX,NDIAG
98  INTEGER, INTENT(IN) :: IKLE(nelmax,2)
99 !
100  DOUBLE PRECISION, INTENT(IN) :: DN(*),D(*)
101  DOUBLE PRECISION, INTENT(INOUT) :: DM(*)
102  DOUBLE PRECISION, INTENT(INOUT) :: XM(nelmax,*)
103  DOUBLE PRECISION, INTENT(IN) :: XN(nelmax,*)
104 !
105  DOUBLE PRECISION, INTENT(IN) :: C
106 !
107  CHARACTER(LEN=8), INTENT(IN) :: OP
108  CHARACTER(LEN=1), INTENT(INOUT) :: TYPDIM,TYPEXM,TYPDIN,TYPEXN
109 !
110 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
111 !
112  INTEGER I,J,IELEM
113 !
114 !-----------------------------------------------------------------------
115 !
116  IF(op(3:8).EQ.'N ') THEN
117 !
118  IF(typdin(1:1).EQ.'Q') THEN
119  CALL ov('X=Y ', x=dm, y=dn, dim1=ndiag)
120  ELSEIF(typdin(1:1).EQ.'I'.OR.typdin(1:1).EQ.'0') THEN
121 ! NOTHING TO DO, ONLY NEEDS TO COPY TYPDIN
122  ELSE
123  WRITE(lu,6) typdin(1:1)
124 6 FORMAT(1x,'OM0101 (BIEF) : TYPDIN UNKNOWN :',a1)
125  CALL plante(1)
126  stop
127  ENDIF
128  typdim(1:1)=typdin(1:1)
129 !
130  IF(typexn(1:1).EQ.'S') THEN
131  CALL ov('X=Y ', x=xm(1,1), y=xn(1,1), dim1=nelem)
132  ELSEIF(typexn(1:1).EQ.'Q') THEN
133  CALL ov('X=Y ', x=xm(1,1), y=xn(1,1), dim1=nelem)
134  CALL ov('X=Y ', x=xm(1,2), y=xn(1,2), dim1=nelem)
135  ELSEIF(typexn(1:1).NE.'0') THEN
136  WRITE(lu,11) typexn(1:1)
137 11 FORMAT(1x,'OM0101 (BIEF) : TYPEXN UNKNOWN :',a1)
138  CALL plante(1)
139  stop
140  ENDIF
141  typexm(1:1)=typexn(1:1)
142 !
143 !-----------------------------------------------------------------------
144 !
145  ELSEIF(op(3:8).EQ.'CN ') THEN
146 !
147  CALL ov('X=CY ', x=dm, y=dn, c=c, dim1=ndiag)
148 !
149  IF(typexn(1:1).EQ.'S') THEN
150  CALL ov('X=CY ', x=xm(1,1), y=xn(1,1), dim1=nelem)
151  ELSEIF(typexn(1:1).EQ.'Q') THEN
152  CALL ov('X=CY ', x=xm(1,1), y=xn(1,1), dim1=nelem)
153  CALL ov('X=CY ', x=xm(1,2), y=xn(1,2), dim1=nelem)
154  ELSEIF(typexn(1:1).NE.'0') THEN
155  WRITE(lu,11) typexn(1:1)
156  CALL plante(1)
157  stop
158  ENDIF
159 !
160  typdim(1:1)=typdin(1:1)
161  typexm(1:1)=typexn(1:1)
162 !
163 !-----------------------------------------------------------------------
164 !
165  ELSEIF(op(3:8).EQ.'M+CN ') THEN
166 !
167  IF(typdin(1:1).EQ.'I') THEN
168  CALL ov('X=X+C ', x=dm, c=c, dim1=ndiag)
169  ELSEIF(typdin(1:1).NE.'0') THEN
170  CALL ov('X=X+CY ', x=dm, y=dn, c=c, dim1=ndiag)
171  ENDIF
172 !
173  IF(typexn(1:1).EQ.'S') THEN
174  CALL ov('X=X+CY ', x=xm(1,1), y=xn(1,1), c=c, dim1=nelem)
175  IF(typexm(1:1).EQ.'Q') THEN
176  CALL ov('X=X+CY ', x=xm(1,2), y=xn(1,1), c=c, dim1=nelem)
177  ENDIF
178  ELSEIF(typexn(1:1).EQ.'Q') THEN
179  IF(typexm(1:1).NE.'Q') THEN
180  WRITE(lu,98) typexm(1:1),op(1:8),typexn(1:1)
181 98 FORMAT(1x,'OM0101 (BIEF) : TYPEXM = ',a1,' DOES NOT GO',
182  & /,1x,'FOR THE OPERATION : ',a8,' WITH TYPEXN = ',a1)
183  CALL plante(1)
184  stop
185  ENDIF
186  CALL ov('X=X+CY ', x=xm(1,1), y=xn(1,1), c=c, dim1=nelem)
187  CALL ov('X=X+CY ', x=xm(1,2), y=xn(1,2), c=c, dim1=nelem)
188  ELSEIF(typexn(1:1).NE.'0') THEN
189  WRITE(lu,11) typexn(1:1)
190  CALL plante(1)
191  stop
192  ENDIF
193 !
194 !-----------------------------------------------------------------------
195 !
196  ELSEIF(op(3:8).EQ.'M+N '.OR.
197  & (op(3:8).EQ.'M+TN ').AND.typexn(1:1).NE.'Q') THEN
198 !
199  CALL ov('X=X+Y ', x=dm, y=dn, dim1=ndiag)
200 !
201  IF(typexn(1:1).EQ.'S') THEN
202  CALL ov('X=X+Y ', x=xm(1,1), y=xn(1,1), dim1=nelem)
203  IF(typexm(1:1).EQ.'Q') THEN
204  CALL ov('X=X+Y ', x=xm(1,2), y=xn(1,1), dim1=nelem)
205  ENDIF
206  ELSEIF(typexn(1:1).EQ.'Q') THEN
207  IF(typexm(1:1).NE.'Q') THEN
208  WRITE(lu,98) typexm(1:1),op(1:8),typexn(1:1)
209  CALL plante(1)
210  stop
211  ENDIF
212  CALL ov('X=X+Y ', x=xm(1,1), y=xn(1,1), dim1=nelem)
213  CALL ov('X=X+Y ', x=xm(1,2), y=xn(1,2), dim1=nelem)
214  ELSEIF(typexn(1:1).NE.'0') THEN
215  WRITE(lu,11) typexn(1:1)
216  CALL plante(1)
217  stop
218  ENDIF
219 !
220 !-----------------------------------------------------------------------
221 !
222  ELSEIF(op(3:8).EQ.'TN ') THEN
223 !
224  CALL ov('X=Y ', x=dm, y=dn, dim1=ndiag)
225 !
226  IF(typexn(1:1).EQ.'S') THEN
227  CALL ov('X=Y ', x=xm(1,1), y=xn(1,1), dim1=nelem)
228  ELSEIF(typexn(1:1).EQ.'Q') THEN
229  CALL ov( 'X=Y ', x=xm(1,1), y=xn(1,2), dim1=nelem)
230  CALL ov( 'X=Y ', x=xm(1,2), y=xn(1,1), dim1=nelem)
231  ELSEIF(typexn(1:1).NE.'0') THEN
232  WRITE(lu,11) typexn(1:1)
233  CALL plante(1)
234  stop
235  ENDIF
236  typdim(1:1)=typdin(1:1)
237  typexm(1:1)=typexn(1:1)
238 !
239 !-----------------------------------------------------------------------
240 !
241  ELSEIF(op(3:8).EQ.'M+TN ') THEN
242 !
243 ! THE CASE WHERE N IS SYMMETRICAL HAS ALREADY BEEN TREATED
244 !
245  CALL ov('X=X+Y ', x=dm, y=dn, dim1=ndiag)
246 !
247  IF(typexm(1:1).EQ.'Q') THEN
248  CALL ov('X=X+Y ', x=xm(1,1), y=xn(1,2), dim1=nelem)
249  CALL ov('X=X+Y ', x=xm(1,2), y=xn(1,1), dim1=nelem)
250  ELSEIF(typexn(1:1).NE.'0') THEN
251  WRITE(lu,11) typexn(1:1)
252  CALL plante(1)
253  stop
254  ENDIF
255  typdim(1:1)=typdin(1:1)
256  typexm(1:1)=typexn(1:1)
257 !
258 !-----------------------------------------------------------------------
259 !
260  ELSEIF(op(3:8).EQ.'MD ') THEN
261 !
262 ! DIAGONAL TERMS
263 !
264  IF(typdim(1:1).EQ.'Q') THEN
265  CALL ov('X=XY ', x=dm, y=d, dim1=ndiag)
266  ELSEIF(typdim(1:1).EQ.'I') THEN
267  CALL ov('X=Y ', x=dm, y=d, dim1=ndiag)
268  typdim(1:1)='Q'
269  ELSEIF(typdim(1:1).NE.'0') THEN
270  WRITE(lu,13) typdim(1:1)
271  CALL plante(1)
272  stop
273  ENDIF
274 !
275 ! EXTRADIAGONAL TERMS
276 !
277  IF(typexm(1:1).EQ.'Q') THEN
278 !
279  DO ielem = 1 , nelem
280 !
281  xm(ielem, 1) = xm(ielem, 1) * d(ikle(ielem,2))
282  xm(ielem, 2) = xm(ielem, 2) * d(ikle(ielem,1))
283 !
284  ENDDO ! IELEM
285 !
286  ELSEIF(typexm(1:1).EQ.'S') THEN
287  WRITE(lu,171)
288 171 FORMAT(1x,
289  & 'OM0101 (BIEF) : M=MD NOT AVAILABLE IF M SYMMETRIC')
290  CALL plante(1)
291  stop
292  ELSEIF(typexm(1:1).NE.'0') THEN
293  WRITE(lu,173) typexm(1:1)
294 173 FORMAT(1x,'OM0101 (BIEF) : TYPEXM NOT AVAILABLE : ',a1)
295  CALL plante(1)
296  stop
297  ENDIF
298 !
299 !-----------------------------------------------------------------------
300 !
301  ELSEIF(op(3:8).EQ.'DM ') THEN
302 !
303 ! DIAGONAL TERMS
304 !
305  IF(typdim(1:1).EQ.'Q') THEN
306  CALL ov('X=XY ', x=dm, y=d, dim1=ndiag)
307  ELSEIF(typdim(1:1).EQ.'I') THEN
308  CALL ov('X=Y ', x=dm, y=d, dim1=ndiag)
309  typdim(1:1)='Q'
310  ELSEIF(typdim(1:1).NE.'0') THEN
311  WRITE(lu,13) typdim(1:1)
312  CALL plante(1)
313  stop
314  ENDIF
315 !
316 ! EXTRADIAGONAL TERMS
317 !
318  IF(typexm(1:1).EQ.'Q') THEN
319 !
320  DO ielem = 1 , nelem
321 !
322  xm(ielem, 1) = xm(ielem, 1) * d(ikle(ielem,1))
323  xm(ielem, 2) = xm(ielem, 2) * d(ikle(ielem,2))
324 !
325  ENDDO ! IELEM
326 !
327  ELSEIF(typexm(1:1).EQ.'S') THEN
328  WRITE(lu,181)
329 181 FORMAT(1x,
330  & 'OM0101 (BIEF) : M=MD NOT AVAILABLE IF M SYMMETRIC')
331  CALL plante(1)
332  stop
333  ELSEIF(typexm(1:1).NE.'0') THEN
334  WRITE(lu,173) typexm(1:1)
335  CALL plante(1)
336  stop
337  ENDIF
338 !
339 !-----------------------------------------------------------------------
340 !
341  ELSEIF(op(3:8).EQ.'DMD ') THEN
342 !
343 ! DIAGONAL TERMS
344 !
345  IF(typdim(1:1).EQ.'Q') THEN
346  CALL ov('X=XY ', x=dm, y=d, dim1=ndiag)
347  CALL ov('X=XY ', x=dm, y=d, dim1=ndiag)
348  ELSEIF(typdim(1:1).EQ.'I') THEN
349  CALL ov('X=YZ ', x=dm, y=d, z=d, dim1=ndiag)
350  typdim(1:1)='Q'
351  ELSEIF(typdim(1:1).NE.'0') THEN
352  WRITE(lu,13) typdim(1:1)
353 13 FORMAT(1x,'OM0101 (BIEF) : TYPDIM UNKNOWN :',a1)
354  CALL plante(1)
355  stop
356  ENDIF
357 !
358 ! EXTRADIAGONAL TERMS
359 !
360  IF(typexm(1:1).EQ.'S') THEN
361 !
362  DO ielem = 1 , nelem
363  xm(ielem,1)=xm(ielem,1)* d(ikle(ielem,2)) * d(ikle(ielem,1))
364  ENDDO ! IELEM
365 !
366  ELSEIF(typexm(1:1).EQ.'Q') THEN
367 !
368  DO ielem = 1 , nelem
369  xm(ielem,1)=xm(ielem,1)* d(ikle(ielem,2)) * d(ikle(ielem,1))
370  xm(ielem,2)=xm(ielem,2)* d(ikle(ielem,1)) * d(ikle(ielem,2))
371  ENDDO ! IELEM
372 !
373  ELSEIF(typexm(1:1).NE.'0') THEN
374  WRITE(lu,21) typexm(1:1)
375 21 FORMAT(1x,'OM0101 (BIEF) : TYPEXM UNKNOWN :',a1)
376  CALL plante(1)
377  stop
378  ENDIF
379 !
380 !-----------------------------------------------------------------------
381 !
382  ELSEIF(op(3:8).EQ.'X(M) ') THEN
383 !
384  IF(typexm(1:1).EQ.'S') THEN
385  CALL ov('X=Y ', x=xm(1,2), y=xm(1,1), dim1=nelem)
386  ELSEIF(typexm(1:1).NE.'0') THEN
387  WRITE(lu,811) typexm(1:1)
388 811 FORMAT(1x,'OM0101 (BIEF) : MATRIX ALREADY NON SYMMETRICAL: ',
389  & a1)
390  CALL plante(1)
391  stop
392  ENDIF
393  typexm(1:1)='Q'
394 !
395 !-----------------------------------------------------------------------
396 !
397  ELSEIF(op(3:8).EQ.'MSK(M)') THEN
398 !
399  IF(typexm(1:1).EQ.'S') THEN
400  j = 1
401  ELSEIF(typexm(1:1).EQ.'Q') THEN
402  j = 2
403  ELSEIF(typexm(1:1).EQ.'0') THEN
404  j = 0
405  ELSE
406 ! J=0 TO AVOID A STUPID COMPILER WARNING
407  j=0
408  WRITE(lu,173) typexm(1:1)
409  CALL plante(1)
410  stop
411  ENDIF
412 !
413  IF(j.GT.0) THEN
414  DO i = 1,j
415  CALL ov('X=XY ', x=xm(1,i), y=d, dim1=nelem)
416  ENDDO
417  ENDIF
418 !
419 !-----------------------------------------------------------------------
420 !
421  ELSE
422 !
423  WRITE(lu,41) op
424 41 FORMAT(1x,'OM0101 (BIEF) : UNKNOWN OPERATION : ',a8)
425  CALL plante(1)
426  stop
427 !
428  ENDIF
429 !
430 !-----------------------------------------------------------------------
431 !
432  RETURN
433  END
ov
subroutine ov(OP, X, Y, Z, C, DIM1)
Definition: ov.f:7
declarations_special
Definition: declarations_special.F:3
om0101
subroutine om0101(OP, DM, TYPDIM, XM, TYPEXM, DN, TYPDIN, XN, TYPEXN, D, C, IKLE, NELEM, NELMAX, NDIAG)
Definition: om0101.f:8
declarations_special::lu
integer lu
Definition: declarations_special.F:45
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:19 for The TELEMAC-MASCARET system by doxygen 1.8.13