The TELEMAC-MASCARET system  trunk
mt12ac.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt12ac
3 ! *****************
4 !
5  &( a11 , a12 , a13 , a14 , a15, a16,
6  & a21 , a22 , a23 , a24 , a25, a26,
7  & a31 , a32 , a33 , a34 , a35, a36,
8  & xmul,sf,su,sv,f,u,v,
9  & xel,yel,surfac,
10  & ikle1,ikle2,ikle3,ikle4,ikle5,ikle6,
11  & nelem,nelmax,icoord)
12 !
13 !***********************************************************************
14 ! BIEF V6P1 21/08/2010
15 !***********************************************************************
16 !
17 !brief COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
18 !code
19 !+ EXAMPLE WITH ICOORD=1
20 !+
21 !+ / D -> -->
22 !+ A(I,J)= XMUL / PSI2(J) * --(F) * U.GRAD(PSI1(I)) D(OMEGA)
23 !+ /OMEGA DX
24 !+
25 !+
26 !+ PSI1 : BASES OF TYPE P1 TRIANGLE
27 !+ PSI2 : BASES OF TYPE P2
28 !+ F : FUNCTION OF TYPE P1 TRIANGLE
29 !+ U : VECTOR OF TYPE P0 OR P2
30 !+
31 !+ IT WOULD BE A DERIVATIVE WRT Y WITH ICOORD=2
32 !
33 !warning THE JACOBIAN MUST BE POSITIVE
34 !
35 !history ALGIANE FROEHLY (MATMECA)
36 !+ 19/06/08
37 !+ V5P9
38 !+
39 !
40 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
41 !+ 13/07/2010
42 !+ V6P0
43 !+ Translation of French comments within the FORTRAN sources into
44 !+ English comments
45 !
46 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
47 !+ 21/08/2010
48 !+ V6P0
49 !+ Creation of DOXYGEN tags for automated documentation and
50 !+ cross-referencing of the FORTRAN sources
51 !
52 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
53 !| A11 |<--| ELEMENTS OF MATRIX
54 !| ... |<--| ELEMENTS OF MATRIX
55 !| A36 |<--| ELEMENTS OF MATRIX
56 !| F |-->| FUNCTION USED IN THE FORMULA
57 !| ICOORD |-->| 1: DERIVATIVE ALONG X, 2: ALONG Y
58 !| IKLE1 |-->| FIRST POINTS OF TRIANGLES
59 !| IKLE2 |-->| SECOND POINTS OF TRIANGLES
60 !| IKLE3 |-->| THIRD POINTS OF TRIANGLES
61 !| IKLE4 |-->| QUASI-BUBBLE POINT
62 !| NELEM |-->| NUMBER OF ELEMENTS
63 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
64 !| SF |-->| STRUCTURE OF FUNCTIONS F
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_MT12AC => MT12AC
76 !
77  USE declarations_special
78  IMPLICIT NONE
79 !
80 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
81 !
82  INTEGER, INTENT(IN) :: NELEM,NELMAX,ICOORD
83  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax)
84  INTEGER, INTENT(IN) :: IKLE3(nelmax),IKLE4(nelmax)
85  INTEGER, INTENT(IN) :: IKLE5(nelmax),IKLE6(nelmax)
86 !
87  DOUBLE PRECISION, INTENT(INOUT) :: A11(*),A12(*),A13(*)
88  DOUBLE PRECISION, INTENT(INOUT) :: A14(*),A15(*),A16(*)
89  DOUBLE PRECISION, INTENT(INOUT) :: A21(*),A22(*),A23(*)
90  DOUBLE PRECISION, INTENT(INOUT) :: A24(*),A25(*),A26(*)
91  DOUBLE PRECISION, INTENT(INOUT) :: A31(*),A32(*),A33(*)
92  DOUBLE PRECISION, INTENT(INOUT) :: A34(*),A35(*),A36(*)
93  DOUBLE PRECISION, INTENT(IN) :: XMUL
94  DOUBLE PRECISION, INTENT(IN) :: F(*),U(*),V(*)
95 !
96 ! STRUCTURES OF F, U, V
97  TYPE(bief_obj), INTENT(IN) :: SF,SU,SV
98 !
99  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,3),YEL(nelmax,3)
100  DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
101 !
102 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
103 !
104  INTEGER IELEM,IELMF,IELMU,IELMV
105  DOUBLE PRECISION X2,X3,Y2,Y3,F1,F2,F3
106  DOUBLE PRECISION U1,U2,U3,U4,U5,U6,V1,V2,V3,V4,V5,V6,UX,UY
107  DOUBLE PRECISION XSUR12 ,XSUR180,XSUR720
108  DOUBLE PRECISION AUX12 ,AUX180 ,AUX720 , UNSURF
109 !
110 !-----------------------------------------------------------------------
111 !
112  ielmf=sf%ELM
113  ielmu=su%ELM
114  ielmv=sv%ELM
115 !
116  xsur12 = xmul / 12.d0
117  xsur720 = xmul / 720.d0
118  xsur180 = xmul / 180.d0
119 !
120 !-----------------------------------------------------------------------
121 ! CASE WHERE F IS OF TYPE P1
122 !-----------------------------------------------------------------------
123 !
124  IF(ielmf.EQ.11) THEN
125 !
126  IF(ielmu.EQ.10.AND.ielmv.EQ.10) THEN
127 !
128 !================================
129 ! DERIVATIVE WRT X =
130 !================================
131 !
132  IF(icoord.EQ.1) THEN
133 !
134 ! LOOP ON THE ELEMENTS
135 !
136  DO ielem = 1 , nelem
137 !
138 ! INITIALISES THE GEOMETRICAL VARIABLES
139 !
140  x2 = xel(ielem,2)
141  x3 = xel(ielem,3)
142  y2 = yel(ielem,2)
143  y3 = yel(ielem,3)
144 !
145  f1 = f(ikle1(ielem))
146  f2 = f(ikle2(ielem)) - f1
147  f3 = f(ikle3(ielem)) - f1
148 !
149  ux = u(ielem)
150  uy = v(ielem)
151 !
152  unsurf= 1.d0 / surfac(ielem)
153  aux12 = xsur12 * unsurf
154 !
155 ! EXTRADIAGONAL TERMS
156 !
157  a12(ielem)= 0.d0
158  a13(ielem)= 0.d0
159  a21(ielem)= 0.d0
160  a23(ielem)= 0.d0
161  a24(ielem)= (y3*ux-uy*x3) * (y3*f2-y2*f3) * aux12
162  a25(ielem)= a24(ielem)
163  a26(ielem)= a24(ielem)
164  a31(ielem)= 0.d0
165  a32(ielem)= 0.d0
166  a34(ielem)= (x2*uy-ux*y2) * (y3*f2-y2*f3) * aux12
167  a35(ielem)= a34(ielem)
168  a36(ielem)= a34(ielem)
169 !
170 ! THE SUM OF EACH COLUMN IS 0
171 !
172  a14(ielem)= - a24(ielem) - a34(ielem)
173  a15(ielem)= a14(ielem)
174  a16(ielem)= a14(ielem)
175 !
176 !C DIAGONAL ELEMENTS
177 !
178  a11(ielem) = 0.d0
179  a22(ielem) = 0.d0
180  a33(ielem) = 0.d0
181 !
182  ENDDO ! IELEM
183 !
184  ELSEIF(icoord.EQ.2) THEN
185 !
186 !================================
187 ! DERIVATIVE WRT Y =
188 !================================
189 !
190  DO ielem = 1 , nelem
191 !
192 ! INITIALISES THE GEOMETRICAL VARIABLES
193 !
194  x2 = xel(ielem,2)
195  x3 = xel(ielem,3)
196  y2 = yel(ielem,2)
197  y3 = yel(ielem,3)
198 !
199  f1 = f(ikle1(ielem))
200  f2 = f(ikle2(ielem)) - f1
201  f3 = f(ikle3(ielem)) - f1
202 !
203  ux = u(ielem)
204  uy = v(ielem)
205 !
206  unsurf = 1.d0 / surfac(ielem)
207  aux12 = xsur12 * unsurf
208 !
209 ! EXTRADIAGONAL TERMS
210 !
211  a12(ielem)= 0.d0
212  a13(ielem)= 0.d0
213  a21(ielem)= 0.d0
214  a23(ielem)= 0.d0
215  a24(ielem)= (x3*f2-x2*f3) * (uy*x3-ux*y3) * aux12
216  a25(ielem)= a24(ielem)
217  a26(ielem)= a24(ielem)
218  a31(ielem)= 0.d0
219  a32(ielem)= 0.d0
220  a34(ielem)= (x3*f2-x2*f3) * (ux*y2-uy*x2) * aux12
221  a35(ielem)= a34(ielem)
222  a36(ielem)= a34(ielem)
223 !
224 ! THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
225 !
226  a14(ielem)= - a24(ielem) - a34(ielem)
227  a15(ielem)= a14(ielem)
228  a16(ielem)= a14(ielem)
229 !
230 !C DIAGONAL TERMS
231 !
232  a11(ielem) = 0.d0
233  a22(ielem) = 0.d0
234  a33(ielem) = 0.d0
235 !
236  ENDDO ! IELEM
237 !
238  ELSE
239 !
240  WRITE(lu,201) icoord
241  CALL plante(1)
242  stop
243  ENDIF
244 !
245  ELSEIF(ielmu.EQ.13) THEN
246 !
247 !================================
248 ! DERIVATIVE WRT X =
249 !================================
250 !
251  IF(icoord.EQ.1) THEN
252 !
253 ! LOOP ON THE ELEMENTS
254 !
255  DO ielem = 1 , nelem
256 !
257 ! INITIALISES THE GEOMETRICAL VARIABLES
258 !
259  x2 = xel(ielem,2)
260  x3 = xel(ielem,3)
261  y2 = yel(ielem,2)
262  y3 = yel(ielem,3)
263 !
264  f1 = f(ikle1(ielem))
265  f2 = f(ikle2(ielem)) - f1
266  f3 = f(ikle3(ielem)) - f1
267 !
268  u1 = u(ikle1(ielem))
269  u2 = u(ikle2(ielem))
270  u3 = u(ikle3(ielem))
271  u4 = u(ikle4(ielem))
272  u5 = u(ikle5(ielem))
273  u6 = u(ikle6(ielem))
274  v1 = v(ikle1(ielem))
275  v2 = v(ikle2(ielem))
276  v3 = v(ikle3(ielem))
277  v4 = v(ikle4(ielem))
278  v5 = v(ikle5(ielem))
279  v6 = v(ikle6(ielem))
280 !
281  unsurf = 1.d0 / surfac(ielem)
282  aux720 = xsur720 * unsurf
283  aux180 = xsur180 * unsurf
284 !
285 ! EXTRADIAGONAL TERMS
286 !
287  a12(ielem)= (y3*f2-y2*f3) * ((u3+u1-6.d0*u2+4.d0*u6) * (y3-y2)
288  & + (v3+v1-6.d0*v2+4.d0*v6) * (x2-x3)) * aux720
289 !
290  a13(ielem)= (y3*f2-y2*f3) * ((4.d0*v4-6.d0*v3+v2+v1) * (x2-x3)
291  & + (6.d0*u3-4.d0*u4-u2-u1) * (y2-y3)) * aux720
292 !
293  a21(ielem)= (y3*f2-y2*f3) * ((4.d0*v5+v3+v2-6.d0*v1) * x3
294  & - (u3-6.d0*u1+4.d0*u5+u2) * y3 ) * aux720
295 !
296  a23(ielem)= (y3*f2-y2*f3) * ((v1-6.d0*v3+4.d0*v4+v2) * x3
297  & + (6.d0*u3-4.d0*u4-u2-u1) * y3 ) * aux720
298 !
299  a24(ielem)= (y3*f2-y2*f3) * ((v3-4.d0*v6-8.d0*v4-4.d0*v5) * x3
300  & - (u3-8.d0*u4-4.d0*u6-4.d0*u5) * y3) * aux180
301 !
302  a25(ielem)= (y3*f2-y2*f3) * ((v1-4.d0*v6-8.d0*v5-4.d0*v4) * x3
303  & + (4.d0*u6-u1+4.d0*u4+8.d0*u5) * y3) * aux180
304 !
305  a26(ielem)= (y3*f2-y2*f3) * ((v2-4.d0*v4-4.d0*v5-8.d0*v6) * x3
306  & + (8.d0*u6+4.d0*u4+4.d0*u5-u2) * y3) * aux180
307 !
308  a31(ielem)= (y3*f2-y2*f3) * ((6.d0*v1-4.d0*v5-v3-v2) * x2
309  & + (u3-6.d0*u1+4.d0*u5+u2) * y2 ) * aux720
310 !
311  a32(ielem)= (y3*f2-y2*f3) * ((-4.d0*v6+6.d0*v2-v1-v3) * x2
312  & + (u3+u1+4.d0*u6-6.d0*u2) * y2 ) * aux720
313 !
314  a34(ielem)= (y3*f2-y2*f3) * ((4.d0*v6+8.d0*v4+4.d0*v5-v3) * x2
315  & + (u3-8.d0*u4-4.d0*u6-4.d0*u5) * y2 ) * aux180
316 !
317  a35(ielem)= (y2*f3-y3*f2) * ((v1-8.d0*v5-4.d0*v4-4.d0*v6) * x2
318  & + (4.d0*u6-u1+4.d0*u4+8.d0*u5) * y2 ) * aux180
319 !
320  a36(ielem)= (y2*f3-y3*f2) * ((v2-4.d0*v4-4.d0*v5-8.d0*v6) * x2
321  & + (8.d0*u6+4.d0*u4+4.d0*u5-u2) * y2 ) * aux180
322 !
323 ! THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
324 !
325  a14(ielem) = - a24(ielem) - a34(ielem)
326  a15(ielem) = - a25(ielem) - a35(ielem)
327  a16(ielem) = - a26(ielem) - a36(ielem)
328 !
329 ! DIAGONAL TERMS
330 ! THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
331 !
332  a11(ielem) = - a21(ielem) - a31(ielem)
333  a22(ielem) = - a12(ielem) - a32(ielem)
334  a33(ielem) = - a13(ielem) - a23(ielem)
335 !
336  ENDDO ! IELEM
337 !
338  ELSEIF(icoord.EQ.2) THEN
339 !
340 !================================
341 ! DERIVATIVE WRT Y =
342 !================================
343 !
344  DO ielem = 1 , nelem
345 !
346 ! INITIALISES THE GEOMETRICAL VARIABLES
347 !
348  x2 = xel(ielem,2)
349  x3 = xel(ielem,3)
350  y2 = yel(ielem,2)
351  y3 = yel(ielem,3)
352 !
353  f1 = f(ikle1(ielem))
354  f2 = f(ikle2(ielem)) - f1
355  f3 = f(ikle3(ielem)) - f1
356 !
357  u1 = u(ikle1(ielem))
358  u2 = u(ikle2(ielem))
359  u3 = u(ikle3(ielem))
360  u4 = u(ikle4(ielem))
361  u5 = u(ikle5(ielem))
362  u6 = u(ikle6(ielem))
363  v1 = v(ikle1(ielem))
364  v2 = v(ikle2(ielem))
365  v3 = v(ikle3(ielem))
366  v4 = v(ikle4(ielem))
367  v5 = v(ikle5(ielem))
368  v6 = v(ikle6(ielem))
369 !
370  unsurf= 1.d0 / surfac(ielem)
371  aux180= xsur180 * unsurf
372  aux720= xsur720 * unsurf
373 !
374 ! EXTRADIAGONAL TERMS
375 !
376  a12(ielem)= (x2*f3-x3*f2) * ((4.d0*v6-6.d0*v2+v3+v1) * (x2-x3)
377  & + (u1+u3+4.d0*u6-6.d0*u2) * (y3-y2)) * aux720
378 !
379  a13(ielem)= (x2*f3-x3*f2) * ((v2-6.d0*v3+4.d0*v4+v1) * (x2-x3)
380  & + (6.d0*u3-u1-u2-4.d0*u4) * (y2-y3)) * aux720
381 !
382  a21(ielem)= (x3*f2-x2*f3) * ((6.d0*v1-v3-v2-4.d0*v5) * x3
383  & + (u3-6.d0*u1+u2+4.d0*u5) * y3 ) * aux720
384 !
385  a23(ielem)= (x2*f3-x3*f2) * ((v2-6.d0*v3+4.d0*v4+v1) * x3
386  & + (6.d0*u3-u1-u2-4.d0*u4) * y3 ) * aux720
387 !
388  a24(ielem)= (x2*f3-x3*f2) * ((v3-4.d0*v6-4.d0*v5-8.d0*v4) * x3
389  & + (4.d0*u5+8.d0*u4-u3+4.d0*u6) * y3 ) * aux180
390 !
391  a25(ielem)= (x2*f3-x3*f2) * ((v1-8.d0*v5-4.d0*v4-4.d0*v6) * x3
392  & + (4.d0*u4-u1+8.d0*u5+4.d0*u6) * y3 ) * aux180
393 !
394  a26(ielem)= (x2*f3-x3*f2) * ((v2-4.d0*v5-4.d0*v4-8.d0*v6) * x3
395  & + (4.d0*u5-u2+8.d0*u6+4.d0*u4) * y3 ) * aux180
396 !
397  a31(ielem)= (x2*f3-x3*f2) * ((6.d0*v1-v3-v2-4.d0*v5) * x2
398  & + (u3-6.d0*u1+u2+4.d0*u5) * y2 ) * aux720
399 !
400  a32(ielem)= (x3*f2-x2*f3) * ((4.d0*v6-6.d0*v2+v3+v1) * x2
401  & + (6.d0*u2-u1-u3-4.d0*u6) * y2 ) * aux720
402 !
403  a34(ielem)= (x3*f2-x2*f3) * ((v3-4.d0*v6-4.d0*v5-8.d0*v4) * x2
404  & + (4.d0*u5+8.d0*u4-u3+4.d0*u6) * y2 ) * aux180
405 !
406  a35(ielem)= (x3*f2-x2*f3) * ((v1-8.d0*v5-4.d0*v4-4.d0*v6) * x2
407  & + (4.d0*u4-u1+8.d0*u5+4.d0*u6) * y2 ) * aux180
408 !
409  a36(ielem)= (x3*f2-x2*f3) * ((v2-4.d0*v5-4.d0*v4-8.d0*v6) * x2
410  & + (4.d0*u5-u2+8.d0*u6+4.d0*u4) * y2 ) * aux180
411 !
412 ! THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
413 !
414  a14(ielem) = - a24(ielem) - a34(ielem)
415  a15(ielem) = - a25(ielem) - a35(ielem)
416  a16(ielem) = - a26(ielem) - a36(ielem)
417 !
418 ! DIAGONAL TERMS
419 ! THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
420 !
421  a11(ielem) = - a21(ielem) - a31(ielem)
422  a22(ielem) = - a12(ielem) - a32(ielem)
423  a33(ielem) = - a13(ielem) - a23(ielem)
424 !
425  ENDDO ! IELEM
426 !
427  ELSE
428  WRITE(lu,201) icoord
429  CALL plante(1)
430  stop
431  ENDIF
432 !
433  ELSE
434 !
435  WRITE(lu,301) ielmu
436  CALL plante(1)
437  stop
438 !
439  ENDIF
440 !
441 !-----------------------------------------------------------------------
442 !
443  ELSE
444  WRITE(lu,101) ielmf
445 101 FORMAT(1x,'MT12AC (BIEF) :',/,
446  & 1x,'DISCRETIZATION OF F : ',1i6,' NOT AVAILABLE')
447  CALL plante(1)
448  stop
449  ENDIF
450 !
451 201 FORMAT(1x,'MT12AC (BIEF) : IMPOSSIBLE COMPONENT ',
452  & 1i6,' CHECK ICOORD')
453 !
454 301 FORMAT(1x,'MT12AC (BIEF) :',/,
455  & 1x,'DISCRETIZATION OF U : ',1i6,' NOT AVAILABLE')
456 !
457 !-----------------------------------------------------------------------
458 !
459  RETURN
460  END
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
mt12ac
subroutine mt12ac(A11, A12, A13, A14, A15, A16, A21, A22, A23, A24, A25, A26, A31, A32, A33, A34, A35, A36, XMUL, SF, SU, SV, F, U, V, XEL, YEL, SURFAC, IKLE1, IKLE2, IKLE3, IKLE4, IKLE5, IKLE6, NELEM, NELMAX, ICOORD)
Definition: mt12ac.f:13
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:12 for The TELEMAC-MASCARET system by doxygen 1.8.13