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