The TELEMAC-MASCARET system  trunk
mt03cc.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt03cc
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,sf,sg,su,sv,f,g,u,v,
12  & xel,yel,ikle1,ikle2,ikle3,ikle4,ikle5,ikle6,
13  & nelem,nelmax)
14 !
15 !***********************************************************************
16 ! BIEF V6P1 21/08/2010
17 !***********************************************************************
18 !
19 !brief COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
20 !code
21 !+ / --> - --> --> --->
22 !+ A(I,J) = XMUL / KEL . GRAD(PSI1(I)) * U . GRAD(PSI2(J)) D(OMEGA)
23 !+ /OMEGA
24 !+ -->
25 !+ KEL CONSTANT VECTOR ON THE ELEMENT, WITH COMPONENTS F AND G
26 !+
27 !+ PSI1 OF P2 DISCRETISATION
28 !+ PSI2 OF P2 DISCRETISATION
29 !
30 !warning THE JACOBIAN MUST BE POSITIVE
31 !
32 !history ALGIANE FROEHLY (MATMECA)
33 !+ 19/06/08
34 !+ V5P9
35 !+
36 !
37 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
38 !+ 13/07/2010
39 !+ V6P0
40 !+ Translation of French comments within the FORTRAN sources into
41 !+ English comments
42 !
43 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
44 !+ 21/08/2010
45 !+ V6P0
46 !+ Creation of DOXYGEN tags for automated documentation and
47 !+ cross-referencing of the FORTRAN sources
48 !
49 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
50 !| A11 |<--| ELEMENTS OF MATRIX
51 !| ... |<--| ELEMENTS OF MATRIX
52 !| A66 |<--| ELEMENTS OF MATRIX
53 !| F |-->| FUNCTION USED IN THE FORMULA
54 !| G |-->| FUNCTION USED IN THE FORMULA
55 !| IKLE1 |-->| FIRST POINTS OF TRIANGLES
56 !| IKLE2 |-->| SECOND POINTS OF TRIANGLES
57 !| IKLE3 |-->| THIRD POINTS OF TRIANGLES
58 !| IKLE4 |-->| FOURTH POINTS OF TRIANGLES (QUADRATIC)
59 !| IKLE5 |-->| FIFTH POINTS OF TRIANGLES (QUADRATIC)
60 !| IKLE6 |-->| SIXTH POINTS OF TRIANGLES (QUADRATIC)
61 !| NELEM |-->| NUMBER OF ELEMENTS
62 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
63 !| SF |-->| STRUCTURE OF FUNCTIONS F
64 !| SG |-->| STRUCTURE OF FUNCTIONS G
65 !| SU |-->| BIEF_OBJ STRUCTURE OF U
66 !| SV |-->| BIEF_OBJ STRUCTURE OF V
67 !| U |-->| FUNCTION U USED IN THE FORMULA
68 !| V |-->| FUNCTION V USED IN THE FORMULA
69 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
70 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
71 !| XMUL |-->| MULTIPLICATION FACTOR
72 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
73 !
74  USE bief!, EX_MT03CC => MT03CC
75 !
76  USE declarations_special
77  IMPLICIT NONE
78 !
79 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
80 !
81  INTEGER, INTENT(IN) :: NELEM,NELMAX
82  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax)
83  INTEGER, INTENT(IN) :: IKLE3(nelmax),IKLE4(nelmax)
84  INTEGER, INTENT(IN) :: IKLE5(nelmax),IKLE6(nelmax)
85 !
86  DOUBLE PRECISION, INTENT(INOUT) :: A11(*),A12(*),A13(*)
87  DOUBLE PRECISION, INTENT(INOUT) :: A14(*),A15(*),A16(*)
88  DOUBLE PRECISION, INTENT(INOUT) :: A21(*),A22(*),A23(*)
89  DOUBLE PRECISION, INTENT(INOUT) :: A24(*),A25(*),A26(*)
90  DOUBLE PRECISION, INTENT(INOUT) :: A31(*),A32(*),A33(*)
91  DOUBLE PRECISION, INTENT(INOUT) :: A34(*),A35(*),A36(*)
92  DOUBLE PRECISION, INTENT(INOUT) :: A41(*),A42(*),A43(*)
93  DOUBLE PRECISION, INTENT(INOUT) :: A44(*),A45(*),A46(*)
94  DOUBLE PRECISION, INTENT(INOUT) :: A51(*),A52(*),A53(*)
95  DOUBLE PRECISION, INTENT(INOUT) :: A54(*),A55(*),A56(*)
96  DOUBLE PRECISION, INTENT(INOUT) :: A61(*),A62(*),A63(*)
97  DOUBLE PRECISION, INTENT(INOUT) :: A64(*),A65(*),A66(*)
98 !
99  DOUBLE PRECISION, INTENT(IN) :: XMUL
100  DOUBLE PRECISION, INTENT(IN) :: F(*),G(*),U(*),V(*)
101 !
102 ! STRUCTURES OF F, G, U, V
103  TYPE(bief_obj), INTENT(IN) :: SF,SG,SU,SV
104 !
105  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,3),YEL(nelmax,3)
106 !
107 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
108 !
109 ! DECLARATIONS SPECIFIC TO THIS SUBROUTINE
110 !
111  INTEGER IELEM,IELMF,IELMG,IELMU,IELMV
112 !
113  DOUBLE PRECISION X2,X3,Y2,Y3,U1,U2,U3,U4,U5,U6,V1,V2,V3,V4,V5,V6
114  DOUBLE PRECISION KXEL,KYEL,ANS1,ANS2
115 !
116 !-----------------------------------------------------------------------
117 !
118  ielmf = sf%ELM
119  ielmg = sg%ELM
120  ielmu = su%ELM
121  ielmv = sv%ELM
122 !
123 !-----------------------------------------------------------------------
124 ! CASE WHERE U IS OF P1 DISCRETISATION
125 !-----------------------------------------------------------------------
126 !
127  IF(ielmf.EQ.10.AND.ielmg.EQ.10.AND.
128  & ielmu.EQ.11.AND.ielmv.EQ.11) THEN
129 !
130 ! LOOP ON THE ELEMENTS
131 !
132  DO ielem = 1 , nelem
133 !
134  x2 = xel(ielem,2)
135  x3 = xel(ielem,3)
136  y2 = yel(ielem,2)
137  y3 = yel(ielem,3)
138 !
139  u1 = u(ikle1(ielem))
140  u2 = u(ikle2(ielem))
141  u3 = u(ikle3(ielem))
142  v1 = v(ikle1(ielem))
143  v2 = v(ikle2(ielem))
144  v3 = v(ikle3(ielem))
145 !
146  kxel = f(ielem)
147  kyel = g(ielem)
148 !
149 ! COMPUTES 30 OF THE 36 TERMS (SELECTED AMONG THE LEAST COMPLEX)
150 !
151  a11(ielem)= ( kxel * (y3-y2) - kyel * (x3-x2) ) *
152  & ((x2-x3)*(v3+3.d0*v1+v2) + (y3-y2)*(u2+3.d0*u1+u3)) *
153  & xmul/((x2*y3-y2*x3)*10.d0)
154 !
155  a12(ielem)= ( kxel * (y3-y2) - kyel * (x3-x2) ) *
156  & ( x3 * (v3+2.d0*(v1+v2)) - y3 * (u3+2.d0*(u1+u2))) *
157  & xmul/((x3*y2-y3*x2)*30.d0)
158 !
159  a13(ielem)= ( kxel * (y3-y2) - kyel * (x3-x2) ) *
160  & ( x2 * (v2+2.d0*(v1+v3)) - y2 * (u2+2.d0*(u1+u3))) *
161  & xmul/((x2*y3-y2*x3)*30.d0)
162 !
163  a14(ielem)= ( kxel * (y3-y2) - kyel * (x3-x2) ) *
164  & ( x2 *( - v3+ 3.d0*v1- 2.d0*v2) +
165  & x3 *(4.d0*v3+11.d0*v1+ 5.d0*v2) +
166  & y2 *( u3+ 2.d0*u2- 3.d0*u1) -
167  & y3 *(11.d0*u1+5.d0*u2+ 4.d0*u3))*
168  & xmul/((x2*y3-y2*x3)*30.d0)
169 !
170  a15(ielem)= ( kxel * (y3-y2) - kyel * (x3-x2) ) *
171  & ( x3 *( v2+2.d0*v3- 3.d0*v1) +
172  & x2 *(3.d0*v1- v3- 2.d0*v2) +
173  & y3 *(3.d0*u1- u2- 2.d0*u3) +
174  & y2 *(-3.d0*u1+2.d0*u2+ u3))*
175  & xmul/((x3*y2-y3*x2)*30.d0)
176 !
177  a21(ielem)= ( kyel * x3 - kxel * y3 ) *
178  & ( (x2-x3) * ( v3 + 2.d0*(v1+v2) ) +
179  & (y3-y2) * ( u3 + 2.d0*(u1+u2) )) *
180  & xmul/((x3*y2-y3*x2)*30.d0)
181 !
182  a22(ielem)= ( kyel * x3 - kxel * y3 ) *
183  & ( x3 * (v3+v1+3.d0*v2) - y3*(u3+u1+3.d0*u2)) *
184  & xmul/((x2*y3-y2*x3)*10.d0)
185 !
186  a23(ielem)= ( kyel * x3 - kxel * y3 ) *
187  & ( x2 * (v1+2.d0*(v2+v3)) - y2 * (u1+2.d0*(u2+u3))) *
188  & xmul/((x2*y3-y2*x3)*30.d0)
189 !
190  a24(ielem)= ( kyel * x3 - kxel * y3 ) *
191  & ( -x3*(4.d0*v3+ 5.d0*v1+11.d0*v2) +
192  & x2*(3.d0*v1+14.d0*v2+ 3.d0*v3) +
193  & y3*(5.d0*u1+11.d0*u2+ 4.d0*u3) -
194  & y2*(3.d0*u1+14.d0*u2+ 3.d0*u3))*
195  & xmul/((x2*y3-y2*x3)*30.d0)
196 !
197  a26(ielem)= ( kxel * y3 - kyel * x3 ) *
198  & (x3*(2.d0*v3 + v1 - 3.d0*v2) +
199  & x2*(v1-v3) + y2*(u3-u1) +
200  & y3*(3.d0*u2 - u1 - 2.d0*u3))*
201  & xmul/((x3*y2-y3*x2)*30.d0)
202 !
203  a31(ielem)= ( kyel * x2 - kxel * y2 ) *
204  & ( (x2-x3) * (2.d0*(v1+v3) + v2 ) +
205  & (y3-y2) * (2.d0*(u1+u3) + u2 )) *
206  & xmul/((x2*y3-y2*x3)*30.d0)
207 !
208  a32(ielem)= ( kyel * x2 - kxel * y2 ) *
209  & ( x3*(2.d0*(v3+v2)+v1)-y3*(u1+2.d0*(u2+u3))) *
210  & xmul/((x2*y3-y2*x3)*30.d0)
211 !
212  a33(ielem)= ( kyel*x2 - kxel * y2 ) *
213  & (x2*(v1+3.d0*v3+v2)-y2*(u1+3.d0*u3+u2)) *
214  & xmul/(10.d0*(x2*y3-y2*x3))
215 !
216  a34(ielem)= ( kyel*x2 - kxel * y2 ) *
217  & (x3*(v1-v2)+x2*(v1+2.d0*v2-3.d0*v3) +
218  & y3*(u2-u1)+y2*(3.d0*u3-u1-2.d0*u2))*
219  & xmul/((x2*y3-y2*x3)*30.d0)
220 !
221  a35(ielem)= ( kxel*y2 - kyel * x2 ) *
222  & (x3*(14.d0*v3+3.d0*(v2+v1)) +
223  & x2*( v1-3.d0*v3+2.d0*v2) -
224  & y3*(3.d0*(u1+u2)+14.d0*u3) +
225  & y2*( 3.d0*u3- u1- 2.d0*u2 ))*
226  & xmul/((x2*y3-y2*x3)*30.d0)
227 !
228  a36(ielem)= ( kyel*x2 - kxel * y2 ) *
229  & ( x3*(-14.d0*v3-3.d0*(v1+v2)) +
230  & x2*( 5.d0*v1+4.d0*v2+11.d0*v3) +
231  & y3*(3.d0*(u1+u2)+14.d0*u3) -
232  & y2*(5.d0*u1+4.d0*u2+11.d0*u3)) *
233  & xmul/((x3*y2-y3*x2)*30.d0)
234 !
235  a42(ielem)=(-5.d0*kxel*y3**2*u1+11.d0*kyel*x3*y3*u2+
236  & 5.d0*kyel*x3*y3*u1+3.d0*kyel*x3*x2*v1+
237  & 5.d0*kxel*y3*x3*v1+3.d0*kxel*y3*y2*u1-
238  & 5.d0*kyel*x3**2*v1-3.d0*kxel*y2*x3*v3-
239  & 3.d0*kyel*x2*y3*u3+4.d0*kyel*x3*y3*u3+
240  & 3.d0*kyel*x3*x2*v3+3.d0*kxel*y3*y2*u3+
241  & 4.d0*kxel*y3*x3*v3-4.d0*kxel*y3**2*u3-
242  & 4.d0*kyel*x3**2*v3-14.d0*kxel*y2*x3*v2-
243  & 14.d0*kyel*x2*y3*u2+14.d0*kxel*y3*y2*u2+
244  & 11.d0*kxel*y3*x3*v2+14.d0*kyel*x3*x2*v2-
245  & 11.d0*kyel*x3**2*v2-11.d0*kxel*y3**2*u2-
246  & 3.d0*kyel*x2*y3*u1-3.d0*kxel*y2*x3*v1)*
247  & xmul/((x2*y3-y2*x3)*30.d0)
248 !
249  a43(ielem)=(kyel*x3*x2*v1-kyel*x3*y2*u1-kxel*y3*x2*v1+
250  & kxel*y3*y2*u1+kyel*x2**2*v1+3.d0*kyel*x2*y2*u3+
251  & 3.d0*kxel*y2*x2*v3-3.d0*kyel*x2**2*v3-
252  & 3.d0*kxel*y2**2*u3-2.d0*kxel*y2*x2*v2+
253  & 2.d0*kyel*x2**2*v2+2.d0*kxel*y2**2*u2-
254  & 2.d0*kyel*x2*y2*u2-kxel*y3*y2*u2+kxel*y3*x2*v2+
255  & kyel*x3*y2*u2-kyel*x3*x2*v2+kxel*y2**2*u1-
256  & kyel*x2*y2*u1-kxel*y2*x2*v1)*
257  & xmul/((x2*y3-y2*x3)*30.d0)
258 !
259  a44(ielem)=(-4.d0*kyel*x3*y3*u2-4.d0*kyel*x3*y3*u1-
260  & 4.d0*kxel*y3*x3*v1-2.d0*kyel*x2*y2*u3+
261  & kxel*y2*x3*v3-2.d0*kxel*y2*x2*v3+kyel*x2*y3*u3+
262  & kyel*x3*y2*u3-2.d0*kyel*x3*x2*v3-
263  & 2.d0*kxel*y3*y2*u3+kxel*y3*x2*v3+
264  & 4.d0*kxel*y2*x3*v2-6.d0*kxel*y2*x2*v2+
265  & 4.d0*kyel*x2*y3*u2+2.d0*kyel*x2**2*v1+
266  & 4.d0*kyel*x3**2*v1+2.d0*kyel*x2**2*v3+
267  & 2.d0*kxel*y2**2*u3+6.d0*kyel*x2**2*v2+
268  & 6.d0*kxel*y2**2*u2+4.d0*kxel*y3**2*u1-
269  & 6.d0*kyel*x2*y2*u2-8.d0*kxel*y3*y2*u2-
270  & 4.d0*kxel*y3*x3*v2+4.d0*kxel*y3*x2*v2+
271  & 4.d0*kyel*x3*y2*u2-8.d0*kyel*x3*x2*v2-
272  & 2.d0*kyel*x2*y2*u1-2.d0*kxel*y2*x2*v1+
273  & 4.d0*kyel*x3**2*v2+2.d0*kxel*y2**2*u1+
274  & 4.d0*kxel*y3**2*u2+2.d0*kxel*y3**2*u3+
275  & 2.d0*kyel*x3**2*v3-2.d0*kyel*x3*y3*u3-
276  & 2.d0*kxel*y3*x3*v3)*
277  & xmul/(15.d0*(x2*y3-y2*x3)/2.d0)
278 !
279  a45(ielem)=(-kyel*x3*y3*u2+kyel*x3*y3*u1-kyel*x3*x2*v1+
280  & kxel*y3*x3*v1-kxel*y3*y2*u1-2.d0*kyel*x2*y2*u3+
281  & 2.d0*kxel*y2*x3*v3-2.d0*kxel*y2*x2*v3+
282  & 2.d0*kyel*x2*y3*u3+kyel*x3*y2*u3-3.d0*kyel*x3*x2*v3-
283  & 3.d0*kxel*y3*y2*u3+kxel*y3*x2*v3+2.d0*kxel*y2*x3*v2-
284  & 6.d0*kxel*y2*x2*v2+2.d0*kyel*x2*y3*u2+
285  & 2.d0*kyel*x2**2*v1-kyel*x3**2*v1+2.d0*kyel*x2**2*v3+
286  & 2.d0*kxel*y2**2*u3+6.d0*kyel*x2**2*v2+
287  & 6.d0*kxel*y2**2*u2-kxel*y3**2*u1-6.d0*kyel*x2*y2*u2-
288  & 6.d0*kxel*y3*y2*u2-kxel*y3*x3*v2+4.d0*kxel*y3*x2*v2+
289  & 4.d0*kyel*x3*y2*u2-6.d0*kyel*x3*x2*v2-
290  & 2.d0*kyel*x2*y2*u1-2.d0*kxel*y2*x2*v1+kyel*x3**2*v2+
291  & 2.d0*kxel*y2**2*u1+kxel*y3**2*u2+kyel*x2*y3*u1+
292  & kxel*y2*x3*v1)*xmul/(15.d0*(x3*y2-y3*x2)/2.d0)
293 !
294  a46(ielem)=(kyel*x3*y3*u2-kyel*x3*y3*u1+4.d0*kyel*x3*x2*v1-
295  & 3.d0*kyel*x3*y2*u1-3.d0*kxel*y3*x2*v1-kxel*y3*x3*v1+
296  & 4.d0*kxel*y3*y2*u1+kyel*x2*y2*u3-
297  & 2.d0*kxel*y2*x3*v3+kxel*y2*x2*v3-
298  & 2.d0*kyel*x2*y3*u3-kyel*x3*y2*u3+
299  & 3.d0*kyel*x3*x2*v3+3.d0*kxel*y3*y2*u3-kxel*y3*x2*v3-
300  & 2.d0*kxel*y2*x3*v2-2.d0*kyel*x2*y3*u2+kyel*x2**2*v1+
301  & kyel*x3**2*v1-kyel*x2**2*v3-kxel*y2**2*u3+
302  & kxel*y3**2*u1+3.d0*kxel*y3*y2*u2+kxel*y3*x3*v2-
303  & kxel*y3*x2*v2-kyel*x3*y2*u2+3.d0*kyel*x3*x2*v2-
304  & kyel*x2*y2*u1-kxel*y2*x2*v1-kyel*x3**2*v2+
305  & kxel*y2**2*u1-kxel*y3**2*u2-kyel*x2*y3*u1-
306  & kxel*y2*x3*v1)*xmul/(15.d0*(x3*y2-y3*x2)/2.d0)
307 !
308  a52(ielem)=(kxel*y3**2*u1+3.d0*kyel*x3*y3*u2-kyel*x3*y3*u1+
309  & 3.d0*kyel*x3*x2*v1-kxel*y3*x3*v1+3.d0*kxel*y3*y2*u1+
310  & kyel*x3**2*v1-3.d0*kxel*y2*x3*v3-3.d0*kyel*x2*y3*u3-
311  & 2.d0*kyel*x3*y3*u3+3.d0*kyel*x3*x2*v3+
312  & 3.d0*kxel*y3*y2*u3-2.d0*kxel*y3*x3*v3+
313  & 2.d0*kxel*y3**2*u3+2.d0*kyel*x3**2*v3-
314  & 14.d0*kxel*y2*x3*v2-14.d0*kyel*x2*y3*u2+
315  & 14.d0*kxel*y3*y2*u2+3.d0*kxel*y3*x3*v2+
316  & 14.d0*kyel*x3*x2*v2-3.d0*kyel*x3**2*v2-
317  & 3.d0*kxel*y3**2*u2-3.d0*kyel*x2*y3*u1-
318  & 3.d0*kxel*y2*x3*v1)*xmul/((x3*y2-y3*x2)*30.d0)
319 !
320  a53(ielem)=(3.d0*kyel*x3*x2*v1-3.d0*kyel*x3*y2*u1-
321  & 3.d0*kxel*y3*x2*v1+3.d0*kxel*y3*y2*u1+
322  & kyel*x2**2*v1+3.d0*kyel*x2*y2*u3+3.d0*kxel*y2*x2*v3-
323  & 3.d0*kyel*x2**2*v3-3.d0*kxel*y2**2*u3-
324  & 14.d0*kyel*x3*y2*u3+14.d0*kyel*x3*x2*v3+
325  & 14.d0*kxel*y3*y2*u3-14.d0*kxel*y3*x2*v3-
326  & 2.d0*kxel*y2*x2*v2+2.d0*kyel*x2**2*v2+
327  & 2.d0*kxel*y2**2*u2-2.d0*kyel*x2*y2*u2+
328  & 3.d0*kxel*y3*y2*u2-3.d0*kxel*y3*x2*v2-
329  & 3.d0*kyel*x3*y2*u2+3.d0*kyel*x3*x2*v2+
330  & kxel*y2**2*u1-kyel*x2*y2*u1-kxel*y2*x2*v1)*
331  & xmul/((x3*y2-y3*x2)*30.d0)
332 !
333  a54(ielem)=(-kyel*x3*y3*u2+kyel*x3*y3*u1-kyel*x3*x2*v1+
334  & kyel*x3*y2*u1+kxel*y3*x2*v1+kxel*y3*x3*v1-
335  & kxel*y3*y2*u1-2.d0*kyel*x2*y2*u3+kxel*y2*x3*v3-
336  & 2.d0*kxel*y2*x2*v3+kyel*x2*y3*u3+2.d0*kyel*x3*y2*u3-
337  & 3.d0*kyel*x3*x2*v3-3.d0*kxel*y3*y2*u3+
338  & 2.d0*kxel*y3*x2*v3+4.d0*kxel*y2*x3*v2-
339  & 6.d0*kxel*y2*x2*v2+4.d0*kyel*x2*y3*u2+
340  & 2.d0*kyel*x2**2*v1-kyel*x3**2*v1+2.d0*kyel*x2**2*v3+
341  & 2.d0*kxel*y2**2*u3+6.d0*kyel*x2**2*v2+
342  & 6.d0*kxel*y2**2*u2-kxel*y3**2*u1-6.d0*kyel*x2*y2*u2-
343  & 6.d0*kxel*y3*y2*u2-kxel*y3*x3*v2+2.d0*kxel*y3*x2*v2+
344  & 2.d0*kyel*x3*y2*u2-6.d0*kyel*x3*x2*v2-
345  & 2.d0*kyel*x2*y2*u1-2.d0*kxel*y2*x2*v1+
346  & kyel*x3**2*v2+2.d0*kxel*y2**2*u1+kxel*y3**2*u2)
347  & *xmul/(15.d0*(x3*y2-y3*x2)/2.d0)
348 !
349  a56(ielem)=(2.d0*kyel*x3*y3*u2+2.d0*kyel*x3*y3*u1+kyel*x3*x2*v1+
350  & 2.d0*kxel*y3*x3*v1+kxel*y3*y2*u1+kyel*x2*y2*u3-
351  & 2.d0*kxel*y2*x3*v3+kxel*y2*x2*v3-2.d0*kyel*x2*y3*u3-
352  & 4.d0*kyel*x3*y2*u3+6.d0*kyel*x3*x2*v3+
353  & 6.d0*kxel*y3*y2*u3-4.d0*kxel*y3*x2*v3-
354  & 2.d0*kxel*y2*x3*v2-2.d0*kyel*x2*y3*u2+kyel*x2**2*v1-
355  & 2.d0*kyel*x3**2*v1-kyel*x2**2*v3-kxel*y2**2*u3-
356  & 2.d0*kxel*y3**2*u1+3.d0*kxel*y3*y2*u2+
357  & 2.d0*kxel*y3*x3*v2-kxel*y3*x2*v2-kyel*x3*y2*u2+
358  & 3.d0*kyel*x3*x2*v2-kyel*x2*y2*u1-kxel*y2*x2*v1-
359  & 2.d0*kyel*x3**2*v2+kxel*y2**2*u1-2.d0*kxel*y3**2*u2-
360  & 6.d0*kxel*y3**2*u3-6.d0*kyel*x3**2*v3+
361  & 6.d0*kyel*x3*y3*u3+6.d0*kxel*y3*x3*v3-
362  & kyel*x2*y3*u1-kxel*y2*x3*v1)*
363  & xmul/(15.d0*(x2*y3-y2*x3)/2.d0)
364 !
365  a62(ielem)=(kxel*y3**2*u1+3.d0*kyel*x3*y3*u2-kyel*x3*y3*u1+
366  & kyel*x3*x2*v1-kxel*y3*x3*v1+kxel*y3*y2*u1+
367  & kyel*x3**2*v1+kxel*y2*x3*v3+kyel*x2*y3*u3-
368  & 2.d0*kyel*x3*y3*u3-kyel*x3*x2*v3-kxel*y3*y2*u3-
369  & 2.d0*kxel*y3*x3*v3+2.d0*kxel*y3**2*u3+
370  & 2.d0*kyel*x3**2*v3+3.d0*kxel*y3*x3*v2-
371  & 3.d0*kyel*x3**2*v2-3.d0*kxel*y3**2*u2-
372  & kyel*x2*y3*u1-kxel*y2*x3*v1)
373  & *xmul/((x2*y3-y2*x3)*30.d0)
374 !
375  a63(ielem)=(-3.d0*kyel*x3*x2*v1+3.d0*kyel*x3*y2*u1+
376  & 3.d0*kxel*y3*x2*v1-3.d0*kxel*y3*y2*u1+
377  & 5.d0*kyel*x2**2*v1-11.d0*kyel*x2*y2*u3-
378  & 11.d0*kxel*y2*x2*v3+11.d0*kyel*x2**2*v3+
379  & 11.d0*kxel*y2**2*u3+14.d0*kyel*x3*y2*u3-
380  & 14.d0*kyel*x3*x2*v3-14.d0*kxel*y3*y2*u3+
381  & 14.d0*kxel*y3*x2*v3-4.d0*kxel*y2*x2*v2+
382  & 4.d0*kyel*x2**2*v2+4.d0*kxel*y2**2*u2-
383  & 4.d0*kyel*x2*y2*u2-3.d0*kxel*y3*y2*u2+
384  & 3.d0*kxel*y3*x2*v2+3.d0*kyel*x3*y2*u2-
385  & 3.d0*kyel*x3*x2*v2+5.d0*kxel*y2**2*u1-
386  & 5.d0*kyel*x2*y2*u1-5.d0*kxel*y2*x2*v1)*
387  & xmul/((x3*y2-y3*x2)*30.d0)
388 !
389  a64(ielem)=(kyel*x3*y3*u2-kyel*x3*y3*u1+4.d0*kyel*x3*x2*v1-
390  & kyel*x3*y2*u1-kxel*y3*x2*v1-kxel*y3*x3*v1+
391  & 4.d0*kxel*y3*y2*u1+kyel*x2*y2*u3-kxel*y2*x3*v3+
392  & kxel*y2*x2*v3-kyel*x2*y3*u3-2.d0*kyel*x3*y2*u3+
393  & 3.d0*kyel*x3*x2*v3+3.d0*kxel*y3*y2*u3-
394  & 2.d0*kxel*y3*x2*v3-kxel*y2*x3*v2-kyel*x2*y3*u2+
395  & kyel*x2**2*v1+kyel*x3**2*v1-kyel*x2**2*v3-
396  & kxel*y2**2*u3+kxel*y3**2*u1+3.d0*kxel*y3*y2*u2+
397  & kxel*y3*x3*v2-2.d0*kxel*y3*x2*v2-
398  & 2.d0*kyel*x3*y2*u2+3.d0*kyel*x3*x2*v2-
399  & kyel*x2*y2*u1-kxel*y2*x2*v1-kyel*x3**2*v2+
400  & kxel*y2**2*u1-kxel*y3**2*u2-3.d0*kyel*x2*y3*u1-
401  & 3.d0*kxel*y2*x3*v1)*
402  & xmul/(15.d0*(x3*y2-y3*x2)/2.d0)
403 !
404  a65(ielem)=(2.d0*kyel*x3*y3*u2+2.d0*kyel*x3*y3*u1+kyel*x3*x2*v1-
405  & kyel*x3*y2*u1-kxel*y3*x2*v1+2.d0*kxel*y3*x3*v1+
406  & kxel*y3*y2*u1+kyel*x2*y2*u3-4.d0*kxel*y2*x3*v3+
407  & kxel*y2*x2*v3-4.d0*kyel*x2*y3*u3-2.d0*kyel*x3*y2*u3+
408  & 6.d0*kyel*x3*x2*v3+6.d0*kxel*y3*y2*u3-
409  & 2.d0*kxel*y3*x2*v3-kxel*y2*x3*v2-kyel*x2*y3*u2+
410  & kyel*x2**2*v1-2.d0*kyel*x3**2*v1-kyel*x2**2*v3-
411  & kxel*y2**2*u3-2.d0*kxel*y3**2*u1+3.d0*kxel*y3*y2*u2+
412  & 2.d0*kxel*y3*x3*v2-2.d0*kxel*y3*x2*v2-
413  & 2.d0*kyel*x3*y2*u2+3.d0*kyel*x3*x2*v2-kyel*x2*y2*u1-
414  & kxel*y2*x2*v1-2.d0*kyel*x3**2*v2+kxel*y2**2*u1-
415  & 2.d0*kxel*y3**2*u2-6.d0*kxel*y3**2*u3-
416  & 6.d0*kyel*x3**2*v3+6.d0*kyel*x3*y3*u3+
417  & 6.d0*kxel*y3*x3*v3)*xmul/(15.d0*(x2*y3-y2*x3)/2.d0)
418 !
419  a66(ielem)=(-2.d0*kyel*x3*y3*u2-2.d0*kyel*x3*y3*u1-
420  & 2.d0*kxel*y3*x3*v1-4.d0*kyel*x2*y2*u3+
421  & 4.d0*kxel*y2*x3*v3-4.d0*kxel*y2*x2*v3+
422  & 4.d0*kyel*x2*y3*u3+4.d0*kyel*x3*y2*u3-
423  & 8.d0*kyel*x3*x2*v3-8.d0*kxel*y3*y2*u3+
424  & 4.d0*kxel*y3*x2*v3+kxel*y2*x3*v2-
425  & 2.d0*kxel*y2*x2*v2+kyel*x2*y3*u2+
426  & 4.d0*kyel*x2**2*v1+2.d0*kyel*x3**2*v1+
427  & 4.d0*kyel*x2**2*v3+4.d0*kxel*y2**2*u3+
428  & 2.d0*kyel*x2**2*v2+2.d0*kxel*y2**2*u2+
429  & 2.d0*kxel*y3**2*u1-2.d0*kyel*x2*y2*u2-
430  & 2.d0*kxel*y3*y2*u2-2.d0*kxel*y3*x3*v2+
431  & kxel*y3*x2*v2+kyel*x3*y2*u2-2.d0*kyel*x3*x2*v2-
432  & 4.d0*kyel*x2*y2*u1-4.d0*kxel*y2*x2*v1+
433  & 2.d0*kyel*x3**2*v2+4.d0*kxel*y2**2*u1+
434  & 2.d0*kxel*y3**2*u2+6.d0*kxel*y3**2*u3+
435  & 6.d0*kyel*x3**2*v3-6.d0*kyel*x3*y3*u3-
436  & 6.d0*kxel*y3*x3*v3)*xmul/(15.d0*(x2*y3-y2*x3)/2.d0)
437 !
438 ! USES HERE THE 'MAGIC SQUARE' PROPERTIES
439 ! (SUM OF EACH LINE = SUM OF EACH COLUMN = 0)
440 !
441  a16(ielem) = - a11(ielem) - a12(ielem) - a13(ielem)
442  & - a14(ielem) - a15(ielem)
443  a25(ielem) = - a21(ielem) - a22(ielem) - a23(ielem)
444  & - a24(ielem) - a26(ielem)
445  a41(ielem) = - a42(ielem) - a43(ielem) - a44(ielem)
446  & - a45(ielem) - a46(ielem)
447  a61(ielem) = - a62(ielem) - a63(ielem) - a64(ielem)
448  & - a65(ielem) - a66(ielem)
449  a51(ielem) = - a11(ielem) - a21(ielem) - a31(ielem)
450  & - a41(ielem) - a61(ielem)
451  a55(ielem) = - a51(ielem) - a52(ielem) - a53(ielem)
452  & - a54(ielem) - a56(ielem)
453 !
454  ENDDO ! IELEM
455 !
456 !-----------------------------------------------------------------------
457 ! CASE WHERE U IS OF P2 DISCRETISATION
458 !-----------------------------------------------------------------------
459 !
460  ELSEIF(ielmf.EQ.10.AND.ielmg.EQ.10.AND.
461  & ielmu.EQ.13.AND.ielmv.EQ.13) THEN
462 !
463 ! LOOP ON THE ELEMENTS
464 !
465  DO ielem = 1 , nelem
466 !
467  x2 = xel(ielem,2)
468  x3 = xel(ielem,3)
469  y2 = yel(ielem,2)
470  y3 = yel(ielem,3)
471 !
472  u1 = u(ikle1(ielem))
473  u2 = u(ikle2(ielem))
474  u3 = u(ikle3(ielem))
475  u4 = u(ikle4(ielem))
476  u5 = u(ikle5(ielem))
477  u6 = u(ikle6(ielem))
478  v1 = v(ikle1(ielem))
479  v2 = v(ikle2(ielem))
480  v3 = v(ikle3(ielem))
481  v4 = v(ikle4(ielem))
482  v5 = v(ikle5(ielem))
483  v6 = v(ikle6(ielem))
484 !
485  kxel = f(ielem)
486  kyel = g(ielem)
487 !
488 ! COMPUTES 30 OF THE 36 TERMS (SELECTED AMONG THE LEAST COMPLEX)
489 !
490  a11(ielem)=(kxel*(y2-y3)+kyel*(x3-x2))*(12.d0*x3*v1+
491  & 12.d0*y2*u1-12.d0*y3*u1-12.d0*x2*v1+2.d0*y3*u3+
492  & 2.d0*x2*v3+15.d0*x3*v6-2.d0*y2*u3+15.d0*y2*u6-
493  & 15.d0*x2*v6-15.d0*y3*u6-2.d0*x3*v3+15.d0*y2*u4-
494  & 7.d0*y3*u5+15.d0*x3*v4-15.d0*x2*v4-15.d0*y3*u4+
495  & 7.d0*y2*u5+7.d0*x3*v5-7.d0*x2*v5-2.d0*y2*u2+
496  & 2.d0*x2*v2-2.d0*x3*v2+2.d0*y3*u2)*
497  & xmul/(x2*y3-y2*x3)/90.d0
498 !
499  a12(ielem)=(kxel*(y3-y2)+kyel*(x2-x3))*(3.d0*x3*v1-
500  & 3.d0*y3*u1+2.d0*y3*u3+5.d0*x3*v6-5.d0*y3*u6-
501  & 2.d0*x3*v3-5.d0*y3*u5+x3*v4-y3*u4+5.d0*x3*v5+
502  & 3.d0*x3*v2-3.d0*y3*u2)*xmul/(-x2*y3+y2*x3)/90.d0
503 !
504  a13(ielem)=(kxel*(y3-y2)+kyel*(x2-x3))*(3.d0*y2*u1-
505  & 3.d0*x2*v1-3.d0*x2*v3+3.d0*y2*u3+y2*u6-x2*v6+
506  & 5.d0*y2*u4-5.d0*x2*v4+5.d0*y2*u5-5.d0*x2*v5-
507  & 2.d0*y2*u2+2.d0*x2*v2)*xmul/(-x2*y3+y2*x3)/90.d0
508 !
509  a14(ielem)=(kxel*(y2-y3)+kyel*(x3-x2))*(15.d0*x3*v1-
510  & 3.d0*y2*u1-15.d0*y3*u1+3.d0*x2*v1+4.d0*y3*u3-
511  & x2*v3+20.d0*x3*v6+y2*u3-4.d0*y2*u6+4.d0*x2*v6-
512  & 20.d0*y3*u6-4.d0*x3*v3-8.d0*y2*u4-12.d0*y3*u5+
513  & 16.d0*x3*v4+8.d0*x2*v4-16.d0*y3*u4+8.d0*y2*u5+
514  & 12.d0*x3*v5-8.d0*x2*v5+6.d0*y2*u2-6.d0*x2*v2+
515  & x3*v2-y3*u2)*xmul/(-x2*y3+y2*x3)/90.d0
516 !
517  a15(ielem)=(kxel*(y2-y3)+kyel*(x3-x2))*(3.d0*x3*v1+
518  & 3.d0*y2*u1-3.d0*y3*u1-3.d0*x2*v1+6.d0*y3*u3+
519  & x2*v3+8.d0*x3*v6-y2*u3+4.d0*y2*u6-4.d0*x2*v6-
520  & 8.d0*y3*u6-6.d0*x3*v3+8.d0*y2*u4+8.d0*y3*u5+
521  & 4.d0*x3*v4-8.d0*x2*v4-4.d0*y3*u4-8.d0*y2*u5-
522  & 8.d0*x3*v5+8.d0*x2*v5-6.d0*y2*u2+6.d0*x2*v2-
523  & x3*v2+y3*u2)*xmul/(-x2*y3+y2*x3)/90.d0
524 !
525  a21(ielem)=(-kxel*y3+kyel*x3)*(3.d0*x3*v1+3.d0*y2*u1-
526  & 3.d0*y3*u1-3.d0*x2*v1+2.d0*y3*u3+2.d0*x2*v3+
527  & 5.d0*x3*v6-2.d0*y2*u3+5.d0*y2*u6-5.d0*x2*v6-
528  & 5.d0*y3*u6-2.d0*x3*v3+y2*u4-5.d0*y3*u5+
529  & x3*v4-x2*v4-y3*u4+5.d0*y2*u5+5.d0*x3*v5-
530  & 5.d0*x2*v5+3.d0*y2*u2-3.d0*x2*v2+3.d0*x3*v2-
531  & 3.d0*y3*u2)*xmul/(x2*y3-y2*x3)/90.d0
532 !
533  a22(ielem)=(-kxel*y3+kyel*x3)*(-2.d0*x3*v1+2.d0*y3*u1+
534  & 2.d0*y3*u3+7.d0*x3*v6-7.d0*y3*u6-2.d0*x3*v3-
535  & 15.d0*y3*u5+15.d0*x3*v4-15.d0*y3*u4+
536  & 15.d0*x3*v5+12.d0*x3*v2-12.d0*y3*u2)
537  & *xmul/(x2*y3-y2*x3)/90.d0
538 !
539  a23(ielem)=(-kxel*y3+kyel*x3)*(2.d0*y2*u1-2.d0*x2*v1+
540  & 3.d0*x2*v3-3.d0*y2*u3-5.d0*y2*u6+5.d0*x2*v6-
541  & 5.d0*y2*u4+5.d0*x2*v4-y2*u5+x2*v5-3.d0*y2*u2+
542  & 3.d0*x2*v2)*xmul/(x2*y3-y2*x3)/90.d0
543 !
544  a24(ielem)=(-kxel*y3+kyel*x3)*(x3*v1-5.d0*y2*u1-y3*u1+
545  & 5.d0*x2*v1+4.d0*y3*u3+5.d0*x2*v3+12.d0*x3*v6-
546  & 5.d0*y2*u3+4.d0*y2*u6-4.d0*x2*v6-12.d0*y3*u6-
547  & 4.d0*x3*v3+24.d0*y2*u4-20.d0*y3*u5+16.d0*x3*v4-
548  & 24.d0*x2*v4-16.d0*y3*u4+24.d0*y2*u5+20.d0*x3*v5-
549  & 24.d0*x2*v5+18.d0*y2*u2-18.d0*x2*v2+15.d0*x3*v2-
550  & 15.d0*y3*u2)*xmul/(-x2*y3+y2*x3)/90.d0
551 !
552  a26(ielem)=(-kxel*y3+kyel*x3)*(-x3*v1+5.d0*y2*u1+y3*u1-
553  & 5.d0*x2*v1+6.d0*y3*u3+5.d0*x2*v3-8.d0*x3*v6-
554  & 5.d0*y2*u3+8.d0*y3*u6-6.d0*x3*v3-4.d0*y2*u4-
555  & 8.d0*y3*u5+4.d0*x3*v4+4.d0*x2*v4-4.d0*y3*u4+
556  & 4.d0*y2*u5+8.d0*x3*v5-4.d0*x2*v5+3.d0*x3*v2-
557  & 3.d0*y3*u2)*xmul/(-x2*y3+y2*x3)/90.d0
558 !
559  a31(ielem)=(kxel*y2-kyel*x2)*(3.d0*x3*v1+3.d0*y2*u1-
560  & 3.d0*y3*u1-3.d0*x2*v1-3.d0*y3*u3-3.d0*x2*v3+
561  & x3*v6+3.d0*y2*u3+y2*u6-x2*v6-y3*u6+3.d0*x3*v3+
562  & 5.d0*y2*u4-5.d0*y3*u5+5.d0*x3*v4-5.d0*x2*v4-
563  & 5.d0*y3*u4+5.d0*y2*u5+5.d0*x3*v5-5.d0*x2*v5-
564  & 2.d0*y2*u2+2.d0*x2*v2-2.d0*x3*v2+2.d0*y3*u2)
565  & *xmul/(x2*y3-y2*x3)/90.d0
566 !
567  a32(ielem)=(kxel*y2-kyel*x2)*(-2.d0*x3*v1+2.d0*y3*u1-
568  & 3.d0*y3*u3+5.d0*x3*v6-5.d0*y3*u6+3.d0*x3*v3-
569  & y3*u5+5.d0*x3*v4-5.d0*y3*u4+x3*v5+3.d0*x3*v2-
570  & 3.d0*y3*u2)*xmul/(-x2*y3+y2*x3)/90.d0
571 !
572  a33(ielem)=(kxel*y2-kyel*x2)*(2.d0*y2*u1-2.d0*x2*v1+
573  & 12.d0*x2*v3-12.d0*y2*u3-15.d0*y2*u6+15.d0*x2*v6-
574  & 7.d0*y2*u4+7.d0*x2*v4-15.d0*y2*u5+15.d0*x2*v5+
575  & 2.d0*y2*u2-2.d0*x2*v2)*xmul/(-x2*y3+y2*x3)/90.d0
576 !
577  a34(ielem)=(kxel*y2-kyel*x2)*(5.d0*x3*v1-y2*u1-5.d0*y3*u1+
578  & x2*v1-3.d0*x2*v3-4.d0*x3*v6+3.d0*y2*u3+
579  & 4.d0*y2*u6-4.d0*x2*v6+4.d0*y3*u6-8.d0*y2*u4-
580  & 4.d0*y3*u5+8.d0*x2*v4+8.d0*y2*u5+4.d0*x3*v5-
581  & 8.d0*x2*v5-6.d0*y2*u2+6.d0*x2*v2-5.d0*x3*v2+
582  & 5.d0*y3*u2)*xmul/(-x2*y3+y2*x3)/90.d0
583 !
584  a35(ielem)=(kxel*y2-kyel*x2)*(-5.d0*x3*v1-y2*u1+5.d0*y3*u1+
585  & x2*v1-18.d0*y3*u3-3.d0*x2*v3+24.d0*x3*v6+
586  & 3.d0*y2*u3+4.d0*y2*u6-4.d0*x2*v6-24.d0*y3*u6+
587  & 18.d0*x3*v3-8.d0*y2*u4-24.d0*y3*u5+4.d0*x3*v4+
588  & 8.d0*x2*v4-4.d0*y3*u4+8.d0*y2*u5+24.d0*x3*v5-
589  & 8.d0*x2*v5-6.d0*y2*u2+6.d0*x2*v2-5.d0*x3*v2+
590  & 5.d0*y3*u2)*xmul/(x2*y3-y2*x3)/90.d0
591 !
592  ans1 = kyel*x3*y2*u2-6.d0*kyel*x2*y2*u2-6.d0*kxel*y2*x2*v2-
593  & 4.d0*kyel*x2**2*v6+15.d0*kxel*y3**2*u1-
594  & 12.d0*kxel*y3*y2*u1+15.d0*kxel*y3*x2*v1+
595  & 6.d0*kyel*x2**2*v2-kxel*y3*x3*v2+16.d0*kxel*y3*x2*v4+
596  & 6.d0*kxel*y2**2*u2+16.d0*kxel*y3**2*u4-
597  & 8.d0*kxel*y3*y2*u4-4.d0*kxel*y2**2*u6+
598  & 15.d0*kyel*x3*y2*u1-15.d0*kyel*x3*y3*u1-
599  & 12.d0*kyel*x3*x2*v1-16.d0*kxel*y3*x3*v4-
600  & 7.d0*kxel*y3*y2*u2+20.d0*kxel*y3**2*u6+
601  & 12.d0*kyel*x3**2*v5+15.d0*kyel*x3**2*v1-
602  & 15.d0*kxel*y3*x3*v1+8.d0*kyel*x2*y3*u5-
603  & 8.d0*kyel*x2*y2*u5+8.d0*kyel*x2*y2*u4+
604  & kxel*y3**2*u2+kyel*x3**2*v2+20.d0*kyel*x3**2*v6-
605  & 4.d0*kyel*x3**2*v3-4.d0*kxel*y3**2*u3-
606  & 8.d0*kxel*y2**2*u4+kyel*x2**2*v3+8.d0*kxel*y2**2*u5+
607  & 8.d0*kyel*x2**2*v5+3.d0*kxel*y2*x2*v1-
608  & 3.d0*kxel*y2*x3*v1-kxel*y2*x2*v3-3.d0*kxel*y2**2*u1+
609  & kxel*y2**2*u3-3.d0*kyel*x2**2*v1-8.d0*kxel*y2*x2*v5+
610  & 8.d0*kxel*y2*x3*v5+8.d0*kxel*y2*x2*v4-kyel*x3*y3*u2
611  ans2 = 4.d0*kyel*x3*y3*u3+3.d0*kyel*x3*x2*v3+
612  & 20.d0*kyel*x3*y2*u6-16.d0*kyel*x3*x2*v6-
613  & 20.d0*kyel*x3*y3*u6+3.d0*kxel*y3*y2*u3-
614  & 20.d0*kxel*y3*x3*v6+4.d0*kxel*y3*x3*v3-
615  & 16.d0*kxel*y3*y2*u6+12.d0*kxel*y3**2*u5+
616  & 16.d0*kyel*x3**2*v4-4.d0*kxel*y2*x3*v6+
617  & 3.d0*kyel*x2*y2*u1-3.d0*kyel*x2*y3*u1-kyel*x2*y2*u3+
618  & kyel*x2*y3*u3-4.d0*kyel*x2*y3*u6+kxel*y2*x3*v3
619  & -8.d0*kyel*x2**2*v4-7.d0*kyel*x3*x2*v2-
620  & 16.d0*kyel*x3*y3*u4+16.d0*kyel*x3*y2*u4+
621  & 12.d0*kxel*y3*x2*v5-20.d0*kxel*y3*y2*u5-
622  & 12.d0*kxel*y3*x3*v5-8.d0*kyel*x3*x2*v4+
623  & 20.d0*kxel*y3*x2*v6-20.d0*kyel*x3*x2*v5-
624  & 12.d0*kyel*x3*y3*u5+12.d0*kyel*x3*y2*u5+
625  & 6.d0*kyel*x2*y3*u2+6.d0*kxel*y2*x3*v2+kxel*y3*x2*v2-
626  & 8.d0*kyel*x2*y3*u4-8.d0*kxel*y2*x3*v4+
627  & 4.d0*kxel*y2*x2*v6-4.d0*kyel*x3*y2*u3+
628  & 4.d0*kyel*x2*y2*u6-4.d0*kxel*y3*x2*v3
629  a41(ielem)= (ans1+ans2)*xmul/(90.d0*(-x2*y3+y2*x3))
630 !
631  ans1 = kxel*y3**2*u1+5.d0*kxel*y3*y2*u1-15.d0*kxel*y3*x3*v2+
632  & 16.d0*kxel*y3**2*u4-24.d0*kxel*y3*y2*u4-kyel*x3*y3*u1+
633  & 5.d0*kyel*x3*x2*v1-16.d0*kxel*y3*x3*v4-
634  & 18.d0*kxel*y3*y2*u2+12.d0*kxel*y3**2*u6+
635  & 20.d0*kyel*x3**2*v5+kyel*x3**2*v1-kxel*y3*x3*v1+
636  & 24.d0*kyel*x2*y3*u5+15.d0*kxel*y3**2*u2+
637  & 15.d0*kyel*x3**2*v2+12.d0*kyel*x3**2*v6-
638  & 4.d0*kyel*x3**2*v3-4.d0*kxel*y3**2*u3-
639  & 5.d0*kxel*y2*x3*v1+24.d0*kxel*y2*x3*v5-
640  & 15.d0*kyel*x3*y3*u2+4.d0*kyel*x3*y3*u3
641  ans2 = 5.d0*kyel*x3*x2*v3-4.d0*kyel*x3*x2*v6-
642  & 12.d0*kyel*x3*y3*u6+5.d0*kxel*y3*y2*u3-
643  & 12.d0*kxel*y3*x3*v6+4.d0*kxel*y3*x3*v3-
644  & 4.d0*kxel*y3*y2*u6+20.d0*kxel*y3**2*u5+
645  & 16.d0*kyel*x3**2*v4+4.d0*kxel*y2*x3*v6-
646  & 5.d0*kyel*x2*y3*u1-5.d0*kyel*x2*y3*u3+
647  & 4.d0*kyel*x2*y3*u6-5.d0*kxel*y2*x3*v3-
648  & 18.d0*kyel*x3*x2*v2-16.d0*kyel*x3*y3*u4-
649  & 24.d0*kxel*y3*y2*u5-20.d0*kxel*y3*x3*v5-
650  & 24.d0*kyel*x3*x2*v4-24.d0*kyel*x3*x2*v5-
651  & 20.d0*kyel*x3*y3*u5+18.d0*kyel*x2*y3*u2+
652  & 18.d0*kxel*y2*x3*v2+24.d0*kyel*x2*y3*u4+
653  & 24.d0*kxel*y2*x3*v4
654  a42(ielem)= (ans1+ans2)*xmul/(90.d0*(-x2*y3+y2*x3))
655 !
656  ans1 = 6.d0*kyel*x3*y2*u2-6.d0*kyel*x2*y2*u2-
657  & 6.d0*kxel*y2*x2*v2+4.d0*kyel*x2**2*v6+
658  & 4.d0*kxel*y3**2.d0*u1+4.d0*kxel*y3*y2*u1-
659  & 2.d0*kxel*y3*x2*v1+6.d0*kyel*x2**2*v2-
660  & 4.d0*kxel*y3*x3*v2+4.d0*kxel*y3*x2*v4+
661  & 6.d0*kxel*y2**2*u2+8.d0*kxel*y3**2.d0*u4-
662  & 8.d0*kxel*y3*y2*u4+4.d0*kxel*y2**2*u6-
663  & 2.d0*kyel*x3*y2*u1-4.d0*kyel*x3*y3*u1+
664  & 4.d0*kyel*x3*x2*v1-8.d0*kxel*y3*x3*v4-
665  & 12.d0*kxel*y3*y2*u2+8.d0*kxel*y3**2*u6+
666  & 8.d0*kyel*x3**2*v5+4.d0*kyel*x3**2*v1-
667  & 4.d0*kxel*y3*x3*v1+8.d0*kyel*x2*y3*u5-
668  & 12.d0*kyel*x2*y2*u5-12.d0*kyel*x2*y2*u4+
669  & 4.d0*kxel*y3**2*u2+4.d0*kyel*x3**2*v2+
670  & 8.d0*kyel*x3**2*v6-2.d0*kyel*x3**2*v3-
671  & 2.d0*kxel*y3**2*u3+12.d0*kxel*y2**2*u4-
672  & 2.d0*kyel*x2**2*v3+12.d0*kxel*y2**2*u5+
673  & 12.d0*kyel*x2**2*v5+2.d0*kxel*y2*x2*v1-
674  & 2.d0*kxel*y2*x3*v1+2.d0*kxel*y2*x2*v3-
675  & 2.d0*kxel*y2**2*u1
676  ans2 = -2.d0*kxel*y2**2*u3-2.d0*kyel*x2**2*v1-
677  & 12.d0*kxel*y2*x2*v5+8.d0*kxel*y2*x3*v5-
678  & 12.d0*kxel*y2*x2*v4-4.d0*kyel*x3*y3*u2+
679  & 2.d0*kyel*x3*y3*u3+2.d0*kyel*x3*x2*v3-
680  & 8.d0*kyel*x3*y3*u6+2.d0*kxel*y3*y2*u3-
681  & 8.d0*kxel*y3*x3*v6+2.d0*kxel*y3*x3*v3+
682  & 8.d0*kxel*y3**2*u5+8.d0*kyel*x3**2*v4+
683  & 2.d0*kyel*x2*y2*u1-2.d0*kyel*x2*y3*u1+
684  & 2.d0*kyel*x2*y2*u3-kyel*x2*y3*u3-kxel*y2*x3*v3+
685  & 12.d0*kyel*x2**2*v4-12.d0*kyel*x3*x2*v2-
686  & 8.d0*kyel*x3*y3*u4+4.d0*kyel*x3*y2*u4+
687  & 8.d0*kxel*y3*x2*v5-16.d0*kxel*y3*y2*u5-
688  & 8.d0*kxel*y3*x3*v5-8.d0*kyel*x3*x2*v4-
689  & 16.d0*kyel*x3*x2*v5-8.d0*kyel*x3*y3*u5+
690  & 8.d0*kyel*x3*y2*u5+6.d0*kyel*x2*y3*u2+
691  & 6.d0*kxel*y2*x3*v2+6.d0*kxel*y3*x2*v2+
692  & 4.d0*kyel*x2*y3*u4+4.d0*kxel*y2*x3*v4-
693  & 4.d0*kxel*y2*x2*v6-kyel*x3*y2*u3-
694  & 4.d0*kyel*x2*y2*u6-kxel*y3*x2*v3
695  a44(ielem)= (ans1+ans2)*2.d0*xmul/(45.d0*(x2*y3-y2*x3))
696 !
697  ans1 = 6.d0*kyel*x3*y2*u2-6.d0*kyel*x2*y2*u2-
698  & 6.d0*kxel*y2*x2*v2+4.d0*kyel*x2**2*v6-
699  & kxel*y3**2*u1+3.d0*kxel*y3*y2*u1-2.d0*kxel*y3*x2*v1+
700  & 6.d0*kyel*x2**2*v2-kxel*y3*x3*v2+4.d0*kxel*y3*x2*v4+
701  & 6.d0*kxel*y2**2*u2-8.d0*kxel*y3*y2*u4+
702  & 4.d0*kxel*y2**2*u6-2.d0*kyel*x3*y2*u1+kyel*x3*y3*u1+
703  & 3.d0*kyel*x3*x2*v1-6.d0*kxel*y3*y2*u2-
704  & 4.d0*kxel*y3**2*u6+4.d0*kyel*x3**2*v5-
705  & kyel*x3**2*v1+kxel*y3*x3*v1+8.d0*kyel*x2*y3*u5-
706  & 12.d0*kyel*x2*y2*u5-12.d0*kyel*x2*y2*u4+
707  & kxel*y3**2*u2+kyel*x3**2*v2-4.d0*kyel*x3**2*v6+
708  & 12.d0*kxel*y2**2*u4-2.d0*kyel*x2**2*v3+
709  & 12.d0*kxel*y2**2*u5+12.d0*kyel*x2**2*v5+
710  & 2.d0*kxel*y2*x2*v1-kxel*y2*x3*v1+2.d0*kxel*y2*x2*v3-
711  & 2.d0*kxel*y2**2*u1
712  ans2 = -2.d0*kxel*y2**2*u3-2.d0*kyel*x2**2*v1-
713  & 12.d0*kxel*y2*x2*v5+8.d0*kxel*y2*x3*v5-
714  & 12.d0*kxel*y2*x2*v4-kyel*x3*y3*u2+
715  & kyel*x3*x2*v3-4.d0*kyel*x3*x2*v6+4.d0*kyel*x3*y3*u6+
716  & kxel*y3*y2*u3+4.d0*kxel*y3*x3*v6-4.d0*kxel*y3*y2*u6+
717  & 4.d0*kxel*y3**2*u5+4.d0*kxel*y2*x3*v6+
718  & 2.d0*kyel*x2*y2*u1-kyel*x2*y3*u1+2.d0*kyel*x2*y2*u3+
719  & 4.d0*kyel*x2*y3*u6+12.d0*kyel*x2**2*v4-
720  & 6.d0*kyel*x3*x2*v2+4.d0*kyel*x3*y2*u4+
721  & 8.d0*kxel*y3*x2*v5-16.d0*kxel*y3*y2*u5-
722  & 4.d0*kxel*y3*x3*v5-8.d0*kyel*x3*x2*v4-
723  & 16.d0*kyel*x3*x2*v5-4.d0*kyel*x3*y3*u5+
724  & 8.d0*kyel*x3*y2*u5+6.d0*kxel*y3*x2*v2+
725  & 4.d0*kyel*x2*y3*u4+4.d0*kxel*y2*x3*v4-
726  & 4.d0*kxel*y2*x2*v6-kyel*x3*y2*u3-
727  & 4.d0*kyel*x2*y2*u6-kxel*y3*x2*v3
728  a45(ielem)= (ans1+ans2)*2.d0*xmul/(45.d0*(-x2*y3+y2*x3))
729 !
730  ans1 = -kyel*x3*y2*u2-kxel*y3**2*u1-4.d0*kxel*y3*y2*u1+
731  & 5.d0*kxel*y3*x2*v1-kxel*y3*x3*v2+4.d0*kxel*y3*x2*v4-
732  & 8.d0*kxel*y3*y2*u4+5.d0*kyel*x3*y2*u1+kyel*x3*y3*u1-
733  & 4.d0*kyel*x3*x2*v1+kxel*y3*y2*u2-4.d0*kxel*y3**2*u6+
734  & 4.d0*kyel*x3**2*v5-kyel*x3**2*v1+kxel*y3*x3*v1+
735  & 8.d0*kyel*x2*y3*u5-4.d0*kyel*x2*y2*u5+
736  & 4.d0*kyel*x2*y2*u4+kxel*y3**2*u2+kyel*x3**2*v2-
737  & 4.d0*kyel*x3**2*v6-4.d0*kxel*y2**2*u4+kyel*x2**2*v3+
738  & 4.d0*kxel*y2**2*u5+4.d0*kyel*x2**2*v5+kxel*y2*x2*v1-
739  & kxel*y2*x3*v1-kxel*y2*x2*v3-kxel*y2**2*u1+
740  & kxel*y2**2*u3-kyel*x2**2*v1-4.d0*kxel*y2*x2*v5
741  ans2 = 8.d0*kxel*y2*x3*v5+4.d0*kxel*y2*x2*v4-kyel*x3*y3*u2+
742  & kyel*x3*x2*v3+4.d0*kyel*x3*y2*u6-8.d0*kyel*x3*x2*v6+
743  & 4.d0*kyel*x3*y3*u6+kxel*y3*y2*u3+4.d0*kxel*y3*x3*v6-
744  & 8.d0*kxel*y3*y2*u6+4.d0*kxel*y3**2*u5+
745  & 4.d0*kxel*y2*x3*v6+kyel*x2*y2*u1-kyel*x2*y3*u1-
746  & kyel*x2*y2*u3+4.d0*kyel*x2*y3*u6-4.d0*kyel*x2**2*v4+
747  & kyel*x3*x2*v2+4.d0*kyel*x3*y2*u4+4.d0*kxel*y3*x2*v5-
748  & 12.d0*kxel*y3*y2*u5-4.d0*kxel*y3*x3*v5-
749  & 8.d0*kyel*x3*x2*v4+4.d0*kxel*y3*x2*v6-
750  & 12.d0*kyel*x3*x2*v5-4.d0*kyel*x3*y3*u5+
751  & 4.d0*kyel*x3*y2*u5-kxel*y3*x2*v2+4.d0*kyel*x2*y3*u4+
752  & 4.d0*kxel*y2*x3*v4-kyel*x3*y2*u3-kxel*y3*x2*v3
753  a46(ielem)= (ans1+ans2)*2.d0*xmul/(45.d0*(x2*y3-y2*x3))
754 !
755  ans1 = -kyel*x3*y2*u2+6.d0*kyel*x2*y2*u2+6.d0*kxel*y2*x2*v2+
756  & 4.d0*kyel*x2**2*v6+3.d0*kxel*y3**2*u1-
757  & 6.d0*kxel*y3*y2*u1+3.d0*kxel*y3*x2*v1-
758  & 6.d0*kyel*x2**2*v2+kxel*y3*x3*v2+4.d0*kxel*y3*x2*v4-
759  & 6.d0*kxel*y2**2*u2+4.d0*kxel*y3**2*u4-
760  & 12.d0*kxel*y3*y2*u4+4.d0*kxel*y2**2*u6+
761  & 3.d0*kyel*x3*y2*u1-3.d0*kyel*x3*y3*u1-
762  & 6.d0*kyel*x3*x2*v1-4.d0*kxel*y3*x3*v4+
763  & 7.d0*kxel*y3*y2*u2+8.d0*kxel*y3**2*u6-
764  & 8.d0*kyel*x3**2*v5+3.d0*kyel*x3**2*v1-
765  & 3.d0*kxel*y3*x3*v1-8.d0*kyel*x2*y3*u5+
766  & 8.d0*kyel*x2*y2*u5-8.d0*kyel*x2*y2*u4-
767  & kxel*y3**2*u2-kyel*x3**2*v2+8.d0*kyel*x3**2*v6-
768  & 6.d0*kyel*x3**2*v3-6.d0*kxel*y3**2*u3+
769  & 8.d0*kxel*y2**2*u4-kyel*x2**2*v3-8.d0*kxel*y2**2*u5-
770  & 8.d0*kyel*x2**2*v5-3.d0*kxel*y2*x2*v1+
771  & 3.d0*kxel*y2*x3*v1+kxel*y2*x2*v3+3.d0*kxel*y2**2*u1-
772  & kxel*y2**2*u3+3.d0*kyel*x2**2*v1+8.d0*kxel*y2*x2*v5
773  ans2 = -8.d0*kxel*y2*x3*v5-8.d0*kxel*y2*x2*v4+
774  & kyel*x3*y3*u2+6.d0*kyel*x3*y3*u3+7.d0*kyel*x3*x2*v3+
775  & 8.d0*kyel*x3*y2*u6-12.d0*kyel*x3*x2*v6-
776  & 8.d0*kyel*x3*y3*u6+7.d0*kxel*y3*y2*u3-
777  & 8.d0*kxel*y3*x3*v6+6.d0*kxel*y3*x3*v3-
778  & 12.d0*kxel*y3*y2*u6-8.d0*kxel*y3**2*u5+
779  & 4.d0*kyel*x3**2*v4+4.d0*kxel*y2*x3*v6-
780  & 3.d0*kyel*x2*y2*u1+3.d0*kyel*x2*y3*u1+
781  & kyel*x2*y2*u3-kyel*x2*y3*u3+4.d0*kyel*x2*y3*u6-
782  & kxel*y2*x3*v3+8.d0*kyel*x2**2*v4+7.d0*kyel*x3*x2*v2-
783  & 4.d0*kyel*x3*y3*u4+4.d0*kyel*x3*y2*u4-
784  & 8.d0*kxel*y3*x2*v5+16.d0*kxel*y3*y2*u5+
785  & 8.d0*kxel*y3*x3*v5-12.d0*kyel*x3*x2*v4+
786  & 8.d0*kxel*y3*x2*v6+16.d0*kyel*x3*x2*v5+
787  & 8.d0*kyel*x3*y3*u5-8.d0*kyel*x3*y2*u5-
788  & 6.d0*kyel*x2*y3*u2-6.d0*kxel*y2*x3*v2-
789  & kxel*y3*x2*v2+8.d0*kyel*x2*y3*u4+8.d0*kxel*y2*x3*v4-
790  & 4.d0*kxel*y2*x2*v6-6.d0*kyel*x3*y2*u3-
791  & 4.d0*kyel*x2*y2*u6-6.d0*kxel*y3*x2*v3
792  a51(ielem) = (ans1+ans2)*xmul/(90.d0*(-x2*y3+y2*x3))
793 !
794  ans1 = -kxel*y3**2*u1+5.d0*kxel*y3*y2*u1-3.d0*kxel*y3*x3*v2+
795  & 4.d0*kxel*y3**2*u4-24.d0*kxel*y3*y2*u4+kyel*x3*y3*u1+
796  & 5.d0*kyel*x3*x2*v1-4.d0*kxel*y3*x3*v4-
797  & 18.d0*kxel*y3*y2*u2-8.d0*kxel*y3**2*u6+
798  & 8.d0*kyel*x3**2*v5-kyel*x3**2*v1+kxel*y3*x3*v1+
799  & 24.d0*kyel*x2*y3*u5+3.d0*kxel*y3**2*u2+
800  & 3.d0*kyel*x3**2*v2-8.d0*kyel*x3**2*v6-
801  & 6.d0*kyel*x3**2*v3-6.d0*kxel*y3**2*u3-
802  & 5.d0*kxel*y2*x3*v1+24.d0*kxel*y2*x3*v5-
803  & 3.d0*kyel*x3*y3*u2+6.d0*kyel*x3*y3*u3
804  ans2 = 5.d0*kyel*x3*x2*v3 -4.d0*kyel*x3*x2*v6+
805  & 8.d0*kyel*x3*y3*u6+5.d0*kxel*y3*y2*u3+
806  & 8.d0*kxel*y3*x3*v6+6.d0*kxel*y3*x3*v3-
807  & 4.d0*kxel*y3*y2*u6+8.d0*kxel*y3**2*u5+
808  & 4.d0*kyel*x3**2*v4+4.d0*kxel*y2*x3*v6-
809  & 5.d0*kyel*x2*y3*u1-5.d0*kyel*x2*y3*u3+
810  & 4.d0*kyel*x2*y3*u6-5.d0*kxel*y2*x3*v3-
811  & 18.d0*kyel*x3*x2*v2-4.d0*kyel*x3*y3*u4-
812  & 24.d0*kxel*y3*y2*u5-8.d0*kxel*y3*x3*v5-
813  & 24.d0*kyel*x3*x2*v4-24.d0*kyel*x3*x2*v5-
814  & 8.d0*kyel*x3*y3*u5+18.d0*kyel*x2*y3*u2+
815  & 18.d0*kxel*y2*x3*v2+24.d0*kyel*x2*y3*u4+
816  & 24.d0*kxel*y2*x3*v4
817  a52(ielem)= (ans1+ans2)*xmul/(90.d0*(x2*y3-y2*x3))
818 !
819  ans1 = 5.d0*kyel*x3*y2*u2-6.d0*kyel*x2*y2*u2-
820  & 6.d0*kxel*y2*x2*v2-4.d0*kyel*x2**2*v6-
821  & 5.d0*kxel*y3*y2*u1+5.d0*kxel*y3*x2*v1+
822  & 6.d0*kyel*x2**2*v2-4.d0*kxel*y3*x2*v4+
823  & 6.d0*kxel*y2**2*u2+4.d0*kxel*y3*y2*u4-
824  & 4.d0*kxel*y2**2*u6+5.d0*kyel*x3*y2*u1-
825  & 5.d0*kyel*x3*x2*v1-5.d0*kxel*y3*y2*u2+
826  & 8.d0*kyel*x2*y2*u5-8.d0*kyel*x2*y2*u4+
827  & 8.d0*kxel*y2**2*u4-3.d0*kyel*x2**2*v3-
828  & 8.d0*kxel*y2**2*u5-8.d0*kyel*x2**2*v5-
829  & kxel*y2*x2*v1+3.d0*kxel*y2*x2*v3+kxel*y2**2*u1-
830  & 3.d0*kxel*y2**2*u3
831  ans2 = kyel*x2**2*v1+8*kxel*y2*x2*v5-8.d0*kxel*y2*x2*v4+
832  & 18.d0*kyel*x3*x2*v3-24.d0*kyel*x3*y2*u6+
833  & 24.d0*kyel*x3*x2*v6+18.d0*kxel*y3*y2*u3+
834  & 24.d0*kxel*y3*y2*u6-kyel*x2*y2*u1+3.d0*kyel*x2*y2*u3+
835  & 8.d0*kyel*x2**2*v4-5.d0*kyel*x3*x2*v2-
836  & 4.d0*kyel*x3*y2*u4-24.d0*kxel*y3*x2*v5+
837  & 24.d0*kxel*y3*y2*u5+4.d0*kyel*x3*x2*v4-
838  & 24.d0*kxel*y3*x2*v6+24.d0*kyel*x3*x2*v5-
839  & 24.d0*kyel*x3*y2*u5+5.d0*kxel*y3*x2*v2+
840  & 4.d0*kxel*y2*x2*v6-18.d0*kyel*x3*y2*u3+
841  & 4.d0*kyel*x2*y2*u6-18.d0*kxel*y3*x2*v3
842  a53(ielem)= (ans1+ans2)*xmul/(90.d0*(-x2*y3+y2*x3))
843 !
844  ans1 = -6.d0*kyel*x2*y2*u2-6.d0*kxel*y2*x2*v2+
845  & 4.d0*kyel*x2**2*v6-kxel*y3**2*u1+3.d0*kxel*y3*y2*u1-
846  & kxel*y3*x2*v1+6.d0*kyel*x2**2*v2-kxel*y3*x3*v2+
847  & 4.d0*kxel*y3*x2*v4+6.d0*kxel*y2**2*u2-
848  & 8.d0*kxel*y3*y2*u4+4.d0*kxel*y2**2*u6-
849  & kyel*x3*y2*u1+kyel*x3*y3*u1+3.d0*kyel*x3*x2*v1-
850  & 6.d0*kxel*y3*y2*u2-4.d0*kxel*y3**2*u6+
851  & 4.d0*kyel*x3**2*v5-kyel*x3**2*v1+kxel*y3*x3*v1+
852  & 8.d0*kyel*x2*y3*u5-12.d0*kyel*x2*y2*u5-
853  & 12.d0*kyel*x2*y2*u4+kxel*y3**2*u2+kyel*x3**2*v2-
854  & 4.d0*kyel*x3**2*v6+12.d0*kxel*y2**2*u4-
855  & 2.d0*kyel*x2**2*v3+12.d0*kxel*y2**2*u5+
856  & 12.d0*kyel*x2**2*v5+2.d0*kxel*y2*x2*v1-
857  & 2.d0*kxel*y2*x3*v1+2.d0*kxel*y2*x2*v3-
858  & 2.d0*kxel*y2**2*u1-2.d0*kxel*y2**2*u3
859  ans2 = -2.d0*kyel*x2**2*v1-12.d0*kxel*y2*x2*v5+
860  & 8.d0*kxel*y2*x3*v5-12.d0*kxel*y2*x2*v4-
861  & kyel*x3*y3*u2+kyel*x3*x2*v3+4.d0*kyel*x3*y2*u6-
862  & 4.d0*kyel*x3*x2*v6+4.d0*kyel*x3*y3*u6+kxel*y3*y2*u3+
863  & 4.d0*kxel*y3*x3*v6-4.d0*kxel*y3*y2*u6+
864  & 4.d0*kxel*y3**2*u5+2.d0*kyel*x2*y2*u1-
865  & 2.d0*kyel*x2*y3*u1+2.d0*kyel*x2*y2*u3-
866  & kyel*x2*y3*u3-kxel*y2*x3*v3+12.d0*kyel*x2**2*v4-
867  & 6.d0*kyel*x3*x2*v2+4.d0*kyel*x3*y2*u4+
868  & 8.d0*kxel*y3*x2*v5-16.d0*kxel*y3*y2*u5-
869  & 4.d0*kxel*y3*x3*v5-8.d0*kyel*x3*x2*v4+
870  & 4.d0*kxel*y3*x2*v6-16.d0*kyel*x3*x2*v5-
871  & 4.d0*kyel*x3*y3*u5+8.d0*kyel*x3*y2*u5+
872  & 6.d0*kyel*x2*y3*u2+6.d0*kxel*y2*x3*v2+
873  & 4.d0*kyel*x2*y3*u4+4.d0*kxel*y2*x3*v4-
874  & 4.d0*kxel*y2*x2*v6-4.d0*kyel*x2*y2*u6
875  a54(ielem)= (ans1+ans2)*2.d0*xmul/(45.d0*(-x2*y3+y2*x3))
876 !
877  ans1 = -kyel*x3*y2*u2-2.d0*kxel*y3**2*u1+3.d0*kxel*y3*y2*u1-
878  & 2.d0*kxel*y3*x2*v1+2.d0*kxel*y3*x3*v2+
879  & 4.d0*kxel*y3**2*u4-4.d0*kxel*y3*y2*u4-
880  & 2.d0*kyel*x3*y2*u1+2.d0*kyel*x3*y3*u1+
881  & 3.d0*kyel*x3*x2*v1-4.d0*kxel*y3*x3*v4+
882  & kxel*y3*y2*u2+12.d0*kxel*y3**2*u6+
883  & 12.d0*kyel*x3**2*v5-2.d0*kyel*x3**2*v1+
884  & 2.d0*kxel*y3*x3*v1+8.d0*kyel*x2*y3*u5-
885  & 4.d0*kyel*x2*y2*u5+4.d0*kyel*x2*y2*u4-
886  & 2.d0*kxel*y3**2*u2-2.d0*kyel*x3**2*v2+
887  & 12.d0*kyel*x3**2*v6+6.d0*kyel*x3**2*v3+
888  & 6.d0*kxel*y3**2*u3-4.d0*kxel*y2**2*u4+
889  & kyel*x2**2*v3+4.d0*kxel*y2**2*u5+4.d0*kyel*x2**2*v5+
890  & kxel*y2*x2*v1-kxel*y2*x3*v1-kxel*y2*x2*v3-
891  & kxel*y2**2*u1+kxel*y2**2*u3-kyel*x2**2*v1-
892  & 4.d0*kxel*y2*x2*v5
893  ans2 = 8.d0*kxel*y2*x3*v5+4.d0*kxel*y2*x2*v4+
894  & 2.d0*kyel*x3*y3*u2-6.d0*kyel*x3*y3*u3-
895  & 6.d0*kyel*x3*x2*v3+4.d0*kyel*x3*y2*u6-
896  & 8.d0*kyel*x3*x2*v6-12.d0*kyel*x3*y3*u6-
897  & 6.d0*kxel*y3*y2*u3-12.d0*kxel*y3*x3*v6-
898  & 6.d0*kxel*y3*x3*v3-8.d0*kxel*y3*y2*u6+
899  & 12.d0*kxel*y3**2*u5+4.d0*kyel*x3**2*v4+
900  & 4.d0*kxel*y2*x3*v6+kyel*x2*y2*u1-kyel*x2*y3*u1-
901  & kyel*x2*y2*u3+4.d0*kyel*x2*y3*u6-4.d0*kyel*x2**2*v4+
902  & kyel*x3*x2*v2-4.d0*kyel*x3*y3*u4+8.d0*kxel*y3*x2*v5-
903  & 16.d0*kxel*y3*y2*u5-12.d0*kxel*y3*x3*v5-
904  & 4.d0*kyel*x3*x2*v4+4.d0*kxel*y3*x2*v6-
905  & 16.d0*kyel*x3*x2*v5-12.d0*kyel*x3*y3*u5+
906  & 8.d0*kyel*x3*y2*u5-kxel*y3*x2*v2+4.d0*kyel*x2*y3*u4+
907  & 4.d0*kxel*y2*x3*v4+6.d0*kyel*x3*y2*u3+6.d0*kxel*y3*x2*v3
908  a56(ielem)= (ans1+ans2)*2.d0*xmul/(45.d0*(-x2*y3+y2*x3))
909 !
910  ans1 = -kyel*x3*y2*u2-4.d0*kyel*x2*y2*u2-4.d0*kxel*y2*x2*v2-
911  & 16.d0*kyel*x2**2*v6+3.d0*kxel*y3**2*u1+
912  & 12.d0*kxel*y3*y2*u1+3.d0*kxel*y3*x2*v1+
913  & 4.d0*kyel*x2**2*v2+kxel*y3*x3*v2+4.d0*kxel*y3*x2*v4+
914  & 4.d0*kxel*y2**2*u2+4.d0*kxel*y3**2*u4+
915  & 16.d0*kxel*y3*y2*u4-16.d0*kxel*y2**2*u6+
916  & 3.d0*kyel*x3*y2*u1-3.d0*kyel*x3*y3*u1+
917  & 12.d0*kyel*x3*x2*v1-4.d0*kxel*y3*x3*v4-
918  & 3.d0*kxel*y3*y2*u2+8.d0*kxel*y3**2*u6-
919  & 8.d0*kyel*x3**2*v5+3.d0*kyel*x3**2*v1-
920  & 3.d0*kxel*y3*x3*v1-12.d0*kyel*x2*y3*u5+
921  & 12.d0*kyel*x2*y2*u5+20.d0*kyel*x2*y2*u4-kxel*y3**2*u2-
922  & kyel*x3**2*v2+8.d0*kyel*x3**2*v6-6.d0*kyel*x3**2*v3-
923  & 6.d0*kxel*y3**2*u3-20.d0*kxel*y2**2*u4-kyel*x2**2*v3-
924  & 12.d0*kxel*y2**2*u5-12.d0*kyel*x2**2*v5+
925  & 15.d0*kxel*y2*x2*v1-15.d0*kxel*y2*x3*v1+
926  & kxel*y2*x2*v3-15.d0*kxel*y2**2*u1-kxel*y2**2*u3-
927  & 15.d0*kyel*x2**2*v1+12.d0*kxel*y2*x2*v5
928  ans2 = -12.d0*kxel*y2*x3*v5+20.d0*kxel*y2*x2*v4+
929  & kyel*x3*y3*u2+6.d0*kyel*x3*y3*u3+7.d0*kyel*x3*x2*v3+
930  & 8.d0*kyel*x3*y2*u6+8.d0*kyel*x3*x2*v6-
931  & 8.d0*kyel*x3*y3*u6+7.d0*kxel*y3*y2*u3-
932  & 8.d0*kxel*y3*x3*v6+6.d0*kxel*y3*x3*v3+
933  & 8.d0*kxel*y3*y2*u6-8.d0*kxel*y3**2*u5+
934  & 4.d0*kyel*x3**2*v4-16.d0*kxel*y2*x3*v6+
935  & 15.d0*kyel*x2*y2*u1-15.d0*kyel*x2*y3*u1+
936  & kyel*x2*y2*u3-kyel*x2*y3*u3-16.d0*kyel*x2*y3*u6-
937  & kxel*y2*x3*v3-20.d0*kyel*x2**2*v4-3.d0*kyel*x3*x2*v2-
938  & 4.d0*kyel*x3*y3*u4+4.d0*kyel*x3*y2*u4-
939  & 8.d0*kxel*y3*x2*v5+20.d0*kxel*y3*y2*u5+
940  & 8.d0*kxel*y3*x3*v5+16.d0*kyel*x3*x2*v4+
941  & 8.d0*kxel*y3*x2*v6+20.d0*kyel*x3*x2*v5+
942  & 8.d0*kyel*x3*y3*u5-8.d0*kyel*x3*y2*u5+
943  & 4.d0*kyel*x2*y3*u2+4.d0*kxel*y2*x3*v2-
944  & kxel*y3*x2*v2-20.d0*kyel*x2*y3*u4-20.d0*kxel*y2*x3*v4+
945  & 16.d0*kxel*y2*x2*v6-6.d0*kyel*x3*y2*u3+16.d0*kyel*x2*
946  & y2*u6-6.d0*kxel*y3*x2*v3
947  a61(ielem)= (ans1+ans2)*xmul/(90.d0*(x2*y3-y2*x3))
948 !
949  ans1 = -kxel*y3**2*u1-5.d0*kxel*y3*y2*u1-
950  & 3.d0*kxel*y3*x3*v2+4.d0*kxel*y3**2*u4+
951  & 4.d0*kxel*y3*y2*u4+kyel*x3*y3*u1-5.d0*kyel*x3*x2*v1-
952  & 4.d0*kxel*y3*x3*v4-8.d0*kxel*y3**2*u6+
953  & 8.d0*kyel*x3**2*v5-kyel*x3**2*v1+kxel*y3*x3*v1+
954  & 4.d0*kyel*x2*y3*u5+3.d0*kxel*y3**2*u2+
955  & 3.d0*kyel*x3**2*v2-8.d0*kyel*x3**2*v6-
956  & 6.d0*kyel*x3**2*v3-6.d0*kxel*y3**2*u3+
957  & 5.d0*kxel*y2*x3*v1+4.d0*kxel*y2*x3*v5-
958  & 3.d0*kyel*x3*y3*u2+6.d0*kyel*x3*y3*u3+
959  & 5.d0*kyel*x3*x2*v3+8.d0*kyel*x3*y3*u6+
960  & 5.d0*kxel*y3*y2*u3+8.d0*kxel*y3*x3*v6+
961  & 6.d0*kxel*y3*x3*v3+8.d0*kxel*y3**2*u5+
962  & 4.d0*kyel*x3**2*v4+5.d0*kyel*x2*y3*u1-
963  & 5.d0*kyel*x2*y3*u3-5.d0*kxel*y2*x3*v3-
964  & 4.d0*kyel*x3*y3*u4-4.d0*kxel*y3*y2*u5-
965  & 8.d0*kxel*y3*x3*v5+4.d0*kyel*x3*x2*v4-
966  & 4.d0*kyel*x3*x2*v5-8.d0*kyel*x3*y3*u5-
967  & 4.d0*kyel*x2*y3*u4-4.d0*kxel*y2*x3*v4
968  a62(ielem)= ans1*xmul/(90.d0*(-x2*y3+y2*x3))
969 !
970  ans1 = 5.d0*kyel*x3*y2*u2-4.d0*kyel*x2*y2*u2-
971  & 4.d0*kxel*y2*x2*v2-16.d0*kyel*x2**2*v6-
972  & 5.d0*kxel*y3*y2*u1+5.d0*kxel*y3*x2*v1+
973  & 4.d0*kyel*x2**2*v2-4.d0*kxel*y3*x2*v4+
974  & 4.d0*kxel*y2**2*u2+4.d0*kxel*y3*y2*u4-
975  & 16.d0*kxel*y2**2*u6+5.d0*kyel*x3*y2*u1-
976  & 5.d0*kyel*x3*x2*v1-5.d0*kxel*y3*y2*u2+
977  & 20.d0*kyel*x2*y2*u5+12.d0*kyel*x2*y2*u4-
978  & 12.d0*kxel*y2**2*u4-15.d0*kyel*x2**2*v3-
979  & 20.d0*kxel*y2**2*u5-20.d0*kyel*x2**2*v5+
980  & kxel*y2*x2*v1+15.d0*kxel*y2*x2*v3-kxel*y2**2*u1-
981  & 15.d0*kxel*y2**2*u3-kyel*x2**2*v1
982  ans2 = 20.d0*kxel*y2*x2*v5+12.d0*kxel*y2*x2*v4+
983  & 18.d0*kyel*x3*x2*v3-24.d0*kyel*x3*y2*u6+
984  & 24.d0*kyel*x3*x2*v6+18.d0*kxel*y3*y2*u3+
985  & 24.d0*kxel*y3*y2*u6+kyel*x2*y2*u1+
986  & 15.d0*kyel*x2*y2*u3-12.d0*kyel*x2**2*v4-
987  & 5.d0*kyel*x3*x2*v2-4.d0*kyel*x3*y2*u4-
988  & 24.d0*kxel*y3*x2*v5+24.d0*kxel*y3*y2*u5+
989  & 4.d0*kyel*x3*x2*v4-24.d0*kxel*y3*x2*v6+
990  & 24.d0*kyel*x3*x2*v5-24.d0*kyel*x3*y2*u5+
991  & 5.d0*kxel*y3*x2*v2+16.d0*kxel*y2*x2*v6-
992  & 18.d0*kyel*x3*y2*u3+16.d0*kyel*x2*y2*u6-
993  & 18.d0*kxel*y3*x2*v3
994  a63(ielem)= (ans1+ans2)*xmul/(90.d0*(x2*y3-y2*x3))
995 !
996  ans1 = -kxel*y3**2*u1-4.d0*kxel*y3*y2*u1-kxel*y3*x2*v1-
997  & kxel*y3*x3*v2+4.d0*kxel*y3*x2*v4-8.d0*kxel*y3*y2*u4-
998  & kyel*x3*y2*u1+kyel*x3*y3*u1-4.d0*kyel*x3*x2*v1+
999  & kxel*y3*y2*u2-4.d0*kxel*y3**2*u6+4.d0*kyel*x3**2*v5-
1000  & kyel*x3**2*v1+kxel*y3*x3*v1+4.d0*kyel*x2*y3*u5-
1001  & 4.d0*kyel*x2*y2*u5+4.d0*kyel*x2*y2*u4+kxel*y3**2*u2+
1002  & kyel*x3**2*v2-4.d0*kyel*x3**2*v6-4.d0*kxel*y2**2*u4+
1003  & kyel*x2**2*v3+4.d0*kxel*y2**2*u5+4.d0*kyel*x2**2*v5+
1004  & kxel*y2*x2*v1+5.d0*kxel*y2*x3*v1-kxel*y2*x2*v3-
1005  & kxel*y2**2*u1+kxel*y2**2*u3-kyel*x2**2*v1-
1006  & 4.d0*kxel*y2*x2*v5+4.d0*kxel*y2*x3*v5+
1007  & 4.d0*kxel*y2*x2*v4-kyel*x3*y3*u2+kyel*x3*x2*v3+
1008  & 4.d0*kyel*x3*y2*u6-8.d0*kyel*x3*x2*v6
1009  ans2 = 4.d0*kyel*x3*y3*u6+kxel*y3*y2*u3+4.d0*kxel*y3*x3*v6-
1010  & 8.d0*kxel*y3*y2*u6+4.d0*kxel*y3**2*u5+
1011  & 4.d0*kxel*y2*x3*v6+kyel*x2*y2*u1+5.d0*kyel*x2*y3*u1-
1012  & kyel*x2*y2*u3-kyel*x2*y3*u3+4.d0*kyel*x2*y3*u6-
1013  & kxel*y2*x3*v3-4.d0*kyel*x2**2*v4+kyel*x3*x2*v2+
1014  & 4.d0*kyel*x3*y2*u4+8.d0*kxel*y3*x2*v5-
1015  & 12.d0*kxel*y3*y2*u5-4.d0*kxel*y3*x3*v5-
1016  & 8.d0*kyel*x3*x2*v4+4.d0*kxel*y3*x2*v6-
1017  & 12.d0*kyel*x3*x2*v5-4.d0*kyel*x3*y3*u5+
1018  & 8.d0*kyel*x3*y2*u5-kyel*x2*y3*u2-kxel*y2*x3*v2+
1019  & 4.d0*kyel*x2*y3*u4+4.d0*kxel*y2*x3*v4
1020  a64(ielem)= (ans1+ans2)*xmul*2.d0/(45.d0*(x2*y3-y2*x3))
1021 !
1022  ans1 = -2.d0*kxel*y3**2*u1+3.d0*kxel*y3*y2*u1-kxel*y3*x2*v1+
1023  & 2.d0*kxel*y3*x3*v2+4.d0*kxel*y3*x2*v4+
1024  & 4.d0*kxel*y3**2*u4-4.d0*kxel*y3*y2*u4-
1025  & kyel*x3*y2*u1+2.d0*kyel*x3*y3*u1+3.d0*kyel*x3*x2*v1-
1026  & 4.d0*kxel*y3*x3*v4+kxel*y3*y2*u2+12.d0*kxel*y3**2*u6+
1027  & 12.d0*kyel*x3**2*v5-2.d0*kyel*x3**2*v1+
1028  & 2.d0*kxel*y3*x3*v1+8.d0*kyel*x2*y3*u5-
1029  & 4.d0*kyel*x2*y2*u5+4.d0*kyel*x2*y2*u4-
1030  & 2.d0*kxel*y3**2*u2-2.d0*kyel*x3**2*v2+
1031  & 12.d0*kyel*x3**2*v6+6.d0*kyel*x3**2*v3+
1032  & 6.d0*kxel*y3**2*u3-4.d0*kxel*y2**2*u4+
1033  & kyel*x2**2*v3+4.d0*kxel*y2**2*u5+4.d0*kyel*x2**2*v5+
1034  & kxel*y2*x2*v1-2.d0*kxel*y2*x3*v1-kxel*y2*x2*v3-
1035  & kxel*y2**2*u1+kxel*y2**2*u3-kyel*x2**2*v1-
1036  & 4.d0*kxel*y2*x2*v5
1037  ans2 = 8.d0*kxel*y2*x3*v5+4.d0*kxel*y2*x2*v4+
1038  & 2.d0*kyel*x3*y3*u2-6.d0*kyel*x3*y3*u3-
1039  & 6.d0*kyel*x3*x2*v3+4.d0*kyel*x3*y2*u6-
1040  & 8.d0*kyel*x3*x2*v6-12.d0*kyel*x3*y3*u6-
1041  & 6.d0*kxel*y3*y2*u3-12.d0*kxel*y3*x3*v6-
1042  & 6.d0*kxel*y3*x3*v3-8.d0*kxel*y3*y2*u6+
1043  & 12.d0*kxel*y3**2*u5+4.d0*kyel*x3**2*v4+
1044  & 4.d0*kxel*y2*x3*v6+kyel*x2*y2*u1-2.d0*kyel*x2*y3*u1-
1045  & kyel*x2*y2*u3+6.d0*kyel*x2*y3*u3+4.d0*kyel*x2*y3*u6+
1046  & 6.d0*kxel*y2*x3*v3-4.d0*kyel*x2**2*v4+kyel*x3*x2*v2-
1047  & 4.d0*kyel*x3*y3*u4+4.d0*kyel*x3*y2*u4+
1048  & 8.d0*kxel*y3*x2*v5-16.d0*kxel*y3*y2*u5-
1049  & 12.d0*kxel*y3*x3*v5-4.d0*kyel*x3*x2*v4+
1050  & 4.d0*kxel*y3*x2*v6-16.d0*kyel*x3*x2*v5-
1051  & 12.d0*kyel*x3*y3*u5+8.d0*kyel*x3*y2*u5-kyel*x2*y3*u2-
1052  & kxel*y2*x3*v2
1053  a65(ielem)= (ans1+ans2)*xmul*2.d0/(45.d0*(-x2*y3+y2*x3))
1054 !
1055 ! USES HERE THE 'MAGIC SQUARE' PROPERTIES
1056 ! (SUM OF EACH LINE = SUM OF EACH COLUMN = 0)
1057 !
1058  a16(ielem) = - a11(ielem) - a12(ielem) - a13(ielem)
1059  & - a14(ielem) - a15(ielem)
1060  a25(ielem) = - a21(ielem) - a22(ielem) - a23(ielem)
1061  & - a24(ielem) - a26(ielem)
1062  a36(ielem) = - a31(ielem) - a32(ielem) - a33(ielem)
1063  & - a34(ielem) - a35(ielem)
1064  a43(ielem) = - a41(ielem) - a42(ielem) - a44(ielem)
1065  & - a45(ielem) - a46(ielem)
1066  a55(ielem) = - a51(ielem) - a52(ielem) - a53(ielem)
1067  & - a54(ielem) - a56(ielem)
1068  a66(ielem) = - a61(ielem) - a62(ielem) - a63(ielem)
1069  & - a64(ielem) - a65(ielem)
1070 !
1071  ENDDO ! IELEM
1072 !
1073 ! OTHER TYPES OF FUNCTIONS F AND G
1074 !
1075 !-----------------------------------------------------------------------
1076 !
1077  ELSE
1078 !
1079  WRITE(lu,101) ielmf,sf%NAME
1080  WRITE(lu,111) ielmg,sg%NAME
1081  WRITE(lu,201) ielmu,su%NAME
1082  WRITE(lu,301)
1083 101 FORMAT(1x,'MT03CC (BIEF) :',/,
1084  & 1x,'DISCRETIZATION OF F:',1i6,
1085  & 1x,'REAL NAME: ',a6)
1086 111 FORMAT(1x,'DISCRETIZATION OF G:',1i6,
1087  & 1x,'REAL NAME: ',a6)
1088 201 FORMAT(1x,'DISCRETIZATION OF U:',1i6,
1089  & 1x,'REAL NAME: ',a6)
1090 301 FORMAT(1x,'CASE NOT IMPLEMENTED')
1091  CALL plante(1)
1092  stop
1093 !
1094  ENDIF
1095 !
1096 !-----------------------------------------------------------------------
1097 !
1098  RETURN
1099  END
mt03cc
subroutine mt03cc(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, SF, SG, SU, SV, F, G, U, V, XEL, YEL, IKLE1, IKLE2, IKLE3, IKLE4, IKLE5, IKLE6, NELEM, NELMAX)
Definition: mt03cc.f:15
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:36:52 for The TELEMAC-MASCARET system by doxygen 1.8.13