The TELEMAC-MASCARET system  trunk
mt02cc.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt02cc
3 ! *****************
4 !
5  &( a11 , a12 , a13 , a14 , a15, a16,
6  & a22 , a23 , a24 , a25, a26,
7  & a33 , a34 , a35, a36,
8  & a44 , a45, a46,
9  & a55, a56,
10  & a66,
11  & xmul,su,u,xel,yel,surfac,ikle1,ikle2,ikle3,nelem,nelmax)
12 !
13 !***********************************************************************
14 ! BIEF V6P1 21/08/2010
15 !***********************************************************************
16 !
17 !brief BUILDS THE DIFFUSION MATRIX FOR P2 TRIANGLES.
18 !+
19 !+ VISCOSITY CAN BE ISOTROPIC, OR NOT ISOTROPIC. IN THIS
20 !+ CASE U IS AN ARRAY WITH SECOND DIMENSION EQUAL TO 3.
21 !
22 !warning THE JACOBIAN MUST BE POSITIVE
23 !
24 !history ALGIANE FROEHLY (MATMECA)
25 !+ 19/06/08
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 !| A12 |<--| ELEMENTS OF MATRIX
44 !| A13 |<--| ELEMENTS OF MATRIX
45 !| A14 |<--| ELEMENTS OF MATRIX
46 !| A15 |<--| ELEMENTS OF MATRIX
47 !| A16 |<--| ELEMENTS OF MATRIX
48 !| A22 |<--| ELEMENTS OF MATRIX
49 !| A23 |<--| ELEMENTS OF MATRIX
50 !| A24 |<--| ELEMENTS OF MATRIX
51 !| A25 |<--| ELEMENTS OF MATRIX
52 !| A26 |<--| ELEMENTS OF MATRIX
53 !| A33 |<--| ELEMENTS OF MATRIX
54 !| A34 |<--| ELEMENTS OF MATRIX
55 !| A35 |<--| ELEMENTS OF MATRIX
56 !| A36 |<--| ELEMENTS OF MATRIX
57 !| A44 |<--| ELEMENTS OF MATRIX
58 !| A45 |<--| ELEMENTS OF MATRIX
59 !| A46 |<--| ELEMENTS OF MATRIX
60 !| FORMUL |-->| FORMULA DESCRIBING THE MATRIX
61 !| IKLE1 |-->| FIRST POINTS OF TRIANGLES
62 !| IKLE2 |-->| SECOND POINTS OF TRIANGLES
63 !| IKLE3 |-->| THIRD POINTS OF TRIANGLES
64 !| NELEM |-->| NUMBER OF ELEMENTS
65 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
66 !| SU |-->| BIEF_OBJ STRUCTURE OF U
67 !| SURFAC |-->| AREA OF TRIANGLES
68 !| SV |-->| BIEF_OBJ STRUCTURE OF V
69 !| U |-->| FUNCTION U USED IN THE FORMULA
70 !| V |-->| FUNCTION V USED IN THE FORMULA
71 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
72 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
73 !| XMUL |-->| MULTIPLICATION FACTOR
74 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
75 !
76  USE bief!, EX_MT02CC => MT02CC
77 !
78  USE declarations_special
79  IMPLICIT NONE
80 !
81 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
82 !
83  INTEGER, INTENT(IN) :: NELEM,NELMAX
84  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax),IKLE3(nelmax)
85  DOUBLE PRECISION, INTENT(INOUT) :: A11(*),A12(*),A13(*)
86  DOUBLE PRECISION, INTENT(INOUT) :: A14(*),A15(*),A16(*)
87  DOUBLE PRECISION, INTENT(INOUT) :: A22(*),A23(*)
88  DOUBLE PRECISION, INTENT(INOUT) :: A24(*),A25(*),A26(*)
89  DOUBLE PRECISION, INTENT(INOUT) :: A33(*),A34(*)
90  DOUBLE PRECISION, INTENT(INOUT) :: A35(*),A36(*)
91  DOUBLE PRECISION, INTENT(INOUT) :: A44(*),A45(*),A46(*)
92  DOUBLE PRECISION, INTENT(INOUT) :: A55(*),A56(*)
93  DOUBLE PRECISION, INTENT(INOUT) :: A66(*)
94  DOUBLE PRECISION, INTENT(IN) :: XMUL,U(*)
95 ! STRUCTURE OF U
96  TYPE(bief_obj), INTENT(IN) :: SU
97 !
98  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,3),YEL(nelmax,3)
99  DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
100 !
101 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
102 !
103 ! DECLARATIONS SPECIFIC TO THIS SUBROUTINE
104 !
105  INTEGER IELMNU,IELEM,ISO,IAD2,IAD3
106 !
107  DOUBLE PRECISION X2,X3,Y2,Y3,AUX1,AUX2,AUX3,AUX4
108  DOUBLE PRECISION NUX1,NUX2,NUX3
109  DOUBLE PRECISION NUY1,NUY2,NUY3
110  DOUBLE PRECISION NUZ1,NUZ2,NUZ3
111 !
112 !=======================================================================
113 !
114 ! EXTRACTS THE TYPE OF ELEMENT FOR VISCOSITY
115 !
116  ielmnu = su%ELM
117  iso = su%DIM2
118 !
119 ! IF(IELMNU.EQ.10.AND.ISO.EQ.1) THEN
120 !
121 !-----------------------------------------------------------------------
122 !
123 ! P0 DISCRETISATION FOR VISCOSITY:
124 !
125 !-----------------------------------------------------------------------
126 !
127  IF(ielmnu.EQ.11.AND.iso.EQ.1) THEN
128 !
129 !-----------------------------------------------------------------------
130 !
131 ! P1 DISCRETISATION FOR ISOTROPIC VISCOSITY:
132 !
133  DO ielem = 1 , nelem
134 !
135 ! INITIALISES THE GEOMETRICAL VARIABLES
136 !
137  x2 = xel(ielem,2)
138  x3 = xel(ielem,3)
139 !
140  y2 = yel(ielem,2)
141  y3 = yel(ielem,3)
142 !
143  nux1 = u(ikle1(ielem))
144  nux2 = u(ikle2(ielem))
145  nux3 = u(ikle3(ielem))
146 !
147  aux1 = xmul/(60.d0*surfac(ielem))
148  aux2 = 3.d0 * aux1
149  aux3 = 4.d0 * aux1
150  aux4 = 8.d0 * aux1
151 !
152 ! EXTRADIAGONAL TERMS
153 !
154  a12(ielem)= -(-y3**2+y3*y2-x3**2+x3*x2) *
155  & (2.d0*nux1+2.d0*nux2+nux3) * aux1
156 !
157  a13(ielem)= (-y3*y2+y2**2-x3*x2+x2**2) *
158  & (2.d0*nux1+nux2+2.d0*nux3) * aux1
159 !
160  a14(ielem)= ((-11.d0*nux1-4.d0*nux3-5.d0*nux2) * (y3**2+x3**2)
161  & + (3.d0*nux1-nux3-2.d0*nux2 ) * (y2**2+x2**2)
162  & + (8.d0*nux1+7.d0*nux2+5.d0*nux3 ) * (x3*x2+y3*y2))
163  & * aux1
164 !
165  a15(ielem)=-((3.d0*nux1-2.d0*nux3-nux2 ) * (y3**2+x3**2)
166  & + (3.d0*nux1-nux3-2.d0*nux2 ) * (y2**2+x2**2)
167  & + (-6.d0*nux1+3.d0*nux2+3.d0*nux3 ) * (x3*x2+y3*y2))
168  & * aux1
169 !
170  a16(ielem)=-((-3.d0*nux1+2.d0*nux3+nux2 ) * (y3**2+x3**2)
171  & + (11.d0*nux1+5.d0*nux3+4.d0*nux2 ) * (y2**2+x2**2)
172  & + (-8.d0*nux1-5.d0*nux2-7.d0*nux3 ) * (x3*x2+y3*y2))
173  & * aux1
174 !
175  a23(ielem)= (y3*y2+x3*x2) * (nux1+2.d0*nux3+2.d0*nux2) * aux1
176 !
177  a24(ielem)= ((-5.d0*nux1-4.d0*nux3-11.d0*nux2) * (y3**2+x3**2)
178  & + (3.d0*nux1+14.d0*nux2+3.d0*nux3 ) * (x3*x2+y3*y2))
179  & * aux1
180 !
181  a25(ielem)=-((nux1+2.d0*nux3-3.d0*nux2 ) * (y3**2+x3**2)
182  & + (3.d0*nux1+3.d0*nux3+14.d0*nux2 ) * (y3*y2+x3*x2))
183  & * aux1
184 !
185  a26(ielem)= ((nux1+2.d0*nux3-3.d0*nux2 ) * (y3**2+x3**2)
186  & + (nux1-nux3) * (x3*x2+y3*y2))
187  & * aux1
188 !
189  a34(ielem)= ((nux1-3.d0*nux3+2.d0*nux2 ) * (y2**2+x2**2)
190  & + (nux1-nux2) * (x3*x2+y3*y2))
191  & * aux1
192 !
193  a35(ielem)=-((nux1-3.d0*nux3+2.d0*nux2 ) * (y2**2+x2**2)
194  & + (3.d0*nux1+3.d0*nux2+14.d0*nux3) * (x3*x2+y3*y2))
195  & * aux1
196 !
197  a36(ielem)=-((5.d0*nux1+11.d0*nux3+4.d0*nux2) * (y2**2+x2**2)
198  & + (-3.d0*(nux1+nux2)-14.d0*nux3 ) * (x3*x2+y3*y2))
199  & * aux1
200 !
201  a45(ielem)=-((-nux1+nux2 ) * (y3**2+x3**2)
202  & + (-nux1-6.d0*nux2-3.d0*nux3 ) * (x3*x2+y3*y2)
203  & + (6.d0*nux2+2.d0*nux1+2.d0*nux3 ) * (y2**2+x2**2))
204  & * aux3
205 !
206  a46(ielem)=-((nux1-nux2 ) * (y3**2+x3**2)
207  & + (nux1-nux3 ) * (y2**2+x2**2)
208  & + (4.d0*nux1+3.d0*nux2+3.d0*nux3 ) * (x3*x2+y3*y2))
209  & * aux3
210 !
211  a56(ielem)= ((-2.d0*nux1-2.d0*nux2-6.d0*nux3) * (y3**2+x3**2)
212  & + (nux1-nux3 ) * (y2**2+x2**2)
213  & + (nux1+3.d0*nux2+6.d0*nux3 ) * (x3*x2+y3*y2))
214  & * aux3
215 !
216 ! DIAGONAL TERMS
217 !
218  a11(ielem)= ((y3-y2)**2+(x3-x2)**2)
219  & * (3.d0*nux1+nux2+nux3) * aux2
220 !
221  a22(ielem)= (y3**2+x3**2) * (nux1+3.d0*nux2+nux3) * aux2
222 !
223  a33(ielem)= (y2**2+x2**2) * (nux1+nux2+3.d0*nux3) * aux2
224 !
225  a44(ielem)= ((2.d0*nux1+nux3+2.d0*nux2) * (y3**2+x3**2)
226  & + (nux1+nux3+3.d0*nux2 ) * (y2**2+x2**2)
227  & + (-4.d0*nux2-nux3 ) * (x3*x2+y3*y2))
228  & * aux4
229 !
230  a55(ielem)= ((nux1+nux2+3.d0*nux3 ) * (y3**2+x3**2)
231  & + (nux1+nux3+3.d0*nux2 ) * (y2**2+x2**2)
232  & + (-nux1-2.d0*nux2-2.d0*nux3) * (x3*x2+y3*y2))
233  & * aux4
234 !
235  a66(ielem) = ((nux1+3.d0*nux3+nux2 ) * (y3**2+x3**2)
236  & + (2.d0*nux1+2.d0*nux3+nux2 ) * (y2**2+x2**2)
237  & + (-nux2-4.d0*nux3 ) * (x3*x2+y3*y2))
238  & * aux4
239 !
240  ENDDO ! IELEM
241 !
242 !-----------------------------------------------------------------------
243 !
244  ELSEIF(ielmnu.EQ.11.AND.iso.EQ.3) THEN
245 !
246 !-----------------------------------------------------------------------
247 !
248 ! P1 DISCRETISATION FOR NONISOTROPIC VISCOSITY:
249 !
250  iad2 = su%MAXDIM1
251  iad3 = 2*iad2
252  DO ielem = 1 , nelem
253 !
254 ! INITIALISES THE GEOMETRICAL VARIABLES
255 !
256  x2 = xel(ielem,2)
257  x3 = xel(ielem,3)
258 !
259  y2 = yel(ielem,2)
260  y3 = yel(ielem,3)
261 !
262  nux1 = u(ikle1(ielem))
263  nux2 = u(ikle2(ielem))
264  nux3 = u(ikle3(ielem))
265  nuy1 = u(ikle1(ielem) + iad2)
266  nuy2 = u(ikle2(ielem) + iad2)
267  nuy3 = u(ikle3(ielem) + iad2)
268  nuz1 = u(ikle1(ielem) + iad3)
269  nuz2 = u(ikle2(ielem) + iad3)
270  nuz3 = u(ikle3(ielem) + iad3)
271 !
272  aux1 = xmul/(60.d0 * surfac(ielem))
273  aux2 = 4.d0 * aux1
274  aux3 = 3.d0 * aux1
275 !
276 ! EXTRADIAGONAL TERMS
277 !
278  a12(ielem) =
279  & ((nuy3+2.d0*nuy1+2.d0*nuy2) * (x3**2-x3*x2 ) +
280  & (nuz3+2.d0*nuz1+2.d0*nuz2) * (y2*x3+y3*x2-2.d0*y3*x3) +
281  & (nux3+2.d0*nux2+2.d0*nux1) * (y3**2-y3*y2 ) ) * aux1
282 !
283  a13(ielem) =
284  & ((2.d0*nuy1+2.d0*nuy3+nuy2) * (x2**2-x3*x2 ) +
285  & (2.d0*nuz1+2.d0*nuz3+nuz2) * (y3*x2+y2*x3-2.d0*y2*x2) +
286  & (2.d0*nux1+nux2+2.d0*nux3) * (y2**2-y3*y2 ) ) * aux1
287 !
288  a14(ielem) =
289  & -((- 3.d0*nux1+ 2.d0*nux2+ nux3) * y2**2 +
290  & (- 8.d0*nux1- 7.d0*nux2- 5.d0*nux3) * y3*y2 +
291  & ( 11.d0*nux1+ 5.d0*nux2+ 4.d0*nux3) * y3**2 +
292  & (- 3.d0*nuy1+ 2.d0*nuy2+ nuy3) * x2**2 +
293  & (- 8.d0*nuy1- 7.d0*nuy2- 5.d0*nuy3) * x3*x2 +
294  & ( 11.d0*nuy1+ 5.d0*nuy2+ 4.d0*nuy3) * x3**2 +
295  & ( 6.d0*nuz1- 4.d0*nuz2- 2.d0*nuz3) * y2*x2 +
296  & ( 8.d0*nuz1+ 7.d0*nuz2+ 5.d0*nuz3) * (y3*x2+y2*x3) +
297  & (-22.d0*nuz1-10.d0*nuz2- 8.d0*nuz3) * y3*x3 ) * aux1
298 !
299  a15(ielem) =
300  & -(( 3.d0*nux1- 2.d0*nux2- nux3) * y2**2 +
301  & (- 6.d0*nux1+ 3.d0*nux2+ 3.d0*nux3) * y3*y2 +
302  & ( 3.d0*nux1- nux2- 2.d0*nux3) * y3**2 +
303  & ( 3.d0*nuy1- 2.d0*nuy2- nuy3) * x2**2 +
304  & (- 6.d0*nuy1+ 3.d0*nuy2+ 3.d0*nuy3) * x3*x2 +
305  & ( 3.d0*nuy1- nuy2- 2.d0*nuy3) * x3**2 +
306  & (- 6.d0*nuz1+ 4.d0*nuz2+ 2.d0*nuz3) * y2*x2 +
307  & ( 6.d0*nuz1- 3.d0*nuz2- 3.d0*nuz3) * (y3*x2+y2*x3) +
308  & (- 6.d0*nuz1+ 2.d0*nuz2+ 4.d0*nuz3) * y3*x3 ) * aux1
309 !
310  a16(ielem) =
311  & ((-11.d0*nux1- 4.d0*nux2- 5.d0*nux3) * y2**2 +
312  & ( 8.d0*nux1+ 5.d0*nux2+ 7.d0*nux3) * y3*y2 +
313  & ( 3.d0*nux1- nux2- 2.d0*nux3) * y3**2 +
314  & (-11.d0*nuy1- 4.d0*nuy2- 5.d0*nuy3) * x2**2 +
315  & (+ 8.d0*nuy1+ 5.d0*nuy2+ 7.d0*nuy3) * x3*x2 +
316  & ( 3.d0*nuy1- nuy2- 2.d0*nuy3) * x3**2 +
317  & ( 22.d0*nuz1+ 8.d0*nuz2+10.d0*nuz3) * y2*x2 +
318  & (- 8.d0*nuz1- 5.d0*nuz2- 7.d0*nuz3) * (y3*x2+y2*x3) +
319  & (- 6.d0*nuz1+ 2.d0*nuz2+ 4.d0*nuz3) * y3*x3 ) * aux1
320 !
321  a23(ielem) =
322  & (( nuy1+2.d0*nuy2+2.d0*nuy3) * x3*x2 +
323  & (-nuz1-2.d0*nuz2-2.d0*nuz3) * (y3*x2+y2*x3) +
324  & ( nux1+2.d0*nux2+2.d0*nux3) * y3*y2 ) * aux1
325 !
326  a24(ielem) =
327  & -((- 3.d0*nux1-14.d0*nux2- 3.d0*nux3) * y3*y2 +
328  & ( 5.d0*nux1+11.d0*nux2+ 4.d0*nux3) * y3**2 +
329  & (- 3.d0*nuy1-14.d0*nuy2- 3.d0*nuy3) * x3*x2 +
330  & ( 5.d0*nuy1+11.d0*nuy2+ 4.d0*nuy3) * x3**2 +
331  & ( 3.d0*nuz1+14.d0*nuz2+ 3.d0*nuz3) * (y3*x2+y2*x3) +
332  & (-10.d0*nuz1-22.d0*nuz2- 8.d0*nuz3) * y3*x3 ) * aux1
333 !
334  a25(ielem) =
335  & -(( 3.d0*nux1+14.d0*nux2+ 3.d0*nux3) * y3*y2 +
336  & ( nux1- 3.d0*nux2+ 2.d0*nux3) * y3**2 +
337  & ( 3.d0*nuy1+14.d0*nuy2+ 3.d0*nuy3) * x3*x2 +
338  & ( nuy1- 3.d0*nuy2+ 2.d0*nuy3) * x3**2 +
339  & (- 3.d0*nuz1-14.d0*nuz2- 3.d0*nuz3) * (y3*x2+y2*x3) +
340  & (- 2.d0*nuz1+ 6.d0*nuz2- 4.d0*nuz3) * y3*x3 ) * aux1
341 !
342  a26(ielem) =
343  & (( nux1-nux3) * y3*y2 + (-nuy3+nuy1) * x3*x2 +
344  & ( nuz3-nuz1) * (y3*x2+x3*y2) +
345  & (2.d0*nuy3+nuy1-3.d0*nuy2) * x3**2 +
346  & (6.d0*nuz2-4.d0*nuz3-2.d0*nuz1) * y3*x3 +
347  & (2.d0*nux3-3.d0*nux2+ nux1) * y3**2 ) * aux1
348 !
349  a34(ielem) =
350  & (( nuy1-3.d0*nuy3+2.d0*nuy2) * x2**2 +
351  & (-2.d0*nuz1+6.d0*nuz3-4.d0*nuz2) * y2*x2 +
352  & ( nux1+2.d0*nux2-3.d0*nux3) * y2**2 +
353  & (nuy1-nuy2) * x3*x2 + (nux1-nux2) * y3*y2 +
354  & (nuz2-nuz1) * (y3*x2+y2*x3) )*aux1
355 !
356  a35(ielem) =
357  & -(( nux1+ 2.d0*nux2- 3.d0*nux3) * y2**2 +
358  & ( 3.d0*nux1+ 3.d0*nux2+14.d0*nux3) * y3*y2 +
359  & ( nuy1+ 2.d0*nuy2- 3.d0*nuy3) * x2**2 +
360  & ( 3.d0*nuy1+ 3.d0*nuy2+14.d0*nuy3) * x3*x2 +
361  & (- 3.d0*nuz1- 3.d0*nuz2-14.d0*nuz3) * (y3*x2+y2*x3) +
362  & (- 2.d0*nuz1- 4.d0*nuz2+ 6.d0*nuz3) * y2*x2 ) * aux1
363 !
364  a36(ielem) =
365  & ((- 5.d0*nux1- 4.d0*nux2-11.d0*nux3) * y2**2 +
366  & ( 3.d0*nux1+3.d0*nux2+ 14.d0*nux3) * y3*y2 +
367  & (- 5.d0*nuy1- 4.d0*nuy2-11.d0*nuy3) * x2**2 +
368  & ( 3.d0*nuy1+ 3.d0*nuy2+14.d0*nuy3) * x3*x2 +
369  & ( 10.d0*nuz1+ 8.d0*nuz2+22.d0*nuz3) * y2*x2 +
370  & (- 3.d0*nuz1- 3.d0*nuz2-14.d0*nuz3) * (y3*x2+y2*x3) ) * aux1
371 !
372  a45(ielem) =
373  & ((- 2.d0*nuy1- 6.d0*nuy2- 2.d0*nuy3) * x2**2 +
374  & ( nuy1+ 6.d0*nuy2+ 3.d0*nuy3) * x3*x2 +
375  & ( 4.d0*nuz1+12.d0*nuz2+ 4.d0*nuz3) * y2*x2 +
376  & (- nuz1- 6.d0*nuz2- 3.d0*nuz3) * (y3*x2+y2*x3) +
377  & (nuy1-nuy2)*x3**2 + (nux1-nux2)*y3**2 +
378  & (- 2.d0*nuz1+ 2.d0*nuz2 ) * y3*x3 +
379  & (- 2.d0*nux1- 2.d0*nux3- 6.d0*nux2) * y2**2 +
380  & ( 6.d0*nux2+ 3.d0*nux3+ nux1) * y3*y2 ) * aux2
381 !
382  a46(ielem) =
383  & -((-nuy3+nuy1)*x2**2 + (nux1-nux2) * y3**2 +
384  & ( nux1-nux3)*y2**2 + (nuy1-nuy2) * x3**2 +
385  & ( 4.d0*nuy1+ 3.d0*nuy3+ 3.d0*nuy2) * x3*x2 +
386  & ( 2.d0*nuz3- 2.d0*nuz1 ) * y2*x2 +
387  & ( -3.d0*nuz3- 4.d0*nuz1- 3.d0*nuz2) * (y3*x2+y2*x3) +
388  & (- 2.d0*nuz1+ 2.d0*nuz2 ) * y3*x3 +
389  & ( 3.d0*nux2+ 4.d0*nux1+ 3.d0*nux3) * y3*y2 ) * aux2
390 !
391  a56(ielem) =
392  & -(( nuy3-nuy1)*x2**2 + (nux3-nux1) * y2**2 +
393  & (- 3.d0*nuy2- 6.d0*nuy3- nuy1) * x3*x2 +
394  & ( 2.d0*nuz1- 2.d0*nuz3 ) * y2*x2 +
395  & ( 3.d0*nuz2+ nuz1+ 6.d0*nuz3) * (y3*x2+y2*x3) +
396  & ( 6.d0*nuy3+ 2.d0*nuy1+ 2.d0*nuy2) * x3**2 +
397  & (-12.d0*nuz3- 4.d0*nuz2- 4.d0*nuz1) * y3*x3 +
398  & (- 6.d0*nux3- 3.d0*nux2- nux1) * y3*y2 +
399  & ( 2.d0*nux2+ 2.d0*nux1+ 6.d0*nux3) * y3**2 ) * aux2
400 !
401 ! THE DIAGONAL TERMS ARE OBTAINED BY MEANS OF THE
402 ! MAGIC SQUARE OR BY DIRECT COMPUTATION:
403 !
404  a11(ielem) = - a12(ielem) - a13(ielem) - a14(ielem)
405  & - a15(ielem) - a16(ielem)
406  a22(ielem) = ((nuy1+nuy3+3.d0*nuy2) * x3**2 +
407  & (-6.d0*nuz2-2.d0*nuz1-2.d0*nuz3) * y3*x3 +
408  & (nux3+nux1+3.d0*nux2) * y3**2 ) * aux3
409  a33(ielem) = ((nuy1+3.d0*nuy3+nuy2) * x2**2 +
410  & (-6.d0*nuz3-2.d0*nuz1-2.d0*nuz2) * y2*x2 +
411  & ( 3.d0*nux3+nux1+nux2 ) * y2**2 ) * aux3
412  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
413  & - a45(ielem) - a46(ielem)
414  a55(ielem) = - a15(ielem) - a25(ielem) - a35(ielem)
415  & - a45(ielem) - a56(ielem)
416  a66(ielem) = - a16(ielem) - a26(ielem) - a36(ielem)
417  & - a46(ielem) - a56(ielem)
418 !
419  ENDDO ! IELEM
420 !
421 !-----------------------------------------------------------------------
422 !
423  ELSE
424 !
425  WRITE(lu,11) ielmnu,iso
426 11 FORMAT(1x,'MT02CC (BIEF) :TYPE OF VISCOSITY NOT AVAILABLE: ',2i6)
427  CALL plante(1)
428  stop
429 !
430  ENDIF
431 !
432 !-----------------------------------------------------------------------
433 !
434  RETURN
435  END
declarations_special
Definition: declarations_special.F:3
mt02cc
subroutine mt02cc(A11, A12, A13, A14, A15, A16, A22, A23, A24, A25, A26, A33, A34, A35, A36, A44, A45, A46, A55, A56, A66, XMUL, SU, U, XEL, YEL, SURFAC, IKLE1, IKLE2, IKLE3, NELEM, NELMAX)
Definition: mt02cc.f:13
declarations_special::lu
integer lu
Definition: declarations_special.F:45
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:36:45 for The TELEMAC-MASCARET system by doxygen 1.8.13