The TELEMAC-MASCARET system  trunk
mt14pp.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt14pp
3 ! *****************
4 !
5  &( t,xm,ppq,lego,xmul,sw,w,h,
6  & surfac,ikle,nelem,nelmax)
7 !
8 !***********************************************************************
9 ! BIEF V6P3 21/08/2010
10 !***********************************************************************
11 !
12 !brief BUILDS COEFFICIENTS LAMBDA(I,J) FOR N-TYPE MURD SCHEME.
13 !+
14 !+ THE ELEMENT IS THE P1 PRISM.
15 !
16 !warning THE JACOBIAN MUST BE POSITIVE
17 !
18 !reference J-M HERVOUET THESIS OR BOOK: MURD SCHEME IN 3 DIMENSIONS
19 !+ CHAPTER 6.
20 !
21 !history J-M JANIN (LNH)
22 !+ 28/11/1994
23 !+
24 !+
25 !
26 !history JMH
27 !+ 16/08/1999
28 !+
29 !+ EPS ADDED FOR THE DIVISION BY 0 TESTS
30 !
31 !history JMH
32 !+ 21/10/2004
33 !+
34 !+ MASS-LUMPING FOR COMPATIBILITY WITH NEW VERSION OF TRIDW2,
35 !
36 !history JMH
37 !+ 04/08/2008
38 !+
39 !+ INVERTED DIMENSIONS OF XM (SEE ALSO MURD3D)
40 !
41 !history J-M HERVOUET / A. DECOENE (LNHE)
42 !+ 27/04/2010
43 !+ V6P0
44 !+
45 !
46 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
47 !+ 13/07/2010
48 !+ V6P0
49 !+ Translation of French comments within the FORTRAN sources into
50 !+ English comments
51 !
52 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
53 !+ 21/08/2010
54 !+ V6P0
55 !+ Creation of DOXYGEN tags for automated documentation and
56 !+ cross-referencing of the FORTRAN sources
57 !
58 !history J-M HERVOUET (EDF R&D, LNHE)
59 !+ 11/01/2013
60 !+ V6P3
61 !+ Arguments X,Y,Z removed.
62 !
63 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
64 !| H |-->| FUNCTION USED IN THE FORMULA
65 !| IKLE |-->| CONNECTIVITY TABLE
66 !| LEGO |-->| LOGICAL. IF YES RESULT ASSEMBLED.
67 !| NELEM |-->| NUMBER OF ELEMENTS
68 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
69 !| PPQ |-->| STORAGE IN XM OF OFF-DIAGONAL ELEMENTS
70 !| SW |-->| BIEF_OBJ STRUCTURE OF W
71 !| SURFAC |-->| AREA OF 2D ELEMENTS
72 !| T |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
73 !| W |-->| FUNCTION USED IN THE FORMULA
74 !| XM |<->| OFF-DIAGONAL TERMS
75 !| XMUL |-->| COEFFICIENT FOR MULTIPLICATION
76 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
77 !
78  USE bief, ex_mt14pp => mt14pp
79 !
81  IMPLICIT NONE
82 !
83 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
84 !
85  INTEGER, INTENT(IN) :: NELEM,NELMAX
86  INTEGER, INTENT(IN) :: IKLE(nelmax,6),PPQ(6,6)
87 !
88  DOUBLE PRECISION, INTENT(INOUT) :: T(nelmax,6),XM(30,nelmax)
89 !
90  DOUBLE PRECISION, INTENT(IN) :: XMUL
91  DOUBLE PRECISION, INTENT(IN) :: W(*)
92  DOUBLE PRECISION, INTENT(IN) :: H(nelmax,6)
93 !
94  LOGICAL, INTENT(IN) :: LEGO
95 !
96 ! STRUCTURES DE U,V,W
97 !
98  TYPE(bief_obj), INTENT(IN) :: SW
99 !
100  DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
101 !
102 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
103 !
104 ! DECLARATIONS SPECIFIC TO THIS SUBROUTINE
105 !
106  INTEGER IELEM,IELMW
107  DOUBLE PRECISION XSUR3
108  DOUBLE PRECISION LI0J0,LI0K0,LJ0I0,LJ0K0,LK0I0,LK0J0
109  DOUBLE PRECISION LI3J3,LI3K3,LJ3I3,LJ3K3,LK3I3,LK3J3
110  DOUBLE PRECISION LI3J0,LI3K0,LJ3I0,LJ3K0,LK3I0,LK3J0
111  DOUBLE PRECISION LI3I0,LJ3J0,LK3K0,LI0I3
112  DOUBLE PRECISION ALFA,ALFAJ,ALFAK,ALFA1,ALFA2,BETA,SOM0,SOM3
113  DOUBLE PRECISION ALFAI0,ALFAJ0,ALFAK0,ALFAI3,ALFAJ3,ALFAK3
114  DOUBLE PRECISION ALFAII,ALFAIJ,ALFAIK,ALFAJI,ALFAJJ,ALFAJK
115  DOUBLE PRECISION ALFAKI,ALFAKJ,ALFAKK
116  DOUBLE PRECISION, PARAMETER :: EPS = 1.d-10
117 !
118  INTEGER :: IPLUS1(6),IPLUS3(6),INDIC(0:7)
119  INTEGER IXM,I0,J0,K0,I3,J3,K3
120 !
121 !-----------------------------------------------------------------------
122 !
123 ! DONNEES
124 !
125  parameter( iplus1 = (/ 2 , 3 , 1 , 5 , 6 , 4 /) )
126  parameter( iplus3 = (/ 4 , 5 , 6 , 1 , 2 , 3 /) )
127  parameter( indic = (/ 4 , 4 , 5 , 3 , 6 , 2 , 1 , 1 /) )
128 !
129 !=======================================================================
130 !
131  ielmw = sw%ELM
132 ! * *
133 ! COMPUTES THE COEFFICIENTS A(I) (INTEGRAL OF H U.GRAD(PSI), ONLY
134 ! CONSIDERING THE HORIZONTAL PART OF U). SEE JMH THESIS
135 !
136 ! NOTE: THIS VC04PP CALL IS ALREADY DONE BY FLUX3D IN TELEMAC-3D
137 ! THROUGH A CALL TO VECTOR. THE EQUIVALENT OF T HAS BEEN
138 ! SAVED IN WHAT IS HERE H
139 !
140 ! FORMUL=' HOR'
141 ! CALL VC04PP(XMUL,SU,SV,SW,U,V,W,SF,SG,SH,F,G,H,X,Y,Z,
142 ! * IKLE(1,1),IKLE(1,2),IKLE(1,3),IKLE(1,4),IKLE(1,5),IKLE(1,6),
143 ! * NELEM,NELMAX,T(1,1),T(1,2),T(1,3),T(1,4),T(1,5),T(1,6),
144 ! * SPECAD,FORMUL,NPLAN)
145 !
146 ! NOW T IS RETRIEVED BY STORAGE DONE INTO FLUX3D
147 ! FUNCTION H MUST CONTAIN THE EQUIVALENT OF T ABOVE, WHICH IS
148 ! THE NON ASSEMBLED VALUE OF FLUINT
149 !
150 ! DO IELEM=1,NELEM
151 ! T(IELEM,1)=H( IELEM)
152 ! T(IELEM,2)=H( NELMAX+IELEM)
153 ! T(IELEM,3)=H(2*NELMAX+IELEM)
154 ! T(IELEM,4)=H(3*NELMAX+IELEM)
155 ! T(IELEM,5)=H(4*NELMAX+IELEM)
156 ! T(IELEM,6)=H(5*NELMAX+IELEM)
157 ! ENDDO
158 !
159 ! FORMUL NORMALLY STARTS WITH VGRADP OR VGRADP2, OR VGRADP22 TO KNOW
160 ! SIGMAG AND SPECAD, USELESS HERE.
161  xsur3 = xmul/3.d0
162 !
163 !-----------------------------------------------------------------------
164 !
165  IF(ielmw.NE.41) THEN
166  WRITE(lu,*) 'MT14PP DISCRETISATION NOT HANDLED'
167  CALL plante(1)
168  stop
169  ENDIF
170 !
171 ! COMPUTES THE COEFFICIENTS B(I) BY ELEMENTS
172 ! AND COMPUTES LAMBDA(I,J)
173 !
174  DO ielem = 1 , nelem
175 !
176  ixm = 0
177  IF (w(ikle(ielem,1)).GT.0.d0) ixm = 1
178  IF (w(ikle(ielem,2)).GT.0.d0) ixm = ixm + 2
179  IF (w(ikle(ielem,3)).GT.0.d0) ixm = ixm + 4
180 !
181 ! NOTE JMH:
182 ! IT SEEMS THAT INDIC IS THE POINT THAT RECEIVES THE VERTICAL
183 ! VELOCITY WHICH IS THE ONLY ONE OF ITS SIGN
184 !
185 ! ALL W < 0 INDIC = 4 (RANDOM, WHY NOT 1 ?)
186 ! ONLY W1 > 0 INDIC = 4
187 ! ONLY W2 > 0 INCIC = 5
188 ! ONLY W3 > 0 INDIC = 6
189 ! ALL W > 0 INDIC = 1 (RANDOM, WHY NOT 4 ?)
190 ! ONLY W1 < 0 INDIC = 1
191 ! ONLY W2 < 0 INDIC = 2
192 ! ONLY W3 < 0 INDIC = 3
193 !
194  i0 = indic(ixm)
195  j0 = iplus1(i0)
196  k0 = iplus1(j0)
197  i3 = iplus3(i0)
198  j3 = iplus1(i3)
199  k3 = iplus1(j3)
200 !
201  alfa = xsur3 * surfac(ielem)
202 !
203  li3i0 = alfa * abs(w(ikle(ielem,min(i0,i3))))
204  li0i3 = 0.d0
205  IF (ixm.GE.1.AND.ixm.LE.6) THEN
206  li0i3 = li3i0
207  li3i0 = 0.d0
208  ENDIF
209  lj3j0 = alfa * abs(w(ikle(ielem,min(j0,j3))))
210  xm(ppq(j0,j3),ielem) = 0.d0
211  lk3k0 = alfa * abs(w(ikle(ielem,min(k0,k3))))
212  xm(ppq(k0,k3),ielem) = 0.d0
213 !
214  alfai0 = -h(ielem,i0)
215  alfaj0 = -h(ielem,j0)
216  alfak0 = -h(ielem,k0)
217  alfai3 = -h(ielem,i3)
218  alfaj3 = -h(ielem,j3)
219  alfak3 = -h(ielem,k3)
220 !
221  lj0i0 = max(0.d0,min(alfai0,-alfaj0))
222  lk0i0 = max(0.d0,min(alfai0,-alfak0))
223  li3j3 = max(0.d0,min(alfaj3,-alfai3))
224  li3k3 = max(0.d0,min(alfak3,-alfai3))
225 !
226  alfaj = min(lj3j0,max(li3j3,lj0i0))
227  alfak = min(lk3k0,max(li3k3,lk0i0))
228  alfa = min(li0i3,alfaj+alfak)
229  IF (alfa.GT.eps) THEN
230 !
231  beta = alfa / (alfaj+alfak)
232  alfaj = alfaj * beta
233  alfak = alfak * beta
234 !
235  li0i3 = li0i3-alfa
236  lj3j0 = lj3j0-alfaj
237  lk3k0 = lk3k0-alfak
238 !
239  alfai3 = alfai3+alfa
240  alfai0 = alfai0-alfa
241  alfaj0 = alfaj0+alfaj
242  alfaj3 = alfaj3-alfaj
243  alfak0 = alfak0+alfak
244  alfak3 = alfak3-alfak
245 !
246  lj0i0 = max(0.d0,min(alfai0,-alfaj0))
247  lk0i0 = max(0.d0,min(alfai0,-alfak0))
248  li3j3 = max(0.d0,min(alfaj3,-alfai3))
249  li3k3 = max(0.d0,min(alfak3,-alfai3))
250 !
251  ENDIF
252 !
253  li0j0 = max(0.d0,min(alfaj0,-alfai0))
254  lk0j0 = max(0.d0,min(alfaj0,-alfak0))
255  li0k0 = max(0.d0,min(alfak0,-alfai0))
256  lj0k0 = max(0.d0,min(alfak0,-alfaj0))
257  lj3i3 = max(0.d0,min(alfai3,-alfaj3))
258  lj3k3 = max(0.d0,min(alfak3,-alfaj3))
259  lk3i3 = max(0.d0,min(alfai3,-alfak3))
260  lk3j3 = max(0.d0,min(alfaj3,-alfak3))
261 !
262  xm(ppq(j0,i3),ielem) = 0.d0
263  xm(ppq(k0,i3),ielem) = 0.d0
264  xm(ppq(i0,j3),ielem) = 0.d0
265  xm(ppq(k0,j3),ielem) = 0.d0
266  xm(ppq(i0,k3),ielem) = 0.d0
267  xm(ppq(j0,k3),ielem) = 0.d0
268 !
269  alfai0 = 0.d0
270  alfaj0 = 0.d0
271  alfak0 = 0.d0
272  alfai3 = 0.d0
273  alfaj3 = 0.d0
274  alfak3 = 0.d0
275 !
276  som0 = lj0i0+lk0i0
277  som3 = li3j3+li3k3
278 !
279  alfa1 = min(li0i3,som0)
280  IF (alfa1.GT.eps) THEN
281 !
282  beta = 1.d0 / som0
283  alfaj0 = lj0i0 * beta
284  alfak0 = lk0i0 * beta
285  alfa = max(0.d0,alfa1-som3)
286 !
287  lj0i0 = lj0i0 - alfaj0*alfa1
288  lk0i0 = lk0i0 - alfak0*alfa1
289  xm(ppq(j0,i3),ielem) = alfaj0*alfa
290  xm(ppq(k0,i3),ielem) = alfak0*alfa
291 !
292  ENDIF
293 !
294  alfa2 = min(li0i3,som3)
295  IF (alfa2.GT.eps) THEN
296 !
297  beta = 1.d0 / som3
298  alfaj3 = li3j3 * beta
299  alfak3 = li3k3 * beta
300  alfa = max(0.d0,alfa2-som0)
301 !
302  li3j3 = li3j3 - alfaj3*alfa2
303  li3k3 = li3k3 - alfak3*alfa2
304  xm(ppq(i0,j3),ielem) = alfaj3*alfa
305  xm(ppq(i0,k3),ielem) = alfak3*alfa
306 !
307  ENDIF
308 !
309  alfa = min(alfa1,alfa2)
310  IF (alfa.GT.0.d0) THEN
311 !
312  alfaj3 = alfaj3 * alfa
313  alfak3 = alfak3 * alfa
314  alfajj = alfaj3 * alfaj0
315  alfakk = alfak3 * alfak0
316  alfajk = alfaj3 * alfak0
317  alfakj = alfak3 * alfaj0
318  alfa = min(alfajk,alfakj)
319 !
320  xm(ppq(j0,j3),ielem) = alfajj + alfa
321  xm(ppq(k0,k3),ielem) = alfakk + alfa
322  xm(ppq(k0,j3),ielem) = alfajk - alfa
323  xm(ppq(j0,k3),ielem) = alfakj - alfa
324 !
325  ENDIF
326 !
327  xm(ppq(i0,i3),ielem) = max(0.d0,li0i3-max(som3,som0))
328 !
329  lj3i0 = 0.d0
330  lk3i0 = 0.d0
331  li3j0 = 0.d0
332  lk3j0 = 0.d0
333  li3k0 = 0.d0
334  lj3k0 = 0.d0
335 !
336  som0 = li0j0+li0k0
337  som3 = lj3i3+lk3i3
338 !
339  alfa1 = min(li3i0,som0)
340  IF (alfa1.GT.eps) THEN
341 !
342  beta = 1.d0 / som0
343  alfaj0 = li0j0 * beta
344  alfak0 = li0k0 * beta
345  alfa = max(0.d0,alfa1-som3)
346 !
347  li0j0 = li0j0 - alfaj0*alfa1
348  li0k0 = li0k0 - alfak0*alfa1
349  li3j0 = li3j0 + alfaj0*alfa
350  li3k0 = li3k0 + alfak0*alfa
351 !
352  ENDIF
353 !
354  xm(ppq(i0,j0),ielem) = li0j0
355  xm(ppq(i0,k0),ielem) = li0k0
356 !
357  alfa2 = min(li3i0,som3)
358  IF (alfa2.GT.eps) THEN
359 !
360  beta = 1.d0 / som3
361  alfaj3 = lj3i3 * beta
362  alfak3 = lk3i3 * beta
363  alfa = max(0.d0,alfa2-som0)
364 !
365  lj3i3 = lj3i3 - alfaj3*alfa2
366  lk3i3 = lk3i3 - alfak3*alfa2
367  lj3i0 = lj3i0 + alfaj3*alfa
368  lk3i0 = lk3i0 + alfak3*alfa
369 !
370  ENDIF
371 !
372  xm(ppq(j3,i3),ielem) = lj3i3
373  xm(ppq(k3,i3),ielem) = lk3i3
374 !
375  alfa = min(alfa1,alfa2)
376  IF (alfa.GT.0.d0) THEN
377 !
378  alfaj0 = alfaj0 * alfa
379  alfak0 = alfak0 * alfa
380  alfajj = alfaj0 * alfaj3
381  alfakk = alfak0 * alfak3
382  alfajk = alfaj0 * alfak3
383  alfakj = alfak0 * alfaj3
384  alfa = min(alfajk,alfakj)
385 !
386  lj3j0 = lj3j0 + alfajj + alfa
387  lk3k0 = lk3k0 + alfakk + alfa
388  lk3j0 = lk3j0 + alfajk - alfa
389  lj3k0 = lj3k0 + alfakj - alfa
390 !
391  ENDIF
392 !
393  li3i0 = max(0.d0,li3i0-max(som0,som3))
394 !
395  som0 = lj0i0+lj0k0
396  som3 = li3j3+lk3j3
397 !
398  alfa1 = min(lj3j0,som0)
399  IF (alfa1.GT.eps) THEN
400 !
401  beta = 1.d0 / som0
402  alfai0 = lj0i0 * beta
403  alfak0 = lj0k0 * beta
404  alfa = max(0.d0,alfa1-som3)
405 !
406  lj0i0 = lj0i0 - alfai0*alfa1
407  lj0k0 = lj0k0 - alfak0*alfa1
408  lj3i0 = lj3i0 + alfai0*alfa
409  lj3k0 = lj3k0 + alfak0*alfa
410 !
411  ENDIF
412 !
413  xm(ppq(j0,i0),ielem) = lj0i0
414  xm(ppq(j0,k0),ielem) = lj0k0
415 !
416  alfa2 = min(lj3j0,som3)
417  IF (alfa2.GT.eps) THEN
418 !
419  beta = 1.d0 / som3
420  alfai3 = li3j3 * beta
421  alfak3 = lk3j3 * beta
422  alfa = max(0.d0,alfa2-som0)
423 !
424  li3j3 = li3j3 - alfai3*alfa2
425  lk3j3 = lk3j3 - alfak3*alfa2
426  li3j0 = li3j0 + alfai3*alfa
427  lk3j0 = lk3j0 + alfak3*alfa
428 !
429  ENDIF
430 !
431  xm(ppq(i3,j3),ielem) = li3j3
432  xm(ppq(k3,j3),ielem) = lk3j3
433 !
434  alfa = min(alfa1,alfa2)
435  IF (alfa.GT.0.d0) THEN
436 !
437  alfai0 = alfai0 * alfa
438  alfak0 = alfak0 * alfa
439  alfaii = alfai0 * alfai3
440  alfakk = alfak0 * alfak3
441  alfaik = alfai0 * alfak3
442  alfaki = alfak0 * alfai3
443  alfa = min(alfaik,alfaki)
444 !
445  li3i0 = li3i0 + alfaii + alfa
446  lk3k0 = lk3k0 + alfakk + alfa
447  lk3i0 = lk3i0 + alfaik - alfa
448  li3k0 = li3k0 + alfaki - alfa
449 !
450  ENDIF
451 !
452  lj3j0 = max(0.d0,lj3j0-max(som0,som3))
453 !
454  som0 = lk0i0+lk0j0
455  som3 = li3k3+lj3k3
456 !
457  alfa1 = min(lk3k0,som0)
458  IF (alfa1.GT.eps) THEN
459 !
460  beta = 1.d0 / som0
461  alfai0 = lk0i0 * beta
462  alfaj0 = lk0j0 * beta
463  alfa = max(0.d0,alfa1-som3)
464 !
465  lk0i0 = lk0i0 - alfai0*alfa1
466  lk0j0 = lk0j0 - alfaj0*alfa1
467  lk3i0 = lk3i0 + alfai0*alfa
468  lk3j0 = lk3j0 + alfaj0*alfa
469 !
470  ENDIF
471 !
472  xm(ppq(k0,i0),ielem) = lk0i0
473  xm(ppq(k0,j0),ielem) = lk0j0
474 !
475  alfa2 = min(lk3k0,som3)
476  IF (alfa2.GT.eps) THEN
477 !
478  beta = 1.d0 / som3
479  alfai3 = li3k3 * beta
480  alfaj3 = lj3k3 * beta
481  alfa = max(0.d0,alfa2-som0)
482 !
483  li3k3 = li3k3 - alfai3*alfa2
484  lj3k3 = lj3k3 - alfaj3*alfa2
485  li3k0 = li3k0 + alfai3*alfa
486  lj3k0 = lj3k0 + alfaj3*alfa
487 !
488  ENDIF
489 !
490  xm(ppq(i3,k3),ielem) = li3k3
491  xm(ppq(j3,k3),ielem) = lj3k3
492 !
493  alfa = min(alfa1,alfa2)
494  IF (alfa.GT.0.d0) THEN
495 !
496  alfai0 = alfai0 * alfa
497  alfaj0 = alfaj0 * alfa
498  alfaii = alfai0 * alfai3
499  alfajj = alfaj0 * alfaj3
500  alfaij = alfai0 * alfaj3
501  alfaji = alfaj0 * alfai3
502  alfa = min(alfaij,alfaji)
503 !
504  li3i0 = li3i0 + alfaii + alfa
505  lj3j0 = lj3j0 + alfajj + alfa
506  lj3i0 = lj3i0 + alfaij - alfa
507  li3j0 = li3j0 + alfaji - alfa
508 !
509  ENDIF
510 !
511  lk3k0 = max(0.d0,lk3k0-max(som0,som3))
512 !
513  xm(ppq(i3,i0),ielem) = li3i0
514  xm(ppq(j3,i0),ielem) = lj3i0
515  xm(ppq(k3,i0),ielem) = lk3i0
516  xm(ppq(i3,j0),ielem) = li3j0
517  xm(ppq(j3,j0),ielem) = lj3j0
518  xm(ppq(k3,j0),ielem) = lk3j0
519  xm(ppq(i3,k0),ielem) = li3k0
520  xm(ppq(j3,k0),ielem) = lj3k0
521  xm(ppq(k3,k0),ielem) = lk3k0
522 !
523  ENDDO ! IELEM
524 !
525 !-----------------------------------------------------------------------
526 !
527 ! COMPUTES THE SUM OF EACH ROW (WITH A - SIGN)
528 ! THE DIAGONAL TERMS ARE 0
529 !
530  IF(lego) THEN
531 !
532  DO ielem = 1,nelem
533 !
534  t(ielem,1) = -xm(01,ielem)-xm(02,ielem)
535  & -xm(03,ielem)-xm(04,ielem)-xm(05,ielem)
536  t(ielem,2) = -xm(16,ielem)-xm(06,ielem)
537  & -xm(07,ielem)-xm(08,ielem)-xm(09,ielem)
538  t(ielem,3) = -xm(17,ielem)-xm(21,ielem)
539  & -xm(10,ielem)-xm(11,ielem)-xm(12,ielem)
540  t(ielem,4) = -xm(18,ielem)-xm(22,ielem)
541  & -xm(25,ielem)-xm(13,ielem)-xm(14,ielem)
542  t(ielem,5) = -xm(19,ielem)-xm(23,ielem)
543  & -xm(26,ielem)-xm(28,ielem)-xm(15,ielem)
544  t(ielem,6) = -xm(20,ielem)-xm(24,ielem)
545  & -xm(27,ielem)-xm(29,ielem)-xm(30,ielem)
546 !
547  ENDDO
548 !
549  ENDIF
550 !
551 !-----------------------------------------------------------------------
552 !
553  RETURN
554  END
subroutine mt14pp(T, XM, PPQ, LEGO, XMUL, SW, W, H, SURFAC, IKLE, NELEM, NELMAX)
Definition: mt14pp.f:8
Definition: bief.f:3