The TELEMAC-MASCARET system  trunk
om1113.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE om1113
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 6 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(1,5)
45 !+ XM(IELEM, 5) ----> M(1,6)
46 !+ XM(IELEM, 6) ----> M(2,1)
47 !+ XM(IELEM, 7) ----> M(2,3)
48 !+ XM(IELEM, 8) ----> M(2,4)
49 !+ XM(IELEM, 9) ----> M(2,5)
50 !+ XM(IELEM, 10) ----> M(2,6)
51 !+ XM(IELEM, 11) ----> M(3,1)
52 !+ XM(IELEM, 12) ----> M(3,2)
53 !+ XM(IELEM, 13) ----> M(3,4)
54 !+ XM(IELEM, 14) ----> M(3,5)
55 !+ XM(IELEM, 15) ----> M(3,6)
56 !
57 !history ALGIANE FROEHLY (MATMECA); J-M HERVOUET (LNHE) ; F LEPEINTRE (LNH)
58 !+ 13/02/08
59 !+ V5P9
60 !+
61 !
62 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
63 !+ 13/07/2010
64 !+ V6P0
65 !+ Translation of French comments within the FORTRAN sources into
66 !+ English comments
67 !
68 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
69 !+ 21/08/2010
70 !+ V6P0
71 !+ Creation of DOXYGEN tags for automated documentation and
72 !+ cross-referencing of the FORTRAN sources
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_OM1113 => OM1113
105 !
106  USE declarations_special
107  IMPLICIT NONE
108 !
109 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
110 !
111  INTEGER, INTENT(IN) :: NELEM,NELMAX,NDIAG
112  INTEGER, INTENT(IN) :: IKLE(nelmax,6)
113  CHARACTER(LEN=8), INTENT(IN) :: OP
114  DOUBLE PRECISION, INTENT(IN) :: DN(*),D(*),XN(nelmax,*)
115  DOUBLE PRECISION, INTENT(INOUT) :: DM(*),XM(nelmax,*)
116  CHARACTER(LEN=1), INTENT(INOUT) :: TYPDIM,TYPEXM,TYPDIN,TYPEXN
117  DOUBLE PRECISION, INTENT(IN) :: C
118 !
119 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
120 !
121  INTEGER I,J,IELEM
122 !
123 !-----------------------------------------------------------------------
124 !
125  IF(op(1:8).EQ.'M=N ') THEN
126 !
127  IF(typdin(1:1).EQ.'Q') THEN
128  CALL ov('X=Y ', x=dm, y=dn, dim1=ndiag)
129  ELSEIF(typdin(1:1).EQ.'I'.OR.typdin(1:1).EQ.'0') THEN
130 ! NOTHING TO DO, ONLY NEEDS TO COPY TYPDIN
131  ELSE
132  WRITE(lu,6) typdin(1:1)
133 6 FORMAT(1x,'OM1113 (BIEF) : TYPDIN UNKNOWN :',a1)
134  CALL plante(1)
135  stop
136  ENDIF
137  typdim(1:1)=typdin(1:1)
138 !
139  IF(typexn(1:1).EQ.'Q') THEN
140  DO i=1,15
141  CALL ov('X=Y ', x=xm(1,i), y=xn(1,i), dim1=nelem)
142  ENDDO ! I
143  ELSEIF(typexn(1:1).NE.'0') THEN
144  WRITE(lu,11) typexn(1:1)
145 11 FORMAT(1x,'OM1113 (BIEF) : TYPEXN UNKNOWN :',a1)
146  CALL plante(1)
147  stop
148  ENDIF
149 !
150  typexm(1:1)=typexn(1:1)
151 !
152 !-----------------------------------------------------------------------
153 !
154  ELSEIF(op(1:8).EQ.'M=CN ') THEN
155 !
156  CALL ov('X=CY ', x=dm, y=dn, c=c, dim1=ndiag)
157 !
158  IF(typexn(1:1).EQ.'Q') THEN
159  DO i=1,15
160  CALL ov('X=CY ', x=xm(1,i), y=xn(1,i), c=c, dim1=nelem)
161  ENDDO ! I
162  ELSEIF(typexn(1:1).NE.'0') THEN
163  WRITE(lu,11) typexn(1:1)
164  CALL plante(1)
165  stop
166  ENDIF
167 !
168  typdim(1:1)=typdin(1:1)
169  typexm(1:1)=typexn(1:1)
170 !
171 !-----------------------------------------------------------------------
172 !
173  ELSEIF(op(1:8).EQ.'M=CM ') THEN
174 !
175  CALL ov('X=CX ', x=dm, c=c, dim1=ndiag)
176 !
177  IF(typexm(1:1).EQ.'Q') THEN
178  DO i=1,15
179  CALL ov('X=CX ', x=xm(1,i), c=c, dim1=nelem)
180  ENDDO ! I
181  ELSEIF(typexm(1:1).NE.'0') THEN
182  WRITE(lu,11) typexm(1:1)
183  CALL plante(1)
184  stop
185  ENDIF
186 !
187 !-----------------------------------------------------------------------
188 !
189  ELSEIF(op(1:8).EQ.'M=M+CN ') THEN
190 !
191  IF(typdin(1:1).EQ.'I') THEN
192  CALL ov('X=X+C ', x=dm, c=c, dim1=ndiag)
193  ELSEIF(typdin(1:1).NE.'0') THEN
194  CALL ov('X=X+CY ', x=dm, y=dn, c=c, dim1=ndiag)
195  ENDIF
196 !
197  IF(typexn(1:1).EQ.'Q') THEN
198  IF(typexm(1:1).NE.'Q') THEN
199  WRITE(lu,98) typexm(1:1),op(1:8),typexn(1:1)
200 98 FORMAT(1x,'OM1113 (BIEF) : TYPEXM = ',a1,' DOES NOT GO ',
201  & /,1x,'FOR THE OPERATION : ',a8,' WITH TYPEXN = ',a1)
202  CALL plante(1)
203  stop
204  ENDIF
205  DO i=1,15
206  CALL ov('X=X+CY ', x=xm(1,i), y=xn(1,i), c=c, dim1=nelem)
207  ENDDO ! I
208  ELSEIF(typexn(1:1).NE.'0') THEN
209  WRITE(lu,11) typexn(1:1)
210  CALL plante(1)
211  stop
212  ENDIF
213 !
214 !-----------------------------------------------------------------------
215 !
216  ELSEIF(op(1:8).EQ.'M=M+N ') THEN
217 !
218  CALL ov('X=X+Y ', x=dm, y=dn, dim1=ndiag)
219 !
220  IF(typexn(1:1).EQ.'Q') THEN
221  IF(typexm(1:1).NE.'Q') THEN
222  WRITE(lu,98) typexm(1:1),op(1:8),typexn(1:1)
223  CALL plante(1)
224  stop
225  ENDIF
226  DO i=1,15
227  CALL ov('X=X+Y ', x=xm(1,i), y=xn(1,i), dim1=nelem)
228  ENDDO ! I
229  ELSEIF(typexn(1:1).NE.'0') THEN
230  WRITE(lu,11) typexn(1:1)
231  CALL plante(1)
232  stop
233  ENDIF
234 !
235 !-----------------------------------------------------------------------
236 !
237  ELSEIF(op(1:8).EQ.'M=MD ') THEN
238 !
239 ! DIAGONAL TERMS (DM AND D DON'T HAVE THE SAME TYPE HERE)
240 !
241  IF(typdim(1:1).EQ.'Q') THEN
242  CALL ov('X=XY ', x=dm, y=d, dim1=ndiag)
243  ELSEIF(typdim(1:1).EQ.'I') THEN
244  CALL ov('X=Y ', x=dm, y=d, dim1=ndiag)
245  typdim(1:1)='Q'
246  ELSEIF(typdim(1:1).NE.'0') THEN
247  WRITE(lu,13) typdim(1:1)
248 13 FORMAT(1x,'OM1113 (BIEF) : TYPDIM UNKNOWN :',a1)
249  CALL plante(1)
250  stop
251  ENDIF
252 !
253 ! EXTRADIAGONAL TERMS
254 !
255  IF(typexm(1:1).EQ.'Q') THEN
256 !
257  DO ielem = 1 , nelem
258 !
259  xm(ielem, 1) = xm(ielem, 1) * d(ikle(ielem,2))
260  xm(ielem, 2) = xm(ielem, 2) * d(ikle(ielem,3))
261  xm(ielem, 3) = xm(ielem, 3) * d(ikle(ielem,4))
262  xm(ielem, 4) = xm(ielem, 4) * d(ikle(ielem,5))
263  xm(ielem, 5) = xm(ielem, 5) * d(ikle(ielem,6))
264  xm(ielem, 6) = xm(ielem, 6) * d(ikle(ielem,1))
265  xm(ielem, 7) = xm(ielem, 7) * d(ikle(ielem,3))
266  xm(ielem, 8) = xm(ielem, 8) * d(ikle(ielem,4))
267  xm(ielem, 9) = xm(ielem, 9) * d(ikle(ielem,5))
268  xm(ielem, 10) = xm(ielem, 10) * d(ikle(ielem,6))
269  xm(ielem, 11) = xm(ielem, 11) * d(ikle(ielem,1))
270  xm(ielem, 12) = xm(ielem, 12) * d(ikle(ielem,2))
271  xm(ielem, 13) = xm(ielem, 13) * d(ikle(ielem,4))
272  xm(ielem, 14) = xm(ielem, 14) * d(ikle(ielem,5))
273  xm(ielem, 15) = xm(ielem, 15) * d(ikle(ielem,6))
274 !
275  ENDDO
276 !
277  ELSEIF(typexm(1:1).NE.'0') THEN
278  WRITE(lu,163)
279  CALL plante(1)
280  stop
281  ENDIF
282 !
283 !-----------------------------------------------------------------------
284 !
285  ELSEIF(op(1:8).EQ.'M=DM ') THEN
286 !
287 ! DIAGONAL TERMS
288 !
289  IF(typdim(1:1).EQ.'Q') THEN
290  CALL ov('X=XY ', x=dm, y=d, dim1=ndiag)
291  ELSEIF(typdim(1:1).EQ.'I') THEN
292  CALL ov('X=Y ', x=dm, y=d, dim1=ndiag)
293  typdim(1:1)='Q'
294  ELSEIF(typdim(1:1).NE.'0') THEN
295  WRITE(lu,13) typdim(1:1)
296  CALL plante(1)
297  stop
298  ENDIF
299 !
300 ! EXTRADIAGONAL TERMS
301 !
302  IF(typexm(1:1).EQ.'Q') THEN
303 !
304  DO ielem = 1 , nelem
305 !
306  xm(ielem, 1) = xm(ielem, 1) * d(ikle(ielem,1))
307  xm(ielem, 2) = xm(ielem, 2) * d(ikle(ielem,1))
308  xm(ielem, 3) = xm(ielem, 3) * d(ikle(ielem,1))
309  xm(ielem, 4) = xm(ielem, 4) * d(ikle(ielem,1))
310  xm(ielem, 5) = xm(ielem, 5) * d(ikle(ielem,1))
311  xm(ielem, 6) = xm(ielem, 6) * d(ikle(ielem,2))
312  xm(ielem, 7) = xm(ielem, 7) * d(ikle(ielem,2))
313  xm(ielem, 8) = xm(ielem, 8) * d(ikle(ielem,2))
314  xm(ielem, 9) = xm(ielem, 9) * d(ikle(ielem,2))
315  xm(ielem, 10) = xm(ielem, 10) * d(ikle(ielem,2))
316  xm(ielem, 11) = xm(ielem, 11) * d(ikle(ielem,3))
317  xm(ielem, 12) = xm(ielem, 12) * d(ikle(ielem,3))
318  xm(ielem, 13) = xm(ielem, 13) * d(ikle(ielem,3))
319  xm(ielem, 14) = xm(ielem, 14) * d(ikle(ielem,3))
320  xm(ielem, 15) = xm(ielem, 15) * d(ikle(ielem,3))
321 !
322  ENDDO
323 !
324  ELSEIF(typexm(1:1).NE.'0') THEN
325  WRITE(lu,163)
326 163 FORMAT(1x,'OM1113 (BIEF) : TYPEXM NOT AVAILABLE : ',a1)
327  CALL plante(1)
328  stop
329  ENDIF
330 !
331 !-----------------------------------------------------------------------
332 !
333  ELSEIF(op(1:8).EQ.'M=TN ') THEN
334 !
335 ! RECOPIES THE DIAGONAL
336 !
337  IF(typdin(1:1).EQ.'Q') THEN
338  CALL ov('X=Y ', x=dm, y=dn, dim1=ndiag)
339  ELSEIF(typdin(1:1).EQ.'I'.OR.typdin(1:1).EQ.'0') THEN
340 ! NOTHING TO DO, ONLY NEEDS TO COPY TYPDIN
341  ELSE
342  WRITE(lu,6) typdin(1:1)
343  CALL plante(1)
344  stop
345  ENDIF
346  typdim(1:1)=typdin(1:1)
347 !
348 ! TRANSPOSES EXTRADIAGONAL TERMS
349 !
350  IF(typexm(1:1).EQ.'Q') THEN
351 !
352  DO ielem = 1 , nelem
353 !
354  xm(ielem, 1) = xn(ielem, 3)
355  xm(ielem, 2) = xn(ielem, 5)
356  xm(ielem, 3) = xn(ielem, 7)
357  xm(ielem, 4) = xn(ielem, 10)
358  xm(ielem, 5) = xn(ielem, 13)
359  xm(ielem, 6) = xn(ielem, 1)
360  xm(ielem, 7) = xn(ielem, 6)
361  xm(ielem, 8) = xn(ielem, 8)
362  xm(ielem, 9) = xn(ielem, 11)
363  xm(ielem, 10) = xn(ielem, 14)
364  xm(ielem, 11) = xn(ielem, 2)
365  xm(ielem, 12) = xn(ielem, 4)
366  xm(ielem, 13) = xn(ielem, 9)
367  xm(ielem, 14) = xn(ielem, 12)
368  xm(ielem, 15) = xn(ielem, 15)
369 !
370  ENDDO
371 !
372  ELSEIF(typexm(1:1).NE.'0') THEN
373  WRITE(lu,163)
374  CALL plante(1)
375  stop
376  ENDIF
377 !
378 !-----------------------------------------------------------------------
379 !
380  ELSEIF(op(1:8).EQ.'M=M+D ') THEN
381 !
382  IF(typdim(1:1).EQ.'Q') THEN
383  CALL ov('X=X+Y ', x=dm, y=d, dim1=ndiag)
384  ELSE
385  WRITE(lu,13) typdim(1:1)
386  CALL plante(1)
387  stop
388  ENDIF
389 !
390 !-----------------------------------------------------------------------
391 !
392  ELSEIF(op(1:8).EQ.'M=0 ') THEN
393 !
394  CALL ov('X=C ', x=dm, c=0.d0, dim1=ndiag)
395 !
396  IF(typexm(1:1).EQ.'Q') THEN
397  DO i=1,15
398  CALL ov('X=C ', x=xm(1,i), c=0.d0, dim1=nelem)
399  ENDDO
400  ELSEIF(typexm(1:1).NE.'0') THEN
401  WRITE(lu,711) typexm(1:1)
402 711 FORMAT(1x,'OM1113 (BIEF) : TYPEXM UNKNOWN :',a1)
403  CALL plante(1)
404  stop
405  ENDIF
406 !
407  typdim(1:1)='0'
408 ! TYPEXM(1:1) IS NOT CHANGED
409 !
410 !-----------------------------------------------------------------------
411 !
412  ELSEIF(op(1:8).EQ.'M=MSK(M)') THEN
413 !
414  IF(typexm(1:1).EQ.'Q') THEN
415  j = 15
416  ELSEIF(typexm(1:1).EQ.'0') THEN
417  j = 0
418  ELSE
419  WRITE(lu,711) typexm
420  j=0
421  CALL plante(1)
422  stop
423  ENDIF
424 !
425  IF(j.GT.0) THEN
426  DO i = 1,j
427  CALL ov('X=XY ', x=xm(1,i), y=d, dim1=nelem)
428  ENDDO
429  ENDIF
430 !
431 !-----------------------------------------------------------------------
432 !
433  ELSE
434 !
435  WRITE(lu,41) op
436 41 FORMAT(1x,'OM1113 (BIEF) : UNKNOWN OPERATION : ',a8)
437  CALL plante(1)
438  stop
439 !
440  ENDIF
441 !
442 !-----------------------------------------------------------------------
443 !
444  RETURN
445  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
om1113
subroutine om1113(OP, DM, TYPDIM, XM, TYPEXM, DN, TYPDIN, XN, TYPEXN, D, C, IKLE, NELEM, NELMAX, NDIAG)
Definition: om1113.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