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