The TELEMAC-MASCARET system  trunk
mt02pt.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt02pt
3 ! *****************
4 !
5  &(t,xm,xmul,sf,sg,sh,f,g,h,x,y,z,ikle,nelem,nelmax,inchyd)
6 !
7 !***********************************************************************
8 ! BIEF V6P2 21/08/2010
9 !***********************************************************************
10 !
11 !brief COMPUTES THE DIFFUSION MATRIX.
12 !+
13 !+ THE FUNCTION DIFFUSION COEFFICIENT IS HERE A P1
14 !+ DIAGONAL TENSOR.
15 !+
16 !+ CASE WHERE THE PRISM IS DIVIDED INTO 3 TETRAHEDRONS.
17 !+ OPTIMISED COEFFICIENT COMPUTATION.
18 !code
19 !+ STORAGE CONVENTION FOR EXTRA-DIAGONAL TERMS:
20 !+
21 !+ XM(IELEM, 1) ----> M(1,2) = M(2,1)
22 !+ XM(IELEM, 2) ----> M(1,3) = M(3,1)
23 !+ XM(IELEM, 3) ----> M(1,4) = M(4,1)
24 !+ XM(IELEM, 4) ----> M(1,5) = M(5,1)
25 !+ XM(IELEM, 5) ----> M(1,6) = M(6,1)
26 !+ XM(IELEM, 6) ----> M(2,3) = M(3,2)
27 !+ XM(IELEM, 7) ----> M(2,4) = M(4,2)
28 !+ XM(IELEM, 8) ----> M(2,5) = M(5,2)
29 !+ XM(IELEM, 9) ----> M(2,6) = M(6,2)
30 !+ XM(IELEM,10) ----> M(3,4) = M(4,3)
31 !+ XM(IELEM,11) ----> M(3,5) = M(5,3)
32 !+ XM(IELEM,12) ----> M(3,6) = M(6,3)
33 !+ XM(IELEM,13) ----> M(4,5) = M(5,4)
34 !+ XM(IELEM,14) ----> M(4,6) = M(6,4)
35 !+ XM(IELEM,15) ----> M(5,6) = M(6,5)
36 !
37 !history J-M HERVOUET (LNH) ; F LEPEINTRE (LNH)
38 !+ 28/11/94
39 !+ V5P3
40 !+
41 !
42 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
43 !+ 13/07/2010
44 !+ V6P0
45 !+ Translation of French comments within the FORTRAN sources into
46 !+ English comments
47 !
48 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
49 !+ 21/08/2010
50 !+ V6P0
51 !+ Creation of DOXYGEN tags for automated documentation and
52 !+ cross-referencing of the FORTRAN sources
53 !
54 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
55 !| F |-->| FUNCTION USED IN THE FORMULA
56 !| FORMUL |-->| FORMULA DESCRIBING THE RESULTING MATRIX
57 !| G |-->| FUNCTION USED IN THE FORMULA
58 !| H |-->| FUNCTION USED IN THE FORMULA
59 !| IKLE |-->| CONNECTIVITY TABLE.
60 !| INCHYD |-->| IF YES, TREATS HYDROSTATIC INCONSISTENCIES
61 !| NELEM |-->| NUMBER OF ELEMENTS
62 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
63 !| NPLAN |-->| NUMBER OF PLANES IN THE MESH OF PRISMS
64 !| SF |-->| STRUCTURE OF FUNCTIONS F
65 !| SG |-->| STRUCTURE OF FUNCTIONS G
66 !| SH |-->| STRUCTURE OF FUNCTIONS H
67 !| SURFAC |-->| AREA OF 2D ELEMENTS
68 !| T |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
69 !| X |-->| ABSCISSAE OF POINTS
70 !| Y |-->| ORDINATES OF POINTS
71 !| Z |-->| ELEVATIONS OF POINTS
72 !| XM |<->| OFF-DIAGONAL TERMS
73 !| XMUL |-->| COEFFICIENT FOR MULTIPLICATION
74 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
75 !
76  USE bief, ex_mt02pt => mt02pt
77  USE declarations_telemac, ONLY : tetra
78 !
79  USE declarations_special
80  IMPLICIT NONE
81 !
82 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
83 !
84  INTEGER, INTENT(IN) :: NELEM,NELMAX
85  INTEGER, INTENT(IN) :: IKLE(nelmax,6)
86  DOUBLE PRECISION, INTENT(INOUT) :: T(nelmax,6),XM(nelmax,15)
87  DOUBLE PRECISION, INTENT(IN) :: XMUL
88  DOUBLE PRECISION, INTENT(IN) :: F(*),G(*),H(*)
89  TYPE(bief_obj), INTENT(IN) :: SF,SG,SH
90  DOUBLE PRECISION, INTENT(IN) :: X(*),Y(*),Z(*)
91  LOGICAL, INTENT(IN) :: INCHYD
92 !
93 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
94 !
95 ! SPECIFIC DECLARATIONS
96 !
97  DOUBLE PRECISION X2,Y2,Z2,X3,Y3,Z3,X4,Y4,Z4
98  DOUBLE PRECISION DIAG1,DIAG2,DIAG3,DIAG4
99  DOUBLE PRECISION EXTR12,EXTR13,EXTR14,EXTR23,EXTR24,EXTR34
100  INTEGER IT1,IT2,IT3,IT4,I,I1,I2,I3,I4,I5,I6,NUM1,NUM2,NUM3,NUM4
101  INTEGER :: IGLOB(6),SENS(3),STO(6,6),IELEM
102 !
103  DOUBLE PRECISION SUR24,HTOT,VTOT,WTOT,COEF
104  DOUBLE PRECISION T1,T3,T5,T7,T9,T11,T13,T15,T17,T19,T21,T23
105  DOUBLE PRECISION T35,T49,T28,T42,T51,T54,SURF
106 !
107 ! PRISM SPLITTING (SEE EXPLANATIONS IN DECLARATIONS_TELEMAC)
108 !
109 ! INTEGER TETRA(2,2,2,3,4)
110 ! DATA TETRA / 0,1,1,1,1,1,1,0,0,4,4,4,4,4,4,0,0,6,4,5,5,4,6,0,
111 ! & 0,2,2,2,2,2,2,0,0,6,6,6,6,6,6,0,0,3,1,2,2,1,3,0,
112 ! & 0,3,3,3,3,3,3,0,0,5,5,5,5,5,5,0,0,2,3,4,1,6,5,0,
113 ! & 0,4,5,4,6,6,5,0,0,2,3,3,1,2,1,0,0,4,5,3,6,2,1,0 /
114 !
115 ! STORAGE CONVENTION FOR THE PRISM EXTRADIAGONAL TERMS
116 !
117  parameter( sto = reshape((/
118  & 00 , 01 , 02 , 03 , 04 , 05 ,
119  & 01 , 00 , 06 , 07 , 08 , 09 ,
120  & 02 , 06 , 00 , 10 , 11 , 12 ,
121  & 03 , 07 , 10 , 00 , 13 , 14 ,
122  & 04 , 08 , 11 , 13 , 00 , 15 ,
123  & 05 , 09 , 12 , 14 , 15 , 00 /), (/6,6/)) )
124 !
125 !***********************************************************************
126 !
127 ! DIFFERENT WAYS TO SPLIT PRISMS (TO ENSURE MATCHING OF TETRAHEDRONS)
128 !
129 ! TETRA(2,2,2,3,4)
130 !
131 ! FIRST 3 DIMENSIONS : TYPE OF FACE
132 ! 1 : CUT RECTANGLE BETWEEN LOW-LEFT AND HIGH-RIGHT
133 ! 2 : CUT RECTANGLE BETWEEN HIGH-LEFT AND LOW-RIGHT
134 !
135 ! 4TH DIMENSION : NUMBER OF TETRAHEDRON
136 ! 5TH DIMENSION : 4 POINTS OF THE TETRAHEDRON (IN LOCAL PRISM NUMBERING)
137 !
138 ! 1 1 2 CUTTING
139 !
140 ! TETRA(1,1,2,1,1)= 1
141 ! TETRA(1,1,2,1,2)= 2
142 ! TETRA(1,1,2,1,3)= 3
143 ! TETRA(1,1,2,1,4)= 6
144 !
145 ! TETRA(1,1,2,2,1)= 4
146 ! TETRA(1,1,2,2,2)= 6
147 ! TETRA(1,1,2,2,3)= 5
148 ! TETRA(1,1,2,2,4)= 1
149 !
150 ! TETRA(1,1,2,3,1)= 5
151 ! TETRA(1,1,2,3,2)= 2
152 ! TETRA(1,1,2,3,3)= 1
153 ! TETRA(1,1,2,3,4)= 6
154 !
155 ! 2 1 1 CUTTING
156 !
157 ! TETRA(2,1,1,1,1)= 1
158 ! TETRA(2,1,1,1,2)= 2
159 ! TETRA(2,1,1,1,3)= 3
160 ! TETRA(2,1,1,1,4)= 4
161 !
162 ! TETRA(2,1,1,2,1)= 4
163 ! TETRA(2,1,1,2,2)= 6
164 ! TETRA(2,1,1,2,3)= 5
165 ! TETRA(2,1,1,2,4)= 2
166 !
167 ! TETRA(2,1,1,3,1)= 6
168 ! TETRA(2,1,1,3,2)= 3
169 ! TETRA(2,1,1,3,3)= 2
170 ! TETRA(2,1,1,3,4)= 4
171 !
172 ! 1 2 1 CUTTING
173 !
174 ! TETRA(1,2,1,1,1)= 1
175 ! TETRA(1,2,1,1,2)= 2
176 ! TETRA(1,2,1,1,3)= 3
177 ! TETRA(1,2,1,1,4)= 5
178 !
179 ! TETRA(1,2,1,2,1)= 4
180 ! TETRA(1,2,1,2,2)= 6
181 ! TETRA(1,2,1,2,3)= 5
182 ! TETRA(1,2,1,2,4)= 3
183 !
184 ! TETRA(1,2,1,3,1)= 4
185 ! TETRA(1,2,1,3,2)= 1
186 ! TETRA(1,2,1,3,3)= 3
187 ! TETRA(1,2,1,3,4)= 5
188 !
189 ! 2 2 1 CUTTING
190 !
191 ! TETRA(2,2,1,1,1)= 1
192 ! TETRA(2,2,1,1,2)= 2
193 ! TETRA(2,2,1,1,3)= 3
194 ! TETRA(2,2,1,1,4)= 4
195 !
196 ! TETRA(2,2,1,2,1)= 4
197 ! TETRA(2,2,1,2,2)= 6
198 ! TETRA(2,2,1,2,3)= 5
199 ! TETRA(2,2,1,2,4)= 3
200 !
201 ! TETRA(2,2,1,3,1)= 5
202 ! TETRA(2,2,1,3,2)= 2
203 ! TETRA(2,2,1,3,3)= 4
204 ! TETRA(2,2,1,3,4)= 3
205 !
206 ! 1 2 2 CUTTING
207 !
208 ! TETRA(1,2,2,1,1)= 1
209 ! TETRA(1,2,2,1,2)= 2
210 ! TETRA(1,2,2,1,3)= 3
211 ! TETRA(1,2,2,1,4)= 5
212 !
213 ! TETRA(1,2,2,2,1)= 4
214 ! TETRA(1,2,2,2,2)= 6
215 ! TETRA(1,2,2,2,3)= 5
216 ! TETRA(1,2,2,2,4)= 1
217 !
218 ! TETRA(1,2,2,3,1)= 6
219 ! TETRA(1,2,2,3,2)= 3
220 ! TETRA(1,2,2,3,3)= 5
221 ! TETRA(1,2,2,3,4)= 1
222 !
223 ! 2 1 2 CUTTING
224 !
225 ! TETRA(2,1,2,1,1)= 1
226 ! TETRA(2,1,2,1,2)= 2
227 ! TETRA(2,1,2,1,3)= 3
228 ! TETRA(2,1,2,1,4)= 6
229 !
230 ! TETRA(2,1,2,2,1)= 4
231 ! TETRA(2,1,2,2,2)= 6
232 ! TETRA(2,1,2,2,3)= 5
233 ! TETRA(2,1,2,2,4)= 2
234 !
235 ! TETRA(2,1,2,3,1)= 4
236 ! TETRA(2,1,2,3,2)= 1
237 ! TETRA(2,1,2,3,3)= 6
238 ! TETRA(2,1,2,3,4)= 2
239 !
240 !-----------------------------------------------------------------------
241 !
242  sur24=1.d0/24.d0
243 !
244  IF(sf%ELM.EQ.41.AND.sg%ELM.EQ.41.AND.sh%ELM.EQ.41) THEN
245 !
246 !-----------------------------------------------------------------------
247 !
248 ! LINEAR DISCRETISATION OF DIFFUSION COEFFICIENTS
249 !
250 ! LOOP ON THE PRISMS
251 !
252  DO ielem=1,nelem
253 !
254  iglob(1)=ikle(ielem,1)
255  iglob(2)=ikle(ielem,2)
256  iglob(3)=ikle(ielem,3)
257  iglob(4)=ikle(ielem,4)
258  iglob(5)=ikle(ielem,5)
259  iglob(6)=ikle(ielem,6)
260 !
261  IF(iglob(1).GT.iglob(2)) THEN
262  sens(1)=1
263  ELSE
264  sens(1)=2
265  ENDIF
266  IF(iglob(2).GT.iglob(3)) THEN
267  sens(2)=1
268  ELSE
269  sens(2)=2
270  ENDIF
271  IF(iglob(3).GT.iglob(1)) THEN
272  sens(3)=1
273  ELSE
274  sens(3)=2
275  ENDIF
276 !
277 ! FOOTPRINT OF THE PRISM
278 !
279  surf=0.5d0*
280  & ((x(iglob(2))-x(iglob(1)))*(y(iglob(3))-y(iglob(1)))
281  & -(x(iglob(3))-x(iglob(1)))*(y(iglob(2))-y(iglob(1))))
282 !
283 ! INITIALISES TO 0.D0
284 !
285  t(ielem,1)=0.d0
286  t(ielem,2)=0.d0
287  t(ielem,3)=0.d0
288  t(ielem,4)=0.d0
289  t(ielem,5)=0.d0
290  t(ielem,6)=0.d0
291 !
292  xm(ielem, 1)= 0.d0
293  xm(ielem, 2)= 0.d0
294  xm(ielem, 3)= 0.d0
295  xm(ielem, 4)= 0.d0
296  xm(ielem, 5)= 0.d0
297 !
298  xm(ielem, 6)= 0.d0
299  xm(ielem, 7)= 0.d0
300  xm(ielem, 8)= 0.d0
301  xm(ielem, 9)= 0.d0
302 !
303  xm(ielem,10)= 0.d0
304  xm(ielem,11)= 0.d0
305  xm(ielem,12)= 0.d0
306 !
307  xm(ielem,13)= 0.d0
308  xm(ielem,14)= 0.d0
309  xm(ielem,15)= 0.d0
310 !
311 !-----------------------------------------------------------------------
312 ! LOOP OVER TI
313 !
314  DO i=1,3
315 !
316 ! TETRAHEDRON POINTS NUMBERS IN THE PRISM NUMBERING
317 !
318  num1=tetra(sens(1),sens(2),sens(3),i,1)
319  num2=tetra(sens(1),sens(2),sens(3),i,2)
320  num3=tetra(sens(1),sens(2),sens(3),i,3)
321  num4=tetra(sens(1),sens(2),sens(3),i,4)
322 !
323 ! GLOBAL NUMBERS OF THE TETRAHEDRON POINTS
324 !
325  it1=iglob(num1)
326  it2=iglob(num2)
327  it3=iglob(num3)
328  it4=iglob(num4)
329 !
330 ! VISCOSITY ALONG X Y AND Z
331 !
332  htot=f(it1)+f(it2)+f(it3)+f(it4)
333  vtot=g(it1)+g(it2)+g(it3)+g(it4)
334  wtot=h(it1)+h(it2)+h(it3)+h(it4)
335 !
336  x2=x(it2)-x(it1)
337  y2=y(it2)-y(it1)
338  z2=z(it2)-z(it1)
339  x3=x(it3)-x(it1)
340  y3=y(it3)-y(it1)
341  z3=z(it3)-z(it1)
342  x4=x(it4)-x(it1)
343  y4=y(it4)-y(it1)
344  z4=z(it4)-z(it1)
345 !
346 !-----------------------------------------------------------------------
347 ! COEF: THANKS MAPLE...
348 !-----------------------------------------------------------------------
349 !
350  t1 = x2*y3
351  t3 = x2*y4
352  t5 = x3*y2
353  t7 = x4*y2
354  t9 = x3*z2
355  t11 = x4*z2
356 ! T13 = 4 TIMES THE TETRAHEDRON VOLUME ?
357  t13 = t1*z4-t3*z3-t5*z4+t7*z3+t9*y4-t11*y3
358 !
359  t15 = -y2*z3+y3*z2
360  t17 = x2*z4-t11
361  t19 = -y3*z4+y4*z3
362  t21 = x2*z3-t9
363  t23 = -y2*z4+y4*z2
364  t35 = x3*z4-x4*z3
365  t49 = x3*y4-x4*y3
366 !
367 ! IF WIDTH MORE THAN 0.01 M
368 !
369  coef=xmul*sur24/max(t13,0.01d0*surf)
370 !
371  t28 = -t19+t23-t15
372  t42 = -t35+t17-t21
373  t51 = t3-t7
374  t54 = t49-t3+t7+t1-t5
375 !
376  diag1 =coef*( htot*t28**2 +vtot*t42**2 +wtot*t54**2)
377  diag2 =coef*( htot*t19**2 +vtot*t35**2 +wtot*t49**2)
378  diag3 =coef*( htot*t23**2 +vtot*t17**2 +wtot*t51**2)
379  extr12 =coef*( htot*t28*t19+vtot*t42*t35-wtot*t54*t49)
380  extr13 =coef*(-htot*t28*t23-vtot*t42*t17+wtot*t54*t51)
381  extr23 =coef*(-htot*t19*t23-vtot*t35*t17-wtot*t49*t51)
382 !
383 ! DEDUCED FROM PROPERTIES OF THE DIFFUSION MATRIX
384 !
385  extr14 = -(extr13+extr12+diag1)
386  extr24 = -(extr23+diag2+extr12)
387  extr34 = -(diag3+extr23+extr13)
388  diag4 = -(extr14+extr24+extr34)
389 !
390 ! ASSEMBLES ON PRISM
391 !
392  t(ielem,num1) = t(ielem,num1)+ diag1
393  t(ielem,num2) = t(ielem,num2)+ diag2
394  t(ielem,num3) = t(ielem,num3)+ diag3
395  t(ielem,num4) = t(ielem,num4)+ diag4
396 !
397  xm(ielem,sto(num1,num2))=xm(ielem,sto(num1,num2))+extr12
398  xm(ielem,sto(num1,num3))=xm(ielem,sto(num1,num3))+extr13
399  xm(ielem,sto(num1,num4))=xm(ielem,sto(num1,num4))+extr14
400  xm(ielem,sto(num2,num3))=xm(ielem,sto(num2,num3))+extr23
401  xm(ielem,sto(num2,num4))=xm(ielem,sto(num2,num4))+extr24
402  xm(ielem,sto(num3,num4))=xm(ielem,sto(num3,num4))+extr34
403 !
404  ENDDO ! I
405 !
406 !---------------------------------------------------------------
407 !
408  ENDDO ! IELEM
409 !
410  ELSEIF(sf%ELM.EQ.40.AND.sg%ELM.EQ.40.AND.sh%ELM.EQ.40) THEN
411 !
412 !
413 !-----------------------------------------------------------------------
414 !
415 ! P0 DISCRETISATION OF DIFFUSION COEFFICIENTS (CONSTANT ON A PRISM)
416 !
417 ! LOOP ON THE PRISMS
418 !
419  DO ielem=1,nelem
420 !
421  iglob(1)=ikle(ielem,1)
422  iglob(2)=ikle(ielem,2)
423  iglob(3)=ikle(ielem,3)
424  iglob(4)=ikle(ielem,4)
425  iglob(5)=ikle(ielem,5)
426  iglob(6)=ikle(ielem,6)
427 !
428  IF(iglob(1).GT.iglob(2)) THEN
429  sens(1)=1
430  ELSE
431  sens(1)=2
432  ENDIF
433  IF(iglob(2).GT.iglob(3)) THEN
434  sens(2)=1
435  ELSE
436  sens(2)=2
437  ENDIF
438  IF(iglob(3).GT.iglob(1)) THEN
439  sens(3)=1
440  ELSE
441  sens(3)=2
442  ENDIF
443 !
444 ! FOOTPRINT OF THE PRISM
445 !
446  surf=0.5d0*
447  & ((x(iglob(2))-x(iglob(1)))*(y(iglob(3))-y(iglob(1)))
448  & -(x(iglob(3))-x(iglob(1)))*(y(iglob(2))-y(iglob(1))))
449 !
450 ! INITIALISES TO 0.D0
451 !
452  t(ielem,1)=0.d0
453  t(ielem,2)=0.d0
454  t(ielem,3)=0.d0
455  t(ielem,4)=0.d0
456  t(ielem,5)=0.d0
457  t(ielem,6)=0.d0
458 !
459  xm(ielem, 1)= 0.d0
460  xm(ielem, 2)= 0.d0
461  xm(ielem, 3)= 0.d0
462  xm(ielem, 4)= 0.d0
463  xm(ielem, 5)= 0.d0
464 !
465  xm(ielem, 6)= 0.d0
466  xm(ielem, 7)= 0.d0
467  xm(ielem, 8)= 0.d0
468  xm(ielem, 9)= 0.d0
469 !
470  xm(ielem,10)= 0.d0
471  xm(ielem,11)= 0.d0
472  xm(ielem,12)= 0.d0
473 !
474  xm(ielem,13)= 0.d0
475  xm(ielem,14)= 0.d0
476  xm(ielem,15)= 0.d0
477 !
478 !-----------------------------------------------------------------------
479 ! LOOP OVER TI
480 !
481  DO i=1,3
482 !
483 ! TETRAHEDRON POINTS NUMBERS IN THE PRISM NUMBERING
484 !
485  num1=tetra(sens(1),sens(2),sens(3),i,1)
486  num2=tetra(sens(1),sens(2),sens(3),i,2)
487  num3=tetra(sens(1),sens(2),sens(3),i,3)
488  num4=tetra(sens(1),sens(2),sens(3),i,4)
489 !
490 ! GLOBAL NUMBERS OF THE TETRAHEDRON POINTS
491 !
492  it1=iglob(num1)
493  it2=iglob(num2)
494  it3=iglob(num3)
495  it4=iglob(num4)
496 !
497 ! VISCOSITY ALONG X Y AND Z
498 !
499  htot=4*f(ielem)
500  vtot=4*g(ielem)
501  wtot=4*h(ielem)
502 !
503  x2=x(it2)-x(it1)
504  y2=y(it2)-y(it1)
505  z2=z(it2)-z(it1)
506  x3=x(it3)-x(it1)
507  y3=y(it3)-y(it1)
508  z3=z(it3)-z(it1)
509  x4=x(it4)-x(it1)
510  y4=y(it4)-y(it1)
511  z4=z(it4)-z(it1)
512 !
513 !-----------------------------------------------------------------------
514 ! COEF: THANKS MAPLE...
515 !-----------------------------------------------------------------------
516 !
517  t1 = x2*y3
518  t3 = x2*y4
519  t5 = x3*y2
520  t7 = x4*y2
521  t9 = x3*z2
522  t11 = x4*z2
523 ! T13 = 4 TIMES THE TETRAHEDRON VOLUME ?
524  t13 = t1*z4-t3*z3-t5*z4+t7*z3+t9*y4-t11*y3
525 !
526  t15 = -y2*z3+y3*z2
527  t17 = x2*z4-t11
528  t19 = -y3*z4+y4*z3
529  t21 = x2*z3-t9
530  t23 = -y2*z4+y4*z2
531  t35 = x3*z4-x4*z3
532  t49 = x3*y4-x4*y3
533 !
534 ! IF WIDTH MORE THAN 0.01 M
535 !
536  coef=xmul*sur24/max(t13,0.01d0*surf)
537 !
538  t28 = -t19+t23-t15
539  t42 = -t35+t17-t21
540  t51 = t3-t7
541  t54 = t49-t3+t7+t1-t5
542 !
543  diag1 =coef*( htot*t28**2 +vtot*t42**2 +wtot*t54**2)
544  diag2 =coef*( htot*t19**2 +vtot*t35**2 +wtot*t49**2)
545  diag3 =coef*( htot*t23**2 +vtot*t17**2 +wtot*t51**2)
546  extr12 =coef*( htot*t28*t19+vtot*t42*t35-wtot*t54*t49)
547  extr13 =coef*(-htot*t28*t23-vtot*t42*t17+wtot*t54*t51)
548  extr23 =coef*(-htot*t19*t23-vtot*t35*t17-wtot*t49*t51)
549 !
550 ! DEDUCED FROM PROPERTIES OF THE DIFFUSION MATRIX
551 !
552  extr14 = -(extr13+extr12+diag1)
553  extr24 = -(extr23+diag2+extr12)
554  extr34 = -(diag3+extr23+extr13)
555  diag4 = -(extr14+extr24+extr34)
556 !
557 ! ASSEMBLES ON PRISM
558 !
559  t(ielem,num1) = t(ielem,num1)+ diag1
560  t(ielem,num2) = t(ielem,num2)+ diag2
561  t(ielem,num3) = t(ielem,num3)+ diag3
562  t(ielem,num4) = t(ielem,num4)+ diag4
563 !
564  xm(ielem,sto(num1,num2))=xm(ielem,sto(num1,num2))+extr12
565  xm(ielem,sto(num1,num3))=xm(ielem,sto(num1,num3))+extr13
566  xm(ielem,sto(num1,num4))=xm(ielem,sto(num1,num4))+extr14
567  xm(ielem,sto(num2,num3))=xm(ielem,sto(num2,num3))+extr23
568  xm(ielem,sto(num2,num4))=xm(ielem,sto(num2,num4))+extr24
569  xm(ielem,sto(num3,num4))=xm(ielem,sto(num3,num4))+extr34
570 !
571  ENDDO ! I
572 !
573 !---------------------------------------------------------------
574 !
575  ENDDO ! IELEM
576 !
577 !-----------------------------------------------------------------------
578 !
579 !
580  ELSE
581 !
582  WRITE(lu,1001) sf%ELM,sg%ELM,sh%ELM
583 1000 FORMAT(1x,'MT02PT (BIEF) : MAUVAIS TYPE DE F,G OU H : ',
584  & i6,1x,i6,1x,i6)
585 1001 FORMAT(1x,'MT02PT (BIEF) : WRONG TYPE OF F,G OR H: ',
586  & i6,1x,i6,1x,i6)
587  CALL plante(1)
588  stop
589 !
590  ENDIF
591 !
592 !-----------------------------------------------------------------------
593 !
594 ! TREATMENT OF HYDROSTATIC INCONSISTENCIES
595 !
596  IF(inchyd) THEN
597 !
598  DO ielem=1,nelem
599 !
600  i1=ikle(ielem,1)
601  i2=ikle(ielem,2)
602  i3=ikle(ielem,3)
603  i4=ikle(ielem,4)
604  i5=ikle(ielem,5)
605  i6=ikle(ielem,6)
606 !
607  IF(max(z(i1),z(i2),z(i3)).GT.min(z(i4),z(i5),z(i6))) THEN
608 !
609  t(ielem,1) =0.d0
610  t(ielem,2) =0.d0
611  t(ielem,3) =0.d0
612  t(ielem,4) =0.d0
613  t(ielem,5) =0.d0
614  t(ielem,6) =0.d0
615  xm(ielem, 1)=0.d0
616  xm(ielem, 2)=0.d0
617  xm(ielem, 3)=0.d0
618  xm(ielem, 4)=0.d0
619  xm(ielem, 5)=0.d0
620  xm(ielem, 6)=0.d0
621  xm(ielem, 7)=0.d0
622  xm(ielem, 8)=0.d0
623  xm(ielem, 9)=0.d0
624  xm(ielem,10)=0.d0
625  xm(ielem,11)=0.d0
626  xm(ielem,12)=0.d0
627  xm(ielem,13)=0.d0
628  xm(ielem,14)=0.d0
629  xm(ielem,15)=0.d0
630 !
631  ENDIF
632 !
633  ENDDO ! IELEM
634 !
635  ENDIF
636 !
637 !-----------------------------------------------------------------------
638 !
639  RETURN
640  END
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
declarations_telemac
Definition: declarations_telemac.f:3
declarations_telemac::tetra
integer, dimension(2, 2, 2, 3, 4) tetra
Definition: declarations_telemac.f:237
mt02pt
subroutine mt02pt(T, XM, XMUL, SF, SG, SH, F, G, H, X, Y, Z, IKLE, NELEM, NELMAX, INCHYD)
Definition: mt02pt.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