The TELEMAC-MASCARET system  trunk
om1111.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE om1111
3 ! *****************
4 !
5  &(op , dm,typdim,xm,typexm, dn,typdin,xn,typexn, d,c,
6  & ikle,nelem,nelmax,ndiag,dm_err, dn_err, d_err)
7 !
8 !***********************************************************************
9 ! BIEF V6P1 21/08/2010
10 !***********************************************************************
11 !
12 !brief OPERATIONS ON MATRICES WITH P1 TRIANGLE.
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=CM ' : MULTIPLIES M BY C
25 !+ OP = 'M=M+CN ' : ADDS CN TO M
26 !+ OP = 'M=TN ' : COPIES TRANSPOSE OF N IN M
27 !+ OP = 'M=M+TN ' : ADDS TRANSPOSE(N) TO M
28 !+ OP = 'M=M+CTN ' : ADDS C TRANSPOSE(N) TO M
29 !+ OP = 'M=M+N ' : ADDS N TO M
30 !+ OP = 'M=MD ' : M X D
31 !+ OP = 'M=DM ' : D X M
32 !+ OP = 'M=M-ND ' : SUBTRACTS N X D TO M
33 !+ OP = 'M=M-DN ' : SUBTRACTS D X N TO M
34 !+ OP = 'M=DMD ' : D X M X D
35 !+ OP = 'M=0 ' : SETS M TO 0
36 !+ OP = 'M=X(M) ' : NOT SYMMETRICAL FORM OF M
37 !+ (OLD MATSNS)
38 !+ OP = 'M=MSK(M)' : MASKS M EXTRADIAGONAL TERMS
39 !+ (OLD MASKEX)
40 !+ THE MASK IS TAKEN FROM D
41 !+ OP = 'M=M+D ' : ADDS D TO M
42 !
43 !code
44 !+ CONVENTION FOR THE STORAGE OF EXTRA-DIAGONAL TERMS:
45 !+
46 !+ XM(IELEM,1) ----> M(1,2)
47 !+ XM(IELEM,2) ----> M(1,3)
48 !+ XM(IELEM,3) ----> M(2,3)
49 !+ XM(IELEM,4) ----> M(2,1)
50 !+ XM(IELEM,5) ----> M(3,1)
51 !+ XM(IELEM,6) ----> M(3,2)
52 !
53 !history J-M HERVOUET (LNHE) ; F LEPEINTRE (LNH)
54 !+ 05/02/91
55 !+ V5P1
56 !+
57 !
58 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
59 !+ 13/07/2010
60 !+ V6P0
61 !+ Translation of French comments within the FORTRAN sources into
62 !+ English comments
63 !
64 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
65 !+ 21/08/2010
66 !+ V6P0
67 !+ Creation of DOXYGEN tags for automated documentation and
68 !+ cross-referencing of the FORTRAN sources
69 !history R.NHEILI (Univerte de Perpignan, DALI)
70 !+ 24/02/2016
71 !+ V7P3
72 !+ ADD MODASS=3
73 !
74 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
75 !| C |-->| A GIVEN CONSTANT USED IN OPERATION OP
76 !| D |-->| A DIAGONAL MATRIX
77 !| DM |<->| DIAGONAL OF M
78 !| DN |-->| DIAGONAL OF N
79 !| IKLE |-->| CONNECTIVITY TABLE.
80 !| NDIAG |-->| NUMBER OF TERMS IN THE DIAGONAL
81 !| NELEM |-->| NUMBER OF ELEMENTS
82 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
83 !| OP |-->| OPERATION TO BE DONE (SEE ABOVE)
84 !| TYPDIM |<->| TYPE OF DIAGONAL OF M:
85 !| | | TYPDIM = 'Q' : ANY VALUE
86 !| | | TYPDIM = 'I' : IDENTITY
87 !| | | TYPDIM = '0' : ZERO
88 !| TYPDIN |<->| TYPE OF DIAGONAL OF N:
89 !| | | TYPDIN = 'Q' : ANY VALUE
90 !| | | TYPDIN = 'I' : IDENTITY
91 !| | | TYPDIN = '0' : ZERO
92 !| TYPEXM |-->| TYPE OF OFF-DIAGONAL TERMS OF M:
93 !| | | TYPEXM = 'Q' : ANY VALUE
94 !| | | TYPEXM = 'S' : SYMMETRIC
95 !| | | TYPEXM = '0' : ZERO
96 !| TYPEXN |-->| TYPE OF OFF-DIAGONAL TERMS OF N:
97 !| | | TYPEXN = 'Q' : ANY VALUE
98 !| | | TYPEXN = 'S' : SYMMETRIC
99 !| | | TYPEXN = '0' : ZERO
100 !| XM |-->| OFF-DIAGONAL TERMS OF M
101 !| XN |-->| OFF-DIAGONAL TERMS OF N
102 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
103 !
104  USE bief, ex_om1111 => om1111
105  USE declarations_telemac, ONLY : modass
106 !
107  USE declarations_special
108  IMPLICIT NONE
109 !
110 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
111 !
112  INTEGER, INTENT(IN) :: NELEM,NELMAX,NDIAG
113  INTEGER, INTENT(IN) :: IKLE(nelmax,3)
114  CHARACTER(LEN=8), INTENT(IN) :: OP
115  DOUBLE PRECISION, INTENT(IN) :: DN(*),D(*),XN(nelmax,*)
116  DOUBLE PRECISION, INTENT(INOUT) :: DM(*),XM(nelmax,*)
117  CHARACTER(LEN=1), INTENT(INOUT) :: TYPDIM,TYPEXM,TYPDIN,TYPEXN
118  DOUBLE PRECISION, INTENT(IN) :: C
119  DOUBLE PRECISION,OPTIONAL, INTENT(INOUT) :: DM_ERR(*)
120  DOUBLE PRECISION,OPTIONAL, INTENT(IN) :: DN_ERR(*),D_ERR(*)
121 !
122 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
123 !
124  INTEGER IELEM,I,J
125 !
126  DOUBLE PRECISION Z(1)
127 !
128 !-----------------------------------------------------------------------
129 !
130  IF(op(3:8).EQ.'N ') THEN
131 !
132  IF(typdin(1:1).EQ.'Q') THEN
133  IF ( modass .EQ.1) THEN
134  CALL ov('X=Y ', x=dm, y=dn, dim1=ndiag)
135  ELSEIF (modass .EQ. 3) THEN
136  CALL ov_comp( 'X=Y ' , dm , dn , z , c , ndiag,
137  & x_err=dm_err, y_err=dn_err )
138  ENDIF
139 !
140  ELSEIF(typdin(1:1).EQ.'I'.OR.typdin(1:1).EQ.'0') THEN
141 ! NOTHING TO DO, ONLY NEEDS TO COPY TYPDIN
142  ELSE
143  WRITE(lu,6) typdin(1:1)
144 6 FORMAT(1x,'OM1111 (BIEF) : TYPDIN UNKNOWN :',a1)
145  CALL plante(1)
146  stop
147  ENDIF
148  typdim(1:1)=typdin(1:1)
149 !
150  IF(typexn(1:1).EQ.'S') THEN
151  CALL ov('X=Y ', x=xm(1,1), y=xn(1,1), dim1=nelem)
152  CALL ov('X=Y ', x=xm(1,2), y=xn(1,2), dim1=nelem)
153  CALL ov('X=Y ', x=xm(1,3), y=xn(1,3), dim1=nelem)
154  ELSEIF(typexn(1:1).EQ.'Q') THEN
155  CALL ov('X=Y ', x=xm(1,1), y=xn(1,1), dim1=nelem)
156  CALL ov('X=Y ', x=xm(1,2), y=xn(1,2), dim1=nelem)
157  CALL ov('X=Y ', x=xm(1,3), y=xn(1,3), dim1=nelem)
158  CALL ov('X=Y ', x=xm(1,4), y=xn(1,4), dim1=nelem)
159  CALL ov('X=Y ', x=xm(1,5), y=xn(1,5), dim1=nelem)
160  CALL ov('X=Y ', x=xm(1,6), y=xn(1,6), dim1=nelem)
161  ELSEIF(typexn(1:1).NE.'0') THEN
162  WRITE(lu,11) typexn(1:1)
163 11 FORMAT(1x,'OM1111 (BIEF) : TYPEXN UNKNOWN :',a1)
164  CALL plante(1)
165  stop
166  ENDIF
167  typexm(1:1)=typexn(1:1)
168 !
169 !-----------------------------------------------------------------------
170 !
171  ELSEIF(op(3:8).EQ.'CN ') THEN
172 !
173  CALL ov('X=CY ', x=dm, y=dn, c=c, dim1=ndiag)
174 !
175  IF(typexn(1:1).EQ.'S') THEN
176  CALL ov('X=CY ', x=xm(1,1), y=xn(1,1), c=c, dim1=nelem)
177  CALL ov('X=CY ', x=xm(1,2), y=xn(1,2), c=c, dim1=nelem)
178  CALL ov('X=CY ', x=xm(1,3), y=xn(1,3), c=c, dim1=nelem)
179  ELSEIF(typexn(1:1).EQ.'Q') THEN
180  CALL ov('X=CY ', x=xm(1,1), y=xn(1,1), c=c, dim1=nelem)
181  CALL ov('X=CY ', x=xm(1,2), y=xn(1,2), c=c, dim1=nelem)
182  CALL ov('X=CY ', x=xm(1,3), y=xn(1,3), c=c, dim1=nelem)
183  CALL ov('X=CY ', x=xm(1,4), y=xn(1,4), c=c, dim1=nelem)
184  CALL ov('X=CY ', x=xm(1,5), y=xn(1,5), c=c, dim1=nelem)
185  CALL ov('X=CY ', x=xm(1,6), y=xn(1,6), c=c, dim1=nelem)
186  ELSEIF(typexn(1:1).NE.'0') THEN
187  WRITE(lu,11) typexn(1:1)
188  CALL plante(1)
189  stop
190  ENDIF
191 !
192  typdim(1:1)=typdin(1:1)
193  typexm(1:1)=typexn(1:1)
194 !
195 !-----------------------------------------------------------------------
196 !
197  ELSEIF(op(3:8).EQ.'CM ') THEN
198 !
199  CALL ov('X=CX ', x=dm, c=c, dim1=ndiag)
200 !
201  IF(typexm(1:1).EQ.'S') THEN
202  CALL ov('X=CX ', x=xm(1,1), c=c, dim1=nelem)
203  CALL ov('X=CX ', x=xm(1,2), c=c, dim1=nelem)
204  CALL ov('X=CX ', x=xm(1,3), c=c, dim1=nelem)
205  ELSEIF(typexm(1:1).EQ.'Q') THEN
206  CALL ov('X=CX ', x=xm(1,1), c=c, dim1=nelem)
207  CALL ov('X=CX ', x=xm(1,2), c=c, dim1=nelem)
208  CALL ov('X=CX ', x=xm(1,3), c=c, dim1=nelem)
209  CALL ov('X=CX ', x=xm(1,4), c=c, dim1=nelem)
210  CALL ov('X=CX ', x=xm(1,5), c=c, dim1=nelem)
211  CALL ov('X=CX ', x=xm(1,6), c=c, dim1=nelem)
212  ELSEIF(typexm(1:1).NE.'0') THEN
213  WRITE(lu,11) typexm(1:1)
214  CALL plante(1)
215  stop
216  ENDIF
217 !
218 !-----------------------------------------------------------------------
219 !
220  ELSEIF(op(3:8).EQ.'M+CN ' .OR.
221  & (op(3:8).EQ.'M+CTN ').AND.typexn(1:1).NE.'Q') THEN
222 !
223  IF(typdin(1:1).EQ.'I') THEN
224  CALL ov('X=X+C ', x=dm, c=c, dim1=ndiag)
225  ELSEIF(typdin(1:1).NE.'0') THEN
226  CALL ov('X=X+CY ', x=dm, y=dn, c=c, dim1=ndiag)
227  ENDIF
228 !
229  IF(typexn(1:1).EQ.'S') THEN
230  CALL ov('X=X+CY ', x=xm(1,1), y=xn(1,1), c=c, dim1=nelem)
231  CALL ov('X=X+CY ', x=xm(1,2), y=xn(1,2), c=c, dim1=nelem)
232  CALL ov('X=X+CY ', x=xm(1,3), y=xn(1,3), c=c, dim1=nelem)
233  IF(typexm(1:1).EQ.'Q') THEN
234  CALL ov('X=X+CY ', x=xm(1,4), y=xn(1,1), c=c, dim1=nelem)
235  CALL ov('X=X+CY ', x=xm(1,5), y=xn(1,2), c=c, dim1=nelem)
236  CALL ov('X=X+CY ', x=xm(1,6), y=xn(1,3), c=c, dim1=nelem)
237  ENDIF
238  ELSEIF(typexn(1:1).EQ.'Q') THEN
239  IF(typexm(1:1).NE.'Q') THEN
240  WRITE(lu,98) typexm(1:1),op(1:8),typexn(1:1)
241 98 FORMAT(1x,'OM1111 (BIEF) : TYPEXM = ',a1,' DOES NOT GO',
242  & /,1x,'FOR THE OPERATION : ',a8,' WITH TYPEXN = ',a1)
243  CALL plante(1)
244  stop
245  ENDIF
246  CALL ov('X=X+CY ', x=xm(1,1), y=xn(1,1), c=c, dim1=nelem)
247  CALL ov('X=X+CY ', x=xm(1,2), y=xn(1,2), c=c, dim1=nelem)
248  CALL ov('X=X+CY ', x=xm(1,3), y=xn(1,3), c=c, dim1=nelem)
249  CALL ov('X=X+CY ', x=xm(1,4), y=xn(1,4), c=c, dim1=nelem)
250  CALL ov('X=X+CY ', x=xm(1,5), y=xn(1,5), c=c, dim1=nelem)
251  CALL ov('X=X+CY ', x=xm(1,6), y=xn(1,6), c=c, dim1=nelem)
252  ELSEIF(typexn(1:1).NE.'0') THEN
253  WRITE(lu,11) typexn(1:1)
254  CALL plante(1)
255  stop
256  ENDIF
257 !
258 !-----------------------------------------------------------------------
259 !
260  ELSEIF(op(3:8).EQ.'M+CTN ') THEN
261 !
262 ! THE CASES WHERE N IS SYMMETRICAL ARE TREATED WITH M=M+CN
263 !
264  CALL ov('X=X+CY ', x=dm, y=dn, c=c, dim1=ndiag)
265 !
266  IF(typexn(1:1).EQ.'Q') THEN
267  IF(typexm(1:1).NE.'Q') THEN
268  WRITE(lu,98) typexm(1:1),op(1:8),typexn(1:1)
269  CALL plante(1)
270  stop
271  ENDIF
272  CALL ov('X=X+CY ', x=xm(1,1), y=xn(1,4), c=c, dim1=nelem)
273  CALL ov('X=X+CY ', x=xm(1,2), y=xn(1,5), c=c, dim1=nelem)
274  CALL ov('X=X+CY ', x=xm(1,3), y=xn(1,6), c=c, dim1=nelem)
275  CALL ov('X=X+CY ', x=xm(1,4), y=xn(1,1), c=c, dim1=nelem)
276  CALL ov('X=X+CY ', x=xm(1,5), y=xn(1,2), c=c, dim1=nelem)
277  CALL ov('X=X+CY ', x=xm(1,6), y=xn(1,3), c=c, dim1=nelem)
278  ELSE
279  WRITE(lu,11) typexn(1:1)
280  CALL plante(1)
281  stop
282  ENDIF
283 !
284 !-----------------------------------------------------------------------
285 !
286  ELSEIF(op(3:8).EQ.'TN ') THEN
287 !
288  CALL ov('X=Y ', x=dm, y=dn, dim1=ndiag)
289 !
290  IF(typexn(1:1).EQ.'S') THEN
291  CALL ov('X=Y ', x=xm(1,1), y=xn(1,1), dim1=nelem)
292  CALL ov('X=Y ', x=xm(1,2), y=xn(1,2), dim1=nelem)
293  CALL ov('X=Y ', x=xm(1,3), y=xn(1,3), dim1=nelem)
294  ELSEIF(typexn(1:1).EQ.'Q') THEN
295  CALL ov('X=Y ', x=xm(1,1), y=xn(1,4), dim1=nelem)
296  CALL ov('X=Y ', x=xm(1,2), y=xn(1,5), dim1=nelem)
297  CALL ov('X=Y ', x=xm(1,3), y=xn(1,6), dim1=nelem)
298  CALL ov('X=Y ', x=xm(1,4), y=xn(1,1), dim1=nelem)
299  CALL ov('X=Y ', x=xm(1,5), y=xn(1,2), dim1=nelem)
300  CALL ov('X=Y ', x=xm(1,6), y=xn(1,3), dim1=nelem)
301  ELSEIF(typexn(1:1).NE.'0') THEN
302  WRITE(lu,11) typexn(1:1)
303  CALL plante(1)
304  stop
305  ENDIF
306  typdim(1:1)=typdin(1:1)
307  typexm(1:1)=typexn(1:1)
308 !
309 !-----------------------------------------------------------------------
310 !
311  ELSEIF(op(3:8).EQ.'M+N '.OR.
312  & (op(3:8).EQ.'M+TN ').AND.typexn(1:1).NE.'Q') THEN
313 !
314  IF ( modass .EQ.1) THEN
315  CALL ov('X=X+Y ', x=dm, y=dn, dim1=ndiag)
316  ELSEIF (modass .EQ. 3) THEN
317  CALL ov_comp( 'X=X+Y ' , dm , dn , z , c , ndiag,
318  & x_err=dm_err, y_err=dn_err )
319  ENDIF
320 !
321  IF(typexn(1:1).EQ.'S') THEN
322  CALL ov('X=X+Y ', x=xm(1,1), y=xn(1,1), dim1=nelem)
323  CALL ov('X=X+Y ', x=xm(1,2), y=xn(1,2), dim1=nelem)
324  CALL ov('X=X+Y ', x=xm(1,3), y=xn(1,3), dim1=nelem)
325  IF(typexm(1:1).EQ.'Q') THEN
326  CALL ov('X=X+Y ', x=xm(1,4), y=xn(1,1), dim1=nelem)
327  CALL ov('X=X+Y ', x=xm(1,5), y=xn(1,2), dim1=nelem)
328  CALL ov('X=X+Y ', x=xm(1,6), y=xn(1,3), dim1=nelem)
329  ENDIF
330  ELSEIF(typexn(1:1).EQ.'Q') THEN
331  IF(typexm(1:1).NE.'Q') THEN
332  WRITE(lu,98) typexm(1:1),op(1:8),typexn(1:1)
333  CALL plante(1)
334  stop
335  ENDIF
336  CALL ov('X=X+Y ', x=xm(1,1), y=xn(1,1), dim1=nelem)
337  CALL ov('X=X+Y ', x=xm(1,2), y=xn(1,2), dim1=nelem)
338  CALL ov('X=X+Y ', x=xm(1,3), y=xn(1,3), dim1=nelem)
339  CALL ov('X=X+Y ', x=xm(1,4), y=xn(1,4), dim1=nelem)
340  CALL ov('X=X+Y ', x=xm(1,5), y=xn(1,5), dim1=nelem)
341  CALL ov('X=X+Y ', x=xm(1,6), y=xn(1,6), dim1=nelem)
342  ELSEIF(typexn(1:1).NE.'0') THEN
343  WRITE(lu,11) typexn(1:1)
344  CALL plante(1)
345  stop
346  ENDIF
347 !
348 !-----------------------------------------------------------------------
349 !
350  ELSEIF(op(3:8).EQ.'M+TN ') THEN
351 !
352 ! THE CASE WHERE N IS SYMMETRICAL HAS ALREADY BEEN TREATED
353 !
354  CALL ov('X=X+Y ', x=dm, y=dn, dim1=ndiag)
355 !
356  IF(typexm(1:1).EQ.'Q') THEN
357  CALL ov('X=X+Y ', x=xm(1,1), y=xn(1,4), dim1=nelem)
358  CALL ov('X=X+Y ', x=xm(1,2), y=xn(1,5), dim1=nelem)
359  CALL ov('X=X+Y ', x=xm(1,3), y=xn(1,6), dim1=nelem)
360  CALL ov('X=X+Y ', x=xm(1,4), y=xn(1,1), dim1=nelem)
361  CALL ov('X=X+Y ', x=xm(1,5), y=xn(1,2), dim1=nelem)
362  CALL ov('X=X+Y ', x=xm(1,6), y=xn(1,3), dim1=nelem)
363  ELSEIF(typexn(1:1).NE.'0') THEN
364  WRITE(lu,11) typexn(1:1)
365  CALL plante(1)
366  stop
367  ENDIF
368  typdim(1:1)=typdin(1:1)
369  typexm(1:1)=typexn(1:1)
370 !
371 !-----------------------------------------------------------------------
372 !
373  ELSEIF(op(3:8).EQ.'MD ') THEN
374 !
375 ! DIAGONAL TERMS
376 !
377  IF(typdim(1:1).EQ.'Q') THEN
378  CALL ov('X=XY ', x=dm, y=d, dim1=ndiag)
379  ELSEIF(typdim(1:1).EQ.'I') THEN
380  CALL ov('X=Y ', x=dm, y=d, dim1=ndiag)
381  typdim(1:1)='Q'
382  ELSEIF(typdim(1:1).NE.'0') THEN
383  WRITE(lu,13) typdim(1:1)
384  CALL plante(1)
385  stop
386  ENDIF
387 !
388 ! EXTRADIAGONAL TERMS
389 !
390  IF(typexm(1:1).EQ.'Q') THEN
391 !
392  DO ielem = 1 , nelem
393 !
394  xm(ielem, 1) = xm(ielem, 1) * d(ikle(ielem,2))
395  xm(ielem, 2) = xm(ielem, 2) * d(ikle(ielem,3))
396  xm(ielem, 3) = xm(ielem, 3) * d(ikle(ielem,3))
397  xm(ielem, 4) = xm(ielem, 4) * d(ikle(ielem,1))
398  xm(ielem, 5) = xm(ielem, 5) * d(ikle(ielem,1))
399  xm(ielem, 6) = xm(ielem, 6) * d(ikle(ielem,2))
400 !
401  ENDDO
402 !
403  ELSEIF(typexm(1:1).EQ.'S') THEN
404  WRITE(lu,171)
405 171 FORMAT(1x,'OM1111 (BIEF) : M=MD , M MUST BE NON-SYMMETRIC')
406  CALL plante(1)
407  stop
408  ELSEIF(typexm(1:1).NE.'0') THEN
409  WRITE(lu,173) typexm(1:1)
410 173 FORMAT(1x,'OM1111 (BIEF) : TYPEXM NOT AVAILABLE : ',a1)
411  CALL plante(1)
412  stop
413  ENDIF
414 !
415 !-----------------------------------------------------------------------
416 !
417  ELSEIF(op(3:8).EQ.'DM ') THEN
418 !
419 ! DIAGONAL TERMS
420 !
421  IF(typdim(1:1).EQ.'Q') THEN
422  CALL ov('X=XY ', x=dm, y=d, dim1=ndiag)
423  ELSEIF(typdim(1:1).EQ.'I') THEN
424  CALL ov('X=Y ', x=dm, y=d, dim1=ndiag)
425  typdim(1:1)='Q'
426  ELSEIF(typdim(1:1).NE.'0') THEN
427  WRITE(lu,13) typdim(1:1)
428  CALL plante(1)
429  stop
430  ENDIF
431 !
432 ! EXTRADIAGONAL TERMS
433 !
434  IF(typexm(1:1).EQ.'Q') THEN
435 !
436  DO ielem = 1 , nelem
437 !
438  xm(ielem, 1) = xm(ielem, 1) * d(ikle(ielem,1))
439  xm(ielem, 2) = xm(ielem, 2) * d(ikle(ielem,1))
440  xm(ielem, 3) = xm(ielem, 3) * d(ikle(ielem,2))
441  xm(ielem, 4) = xm(ielem, 4) * d(ikle(ielem,2))
442  xm(ielem, 5) = xm(ielem, 5) * d(ikle(ielem,3))
443  xm(ielem, 6) = xm(ielem, 6) * d(ikle(ielem,3))
444 !
445  ENDDO
446 !
447  ELSEIF(typexm(1:1).EQ.'S') THEN
448  WRITE(lu,181)
449 181 FORMAT(1x,
450  & 'OM1111 (BIEF) : M=MD NOT AVAILABLE IF M SYMMETRIC')
451  CALL plante(1)
452  stop
453  ELSEIF(typexm(1:1).NE.'0') THEN
454  WRITE(lu,173) typexm(1:1)
455  CALL plante(1)
456  stop
457  ENDIF
458 !
459 !-----------------------------------------------------------------------
460 !
461  ELSEIF(op(3:8).EQ.'M-DN ') THEN
462 !
463 ! DIAGONAL TERMS
464 !
465  IF(typdim(1:1).EQ.'Q') THEN
466  CALL ov('X=X-YZ ', x=dm, y=dn, z=d, dim1=ndiag)
467  ELSEIF(typdim(1:1).NE.'0') THEN
468  WRITE(lu,13) typdim(1:1)
469  CALL plante(1)
470  stop
471  ENDIF
472 !
473 ! EXTRADIAGONAL TERMS
474 !
475  IF(typexm(1:1).EQ.'Q') THEN
476  IF(typexn(1:1).EQ.'Q') THEN
477  DO ielem = 1 , nelem
478  xm(ielem, 1) = xm(ielem,1) - xn(ielem, 1) * d(ikle(ielem,1))
479  xm(ielem, 2) = xm(ielem,2) - xn(ielem, 2) * d(ikle(ielem,1))
480  xm(ielem, 3) = xm(ielem,3) - xn(ielem, 3) * d(ikle(ielem,2))
481  xm(ielem, 4) = xm(ielem,4) - xn(ielem, 4) * d(ikle(ielem,2))
482  xm(ielem, 5) = xm(ielem,5) - xn(ielem, 5) * d(ikle(ielem,3))
483  xm(ielem, 6) = xm(ielem,6) - xn(ielem, 6) * d(ikle(ielem,3))
484  ENDDO ! IELEM
485  ELSEIF(typexn(1:1).EQ.'S') THEN
486  DO ielem = 1 , nelem
487  xm(ielem, 1) = xm(ielem,1) - xn(ielem, 1) * d(ikle(ielem,1))
488  xm(ielem, 2) = xm(ielem,2) - xn(ielem, 2) * d(ikle(ielem,1))
489  xm(ielem, 3) = xm(ielem,3) - xn(ielem, 3) * d(ikle(ielem,2))
490  xm(ielem, 4) = xm(ielem,4) - xn(ielem, 1) * d(ikle(ielem,2))
491  xm(ielem, 5) = xm(ielem,5) - xn(ielem, 2) * d(ikle(ielem,3))
492  xm(ielem, 6) = xm(ielem,6) - xn(ielem, 3) * d(ikle(ielem,3))
493  ENDDO ! IELEM
494  ELSEIF(typexn(1:1).NE.'0') THEN
495  WRITE(lu,11) typexn(1:1)
496  CALL plante(1)
497  stop
498  ENDIF
499  ELSE
500  WRITE(lu,173) typexm(1:1)
501  CALL plante(1)
502  stop
503  ENDIF
504 !
505 !-----------------------------------------------------------------------
506 !
507  ELSEIF(op(3:8).EQ.'M-ND ') THEN
508 !
509 ! DIAGONAL TERMS
510 !
511  IF(typdim(1:1).EQ.'Q') THEN
512  CALL ov('X=X-YZ ', x=dm, y=dn, z=d, dim1=ndiag)
513  ELSEIF(typdim(1:1).NE.'0') THEN
514  WRITE(lu,13) typdim(1:1)
515  CALL plante(1)
516  stop
517  ENDIF
518 !
519 ! EXTRADIAGONAL TERMS
520 !
521  IF(typexm(1:1).EQ.'Q') THEN
522  IF(typexn(1:1).EQ.'Q') THEN
523  DO ielem = 1 , nelem
524  xm(ielem, 1) = xm(ielem,1) - xn(ielem, 1) * d(ikle(ielem,2))
525  xm(ielem, 2) = xm(ielem,2) - xn(ielem, 2) * d(ikle(ielem,3))
526  xm(ielem, 3) = xm(ielem,3) - xn(ielem, 3) * d(ikle(ielem,3))
527  xm(ielem, 4) = xm(ielem,4) - xn(ielem, 4) * d(ikle(ielem,1))
528  xm(ielem, 5) = xm(ielem,5) - xn(ielem, 5) * d(ikle(ielem,1))
529  xm(ielem, 6) = xm(ielem,6) - xn(ielem, 6) * d(ikle(ielem,2))
530  ENDDO ! IELEM
531  ELSEIF(typexn(1:1).EQ.'S') THEN
532  DO ielem = 1 , nelem
533  xm(ielem, 1) = xm(ielem,1) - xn(ielem, 1) * d(ikle(ielem,2))
534  xm(ielem, 2) = xm(ielem,2) - xn(ielem, 2) * d(ikle(ielem,3))
535  xm(ielem, 3) = xm(ielem,3) - xn(ielem, 3) * d(ikle(ielem,3))
536  xm(ielem, 4) = xm(ielem,4) - xn(ielem, 1) * d(ikle(ielem,1))
537  xm(ielem, 5) = xm(ielem,5) - xn(ielem, 2) * d(ikle(ielem,1))
538  xm(ielem, 6) = xm(ielem,6) - xn(ielem, 3) * d(ikle(ielem,2))
539  ENDDO ! IELEM
540  ELSEIF(typexn(1:1).NE.'0') THEN
541  WRITE(lu,11) typexn(1:1)
542  CALL plante(1)
543  stop
544  ENDIF
545  ELSE
546  WRITE(lu,173) typexm(1:1)
547  CALL plante(1)
548  stop
549  ENDIF
550 !
551 !-----------------------------------------------------------------------
552 !
553  ELSEIF(op(3:8).EQ.'DMD ') THEN
554 !
555 ! DIAGONAL TERMS
556 !
557  IF(typdim(1:1).EQ.'Q') THEN
558  IF (modass.EQ. 1)THEN
559  CALL ov('X=XY ', x=dm, y=d, dim1=ndiag)
560  CALL ov('X=XY ', x=dm, y=d, dim1=ndiag)
561  ELSEIF (modass.EQ. 3)THEN
562  CALL ov_comp( 'X=XY ' , dm , d , z , c , ndiag,
563  & x_err=dm_err, y_err= d_err )
564  CALL ov_comp( 'X=XY ' , dm , d , z , c , ndiag,
565  & x_err=dm_err, y_err= d_err )
566  ENDIF
567  ELSEIF(typdim(1:1).EQ.'I') THEN
568  CALL ov('X=YZ ', x=dm, y=d, z=d, dim1=ndiag)
569  typdim(1:1)='Q'
570  ELSEIF(typdim(1:1).NE.'0') THEN
571  WRITE(lu,13) typdim(1:1)
572 13 FORMAT(1x,'OM1111 (BIEF) : TYPDIM UNKNOWN :',a1)
573  CALL plante(1)
574  stop
575  ENDIF
576 !
577 ! EXTRADIAGONAL TERMS
578 !
579  IF(typexm(1:1).EQ.'S') THEN
580 !
581  DO ielem = 1 , nelem
582  xm(ielem,1)=xm(ielem,1)* d(ikle(ielem,2)) * d(ikle(ielem,1))
583  xm(ielem,2)=xm(ielem,2)* d(ikle(ielem,3)) * d(ikle(ielem,1))
584  xm(ielem,3)=xm(ielem,3)* d(ikle(ielem,3)) * d(ikle(ielem,2))
585  ENDDO
586 !
587  ELSEIF(typexm(1:1).EQ.'Q') THEN
588 !
589  DO ielem = 1 , nelem
590  xm(ielem,1)=xm(ielem,1)* d(ikle(ielem,2)) * d(ikle(ielem,1))
591  xm(ielem,2)=xm(ielem,2)* d(ikle(ielem,3)) * d(ikle(ielem,1))
592  xm(ielem,3)=xm(ielem,3)* d(ikle(ielem,3)) * d(ikle(ielem,2))
593  xm(ielem,4)=xm(ielem,4)* d(ikle(ielem,2)) * d(ikle(ielem,1))
594  xm(ielem,5)=xm(ielem,5)* d(ikle(ielem,3)) * d(ikle(ielem,1))
595  xm(ielem,6)=xm(ielem,6)* d(ikle(ielem,3)) * d(ikle(ielem,2))
596  ENDDO
597 !
598  ELSEIF(typexm(1:1).NE.'0') THEN
599  WRITE(lu,21) typexm(1:1)
600 21 FORMAT(1x,'OM1111 (BIEF) : TYPEXM UNKNOWN :',a1)
601  CALL plante(1)
602  stop
603  ENDIF
604 !
605 !-----------------------------------------------------------------------
606 !
607  ELSEIF(op(3:8).EQ.'M+D ') THEN
608 !
609  IF ( modass .EQ.1) THEN
610  CALL ov('X=X+Y ', x=dm, y=d, dim1=ndiag)
611  ELSEIF (modass .EQ. 3) THEN
612  CALL ov_comp( 'X=X+Y ' , dm , d , z , 0.d0, ndiag,
613  & x_err=dm_err, y_err= d_err )
614  ENDIF
615 ! HERE THERE IS A DOUBT ABOUT TYPDIM
616  typdim(1:1)='Q'
617 !
618 !-----------------------------------------------------------------------
619 !
620  ELSEIF(op(3:8).EQ.'0 ') THEN
621 !
622  CALL ov('X=C ', x=dm, c=0.d0, dim1=ndiag)
623 !
624  IF(typexm(1:1).EQ.'S') THEN
625  CALL ov('X=C ', x=xm(1,1), c=0.d0, dim1=nelem)
626  CALL ov('X=C ', x=xm(1,2), c=0.d0, dim1=nelem)
627  CALL ov('X=C ', x=xm(1,3), c=0.d0, dim1=nelem)
628  ELSEIF(typexm(1:1).EQ.'Q') THEN
629  CALL ov('X=C ', x=xm(1,1), c=0.d0, dim1=nelem)
630  CALL ov('X=C ', x=xm(1,2), c=0.d0, dim1=nelem)
631  CALL ov('X=C ', x=xm(1,3), c=0.d0, dim1=nelem)
632  CALL ov('X=C ', x=xm(1,4), c=0.d0, dim1=nelem)
633  CALL ov('X=C ', x=xm(1,5), c=0.d0, dim1=nelem)
634  CALL ov('X=C ', x=xm(1,6), c=0.d0, dim1=nelem)
635  ELSEIF(typexm(1:1).NE.'0') THEN
636  WRITE(lu,711) typexm(1:1)
637 711 FORMAT(1x,'OM1111 (BIEF) : TYPEXM UNKNOWN :',a1)
638  CALL plante(1)
639  stop
640  ENDIF
641 ! TYPDIM IS NOT CHANGED
642 ! TYPDIM(1:1)='0'
643 ! TYPEXM IS NOT CHANGED
644 ! TYPEXM(1:1)='0'
645 !-----------------------------------------------------------------------
646 !
647  ELSEIF(op(3:8).EQ.'X(M) ') THEN
648 !
649  IF(typexm(1:1).EQ.'S') THEN
650  CALL ov('X=Y ', x=xm(1,4), y=xm(1,1), dim1=nelem)
651  CALL ov('X=Y ', x=xm(1,5), y=xm(1,2), dim1=nelem)
652  CALL ov('X=Y ', x=xm(1,6), y=xm(1,3), dim1=nelem)
653  ELSEIF(typexm(1:1).NE.'0') THEN
654  WRITE(lu,811) typexm(1:1)
655 811 FORMAT(1x,'OM1111 (BIEF) : MATRIX ALREADY NON SYMMETRICAL: ',
656  & a1)
657  CALL plante(1)
658  stop
659  ENDIF
660  typexm(1:1)='Q'
661 !
662 !-----------------------------------------------------------------------
663 !
664  ELSEIF(op(3:8).EQ.'MSK(M)') THEN
665 !
666  IF(typexm(1:1).EQ.'S') THEN
667  j = 3
668  ELSEIF(typexm(1:1).EQ.'Q') THEN
669  j = 6
670  ELSEIF(typexm(1:1).EQ.'0') THEN
671  j = 0
672  ELSE
673  WRITE(lu,711) typexm
674  j=0
675  CALL plante(1)
676  stop
677  ENDIF
678 !
679  IF(j.GT.0) THEN
680  DO i = 1,j
681  CALL ov('X=XY ', x=xm(1,i), y=d, dim1=nelem)
682  ENDDO
683  ENDIF
684 !
685 !-----------------------------------------------------------------------
686 !
687  ELSE
688 !
689  WRITE(lu,41) op
690 41 FORMAT(1x,'OM1111 (BIEF) : UNKNOWN OPERATION : ',a8)
691  CALL plante(1)
692  stop
693 !
694  ENDIF
695 !
696 !-----------------------------------------------------------------------
697 !
698  RETURN
699  END
ov_comp
subroutine ov_comp(OP, X, Y, Z, C, NPOIN, X_ERR, Y_ERR, Z_ERR)
Definition: ov_comp.f:8
ov
subroutine ov(OP, X, Y, Z, C, DIM1)
Definition: ov.f:7
declarations_special
Definition: declarations_special.F:3
om1111
subroutine om1111(OP, DM, TYPDIM, XM, TYPEXM, DN, TYPDIN, XN, TYPEXN, D, C, IKLE, NELEM, NELMAX, NDIAG, DM_ERR, DN_ERR, D_ERR)
Definition: om1111.f:8
declarations_special::lu
integer lu
Definition: declarations_special.F:45
declarations_telemac
Definition: declarations_telemac.f:3
declarations_telemac::modass
integer modass
Definition: declarations_telemac.f:265
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:19 for The TELEMAC-MASCARET system by doxygen 1.8.13