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