The TELEMAC-MASCARET system  trunk
mt05cc.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt05cc
3 ! *****************
4 !
5  &( a11 , a12 , a13 , a14 , a15 , a16 ,
6  & a21 , a22 , a23 , a24 , a25 , a26 ,
7  & a31 , a32 , a33 , a34 , a35 , a36 ,
8  & a41 , a42 , a43 , a44 , a45 , a46 ,
9  & a51 , a52 , a53 , a54 , a55 , a56 ,
10  & a61 , a62 , a63 , a64 , a65 , a66 ,
11  & xmul,su,sv,u,v,
12  & xel,yel,ikle1,ikle2,ikle3,ikle4,ikle5,ikle6,
13  & nelem,nelmax,formul)
14 !
15 !***********************************************************************
16 ! BIEF V6P1 21/08/2010
17 !***********************************************************************
18 !
19 !brief BUILDS THE FOLLOWING MATRIX FOR P2 TRIANGLES:
20 !code
21 !+ ->--->
22 !+ U.GRAD
23 !
24 !history ALGIANE FROEHLY; C MOULIN (LNH)
25 !+ 09/07/2008
26 !+ V5P9
27 !+
28 !
29 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
30 !+ 13/07/2010
31 !+ V6P0
32 !+ Translation of French comments within the FORTRAN sources into
33 !+ English comments
34 !
35 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
36 !+ 21/08/2010
37 !+ V6P0
38 !+ Creation of DOXYGEN tags for automated documentation and
39 !+ cross-referencing of the FORTRAN sources
40 !
41 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
42 !| A11 |<--| ELEMENTS OF MATRIX
43 !| ... |<--| ELEMENTS OF MATRIX
44 !| A66 |<--| ELEMENTS OF MATRIX
45 !| FORMUL |-->| FORMULA DESCRIBING THE RESULTING MATRIX
46 !| IKLE1 |-->| FIRST POINTS OF TRIANGLES
47 !| IKLE2 |-->| SECOND POINTS OF TRIANGLES
48 !| IKLE3 |-->| THIRD POINTS OF TRIANGLES
49 !| IKLE4 |-->| FOURTH POINTS OF TRIANGLES (QUADRATIC)
50 !| IKLE5 |-->| FIFTH POINTS OF TRIANGLES (QUADRATIC)
51 !| IKLE6 |-->| SIXTH POINTS OF TRIANGLES (QUADRATIC)
52 !| NELEM |-->| NUMBER OF ELEMENTS
53 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
54 !| SU |-->| BIEF_OBJ STRUCTURE OF U
55 !| SURFAC |-->| AREA OF TRIANGLES
56 !| SV |-->| BIEF_OBJ STRUCTURE OF V
57 !| U |-->| FUNCTION U USED IN THE FORMULA
58 !| V |-->| FUNCTION V USED IN THE FORMULA
59 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
60 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
61 !| XMUL |-->| MULTIPLICATION FACTOR
62 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
63 !
64  USE bief!, EX_MT05CC => MT05CC
65 !
66  USE declarations_special
67  IMPLICIT NONE
68 !
69 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
70 !
71  INTEGER, INTENT(IN) :: NELEM,NELMAX
72  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax)
73  INTEGER, INTENT(IN) :: IKLE3(nelmax),IKLE4(nelmax)
74  INTEGER, INTENT(IN) :: IKLE5(nelmax),IKLE6(nelmax)
75 !
76  DOUBLE PRECISION, INTENT(INOUT) :: A11(*),A12(*),A13(*)
77  DOUBLE PRECISION, INTENT(INOUT) :: A14(*),A15(*),A16(*)
78  DOUBLE PRECISION, INTENT(INOUT) :: A21(*),A22(*),A23(*)
79  DOUBLE PRECISION, INTENT(INOUT) :: A24(*),A25(*),A26(*)
80  DOUBLE PRECISION, INTENT(INOUT) :: A31(*),A32(*),A33(*)
81  DOUBLE PRECISION, INTENT(INOUT) :: A34(*),A35(*),A36(*)
82  DOUBLE PRECISION, INTENT(INOUT) :: A41(*),A42(*),A43(*)
83  DOUBLE PRECISION, INTENT(INOUT) :: A44(*),A45(*),A46(*)
84  DOUBLE PRECISION, INTENT(INOUT) :: A51(*),A52(*),A53(*)
85  DOUBLE PRECISION, INTENT(INOUT) :: A54(*),A55(*),A56(*)
86  DOUBLE PRECISION, INTENT(INOUT) :: A61(*),A62(*),A63(*)
87  DOUBLE PRECISION, INTENT(INOUT) :: A64(*),A65(*),A66(*)
88 !
89  DOUBLE PRECISION, INTENT(IN) :: XMUL
90  DOUBLE PRECISION, INTENT(IN) :: U(*),V(*)
91 !
92 !
93 ! STRUCTURES OF U, V
94  TYPE(bief_obj), INTENT(IN) :: SU,SV
95 !
96  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,3),YEL(nelmax,3)
97 !
98  CHARACTER(LEN=16) :: FORMUL
99 !
100 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
101 !
102 !
103 ! DECLARATIONS SPECIFIC TO THIS SUBROUTINE
104 !
105  INTEGER IELEM,IELMU,IELMV
106 !
107  DOUBLE PRECISION X2,X3
108  DOUBLE PRECISION Y2,Y3
109  DOUBLE PRECISION U1,U2,U3,U4,U5,U6
110  DOUBLE PRECISION V1,V2,V3,V4,V5,V6
111  DOUBLE PRECISION XSU360,XSUR90,XSUR45
112  DOUBLE PRECISION XSUR2520,XSUR630
113 !=======================================================================
114 !
115 ! EXTRACTS THE TYPE OF ELEMENT FOR VELOCITY
116 !
117  ielmu = su%ELM
118  ielmv = sv%ELM
119 !
120 !-----------------------------------------------------------------------
121 !
122  xsu360 = xmul / 360.d0
123  xsur45 = xmul / 45.d0
124  xsur90 = xmul / 90.d0
125  xsur2520 = xmul / 2520.d0
126  xsur630 = xmul / 630.d0
127 !
128 !-----------------------------------------------------------------------
129 !
130 ! CASE WHERE U AND V ARE CONSTANT BY ELEMENT
131 !
132  IF(ielmu.EQ.10.AND.ielmv.EQ.10) THEN
133 !
134 !-----------------------------------------------------------------------
135 !
136 ! P1 DISCRETISATION OF THE VELOCITY:
137 !
138  DO ielem = 1 , nelem
139 !
140 ! INITIALISES THE GEOMETRICAL VARIABLES
141 !
142  x2 = xel(ielem,2)
143  x3 = xel(ielem,3)
144 !
145  y2 = yel(ielem,2)
146  y3 = yel(ielem,3)
147 !
148  u1 = u(ielem)
149  u2 = u(ielem)
150  u3 = u(ielem)
151  v1 = v(ielem)
152  v2 = v(ielem)
153  v3 = v(ielem)
154 !
155 ! EXTRADIAGONAL TERMS
156 !
157  a12(ielem)=(( 5.d0*v2+v3+v1*6.d0)*x3 +
158  & ( -u1*6.d0-u3-u2*5.d0)*y3) * xsu360
159  a13(ielem)=((-v2-v3*5.d0-v1*6.d0)*x2 +
160  & ( u1*6.d0+u3*5.d0+u2)*y2) * xsu360
161  a14(ielem)=(( v2*2.d0+v3)*x2+(-v2*2.d0-v3-v1*6.d0)*x3 +
162  & ( -u2*2.d0-u3)*y2+(u1*6.d0+u3+u2*2.d0)*y3) * xsur90
163  a15(ielem)=((-v3-v2*2.d0)*x2+(v2+v3*2.d0)*x3 +
164  & (u2*2.d0+u3)*y2+(-u2-u3*2.d0)*y3) * xsur90
165  a16(ielem)=((v2+v3*2.d0+v1*6.d0)*x2+(-v2-v3*2.d0)*x3 +
166  & (-u2-u3*2.d0-u1*6.d0)*y2+(u2+u3*2.d0)*y3) * xsur90
167 !
168  a21(ielem)=((v2*6.d0+v3+v1*5.d0)*x2 +
169  & (-v2*6.d0-v3-v1*5.d0)*x3 +
170  & (-u2*6.d0-u3-u1*5.d0)*y2 +
171  & (u1*5.d0+u3+u2*6.d0)*y3) * xsu360
172  a23(ielem)=((-v2*6.d0-v3*5.d0-v1)*x2 +
173  & ( u1+u3*5.d0+u2*6.d0)*y2) * xsu360
174  a24(ielem)=(-v2*x2*6.d0+(v2*6.d0+v3+v1*2.d0)*x3+u2*y2*6.d0 +
175  & (-u2*6.d0-u3-u1*2.d0)*y3) * xsur90
176  a25(ielem)=(v2*x2*6.d0+(v3*2.d0+v1)*x3-u2*y2*6.d0 +
177  & (-u3*2.d0-u1)*y3) * xsur90
178  a26(ielem)=((v3-v1)*x2+(-v3*2.d0-v1)*x3+(-u3+u1)*y2 +
179  & (u1+u3*2.d0)*y3) * xsur90
180 !
181  a31(ielem)=((v2+v3*6.d0+v1*5.d0)*x2 +
182  & (-v2-v3*6.d0-v1*5.d0)*x3 +
183  & (-u2-u3*6.d0-u1*5.d0)*y2 +
184  & (u1*5.d0+u3*6.d0+u2)*y3) * xsu360
185  a32(ielem)=((v2*5.d0+v3*6.d0+v1)*x3 +
186  & (-u1-u3*6.d0-u2*5.d0)*y3) * xsu360
187  a34(ielem)=((v2*2.d0+v1)*x2+(-v2+v1)*x3+(-u1-u2*2.d0)*y2+
188  & (u2-u1)*y3)*xsur90
189  a35(ielem)=((-v2*2.d0-v1)*x2-v3*x3*6.d0 +
190  & ( u2*2.d0+u1)*y2+u3*y3*6.d0) * xsur90
191  a36(ielem)=((-v3*6.d0-v2-v1*2.d0)*x2+v3*x3*6.d0 +
192  & ( u2+u3*6.d0+u1*2.d0)*y2-u3*y3*6.d0) * xsur90
193 !
194  a41(ielem)=((-v2*2.d0-v3-v1*6.d0)*x2 +
195  & ( v2*2.d0+v3+v1*6.d0)*x3 +
196  & ( u1*6.d0+u3+u2*2.d0)*y2 +
197  & (-u1*6.d0-u3-u2*2.d0)*y3) * xsur90
198  a42(ielem)=((-v2*6.d0-v3-v1*2.d0)*x3 +
199  & ( u1*2.d0+u3+u2*6.d0)*y3) * xsur90
200  a43(ielem)=((-v2*2.d0+v3-v1*2.d0)*x2 +
201  & ( u1*2.d0-u3+u2*2.d0)*y2) * xsur90
202  a45(ielem)=(( 6.d0*v2+2.d0*v3+4.d0*v1)*x2 +
203  & (-2.d0*v2-2.d0*v3-2.d0*v1)*x3 +
204  & (-6.d0*u2-2.d0*u3-4.d0*u1)*y2 +
205  & (2.d0*u1+2.d0*u3+2.d0*u2)*y3) * xsur45
206  a46(ielem)=((2.d0*v2+4.d0*v1)*x2 +
207  & (2.d0*v3+2.d0*v2+2.d0*v1)*x3 +
208  & (-4.d0*u1-2.d0*u2)*y2 +
209  & (-2.d0*u2-2.d0*u3-2.d0*u1)*y3) * xsur45
210 !
211  a51(ielem)=((v2*2.d0+v3*2.d0-v1)*x2 +
212  & (-v2*2.d0-v3*2.d0+v1)*x3 +
213  & (-u2*2.d0-u3*2.d0+u1)*y2 +
214  & (-u1+u3*2.d0+u2*2.d0)*y3) * xsur90
215  a52(ielem)=((-v2*6.d0-v3*2.d0-v1)*x3 +
216  & (u1+u3*2.d0+u2*6.d0)*y3) * xsur90
217  a53(ielem)=((v2*2.d0+v3*6.d0+v1)*x2 +
218  & (-u1-u3*6.d0-u2*2.d0)*y2) * xsur90
219  a54(ielem)=((-6.d0*v2-4.d0*v3-2.d0*v1)*x2 +
220  & (4.d0*v2+2.d0*v3)*x3 +
221  & (6.d0*u2+4.d0*u3+2.d0*u1)*y2 +
222  & (-4.d0*u2-2.d0*u3)*y3) * xsur45
223  a56(ielem)=((-2.d0*v2-4.d0*v3)*x2 +
224  & (4.d0*v2+6.d0*v3+2.d0*v1)*x3 +
225  & (2.d0*u2+4.d0*u3)*y2 +
226  & (-2.d0*u1-6.d0*u3-4.d0*u2)*y3) * xsur45
227 !
228  a61(ielem)=((-v2-v3*2.d0-v1*6.d0)*x2 +
229  & (v2+v3*2.d0+v1*6.d0)*x3 +
230  & (u2+u3*2.d0+u1*6.d0)*y2 +
231  & (-u2-u3*2.d0-u1*6.d0)*y3) * xsur90
232  a62(ielem)=((-v2+v3*2.d0+v1*2.d0)*x3 +
233  & (-u1*2.d0-u3*2.d0+u2)*y3) * xsur90
234  a63(ielem)=((v2+v3*6.d0+v1*2.d0)*x2 +
235  & (-u1*2.d0-u3*6.d0-u2)*y2) * xsur90
236  a64(ielem)=((-2.d0*v2-2.d0*v3-2.d0*v1)*x2 +
237  & (-4.d0*v1-2.d0*v3)*x3 +
238  & (2.d0*u1+2.d0*u3+2.d0*u2)*y2 +
239  & (4.d0*u1+2.d0*u3)*y3) * xsur45
240  a65(ielem)=((2.d0*v3+2.d0*v2+2.d0*v1)*x2 +
241  & (-2.d0*v2-6.d0*v3-4.d0*v1)*x3 +
242  & (-2.d0*u2-2.d0*u3-2.d0*u1)*y2 +
243  & (4.d0*u1+6.d0*u3+2.d0*u2)*y3) * xsur45
244 !
245 ! DIAGONAL TERMS:
246 !
247  a11(ielem) = - a12(ielem) - a13(ielem) - a14(ielem)
248  & - a15(ielem) - a16(ielem)
249  a22(ielem) = - a21(ielem) - a23(ielem) - a24(ielem)
250  & - a25(ielem) - a26(ielem)
251  a33(ielem) = - a31(ielem) - a32(ielem) - a34(ielem)
252  & - a35(ielem) - a36(ielem)
253  a44(ielem) = - a41(ielem) - a42(ielem) - a43(ielem)
254  & - a45(ielem) - a46(ielem)
255  a55(ielem) = - a51(ielem) - a52(ielem) - a53(ielem)
256  & - a54(ielem) - a56(ielem)
257  a66(ielem) = - a61(ielem) - a62(ielem) - a63(ielem)
258  & - a64(ielem) - a65(ielem)
259 !
260  ENDDO ! IELEM
261 !
262 !-----------------------------------------------------------------------
263 !
264 ! CASE WHERE U AND V ARE LINEAR
265 !
266  ELSEIF(ielmu.EQ.11.AND.ielmv.EQ.11) THEN
267 !
268 !-----------------------------------------------------------------------
269 !
270 ! P1 DISCRETISATION OF THE VELOCITY:
271 !
272  DO ielem = 1 , nelem
273 !
274 ! INITIALISES THE GEOMETRICAL VARIABLES
275 !
276  x2 = xel(ielem,2)
277  x3 = xel(ielem,3)
278 !
279  y2 = yel(ielem,2)
280  y3 = yel(ielem,3)
281 !
282  u1 = u(ikle1(ielem))
283  u2 = u(ikle2(ielem))
284  u3 = u(ikle3(ielem))
285  v1 = v(ikle1(ielem))
286  v2 = v(ikle2(ielem))
287  v3 = v(ikle3(ielem))
288 !
289 ! EXTRADIAGONAL TERMS
290 !
291  a12(ielem)=(( 5.d0*v2+v3+v1*6.d0)*x3 +
292  & ( -u1*6.d0-u3-u2*5.d0)*y3) * xsu360
293  a13(ielem)=((-v2-v3*5.d0-v1*6.d0)*x2 +
294  & ( u1*6.d0+u3*5.d0+u2)*y2) * xsu360
295  a14(ielem)=(( v2*2.d0+v3)*x2+(-v2*2.d0-v3-v1*6.d0)*x3 +
296  & ( -u2*2.d0-u3)*y2+(u1*6.d0+u3+u2*2.d0)*y3) * xsur90
297  a15(ielem)=((-v3-v2*2.d0)*x2+(v2+v3*2.d0)*x3 +
298  & (u2*2.d0+u3)*y2+(-u2-u3*2.d0)*y3) * xsur90
299  a16(ielem)=((v2+v3*2.d0+v1*6.d0)*x2+(-v2-v3*2.d0)*x3 +
300  & (-u2-u3*2.d0-u1*6.d0)*y2+(u2+u3*2.d0)*y3) * xsur90
301 !
302  a21(ielem)=((v2*6.d0+v3+v1*5.d0)*x2 +
303  & (-v2*6.d0-v3-v1*5.d0)*x3 +
304  & (-u2*6.d0-u3-u1*5.d0)*y2 +
305  & (u1*5.d0+u3+u2*6.d0)*y3) * xsu360
306  a23(ielem)=((-v2*6.d0-v3*5.d0-v1)*x2 +
307  & ( u1+u3*5.d0+u2*6.d0)*y2) * xsu360
308  a24(ielem)=(-v2*x2*6.d0+(v2*6.d0+v3+v1*2.d0)*x3+u2*y2*6.d0 +
309  & (-u2*6.d0-u3-u1*2.d0)*y3) * xsur90
310  a25(ielem)=(v2*x2*6.d0+(v3*2.d0+v1)*x3-u2*y2*6.d0 +
311  & (-u3*2.d0-u1)*y3) * xsur90
312  a26(ielem)=((v3-v1)*x2+(-v3*2.d0-v1)*x3+(-u3+u1)*y2 +
313  & (u1+u3*2.d0)*y3) * xsur90
314 !
315  a31(ielem)=((v2+v3*6.d0+v1*5.d0)*x2 +
316  & (-v2-v3*6.d0-v1*5.d0)*x3 +
317  & (-u2-u3*6.d0-u1*5.d0)*y2 +
318  & (u1*5.d0+u3*6.d0+u2)*y3) * xsu360
319  a32(ielem)=((v2*5.d0+v3*6.d0+v1)*x3 +
320  & (-u1-u3*6.d0-u2*5.d0)*y3) * xsu360
321  a34(ielem)=((v2*2.d0+v1)*x2+(-v2+v1)*x3+(-u1-u2*2.d0)*y2+
322  & (u2-u1)*y3)*xsur90
323  a35(ielem)=((-v2*2.d0-v1)*x2-v3*x3*6.d0 +
324  & ( u2*2.d0+u1)*y2+u3*y3*6.d0) * xsur90
325  a36(ielem)=((-v3*6.d0-v2-v1*2.d0)*x2+v3*x3*6.d0 +
326  & ( u2+u3*6.d0+u1*2.d0)*y2-u3*y3*6.d0) * xsur90
327 !
328  a41(ielem)=((-v2*2.d0-v3-v1*6.d0)*x2 +
329  & ( v2*2.d0+v3+v1*6.d0)*x3 +
330  & ( u1*6.d0+u3+u2*2.d0)*y2 +
331  & (-u1*6.d0-u3-u2*2.d0)*y3) * xsur90
332  a42(ielem)=((-v2*6.d0-v3-v1*2.d0)*x3 +
333  & ( u1*2.d0+u3+u2*6.d0)*y3) * xsur90
334  a43(ielem)=((-v2*2.d0+v3-v1*2.d0)*x2 +
335  & ( u1*2.d0-u3+u2*2.d0)*y2) * xsur90
336  a45(ielem)=(( 6.d0*v2+2.d0*v3+4.d0*v1)*x2 +
337  & (-2.d0*v2-2.d0*v3-2.d0*v1)*x3 +
338  & (-6.d0*u2-2.d0*u3-4.d0*u1)*y2 +
339  & (2.d0*u1+2.d0*u3+2.d0*u2)*y3) * xsur45
340  a46(ielem)=((2.d0*v2+4.d0*v1)*x2 +
341  & (2.d0*v3+2.d0*v2+2.d0*v1)*x3 +
342  & (-4.d0*u1-2.d0*u2)*y2 +
343  & (-2.d0*u2-2.d0*u3-2.d0*u1)*y3) * xsur45
344 !
345  a51(ielem)=((v2*2.d0+v3*2.d0-v1)*x2 +
346  & (-v2*2.d0-v3*2.d0+v1)*x3 +
347  & (-u2*2.d0-u3*2.d0+u1)*y2 +
348  & (-u1+u3*2.d0+u2*2.d0)*y3) * xsur90
349  a52(ielem)=((-v2*6.d0-v3*2.d0-v1)*x3 +
350  & (u1+u3*2.d0+u2*6.d0)*y3) * xsur90
351  a53(ielem)=((v2*2.d0+v3*6.d0+v1)*x2 +
352  & (-u1-u3*6.d0-u2*2.d0)*y2) * xsur90
353  a54(ielem)=((-6.d0*v2-4.d0*v3-2.d0*v1)*x2 +
354  & (4.d0*v2+2.d0*v3)*x3 +
355  & (6.d0*u2+4.d0*u3+2.d0*u1)*y2 +
356  & (-4.d0*u2-2.d0*u3)*y3) * xsur45
357  a56(ielem)=((-2.d0*v2-4.d0*v3)*x2 +
358  & (4.d0*v2+6.d0*v3+2.d0*v1)*x3 +
359  & (2.d0*u2+4.d0*u3)*y2 +
360  & (-2.d0*u1-6.d0*u3-4.d0*u2)*y3) * xsur45
361 !
362  a61(ielem)=((-v2-v3*2.d0-v1*6.d0)*x2 +
363  & (v2+v3*2.d0+v1*6.d0)*x3 +
364  & (u2+u3*2.d0+u1*6.d0)*y2 +
365  & (-u2-u3*2.d0-u1*6.d0)*y3) * xsur90
366  a62(ielem)=((-v2+v3*2.d0+v1*2.d0)*x3 +
367  & (-u1*2.d0-u3*2.d0+u2)*y3) * xsur90
368  a63(ielem)=((v2+v3*6.d0+v1*2.d0)*x2 +
369  & (-u1*2.d0-u3*6.d0-u2)*y2) * xsur90
370  a64(ielem)=((-2.d0*v2-2.d0*v3-2.d0*v1)*x2 +
371  & (-4.d0*v1-2.d0*v3)*x3 +
372  & (2.d0*u1+2.d0*u3+2.d0*u2)*y2 +
373  & (4.d0*u1+2.d0*u3)*y3) * xsur45
374  a65(ielem)=((2.d0*v3+2.d0*v2+2.d0*v1)*x2 +
375  & (-2.d0*v2-6.d0*v3-4.d0*v1)*x3 +
376  & (-2.d0*u2-2.d0*u3-2.d0*u1)*y2 +
377  & (4.d0*u1+6.d0*u3+2.d0*u2)*y3) * xsur45
378 !
379 ! DIAGONAL TERMS:
380 !
381  a11(ielem) = - a12(ielem) - a13(ielem) - a14(ielem)
382  & - a15(ielem) - a16(ielem)
383  a22(ielem) = - a21(ielem) - a23(ielem) - a24(ielem)
384  & - a25(ielem) - a26(ielem)
385  a33(ielem) = - a31(ielem) - a32(ielem) - a34(ielem)
386  & - a35(ielem) - a36(ielem)
387  a44(ielem) = - a41(ielem) - a42(ielem) - a43(ielem)
388  & - a45(ielem) - a46(ielem)
389  a55(ielem) = - a51(ielem) - a52(ielem) - a53(ielem)
390  & - a54(ielem) - a56(ielem)
391  a66(ielem) = - a61(ielem) - a62(ielem) - a63(ielem)
392  & - a64(ielem) - a65(ielem)
393  ENDDO ! IELEM
394 !
395 !-----------------------------------------------------------------------
396 !
397  ELSEIF(ielmu.EQ.13.AND.ielmv.EQ.13) THEN
398 !
399 !-----------------------------------------------------------------------
400 !
401 ! P2 DISCRETISATION OF THE VELOCITY:
402 !
403  IF(formul(16:16).EQ.'N') THEN
404 !
405 ! N SCHEME
406 !
407  WRITE(lu,13)
408  CALL plante(1)
409  ELSE
410 !
411  DO ielem = 1 , nelem
412 !
413 ! TRADITIONAL METHOD
414 !
415 !
416 ! INITIALISES THE GEOMETRICAL VARIABLES
417 !
418  x2 = xel(ielem,2)
419  x3 = xel(ielem,3)
420 !
421  y2 = yel(ielem,2)
422  y3 = yel(ielem,3)
423 !
424  u1 = u(ikle1(ielem))
425  u2 = u(ikle2(ielem))
426  u3 = u(ikle3(ielem))
427  u4 = u(ikle4(ielem))
428  u5 = u(ikle5(ielem))
429  u6 = u(ikle6(ielem))
430  v1 = v(ikle1(ielem))
431  v2 = v(ikle2(ielem))
432  v3 = v(ikle3(ielem))
433  v4 = v(ikle4(ielem))
434  v5 = v(ikle5(ielem))
435  v6 = v(ikle6(ielem))
436 !
437 ! EXTRADIAGONAL TERMS
438 !
439  a12(ielem)= ((v5*20.d0+4.d0*8.d0*v4-11.d0*v3+16.d0*v6 +
440  & v1*18.d0+v2*9.d0)*x3+(-16.d0*u6-u5*20.d0+
441  & 11.d0*u3-32.d0*u4-u1*18.d0-u2*9.d0)*y3) * xsur2520
442  a13(ielem)= ((-16.d0*v4-v5*20.d0-4.d0*8.d0*v6-v3*9.d0 -
443  & v1*18.d0+11.d0*v2)*x2+(u1*18.d0+u5*20.d0 +
444  & u3*9.d0+32.d0*u6+16.d0*u4-11.d0*u2)*y2) * xsur2520
445  a14(ielem)= ((-v3+8.d0*v4+4.d0*v6-v1*6.d0+12.d0*v5 +
446  & 2.d0*2.d0*v2)*x2+(-20.d0*v4+v3*5.d0-16.d0*v6-
447  & 24.d0*v1-8.d0*v5)*x3+(-4.d0*u2-8.d0*u4-12.d0*u5+
448  & u3+u1*6.d0-4.d0*u6)*y2+(16.d0*u6+20.d0*u4+
449  & 24.d0*u1-u3*5.d0+4.d0*2.d0*u5)*y3) * xsur630
450  a15(ielem)= ((v3-12.d0*v5-4.d0*v6+v1*6.d0-8.d0*v4-4.d0*v2)*x2 +
451  & (4.d0*v4-v1*6.d0-v2+12.d0*v5+4.d0*v3+8.d0*v6)*x3 +
452  & (4.d0*u2+12.d0*u5-u1*6.d0+8.d0*u4+4.d0*u6-u3)*y2 +
453  & (-8.d0*u6-4.d0*u4+u1*6.d0+u2-4.d0*u3-12.d0*u5)*y3)
454  & * xsur630
455  a16(ielem)= (( 16.d0*v4+20.d0*v6+24.d0*v1+8.d0*v5-v2*5.d0)*x2 +
456  & (-8.d0*v6-4.d0*v4+v1*6.d0-4.d0*v3-12.d0*v5+v2)*x3+
457  & (-24.d0*u1-8.d0*u5+u2*5.d0-20.d0*u6-16.d0*u4)*y2 +
458  & (8.d0*u6+4.d0*u4-u1*6.d0+4.d0*u3+12.d0*u5-u2)*y3)
459  & * xsur630
460 !
461  a21(ielem)= ((-11.d0*v3+2.d0*8.d0*v5+v6*20.d0+v1*9.d0+32.d0*v4+
462  & v2*18.d0)*x2+(-4.d0*8.d0*v4-v1*9.d0-v2*18.d0-
463  & 2.d0*8.d0*v5+11.d0*v3-v6*20.d0)*x3+(-u2*18.d0-
464  & 2.d0*8.d0*u5-u1*9.d0-32.d0*u4-u6*20.d0+11.d0*u3)*y2+
465  & (u6*20.d0+4.d0*8.d0*u4+u1*9.d0+u2*18.d0-11.d0*u3+
466  & 2.d0*8.d0*u5)*y3) * xsur2520
467  a23(ielem)= ((-2.d0*8.d0*v4-4.d0*8.d0*v5-v6*20.d0-v3*9.d0 +
468  & 11.d0*v1-v2*18.d0)*x2+(-11.d0*u1+4.d0*8.d0*u5 +
469  & u3*9.d0+u6*20.d0+2.d0*8.d0*u4+u2*18.d0)*y2)
470  & *xsur2520
471  a24(ielem)=((4.d0*v3-12.d0*v4+4.d0*v1+4.d0*v6-12.d0*v5-
472  & v2*30.d0)*x2+(20.d0*v4-v3*5.d0+8.d0*v6+24.d0*v2+
473  & 16.d0*v5)*x3+(u2*30.d0+12.d0*u4+12.d0*u5-4.d0*u3-
474  & 4.d0*u1-4.d0*u6)*y2+(-4.d0*2.d0*u6-20.d0*u4-24.d0*u2+
475  & u3*5.d0-16.d0*u5)*y3)*xsur630
476  a25(ielem)=((-4.d0*v3+12.d0*v5-4.d0*v6-4.d0*v1+12.d0*v4+
477  & v2*30.d0)*x2+(4.d0*v4-v1-v2*6.d0+8.d0*v5+4.d0*v3+
478  & 12.d0*v6)*x3+(-u2*30.d0-12.d0*u5+4.d0*u1-12.d0*u4+
479  & 4.d0*u6+4.d0*u3)*y2+(-12.d0*u6-4.d0*u4+u1+u2*6.d0-
480  & 4.d0*u3-8.d0*u5)*y3)*xsur630
481  a26(ielem)=((4.d0*v5-4.d0*v4-v1*5.d0+v3*5.d0)*x2+(-12.d0*v6+
482  & v2*6.d0-4.d0*v4+v1-4.d0*v3-8.d0*v5)*x3+(-4.d0*u5-
483  & u3*5.d0+u1*5.d0+4.d0*u4)*y2+(12.d0*u6+4.d0*u4-u2*6.d0+
484  & 8.d0*u5+4.d0*u3-u1)*y3)*xsur630
485 !
486  a31(ielem)=((v3*18.d0+16.d0*v5+32.d0*v6+v1*9.d0+v4*20.d0-
487  & 11.d0*v2)*x2+(-v4*20.d0-v1*9.d0+11.d0*v2-2.d0*8.d0*v5-
488  & v3*18.d0-4.d0*8.d0*v6)*x3+(11.d0*u2-2.d0*8.d0*u5-
489  & u1*9.d0-u4*20.d0-4.d0*8.d0*u6-u3*18.d0)*y2+
490  & (4.d0*8.d0*u6+u4*20.d0+u1*9.d0-11.d0*u2+u3*18.d0+
491  & 16.d0*u5)*y3)*xsur2520
492  a32(ielem)=((32.d0*v5+v4*20.d0+v3*18.d0+16.d0*v6-11.d0*v1+
493  & v2*9.d0)*x3+(-16.d0*u6-32.d0*u5-u3*18.d0-u4*20.d0+
494  & 11.d0*u1-u2*9.d0)*y3)*xsur2520
495  a34(ielem)=((-v3*6.d0+8.d0*v5+4.d0*v6+12.d0*v4-v1+4.d0*v2)*x2+
496  & (-v2*5.d0-4.d0*v5+v1*5.d0+4.d0*v6)*x3+(u1-8.d0*u5+
497  & u3*6.d0-4.d0*u6-4.d0*u2-12.d0*u4)*y2+(-4.d0*u6+u2*5.d0+
498  & 4.d0*u5-u1*5.d0)*y3)*xsur630
499  a35(ielem)=((v3*6.d0-8.d0*v5-4.d0*v6+v1-12.d0*v4-4.d0*v2)*x2+
500  & (4.d0*v4+4.d0*v1+4.d0*v2-12.d0*v5-v3*30.d0-
501  & 12.d0*v6)*x3+(4.d0*u2+8.d0*u5-u1+12.d0*u4+4.d0*u6-
502  & u3*6.d0)*y2+(12.d0*u6-4.d0*u4-4.d0*u1-4.d0*u2+u3*30.d0+
503  & 12.d0*u5)*y3)*xsur630
504  a36(ielem)=((-24.d0*v3-8.d0*v4-20.d0*v6-16.d0*v5+v2*5.d0)*x2+
505  & (-4.d0*v4+v3*30.d0+12.d0*v6-4.d0*v2+12.d0*v5-
506  & 4.d0*v1)*x3+(-u2*5.d0+16.d0*u5+24.d0*u3+8.d0*u4+
507  & 20.d0*u6)*y2+(-12.d0*u6+4.d0*u4+4.d0*u2-u3*30.d0+
508  & 4.d0*u1-12.d0*u5)*y3)*xsur630
509 !
510  a41(ielem)=(-u3*y2*5.d0-20.d0*v6*x2-v3*x3*5.d0+12.*v1*x3-
511  & 12.d0*v1*x2-12.d0*u1*y3+12.d0*u1*y2-20.d0*u6*y3+
512  & v3*x2*5.d0+u3*y3*5.d0+20.d0*v6*x3+20.*u6*y2-
513  & 40.d0*u4*y3+40.*v4*x3-4.d0*u5*y3+4.d0*v5*x3+
514  & 40.d0*u4*y2+4.d0*u5*y2-4.d0*v5*x2-40.*v4*x2-
515  & 8.d0*v2*x3+8.d0*u2*y3+8.d0*v2*x2-
516  & 8.d0*u2*y2)*xsur630
517  a42(ielem)=((-20.d0*v5-40.d0*v4+5.d0*v3-4.d0*v6+8.d0*v1-
518  & 12.d0*v2)*x3+(4.d0*u6+20.d0*u5-5.d0*u3+40.d0*u4-
519  & 8.d0*u1+12.d0*u2)*y3) * xsur630
520  a43(ielem)=((-24.d0*v4+4.d0*v5+4.d0*v6+3.d0*v3-4.d0*v1-
521  & 4.d0*v2)*x2+(4.d0*u1-4.d0*u5-3.d0*u3-4.d0*u6+
522  & 24.d0*u4+4.d0*u2)*y2) * xsur630
523  a45(ielem)=(12.d0*u3*y2+32.d0*v6*x2+4.d0*v3*x3+4.d0*v1*x3-
524  & 8.d0*v1*x2-4.d0*u1*y3+8.d0*u1*y2+32.d0*u6*y3-
525  & 12.d0*v3*x2-4.d0*u3*y3-32.d0*v6*x3-32.d0*u6*y2+
526  & 32.d0*u4*y3-32.d0*v4*x3+32.d0*u5*y3-32.d0*v5*x3-
527  & 96.d0*u4*y2-48.d0*u5*y2+48.d0*v5*x2+96.d0*v4*x2+
528  & 4.d0*v2*x3-4.d0*u2*y3+12.d0*v2*x2-12.d0*u2*y2)
529  & * xsur630
530  a46(ielem)=((-8.d0*v3+16.d0*v6+16.d0*v1+64.d0*v4-4.d0*v2)*x2+
531  & (-4.d0*v1+32.d0*v4-4.d0*v3+32.d0*v6+32.d0*v5-
532  & 4.d0*v2)*x3+(8.d0*u3-16.d0*u1+4.d0*u2-16.d0*u6-
533  & 64.d0*u4)*y2+(-32.d0*u6+4.d0*u1-32.d0*u5+4.d0*u3-
534  & 32.d0*u4+4.d0*u2)*y3)*xsur630
535 !
536  a51(ielem)=((4.d0*v3+24.d0*v5-4.d0*v6-v1*3.d0-4.d0*v4+
537  & 4.d0*v2)*x2+(4.d0*v4+v1*3.d0-4.d0*v2-24.d0*v5-
538  & 4.d0*v3+4.d0*v6)*x3+(-4.d0*u2-24.d0*u5+u1*3.d0+
539  & 4.d0*u4+4.d0*u6-4.d0*u3)*y2+(-4.d0*u6-4.d0*u4-u1*3.d0+
540  & 4.d0*u2+4.d0*u3+24.d0*u5)*y3)*xsur630
541  a52(ielem)=((-40.d0*v5-20.d0*v4+8.d0*v3-4.d0*v6+v1*5.d0-
542  & 12.d0*v2)*x3+(4.d0*u6+40.d0*u5-8.d0*u3+
543  & 20.d0*u4-u1*5.d0+12.d0*u2)*y3)*xsur630
544  a53(ielem)=((4.d0*v4+40.d0*v5+20.d0*v6+2.d0*6.d0*v3-v1*5.d0-
545  & 8.d0*v2)*x2+(u1*5.d0-40.d0*u5-12.d0*u3-20.d0*u6-
546  & 4.d0*u4+8.d0*u2)*y2)*xsur630
547  a54(ielem)=((8.d0*v3-96.d0*v5-32.d0*v6+12.d0*v1-48.d0*v4-
548  & 12.d0*v2)*x2+(16.d0*v4-8.d0*v1+16.d0*v2+64.d0*v5-
549  & 4.d0*v3)*x3+(-12.d0*u1+96.d0*u5-8.d0*u3+12.d0*u2+
550  & 32.d0*u6+48.d0*u4)*y2+(8.d0*u1+4.d0*u3-16.d0*u4-
551  & 64.d0*u5-16.d0*u2)*y3)*xsur630
552  a56(ielem)=((-16.d0*v3-16.d0*v6+8.d0*v1-64.d0*v5+4.d0*v2)*x2+
553  & (32.d0*v4+12.d0*v3-12.d0*v1+48.d0*v6+96.d0*v5-
554  & 8.d0*v2)*x3+(-8.d0*u1-4.d0*u2+16.d0*u3+16.d0*u6+
555  & 64.d0*u5)*y2+(-48.d0*u6+12.d0*u1+8.d0*u2-12.d0*u3-
556  & 32.d0*u4-96.d0*u5)*y3)*xsur630
557 !
558  a61(ielem)=((8.d0*v3-4.d0*v5-40.d0*v6-12.d0*v1-20.d0*v4+
559  & v2*5.d0)*x2+(20.d0*v4+12.d0*v1-v2*5.d0+4.d0*v5-
560  & 8.d0*v3+40.d0*v6)*x3+(-u2*5.d0+4.d0*u5+12.d0*u1+
561  & 20.d0*u4+40.d0*u6-8.d0*u3)*y2+(-40.d0*u6-20.d0*u4-
562  & 12.d0*u1+u2*5.d0+8.d0*u3-4.d0*u5)*y3)*xsur630
563  a62(ielem)=((-4.d0*v5-4.d0*v4+4.d0*v3+24.d0*v6+4.d0*v1-
564  & v2*3.d0)*x3+(-24.d0*u6+4.d0*u5-4.d0*u3+4.d0*u4-
565  & 4.d0*u1+u2*3.d0)*y3)*xsur630
566  a63(ielem)=((4.d0*v4+20.d0*v5+40.d0*v6+12.d0*v3-8.d0*v1-
567  & v2*5.d0)*x2+(8.d0*u1-20.d0*u5-12.d0*u3-40.d0*u6-
568  & 4.d0*u4+u2*5.d0)*y2)*xsur630
569  a64(ielem)=((-32.d0*v4+4.d0*v1+4.d0*v2-32.d0*v5+4.d0*v3-
570  & 32.d0*v6)*x2+(-16.d0*v4-16.d0*v1+8.d0*v2+4.d0*v3-
571  & 64.d0*v6)*x3+(32.d0*u6+32.d0*u4-4.d0*u1-4.d0*u2-
572  & 4.d0*u3+32.d0*u5)*y2+(64.d0*u6-8.d0*u2+16.d0*u1-
573  & 4.d0*u3+16.d0*u4)*y3)*xsur630
574  a65(ielem)=((-4.d0*v1+32.d0*v4-4.d0*v3+32.d0*v6+32.d0*v5-
575  & 4.d0*v2)*x2+(-32.d0*v4+8.d0*v1+12.d0*v2-48.d0*v5-
576  & 12.d0*v3-96.d0*v6)*x3+(-32.d0*u6+4.d0*u1-32.d0*u5+
577  & 4.d0*u3-32.d0*u4+4.d0*u2)*y2+(96.d0*u6+32.d0*u4-
578  & 8.d0*u1-12.d0*u2+12.d0*u3+48.d0*u5)*y3)*xsur630
579 !
580 ! THE DIAGONAL TERMS ARE OBTAINED BY MEANS OF THE 'MAGIC SQUARE':
581 !
582  a11(ielem) = - a12(ielem) - a13(ielem) - a14(ielem)
583  & - a15(ielem) - a16(ielem)
584  a22(ielem) = - a21(ielem) - a23(ielem) - a24(ielem)
585  & - a25(ielem) - a26(ielem)
586  a33(ielem) = - a31(ielem) - a32(ielem) - a34(ielem)
587  & - a35(ielem) - a36(ielem)
588  a44(ielem) = - a41(ielem) - a42(ielem) - a43(ielem)
589  & - a45(ielem) - a46(ielem)
590  a55(ielem) = - a51(ielem) - a52(ielem) - a53(ielem)
591  & - a54(ielem) - a56(ielem)
592  a66(ielem) = - a61(ielem) - a62(ielem) - a63(ielem)
593  & - a64(ielem) - a65(ielem)
594 !
595  ENDDO ! IELEM
596 !
597  ENDIF
598 !
599 !-----------------------------------------------------------------------
600 !
601  ELSE
602 !
603  WRITE(lu,11) ielmu,ielmv
604 11 FORMAT(1x,
605  & 'MT05CC (BIEF) : TYPES OF VELOCITIES NOT AVAILABLE : ',2i6)
606 13 FORMAT(1x,
607  & 'MT05CC (BIEF) : N SCHEMES NOT AVAILABLE ')
608  CALL plante(1)
609  stop
610 !
611  ENDIF
612 !
613 !-----------------------------------------------------------------------
614 !
615  RETURN
616  END
declarations_special
Definition: declarations_special.F:3
mt05cc
subroutine mt05cc(A11, A12, A13, A14, A15, A16, A21, A22, A23, A24, A25, A26, A31, A32, A33, A34, A35, A36, A41, A42, A43, A44, A45, A46, A51, A52, A53, A54, A55, A56, A61, A62, A63, A64, A65, A66, XMUL, SU, SV, U, V, XEL, YEL, IKLE1, IKLE2, IKLE3, IKLE4, IKLE5, IKLE6, NELEM, NELMAX, FORMUL)
Definition: mt05cc.f:15
declarations_special::lu
integer lu
Definition: declarations_special.F:45
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:07 for The TELEMAC-MASCARET system by doxygen 1.8.13