The TELEMAC-MASCARET system  trunk
om3181.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE om3181
3 ! *****************
4 !
5  &(op , dm,typdim,xm,typexm, dn,typdin,xn,typexn, c,
6  & nulone,nelbor,nbor,nelmax,sizdn,neleb)
7 !
8 !***********************************************************************
9 ! BIEF V6P1 21/08/2010
10 !***********************************************************************
11 !
12 !brief OPERATIONS ON MATRICES.
13 !code
14 !+ M: T1 TETRAHEDRON
15 !+ N: BOUNDARY MATRIX, T1 TRIANGLE
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=M+N ' : ADDS N TO M
25 !+ OP = 'M=M+TN ' : ADDS TRANSPOSE(N) TO M
26 !
27 !code
28 !+ CONVENTION FOR THE STORAGE OF EXTRA-DIAGONAL TERMS:
29 !+
30 !+ XM( ,1) ----> M(1,2)
31 !+ XM( ,2) ----> M(1,3)
32 !+ XM( ,3) ----> M(2,3)
33 !+ XM( ,4) ----> M(2,1)
34 !+ XM( ,5) ----> M(3,1)
35 !+ XM( ,6) ----> M(3,2)
36 !
37 !history J-M HERVOUET (LNHE)
38 !+ 23/06/2008
39 !+ V5P9
40 !+
41 !
42 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
43 !+ 13/07/2010
44 !+ V6P0
45 !+ Translation of French comments within the FORTRAN sources into
46 !+ English comments
47 !
48 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
49 !+ 21/08/2010
50 !+ V6P0
51 !+ Creation of DOXYGEN tags for automated documentation and
52 !+ cross-referencing of the FORTRAN sources
53 !
54 !history S.E.BOURBAN (HRW)
55 !+ 21/03/2017
56 !+ V7P3
57 !+ Replacement of the DATA declarations by the PARAMETER associates
58 !
59 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
60 !| C |-->| A GIVEN CONSTANT USED IN OPERATION OP
61 !| DM |<->| DIAGONAL OF M
62 !| DN |-->| DIAGONAL OF N
63 !| NBOR |-->| GLOBAL NUMBER OF BOUNDARY POINTS
64 !| NELBOR |-->| FOR THE KTH BOUNDARY EDGE, GIVES THE CORRESPONDING
65 !| | | ELEMENT.
66 !| NELEB |-->| NUMBER OF BOUNDARY ELEMENTS
67 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
68 !| NULONE |-->| GOES WITH ARRAY NELBOR. NELBOR GIVES THE
69 !| | | ADJACENT ELEMENT, NULONE GIVES THE LOCAL
70 !| | | NUMBER OF THE FIRST NODE OF THE BOUNDARY EDGE
71 !| | | I.E. 1, 2 OR 3 FOR TRIANGLES.
72 !| OP |-->| OPERATION TO BE DONE (SEE ABOVE)
73 !| SIZDN |-->| SIZE OF DIAGONAL DN
74 !| TYPDIM |<->| TYPE OF DIAGONAL OF M:
75 !| | | TYPDIM = 'Q' : ANY VALUE
76 !| | | TYPDIM = 'I' : IDENTITY
77 !| | | TYPDIM = '0' : ZERO
78 !| TYPDIN |<->| TYPE OF DIAGONAL OF N:
79 !| | | TYPDIN = 'Q' : ANY VALUE
80 !| | | TYPDIN = 'I' : IDENTITY
81 !| | | TYPDIN = '0' : ZERO
82 !| TYPEXM |-->| TYPE OF OFF-DIAGONAL TERMS OF M:
83 !| | | TYPEXM = 'Q' : ANY VALUE
84 !| | | TYPEXM = 'S' : SYMMETRIC
85 !| | | TYPEXM = '0' : ZERO
86 !| TYPEXN |-->| TYPE OF OFF-DIAGONAL TERMS OF N:
87 !| | | TYPEXN = 'Q' : ANY VALUE
88 !| | | TYPEXN = 'S' : SYMMETRIC
89 !| | | TYPEXN = '0' : ZERO
90 !| XM |-->| OFF-DIAGONAL TERMS OF M
91 !| XN |-->| OFF-DIAGONAL TERMS OF N
92 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
93 !
94  USE bief
95 !
96  USE declarations_special
97  IMPLICIT NONE
98 !
99 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
100 !
101  INTEGER, INTENT(IN) :: NELMAX,SIZDN,NELEB
102  CHARACTER(LEN=8), INTENT(IN) :: OP
103  INTEGER, INTENT(IN) :: NULONE(3*neleb)
104  INTEGER, INTENT(IN) :: NELBOR(*),NBOR(*)
105  DOUBLE PRECISION, INTENT(IN) :: DN(*),XN(neleb,*)
106  DOUBLE PRECISION, INTENT(INOUT) :: DM(*),XM(nelmax,*)
107  CHARACTER(LEN=1), INTENT(INOUT) :: TYPDIM,TYPEXM,TYPDIN,TYPEXN
108  DOUBLE PRECISION, INTENT(IN) :: C
109 !
110 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
111 !
112  INTEGER K,IEL,NUL1,NUL2,NUL3
113 !
114  DOUBLE PRECISION Z(1)
115 !
116 !-----------------------------------------------------------------------
117 !
118 ! BEWARE: IN FORTRAN THE FIRST INDEX VARIES MOST QUICKLY
119 ! 123456789 SHOULD NOT BE USED
120 !
121  INTEGER :: CONVNSY(4,4)
122  parameter( convnsy = reshape( (/
123  & 123456789, 7 , 8 , 9 ,
124  & 1 , 123456789, 10 , 11 ,
125  & 2 , 4 , 123456789, 12 ,
126  & 3 , 5 , 6 , 123456789 /),
127  & shape=(/ 4,4 /) ) )
128  INTEGER :: CONVSYM(4,4)
129  parameter( convsym = reshape( (/
130  & 123456789, 1 , 2 , 3 ,
131  & 1 , 123456789, 4 , 5 ,
132  & 2 , 4 , 123456789, 6 ,
133  & 3 , 5 , 6 , 123456789 /),
134  & shape=(/ 4,4 /) ) )
135 !
136 !***********************************************************************
137 !
138  IF(op(1:8).EQ.'M=M+N ') THEN
139  IF(typdim.EQ.'Q'.AND.typdin.EQ.'Q') THEN
140  CALL ovdb( 'X=X+Y ' , dm , dn , z , c , nbor , sizdn )
141  ELSE
142  WRITE(lu,199) typdim(1:1),op(1:8),typdin(1:1)
143 199 FORMAT(1x,'OM5161 (BIEF) : TYPDIM = ',a1,' NOT IMPLEMENTED',
144  & /,1x,'FOR THE OPERATION : ',a8,' WITH TYPDIN = ',a1)
145  CALL plante(1)
146  stop
147  ENDIF
148 !
149  IF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'Q') THEN
150 !
151 ! CASE WHERE BOTH MATRICES ARE NONSYMMETRICAL
152 !
153  IF(ncsize.GT.1) THEN
154  DO k = 1 , neleb
155  iel = nelbor(k)
156  IF(iel.GT.0) THEN
157  nul1=nulone(k)
158  nul2=nulone(k+neleb)
159  nul3=nulone(k+2*neleb)
160  xm( iel , convnsy(nul1,nul2) ) =
161  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 1)
162  xm( iel , convnsy(nul1,nul3) ) =
163  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 2)
164  xm( iel , convnsy(nul2,nul3) ) =
165  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 3)
166  xm( iel , convnsy(nul2,nul1) ) =
167  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 4)
168  xm( iel , convnsy(nul3,nul1) ) =
169  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 5)
170  xm( iel , convnsy(nul3,nul2) ) =
171  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 6)
172  ENDIF
173  ENDDO
174  ELSE
175  DO k = 1 , neleb
176  iel = nelbor(k)
177  nul1=nulone(k)
178  nul2=nulone(k+neleb)
179  nul3=nulone(k+2*neleb)
180  xm( iel , convnsy(nul1,nul2) ) =
181  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 1)
182  xm( iel , convnsy(nul1,nul3) ) =
183  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 2)
184  xm( iel , convnsy(nul2,nul3) ) =
185  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 3)
186  xm( iel , convnsy(nul2,nul1) ) =
187  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 4)
188  xm( iel , convnsy(nul3,nul1) ) =
189  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 5)
190  xm( iel , convnsy(nul3,nul2) ) =
191  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 6)
192  ENDDO
193  ENDIF
194 !
195  ELSEIF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'S') THEN
196 !
197 ! CASE WHERE M CAN BE ANYTHING AND N IS SYMMETRICAL
198 !
199  IF(ncsize.GT.1) THEN
200  DO k = 1 , neleb
201  iel = nelbor(k)
202  IF(iel.GT.0) THEN
203  nul1=nulone(k)
204  nul2=nulone(k+neleb)
205  nul3=nulone(k+2*neleb)
206  xm( iel , convnsy(nul1,nul2) ) =
207  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 1)
208  xm( iel , convnsy(nul1,nul3) ) =
209  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 2)
210  xm( iel , convnsy(nul2,nul3) ) =
211  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 3)
212  xm( iel , convnsy(nul2,nul1) ) =
213  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 1)
214  xm( iel , convnsy(nul3,nul1) ) =
215  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 2)
216  xm( iel , convnsy(nul3,nul2) ) =
217  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 3)
218  ENDIF
219  ENDDO
220  ELSE
221  DO k = 1 , neleb
222  iel = nelbor(k)
223  nul1=nulone(k)
224  nul2=nulone(k+neleb)
225  nul3=nulone(k+2*neleb)
226  xm( iel , convnsy(nul1,nul2) ) =
227  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 1)
228  xm( iel , convnsy(nul1,nul3) ) =
229  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 2)
230  xm( iel , convnsy(nul2,nul3) ) =
231  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 3)
232  xm( iel , convnsy(nul2,nul1) ) =
233  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 1)
234  xm( iel , convnsy(nul3,nul1) ) =
235  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 2)
236  xm( iel , convnsy(nul3,nul2) ) =
237  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 3)
238  ENDDO
239  ENDIF
240 !
241  ELSEIF(typexm(1:1).EQ.'S'.AND.typexn(1:1).EQ.'S') THEN
242 !
243 ! CASE WHERE BOTH MATRICES ARE SYMMETRICAL
244 !
245  IF(ncsize.GT.1) THEN
246  DO k = 1 , neleb
247  iel = nelbor(k)
248  IF(iel.GT.0) THEN
249  nul1=nulone(k)
250  nul2=nulone(k+neleb)
251  nul3=nulone(k+2*neleb)
252  xm( iel , convsym(nul1,nul2) ) =
253  & xm( iel , convsym(nul1,nul2) ) + xn(k, 1)
254  xm( iel , convsym(nul1,nul3) ) =
255  & xm( iel , convsym(nul1,nul3) ) + xn(k, 2)
256  xm( iel , convsym(nul2,nul3) ) =
257  & xm( iel , convsym(nul2,nul3) ) + xn(k, 3)
258  ENDIF
259  ENDDO
260  ELSE
261  DO k = 1 , neleb
262  iel = nelbor(k)
263  nul1=nulone(k)
264  nul2=nulone(k+neleb)
265  nul3=nulone(k+2*neleb)
266  xm( iel , convsym(nul1,nul2) ) =
267  & xm( iel , convsym(nul1,nul2) ) + xn(k, 1)
268  xm( iel , convsym(nul1,nul3) ) =
269  & xm( iel , convsym(nul1,nul3) ) + xn(k, 2)
270  xm( iel , convsym(nul2,nul3) ) =
271  & xm( iel , convsym(nul2,nul3) ) + xn(k, 3)
272  ENDDO
273  ENDIF
274 !
275  ELSE
276  WRITE(lu,99) typexm(1:1),op(1:8),typexn(1:1)
277 99 FORMAT(1x,'OM3181 (BIEF) : TYPEXM = ',a1,' DOES NOT GO',
278  & /,1x,'FOR THE OPERATION : ',a8,' WITH TYPEXN = ',a1)
279  CALL plante(1)
280  stop
281  ENDIF
282 !
283 !-----------------------------------------------------------------------
284 !
285  ELSEIF(op(1:8).EQ.'M=M+TN ') THEN
286 !
287  CALL ovdb( 'X=X+Y ' , dm , dn , z , c , nbor , neleb )
288 !
289  IF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'Q') THEN
290 !
291 ! CASE WHERE BOTH MATRICES ARE NONSYMMETRICAL
292 !
293  IF(ncsize.GT.1) THEN
294  DO k = 1 , neleb
295  iel = nelbor(k)
296  IF(iel.GT.0) THEN
297  nul1=nulone(k)
298  nul2=nulone(k+neleb)
299  nul3=nulone(k+2*neleb)
300  xm( iel , convnsy(nul1,nul2) ) =
301  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 4)
302  xm( iel , convnsy(nul1,nul3) ) =
303  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 5)
304  xm( iel , convnsy(nul2,nul3) ) =
305  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 6)
306  xm( iel , convnsy(nul2,nul1) ) =
307  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 1)
308  xm( iel , convnsy(nul3,nul1) ) =
309  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 2)
310  xm( iel , convnsy(nul3,nul2) ) =
311  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 3)
312  ENDIF
313  ENDDO
314  ELSE
315  DO k = 1 , neleb
316  iel = nelbor(k)
317  nul1=nulone(k)
318  nul2=nulone(k+neleb)
319  nul3=nulone(k+2*neleb)
320  xm( iel , convnsy(nul1,nul2) ) =
321  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 4)
322  xm( iel , convnsy(nul1,nul3) ) =
323  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 5)
324  xm( iel , convnsy(nul2,nul3) ) =
325  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 6)
326  xm( iel , convnsy(nul2,nul1) ) =
327  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 1)
328  xm( iel , convnsy(nul3,nul1) ) =
329  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 2)
330  xm( iel , convnsy(nul3,nul2) ) =
331  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 3)
332  ENDDO
333  ENDIF
334 !
335  ELSEIF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'S') THEN
336 !
337 ! CASE WHERE M CAN BE ANYTHING AND N IS SYMMETRICAL
338 !
339  IF(ncsize.GT.1) THEN
340  DO k = 1 , neleb
341  iel = nelbor(k)
342  IF(iel.GT.0) THEN
343  nul1=nulone(k)
344  nul2=nulone(k+neleb)
345  nul3=nulone(k+2*neleb)
346  xm( iel , convnsy(nul1,nul2) ) =
347  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 1)
348  xm( iel , convnsy(nul1,nul3) ) =
349  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 2)
350  xm( iel , convnsy(nul2,nul3) ) =
351  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 3)
352  xm( iel , convnsy(nul2,nul1) ) =
353  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 1)
354  xm( iel , convnsy(nul3,nul1) ) =
355  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 2)
356  xm( iel , convnsy(nul3,nul2) ) =
357  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 3)
358  ENDIF
359  ENDDO
360  ELSE
361  DO k = 1 , neleb
362  iel = nelbor(k)
363  nul1=nulone(k)
364  nul2=nulone(k+neleb)
365  nul3=nulone(k+2*neleb)
366  xm( iel , convnsy(nul1,nul2) ) =
367  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 1)
368  xm( iel , convnsy(nul1,nul3) ) =
369  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 2)
370  xm( iel , convnsy(nul2,nul3) ) =
371  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 3)
372  xm( iel , convnsy(nul2,nul1) ) =
373  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 1)
374  xm( iel , convnsy(nul3,nul1) ) =
375  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 2)
376  xm( iel , convnsy(nul3,nul2) ) =
377  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 3)
378  ENDDO
379  ENDIF
380 !
381  ELSEIF(typexm(1:1).EQ.'S'.AND.typexn(1:1).EQ.'S') THEN
382 !
383 ! CASE WHERE BOTH MATRICES ARE SYMMETRICAL
384 !
385  IF(ncsize.GT.1) THEN
386  DO k = 1 , neleb
387  iel = nelbor(k)
388  IF(iel.GT.0) THEN
389  nul1=nulone(k)
390  nul2=nulone(k+neleb)
391  nul3=nulone(k+2*neleb)
392  xm( iel , convsym(nul1,nul2) ) =
393  & xm( iel , convsym(nul1,nul2) ) + xn(k, 1)
394  xm( iel , convsym(nul1,nul3) ) =
395  & xm( iel , convsym(nul1,nul3) ) + xn(k, 2)
396  xm( iel , convsym(nul2,nul3) ) =
397  & xm( iel , convsym(nul2,nul3) ) + xn(k, 3)
398  ENDIF
399  ENDDO
400  ELSE
401  DO k = 1 , neleb
402  iel = nelbor(k)
403  nul1=nulone(k)
404  nul2=nulone(k+neleb)
405  nul3=nulone(k+2*neleb)
406  xm( iel , convsym(nul1,nul2) ) =
407  & xm( iel , convsym(nul1,nul2) ) + xn(k, 1)
408  xm( iel , convsym(nul1,nul3) ) =
409  & xm( iel , convsym(nul1,nul3) ) + xn(k, 2)
410  xm( iel , convsym(nul2,nul3) ) =
411  & xm( iel , convsym(nul2,nul3) ) + xn(k, 3)
412  ENDDO
413  ENDIF
414 !
415  ELSE
416  WRITE(lu,99) typexm(1:1),op(1:8),typexn(1:1)
417  CALL plante(1)
418  stop
419  ENDIF
420 !
421 !-----------------------------------------------------------------------
422 !
423  ELSE
424 !
425  WRITE(lu,71) op
426 71 FORMAT(1x,'OM3181 (BIEF) : UNKNOWN OPERATION : ',a8)
427  CALL plante(1)
428  stop
429 !
430  ENDIF
431 !
432 !-----------------------------------------------------------------------
433 !
434  RETURN
435  END
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
ovdb
subroutine ovdb(OP, X, Y, Z, C, NBOR, NPTFR)
Definition: ovdb.f:7
om3181
subroutine om3181(OP, DM, TYPDIM, XM, TYPEXM, DN, TYPDIN, XN, TYPEXN, C, NULONE, NELBOR, NBOR, NELMAX, SIZDN, NELEB)
Definition: om3181.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