The TELEMAC-MASCARET system  trunk
mt12ab.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt12ab
3 ! *****************
4 !
5  &( a11 , a12 , a13 , a14 ,
6  & a21 , a22 , a23 , a24 ,
7  & a31 , a32 , a33 , a34 ,
8  & xmul,sf,su,sv,f,u,v,
9  & xel,yel,surfac,
10  & ikle1,ikle2,ikle3,ikle4,
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 QUASI-BUBBLE
28 !+ F : FUNCTION OF TYPE P1 TRIANGLE
29 !+ U : VECTOR OF TYPE P0 OR QUASI-BUBBLE
30 !+
31 !+ IT WOULD BE A DERIVATIVE WRT Y WITH ICOORD=2
32 !
33 !warning THE JACOBIAN MUST BE POSITIVE
34 !
35 !history J-M HERVOUET (LNH) ; C MOULIN (LNH)
36 !+ 09/12/94
37 !+ V5P1
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 !| A12 |<--| ELEMENTS OF MATRIX
55 !| A13 |<--| ELEMENTS OF MATRIX
56 !| A14 |<--| ELEMENTS OF MATRIX
57 !| A21 |<--| ELEMENTS OF MATRIX
58 !| A22 |<--| ELEMENTS OF MATRIX
59 !| A23 |<--| ELEMENTS OF MATRIX
60 !| A24 |<--| ELEMENTS OF MATRIX
61 !| A31 |<--| ELEMENTS OF MATRIX
62 !| A32 |<--| ELEMENTS OF MATRIX
63 !| A33 |<--| ELEMENTS OF MATRIX
64 !| A34 |<--| ELEMENTS OF MATRIX
65 !| F |-->| FUNCTION USED IN THE FORMULA
66 !| ICOORD |-->| 1: DERIVATIVE ALONG X, 2: ALONG Y
67 !| IKLE1 |-->| FIRST POINTS OF TRIANGLES
68 !| IKLE2 |-->| SECOND POINTS OF TRIANGLES
69 !| IKLE3 |-->| THIRD POINTS OF TRIANGLES
70 !| IKLE4 |-->| QUASI-BUBBLE POINT
71 !| NELEM |-->| NUMBER OF ELEMENTS
72 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
73 !| SF |-->| STRUCTURE OF FUNCTIONS F
74 !| SU |-->| BIEF_OBJ STRUCTURE OF U
75 !| SURFAC |-->| AREA OF TRIANGLES
76 !| SV |-->| BIEF_OBJ STRUCTURE OF V
77 !| U |-->| FUNCTION U USED IN THE FORMULA
78 !| V |-->| FUNCTION V USED IN THE FORMULA
79 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
80 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
81 !| XMUL |-->| MULTIPLICATION FACTOR
82 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
83 !
84  USE bief, ex_mt12ab => mt12ab
85 !
86  USE declarations_special
87  IMPLICIT NONE
88 !
89 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
90 !
91  INTEGER, INTENT(IN) :: NELEM,NELMAX,ICOORD
92  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax)
93  INTEGER, INTENT(IN) :: IKLE3(nelmax),IKLE4(nelmax)
94 !
95  DOUBLE PRECISION, INTENT(INOUT) :: A11(*),A12(*),A13(*),A14(*)
96  DOUBLE PRECISION, INTENT(INOUT) :: A21(*),A22(*),A23(*),A24(*)
97  DOUBLE PRECISION, INTENT(INOUT) :: A31(*),A32(*),A33(*),A34(*)
98 !
99  DOUBLE PRECISION, INTENT(IN) :: XMUL
100  DOUBLE PRECISION, INTENT(IN) :: F(*),U(*),V(*)
101 !
102 ! STRUCTURES OF F, U, V
103  TYPE(bief_obj), INTENT(IN) :: SF,SU,SV
104 !
105  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,3),YEL(nelmax,3)
106  DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
107 !
108 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
109 !
110  INTEGER IELEM,IELMF,IELMU,IELMV
111  DOUBLE PRECISION X2,X3,Y2,Y3,F1,F2,F3
112  DOUBLE PRECISION U1,U2,U3,U4,V1,V2,V3,V4,UX,UY
113  DOUBLE PRECISION XSUR18,XSUR12,XSUR72,XSU144
114  DOUBLE PRECISION AUX18,AUX12,AUX144,AUX72,UNSURF
115 !
116 !-----------------------------------------------------------------------
117 !
118  ielmf=sf%ELM
119  ielmu=su%ELM
120  ielmv=sv%ELM
121 !
122  xsur18 = xmul / 18.d0
123  xsur12 = xmul / 12.d0
124  xsur72 = xmul / 72.d0
125  xsu144 = xmul / 144.d0
126 !
127 !-----------------------------------------------------------------------
128 ! CASE WHERE F IS OF TYPE P1
129 !-----------------------------------------------------------------------
130 !
131  IF(ielmf.EQ.11) THEN
132 !
133  IF(ielmu.EQ.10.AND.ielmv.EQ.10) THEN
134 !
135 !================================
136 ! DERIVATIVE WRT X =
137 !================================
138 !
139  IF(icoord.EQ.1) THEN
140 !
141 ! LOOP ON THE ELEMENTS
142 !
143  DO ielem = 1 , nelem
144 !
145 ! INITIALISES THE GEOMETRICAL VARIABLES
146 !
147  x2 = xel(ielem,2)
148  x3 = xel(ielem,3)
149  y2 = yel(ielem,2)
150  y3 = yel(ielem,3)
151 !
152  f1 = f(ikle1(ielem))
153  f2 = f(ikle2(ielem)) - f1
154  f3 = f(ikle3(ielem)) - f1
155 !
156  ux = u(ielem)
157  uy = v(ielem)
158 !
159  unsurf= 1.d0 / surfac(ielem)
160  aux12 = xsur12 * unsurf
161  aux18 = xsur18 * unsurf
162 !
163 ! EXTRADIAGONAL TERMS
164 !
165  a12(ielem)=-( (y3-y2)*ux+x2*uy-x3*uy )*(y3*f2-y2*f3)*aux18
166  a14(ielem)=-( (y3-y2)*ux+x2*uy-x3*uy )*(y3*f2-y2*f3)*aux12
167  a21(ielem)=-(x3*uy-ux*y3)*(y3*f2-y2*f3)*aux18
168  a24(ielem)=-(x3*uy-ux*y3)*(y3*f2-y2*f3)*aux12
169  a31(ielem)= (x2*uy-ux*y2)*(y3*f2-y2*f3)*aux18
170  a34(ielem)= (x2*uy-ux*y2)*(y3*f2-y2*f3)*aux12
171  a13(ielem)= a12(ielem)
172  a23(ielem)= a21(ielem)
173  a32(ielem)= a31(ielem)
174 !
175 ! DIAGONAL TERMS
176 ! THE SUM OF EACH COLUMN IS 0
177 !
178  a11(ielem) = - a21(ielem) - a31(ielem)
179  a22(ielem) = - a12(ielem) - a32(ielem)
180  a33(ielem) = - a13(ielem) - a23(ielem)
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  aux18 = xsur18 * unsurf
209 !
210 ! EXTRADIAGONAL TERMS
211 !
212  a12(ielem)=-(x2*uy-x3*uy+ux*y3-ux*y2)*(x2*f3-x3*f2)*aux18
213  a14(ielem)=-(x2*uy-x3*uy+ux*y3-ux*y2)*(x2*f3-x3*f2)*aux12
214  a21(ielem)=-(x2*f3-x3*f2)*(x3*uy-ux*y3)*aux18
215  a24(ielem)=-(x2*f3-x3*f2)*(x3*uy-ux*y3)*aux12
216  a31(ielem)= (x2*uy-ux*y2)*(x2*f3-x3*f2)*aux18
217  a34(ielem)= (x2*uy-ux*y2)*(x2*f3-x3*f2)*aux12
218  a13(ielem)= a12(ielem)
219  a23(ielem)= a21(ielem)
220  a32(ielem)= a31(ielem)
221 !
222 ! DIAGONAL TERMS
223 ! THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
224 !
225  a11(ielem) = - a21(ielem) - a31(ielem)
226  a22(ielem) = - a12(ielem) - a32(ielem)
227  a33(ielem) = - a13(ielem) - a23(ielem)
228 !
229  ENDDO ! IELEM
230 !
231  ELSE
232 !
233  WRITE(lu,201) icoord
234  CALL plante(0)
235  stop
236  ENDIF
237 !
238 !
239  ELSEIF(ielmu.EQ.12) THEN
240 !
241 !================================
242 ! DERIVATIVE WRT X =
243 !================================
244 !
245  IF(icoord.EQ.1) THEN
246 !
247 ! LOOP ON THE ELEMENTS
248 !
249  DO ielem = 1 , nelem
250 !
251 ! INITIALISES THE GEOMETRICAL VARIABLES
252 !
253  x2 = xel(ielem,2)
254  x3 = xel(ielem,3)
255  y2 = yel(ielem,2)
256  y3 = yel(ielem,3)
257 !
258  f1 = f(ikle1(ielem))
259  f2 = f(ikle2(ielem)) - f1
260  f3 = f(ikle3(ielem)) - f1
261 !
262  u1 = u(ikle1(ielem))
263  u2 = u(ikle2(ielem))
264  u3 = u(ikle3(ielem))
265  u4 = u(ikle4(ielem))
266  v1 = v(ikle1(ielem))
267  v2 = v(ikle2(ielem))
268  v3 = v(ikle3(ielem))
269  v4 = v(ikle4(ielem))
270 !
271  unsurf= 1.d0 / surfac(ielem)
272  aux144= xsu144 * unsurf
273  aux72 = xsur72 * unsurf
274 !
275 ! EXTRADIAGONAL TERMS
276 !
277  a12(ielem)=(-((v3+v4+2*v2)*x2-(v3+v4+2*v2)*x3+(v4+2*v2+
278  & v1)*x2-(v4+2*v2+v1)*x3+(y3-y2)*u3+2*(y3-y2)*u4+4*(y3-
279  & y2)*u2+(y3-y2)*u1)*(y3*f2-y2*f3))*aux144
280  a13(ielem)=(-((2*v3+v4+v2)*x2-(2*v3+v4+v2)*x3+(2*v3+v4+
281  & v1)*x2-(2*v3+v4+v1)*x3+4*(y3-y2)*u3+2*(y3-y2)*u4+(y3-
282  & y2)*u2+(y3-y2)*u1)*(y3*f2-y2*f3))*aux144
283  a14(ielem)=(-(x2*v3+3*x2*v4+x2*v2+x2*v1-x3*v3-3*x3*v4-x3
284  & *v2-x3*v1+u3*y3-u3*y2+3*u4*y3-3*u4*y2+u2*y3-u2*y2+u1*y3
285  & -u1*y2)*(y3*f2-y2*f3))*aux72
286  a21(ielem)=(-(x3*v3+2*x3*v4+x3*v2+4*x3*v1-u3*y3-2*u4*y3
287  & -u2*y3-4*u1*y3)*(y3*f2-y2*f3))*aux144
288  a23(ielem)=(-(4*x3*v3+2*x3*v4+x3*v2+x3*v1-4*u3*y3-2*u4
289  & *y3-u2*y3-u1*y3)*(y3*f2-y2*f3))*aux144
290  a24(ielem)=(-(x3*v3+3*x3*v4+x3*v2+x3*v1-u3*y3-3*u4*y3-u2
291  & *y3-u1*y3)*(y3*f2-y2*f3))*aux72
292  a31(ielem)=((x2*v3+2*x2*v4+x2*v2+4*x2*v1-u3*y2-2*u4*y2-
293  & u2*y2-4*u1*y2)*(y3*f2-y2*f3))*aux144
294  a32(ielem)=((x2*v3+2*x2*v4+4*x2*v2+x2*v1-u3*y2-2*u4*y2-
295  & 4*u2*y2-u1*y2)*(y3*f2-y2*f3))*aux144
296  a34(ielem)=((x2*v3+3*x2*v4+x2*v2+x2*v1-u3*y2-3*u4*y2-u2*
297  & y2-u1*y2)*(y3*f2-y2*f3))*aux72
298 !
299 ! DIAGONAL TERMS
300 ! THE SUM OF EACH COLUMN IS 0
301 !
302  a11(ielem) = - a21(ielem) - a31(ielem)
303  a22(ielem) = - a12(ielem) - a32(ielem)
304  a33(ielem) = - a13(ielem) - a23(ielem)
305 !
306  ENDDO ! IELEM
307 !
308  ELSEIF(icoord.EQ.2) THEN
309 !
310 !================================
311 ! DERIVATIVE WRT Y =
312 !================================
313 !
314  DO ielem = 1 , nelem
315 !
316 ! INITIALISES THE GEOMETRICAL VARIABLES
317 !
318  x2 = xel(ielem,2)
319  x3 = xel(ielem,3)
320  y2 = yel(ielem,2)
321  y3 = yel(ielem,3)
322 !
323  f1 = f(ikle1(ielem))
324  f2 = f(ikle2(ielem)) - f1
325  f3 = f(ikle3(ielem)) - f1
326 !
327  u1 = u(ikle1(ielem))
328  u2 = u(ikle2(ielem))
329  u3 = u(ikle3(ielem))
330  u4 = u(ikle4(ielem))
331  v1 = v(ikle1(ielem))
332  v2 = v(ikle2(ielem))
333  v3 = v(ikle3(ielem))
334  v4 = v(ikle4(ielem))
335 !
336  unsurf= 1.d0 / surfac(ielem)
337  aux144= xsu144 * unsurf
338  aux72= xsur72 * unsurf
339 !
340 ! EXTRADIAGONAL TERMS
341 !
342  a12(ielem)=(-(x2*v3+2*x2*v4+4*x2*v2+x2*v1-x3*v3-2*x3*v4
343  & -4*x3*v2-x3*v1+u3*y3-u3*y2+2*u4*y3-2*u4*y2+4*u2*y3-4
344  & *u2*y2+u1*y3-u1*y2)*(x2*f3-x3*f2))*aux144
345  a13(ielem)=(-(4*x2*v3+2*x2*v4+x2*v2+x2*v1-4*x3*v3-2*x3
346  & *v4-x3*v2-x3*v1+4*u3*y3-4*u3*y2+2*u4*y3-2*u4*y2+u2*y3
347  & -u2*y2+u1*y3-u1*y2)*(x2*f3-x3*f2))*aux144
348  a14(ielem)=(-(x2*v3+3*x2*v4+x2*v2+x2*v1-x3*v3-3*x3*v4-x3
349  & *v2-x3*v1+u3*y3-u3*y2+3*u4*y3-3*u4*y2+u2*y3-u2*y2+u1*y3
350  & -u1*y2)*(x2*f3-x3*f2))*aux72
351  a21(ielem)=(-(x2*f3-x3*f2)*(x3*v3+2*x3*v4+x3*v2+4*x3*v1-
352  & u3*y3-2*u4*y3-u2*y3-4*u1*y3))*aux144
353  a23(ielem)=(-(x2*f3-x3*f2)*(4*x3*v3+2*x3*v4+x3*v2+x3*v1-
354  & 4*u3*y3-2*u4*y3-u2*y3-u1*y3))*aux144
355  a24(ielem)=(-(x2*f3-x3*f2)*(x3*v3+3*x3*v4+x3*v2+x3*v1-u3*
356  & y3-3*u4*y3-u2*y3-u1*y3))*aux72
357  a31(ielem)=((x2*v3+2*x2*v4+x2*v2+4*x2*v1-u3*y2-2*u4*y2-
358  & u2*y2-4*u1*y2)*(x2*f3-x3*f2))*aux144
359  a32(ielem)=((x2*v3+2*x2*v4+4*x2*v2+x2*v1-u3*y2-2*u4*y2-
360  & 4*u2*y2-u1*y2)*(x2*f3-x3*f2))*aux144
361  a34(ielem)=((x2*v3+3*x2*v4+x2*v2+x2*v1-u3*y2-3*u4*y2-u2*
362  & y2-u1*y2)*(x2*f3-x3*f2))*aux72
363 !
364 ! DIAGONAL TERMS
365 ! THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
366 !
367  a11(ielem) = - a21(ielem) - a31(ielem)
368  a22(ielem) = - a12(ielem) - a32(ielem)
369  a33(ielem) = - a13(ielem) - a23(ielem)
370 !
371  ENDDO ! IELEM
372 !
373  ELSE
374 !
375  WRITE(lu,201) icoord
376  CALL plante(0)
377  stop
378  ENDIF
379 !
380  ELSE
381 !
382  WRITE(lu,301) ielmu
383  CALL plante(0)
384  stop
385 !
386  ENDIF
387 !
388 !-----------------------------------------------------------------------
389 !
390  ELSE
391  WRITE(lu,101) ielmf
392 101 FORMAT(1x,'MT12AB (BIEF) :',/,
393  & 1x,'DISCRETIZATION OF F : ',1i6,' NOT AVAILABLE')
394  CALL plante(0)
395  stop
396  ENDIF
397 !
398 201 FORMAT(1x,'MT12AB (BIEF) : IMPOSSIBLE COMPONENT ',
399  & 1i6,' CHECK ICOORD')
400 !
401 301 FORMAT(1x,'MT12AB (BIEF) :',/,
402  & 1x,'DISCRETIZATION OF U : ',1i6,' NOT AVAILABLE')
403 !
404 !-----------------------------------------------------------------------
405 !
406  RETURN
407  END
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
mt12ab
subroutine mt12ab(A11, A12, A13, A14, A21, A22, A23, A24, A31, A32, A33, A34, XMUL, SF, SU, SV, F, U, V, XEL, YEL, SURFAC, IKLE1, IKLE2, IKLE3, IKLE4, NELEM, NELMAX, ICOORD)
Definition: mt12ab.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