The TELEMAC-MASCARET system  trunk
om5111.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE om5111
3 ! *****************
4 !
5  &(op , dm,typdim,xm,typexm, dn,typdin,xn,typexn,
6  & sizdn,szmdn,sizxn,netage, nelmax3d)
7 !
8 !***********************************************************************
9 ! BIEF V6P1 21/08/2010
10 !***********************************************************************
11 !
12 !brief OPERATIONS ON MATRICES.
13 !code
14 !+ M: TETRAHEDRON
15 !+ N: TRIANGLE MATRIX (BOTTOM OR SURFACE)
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 !+ 28/08/02
39 !+ V5P3
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 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
55 !| DM |<->| DIAGONAL OF M
56 !| DN |-->| DIAGONAL OF N
57 !| NELMAX3D |-->| MAXIMUM NUMBER OF 3D ELEMENTS
58 !| NETAGE |-->| NUMBER OF PLANES - 1
59 !| OP |-->| OPERATION TO BE DONE (SEE ABOVE)
60 !| SIZDN |-->| SIZE OF DIAGONAL DN
61 !| SIZXN |-->| SIZE OF OFF-DIAGONAL TERMS XN
62 !| SZMDN |-->| MAXIMUM SIZE OF DIAGONAL DN
63 !| TYPDIM |<->| TYPE OF DIAGONAL OF M:
64 !| | | TYPDIM = 'Q' : ANY VALUE
65 !| | | TYPDIM = 'I' : IDENTITY
66 !| | | TYPDIM = '0' : ZERO
67 !| TYPDIN |<->| TYPE OF DIAGONAL OF N:
68 !| | | TYPDIN = 'Q' : ANY VALUE
69 !| | | TYPDIN = 'I' : IDENTITY
70 !| | | TYPDIN = '0' : ZERO
71 !| TYPEXM |-->| TYPE OF OFF-DIAGONAL TERMS OF M:
72 !| | | TYPEXM = 'Q' : ANY VALUE
73 !| | | TYPEXM = 'S' : SYMMETRIC
74 !| | | TYPEXM = '0' : ZERO
75 !| TYPEXN |-->| TYPE OF OFF-DIAGONAL TERMS OF N:
76 !| | | TYPEXN = 'Q' : ANY VALUE
77 !| | | TYPEXN = 'S' : SYMMETRIC
78 !| | | TYPEXN = '0' : ZERO
79 !| XM |-->| OFF-DIAGONAL TERMS OF M
80 !| XN |-->| OFF-DIAGONAL TERMS OF N
81 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
82 !
83  USE bief
84 !
86  IMPLICIT NONE
87 !
88 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
89 !
90  INTEGER, INTENT(IN) :: NETAGE,SIZDN,SZMDN,SIZXN
91  INTEGER, INTENT(IN) :: NELMAX3D
92  CHARACTER(LEN=8), INTENT(IN) :: OP
93  DOUBLE PRECISION, INTENT(IN) :: DN(*),XN(nelmax3d/(3*netage),*)
94  DOUBLE PRECISION, INTENT(INOUT) :: DM(szmdn,*)
95  DOUBLE PRECISION,INTENT(INOUT)::XM(nelmax3d/(3*netage),3,netage,*)
96  CHARACTER(LEN=1), INTENT(INOUT) :: TYPDIM,TYPEXM,TYPDIN,TYPEXN
97 !
98 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
99 !
100  INTEGER K
101 !
102 !-----------------------------------------------------------------------
103 !
104  IF(op(1:8).EQ.'M=M+NF ') THEN
105 !
106  IF(typdim.EQ.'Q'.AND.typdin.EQ.'Q') THEN
107  CALL ov('X=X+Y ', x=dm, y=dn, dim1=sizdn)
108  ELSE
109  WRITE(lu,199) typdim(1:1),op(1:8),typdin(1:1)
110 199 FORMAT(1x,'OM5111 (BIEF) : TYPDIM = ',a1,' NOT IMPLEMENTED',
111  & /,1x,'FOR THE OPERATION : ',a8,' WITH TYPDIN = ',a1)
112  CALL plante(1)
113  stop
114  ENDIF
115 !
116  IF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'Q') THEN
117 !
118 ! CASE WHERE BOTH MATRICES ARE NONSYMMETRICAL
119 !
120 ! XM(K,1,1, ) K : TRIANGLE NUMBER
121 ! 1 : T1 TETRAHEDRON, ONE SIDE OF WHICH IS AT THE BOTTOM
122 ! 1 : 1ST LAYER, THAT AT THE BOTTOM
123 !
124  DO k = 1 , sizxn
125  xm(k,1,1,01) = xm(k,1,1,01) + xn(k,1)
126  xm(k,1,1,02) = xm(k,1,1,02) + xn(k,2)
127  xm(k,1,1,04) = xm(k,1,1,04) + xn(k,3)
128  xm(k,1,1,07) = xm(k,1,1,07) + xn(k,4)
129  xm(k,1,1,08) = xm(k,1,1,08) + xn(k,5)
130  xm(k,1,1,10) = xm(k,1,1,10) + xn(k,6)
131  ENDDO
132 !
133  ELSEIF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'S') THEN
134 !
135 ! CASE WHERE M CAN BE ANYTHING AND N IS SYMMETRICAL
136 !
137  DO k = 1 , sizxn
138  xm(k,1,1,01) = xm(k,1,1,01) + xn(k,1)
139  xm(k,1,1,02) = xm(k,1,1,02) + xn(k,2)
140  xm(k,1,1,04) = xm(k,1,1,04) + xn(k,3)
141  xm(k,1,1,07) = xm(k,1,1,07) + xn(k,1)
142  xm(k,1,1,08) = xm(k,1,1,08) + xn(k,2)
143  xm(k,1,1,10) = xm(k,1,1,10) + xn(k,3)
144  ENDDO
145 !
146  ELSEIF(typexm(1:1).EQ.'S'.AND.typexn(1:1).EQ.'S') THEN
147 !
148 ! CASE WHERE BOTH MATRICES ARE SYMMETRICAL
149 !
150  DO k = 1 , sizxn
151  xm(k,1,1,01) = xm(k,1,1,01) + xn(k,1)
152  xm(k,1,1,02) = xm(k,1,1,02) + xn(k,2)
153  xm(k,1,1,04) = xm(k,1,1,04) + xn(k,3)
154  ENDDO
155 !
156  ELSE
157  WRITE(lu,99) typexm(1:1),op(1:8),typexn(1:1)
158 99 FORMAT(1x,'OM5111 (BIEF) : TYPEXM = ',a1,' DOES NOT GO',
159  & /,1x,'FOR THE OPERATION : ',a8,' WITH TYPEXN = ',a1)
160  CALL plante(1)
161  stop
162  ENDIF
163 !
164 !-----------------------------------------------------------------------
165 !
166  ELSEIF(op(1:8).EQ.'M=M+TNF ') THEN
167 !
168  CALL ov('X=X+Y ', x=dm, y=dn, dim1=sizdn)
169 !
170  IF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'Q') THEN
171 !
172 ! CASE WHERE BOTH MATRICES ARE NONSYMMETRICAL
173 !
174  DO k = 1 , sizxn
175  xm(k,1,1,01) = xm(k,1,1,01) + xn(k,4)
176  xm(k,1,1,02) = xm(k,1,1,02) + xn(k,5)
177  xm(k,1,1,04) = xm(k,1,1,04) + xn(k,6)
178  xm(k,1,1,07) = xm(k,1,1,07) + xn(k,1)
179  xm(k,1,1,08) = xm(k,1,1,08) + xn(k,2)
180  xm(k,1,1,10) = xm(k,1,1,10) + xn(k,3)
181  ENDDO
182 !
183  ELSEIF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'S') THEN
184 !
185 ! CASE WHERE M CAN BE ANYTHING AND N IS SYMMETRICAL
186 !
187  DO k = 1 , sizxn
188  xm(k,1,1,01) = xm(k,1,1,01) + xn(k,1)
189  xm(k,1,1,02) = xm(k,1,1,02) + xn(k,2)
190  xm(k,1,1,04) = xm(k,1,1,04) + xn(k,3)
191  xm(k,1,1,07) = xm(k,1,1,07) + xn(k,1)
192  xm(k,1,1,08) = xm(k,1,1,08) + xn(k,2)
193  xm(k,1,1,10) = xm(k,1,1,10) + xn(k,3)
194  ENDDO
195 !
196  ELSEIF(typexm(1:1).EQ.'S'.AND.typexn(1:1).EQ.'S') THEN
197 !
198 ! CASE WHERE BOTH MATRICES ARE SYMMETRICAL
199 !
200  DO k = 1 , sizxn
201  xm(k,1,1,01) = xm(k,1,1,01) + xn(k,1)
202  xm(k,1,1,02) = xm(k,1,1,02) + xn(k,2)
203  xm(k,1,1,04) = xm(k,1,1,04) + xn(k,3)
204  ENDDO
205 !
206  ELSE
207  WRITE(lu,99) typexm(1:1),op(1:8),typexn(1:1)
208  CALL plante(1)
209  stop
210  ENDIF
211 !
212 !-----------------------------------------------------------------------
213 !
214  ELSEIF(op(1:8).EQ.'M=M+NS ') THEN
215 !
216  CALL ov('X=X+Y ', x=dm(1,netage+1), y=dn, dim1=sizdn)
217 !
218  IF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'Q') THEN
219 !
220 ! CASE WHERE BOTH MATRICES ARE NONSYMMETRICAL
221 !
222 ! XM(K,1,1, ) K : TRIANGLE NUMBER
223 ! 2 : T2 TETRAHEDRON, ONE SIDE OF WHICH IS AT THE SURFACE
224 ! NETAGE : LAST LAYER, THAT AT THE SURFACE
225 !
226  DO k = 1 , sizxn
227  xm(k,2,netage,02) = xm(k,2,netage,02) + xn(k,1)
228  xm(k,2,netage,01) = xm(k,2,netage,01) + xn(k,2)
229  xm(k,2,netage,10) = xm(k,2,netage,10) + xn(k,3)
230  xm(k,2,netage,08) = xm(k,2,netage,08) + xn(k,4)
231  xm(k,2,netage,07) = xm(k,2,netage,07) + xn(k,5)
232  xm(k,2,netage,04) = xm(k,2,netage,04) + xn(k,6)
233  ENDDO
234 !
235  ELSEIF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'S') THEN
236 !
237 ! CASE WHERE M CAN BE ANYTHING AND N IS SYMMETRICAL
238 !
239  DO k = 1 , sizxn
240  xm(k,2,netage,02) = xm(k,2,netage,02) + xn(k,1)
241  xm(k,2,netage,01) = xm(k,2,netage,01) + xn(k,2)
242  xm(k,2,netage,10) = xm(k,2,netage,10) + xn(k,3)
243  xm(k,2,netage,08) = xm(k,2,netage,08) + xn(k,1)
244  xm(k,2,netage,07) = xm(k,2,netage,07) + xn(k,2)
245  xm(k,2,netage,04) = xm(k,2,netage,04) + xn(k,3)
246  ENDDO
247 !
248  ELSEIF(typexm(1:1).EQ.'S'.AND.typexn(1:1).EQ.'S') THEN
249 !
250 ! CASE WHERE BOTH MATRICES ARE SYMMETRICAL
251 !
252  DO k = 1 , sizxn
253  xm(k,2,netage,02) = xm(k,2,netage,02) + xn(k,1)
254  xm(k,2,netage,01) = xm(k,2,netage,01) + xn(k,2)
255  xm(k,2,netage,04) = xm(k,2,netage,04) + xn(k,3)
256  ENDDO
257 !
258  ELSE
259  WRITE(lu,99) typexm(1:1),op(1:8),typexn(1:1)
260  CALL plante(1)
261  stop
262  ENDIF
263 !
264 !-----------------------------------------------------------------------
265 !
266  ELSEIF(op(1:8).EQ.'M=M+TNS ') THEN
267 !
268  CALL ov('X=X+Y ', x=dm(1,netage+1), y=dn, dim1=sizdn)
269 !
270  IF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'Q') THEN
271 !
272 ! CASE WHERE BOTH MATRICES ARE NONSYMMETRICAL
273 !
274  DO k = 1 , sizxn
275  xm(k,2,netage,02) = xm(k,2,netage,02) + xn(k,4)
276  xm(k,2,netage,01) = xm(k,2,netage,01) + xn(k,5)
277  xm(k,2,netage,10) = xm(k,2,netage,10) + xn(k,6)
278  xm(k,2,netage,08) = xm(k,2,netage,08) + xn(k,1)
279  xm(k,2,netage,07) = xm(k,2,netage,07) + xn(k,2)
280  xm(k,2,netage,04) = xm(k,2,netage,04) + xn(k,3)
281  ENDDO
282 !
283  ELSEIF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'S') THEN
284 !
285 ! CASE WHERE M CAN BE ANYTHING AND N IS SYMMETRICAL
286 !
287  DO k = 1 , sizxn
288  xm(k,2,netage,02) = xm(k,2,netage,02) + xn(k,1)
289  xm(k,2,netage,01) = xm(k,2,netage,01) + xn(k,2)
290  xm(k,2,netage,10) = xm(k,2,netage,10) + xn(k,3)
291  xm(k,2,netage,08) = xm(k,2,netage,08) + xn(k,1)
292  xm(k,2,netage,07) = xm(k,2,netage,07) + xn(k,2)
293  xm(k,2,netage,04) = xm(k,2,netage,04) + xn(k,3)
294  ENDDO
295 !
296  ELSEIF(typexm(1:1).EQ.'S'.AND.typexn(1:1).EQ.'S') THEN
297 !
298 ! CASE WHERE BOTH MATRICES ARE SYMMETRICAL
299 !
300  DO k = 1 , sizxn
301  xm(k,2,netage,02) = xm(k,2,netage,02) + xn(k,1)
302  xm(k,2,netage,01) = xm(k,2,netage,01) + xn(k,2)
303  xm(k,2,netage,04) = xm(k,2,netage,04) + xn(k,3)
304  ENDDO
305 !
306  ELSE
307  WRITE(lu,99) typexm(1:1),op(1:8),typexn(1:1)
308  CALL plante(1)
309  stop
310  ENDIF
311 !
312 !-----------------------------------------------------------------------
313 !
314  ELSE
315 !
316  WRITE(lu,71) op
317 71 FORMAT(1x,'OM5111 (BIEF) : UNKNOWN OPERATION : ',a8)
318  CALL plante(1)
319  stop
320 !
321  ENDIF
322 !
323 !-----------------------------------------------------------------------
324 !
325  RETURN
326  END
subroutine ov(OP, X, Y, Z, C, DIM1)
Definition: ov.f:7
subroutine om5111(OP, DM, TYPDIM, XM, TYPEXM, DN, TYPDIN, XN, TYPEXN, SIZDN, SZMDN, SIZXN, NETAGE, NELMAX3D)
Definition: om5111.f:8
Definition: bief.f:3