The TELEMAC-MASCARET system  trunk
mt02tt.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt02tt
3 ! *****************
4 !
5  &(t,xm,xmul,sf,sg,sh,f,g,h,x,y,z,ikle,nelem,nelmax,npoin2)
6 !
7 !***********************************************************************
8 ! BIEF V6P2 21/08/2010
9 !***********************************************************************
10 !
11 !brief COMPUTES THE DIFFUSION MATRIX FOR TETRAHEDRONS.
12 !+
13 !+ THE FUNCTION DIFFUSION COEFFICIENT IS HERE A P1
14 !+ DIAGONAL TENSOR.
15 !
16 !code
17 !+ STORAGE CONVENTION FOR EXTRA-DIAGONAL TERMS:
18 !+
19 !+ XM(IELEM, 1) ----> M(1,2) = M(2,1)
20 !+ XM(IELEM, 2) ----> M(1,3) = M(3,1)
21 !+ XM(IELEM, 3) ----> M(1,4) = M(4,1)
22 !+ XM(IELEM, 4) ----> M(2,3) = M(3,2)
23 !+ XM(IELEM, 5) ----> M(2,4) = M(4,2)
24 !+ XM(IELEM, 6) ----> M(3,4) = M(4,3)
25 !
26 !history J-M HERVOUET (LNH) ; F LEPEINTRE (LNH)
27 !+ 28/11/94
28 !+ V5P3
29 !+
30 !
31 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
32 !+ 13/07/2010
33 !+ V6P0
34 !+ Translation of French comments within the FORTRAN sources into
35 !+ English comments
36 !
37 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
38 !+ 21/08/2010
39 !+ V6P0
40 !+ Creation of DOXYGEN tags for automated documentation and
41 !+ cross-referencing of the FORTRAN sources
42 !
43 !history J-M HERVOUET (LNHE)
44 !+ 09/12/2011
45 !+ V6P2
46 !+ Element 51 now treated separately.
47 !
48 !history U.H. Merkel
49 !+ 18/07/2012
50 !+ V6P2
51 !+ Replaced EPSILON with CHOUIA due to nag compiler problems
52 !
53 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
54 !| F |-->| FUNCTION USED IN THE FORMULA
55 !| FORMUL |-->| FORMULA DESCRIBING THE RESULTING MATRIX
56 !| G |-->| FUNCTION USED IN THE FORMULA
57 !| H |-->| FUNCTION USED IN THE FORMULA
58 !| IKLE |-->| CONNECTIVITY TABLE.
59 !| INCHYD |-->| IF YES, TREATS HYDROSTATIC INCONSISTENCIES
60 !| NELEM |-->| NUMBER OF ELEMENTS
61 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
62 !| NPLAN |-->| NUMBER OF PLANES IN THE MESH OF PRISMS
63 !| NPOIN2 |-->| NUMBER OF POINTS IN UNDERLYING 2D MESH
64 !| | | (CASE OF ELEMENT 51, PRISMS CUT INTO TETRAHEDRA)
65 !| SF |-->| STRUCTURE OF FUNCTIONS F
66 !| SG |-->| STRUCTURE OF FUNCTIONS G
67 !| SH |-->| STRUCTURE OF FUNCTIONS H
68 !| SURFAC |-->| AREA OF 2D ELEMENTS
69 !| T |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
70 !| X |-->| ABSCISSAE OF POINTS
71 !| Y |-->| ORDINATES OF POINTS
72 !| Z |-->| ELEVATIONS OF POINTS
73 !| XM |<->| OFF-DIAGONAL TERMS
74 !| XMUL |-->| COEFFICIENT FOR MULTIPLICATION
75 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
76 !
77  USE bief, ex_mt02tt => mt02tt
78 !
79  USE declarations_special
80  IMPLICIT NONE
81 !
82 !
83 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
84 !
85  INTEGER, INTENT(IN) :: NELEM,NELMAX,NPOIN2
86  INTEGER, INTENT(IN) :: IKLE(nelmax,4)
87 !
88  DOUBLE PRECISION, INTENT(INOUT) :: T(nelmax,4),XM(nelmax,6)
89 !
90  DOUBLE PRECISION, INTENT(IN) :: XMUL
91  DOUBLE PRECISION, INTENT(IN) :: F(*),G(*),H(*)
92 !
93 ! STRUCTURES OF F, G, H
94 !
95  TYPE(bief_obj), INTENT(IN) :: SF,SG,SH
96 !
97  DOUBLE PRECISION, INTENT(IN) :: X(*),Y(*),Z(*)
98 !
99 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
100 !
101 ! SPECIFIC DECLARATIONS
102 !
103  DOUBLE PRECISION X2,Y2,Z2,X3,Y3,Z3,X4,Y4,Z4,VOL6
104  INTEGER I1,I2,I3,I4,IELEM,ISO,IL1,IL2,IL3,IL4,ILOW
105 !
106  DOUBLE PRECISION COEF,SUR24,DZ1,DZ2,DZ3,DZ4
107  DOUBLE PRECISION FTOT,GTOT,HTOT,VTOT,WTOT,FGTOT,GHTOT,FHTOT
108  DOUBLE PRECISION T1,T3,T5,T7,T9,T11,T13,T15,T17,T19,T21,T23
109  DOUBLE PRECISION T35,T49,T28,T42,T51,T54
110  DOUBLE PRECISION AUX,AUXX,AUXXX
111  DOUBLE PRECISION AUX1,AUX2,AUX3,AUX4,AUX5,AUX6,AUX7,AUX8,AUX9
112 !
113  DOUBLE PRECISION, PARAMETER :: CHOUIA = 1.d-4
114 !
115 !***********************************************************************
116 !
117  sur24=1.d0/24.d0
118  iso = sf%DIM2
119 !
120 ! SEE VISCLM OF TELEMAC-3D
121 ! FOR THE TREATMENT OF P0 VERTICAL VISCOSITY ON THE VERTICAL
122 !
123  IF(sh%DIMDISC.EQ.0) THEN
124 ! P1 VERTICAL VISCOSITY
125 ! NPOU0=SH%DIM1/NPLAN
126 ! ELSEIF(SH%DIMDISC.EQ.4111) THEN
127 ! P0 VERTICAL VISCOSITY ON THE VERTICAL (SEE II4,5,6)
128 ! NPOU0=0
129  ELSE
130  WRITE(lu,4001) sh%DIMDISC
131 4001 FORMAT(1x,'MT02TT (BIEF): DIMDISC OF H NOT IMPLEMENTED: ',i6)
132  CALL plante(1)
133  stop
134  ENDIF
135 !
136 ! TRUE TETRAHEDRA
137 !
138  IF(sf%ELM.EQ.31.AND.sg%ELM.EQ.31.AND.
139  & sh%ELM.EQ.31.AND.iso.EQ.1) THEN
140 !
141 !-----------------------------------------------------------------------
142 !
143 ! LINEAR DISCRETISATION OF DIFFUSION COEFFICIENTS
144 !
145 ! LOOP ON THE TETRAHEDRONS
146 !
147  DO ielem=1,nelem
148 !
149  i1=ikle(ielem,1)
150  i2=ikle(ielem,2)
151  i3=ikle(ielem,3)
152  i4=ikle(ielem,4)
153 !
154 !-----------------------------------------------------------------------
155 !
156 ! VISCOSITY ALONG X Y AND Z
157 !
158  htot=f(i1)+f(i2)+f(i3)+f(i4)
159  vtot=g(i1)+g(i2)+g(i3)+g(i4)
160  wtot=h(i1)+h(i2)+h(i3)+h(i4)
161 !
162  x2=x(i2)-x(i1)
163  y2=y(i2)-y(i1)
164  z2=z(i2)-z(i1)
165  x3=x(i3)-x(i1)
166  y3=y(i3)-y(i1)
167  z3=z(i3)-z(i1)
168  x4=x(i4)-x(i1)
169  y4=y(i4)-y(i1)
170  z4=z(i4)-z(i1)
171 !
172 !-----------------------------------------------------------------------
173 ! COEF: THANKS MAPLE...
174 !-----------------------------------------------------------------------
175 !
176  t1 = x2*y3
177  t3 = x2*y4
178  t5 = x3*y2
179  t7 = x4*y2
180  t9 = x3*z2
181  t11 = x4*z2
182 ! T13 = 4 TIMES THE TETRAHEDRON VOLUME ?
183  t13 = t1*z4-t3*z3-t5*z4+t7*z3+t9*y4-t11*y3
184 !
185  t15 = -y2*z3+y3*z2
186  t17 = x2*z4-t11
187  t19 = -y3*z4+y4*z3
188  t21 = x2*z3-t9
189  t23 = -y2*z4+y4*z2
190  t35 = x3*z4-x4*z3
191  t49 = x3*y4-x4*y3
192 !
193 ! BEWARE, PROBLEM WITH TIDAL FLATS HERE
194 !
195  coef=xmul*sur24/max(t13,1.d-10)
196 ! COEF=XMUL*SUR24/T13
197 !
198  t28 = -t19+t23-t15
199  t42 = -t35+t17-t21
200  t51 = t3-t7
201  t54 = t49-t3+t7+t1-t5
202 !
203  t(ielem,1) =coef*(htot*t28**2+vtot*t42**2+wtot*t54**2)
204  t(ielem,2) =coef*(htot*t19**2+vtot*t35**2+wtot*t49**2)
205  t(ielem,3) =coef*(htot*t23**2+vtot*t17**2+wtot*t51**2)
206  xm(ielem,1)=coef*(htot*t28*t19+vtot*t42*t35-wtot*t54*t49)
207  xm(ielem,2)=coef*(-htot*t28*t23-vtot*t42*t17+wtot*t54*t51)
208  xm(ielem,3)=-(xm(ielem,2)+xm(ielem,1)+t(ielem,1))
209  xm(ielem,4)=coef*(-htot*t19*t23-vtot*t35*t17-wtot*t49*t51)
210  xm(ielem,5)= -(xm(ielem,4)+t(ielem,2)+xm(ielem,1))
211  xm(ielem,6)= -(t(ielem,3)+xm(ielem,4)+xm(ielem,2))
212  t(ielem,4) = -(xm(ielem,3)+xm(ielem,5)+xm(ielem,6))
213 !
214 !-----------------------------------------------------------------------
215 !
216  ENDDO
217 !
218 ! TETRAHEDRA FORMING PRISMS
219 !
220  ELSEIF(sf%ELM.EQ.51.AND.sg%ELM.EQ.51.AND.sh%ELM.EQ.51.AND.
221  & iso.EQ.1) THEN
222 !
223 !-----------------------------------------------------------------------
224 !
225 ! LINEAR DISCRETISATION OF DIFFUSION COEFFICIENTS
226 !
227 ! LOOP ON THE TETRAHEDRONS
228 !
229  DO ielem=1,nelem
230 !
231  i1=ikle(ielem,1)
232  i2=ikle(ielem,2)
233  i3=ikle(ielem,3)
234  i4=ikle(ielem,4)
235 !
236 !-----------------------------------------------------------------------
237 !
238 ! VISCOSITY ALONG X, Y, AND Z
239 !
240 ! IL1: PLANE NUMBER OF POINT 1 MINUS 1 (0 FOR FIRST PLANE), ETC.
241  il1=(i1-1)/npoin2
242  il2=(i2-1)/npoin2
243  il3=(i3-1)/npoin2
244  il4=(i4-1)/npoin2
245 ! RANK OF LOWER PLANE MINUS 1
246  ilow=min(il1,il2,il3,il4)
247 ! NOW IL1=0 IF POINT AT THE BACK OF ORIGINAL PRISM
248 ! IL1=1 IF POINT AT THE TOP OF ORIGINAL PRISM ,ETC. FOR IL2...
249  il1=il1-ilow
250  il2=il2-ilow
251  il3=il3-ilow
252  il4=il4-ilow
253 ! DZ1, DZ2, DZ3, DZ4 WILL BE THE DELTA(Z) ABOVE OR BELOW EVERY POINT IN
254 ! THE ORIGINAL PRISM
255  IF(il1.EQ.0) THEN
256  dz1=z(i1+npoin2)-z(i1)
257  ELSE
258  dz1=z(i1)-z(i1-npoin2)
259  ENDIF
260  IF(il2.EQ.0) THEN
261  dz2=z(i2+npoin2)-z(i2)
262  ELSE
263  dz2=z(i2)-z(i2-npoin2)
264  ENDIF
265  IF(il3.EQ.0) THEN
266  dz3=z(i3+npoin2)-z(i3)
267  ELSE
268  dz3=z(i3)-z(i3-npoin2)
269  ENDIF
270  IF(il4.EQ.0) THEN
271  dz4=z(i4+npoin2)-z(i4)
272  ELSE
273  dz4=z(i4)-z(i4-npoin2)
274  ENDIF
275 !
276 ! SEE MT02PP FOR A SIMILAR LIMITATION
277 ! DIFFUSION CONNECTIONS CUT IF ONE POINT IS CRUSHED IN THE ELEMENT
278 ! THIS IS NECESSARY TO AVOID A MASS ERROR
279 !
280  IF(dz1.LT.chouia.OR.dz2.LT.chouia.OR.
281  & dz3.LT.chouia.OR.dz4.LT.chouia ) THEN
282  htot=0.d0
283  vtot=0.d0
284  wtot=0.d0
285  ELSE
286  htot=f(i1)+f(i2)+f(i3)+f(i4)
287  vtot=g(i1)+g(i2)+g(i3)+g(i4)
288  wtot=h(i1)+h(i2)+h(i3)+h(i4)
289  ENDIF
290 !
291  x2=x(i2)-x(i1)
292  y2=y(i2)-y(i1)
293  z2=z(i2)-z(i1)
294  x3=x(i3)-x(i1)
295  y3=y(i3)-y(i1)
296  z3=z(i3)-z(i1)
297  x4=x(i4)-x(i1)
298  y4=y(i4)-y(i1)
299  z4=z(i4)-z(i1)
300 !
301 !-----------------------------------------------------------------------
302 !
303  vol6 = x2*(y3*z4-y4*z3)+y2*(x4*z3-x3*z4)+z2*(x3*y4-x4*y3)
304 !
305  t15 = -y2*z3+y3*z2
306  t17 = x2*z4-x4*z2
307  t19 = -y3*z4+y4*z3
308  t21 = x2*z3-x3*z2
309  t23 = -y2*z4+y4*z2
310  t35 = x3*z4-x4*z3
311  t49 = x3*y4-x4*y3
312 !
313 ! BEWARE, PROBLEM WITH TIDAL FLATS HERE
314 !
315  coef=xmul*sur24/max(vol6,1.d-10)
316 !
317  t28 = -t19+t23-t15
318  t42 = -t35+t17-t21
319  t51 = x2*y4-x4*y2
320  t54 = t49-x2*y4+x4*y2+x2*y3-x3*y2
321 !
322  t(ielem,1) =coef*( htot*t28**2+vtot*t42**2+wtot*t54**2)
323  t(ielem,2) =coef*( htot*t19**2+vtot*t35**2+wtot*t49**2)
324  t(ielem,3) =coef*( htot*t23**2+vtot*t17**2+wtot*t51**2)
325  xm(ielem,1)=coef*( htot*t28*t19+vtot*t42*t35-wtot*t54*t49)
326  xm(ielem,2)=coef*(-htot*t28*t23-vtot*t42*t17+wtot*t54*t51)
327  xm(ielem,3)=-(xm(ielem,2)+xm(ielem,1)+t(ielem,1))
328  xm(ielem,4)=coef*(-htot*t19*t23-vtot*t35*t17-wtot*t49*t51)
329  xm(ielem,5)= -(xm(ielem,4)+t(ielem,2)+xm(ielem,1))
330  xm(ielem,6)= -(t(ielem,3)+xm(ielem,4)+xm(ielem,2))
331  t(ielem,4) = -(xm(ielem,3)+xm(ielem,5)+xm(ielem,6))
332 !
333 !-----------------------------------------------------------------------
334 !
335  ENDDO
336 !
337  ELSEIF(sf%ELM.EQ.30.AND.sg%ELM.EQ.30.AND.sh%ELM.EQ.30.AND.
338  & iso.EQ.1) THEN
339 !
340 !
341 !-----------------------------------------------------------------------
342 !
343 ! P0 DISCRETISATION OF DIFFUSION COEFFICIENTS (CONSTANT ON AN ELEMENT)
344 !
345 ! LOOP ON THE PRISMS
346 !
347  DO ielem=1,nelem
348 !
349  i1=ikle(ielem,1)
350  i2=ikle(ielem,2)
351  i3=ikle(ielem,3)
352  i4=ikle(ielem,4)
353 !
354 ! ONLY DIFFERENCE WITH LOOP 20
355  htot=4*f(ielem)
356  vtot=4*g(ielem)
357  wtot=4*h(ielem)
358 ! END OF ONLY DIFFERENCE WITH LOOP 20
359 !
360 !-----------------------------------------------------------------------
361 !
362  x2=x(i2)-x(i1)
363  y2=y(i2)-y(i1)
364  z2=z(i2)-z(i1)
365  x3=x(i3)-x(i1)
366  y3=y(i3)-y(i1)
367  z3=z(i3)-z(i1)
368  x4=x(i4)-x(i1)
369  y4=y(i4)-y(i1)
370  z4=z(i4)-z(i1)
371 !
372 !-----------------------------------------------------------------------
373 ! COEF: THANKS MAPLE...
374 !-----------------------------------------------------------------------
375 !
376  t1 = x2*y3
377  t3 = x2*y4
378  t5 = x3*y2
379  t7 = x4*y2
380  t9 = x3*z2
381  t11 = x4*z2
382 ! T13 = 4 TIMES THE TETRAHEDRON VOLUME ?
383  t13 = t1*z4-t3*z3-t5*z4+t7*z3+t9*y4-t11*y3
384 !
385  t15 = -y2*z3+y3*z2
386  t17 = x2*z4-t11
387  t19 = -y3*z4+y4*z3
388  t21 = x2*z3-t9
389  t23 = -y2*z4+y4*z2
390  t35 = x3*z4-x4*z3
391  t49 = x3*y4-x4*y3
392 !
393 ! IF WIDTH MORE THAN 0.01 M
394 !
395 ! BEWARE, PROBLEM WITH TIDAL FLAT HERE
396 !
397 ! COEF=XMUL*SUR24/MAX(T13,0.01D0*SURF)
398 ! COEF=XMUL*SUR24/MAX(T13,1.D-3)
399  coef=xmul*sur24/max(t13,1.d-10)
400 ! COEF=XMUL*SUR24/T13
401 !
402  t28 = -t19+t23-t15
403  t42 = -t35+t17-t21
404  t51 = t3-t7
405  t54 = t49-t3+t7+t1-t5
406 !
407  t(ielem,1) =coef*(htot*t28**2+vtot*t42**2+wtot*t54**2)
408  t(ielem,2) =coef*(htot*t19**2+vtot*t35**2+wtot*t49**2)
409  t(ielem,3) =coef*(htot*t23**2+vtot*t17**2+wtot*t51**2)
410  xm(ielem,1)=coef*(htot*t28*t19+vtot*t42*t35-wtot*t54*t49)
411  xm(ielem,2)=coef*(-htot*t28*t23-vtot*t42*t17+wtot*t54*t51)
412  xm(ielem,3)=-(xm(ielem,2)+xm(ielem,1)+t(ielem,1))
413  xm(ielem,4)=coef*(-htot*t19*t23-vtot*t35*t17-wtot*t49*t51)
414  xm(ielem,5)= -(xm(ielem,4)+t(ielem,2)+xm(ielem,1))
415  xm(ielem,6)= -(t(ielem,3)+xm(ielem,4)+xm(ielem,2))
416  t(ielem,4) = -(xm(ielem,3)+xm(ielem,5)+xm(ielem,6))
417 !
418 !---------------------------------------------------------------
419 !
420  ENDDO ! IELEM
421 !
422  ELSEIF(sf%ELM.EQ.30.AND.iso.EQ.6) THEN
423 !
424 ! P0 DISCRETISATION OF DIFFUSION COEFFICIENTS (CONSTANT ON AN ELEMENT)
425 ! NONISOTROPIC VISCOSITY ==> 6 COMPONENTS
426 ! THE FUNCTION DIFFUSION COEFFICIENT IS HERE A P0 NON DIAGONAL
427 ! SYMMETRICAL TENSOR
428 !
429 ! / DXX DXY DXZ \ / F FG FH \
430 ! D = | DYX DYY DYZ | = | X G GH |
431 ! \ DZX DZY DZZ / \ X XX H /
432 !-----------------------------------------------------------------------
433 !
434 ! LOOP ON THE TETRAHEDRONS
435 !
436  DO ielem=1,nelem
437 !
438  i1=ikle(ielem,1)
439  i2=ikle(ielem,2)
440  i3=ikle(ielem,3)
441  i4=ikle(ielem,4)
442 !
443  ftot = 4.d0*f(ielem )
444  gtot = 4.d0*f(ielem + nelem)
445  htot = 4.d0*f(ielem + 2*nelem)
446  ghtot= 4.d0*f(ielem + 3*nelem)
447  fhtot= 4.d0*f(ielem + 4*nelem)
448  fgtot= 4.d0*f(ielem + 5*nelem)
449 !
450 !-----------------------------------------------------------------------
451 !
452  x2=x(i2)-x(i1)
453  y2=y(i2)-y(i1)
454  z2=z(i2)-z(i1)
455  x3=x(i3)-x(i1)
456  y3=y(i3)-y(i1)
457  z3=z(i3)-z(i1)
458  x4=x(i4)-x(i1)
459  y4=y(i4)-y(i1)
460  z4=z(i4)-z(i1)
461 !
462 !-----------------------------------------------------------------------
463 ! COEF: THANKS MAPLE...
464 !-----------------------------------------------------------------------
465 !
466  t13 = x2*y3*z4-x2*y4*z3-y2*x3*z4+y2*x4*z3+z2*x3*y4-z2*x4*y3
467 !
468  aux = y3*z4-y4*z3-y2*z4+z2*y4+y2*z3-z2*y3
469  auxx = x3*z4-x4*z3-x2*z4+z2*x4+x2*z3-z2*x3
470  auxxx = x3*y4-x4*y3-x2*y4+y2*x4+x2*y3-y2*x3
471  aux1 = ftot*(y3*z4-y4*z3-y2*z4+z2*y4+y2*z3-z2*y3)
472  aux2 = fgtot*(-x3*z4+x4*z3+x2*z4-z2*x4-x2*z3+z2*x3)
473  aux3 = fhtot*(x3*y4-x4*y3-x2*y4+y2*x4+x2*y3-y2*x3)
474  aux4 = gtot*(-x3*z4+x4*z3+x2*z4-z2*x4-x2*z3+z2*x3)
475  aux5 = ghtot*(x3*y4-x4*y3-x2*y4+y2*x4+x2*y3-y2*x3)
476  aux6 = fgtot*(y3*z4-y4*z3-y2*z4+z2*y4+y2*z3-z2*y3)
477  aux7 = htot*(x3*y4-x4*y3-x2*y4+y2*x4+x2*y3-y2*x3)
478  aux8 = fhtot*(y3*z4-y4*z3-y2*z4+z2*y4+y2*z3-z2*y3)
479  aux9 = ghtot*(-x3*z4+x4*z3+x2*z4-z2*x4-x2*z3+z2*x3)
480 !
481  coef=xmul*sur24/max(t13,1.d-14)
482 !
483  t(ielem,1) =coef*(aux*(aux1+aux2+aux3)
484  & -auxx*(aux4+aux5+aux6)+auxxx*(aux7+aux8+aux9))
485  t(ielem,2) =coef*((y3*z4-y4*z3)*
486  & (ftot*(y3*z4-y4*z3)+fgtot*(-x3*z4+x4*z3)+fhtot*(x3*y4-x4*y3))
487  & -(x3*z4-x4*z3)*(gtot*(-x3*z4+x4*z3)+fgtot*(y3*z4-y4*z3)
488  & +ghtot*(x3*y4-x4*y3))+(x3*y4-x4*y3)*(htot*(x3*y4-x4*y3)
489  & +fhtot*(y3*z4-y4*z3)+ghtot*(-x3*z4+x4*z3)))
490  t(ielem,3) =coef*((y2*z4-z2*y4)*(ftot*(y2*z4-z2*y4)+
491  & fgtot*(-x2*z4+z2*x4)+fhtot*(x2*y4-y2*x4))-
492  & (x2*z4-z2*x4)*(gtot*(-x2*z4+z2*x4)+fgtot*(y2*z4-z2*y4)+
493  & ghtot*(x2*y4-y2*x4))+(x2*y4-y2*x4)*(htot*(x2*y4-y2*x4)+
494  & fhtot*(y2*z4-z2*y4)+ghtot*(-x2*z4+z2*x4)))
495  t(ielem,4) =coef*((y2*z3-z2*y3)*(ftot*(y2*z3-z2*y3)+
496  & fgtot*(-x2*z3+z2*x3)+fhtot*(x2*y3-y2*x3))-
497  & (x2*z3-z2*x3)*(gtot*(-x2*z3+z2*x3)+fgtot*(y2*z3-z2*y3)+
498  & ghtot*(x2*y3-y2*x3))+(x2*y3-y2*x3)*(htot*(x2*y3-y2*x3)+
499  & fhtot*(y2*z3-z2*y3)+ghtot*(-x2*z3+z2*x3)))
500  xm(ielem,1)=coef*((x3*z4-x4*z3)*(aux4+aux5+aux6)-
501  & (y3*z4-y4*z3)*(aux1+aux2+aux3)-(x3*y4-x4*y3)*(aux7+aux8+aux9))
502  xm(ielem,2)=coef*((y2*z4-z2*y4)*(aux1+aux2+aux3)-
503  & (x2*z4-z2*x4)*(aux4+aux5+aux6)+(x2*y4-y2*x4)*(aux7+aux8+aux9))
504  xm(ielem,3)=coef*((x2*z3-z2*x3)*(aux4+aux5+aux6)-
505  & (y2*z3-z2*y3)*(aux1+aux2+aux3)-(x2*y3-y2*x3)*(aux7+aux8+aux9))
506  xm(ielem,4)=coef*((x2*z4-z2*x4)*(gtot*(-x3*z4+x4*z3)+
507  & fgtot*(y3*z4-y4*z3)+ghtot*(x3*y4-x4*y3))-
508  & (y2*z4-z2*y4)*(ftot*(y3*z4-y4*z3)+fgtot*(-x3*z4+x4*z3)+
509  & fhtot*(x3*y4-x4*y3))-(x2*y4-y2*x4)*(htot*(x3*y4-x4*y3)+
510  & fhtot*(y3*z4-y4*z3)+ghtot*(-x3*z4+x4*z3)))
511  xm(ielem,5)=coef*((y2*z3-z2*y3)*(ftot*(y3*z4-y4*z3)+
512  & fgtot*(-x3*z4+x4*z3)+fhtot*(x3*y4-x4*y3))-
513  & (x2*z3-z2*x3)*(gtot*(-x3*z4+x4*z3)+fgtot*(y3*z4-y4*z3)
514  & +ghtot*(x3*y4-x4*y3))+(x2*y3-y2*x3)*(htot*(x3*y4-x4*y3)
515  & +fhtot*(y3*z4-y4*z3)+ghtot*(-x3*z4+x4*z3)))
516  xm(ielem,6)=coef*((x2*z3-z2*x3)*(gtot*(-x2*z4+z2*x4)+
517  & fgtot*(y2*z4-z2*y4)+ghtot*(x2*y4-y2*x4))-
518  & (y2*z3-z2*y3)*(ftot*(y2*z4-z2*y4)+fgtot*(-x2*z4+z2*x4)+
519  & fhtot*(x2*y4-y2*x4))-(x2*y3-y2*x3)*(htot*(x2*y4-y2*x4)+
520  & fhtot*(y2*z4-z2*y4)+ghtot*(-x2*z4+z2*x4)))
521  ENDDO ! IELEM
522 !
523  ELSE
524 !
525  WRITE(lu,1001) sf%ELM,sg%ELM,sh%ELM
526 1001 FORMAT(1x,'MT02TT (BIEF) : WRONG TYPE OF F,G OR H: ',
527  & i6,1x,i6,1x,i6)
528  CALL plante(1)
529  stop
530 !
531  ENDIF
532 !
533 !-----------------------------------------------------------------------
534 !
535  RETURN
536  END
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
mt02tt
subroutine mt02tt(T, XM, XMUL, SF, SG, SH, F, G, H, X, Y, Z, IKLE, NELEM, NELMAX, NPOIN2)
Definition: mt02tt.f:7
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:36:46 for The TELEMAC-MASCARET system by doxygen 1.8.13