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