The TELEMAC-MASCARET system  trunk
om5161.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE om5161
3 ! *****************
4 !
5  &(op , dm,typdim,xm,typexm, dn,typdin,xn,typexn, c,
6  & nulone,nelbor,nbor,nelmax,sizdn,sizxn,szmxn)
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 !| NELBOR |-->| FOR THE KTH BOUNDARY EDGE, GIVES THE CORRESPONDING
64 !| | | ELEMENT.
65 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
66 !| NULONE |-->| GOES WITH ARRAY NELBOR. NELBOR GIVES THE
67 !| | | ADJACENT ELEMENT, NULONE GIVES THE LOCAL
68 !| | | NUMBER OF THE FIRST NODE OF THE BOUNDARY EDGE
69 !| | | I.E. 1, 2 OR 3 FOR TRIANGLES.
70 !| OP |-->| OPERATION TO BE DONE (SEE ABOVE)
71 !| SIZDN |-->| SIZE OF DIAGONAL DN
72 !| SIZXN |-->| SIZE OF OFF-DIAGONAL TERMS XN
73 !| SZMXN |-->| MAXIMUM SIZE OF OFF-DIAGONAL TERMS XN
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 !
95  USE bief, ex_om5161 => om5161
96 !
97  USE declarations_special
98  IMPLICIT NONE
99 !
100 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
101 !
102  INTEGER, INTENT(IN) :: NELMAX,SIZDN,SIZXN,SZMXN
103  CHARACTER(LEN=8), INTENT(IN) :: OP
104  INTEGER, INTENT(IN) :: NULONE(szmxn,3)
105  INTEGER, INTENT(IN) :: NELBOR(*),NBOR(*)
106  DOUBLE PRECISION, INTENT(IN) :: DN(*),XN(szmxn,*)
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 K,IEL,NUL1,NUL2,NUL3
114 !
115  DOUBLE PRECISION Z(1)
116 !
117 !-----------------------------------------------------------------------
118 !
119 ! BEWARE: IN FORTRAN THE FIRST INDEX VARIES MOST QUICKLY
120 ! 123456789 SHOULD NOT BE USED
121 !
122  INTEGER :: CONVNSY(4,4)
123  parameter( convnsy = reshape( (/
124  & 123456789, 7 , 8 , 9 ,
125  & 1 , 123456789, 10 , 11 ,
126  & 2 , 4 , 123456789, 12 ,
127  & 3 , 5 , 6 , 123456789 /),
128  & shape=(/ 4,4 /) ) )
129  INTEGER :: CONVSYM(4,4)
130  parameter( convsym = reshape( (/
131  & 123456789, 1 , 2 , 3 ,
132  & 1 , 123456789, 4 , 5 ,
133  & 2 , 4 , 123456789, 6 ,
134  & 3 , 5 , 6 , 123456789 /),
135  & shape=(/ 4,4 /) ) )
136 !
137 !***********************************************************************
138 !
139  IF(op(1:8).EQ.'M=M+N ') THEN
140 !
141  IF(typdim.EQ.'Q'.AND.typdin.EQ.'Q') THEN
142  CALL ovdb( 'X=X+Y ' , dm , dn , z , c , nbor , sizdn )
143  ELSE
144  WRITE(lu,199) typdim(1:1),op(1:8),typdin(1:1)
145 199 FORMAT(1x,'OM5161 (BIEF) : TYPDIM = ',a1,' NOT IMPLEMENTED',
146  & /,1x,'FOR THE OPERATION : ',a8,' WITH TYPDIN = ',a1)
147  CALL plante(1)
148  stop
149  ENDIF
150 !
151  IF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'Q') THEN
152 !
153 ! CASE WHERE BOTH MATRICES ARE NONSYMMETRICAL
154 !
155  IF(ncsize.GT.1) THEN
156  DO k = 1 , sizxn
157  iel = nelbor(k)
158  IF(iel.GT.0) THEN
159  nul1=nulone(k,1)
160  nul2=nulone(k,2)
161  nul3=nulone(k,3)
162  xm( iel , convnsy(nul1,nul2) ) =
163  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 1)
164  xm( iel , convnsy(nul1,nul3) ) =
165  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 2)
166  xm( iel , convnsy(nul2,nul3) ) =
167  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 3)
168  xm( iel , convnsy(nul2,nul1) ) =
169  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 4)
170  xm( iel , convnsy(nul3,nul1) ) =
171  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 5)
172  xm( iel , convnsy(nul3,nul2) ) =
173  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 6)
174  ENDIF
175  ENDDO
176  ELSE
177  DO k = 1 , sizxn
178  iel = nelbor(k)
179  nul1=nulone(k,1)
180  nul2=nulone(k,2)
181  nul3=nulone(k,3)
182  xm( iel , convnsy(nul1,nul2) ) =
183  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 1)
184  xm( iel , convnsy(nul1,nul3) ) =
185  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 2)
186  xm( iel , convnsy(nul2,nul3) ) =
187  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 3)
188  xm( iel , convnsy(nul2,nul1) ) =
189  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 4)
190  xm( iel , convnsy(nul3,nul1) ) =
191  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 5)
192  xm( iel , convnsy(nul3,nul2) ) =
193  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 6)
194  ENDDO
195  ENDIF
196 !
197  ELSEIF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'S') THEN
198 !
199 ! CASE WHERE M CAN BE ANYTHING AND N IS SYMMETRICAL
200 !
201  IF(ncsize.GT.1) THEN
202  DO k = 1 , sizxn
203  iel = nelbor(k)
204  IF(iel.GT.0) THEN
205  nul1=nulone(k,1)
206  nul2=nulone(k,2)
207  nul3=nulone(k,3)
208  xm( iel , convnsy(nul1,nul2) ) =
209  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 1)
210  xm( iel , convnsy(nul1,nul3) ) =
211  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 2)
212  xm( iel , convnsy(nul2,nul3) ) =
213  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 3)
214  xm( iel , convnsy(nul2,nul1) ) =
215  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 1)
216  xm( iel , convnsy(nul3,nul1) ) =
217  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 2)
218  xm( iel , convnsy(nul3,nul2) ) =
219  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 3)
220  ENDIF
221  ENDDO
222  ELSE
223  DO k = 1 , sizxn
224  iel = nelbor(k)
225  nul1=nulone(k,1)
226  nul2=nulone(k,2)
227  nul3=nulone(k,3)
228  xm( iel , convnsy(nul1,nul2) ) =
229  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 1)
230  xm( iel , convnsy(nul1,nul3) ) =
231  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 2)
232  xm( iel , convnsy(nul2,nul3) ) =
233  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 3)
234  xm( iel , convnsy(nul2,nul1) ) =
235  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 1)
236  xm( iel , convnsy(nul3,nul1) ) =
237  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 2)
238  xm( iel , convnsy(nul3,nul2) ) =
239  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 3)
240  ENDDO
241  ENDIF
242 !
243  ELSEIF(typexm(1:1).EQ.'S'.AND.typexn(1:1).EQ.'S') THEN
244 !
245 ! CASE WHERE BOTH MATRICES ARE SYMMETRICAL
246 !
247  IF(ncsize.GT.1) THEN
248  DO k = 1 , sizxn
249  iel = nelbor(k)
250  IF(iel.GT.0) THEN
251  nul1=nulone(k,1)
252  nul2=nulone(k,2)
253  nul3=nulone(k,3)
254  xm( iel , convsym(nul1,nul2) ) =
255  & xm( iel , convsym(nul1,nul2) ) + xn(k, 1)
256  xm( iel , convsym(nul1,nul3) ) =
257  & xm( iel , convsym(nul1,nul3) ) + xn(k, 2)
258  xm( iel , convsym(nul2,nul3) ) =
259  & xm( iel , convsym(nul2,nul3) ) + xn(k, 3)
260  ENDIF
261  ENDDO
262  ELSE
263  DO k = 1 , sizxn
264  iel = nelbor(k)
265  nul1=nulone(k,1)
266  nul2=nulone(k,2)
267  nul3=nulone(k,3)
268  xm( iel , convsym(nul1,nul2) ) =
269  & xm( iel , convsym(nul1,nul2) ) + xn(k, 1)
270  xm( iel , convsym(nul1,nul3) ) =
271  & xm( iel , convsym(nul1,nul3) ) + xn(k, 2)
272  xm( iel , convsym(nul2,nul3) ) =
273  & xm( iel , convsym(nul2,nul3) ) + xn(k, 3)
274  ENDDO
275  ENDIF
276 !
277  ELSE
278  WRITE(lu,99) typexm(1:1),op(1:8),typexn(1:1)
279 99 FORMAT(1x,'OM5161 (BIEF) : TYPEXM = ',a1,' DOES NOT GO',
280  & /,1x,'FOR THE OPERATION : ',a8,' WITH TYPEXN = ',a1)
281  CALL plante(1)
282  stop
283  ENDIF
284 !
285 !-----------------------------------------------------------------------
286 !
287  ELSEIF(op(1:8).EQ.'M=M+TN ') THEN
288 !
289  CALL ovdb( 'X=X+Y ' , dm , dn , z , c , nbor , sizxn )
290 !
291  IF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'Q') THEN
292 !
293 ! CASE WHERE BOTH MATRICES ARE NONSYMMETRICAL
294 !
295  IF(ncsize.GT.1) THEN
296  DO k = 1 , sizxn
297  iel = nelbor(k)
298  IF(iel.GT.0) THEN
299  nul1=nulone(k,1)
300  nul2=nulone(k,2)
301  nul3=nulone(k,3)
302  xm( iel , convnsy(nul1,nul2) ) =
303  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 4)
304  xm( iel , convnsy(nul1,nul3) ) =
305  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 5)
306  xm( iel , convnsy(nul2,nul3) ) =
307  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 6)
308  xm( iel , convnsy(nul2,nul1) ) =
309  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 1)
310  xm( iel , convnsy(nul3,nul1) ) =
311  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 2)
312  xm( iel , convnsy(nul3,nul2) ) =
313  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 3)
314  ENDIF
315  ENDDO
316  ELSE
317  DO k = 1 , sizxn
318  iel = nelbor(k)
319  nul1=nulone(k,1)
320  nul2=nulone(k,2)
321  nul3=nulone(k,3)
322  xm( iel , convnsy(nul1,nul2) ) =
323  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 4)
324  xm( iel , convnsy(nul1,nul3) ) =
325  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 5)
326  xm( iel , convnsy(nul2,nul3) ) =
327  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 6)
328  xm( iel , convnsy(nul2,nul1) ) =
329  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 1)
330  xm( iel , convnsy(nul3,nul1) ) =
331  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 2)
332  xm( iel , convnsy(nul3,nul2) ) =
333  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 3)
334  ENDDO
335  ENDIF
336 !
337  ELSEIF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'S') THEN
338 !
339 ! CASE WHERE M CAN BE ANYTHING AND N IS SYMMETRICAL
340 !
341  IF(ncsize.GT.1) THEN
342  DO k = 1 , sizxn
343  iel = nelbor(k)
344  IF(iel.GT.0) THEN
345  nul1=nulone(k,1)
346  nul2=nulone(k,2)
347  nul3=nulone(k,3)
348  xm( iel , convnsy(nul1,nul2) ) =
349  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 1)
350  xm( iel , convnsy(nul1,nul3) ) =
351  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 2)
352  xm( iel , convnsy(nul2,nul3) ) =
353  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 3)
354  xm( iel , convnsy(nul2,nul1) ) =
355  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 1)
356  xm( iel , convnsy(nul3,nul1) ) =
357  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 2)
358  xm( iel , convnsy(nul3,nul2) ) =
359  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 3)
360  ENDIF
361  ENDDO
362  ELSE
363  DO k = 1 , sizxn
364  iel = nelbor(k)
365  nul1=nulone(k,1)
366  nul2=nulone(k,2)
367  nul3=nulone(k,3)
368  xm( iel , convnsy(nul1,nul2) ) =
369  & xm( iel , convnsy(nul1,nul2) ) + xn(k, 1)
370  xm( iel , convnsy(nul1,nul3) ) =
371  & xm( iel , convnsy(nul1,nul3) ) + xn(k, 2)
372  xm( iel , convnsy(nul2,nul3) ) =
373  & xm( iel , convnsy(nul2,nul3) ) + xn(k, 3)
374  xm( iel , convnsy(nul2,nul1) ) =
375  & xm( iel , convnsy(nul2,nul1) ) + xn(k, 1)
376  xm( iel , convnsy(nul3,nul1) ) =
377  & xm( iel , convnsy(nul3,nul1) ) + xn(k, 2)
378  xm( iel , convnsy(nul3,nul2) ) =
379  & xm( iel , convnsy(nul3,nul2) ) + xn(k, 3)
380  ENDDO
381  ENDIF
382 !
383  ELSEIF(typexm(1:1).EQ.'S'.AND.typexn(1:1).EQ.'S') THEN
384 !
385 ! CASE WHERE BOTH MATRICES ARE SYMMETRICAL
386 !
387  IF(ncsize.GT.1) THEN
388  DO k = 1 , sizxn
389  iel = nelbor(k)
390  IF(iel.GT.0) THEN
391  nul1=nulone(k,1)
392  nul2=nulone(k,2)
393  nul3=nulone(k,3)
394  xm( iel , convsym(nul1,nul2) ) =
395  & xm( iel , convsym(nul1,nul2) ) + xn(k, 1)
396  xm( iel , convsym(nul1,nul3) ) =
397  & xm( iel , convsym(nul1,nul3) ) + xn(k, 2)
398  xm( iel , convsym(nul2,nul3) ) =
399  & xm( iel , convsym(nul2,nul3) ) + xn(k, 3)
400  ENDIF
401  ENDDO
402  ELSE
403  DO k = 1 , sizxn
404  iel = nelbor(k)
405  nul1=nulone(k,1)
406  nul2=nulone(k,2)
407  nul3=nulone(k,3)
408  xm( iel , convsym(nul1,nul2) ) =
409  & xm( iel , convsym(nul1,nul2) ) + xn(k, 1)
410  xm( iel , convsym(nul1,nul3) ) =
411  & xm( iel , convsym(nul1,nul3) ) + xn(k, 2)
412  xm( iel , convsym(nul2,nul3) ) =
413  & xm( iel , convsym(nul2,nul3) ) + xn(k, 3)
414  ENDDO
415  ENDIF
416 !
417  ELSE
418  WRITE(lu,99) typexm(1:1),op(1:8),typexn(1:1)
419  CALL plante(1)
420  stop
421  ENDIF
422 !
423 !-----------------------------------------------------------------------
424 !
425  ELSE
426 !
427  WRITE(lu,71) op
428 71 FORMAT(1x,'OM5161 (BIEF) : UNKNOWN OPERATION : ',a8)
429  CALL plante(1)
430  stop
431 !
432  ENDIF
433 !
434 !-----------------------------------------------------------------------
435 !
436  RETURN
437  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
om5161
subroutine om5161(OP, DM, TYPDIM, XM, TYPEXM, DN, TYPDIN, XN, TYPEXN, C, NULONE, NELBOR, NBOR, NELMAX, SIZDN, SIZXN, SZMXN)
Definition: om5161.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