The TELEMAC-MASCARET system  trunk
mt05aa.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt05aa
3 ! *****************
4 !
5  &( a11 , a12 , a13 ,
6  & a21 , a22 , a23 ,
7  & a31 , a32 , a33 ,
8  & xmul,su,sv,u,v,
9  & xel,yel,ikle,nelem,nelmax,formul)
10 !
11 !***********************************************************************
12 ! BIEF V6P1 21/08/2010
13 !***********************************************************************
14 !
15 !brief COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
16 !code
17 !+ / -> --->
18 !+ A(I,J) = XMUL / PSI1(I) * U . GRAD(PSI2(J)) D(OMEGA)
19 !+ /OMEGA
20 !+
21 !+ PSI1: BASES OF TYPE P1 TRIANGLE
22 !+ PSI2: BASES OF TYPE P1 TRIANGLE
23 !
24 !warning THE JACOBIAN MUST BE POSITIVE
25 !
26 !history J-M HERVOUET (LNH)
27 !+ 05/02/91
28 !+ V5P1
29 !+
30 !
31 !history A FROEHLY (MATMECA)
32 !+ 01/07/08
33 !+ V5P9
34 !+
35 !
36 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
37 !+ 13/07/2010
38 !+ V6P0
39 !+ Translation of French comments within the FORTRAN sources into
40 !+ English comments
41 !
42 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
43 !+ 21/08/2010
44 !+ V6P0
45 !+ Creation of DOXYGEN tags for automated documentation and
46 !+ cross-referencing of the FORTRAN sources
47 !
48 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
49 !| A11 |<--| ELEMENTS OF MATRIX
50 !| A12 |<--| ELEMENTS OF MATRIX
51 !| A13 |<--| ELEMENTS OF MATRIX
52 !| A21 |<--| ELEMENTS OF MATRIX
53 !| A22 |<--| ELEMENTS OF MATRIX
54 !| A23 |<--| ELEMENTS OF MATRIX
55 !| A31 |<--| ELEMENTS OF MATRIX
56 !| A32 |<--| ELEMENTS OF MATRIX
57 !| A33 |<--| ELEMENTS OF MATRIX
58 !| FORMUL |-->| FORMULA DESCRIBING THE RESULTING MATRIX
59 !| IKLE1 |-->| FIRST POINTS OF TRIANGLES
60 !| IKLE2 |-->| SECOND POINTS OF TRIANGLES
61 !| IKLE3 |-->| THIRD POINTS OF TRIANGLES
62 !| IKLE4 |-->| QUASI-BUBBLE POINT
63 !| NELEM |-->| NUMBER OF ELEMENTS
64 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
65 !| SU |-->| BIEF_OBJ STRUCTURE OF U
66 !| SURFAC |-->| AREA OF TRIANGLES
67 !| SV |-->| BIEF_OBJ STRUCTURE OF V
68 !| U |-->| FUNCTION U USED IN THE FORMULA
69 !| V |-->| FUNCTION V USED IN THE FORMULA
70 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
71 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
72 !| XMUL |-->| MULTIPLICATION FACTOR
73 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
74 !
75  USE bief, ex_mt05aa => mt05aa
76 !
77  USE declarations_special
78  IMPLICIT NONE
79 !
80 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
81 !
82  INTEGER, INTENT(IN) :: NELEM,NELMAX
83  INTEGER, INTENT(IN) :: IKLE(nelmax,*)
84  DOUBLE PRECISION, INTENT(INOUT) :: A11(*),A12(*),A13(*)
85  DOUBLE PRECISION, INTENT(INOUT) :: A21(*),A22(*),A23(*)
86  DOUBLE PRECISION, INTENT(INOUT) :: A31(*),A32(*),A33(*)
87  DOUBLE PRECISION, INTENT(IN) :: XMUL,U(*),V(*)
88  TYPE(bief_obj), INTENT(IN) :: SU,SV
89  CHARACTER(LEN=16) :: FORMUL
90  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,3),YEL(nelmax,3)
91 !
92 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
93 !
94 ! DECLARATIONS SPECIFIC TO THIS SUBROUTINE
95 !
96  INTEGER IELMU,IELMV,IELEM
97 !
98  DOUBLE PRECISION SUR24,SUR120,SUR216
99  DOUBLE PRECISION X2,X3,Y2,Y3
100  DOUBLE PRECISION U1,U2,U3,U4,U5,U6,V1,V2,V3,V4,V5,V6
101  DOUBLE PRECISION U123,V123,COU1,COV1,COU2,COV2,COU3,COV3
102  DOUBLE PRECISION QUATRU,QUATRV,SUR6,USUR2,VSUR2
103  DOUBLE PRECISION K1,K2,K3,L12,L13,L21,L23,L31,L32
104 !
105  INTRINSIC max,min
106 !
107 !-----------------------------------------------------------------------
108 !
109  sur6 = xmul/ 6.d0
110  sur24 = xmul/ 24.d0
111  sur120 = xmul/120.d0
112  sur216 = xmul/216.d0
113 !
114 !-----------------------------------------------------------------------
115 !
116  ielmu = su%ELM
117  ielmv = sv%ELM
118 !
119 !-----------------------------------------------------------------------
120 !
121 ! CASE WHERE U AND V ARE CONSTANT BY ELEMENT
122 !
123  IF(ielmu.EQ.10.AND.ielmv.EQ.10) THEN
124 !
125 ! LOOP ON THE ELEMENTS
126 !
127  DO ielem = 1 , nelem
128 !
129  x2 = xel(ielem,2) * sur24
130  x3 = xel(ielem,3) * sur24
131  y2 = yel(ielem,2) * sur24
132  y3 = yel(ielem,3) * sur24
133 !
134  quatru = 4 * u(ielem)
135  quatrv = 4 * v(ielem)
136 !
137 ! DIAGONAL TERMS
138 !
139  a11(ielem) = (x3-x2) * quatrv + (y2-y3) * quatru
140  a22(ielem) = -x3 * quatrv +y3 * quatru
141  a33(ielem) = x2 * quatrv - y2 * quatru
142 !
143 ! EXTRADIAGONAL TERMS
144 !
145  a12(ielem) = -x3 * quatrv + y3 * quatru
146  a13(ielem) = x2 * quatrv - y2 * quatru
147  a23(ielem) = x2 * quatrv - y2 * quatru
148  a21(ielem) = (x3-x2) * quatrv + (y2-y3) * quatru
149  a31(ielem) = (x3-x2) * quatrv + (y2-y3) * quatru
150  a32(ielem) = -x3 * quatrv + y3 * quatru
151 !
152  ENDDO ! IELEM
153 !
154 !-----------------------------------------------------------------------
155 !
156 ! CASE WHERE U AND V ARE LINEAR, QUASI-BUBBLE OR QUADRATIC AND N SCHEME
157 !
158  ELSEIF(formul(16:16).EQ.'N' .AND.
159  & ( (ielmu.EQ.11.AND.ielmv.EQ.11).OR.
160  & (ielmu.EQ.12.AND.ielmv.EQ.12).OR.
161  & (ielmu.EQ.13.AND.ielmv.EQ.13) ) ) THEN
162 !
163 ! N SCHEME: U AND V ARE TREATED AS IF LINEAR
164 !
165  DO ielem = 1 , nelem
166 !
167  x2 = xel(ielem,2)
168  x3 = xel(ielem,3)
169  y2 = yel(ielem,2)
170  y3 = yel(ielem,3)
171 !
172  u1 = u(ikle(ielem,1))
173  u2 = u(ikle(ielem,2))
174  u3 = u(ikle(ielem,3))
175  v1 = v(ikle(ielem,1))
176  v2 = v(ikle(ielem,2))
177  v3 = v(ikle(ielem,3))
178 !
179  usur2 = (u1+u2+u3)*sur6
180  vsur2 = (v1+v2+v3)*sur6
181 !
182  k1 = usur2 * (y2-y3) - vsur2 * (x2-x3)
183  k2 = usur2 * (y3 ) - vsur2 * (x3 )
184  k3 = usur2 * ( -y2) - vsur2 * ( -x2)
185 !
186  l12 = max( min(k1,-k2) , 0.d0 )
187  l13 = max( min(k1,-k3) , 0.d0 )
188  l21 = max( min(k2,-k1) , 0.d0 )
189  l23 = max( min(k2,-k3) , 0.d0 )
190  l31 = max( min(k3,-k1) , 0.d0 )
191  l32 = max( min(k3,-k2) , 0.d0 )
192 !
193 ! DIAGONAL TERMS
194 !
195  a11(ielem) = l12 + l13
196  a22(ielem) = l21 + l23
197  a33(ielem) = l31 + l32
198 !
199 ! EXTRADIAGONAL TERMS
200 !
201  a12(ielem) = - l12
202  a13(ielem) = - l13
203  a23(ielem) = - l23
204  a21(ielem) = - l21
205  a31(ielem) = - l31
206  a32(ielem) = - l32
207 !
208  ENDDO ! IELEM
209 !
210  ELSEIF(ielmu.EQ.11.AND.ielmv.EQ.11) THEN
211 !
212 ! TRADITIONAL SCHEME, U AND V LINEAR
213 !
214 ! LOOP ON THE ELEMENTS
215 !
216  DO ielem = 1 , nelem
217 !
218  x2 = xel(ielem,2) * sur24
219  x3 = xel(ielem,3) * sur24
220  y2 = yel(ielem,2) * sur24
221  y3 = yel(ielem,3) * sur24
222 !
223  u1 = u(ikle(ielem,1))
224  u2 = u(ikle(ielem,2))
225  u3 = u(ikle(ielem,3))
226  v1 = v(ikle(ielem,1))
227  v2 = v(ikle(ielem,2))
228  v3 = v(ikle(ielem,3))
229 !
230  u123 = u1 + u2 + u3
231  v123 = v1 + v2 + v3
232 !
233 ! DIAGONAL TERMS
234 !
235  a11(ielem) = (x3-x2) * (v123+v1) + (y2-y3) * (u123+u1)
236  a22(ielem) = -x3 * (v123+v2) +y3 * (u123+u2)
237  a33(ielem) = x2 * (v123+v3) - y2 * (u123+u3)
238 !
239 ! EXTRADIAGONAL TERMS
240 !
241  a12(ielem) = -x3 * (v123+v1) + y3 * (u123+u1)
242  a13(ielem) = x2 * (v123+v1) - y2 * (u123+u1)
243  a23(ielem) = x2 * (v123+v2) - y2 * (u123+u2)
244  a21(ielem) = (x3-x2) * (v123+v2) + (y2-y3) * (u123+u2)
245  a31(ielem) = (x3-x2) * (v123+v3) + (y2-y3) * (u123+u3)
246  a32(ielem) = -x3 * (v123+v3) + y3 * (u123+u3)
247 !
248  ENDDO ! IELEM
249 !
250 !-----------------------------------------------------------------------
251 !
252  ELSEIF(ielmu.EQ.12.AND.ielmv.EQ.12) THEN
253 !
254 ! TRADITIONAL SCHEME, U AND V QUASI-BUBBLE
255 !
256 ! LOOP ON THE ELEMENTS
257 !
258  DO ielem = 1 , nelem
259 !
260  x2 = xel(ielem,2)
261  x3 = xel(ielem,3)
262  y2 = yel(ielem,2)
263  y3 = yel(ielem,3)
264 !
265  u1 = u(ikle(ielem,1))
266  u2 = u(ikle(ielem,2))
267  u3 = u(ikle(ielem,3))
268  u4 = u(ikle(ielem,4))
269  v1 = v(ikle(ielem,1))
270  v2 = v(ikle(ielem,2))
271  v3 = v(ikle(ielem,3))
272  v4 = v(ikle(ielem,4))
273 !
274  cov1 = 5*v3+12*v4+5*v2+14*v1
275  cou1 = 5*u3+12*u4+5*u2+14*u1
276  cov2 = 5*v3+12*v4+14*v2+5*v1
277  cou2 = 5*u3+12*u4+14*u2+5*u1
278  cov3 = 14*v3+12*v4+5*v2+5*v1
279  cou3 = 14*u3+12*u4+5*u2+5*u1
280 !
281 ! EXTRADIAGONAL TERMS
282 !
283  a12(ielem) = ( -x3*cov1 + y3*cou1 )*sur216
284  a13(ielem) = ( x2*cov1 - y2*cou1 )*sur216
285  a21(ielem) = ( (x3-x2)*cov2 + (y2-y3)*cou2 )*sur216
286  a23(ielem) = ( x2*cov2 - y2*cou2 )*sur216
287  a31(ielem) = ( (x3-x2)*cov3 + (y2-y3)*cou3 )*sur216
288  a32(ielem) = ( -x3*cov3 + y3*cou3 )*sur216
289 !
290 ! DIAGONAL TERMS (SUM OF EACH COLUMN = 0)
291 !
292  a11(ielem) = - a12(ielem) - a13(ielem)
293  a22(ielem) = - a21(ielem) - a23(ielem)
294  a33(ielem) = - a31(ielem) - a32(ielem)
295 !
296  ENDDO ! IELEM
297 !
298 !-----------------------------------------------------------------------
299 !
300  ELSEIF(ielmu.EQ.13.AND.ielmv.EQ.13) THEN
301 !
302 ! TRADITIONAL SCHEME, U AND V P2
303 !
304 ! LOOP ON THE ELEMENTS
305 !
306  DO ielem = 1 , nelem
307 !
308  x2 = xel(ielem,2)
309  x3 = xel(ielem,3)
310  y2 = yel(ielem,2)
311  y3 = yel(ielem,3)
312 !
313  u1 = u(ikle(ielem,1))
314  u2 = u(ikle(ielem,2))
315  u3 = u(ikle(ielem,3))
316  u4 = u(ikle(ielem,4))
317  u5 = u(ikle(ielem,5))
318  u6 = u(ikle(ielem,6))
319  v1 = v(ikle(ielem,1))
320  v2 = v(ikle(ielem,2))
321  v3 = v(ikle(ielem,3))
322  v4 = v(ikle(ielem,4))
323  v5 = v(ikle(ielem,5))
324  v6 = v(ikle(ielem,6))
325 !
326 ! EXTRADIAGONAL TERMS
327 !
328  a12(ielem) = ((v3+v2-2.d0*v1-4.d0*v5-8.d0*(v4+v6)) * x3
329  & - (u2-4.d0*u5-8.d0*(u4+u6)-2.d0*u1+u3) * y3) * sur120
330  a13(ielem) = ((2.d0*v1-v3-v2+4.d0*v5+8.d0*(v4+v6)) * x2
331  & + (u2-4.d0*u5-8.d0*(u4+u6)-2.d0*u1+u3) * y2) * sur120
332  a21(ielem) = ((v1-8.d0*(v4+v5)-4.d0*v6+v3-2.d0*v2) * (x2-x3)
333  & + (2.d0*u2-u3-u1+8.d0*(u5+u4)+4.d0*u6) * (y2-y3))
334  & * sur120
335  a23(ielem) = ((8.d0*(v4+v5)-v1+4.d0*v6-v3+2.d0*v2) * x2
336  & - (2.d0*u2-u3-u1+8.d0*(u5+u4)+4.d0*u6) * y2) * sur120
337  a31(ielem) = ((v1+v2-4.d0*v4-2.d0*v3-8.d0*(v6+v5)) * (x2-x3)
338  & + (4.d0*u4-u2+8.d0*(u6+u5)+2.d0*u3-u1) * (y2-y3))
339  & * sur120
340  a32(ielem) = ((v1+v2-4.d0*v4-2.d0*v3-8.d0*(v6+v5)) * x3
341  & + (4.d0*u4-u2+8.d0*(u6+u5)+2.d0*u3-u1) * y3) * sur120
342 !
343 ! DIAGONAL TERMS (SUM OF EACH COLUMN = 0)
344 !
345  a11(ielem) = - a12(ielem) - a13(ielem)
346  a22(ielem) = - a21(ielem) - a23(ielem)
347  a33(ielem) = - a31(ielem) - a32(ielem)
348 !
349  ENDDO ! IELEM
350 !
351 ! OTHER TYPES OF U AND V DISCRETISATION
352 !
353 !-----------------------------------------------------------------------
354 !
355  ELSE
356 !
357  IF(ielmu.EQ.ielmv) THEN
358  WRITE(lu,101) ielmu
359 101 FORMAT(1x,'MT05AA (BIEF) :',/,
360  & 1x,'DISCRETIZATION OF U AND V : ',1i6,' NOT AVAILABLE')
361  ELSE
362  WRITE(lu,201) ielmu,ielmv
363 201 FORMAT(1x,'MT05AA (BIEF) :',/,
364  & 1x,'U AND V OF A DIFFERENT DISCRETISATION:',1i6,3x,1i6)
365  ENDIF
366 !
367  CALL plante(1)
368  stop
369 !
370  ENDIF
371 !
372 !-----------------------------------------------------------------------
373 !
374  RETURN
375  END
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
mt05aa
subroutine mt05aa(A11, A12, A13, A21, A22, A23, A31, A32, A33, XMUL, SU, SV, U, V, XEL, YEL, IKLE, NELEM, NELMAX, FORMUL)
Definition: mt05aa.f:11
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:06 for The TELEMAC-MASCARET system by doxygen 1.8.13