The TELEMAC-MASCARET system  trunk
mt12bb.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt12bb
3 ! *****************
4 !
5  &( a11 , a12 , a13 , a14 ,
6  & a21 , a22 , a23 , a24 ,
7  & a31 , a32 , a33 , a34 ,
8  & a41 , a42 , a43 , a44 ,
9  & xmul,sf,su,sv,f,u,v,
10  & xel,yel,
11  & ikle1,ikle2,ikle3,ikle4,
12  & nelem,nelmax,icoord)
13 !
14 !***********************************************************************
15 ! BIEF V6P0 21/08/2010
16 !***********************************************************************
17 !
18 !brief COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
19 !code
20 !+ EXAMPLE WITH ICOORD=1
21 !+
22 !+ / D -> -->
23 !+ A(I,J)= XMUL / PSI2(J) * --(F) * U.GRAD(PSI1(I)) D(OMEGA)
24 !+ /OMEGA DX
25 !+
26 !+
27 !+ PSI1: BASES OF TYPE P1 TRIANGLE
28 !+ PSI2: BASES OF TYPE IELM2
29 !+
30 !+ IT WOULD BE A DERIVATIVE WRT Y WITH ICOORD=2
31 !
32 !warning THE JACOBIAN MUST BE POSITIVE
33 !
34 !history J-M HERVOUET (LNH) ; C MOULIN (LNH)
35 !+ 09/12/94
36 !+ V5P1
37 !+
38 !
39 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
40 !+ 13/07/2010
41 !+ V6P0
42 !+ Translation of French comments within the FORTRAN sources into
43 !+ English comments
44 !
45 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
46 !+ 21/08/2010
47 !+ V6P0
48 !+ Creation of DOXYGEN tags for automated documentation and
49 !+ cross-referencing of the FORTRAN sources
50 !
51 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
52 !| A11 |<--| ELEMENTS OF MATRIX
53 !| A12 |<--| ELEMENTS OF MATRIX
54 !| A13 |<--| ELEMENTS OF MATRIX
55 !| A14 |<--| ELEMENTS OF MATRIX
56 !| A21 |<--| ELEMENTS OF MATRIX
57 !| A22 |<--| ELEMENTS OF MATRIX
58 !| A23 |<--| ELEMENTS OF MATRIX
59 !| A24 |<--| ELEMENTS OF MATRIX
60 !| A31 |<--| ELEMENTS OF MATRIX
61 !| A32 |<--| ELEMENTS OF MATRIX
62 !| A33 |<--| ELEMENTS OF MATRIX
63 !| A34 |<--| ELEMENTS OF MATRIX
64 !| A41 |<--| ELEMENTS OF MATRIX
65 !| A42 |<--| ELEMENTS OF MATRIX
66 !| A43 |<--| ELEMENTS OF MATRIX
67 !| A44 |<--| ELEMENTS OF MATRIX
68 !| F |-->| FUNCTION USED IN THE FORMULA
69 !| ICOORD |-->| 1: DERIVATIVE ALONG X, 2: ALONG Y
70 !| IKLE1 |-->| FIRST POINTS OF TRIANGLES
71 !| IKLE2 |-->| SECOND POINTS OF TRIANGLES
72 !| IKLE3 |-->| THIRD POINTS OF TRIANGLES
73 !| IKLE4 |-->| FOURTH POINTS OF TRIANGLES (QUADRATIC)
74 !| IKLE5 |-->| FIFTH POINTS OF TRIANGLES (QUADRATIC)
75 !| IKLE6 |-->| SIXTH POINTS OF TRIANGLES (QUADRATIC)
76 !| NELEM |-->| NUMBER OF ELEMENTS
77 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
78 !| SF |-->| STRUCTURE OF FUNCTIONS F
79 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
80 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
81 !| XMUL |-->| MULTIPLICATION FACTOR
82 !| F |-->| FUNCTION USED IN THE FORMULA
83 !| ICOORD |-->| 1: DERIVATIVE ALONG X, 2: ALONG Y
84 !| IKLE1 |-->| FIRST POINTS OF TRIANGLES
85 !| IKLE2 |-->| SECOND POINTS OF TRIANGLES
86 !| IKLE3 |-->| THIRD POINTS OF TRIANGLES
87 !| IKLE4 |-->| QUASI-BUBBLE POINT
88 !| NELEM |-->| NUMBER OF ELEMENTS
89 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
90 !| SF |-->| STRUCTURE OF FUNCTIONS F
91 !| SU |-->| BIEF_OBJ STRUCTURE OF U
92 !| SURFAC |-->| AREA OF TRIANGLES
93 !| SV |-->| BIEF_OBJ STRUCTURE OF V
94 !| U |-->| FUNCTION U USED IN THE FORMULA
95 !| V |-->| FUNCTION V USED IN THE FORMULA
96 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
97 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
98 !| XMUL |-->| MULTIPLICATION FACTOR
99 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
100 !
101  USE bief, ex_mt12bb => mt12bb
102 !
103  USE declarations_special
104  IMPLICIT NONE
105 !
106 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
107 !
108  INTEGER, INTENT(IN) :: NELEM,NELMAX,ICOORD
109  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax)
110  INTEGER, INTENT(IN) :: IKLE3(nelmax),IKLE4(nelmax)
111 !
112  DOUBLE PRECISION, INTENT(INOUT) :: A11(*),A12(*),A13(*),A14(*)
113  DOUBLE PRECISION, INTENT(INOUT) :: A21(*),A22(*),A23(*),A24(*)
114  DOUBLE PRECISION, INTENT(INOUT) :: A31(*),A32(*),A33(*),A34(*)
115  DOUBLE PRECISION, INTENT(INOUT) :: A41(*),A42(*),A43(*),A44(*)
116 !
117  DOUBLE PRECISION, INTENT(IN) :: XMUL
118  DOUBLE PRECISION, INTENT(IN) :: F(*),U(*),V(*)
119 !
120 ! STRUCTURES OF F, U, V
121  TYPE(bief_obj), INTENT(IN) :: SF,SU,SV
122 !
123  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,3),YEL(nelmax,3)
124 !
125 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
126 !
127  INTEGER IELEM,IELMF,IELMU,IELMV
128  DOUBLE PRECISION X2,X3,Y2,Y3,F1,F2,F3,F4
129  DOUBLE PRECISION U1,U2,U3,U4,V1,V2,V3,V4,UX,UY
130 !
131 !-----------------------------------------------------------------------
132 !
133  ielmf=sf%ELM
134  ielmu=su%ELM
135  ielmv=sv%ELM
136 !
137 !-----------------------------------------------------------------------
138 ! CASE WHERE F IS OF TYPE P1
139 !-----------------------------------------------------------------------
140 !
141  IF(ielmf.EQ.11) THEN
142 !
143  IF(ielmu.EQ.10.AND.ielmv.EQ.10) THEN
144 !
145 !================================
146 ! DERIVATIVE WRT X =
147 !================================
148 !
149  IF(icoord.EQ.1) THEN
150 !
151 ! LOOP ON THE ELEMENTS
152 !
153  DO ielem = 1 , nelem
154 !
155 ! INITIALISES THE GEOMETRICAL VARIABLES
156 !
157  x2 = xel(ielem,2)
158  x3 = xel(ielem,3)
159  y2 = yel(ielem,2)
160  y3 = yel(ielem,3)
161 !
162  f1 = f(ikle1(ielem))
163  f2 = f(ikle2(ielem)) - f1
164  f3 = f(ikle3(ielem)) - f1
165 !
166  ux = u(ielem)
167  uy = v(ielem)
168 !
169 ! EXTRADIAGONAL TERMS
170 !
171  a12(ielem) = ((ux*y3-2*ux*y2+2*uy*x2-uy*x3)*(f3*y2-f2*y3))
172  & *xmul/(18*(x2*y3-x3*y2))
173 !
174  a13(ielem) = ((2*ux*y3-ux*y2+uy*x2-2*uy*x3)*(f3*y2-f2*y3))
175  & *xmul/(18*(x2*y3-x3*y2))
176 !
177  a14(ielem) = ((ux*y3-ux*y2+uy*x2-uy*x3)*(f3*y2-f2*y3))
178  & *xmul/(6*(x2* y3-x3*y2))
179 !
180  a21(ielem) = (-(ux*y3+ux*y2-uy*x2-uy*x3)*(f3*y2-f2*y3))
181  & *xmul/(18*(x2*y3-x3*y2))
182 !
183  a23(ielem) = (-(2*ux*y3-ux*y2+uy*x2-2*uy*x3)*(f3*y2-f2*y3))
184  & *xmul/(18*(x2*y3-x3*y2))
185 !
186  a24(ielem) = (-(ux*y3-uy*x3)*(f3*y2-f2*y3))
187  & *xmul/(6*(x2*y3-x3*y2))
188 !
189  a31(ielem) = ((ux*y3+ux*y2-uy*x2-uy*x3)*(f3*y2-f2*y3))
190  & *xmul/(18*(x2*y3-x3*y2))
191 !
192  a32(ielem) = (-(ux*y3-2*ux*y2+2*uy*x2-uy*x3)*(f3*y2-f2*y3))
193  & *xmul/(18*(x2*y3-x3*y2))
194 !
195  a34(ielem) = ((ux*y2-uy*x2)*(f3*y2-f2*y3))
196  & *xmul/(6*(x2*y3-x3*y2))
197 !
198  a41(ielem) = (-(ux*y3-ux*y2+uy*x2-uy*x3)*(f3*y2-f2*y3))
199  & *xmul/(6*(x2*y3-x3*y2))
200 !
201  a42(ielem) = ((ux*y3-uy*x3)*(f3*y2-f2*y3))
202  & *xmul/(6*(x2*y3-x3*y2))
203 !
204  a43(ielem) = (-(ux*y2-uy*x2)*(f3*y2-f2*y3))
205  & *xmul/(6*(x2*y3-x3*y2))
206 !
207 ! DIAGONAL TERMS
208 ! THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
209 !
210  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
211  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
212  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
213  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
214 !
215  ENDDO ! IELEM
216 !
217  ELSEIF(icoord.EQ.2) THEN
218 !
219 !================================
220 ! DERIVATIVE WRT Y =
221 !================================
222 !
223  DO ielem = 1 , nelem
224 !
225 ! INITIALISES THE GEOMETRICAL VARIABLES
226 !
227  x2 = xel(ielem,2)
228  x3 = xel(ielem,3)
229  y2 = yel(ielem,2)
230  y3 = yel(ielem,3)
231 !
232  f1 = f(ikle1(ielem))
233  f2 = f(ikle2(ielem)) - f1
234  f3 = f(ikle3(ielem)) - f1
235 !
236  ux = u(ielem)
237  uy = v(ielem)
238 !
239 ! EXTRADIAGONAL TERMS
240 !
241  a12(ielem) = (-(ux*y3-2*ux*y2+2*uy*x2-uy*x3)*(x2*f3-x3*f2))
242  & *xmul/(18*(x2*y3-x3*y2))
243 !
244  a13(ielem) = (-(2*ux*y3-ux*y2+uy*x2-2*uy*x3)*(x2*f3-x3*f2))
245  & *xmul/(18*(x2*y3-x3*y2))
246 !
247  a14(ielem) = (-(ux*y3-ux*y2+uy*x2-uy*x3)*(x2*f3-x3*f2))
248  & *xmul/(6*(x2*y3-x3*y2))
249 !
250  a21(ielem) = ((ux*y3+ux*y2-uy*x2-uy*x3)*(x2*f3-x3*f2))
251  & *xmul/(18*(x2*y3-x3*y2))
252 !
253  a23(ielem) = ((2*ux*y3-ux*y2+uy*x2-2*uy*x3)*(x2*f3-x3*f2))
254  & *xmul/(18*(x2*y3-x3*y2))
255 !
256  a24(ielem) = ((ux*y3-uy*x3)*(x2*f3-x3*f2))
257  & *xmul/(6*(x2*y3-x3*y2))
258 !
259  a31(ielem) = (-(ux*y3+ux*y2-uy*x2-uy*x3)*(x2*f3-x3*f2))
260  & *xmul/(18*(x2*y3-x3*y2))
261 !
262  a32(ielem) = ((ux*y3-2*ux*y2+2*uy*x2-uy*x3)*(x2*f3-x3*f2))
263  & *xmul/(18*(x2*y3-x3*y2))
264 !
265  a34(ielem) = (-(ux*y2-uy*x2)*(x2*f3-x3*f2))
266  & *xmul/(6*(x2*y3-x3*y2))
267 !
268  a41(ielem) = ((ux*y3-ux*y2+uy*x2-uy*x3)*(x2*f3-x3*f2))
269  & *xmul/(6*(x2*y3-x3*y2))
270 !
271  a42(ielem) = (-(ux*y3-uy*x3)*(x2*f3-x3*f2))
272  & *xmul/(6*(x2*y3-x3*y2))
273 !
274  a43(ielem) = ((ux*y2-uy*x2)*(x2*f3-x3*f2))
275  & *xmul/(6*(x2*y3-x3*y2))
276 !
277 ! DIAGONAL TERMS
278 ! THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
279 !
280  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
281  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
282  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
283  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
284 !
285  ENDDO ! IELEM
286 !
287  ELSE
288 !
289  WRITE(lu,201) icoord
290  CALL plante(0)
291  stop
292  ENDIF
293 !
294 !
295  ELSEIF(ielmu.EQ.12) THEN
296 !
297 !================================
298 ! DERIVATIVE WRT X =
299 !================================
300 !
301  IF(icoord.EQ.1) THEN
302 !
303 ! LOOP ON THE ELEMENTS
304 !
305  DO ielem = 1 , nelem
306 !
307 ! INITIALISES THE GEOMETRICAL VARIABLES
308 !
309  x2 = xel(ielem,2)
310  x3 = xel(ielem,3)
311  y2 = yel(ielem,2)
312  y3 = yel(ielem,3)
313 !
314  f1 = f(ikle1(ielem))
315  f2 = f(ikle2(ielem)) - f1
316  f3 = f(ikle3(ielem)) - f1
317 !
318  u1 = u(ikle1(ielem))
319  u2 = u(ikle2(ielem))
320  u3 = u(ikle3(ielem))
321  u4 = u(ikle4(ielem))
322  v1 = v(ikle1(ielem))
323  v2 = v(ikle2(ielem))
324  v3 = v(ikle3(ielem))
325  v4 = v(ikle4(ielem))
326 !
327 ! EXTRADIAGONAL TERMS
328 !
329  a12(ielem) = (-(2*x2*v4+4*x2*v2+2*x2*v1-x3*v4-2*x3*v2-x3
330  & *v1+u4*y3-2*u4*y2+2*u2*y3-4*u2*y2+u1*y3-2*u1*y2)*(f3+
331  & f2)*(y3+y2))*xmul/(216*(x2*y3-x3*y2))
332 !
333  a13(ielem) = ((2*x2*v3+x2*v4+x2*v1-4*x3*v3-2*x3*v4-2*x3*
334  & v1+4*u3*y3-2*u3*y2+2*u4*y3-u4*y2+2*u1*y3-u1*y2)*(f3*
335  & y2-f2*y3))*xmul/(72*(x2*y3-x3*y2))
336 !
337  a14(ielem) = (3*(x2*v3+2*x2*v4+x2*v1-2*x3*v3-4*x3*v4-2*
338  & x3*v1+2*u3*y3-u3*y2+4*u4*y3-2*u4*y2+2*u1*y3-u1*y2)*(
339  & f3*y2-f2*y3)-(4*x2*v4+2*x2*v2+2*x2*v1-2*x3*v4-x3*v2-
340  & x3*v1+2*u4*y3-4*u4*y2+u2*y3-2*u2*y2+u1*y3-2*u1*y2)*(
341  & f3+f2)*(y3+y2))*xmul/(216*(x2*y3-x3*y2))
342 !
343  a21(ielem) = (-(x2*v4+x2*v2+2*x2*v1+x3*v4+x3*v2+2*x3*v1-u4
344  & *y3-u4*y2-u2*y3-u2*y2-2*u1*y3-2*u1*y2)*(f3+f2)*(y3+y2))
345  & *xmul/(216*(x2*y3-x3*y2))
346 !
347  a23(ielem) = (-(2*x2*v3+x2*v4+x2*v2-4*x3*v3-2*x3*v4-2*x3
348  & *v2+4*u3*y3-2*u3*y2+2*u4*y3-u4*y2+2*u2*y3-u2*y2)*(f3*
349  & y2-f2*y3))*xmul/(72*(x2*y3-x3*y2))
350 !
351  a24(ielem) = (-(3*(x2*v3+2*x2*v4+x2*v2-2*x3*v3-4*x3*v4-
352  & 2*x3*v2+2*u3*y3-u3*y2+4*u4*y3-2*u4*y2+2*u2*y3-u2*y2)*
353  & (f3*y2-f2*y3)+(2*x2*v4+x2*v2+x2*v1+2*x3*v4+x3*v2+x3*v1-
354  & 2*u4*y3-2*u4*y2-u2*y3-u2*y2-u1*y3-u1*y2)*(f3+f2)*(y3+y2)
355  & ))*xmul/(216*(x2*y3-x3*y2))
356 !
357  a31(ielem) = (-(x2*v3+x2*v4+2*x2*v1+x3*v3+x3*v4+2*x3*v1-u3
358  & *y3-u3*y2-u4*y3-u4*y2-2*u1*y3-2*u1*y2)*(f3*y2-f2*y3))*xmul/(
359  & 72*(x2*y3-x3*y2))
360 !
361  a32(ielem) = (-(2*x2*v3+2*x2*v4+4*x2*v2-x3*v3-x3*v4-2*x3
362  & *v2+u3*y3-2*u3*y2+u4*y3-2*u4*y2+2*u2*y3-4*u2*y2)*(f3*
363  & y2-f2*y3))*xmul/(72*(x2*y3-x3*y2))
364 !
365  a34(ielem) = (-(3*x2*v3+6*x2*v4+2*x2*v2+x2*v1-x3*v2+x3*v1
366  & -3*u3*y2-6*u4*y2+u2*y3-2*u2*y2-u1*y3-u1*y2)*(f3*y2-f2*
367  & y3))*xmul/(72*(x2*y3-x3*y2))
368 !
369  a41(ielem) = ((x2*v4+x2*v2+2*x2*v1-u4*y2-u2*y2-2*u1*y2)*(
370  & f3+f2)*(y3+y2)+3*(x3*v3+x3*v4+2*x3*v1-u3*y3-u4*y3-2*u1
371  & *y3)*(f3*y2-f2*y3))*xmul/(72*(x2*y3-x3*y2))
372 !
373  a42(ielem) = (3*(x2*v3+x2*v4+2*x2*v2-x3*v3-x3*v4-2*x3*v2+
374  & u3*y3-u3*y2+u4*y3-u4*y2+2*u2*y3-2*u2*y2)*(f3*y2-f2*y3)+
375  & (x2*v4+2*x2*v2+x2*v1-u4*y2-2*u2*y2-u1*y2)*(f3+f2)*(y3+
376  & y2))*xmul/(72*(x2*y3-x3*y2))
377 !
378  a43(ielem) = ((2*x2*v3+x2*v4+x2*v2-x3*v2+x3*v1-2*u3*y2-u4*
379  & y2+u2*y3-u2*y2-u1*y3)*(f3*y2-f2*y3))*xmul/(24*(x2*y3-x3*y2))
380 !
381 ! DIAGONAL TERMS
382 ! THE SUM OF EACH COLUMN IS 0
383 !
384  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
385  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
386  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
387  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
388 !
389  ENDDO ! IELEM
390 !
391  ELSEIF(icoord.EQ.2) THEN
392 !
393 !================================
394 ! DERIVATIVE WRT Y =
395 !================================
396 !
397  DO ielem = 1 , nelem
398 !
399 ! INITIALISES THE GEOMETRICAL VARIABLES
400 !
401  x2 = xel(ielem,2)
402  x3 = xel(ielem,3)
403  y2 = yel(ielem,2)
404  y3 = yel(ielem,3)
405 !
406  f1 = f(ikle1(ielem))
407  f2 = f(ikle2(ielem)) - f1
408  f3 = f(ikle3(ielem)) - f1
409 !
410  u1 = u(ikle1(ielem))
411  u2 = u(ikle2(ielem))
412  u3 = u(ikle3(ielem))
413  u4 = u(ikle4(ielem))
414  v1 = v(ikle1(ielem))
415  v2 = v(ikle2(ielem))
416  v3 = v(ikle3(ielem))
417  v4 = v(ikle4(ielem))
418 !
419 ! EXTRADIAGONAL TERMS
420 !
421  a12(ielem) = (2*x2*v4+4*x2*v2+2*x2*v1-x3*v4-2*x3*v2-x3*
422  & v1+u4*y3-2*u4*y2+2*u2*y3-4*u2*y2+u1*y3-2*u1*y2)*(x2+
423  & x3)*(f3+f2)*xmul/(216*(x2*y3-x3*y2))
424 !
425  a13(ielem) = -(2*x2*v3+x2*v4+x2*v1-4*x3*v3-2*x3*v4-2*x3
426  & *v1+4*u3*y3-2*u3*y2+2*u4*y3-u4*y2+2*u1*y3-u1*y2)*(x2*
427  & f3-x3*f2)*xmul/(72*(x2*y3-x3*y2))
428 !
429  a14(ielem) = -(3*(x2*v3+2*x2*v4+x2*v1-2*x3*v3-4*x3*v4-
430  & 2*x3*v1+2*u3*y3-u3*y2+4*u4*y3-2*u4*y2+2*u1*y3-u1*y2)*
431  & (x2*f3-x3*f2)-(4*x2*v4+2*x2*v2+2*x2*v1-2*x3*v4-x3*v2-
432  & x3*v1+2*u4*y3-4*u4*y2+u2*y3-2*u2*y2+u1*y3-2*u1*y2)*(
433  & x2+x3)*(f3+f2))*xmul/(216*(x2*y3-x3*y2))
434 !
435  a21(ielem) = (x2*v4+x2*v2+2*x2*v1+x3*v4+x3*v2+2*x3*v1-u4*
436  & y3-u4*y2-u2*y3-u2*y2-2*u1*y3-2*u1*y2)*(x2+x3)*(f3+f2)*xmul/
437  & (216*(x2*y3-x3*y2))
438 !
439  a23(ielem) = (2*x2*v3+x2*v4+x2*v2-4*x3*v3-2*x3*v4-2*x3*
440  & v2+4*u3*y3-2*u3*y2+2*u4*y3-u4*y2+2*u2*y3-u2*y2)*(x2*
441  & f3-x3*f2)*xmul/(72*(x2*y3-x3*y2))
442 !
443  a24(ielem) = (3*(x2*v3+2*x2*v4+x2*v2-2*x3*v3-4*x3*v4-2*
444  & x3*v2+2*u3*y3-u3*y2+4*u4*y3-2*u4*y2+2*u2*y3-u2*y2)*(
445  & x2*f3-x3*f2)+(2*x2*v4+x2*v2+x2*v1+2*x3*v4+x3*v2+x3*v1-
446  & 2*u4*y3-2*u4*y2-u2*y3-u2*y2-u1*y3-u1*y2)*(x2+x3)*(f3+f2)
447  & )*xmul/(216*(x2*y3-x3*y2))
448 !
449  a31(ielem) = (x2*v3+x2*v4+2*x2*v1+x3*v3+x3*v4+2*x3*v1-u3*
450  & y3-u3*y2-u4*y3-u4*y2-2*u1*y3-2*u1*y2)*(x2*f3-x3*f2)*xmul/(
451  & 72*(x2*y3-x3*y2))
452 !
453  a32(ielem) = (2*x2*v3+2*x2*v4+4*x2*v2-x3*v3-x3*v4-2*x3*
454  & v2+u3*y3-2*u3*y2+u4*y3-2*u4*y2+2*u2*y3-4*u2*y2)*(x2*
455  & f3-x3*f2)*xmul/(72*(x2*y3-x3*y2))
456 !
457  a34(ielem) = (3*x2*v3+6*x2*v4+2*x2*v2+x2*v1-x3*v2+x3*v1-
458  & 3*u3*y2-6*u4*y2+u2*y3-2*u2*y2-u1*y3-u1*y2)*(x2*f3-x3*f2
459  & )*xmul/(72*(x2*y3-x3*y2))
460 !
461  a41(ielem) = -((x2*v4+x2*v2+2*x2*v1-u4*y2-u2*y2-2*u1*y2)*
462  & (x2+x3)*(f3+f2)+3*(x2*f3-x3*f2)*(x3*v3+x3*v4+2*x3*v1-u3
463  & *y3-u4*y3-2*u1*y3))*xmul/(72*(x2*y3-x3*y2))
464 !
465  a42(ielem) = -(3*(x2*v3+x2*v4+2*x2*v2-x3*v3-x3*v4-2*x3*
466  & v2+u3*y3-u3*y2+u4*y3-u4*y2+2*u2*y3-2*u2*y2)*(x2*f3-x3*
467  & f2)+(x2*v4+2*x2*v2+x2*v1-u4*y2-2*u2*y2-u1*y2)*(x2+x3)*(
468  & f3+f2))*xmul/(72*(x2*y3-x3*y2))
469 !
470  a43(ielem) = -(2*x2*v3+x2*v4+x2*v2-x3*v2+x3*v1-2*u3*y2-u4
471  & *y2+u2*y3-u2*y2-u1*y3)*(x2*f3-x3*f2)*xmul/(24*(x2*y3-x3*y2))
472 !
473 ! DIAGONAL TERMS
474 ! THE SUM OF EACH COLUMN IS 0
475 !
476  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
477  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
478  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
479  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
480 !
481  ENDDO ! IELEM
482 !
483  ELSE
484 !
485  WRITE(lu,201) icoord
486  CALL plante(0)
487  stop
488  ENDIF
489 !
490  ELSE
491 !
492  WRITE(lu,301) ielmu
493  CALL plante(0)
494  stop
495 !
496  ENDIF
497 !
498 !-----------------------------------------------------------------------
499 ! CASE WHERE F IS OF TYPE QUASI-BUBBLE
500 !-----------------------------------------------------------------------
501 !
502  ELSEIF(ielmf.EQ.12) THEN
503 !
504  IF (ielmu.EQ.10) THEN
505 !
506 !================================
507 ! DERIVATIVE WRT X =
508 !================================
509 !
510  IF(icoord.EQ.1) THEN
511 !
512 ! LOOP ON THE ELEMENTS
513 !
514  DO ielem = 1 , nelem
515 !
516 ! INITIALISES THE GEOMETRICAL VARIABLES
517 !
518  x2 = xel(ielem,2)
519  x3 = xel(ielem,3)
520  y2 = yel(ielem,2)
521  y3 = yel(ielem,3)
522 !
523  f1 = f(ikle1(ielem))
524  f2 = f(ikle2(ielem)) - f1
525  f3 = f(ikle3(ielem)) - f1
526  f4 = f(ikle4(ielem)) - f1
527 !
528  ux = u(ielem)
529  uy = v(ielem)
530 !
531 ! EXTRADIAGONAL TERMS
532 !
533  a12(ielem) = ((ux*y3-2*ux*y2+2*uy*x2-uy*x3)*(3*f4*y2-f2*y3-f2
534  & *y2))*xmul/(18*(x2*y3-x3*y2))
535 !
536  a13(ielem) = ((2*ux*y3-ux*y2+uy*x2-2*uy*x3)*(f3*y3+f3*y2-3*f4
537  & *y3))*xmul/(18*(x2*y3-x3*y2))
538 !
539  a14(ielem) = ((2*ux*y3-ux*y2+uy*x2-2*uy*x3)*(f3*y3+f3*y2-3*f4
540  & *y3)+(ux*y3-2*ux*y2+2*uy*x2-uy*x3)*(3*f4*y2-f2*y3-f2*y2))*xmul/
541  & (18*(x2*y3-x3*y2))
542 !
543  a21(ielem) = (-(ux*y3+ux*y2-uy*x2-uy*x3)*(3*f4*y2-f2*y3-f2*y2))
544  & *xmul/(18*(x2*y3-x3*y2))
545 !
546  a23(ielem) = ((2*ux*y3-ux*y2+uy*x2-2*uy*x3)*(f3*y3-2*f3*y2-3
547  & *f4*y3+3*f4*y2+2*f2*y3-f2*y2))*xmul/(18*(x2*y3-x3*y2))
548 !
549  a24(ielem) = ((2*ux*y3-ux*y2+uy*x2-2*uy*x3)*(f3*y3-2*f3*y2-3
550  & *f4*y3+3*f4*y2+2*f2*y3-f2*y2)-(ux*y3+ux*y2-uy*x2-uy*x3)*(3
551  & *f4*y2-f2*y3-f2*y2))*xmul/(18*(x2*y3-x3*y2))
552 !
553  a31(ielem) = ((ux*y3+ux*y2-uy*x2-uy*x3)*(f3*y3+f3*y2-3*f4*y3))
554  & *xmul/(18*(x2*y3-x3*y2))
555 !
556  a32(ielem) = ((ux*y3-2*ux*y2+2*uy*x2-uy*x3)*(f3*y3-2*f3*y2-3
557  & *f4*y3+3*f4*y2+2*f2*y3-f2*y2))*xmul/(18*(x2*y3-x3*y2))
558 !
559  a34(ielem) = ((ux*y3+ux*y2-uy*x2-uy*x3)*(f3*y3+f3*y2-3*f4*y3)+(
560  & ux*y3-2*ux*y2+2*uy*x2-uy*x3)*(f3*y3-2*f3*y2-3*f4*y3+3*f4
561  & *y2+2*f2*y3-f2*y2))*xmul/(18*(x2*y3-x3*y2))
562 !
563  a41(ielem) = (-((ux*y3-uy*x3)*(f3*y3+f3*y2-3*f4*y3)-(ux*y2-uy*
564  & x2)*(3*f4*y2-f2*y3-f2*y2)))*xmul/(6*(x2*y3-x3*y2))
565 !
566  a42(ielem) = (-((ux*y3-ux*y2+uy*x2-uy*x3)*(f3*y3-2*f3*y2-3*f4*
567  & y3+3*f4*y2+2*f2*y3-f2*y2)-(ux*y2-uy*x2)*(3*f4*y2-f2*y3-
568  & f2*y2)))*xmul/(6*(x2*y3-x3*y2))
569 !
570  a43(ielem) = (-((ux*y3-ux*y2+uy*x2-uy*x3)*(f3*y3-2*f3*y2-3*f4*
571  & y3+3*f4*y2+2*f2*y3-f2*y2)+(ux*y3-uy*x3)*(f3*y3+f3*y2-3*
572  & f4*y3)))*xmul/(6*(x2*y3-x3*y2))
573 !
574 !
575 ! DIAGONAL TERMS
576 ! THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
577 !
578  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
579  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
580  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
581  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
582 !
583  ENDDO ! IELEM
584 !
585  ELSEIF(icoord.EQ.2) THEN
586 !
587 !================================
588 ! DERIVATIVE WRT Y =
589 !================================
590 !
591  DO ielem = 1 , nelem
592 !
593 ! INITIALISES THE GEOMETRICAL VARIABLES
594 !
595  x2 = xel(ielem,2)
596  x3 = xel(ielem,3)
597  y2 = yel(ielem,2)
598  y3 = yel(ielem,3)
599 !
600  f1 = f(ikle1(ielem))
601  f2 = f(ikle2(ielem)) - f1
602  f3 = f(ikle3(ielem)) - f1
603  f4 = f(ikle4(ielem)) - f1
604 !
605  ux = u(ielem)
606  uy = v(ielem)
607 !
608 ! EXTRADIAGONAL TERMS
609 !
610  a12(ielem) = (-(ux*y3-2*ux*y2+2*uy*x2-uy*x3)*(3*x2*f4-x2*f2-
611  & x3*f2))*xmul/(18*(x2*y3-x3*y2))
612 !
613  a13(ielem) = (-(2*ux*y3-ux*y2+uy*x2-2*uy*x3)*(x2*f3+x3*f3-3*
614  & x3*f4))*xmul/(18*(x2*y3-x3*y2))
615 !
616  a14(ielem) = (-((2*ux*y3-ux*y2+uy*x2-2*uy*x3)*(x2*f3+x3*f3-3*
617  & x3*f4)+(ux*y3-2*ux*y2+2*uy*x2-uy*x3)*(3*x2*f4-x2*f2-x3*f2)
618  & ))*xmul/(18*(x2*y3-x3*y2))
619 !
620  a21(ielem) = ((ux*y3+ux*y2-uy*x2-uy*x3)*(3*x2*f4-x2*f2-x3*f2))
621  & *xmul/(18*(x2*y3-x3*y2))
622 !
623  a23(ielem) = ((2*ux*y3-ux*y2+uy*x2-2*uy*x3)*(2*x2*f3-3*x2*f4
624  & +x2*f2-x3*f3+3*x3*f4-2*x3*f2))*xmul/(18*(x2*y3-x3*y2))
625 !
626  a24(ielem) = ((2*ux*y3-ux*y2+uy*x2-2*uy*x3)*(2*x2*f3-3*x2*f4
627  & +x2*f2-x3*f3+3*x3*f4-2*x3*f2)+(ux*y3+ux*y2-uy*x2-uy*x3)*(3
628  & *x2*f4-x2*f2-x3*f2))*xmul/(18*(x2*y3-x3*y2))
629 !
630  a31(ielem) = (-(ux*y3+ux*y2-uy*x2-uy*x3)*(x2*f3+x3*f3-3*x3*f4))
631  & *xmul/(18*(x2*y3-x3*y2))
632 !
633  a32(ielem) = ((ux*y3-2*ux*y2+2*uy*x2-uy*x3)*(2*x2*f3-3*x2*f4
634  & +x2*f2-x3*f3+3*x3*f4-2*x3*f2))*xmul/(18*(x2*y3-x3*y2))
635 !
636  a34(ielem) = (-((ux*y3+ux*y2-uy*x2-uy*x3)*(x2*f3+x3*f3-3*x3*f4)
637  & -(ux*y3-2*ux*y2+2*uy*x2-uy*x3)*(2*x2*f3-3*x2*f4+x2*f2-x3*
638  & f3+3*x3*f4-2*x3*f2)))*xmul/(18*(x2*y3-x3*y2))
639 !
640  a41(ielem) = ((ux*y3-uy*x3)*(x2*f3+x3*f3-3*x3*f4)-(ux*y2-uy*x2)
641  & *(3*x2*f4-x2*f2-x3*f2))*xmul/(6*(x2*y3-x3*y2))
642 !
643  a42(ielem) = (-((ux*y3-ux*y2+uy*x2-uy*x3)*(2*x2*f3-3*x2*f4+x2*
644  & f2-x3*f3+3*x3*f4-2*x3*f2)+(ux*y2-uy*x2)*(3*x2*f4-x2*f2-
645  & x3*f2)))*xmul/(6*(x2*y3-x3*y2))
646 !
647  a43(ielem) = (-((ux*y3-ux*y2+uy*x2-uy*x3)*(2*x2*f3-3*x2*f4+x2*
648  & f2-x3*f3+3*x3*f4-2*x3*f2)-(ux*y3-uy*x3)*(x2*f3+x3*f3-3*
649  & x3*f4)))*xmul/(6*(x2*y3-x3*y2))
650 !
651 !
652 ! DIAGONAL TERMS
653 ! THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
654 !
655  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
656  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
657  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
658  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
659 !
660  ENDDO ! IELEM
661 !
662  ELSE
663 !
664  WRITE(lu,201) icoord
665  CALL plante(0)
666  stop
667 !
668  ENDIF
669 !
670 !
671  ELSEIF(ielmu.EQ.12) THEN
672 !
673 !================================
674 ! DERIVATIVE WRT X =
675 !================================
676 !
677  IF(icoord.EQ.1) THEN
678 !
679 ! LOOP ON THE ELEMENTS
680 !
681  DO ielem = 1 , nelem
682 !
683 ! INITIALISES THE GEOMETRICAL VARIABLES
684 !
685  x2 = xel(ielem,2)
686  x3 = xel(ielem,3)
687  y2 = yel(ielem,2)
688  y3 = yel(ielem,3)
689 !
690  f1 = f(ikle1(ielem))
691  f2 = f(ikle2(ielem)) - f1
692  f3 = f(ikle3(ielem)) - f1
693  f4 = f(ikle4(ielem)) - f1
694 !
695  u1 = u(ikle1(ielem))
696  u2 = u(ikle2(ielem))
697  u3 = u(ikle3(ielem))
698  u4 = u(ikle4(ielem))
699  v1 = v(ikle1(ielem))
700  v2 = v(ikle2(ielem))
701  v3 = v(ikle3(ielem))
702  v4 = v(ikle4(ielem))
703 !
704 ! EXTRADIAGONAL TERMS
705 !
706  a12(ielem) = (-(2*x2*v4+4*x2*v2+2*x2*v1-x3*v4-2*x3*v2-x3
707  & *v1+u4*y3-2*u4*y2+2*u2*y3-4*u2*y2+u1*y3-2*u1*y2)*(y3+
708  & y2)*f4)*xmul/(72*(x2*y3-x3*y2))
709 !
710  a13(ielem) = ((2*x2*v3+x2*v4+x2*v1-4*x3*v3-2*x3*v4-2*x3*
711  & v1+4*u3*y3-2*u3*y2+2*u4*y3-u4*y2+2*u1*y3-u1*y2)*(f3*
712  & y3+f3*y2-3*f4*y3))*xmul/(72*(x2*y3-x3*y2))
713 !
714  a14(ielem) = ((x2*v3+2*x2*v4+x2*v1-2*x3*v3-4*x3*v4-2*x3*
715  & v1+2*u3*y3-u3*y2+4*u4*y3-2*u4*y2+2*u1*y3-u1*y2)*(f3*
716  & y3+f3*y2-3*f4*y3)-(4*x2*v4+2*x2*v2+2*x2*v1-2*x3*v4-
717  & x3*v2-x3*v1+2*u4*y3-4*u4*y2+u2*y3-2*u2*y2+u1*y3-2*u1*
718  & y2)*(y3+y2)*f4)*xmul/(72*(x2*y3-x3*y2))
719 !
720  a21(ielem) = (-(x2*v4+x2*v2+2*x2*v1+x3*v4+x3*v2+2*x3*v1-u4
721  & *y3-u4*y2-u2*y3-u2*y2-2*u1*y3-2*u1*y2)*(y3+y2)*f4)*xmul/(72
722  & *(x2*y3-x3*y2))
723 !
724  a23(ielem) = ((2*x2*v3+x2*v4+x2*v2-4*x3*v3-2*x3*v4-2*x3*
725  & v2+4*u3*y3-2*u3*y2+2*u4*y3-u4*y2+2*u2*y3-u2*y2)*(f3*
726  & y3-2*f3*y2-3*f4*y3+3*f4*y2+2*f2*y3-f2*y2))*xmul/(72*(x2*
727  & y3-x3*y2))
728 !
729  a24(ielem) = ((x2*v3+2*x2*v4+x2*v2-2*x3*v3-4*x3*v4-2*x3*
730  & v2+2*u3*y3-u3*y2+4*u4*y3-2*u4*y2+2*u2*y3-u2*y2)*(f3*
731  & y3-2*f3*y2-3*f4*y3+3*f4*y2+2*f2*y3-f2*y2)-(2*x2*v4+
732  & x2*v2+x2*v1+2*x3*v4+x3*v2+x3*v1-2*u4*y3-2*u4*y2-u2*y3-
733  & u2*y2-u1*y3-u1*y2)*(y3+y2)*f4)*xmul/(72*(x2*y3-x3*y2))
734 !
735  a31(ielem) = (-(x2*v3+x2*v4+2*x2*v1+x3*v3+x3*v4+2*x3*v1-u3
736  & *y3-u3*y2-u4*y3-u4*y2-2*u1*y3-2*u1*y2)*(f3*y3+f3*y2-3*
737  & f4*y3))*xmul/(72*(x2*y3-x3*y2))
738 !
739  a32(ielem) = ((2*x2*v3+2*x2*v4+4*x2*v2-x3*v3-x3*v4-2*x3*
740  & v2+u3*y3-2*u3*y2+u4*y3-2*u4*y2+2*u2*y3-4*u2*y2)*(f3*
741  & y3-2*f3*y2-3*f4*y3+3*f4*y2+2*f2*y3-f2*y2))*xmul/(72*(x2*
742  & y3-x3*y2))
743 !
744  a34(ielem) = ((2*x2*v3+4*x2*v4+2*x2*v2-x3*v3-2*x3*v4-x3*
745  & v2+u3*y3-2*u3*y2+2*u4*y3-4*u4*y2+u2*y3-2*u2*y2)*(f3*
746  & y3-2*f3*y2-3*f4*y3+3*f4*y2+2*f2*y3-f2*y2)-(x2*v3+2*
747  & x2*v4+x2*v1+x3*v3+2*x3*v4+x3*v1-u3*y3-u3*y2-2*u4*y3-2*
748  & u4*y2-u1*y3-u1*y2)*(f3*y3+f3*y2-3*f4*y3))*xmul/(72*(x2*y3-x3
749  & *y2))
750 !
751  a41(ielem) = ((x2*v4+x2*v2+2*x2*v1-u4*y2-u2*y2-2*u1*y2)*(
752  & y3+y2)*f4+(x3*v3+x3*v4+2*x3*v1-u3*y3-u4*y3-2*u1*y3)*(f3
753  & *y3+f3*y2-3*f4*y3))*xmul/(24*(x2*y3-x3*y2))
754 !
755  a42(ielem) = (-((x2*v3+x2*v4+2*x2*v2-x3*v3-x3*v4-2*x3*v2+
756  & u3*y3-u3*y2+u4*y3-u4*y2+2*u2*y3-2*u2*y2)*(f3*y3-2*f3*
757  & y2-3*f4*y3+3*f4*y2+2*f2*y3-f2*y2)-(x2*v4+2*x2*v2+x2*
758  & v1-u4*y2-2*u2*y2-u1*y2)*(y3+y2)*f4))*xmul/(24*(x2*y3-x3*y2))
759 !
760  a43(ielem) = (-((2*x2*v3+x2*v4+x2*v2-2*x3*v3-x3*v4-x3*v2+
761  & 2*u3*y3-2*u3*y2+u4*y3-u4*y2+u2*y3-u2*y2)*(f3*y3-2*f3*y2
762  & -3*f4*y3+3*f4*y2+2*f2*y3-f2*y2)-(2*x3*v3+x3*v4+x3*v1-
763  & 2*u3*y3-u4*y3-u1*y3)*(f3*y3+f3*y2-3*f4*y3)))*xmul/(24*(x2*y3
764  & -x3*y2))
765 !
766 !
767 ! DIAGONAL TERMS
768 ! THE SUM OF EACH COLUMN IS 0
769 !
770  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
771  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
772  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
773  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
774 !
775  ENDDO ! IELEM
776 !
777  ELSEIF(icoord.EQ.2) THEN
778 !
779 !================================
780 ! DERIVATIVE WRT Y =
781 !================================
782 !
783  DO ielem = 1 , nelem
784 !
785 ! INITIALISES THE GEOMETRICAL VARIABLES
786 !
787  x2 = xel(ielem,2)
788  x3 = xel(ielem,3)
789  y2 = yel(ielem,2)
790  y3 = yel(ielem,3)
791 !
792  f1 = f(ikle1(ielem))
793  f2 = f(ikle2(ielem)) - f1
794  f3 = f(ikle3(ielem)) - f1
795  f4 = f(ikle4(ielem)) - f1
796 !
797  u1 = u(ikle1(ielem))
798  u2 = u(ikle2(ielem))
799  u3 = u(ikle3(ielem))
800  u4 = u(ikle4(ielem))
801  v1 = v(ikle1(ielem))
802  v2 = v(ikle2(ielem))
803  v3 = v(ikle3(ielem))
804  v4 = v(ikle4(ielem))
805 !
806 ! EXTRADIAGONAL TERMS
807 !
808  a12(ielem) = ((2*x2*v4+4*x2*v2+2*x2*v1-x3*v4-2*x3*v2-x3*
809  & v1+u4*y3-2*u4*y2+2*u2*y3-4*u2*y2+u1*y3-2*u1*y2)*(x2+
810  & x3)*f4)*xmul/(72*(x2*y3-x3*y2))
811 !
812  a13(ielem) = (-(x2*f3+x3*f3-3*x3*f4)*(2*x2*v3+x2*v4+x2*v1-
813  & 4*x3*v3-2*x3*v4-2*x3*v1+4*u3*y3-2*u3*y2+2*u4*y3-u4*
814  & y2+2*u1*y3-u1*y2))*xmul/(72*(x2*y3-x3*y2))
815 !
816  a14(ielem) = (-((x2*f3+x3*f3-3*x3*f4)*(x2*v3+2*x2*v4+x2*v1
817  & -2*x3*v3-4*x3*v4-2*x3*v1+2*u3*y3-u3*y2+4*u4*y3-2*u4
818  & *y2+2*u1*y3-u1*y2)-(4*x2*v4+2*x2*v2+2*x2*v1-2*x3*v4-
819  & x3*v2-x3*v1+2*u4*y3-4*u4*y2+u2*y3-2*u2*y2+u1*y3-2*u1*
820  & y2)*(x2+x3)*f4))*xmul/(72*(x2*y3-x3*y2))
821 !
822  a21(ielem) = ((x2*v4+x2*v2+2*x2*v1+x3*v4+x3*v2+2*x3*v1-u4*
823  & y3-u4*y2-u2*y3-u2*y2-2*u1*y3-2*u1*y2)*(x2+x3)*f4)*xmul/(72*
824  & (x2*y3-x3*y2))
825 !
826  a23(ielem) = ((2*x2*f3-3*x2*f4+x2*f2-x3*f3+3*x3*f4-2*x3*
827  & f2)*(2*x2*v3+x2*v4+x2*v2-4*x3*v3-2*x3*v4-2*x3*v2+4*
828  & u3*y3-2*u3*y2+2*u4*y3-u4*y2+2*u2*y3-u2*y2))*xmul/(72*(x2*
829  & y3-x3*y2))
830 !
831  a24(ielem) = ((2*x2*f3-3*x2*f4+x2*f2-x3*f3+3*x3*f4-2*x3*
832  & f2)*(x2*v3+2*x2*v4+x2*v2-2*x3*v3-4*x3*v4-2*x3*v2+2*
833  & u3*y3-u3*y2+4*u4*y3-2*u4*y2+2*u2*y3-u2*y2)+(2*x2*v4+
834  & x2*v2+x2*v1+2*x3*v4+x3*v2+x3*v1-2*u4*y3-2*u4*y2-u2*y3-
835  & u2*y2-u1*y3-u1*y2)*(x2+x3)*f4)*xmul/(72*(x2*y3-x3*y2))
836 !
837  a31(ielem) = ((x2*f3+x3*f3-3*x3*f4)*(x2*v3+x2*v4+2*x2*v1+
838  & x3*v3+x3*v4+2*x3*v1-u3*y3-u3*y2-u4*y3-u4*y2-2*u1*y3-2*
839  & u1*y2))*xmul/(72*(x2*y3-x3*y2))
840 !
841  a32(ielem) = ((2*x2*f3-3*x2*f4+x2*f2-x3*f3+3*x3*f4-2*x3*
842  & f2)*(2*x2*v3+2*x2*v4+4*x2*v2-x3*v3-x3*v4-2*x3*v2+u3*
843  & y3-2*u3*y2+u4*y3-2*u4*y2+2*u2*y3-4*u2*y2))*xmul/(72*(x2*
844  & y3-x3*y2))
845 !
846  a34(ielem) = ((2*x2*f3-3*x2*f4+x2*f2-x3*f3+3*x3*f4-2*x3*
847  & f2)*(2*x2*v3+4*x2*v4+2*x2*v2-x3*v3-2*x3*v4-x3*v2+u3*
848  & y3-2*u3*y2+2*u4*y3-4*u4*y2+u2*y3-2*u2*y2)+(x2*f3+x3*
849  & f3-3*x3*f4)*(x2*v3+2*x2*v4+x2*v1+x3*v3+2*x3*v4+x3*v1-
850  & u3*y3-u3*y2-2*u4*y3-2*u4*y2-u1*y3-u1*y2))*xmul/(72*(x2*y3-
851  & x3*y2))
852 !
853  a41(ielem) = (-((x2*f3+x3*f3-3*x3*f4)*(x3*v3+x3*v4+2*x3*v1
854  & -u3*y3-u4*y3-2*u1*y3)+(x2*v4+x2*v2+2*x2*v1-u4*y2-u2*y2-
855  & 2*u1*y2)*(x2+x3)*f4))*xmul/(24*(x2*y3-x3*y2))
856 !
857  a42(ielem) = (-((2*x2*f3-3*x2*f4+x2*f2-x3*f3+3*x3*f4-2*
858  & x3*f2)*(x2*v3+x2*v4+2*x2*v2-x3*v3-x3*v4-2*x3*v2+u3*y3-
859  & u3*y2+u4*y3-u4*y2+2*u2*y3-2*u2*y2)+(x2*v4+2*x2*v2+x2*
860  & v1-u4*y2-2*u2*y2-u1*y2)*(x2+x3)*f4))*xmul/(24*(x2*y3-x3*y2))
861 !
862  a43(ielem) = (-((2*x2*f3-3*x2*f4+x2*f2-x3*f3+3*x3*f4-2*
863  & x3*f2)*(2*x2*v3+x2*v4+x2*v2-2*x3*v3-x3*v4-x3*v2+2*u3*
864  & y3-2*u3*y2+u4*y3-u4*y2+u2*y3-u2*y2)+(x2*f3+x3*f3-3*x3*
865  & f4)*(2*x3*v3+x3*v4+x3*v1-2*u3*y3-u4*y3-u1*y3)))*xmul/(24*(
866  & x2*y3-x3*y2))
867 !
868 !
869 ! DIAGONAL TERMS
870 ! THE SUM OF EACH COLUMN IS 0
871 !
872  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
873  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
874  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
875  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
876 !
877  ENDDO ! IELEM
878 !
879  ELSE
880 !
881  WRITE(lu,201) icoord
882  CALL plante(0)
883  stop
884  ENDIF
885 !
886  ELSE
887 !
888  WRITE(lu,301) ielmu
889  CALL plante(0)
890  stop
891  ENDIF
892 !
893 !-----------------------------------------------------------------------
894 !
895  ELSE
896  WRITE(lu,101) ielmf
897 101 FORMAT(1x,'MT12BB (BIEF) :',/,
898  & 1x,'DISCRETIZATION OF F : ',1i6,' NOT AVAILABLE')
899  CALL plante(0)
900  stop
901  ENDIF
902 !
903 201 FORMAT(1x,'MT12BB (BIEF) : IMPOSSIBLE COMPONENT ',
904  & 1i6,' CHECK ICOORD')
905 !
906 301 FORMAT(1x,'MT12BB (BIEF) :',/,
907  & 1x,'DISCRETIZATION OF U : ',1i6,' NOT AVAILABLE')
908 !
909 !-----------------------------------------------------------------------
910 !
911  RETURN
912  END
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
mt12bb
subroutine mt12bb(A11, A12, A13, A14, A21, A22, A23, A24, A31, A32, A33, A34, A41, A42, A43, A44, XMUL, SF, SU, SV, F, U, V, XEL, YEL, IKLE1, IKLE2, IKLE3, IKLE4, NELEM, NELMAX, ICOORD)
Definition: mt12bb.f:14
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:13 for The TELEMAC-MASCARET system by doxygen 1.8.13