The TELEMAC-MASCARET system  trunk
mt04tt.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt04tt
3 ! *****************
4 !
5  &( t,xm,xmul,su,sv,sw,u,v,w,x,y,z,ikle,nelem,nelmax,formul)
6 !
7 !***********************************************************************
8 ! BIEF V6P2 21/08/2010
9 !***********************************************************************
10 !
11 !brief COMPUTES A MATRIX FOR THE SUPG METHOD.
12 !code
13 !+ COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
14 !+
15 !+ / -> ->
16 !+ A = XMUL / U . GRAD(P ) * U . GRAD(P ) * J(X,Y) DXDY
17 !+ I J /S I J
18 !+
19 !+ BY ELEMENTARY CELL; THE ELEMENT IS THE TETRAHEDRON WITH LINEAR INTERPOLATION
20 !+
21 !+ J(X,Y): JACOBIAN OF THE ISOPARAMETRIC TRANSFORMATION
22 !
23 !warning THE VERTICAL COMPONENT IS HERE NEGLECTED,
24 !+ IN ACCORDANCE WITH THE PRINCIPLE NOTE,
25 !+ AND CONTRARY TO THE PRISMS
26 !code
27 !+ STORAGE CONVENTION FOR EXTRA-DIAGONAL TERMS:
28 !+
29 !+ XM(IELEM, 1) ----> M(1,2) = M(2,1)
30 !+ XM(IELEM, 2) ----> M(1,3) = M(3,1)
31 !+ XM(IELEM, 3) ----> M(1,4) = M(4,1)
32 !+ XM(IELEM, 4) ----> M(2,3) = M(3,2)
33 !+ XM(IELEM, 5) ----> M(2,4) = M(4,2)
34 !+ XM(IELEM, 6) ----> M(3,4) = M(4,3)
35 !
36 !history J-M HERVOUET (LNHE) ; F LEPEINTRE (LNH)
37 !+ 22/08/05
38 !+ V5P6
39 !+
40 !
41 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
42 !+ 13/07/2010
43 !+ V6P0
44 !+ Translation of French comments within the FORTRAN sources into
45 !+ English comments
46 !
47 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
48 !+ 21/08/2010
49 !+ V6P0
50 !+ Creation of DOXYGEN tags for automated documentation and
51 !+ cross-referencing of the FORTRAN sources
52 !
53 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
54 !| IKLE |-->| CONNECTIVITY TABLE.
55 !| NELEM |-->| NUMBER OF ELEMENTS
56 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
57 !| SU |-->| STRUCTURE OF FUNCTIONS U
58 !| SV |-->| STRUCTURE OF FUNCTIONS V
59 !| SW |-->| STRUCTURE OF FUNCTIONS W
60 !| SURFAC |-->| AREA OF 2D ELEMENTS
61 !| T |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
62 !| U |-->| FUNCTION USED IN THE FORMULA
63 !| V |-->| FUNCTION USED IN THE FORMULA
64 !| W |-->| FUNCTION USED IN THE FORMULA
65 !| X |-->| ABSCISSAE OF POINTS
66 !| Y |-->| ORDINATES OF POINTS
67 !| Z |-->| ELEVATIONS OF POINTS
68 !| XM |<->| OFF-DIAGONAL TERMS
69 !| XMUL |-->| COEFFICIENT FOR MULTIPLICATION
70 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
71 !
72  USE bief, ex_mt04tt => mt04tt
74  IMPLICIT NONE
75 !
76 !
77 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
78 !
79  INTEGER, INTENT(IN) :: NELEM,NELMAX
80  INTEGER, INTENT(IN) :: IKLE(nelmax,4)
81 !
82  DOUBLE PRECISION, INTENT(INOUT) :: T(nelmax,4),XM(nelmax,6)
83 !
84  DOUBLE PRECISION, INTENT(IN) :: XMUL
85  DOUBLE PRECISION, INTENT(IN) :: U(*),V(*),W(*)
86 !
87 ! STRUCTURES OF U, V, W
88 !
89  TYPE(bief_obj), INTENT(IN) :: SU,SV,SW
90 !
91  DOUBLE PRECISION, INTENT(IN) :: X(*),Y(*),Z(*)
92 !
93  CHARACTER(LEN=16), INTENT(IN) :: FORMUL
94 !
95 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
96 !
97 ! SPECIFIC DECLARATIONS
98 !
99  DOUBLE PRECISION X2,Y2,Z2,X3,Y3,Z3,X4,Y4,Z4
100  DOUBLE PRECISION U1,U2,U3,U4,V1,V2,V3,V4,W1,W2,W3,W4,XJAC
101  DOUBLE PRECISION COEF,SUR120,AUX1,AUX2,AUX4,AUX5,AUX6,AUX7
102  DOUBLE PRECISION A1,A2,A3
103  INTEGER I1,I2,I3,I4,IELEM
104 !
105 !***********************************************************************
106 !
107  sur120=xmul/120.d0
108 !
109  IF(formul(1:7).EQ.'MAUGUG2') THEN
110 !
111  IF((su%ELM.EQ.31.AND.sv%ELM.EQ.31.AND.sw%ELM.EQ.31).OR.
112  & (su%ELM.EQ.51.AND.sv%ELM.EQ.51.AND.sw%ELM.EQ.51) ) THEN
113 !
114 !-----------------------------------------------------------------------
115 !
116 ! LINEAR DISCRETISATION OF DIFFUSION COEFFICIENTS
117 !
118 ! LOOP ON THE TETRAHEDRONS
119 !
120  DO ielem=1,nelem
121 !
122  i1=ikle(ielem,1)
123  i2=ikle(ielem,2)
124  i3=ikle(ielem,3)
125  i4=ikle(ielem,4)
126 !
127 !-----------------------------------------------------------------------
128 !
129  x2=x(i2)-x(i1)
130  y2=y(i2)-y(i1)
131  z2=z(i2)-z(i1)
132  x3=x(i3)-x(i1)
133  y3=y(i3)-y(i1)
134  z3=z(i3)-z(i1)
135  x4=x(i4)-x(i1)
136  y4=y(i4)-y(i1)
137  z4=z(i4)-z(i1)
138 !
139 !-----------------------------------------------------------------------
140 !
141  u1 = u(ikle(ielem,1))
142  u2 = u(ikle(ielem,2))
143  u3 = u(ikle(ielem,3))
144  u4 = u(ikle(ielem,4))
145  v1 = v(ikle(ielem,1))
146  v2 = v(ikle(ielem,2))
147  v3 = v(ikle(ielem,3))
148  v4 = v(ikle(ielem,4))
149  w1 = w(ikle(ielem,1))
150  w2 = w(ikle(ielem,2))
151  w3 = w(ikle(ielem,3))
152  w4 = w(ikle(ielem,4))
153 !
154  xjac = x2*(y3*z4-y4*z3)+y2*(x4*z3-x3*z4)+z2*(x3*y4-x4*y3)
155 !
156  coef = sur120/max(xjac,1.d-8)
157 !
158  a1 = -y3*z4+y4*z3+y2*z4-z2*y4-y2*z3+z2*y3
159  a2 = -x3*z4+x4*z3+x2*z4-z2*x4-x2*z3+z2*x3
160  a3 = -x3*y4+x4*y3+x2*y4-y2*x4-x2*y3+y2*x3
161 !
162  aux1 = u3*u2+u1*u3+u1*u2+u4**2+u1*u4+u4*u2+u3**2+u2**2+u1**2+u3*u4
163 !
164  aux2 = v1*u4+2*v1*u1+v1*u2+v3*u1+v3*u2+v1*u3+v2*u1+v2*u3
165  & +v4*u1+2*v2*u2+2*v3*u3+v4*u2+v2*u4+v4*u3+v3*u4+2*v4*u4
166 !
167  aux4 = v1**2+v1*v3+v3**2+v2**2+v3*v2+v1*v2+v4**2+v1*v4+v3*v4+v2*v4
168  aux5 = v3*w2+v3*w1+v1*w2+2*v1*w1+v2*w3+v2*w1+v1*w3+2*v3*w3+2*v2*w2
169  &+v4*w1+v4*w3+v4*w2+2*v4*w4+v3*w4+v2*w4+v1*w4
170  aux6 = w2**2+w3**2+w2*w3+w1**2+w1*w3+w4**2+w1*w2+w1*w4+w4*w2+w4*w3
171  aux7 = 2*u2*w2+2*u3*w3+u3*w2+u3*w1+u1*w2+2*u1*w1+u2*w3+u2*w1+
172  &u1*w3+u4*w2+u4*w1+2*u4*w4+u1*w4+u2*w4+u4*w3+u3*w4
173 !
174  t(ielem,1)= (a1**2*aux1-a1*a2*aux2+a1*a3*aux7+a2**2*aux4-a2*
175  &a3*aux5+a3**2*aux6)*2*coef
176 !
177  xm(ielem,1) = (2*(y3*z4-y4*z3)*a1*aux1-(y3*z4-y4*z3)*a2*aux2+
178  &(y3*z4-y4*z3)*a3*aux7-(x3*z4-x4*z3)*a1*aux2+
179  &2*(x3*z4-x4*z3)*a2*aux4-(x3*z4-x4*z3)*a3*aux5-
180  &(-x3*y4+x4*y3)*a1*aux7+(-x3*y4+x4*y3)*a2*aux5-
181  &2*(-x3*y4+x4*y3)*a3*aux6)*coef
182 !
183  xm(ielem,2) = (2*(-y2*z4+z2*y4)*a1*aux1-(-y2*z4+z2*y4)*a2*aux2+
184  &(-y2*z4+z2*y4)*a3*aux7-(-x2*z4+z2*x4)*a1*aux2+
185  &2*(-x2*z4+z2*x4)*a2*aux4-(-x2*z4+z2*x4)*a3*aux5-
186  &(x2*y4-y2*x4)*a1*aux7+(x2*y4-y2*x4)*a2*aux5-
187  &2*(x2*y4-y2*x4)*a3*aux6)*coef
188 !
189  xm(ielem,3) = (-2*(-y2*z3+z2*y3)*a1*aux1+(-y2*z3+z2*y3)*a2*aux2-
190  &(-y2*z3+z2*y3)*a3*aux7+(-x2*z3+z2*x3)*a1*aux2-
191  &2*(-x2*z3+z2*x3)*a2*aux4+(-x2*z3+z2*x3)*a3*aux5+
192  &(x2*y3-y2*x3)*a1*aux7-(x2*y3-y2*x3)*a2*aux5+
193  &2*(x2*y3-y2*x3)*a3*aux6)*coef
194 !
195  t(ielem,2) =((y3*z4-y4*z3)**2*aux1-
196  &(y3*z4-y4*z3)*(x3*z4-x4*z3)*aux2
197  &-(y3*z4-y4*z3)*(-x3*y4+x4*y3)*aux7+
198  &(x3*z4-x4*z3)**2*aux4+
199  &(x3*z4-x4*z3)*(-x3*y4+x4*y3)*aux5+
200  &(-x3*y4+x4*y3)**2*aux6)*2*coef
201 !
202  xm(ielem,4) = (2*(-y2*z4+z2*y4)*(y3*z4-y4*z3)*aux1-
203  &(-y2*z4+z2*y4)*(x3*z4-x4*z3)*aux2-
204  &(-y2*z4+z2*y4)*(-x3*y4+x4*y3)*aux7
205  &-(y3*z4-y4*z3)*(-x2*z4+z2*x4)*aux2+
206  &2*(-x2*z4+z2*x4)*(x3*z4-x4*z3)*aux4+
207  &(-x2*z4+z2*x4)*(-x3*y4+x4*y3)*aux5-
208  &(y3*z4-y4*z3)*(x2*y4-y2*x4)*aux7+
209  &(x3*z4-x4*z3)*(x2*y4-y2*x4)*aux5+
210  &2*(x2*y4-y2*x4)*(-x3*y4+x4*y3)*aux6)*coef
211 !
212  xm(ielem,5) = (-2*(-y2*z3+z2*y3)*(y3*z4-y4*z3)*aux1+
213  &(-y2*z3+z2*y3)*(x3*z4-x4*z3)*aux2+
214  &(-y2*z3+z2*y3)*(-x3*y4+x4*y3)*aux7+
215  &(y3*z4-y4*z3)*(-x2*z3+z2*x3)*aux2-
216  &2*(-x2*z3+z2*x3)*(x3*z4-x4*z3)*aux4-
217  &(-x2*z3+z2*x3)*(-x3*y4+x4*y3)*aux5+
218  &(y3*z4-y4*z3)*(x2*y3-y2*x3)*aux7-
219  &(x3*z4-x4*z3)*(x2*y3-y2*x3)*aux5-
220  &2*(x2*y3-y2*x3)*(-x3*y4+x4*y3)*aux6)*coef
221 !
222  t(ielem,3) = ((-y2*z4+z2*y4)**2*aux1-
223  &(-y2*z4+z2*y4)*(-x2*z4+z2*x4)*aux2-
224  &(-y2*z4+z2*y4)*(x2*y4-y2*x4)*aux7+
225  &(-x2*z4+z2*x4)**2*aux4+(-x2*z4+z2*x4)*(x2*y4-y2*x4)*aux5+
226  &(x2*y4-y2*x4)**2*aux6)*2*coef
227 !
228  xm(ielem,6) = (-2*(-y2*z3+z2*y3)*(-y2*z4+z2*y4)*aux1+
229  &(-y2*z3+z2*y3)*(-x2*z4+z2*x4)*aux2+
230  &(-y2*z3+z2*y3)*(x2*y4-y2*x4)*aux7+
231  &(-y2*z4+z2*y4)*(-x2*z3+z2*x3)*aux2-
232  &2*(-x2*z3+z2*x3)*(-x2*z4+z2*x4)*aux4-
233  &(-x2*z3+z2*x3)*(x2*y4-y2*x4)*aux5+
234  &(-y2*z4+z2*y4)*(x2*y3-y2*x3)*aux7-
235  &(-x2*z4+z2*x4)*(x2*y3-y2*x3)*aux5-
236  &2*(x2*y3-y2*x3)*(x2*y4-y2*x4)*aux6)*coef
237 !
238  t(ielem,4) = ((-y2*z3+z2*y3)**2*aux1-
239  &(-y2*z3+z2*y3)*(-x2*z3+z2*x3)*aux2-
240  &(-y2*z3+z2*y3)*(x2*y3-y2*x3)*aux7
241  &+(-x2*z3+z2*x3)**2*aux4+(-x2*z3+z2*x3)*(x2*y3-y2*x3)*aux5+
242  &(x2*y3-y2*x3)**2*aux6)*2*coef
243 !
244 !
245 !-----------------------------------------------------------------------
246 !
247  ENDDO ! IELEM
248 !
249 !---------------------------------------------------------------
250 !
251  ELSE IF((su%ELM.EQ.30.AND.sv%ELM.EQ.30.AND.sw%ELM.EQ.30).OR.
252  & (su%ELM.EQ.50.AND.sv%ELM.EQ.50.AND.sw%ELM.EQ.50) ) THEN
253 !
254 ! P0 DISCRETISATION OF DIFFUSION COEFFICIENTS
255 !
256 ! LOOP ON THE TETRAHEDRONS
257 !
258  DO ielem=1,nelem
259 !
260  i1=ikle(ielem,1)
261  i2=ikle(ielem,2)
262  i3=ikle(ielem,3)
263  i4=ikle(ielem,4)
264 !
265 !-----------------------------------------------------------------------
266 !
267  x2=x(i2)-x(i1)
268  y2=y(i2)-y(i1)
269  z2=z(i2)-z(i1)
270  x3=x(i3)-x(i1)
271  y3=y(i3)-y(i1)
272  z3=z(i3)-z(i1)
273  x4=x(i4)-x(i1)
274  y4=y(i4)-y(i1)
275  z4=z(i4)-z(i1)
276 !
277 !-----------------------------------------------------------------------
278 !
279  u1 = u(ielem)
280  u2 = u1
281  u3 = u1
282  u4 = u1
283  v1 = v(ielem)
284  v2 = v1
285  v3 = v1
286  v4 = v1
287  w1 = w(ielem)
288  w2 = w1
289  w3 = w1
290  w4 = w1
291 !
292  xjac = x2*(y3*z4-y4*z3)+y2*(x4*z3-x3*z4)+z2*(x3*y4-x4*y3)
293 !
294  coef = sur120/max(xjac,1.d-8)
295 !
296  a1 = -y3*z4+y4*z3+y2*z4-z2*y4-y2*z3+z2*y3
297  a2 = -x3*z4+x4*z3+x2*z4-z2*x4-x2*z3+z2*x3
298  a3 = -x3*y4+x4*y3+x2*y4-y2*x4-x2*y3+y2*x3
299 !
300 ! TODO:SIMPLIFICATIONS TO BE DONE HERE !!!!!!!!!!!!!!
301 !
302  aux1 = u3*u2+u1*u3+u1*u2+u4**2+u1*u4+u4*u2+u3**2+u2**2+u1**2+u3*u4
303 !
304  aux2 = v1*u4+2*v1*u1+v1*u2+v3*u1+v3*u2+v1*u3+v2*u1+v2*u3
305  & +v4*u1+2*v2*u2+2*v3*u3+v4*u2+v2*u4+v4*u3+v3*u4+2*v4*u4
306 !
307  aux4 = v1**2+v1*v3+v3**2+v2**2+v3*v2+v1*v2+v4**2+v1*v4+v3*v4+v2*v4
308  aux5 = v3*w2+v3*w1+v1*w2+2*v1*w1+v2*w3+v2*w1+v1*w3+2*v3*w3+2*v2*w2
309  &+v4*w1+v4*w3+v4*w2+2*v4*w4+v3*w4+v2*w4+v1*w4
310  aux6 = w2**2+w3**2+w2*w3+w1**2+w1*w3+w4**2+w1*w2+w1*w4+w4*w2+w4*w3
311  aux7 = 2*u2*w2+2*u3*w3+u3*w2+u3*w1+u1*w2+2*u1*w1+u2*w3+u2*w1+
312  &u1*w3+u4*w2+u4*w1+2*u4*w4+u1*w4+u2*w4+u4*w3+u3*w4
313 !
314  t(ielem,1)= (a1**2*aux1-a1*a2*aux2+a1*a3*aux7+a2**2*aux4-a2*
315  &a3*aux5+a3**2*aux6)*2*coef
316 !
317  xm(ielem,1) = (2*(y3*z4-y4*z3)*a1*aux1-
318  &(y3*z4-y4*z3)*a2*aux2+
319  &(y3*z4-y4*z3)*a3*aux7-
320  &(x3*z4-x4*z3)*a1*aux2+
321  &2*(x3*z4-x4*z3)*a2*aux4-
322  &(x3*z4-x4*z3)*a3*aux5-
323  &(-x3*y4+x4*y3)*a1*aux7+
324  &(-x3*y4+x4*y3)*a2*aux5-
325  &2*(-x3*y4+x4*y3)*a3*aux6)*coef
326 !
327  xm(ielem,2) = (2*(-y2*z4+z2*y4)*
328  &a1*aux1-(-y2*z4+z2*y4)*a2*aux2+(-y2*z4+z2*y4)*a3*aux7-
329  &(-x2*z4+z2*x4)*a1*aux2+2*(-x2*z4+z2*x4)*a2*aux4-
330  &(-x2*z4+z2*x4)*a3*aux5-(x2*y4-y2*x4)*a1*aux7+
331  &(x2*y4-y2*x4)*a2*aux5-2*(x2*y4-y2*x4)*a3*aux6)*coef
332 !
333  xm(ielem,3) = (-2*(-y2*z3+z2*y3)*a1*aux1+(-y2*z3+z2*y3)*a2*aux2-
334  &(-y2*z3+z2*y3)*a3*aux7+(-x2*z3+z2*x3)*a1*aux2-
335  &2*(-x2*z3+z2*x3)*a2*aux4+(-x2*z3+z2*x3)*a3*aux5+
336  &(x2*y3-y2*x3)*a1*aux7-(x2*y3-y2*x3)*a2*aux5+
337  &2*(x2*y3-y2*x3)*a3*aux6)*coef
338 !
339  t(ielem,2) =((y3*z4-y4*z3)**2*aux1-
340  &(y3*z4-y4*z3)*(x3*z4-x4*z3)*aux2
341  &-(y3*z4-y4*z3)*(-x3*y4+x4*y3)*aux7+
342  &(x3*z4-x4*z3)**2*aux4+
343  &(x3*z4-x4*z3)*(-x3*y4+x4*y3)*aux5+
344  &(-x3*y4+x4*y3)**2*aux6)*2*coef
345 !
346  xm(ielem,4) = (2*(-y2*z4+z2*y4)*(y3*z4-y4*z3)*aux1-
347  &(-y2*z4+z2*y4)*(x3*z4-x4*z3)*aux2-
348  &(-y2*z4+z2*y4)*(-x3*y4+x4*y3)*aux7
349  &-(y3*z4-y4*z3)*(-x2*z4+z2*x4)*aux2+
350  &2*(-x2*z4+z2*x4)*(x3*z4-x4*z3)*aux4+
351  &(-x2*z4+z2*x4)*(-x3*y4+x4*y3)*aux5-
352  &(y3*z4-y4*z3)*(x2*y4-y2*x4)*aux7+
353  &(x3*z4-x4*z3)*(x2*y4-y2*x4)*aux5+
354  &2*(x2*y4-y2*x4)*(-x3*y4+x4*y3)*aux6)*coef
355 !
356  xm(ielem,5) = (-2*(-y2*z3+z2*y3)*(y3*z4-y4*z3)*aux1+
357  &(-y2*z3+z2*y3)*(x3*z4-x4*z3)*aux2+
358  &(-y2*z3+z2*y3)*(-x3*y4+x4*y3)*aux7+
359  &(y3*z4-y4*z3)*(-x2*z3+z2*x3)*aux2-
360  &2*(-x2*z3+z2*x3)*(x3*z4-x4*z3)*aux4-
361  &(-x2*z3+z2*x3)*(-x3*y4+x4*y3)*aux5+
362  &(y3*z4-y4*z3)*(x2*y3-y2*x3)*aux7-
363  &(x3*z4-x4*z3)*(x2*y3-y2*x3)*aux5-
364  &2*(x2*y3-y2*x3)*(-x3*y4+x4*y3)*aux6)*coef
365 !
366  t(ielem,3) = ((-y2*z4+z2*y4)**2*aux1-
367  &(-y2*z4+z2*y4)*(-x2*z4+z2*x4)*aux2-
368  &(-y2*z4+z2*y4)*(x2*y4-y2*x4)*aux7+
369  &(-x2*z4+z2*x4)**2*aux4+(-x2*z4+z2*x4)*(x2*y4-y2*x4)*aux5+
370  &(x2*y4-y2*x4)**2*aux6)*2*coef
371 !
372  xm(ielem,6) = (-2*(-y2*z3+z2*y3)*(-y2*z4+z2*y4)*aux1+
373  &(-y2*z3+z2*y3)*(-x2*z4+z2*x4)*aux2+
374  &(-y2*z3+z2*y3)*(x2*y4-y2*x4)*aux7+
375  &(-y2*z4+z2*y4)*(-x2*z3+z2*x3)*aux2-
376  &2*(-x2*z3+z2*x3)*(-x2*z4+z2*x4)*aux4-
377  &(-x2*z3+z2*x3)*(x2*y4-y2*x4)*aux5+
378  &(-y2*z4+z2*y4)*(x2*y3-y2*x3)*aux7-
379  &(-x2*z4+z2*x4)*(x2*y3-y2*x3)*aux5-
380  &2*(x2*y3-y2*x3)*(x2*y4-y2*x4)*aux6)*coef
381 !
382  t(ielem,4) = ((-y2*z3+z2*y3)**2*aux1-
383  &(-y2*z3+z2*y3)*(-x2*z3+z2*x3)*aux2-
384  &(-y2*z3+z2*y3)*(x2*y3-y2*x3)*aux7
385  &+(-x2*z3+z2*x3)**2*aux4+(-x2*z3+z2*x3)*(x2*y3-y2*x3)*aux5+
386  &(x2*y3-y2*x3)**2*aux6)*2*coef
387 !
388 !-----------------------------------------------------------------------
389 !
390  ENDDO ! IELEM
391 !
392 !---------------------------------------------------------------
393  ELSE
394 !
395  WRITE(lu,1001) su%ELM,sv%ELM,sw%ELM
396 1001 FORMAT(1x,'MT04TT (BIEF) : WRONG TYPE OF U,V OR W: ',
397  & i6,1x,i6,1x,i6)
398  CALL plante(1)
399  stop
400 !
401  ENDIF
402 !
403 ! ELSEIF(FORMUL(1:7).EQ.'MAUGUG1') THEN
404 !
405 !
406 !
407  ELSE
408  WRITE(lu,4001) formul
409 4001 FORMAT(1x,'MT04TT (BIEF) : UNEXPECTED FORMULA: ',a16)
410  CALL plante(1)
411  stop
412  ENDIF
413 !
414 !-----------------------------------------------------------------------
415 !
416  RETURN
417  END
subroutine mt04tt(T, XM, XMUL, SU, SV, SW, U, V, W, X, Y, Z, IKLE, NELEM, NELMAX, FORMUL)
Definition: mt04tt.f:7
Definition: bief.f:3