The TELEMAC-MASCARET system  trunk
om1112.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE om1112
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.
15 !+
16 !+ (LINEAR,QUASIBUBBLE FOR EXAMPLE).
17 !code
18 !+ D: DIAGONAL MATRIX
19 !+ C: CONSTANT
20 !+
21 !+ OP IS A STRING OF 8 CHARACTERS, WHICH INDICATES THE OPERATION TO BE
22 !+ PERFORMED ON MATRICES M AND N, D AND C.
23 !+
24 !+ THE RESULT IS MATRIX M.
25 !+
26 !+ OP = 'M=N ' : COPIES N IN M
27 !+ OP = 'M=CN ' : MULTIPLIES N BY C
28 !+ OP = 'M=CM ' : MULTIPLIES M BY C
29 !+ OP = 'M=M+CN ' : ADDS CN TO M
30 !+ OP = 'M=M+N ' : ADDS N TO M
31 !+ OP = 'M=MD ' : M X D
32 !+ OP = 'M=DM ' : D X M
33 !+ OP = 'M=M+D ' : ADDS D TO M
34 !+ OP = 'M=MSK(M)' : MASKS M EXTRADIAGONAL TERMS
35 !+ (OLD MASKEX)
36 !+ THE MASK IS TAKEN FROM D
37 !
38 !code
39 !+ CONVENTION FOR THE STORAGE OF EXTRA-DIAGONAL TERMS:
40 !+
41 !+ XM(IELEM, 1) ----> M(1,2)
42 !+ XM(IELEM, 2) ----> M(1,3)
43 !+ XM(IELEM, 3) ----> M(1,4)
44 !+ XM(IELEM, 4) ----> M(2,1)
45 !+ XM(IELEM, 5) ----> M(2,3)
46 !+ XM(IELEM, 6) ----> M(2,4)
47 !+ XM(IELEM, 7) ----> M(3,1)
48 !+ XM(IELEM, 8) ----> M(3,2)
49 !+ XM(IELEM, 9) ----> M(3,4)
50 !
51 !history J-M HERVOUET (LNHE) ; F LEPEINTRE (LNH)
52 !+ 21/09/93
53 !+ V5P1
54 !+
55 !
56 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
57 !+ 13/07/2010
58 !+ V6P0
59 !+ Translation of French comments within the FORTRAN sources into
60 !+ English comments
61 !
62 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
63 !+ 21/08/2010
64 !+ V6P0
65 !+ Creation of DOXYGEN tags for automated documentation and
66 !+ cross-referencing of the FORTRAN sources
67 !
68 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
69 !| C |-->| A GIVEN CONSTANT USED IN OPERATION OP
70 !| D |-->| A DIAGONAL MATRIX
71 !| DM |<->| DIAGONAL OF M
72 !| DN |-->| DIAGONAL OF N
73 !| IKLE |-->| CONNECTIVITY TABLE.
74 !| NDIAG |-->| NUMBER OF TERMS IN THE DIAGONAL
75 !| NELEM |-->| NUMBER OF ELEMENTS
76 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
77 !| OP |-->| OPERATION TO BE DONE (SEE ABOVE)
78 !| TYPDIM |<->| TYPE OF DIAGONAL OF M:
79 !| | | TYPDIM = 'Q' : ANY VALUE
80 !| | | TYPDIM = 'I' : IDENTITY
81 !| | | TYPDIM = '0' : ZERO
82 !| TYPDIN |<->| TYPE OF DIAGONAL OF N:
83 !| | | TYPDIN = 'Q' : ANY VALUE
84 !| | | TYPDIN = 'I' : IDENTITY
85 !| | | TYPDIN = '0' : ZERO
86 !| TYPEXM |-->| TYPE OF OFF-DIAGONAL TERMS OF M:
87 !| | | TYPEXM = 'Q' : ANY VALUE
88 !| | | TYPEXM = 'S' : SYMMETRIC
89 !| | | TYPEXM = '0' : ZERO
90 !| TYPEXN |-->| TYPE OF OFF-DIAGONAL TERMS OF N:
91 !| | | TYPEXN = 'Q' : ANY VALUE
92 !| | | TYPEXN = 'S' : SYMMETRIC
93 !| | | TYPEXN = '0' : ZERO
94 !| XM |-->| OFF-DIAGONAL TERMS OF M
95 !| XN |-->| OFF-DIAGONAL TERMS OF N
96 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
97 !
98  USE bief, ex_om1112 => om1112
99 !
100  USE declarations_special
101  IMPLICIT NONE
102 !
103 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
104 !
105  INTEGER, INTENT(IN) :: NELEM,NELMAX,NDIAG
106  INTEGER, INTENT(IN) :: IKLE(nelmax,4)
107  CHARACTER(LEN=8), INTENT(IN) :: OP
108  DOUBLE PRECISION, INTENT(IN) :: DN(*),D(*),XN(nelmax,*)
109  DOUBLE PRECISION, INTENT(INOUT) :: DM(*),XM(nelmax,*)
110  CHARACTER(LEN=1), INTENT(INOUT) :: TYPDIM,TYPEXM,TYPDIN,TYPEXN
111  DOUBLE PRECISION, INTENT(IN) :: C
112 !
113 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
114 !
115  INTEGER I,J,IELEM
116 !
117 !-----------------------------------------------------------------------
118 !
119  IF(op(1:8).EQ.'M=N ') THEN
120 !
121  IF(typdin(1:1).EQ.'Q') THEN
122  CALL ov('X=Y ', x=dm, y=dn, dim1=ndiag)
123  ELSEIF(typdin(1:1).EQ.'I'.OR.typdin(1:1).EQ.'0') THEN
124 ! NOTHING TO DO, ONLY NEEDS TO COPY TYPDIN
125  ELSE
126  WRITE(lu,6) typdin(1:1)
127 6 FORMAT(1x,'OM1112 (BIEF) : TYPDIN UNKNOWN :',a1)
128  CALL plante(1)
129  stop
130  ENDIF
131  typdim(1:1)=typdin(1:1)
132 !
133  IF(typexn(1:1).EQ.'Q') THEN
134  DO i=1,9
135  CALL ov('X=Y ', x=xm(1,i), y=xn(1,i), dim1=nelem)
136  ENDDO ! I
137  ELSEIF(typexn(1:1).NE.'0') THEN
138  WRITE(lu,11) typexn(1:1)
139 11 FORMAT(1x,'OM1112 (BIEF) : TYPEXN UNKNOWN :',a1)
140  CALL plante(1)
141  stop
142  ENDIF
143 !
144  typexm(1:1)=typexn(1:1)
145 !
146 !-----------------------------------------------------------------------
147 !
148  ELSEIF(op(1:8).EQ.'M=CN ') THEN
149 !
150  CALL ov('X=CY ', x=dm, y=dn, c=c, dim1=ndiag)
151 !
152  IF(typexn(1:1).EQ.'Q') THEN
153  DO i=1,9
154  CALL ov('X=CY ', x=xm(1,i), y=xn(1,i), c=c, dim1=nelem)
155  ENDDO ! I
156  ELSEIF(typexn(1:1).NE.'0') THEN
157  WRITE(lu,11) typexn(1:1)
158  CALL plante(1)
159  stop
160  ENDIF
161 !
162  typdim(1:1)=typdin(1:1)
163  typexm(1:1)=typexn(1:1)
164 !
165 !-----------------------------------------------------------------------
166 !
167  ELSEIF(op(1:8).EQ.'M=CM ') THEN
168 !
169  CALL ov('X=CX ', x=dm, c=c, dim1=ndiag)
170 !
171  IF(typexm(1:1).EQ.'Q') THEN
172  DO i=1,9
173  CALL ov('X=CX ', x=xm(1,i), c=c, dim1=nelem)
174  ENDDO ! I
175  ELSEIF(typexm(1:1).NE.'0') THEN
176  WRITE(lu,11) typexm(1:1)
177  CALL plante(1)
178  stop
179  ENDIF
180 !
181 !-----------------------------------------------------------------------
182 !
183  ELSEIF(op(1:8).EQ.'M=M+CN ') THEN
184 !
185  IF(typdin(1:1).EQ.'I') THEN
186  CALL ov('X=X+C ', x=dm, c=c, dim1=ndiag)
187  ELSEIF(typdin(1:1).NE.'0') THEN
188  CALL ov('X=X+CY ', x=dm, y=dn, c=c, dim1=ndiag)
189  ENDIF
190 !
191  IF(typexn(1:1).EQ.'Q') THEN
192  IF(typexm(1:1).NE.'Q') THEN
193  WRITE(lu,98) typexm(1:1),op(1:8),typexn(1:1)
194 98 FORMAT(1x,'OM1112 (BIEF) : TYPEXM = ',a1,' DOES NOT GO ',
195  & /,1x,'FOR THE OPERATION : ',a8,' WITH TYPEXN = ',a1)
196  CALL plante(1)
197  stop
198  ENDIF
199  DO i=1,9
200  CALL ov('X=X+CY ', x=xm(1,i), y=xn(1,i), c=c, dim1=nelem)
201  ENDDO ! I
202  ELSEIF(typexn(1:1).NE.'0') THEN
203  WRITE(lu,11) typexn(1:1)
204  CALL plante(1)
205  stop
206  ENDIF
207 !
208 !-----------------------------------------------------------------------
209 !
210  ELSEIF(op(1:8).EQ.'M=M+N ') THEN
211 !
212  CALL ov('X=X+Y ', x=dm, y=dn, dim1=ndiag)
213 !
214  IF(typexn(1:1).EQ.'Q') THEN
215  IF(typexm(1:1).NE.'Q') THEN
216  WRITE(lu,98) typexm(1:1),op(1:8),typexn(1:1)
217  CALL plante(1)
218  stop
219  ENDIF
220  DO i=1,9
221  CALL ov('X=X+Y ', x=xm(1,i), y=xn(1,i), dim1=nelem)
222  ENDDO ! I
223  ELSEIF(typexn(1:1).NE.'0') THEN
224  WRITE(lu,11) typexn(1:1)
225  CALL plante(1)
226  stop
227  ENDIF
228 !
229 !-----------------------------------------------------------------------
230 !
231  ELSEIF(op(1:8).EQ.'M=MD ') THEN
232 !
233 ! DIAGONAL TERMS (DM AND D DON'T HAVE THE SAME TYPE HERE)
234 !
235  IF(typdim(1:1).EQ.'Q') THEN
236  CALL ov('X=XY ', x=dm, y=d, dim1=ndiag)
237  ELSEIF(typdim(1:1).EQ.'I') THEN
238  CALL ov('X=Y ', x=dm, y=d, dim1=ndiag)
239  typdim(1:1)='Q'
240  ELSEIF(typdim(1:1).NE.'0') THEN
241  WRITE(lu,13) typdim(1:1)
242 13 FORMAT(1x,'OM1112 (BIEF) : TYPDIM UNKNOWN :',a1)
243  CALL plante(1)
244  stop
245  ENDIF
246 !
247 ! EXTRADIAGONAL TERMS
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,4))
256  xm(ielem, 4) = xm(ielem, 4) * d(ikle(ielem,1))
257  xm(ielem, 5) = xm(ielem, 5) * d(ikle(ielem,3))
258  xm(ielem, 6) = xm(ielem, 6) * d(ikle(ielem,4))
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,4))
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=DM ') THEN
274 !
275 ! DIAGONAL TERMS
276 !
277  IF(typdim(1:1).EQ.'Q') THEN
278  CALL ov('X=XY ', x=dm, y=d, dim1=ndiag)
279  ELSEIF(typdim(1:1).EQ.'I') THEN
280  CALL ov('X=Y ', x=dm, y=d, dim1=ndiag)
281  typdim(1:1)='Q'
282  ELSEIF(typdim(1:1).NE.'0') THEN
283  WRITE(lu,13) typdim(1:1)
284  CALL plante(1)
285  stop
286  ENDIF
287 !
288 ! EXTRADIAGONAL TERMS
289 !
290  IF(typexm(1:1).EQ.'Q') THEN
291 !
292  DO ielem = 1 , nelem
293 !
294  xm(ielem, 1) = xm(ielem, 1) * d(ikle(ielem,1))
295  xm(ielem, 2) = xm(ielem, 2) * d(ikle(ielem,1))
296  xm(ielem, 3) = xm(ielem, 3) * d(ikle(ielem,1))
297  xm(ielem, 4) = xm(ielem, 4) * d(ikle(ielem,2))
298  xm(ielem, 5) = xm(ielem, 5) * d(ikle(ielem,2))
299  xm(ielem, 6) = xm(ielem, 6) * d(ikle(ielem,2))
300  xm(ielem, 7) = xm(ielem, 7) * d(ikle(ielem,3))
301  xm(ielem, 8) = xm(ielem, 8) * d(ikle(ielem,3))
302  xm(ielem, 9) = xm(ielem, 9) * d(ikle(ielem,3))
303 !
304  ENDDO
305 !
306  ELSEIF(typexm(1:1).NE.'0') THEN
307  WRITE(lu,163)
308 163 FORMAT(1x,'OM1112 (BIEF) : TYPEXM NOT AVAILABLE : ',a1)
309  CALL plante(1)
310  stop
311  ENDIF
312 !
313 !-----------------------------------------------------------------------
314 !
315  ELSEIF(op(1:8).EQ.'M=TN ') THEN
316 !
317 ! RECOPIES THE DIAGONAL
318 !
319  IF(typdin(1:1).EQ.'Q') THEN
320  CALL ov('X=Y ', x=dm, y=dn, dim1=ndiag)
321  ELSEIF(typdin(1:1).EQ.'I'.OR.typdin(1:1).EQ.'0') THEN
322 ! NOTHING TO DO, ONLY NEEDS TO COPY TYPDIN
323  ELSE
324  WRITE(lu,6) typdin(1:1)
325  CALL plante(1)
326  stop
327  ENDIF
328  typdim(1:1)=typdin(1:1)
329 !
330 ! TRANSPOSES EXTRADIAGONAL TERMS
331 !
332  IF(typexm(1:1).EQ.'Q') THEN
333 !
334  DO ielem = 1 , nelem
335 !
336  xm(ielem, 1) = xn(ielem, 3)
337  xm(ielem, 2) = xn(ielem, 5)
338  xm(ielem, 3) = xn(ielem, 7)
339  xm(ielem, 4) = xn(ielem, 1)
340  xm(ielem, 5) = xn(ielem, 6)
341  xm(ielem, 6) = xn(ielem, 8)
342  xm(ielem, 7) = xn(ielem, 2)
343  xm(ielem, 8) = xn(ielem, 4)
344  xm(ielem, 9) = xn(ielem, 9)
345 !
346  ENDDO
347 !
348  ELSEIF(typexm(1:1).NE.'0') THEN
349  WRITE(lu,163)
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,'OM1112 (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,'OM1112 (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
om1112
subroutine om1112(OP, DM, TYPDIM, XM, TYPEXM, DN, TYPDIN, XN, TYPEXN, D, C, IKLE, NELEM, NELMAX, NDIAG)
Definition: om1112.f:8
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:19 for The TELEMAC-MASCARET system by doxygen 1.8.13