The TELEMAC-MASCARET system  trunk
om1311.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE om1311
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 (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 !+ XM(IELEM, 10) ----> M(5,1)
49 !+ XM(IELEM, 11) ----> M(5,2)
50 !+ XM(IELEM, 12) ----> M(5,3)
51 !+ XM(IELEM, 13) ----> M(6,1)
52 !+ XM(IELEM, 14) ----> M(6,2)
53 !+ XM(IELEM, 15) ----> M(6,3)
54 !
55 !history ALGIANE FROEHLY (MATMECA); J-M HERVOUET (LNHE) ; F LEPEINTRE (LNH)
56 !+ 13/02/2008
57 !+ V5P9
58 !+
59 !
60 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
61 !+ 13/07/2010
62 !+ V6P0
63 !+ Translation of French comments within the FORTRAN sources into
64 !+ English comments
65 !
66 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
67 !+ 21/08/2010
68 !+ V6P0
69 !+ Creation of DOXYGEN tags for automated documentation and
70 !+ cross-referencing of the FORTRAN sources
71 !
72 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
73 !| C |-->| A GIVEN CONSTANT USED IN OPERATION OP
74 !| D |-->| A DIAGONAL MATRIX
75 !| DM |<->| DIAGONAL OF M
76 !| DN |-->| DIAGONAL OF N
77 !| IKLE |-->| CONNECTIVITY TABLE.
78 !| NDIAG |-->| NUMBER OF TERMS IN THE DIAGONAL
79 !| NELEM |-->| NUMBER OF ELEMENTS
80 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
81 !| OP |-->| OPERATION TO BE DONE (SEE ABOVE)
82 !| TYPDIM |<->| TYPE OF DIAGONAL OF M:
83 !| | | TYPDIM = 'Q' : ANY VALUE
84 !| | | TYPDIM = 'I' : IDENTITY
85 !| | | TYPDIM = '0' : ZERO
86 !| TYPDIN |<->| TYPE OF DIAGONAL OF N:
87 !| | | TYPDIN = 'Q' : ANY VALUE
88 !| | | TYPDIN = 'I' : IDENTITY
89 !| | | TYPDIN = '0' : ZERO
90 !| TYPEXM |-->| TYPE OF OFF-DIAGONAL TERMS OF M:
91 !| | | TYPEXM = 'Q' : ANY VALUE
92 !| | | TYPEXM = 'S' : SYMMETRIC
93 !| | | TYPEXM = '0' : ZERO
94 !| TYPEXN |-->| TYPE OF OFF-DIAGONAL TERMS OF N:
95 !| | | TYPEXN = 'Q' : ANY VALUE
96 !| | | TYPEXN = 'S' : SYMMETRIC
97 !| | | TYPEXN = '0' : ZERO
98 !| XM |-->| OFF-DIAGONAL TERMS OF M
99 !| XN |-->| OFF-DIAGONAL TERMS OF N
100 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
101 !
102  USE bief !, EX_OM1311 => OM1311
103 !
104  USE declarations_special
105  IMPLICIT NONE
106 !
107 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
108 !
109  INTEGER, INTENT(IN) :: NELEM,NELMAX,NDIAG
110  INTEGER, INTENT(IN) :: IKLE(nelmax,6)
111  CHARACTER(LEN=8), INTENT(IN) :: OP
112  DOUBLE PRECISION, INTENT(IN) :: DN(*),D(*),XN(nelmax,*)
113  DOUBLE PRECISION, INTENT(INOUT) :: DM(*),XM(nelmax,*)
114  CHARACTER(LEN=1), INTENT(INOUT) :: TYPDIM,TYPEXM,TYPDIN,TYPEXN
115  DOUBLE PRECISION, INTENT(IN) :: C
116 !
117 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
118 !
119  INTEGER I,J,IELEM
120 !
121 !-----------------------------------------------------------------------
122 !
123  IF(op(1:8).EQ.'M=N ') THEN
124 !
125  IF(typdin(1:1).EQ.'Q') THEN
126  CALL ov('X=Y ', x=dm, y=dn, dim1=ndiag)
127  ELSEIF(typdin(1:1).EQ.'I'.OR.typdin(1:1).EQ.'0') THEN
128 ! NOTHING TO DO, ONLY NEEDS TO COPY TYPDIN
129  ELSE
130  WRITE(lu,6) typdin(1:1)
131 6 FORMAT(1x,'OM1311 (BIEF) : TYPDIN UNKNOWN :',a1)
132  CALL plante(1)
133  stop
134  ENDIF
135  typdim(1:1)=typdin(1:1)
136 !
137  IF(typexn(1:1).EQ.'Q') THEN
138  DO i=1,15
139  CALL ov('X=Y ', x=xm(1,i), y=xn(1,i), dim1=nelem)
140  ENDDO ! I
141  ELSEIF(typexn(1:1).NE.'0') THEN
142  WRITE(lu,11) typexn(1:1)
143 11 FORMAT(1x,'OM1311 (BIEF) : TYPEXN UNKNOWN :',a1)
144  CALL plante(1)
145  stop
146  ENDIF
147 !
148  typexm(1:1)=typexn(1:1)
149 !
150 !-----------------------------------------------------------------------
151 !
152  ELSEIF(op(1:8).EQ.'M=CN ') THEN
153 !
154  CALL ov('X=CY ', x=dm, y=dn, c=c, dim1=ndiag)
155 !
156  IF(typexn(1:1).EQ.'Q') THEN
157  DO i=1,15
158  CALL ov('X=CY ', x=xm(1,i), y=xn(1,i), c=c, dim1=nelem)
159  ENDDO ! I
160  ELSEIF(typexn(1:1).NE.'0') THEN
161  WRITE(lu,11) typexn(1:1)
162  CALL plante(1)
163  stop
164  ENDIF
165 !
166  typdim(1:1)=typdin(1:1)
167  typexm(1:1)=typexn(1:1)
168 !
169 !-----------------------------------------------------------------------
170 !
171  ELSEIF(op(1:8).EQ.'M=CM ') THEN
172 !
173  CALL ov('X=CX ', x=dm, c=c, dim1=ndiag)
174 !
175  IF(typexm(1:1).EQ.'Q') THEN
176  DO i=1,15
177  CALL ov('X=CX ', x=xm(1,i), c=c, dim1=nelem)
178  ENDDO ! I
179  ELSEIF(typexm(1:1).NE.'0') THEN
180  WRITE(lu,11) typexm(1:1)
181  CALL plante(1)
182  stop
183  ENDIF
184 !
185 !-----------------------------------------------------------------------
186 !
187  ELSEIF(op(1:8).EQ.'M=M+CN ') THEN
188 !
189  IF(typdin(1:1).EQ.'I') THEN
190  CALL ov('X=X+C ', x=dm, c=c, dim1=ndiag)
191  ELSEIF(typdin(1:1).NE.'0') THEN
192  CALL ov('X=X+CY ', x=dm, y=dn, c=c, dim1=ndiag)
193  ENDIF
194 !
195  IF(typexn(1:1).EQ.'Q') THEN
196  IF(typexm(1:1).NE.'Q') THEN
197  WRITE(lu,98) typexm(1:1),op(1:8),typexn(1:1)
198 98 FORMAT(1x,'OM1311 (BIEF) : TYPEXM = ',a1,' DOES NOT GO ',
199  & /,1x,'FOR THE OPERATION : ',a8,' WITH TYPEXN = ',a1)
200  CALL plante(1)
201  stop
202  ENDIF
203  DO i=1,15
204  CALL ov('X=X+CY ', x=xm(1,i), y=xn(1,i), c=c, dim1=nelem)
205  ENDDO ! I
206  ELSEIF(typexn(1:1).NE.'0') THEN
207  WRITE(lu,11) typexn(1:1)
208  CALL plante(1)
209  stop
210  ENDIF
211 !
212 !-----------------------------------------------------------------------
213 !
214  ELSEIF(op(1:8).EQ.'M=M+N ') THEN
215 !
216  CALL ov('X=X+Y ', x=dm, y=dn, dim1=ndiag)
217 !
218  IF(typexn(1:1).EQ.'Q') THEN
219  IF(typexm(1:1).NE.'Q') THEN
220  WRITE(lu,98) typexm(1:1),op(1:8),typexn(1:1)
221  CALL plante(1)
222  stop
223  ENDIF
224  DO i=1,15
225  CALL ov('X=X+Y ', x=xm(1,i), y=xn(1,i), dim1=nelem)
226  ENDDO
227  ELSEIF(typexn(1:1).NE.'0') THEN
228  WRITE(lu,11) typexn(1:1)
229  CALL plante(1)
230  stop
231  ENDIF
232 !
233 !-----------------------------------------------------------------------
234 !
235  ELSEIF(op(1:8).EQ.'M=MD ') THEN
236 !
237 ! DIAGONAL TERMS (DM AND D CAN HAVE DIFFERENT TYPES)
238 !
239 !
240  IF(typdim(1:1).EQ.'Q') THEN
241  CALL ov('X=XY ', x=dm, y=d, dim1=ndiag)
242  ELSEIF(typdim(1:1).EQ.'I') THEN
243  CALL ov('X=Y ', x=dm, y=d, dim1=ndiag)
244  typdim(1:1)='Q'
245  ELSEIF(typdim(1:1).NE.'0') THEN
246  WRITE(lu,13) typdim(1:1)
247 13 FORMAT(1x,'OM1311 (BIEF) : TYPDIM UNKNOWN :',a1)
248  CALL plante(1)
249  stop
250  ENDIF
251 !
252 ! EXTRADIAGONAL TERMS
253 !
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,1))
262  xm(ielem, 4) = xm(ielem, 4) * d(ikle(ielem,3))
263  xm(ielem, 5) = xm(ielem, 5) * d(ikle(ielem,1))
264  xm(ielem, 6) = xm(ielem, 6) * d(ikle(ielem,2))
265  xm(ielem, 7) = xm(ielem, 7) * d(ikle(ielem,1))
266  xm(ielem, 8) = xm(ielem, 8) * d(ikle(ielem,2))
267  xm(ielem, 9) = xm(ielem, 9) * d(ikle(ielem,3))
268  xm(ielem, 10) = xm(ielem, 10) * d(ikle(ielem,1))
269  xm(ielem, 11) = xm(ielem, 11) * d(ikle(ielem,2))
270  xm(ielem, 12) = xm(ielem, 12) * d(ikle(ielem,3))
271  xm(ielem, 13) = xm(ielem, 13) * d(ikle(ielem,1))
272  xm(ielem, 14) = xm(ielem, 14) * d(ikle(ielem,2))
273  xm(ielem, 15) = xm(ielem, 15) * d(ikle(ielem,3))
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=TN ') THEN
286 !
287 ! RECOPIES EXTRADIAGONAL TERMS
288 !
289  IF(typdin(1:1).EQ.'Q') THEN
290  CALL ov('X=Y ', x=dm, y=dn, dim1=ndiag)
291  ELSEIF(typdin(1:1).EQ.'I'.OR.typdin(1:1).EQ.'0') THEN
292 ! NOTHING TO DO, ONLY NEEDS TO COPY TYPDIN
293  ELSE
294  WRITE(lu,6) typdin(1:1)
295  CALL plante(1)
296  stop
297  ENDIF
298  typdim(1:1)=typdin(1:1)
299 !
300 ! TRANSPOSES EXTRADIAGONAL TERMS
301 !
302  IF(typexm(1:1).EQ.'Q') THEN
303 !
304  DO ielem = 1 , nelem
305 !
306  xm(ielem, 1) = xn(ielem, 6)
307  xm(ielem, 2) = xn(ielem, 11)
308  xm(ielem, 3) = xn(ielem, 1)
309  xm(ielem, 4) = xn(ielem, 12)
310  xm(ielem, 5) = xn(ielem, 2)
311  xm(ielem, 6) = xn(ielem, 7)
312  xm(ielem, 7) = xn(ielem, 3)
313  xm(ielem, 8) = xn(ielem, 8)
314  xm(ielem, 9) = xn(ielem, 13)
315  xm(ielem, 10) = xn(ielem, 4)
316  xm(ielem, 11) = xn(ielem, 9)
317  xm(ielem, 12) = xn(ielem, 14)
318  xm(ielem, 13) = xn(ielem, 5)
319  xm(ielem, 14) = xn(ielem, 10)
320  xm(ielem, 15) = xn(ielem, 15)
321 !
322  ENDDO
323 !
324  ELSEIF(typexm(1:1).NE.'0') THEN
325  WRITE(lu,163)
326  CALL plante(1)
327  stop
328  ENDIF
329 !
330 !-----------------------------------------------------------------------
331 !
332  ELSEIF(op(1:8).EQ.'M=DM ') THEN
333 !
334 ! DIAGONAL TERMS
335 !
336  IF(typdim(1:1).EQ.'Q') THEN
337  CALL ov('X=XY ', x=dm, y=d, dim1=ndiag)
338  ELSEIF(typdim(1:1).EQ.'I') THEN
339  CALL ov('X=Y ', x=dm, y=d, dim1=ndiag)
340  typdim(1:1)='Q'
341  ELSEIF(typdim(1:1).NE.'0') THEN
342  WRITE(lu,13) typdim(1:1)
343  CALL plante(1)
344  stop
345  ENDIF
346 !
347 ! EXTRADIAGONAL TERMS
348 !
349  IF(typexm(1:1).EQ.'Q') THEN
350 !
351  DO ielem = 1 , nelem
352 !
353  xm(ielem, 1) = xm(ielem, 1) * d(ikle(ielem,1))
354  xm(ielem, 2) = xm(ielem, 2) * d(ikle(ielem,1))
355  xm(ielem, 3) = xm(ielem, 3) * d(ikle(ielem,2))
356  xm(ielem, 4) = xm(ielem, 4) * d(ikle(ielem,2))
357  xm(ielem, 5) = xm(ielem, 5) * d(ikle(ielem,3))
358  xm(ielem, 6) = xm(ielem, 6) * d(ikle(ielem,3))
359  xm(ielem, 7) = xm(ielem, 7) * d(ikle(ielem,4))
360  xm(ielem, 8) = xm(ielem, 8) * d(ikle(ielem,4))
361  xm(ielem, 9) = xm(ielem, 9) * d(ikle(ielem,4))
362  xm(ielem, 10) = xm(ielem, 10) * d(ikle(ielem,5))
363  xm(ielem, 11) = xm(ielem, 11) * d(ikle(ielem,5))
364  xm(ielem, 12) = xm(ielem, 12) * d(ikle(ielem,5))
365  xm(ielem, 13) = xm(ielem, 13) * d(ikle(ielem,6))
366  xm(ielem, 14) = xm(ielem, 14) * d(ikle(ielem,6))
367  xm(ielem, 15) = xm(ielem, 15) * d(ikle(ielem,6))
368 !
369  ENDDO
370 !
371  ELSEIF(typexm(1:1).NE.'0') THEN
372  WRITE(lu,163)
373 163 FORMAT(1x,'OM1311 (BIEF) : TYPEXM NOT AVAILABLE : ',a1)
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,'OM11311 (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,'OM1311 (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
om1311
subroutine om1311(OP, DM, TYPDIM, XM, TYPEXM, DN, TYPDIN, XN, TYPEXN, D, C, IKLE, NELEM, NELMAX, NDIAG)
Definition: om1311.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