The TELEMAC-MASCARET system  trunk
vc02pp_star.f
Go to the documentation of this file.
1 ! **********************
2  SUBROUTINE vc02pp_star
3 ! **********************
4 !
5  &( xmul,sf,sg,sh,su,f,g,h,u,x,y,z,surfac,
6  & ikle1,ikle2,ikle3,ikle4,ikle5,ikle6,nelem,nelmax,
7  & w1,w2,w3,w4,w5,w6,formul)
8 !
9 !***********************************************************************
10 ! BIEF V6P2 23/06/2011
11 !***********************************************************************
12 !
13 !brief COMPUTES THE PRODUCT OF THE DIFFUSION MATRIX BY FUNCTION U
14 ! CORRESPONDS TO MATRIX COMPUTED IN MT02PP_STAR
15 ! F, G AND H ARE THE DIFFUSION COEFFICIENTS ALONG X, Y AND Z
16 !
17 !history J-M HERVOUET (LNHE)
18 !+ 23/06/2011
19 !+ V6P1
20 !+ First version
21 !
22 !history U.H. MERKEL
23 !+ 18/07/2012
24 !+ V6P2
25 !+ Replaced EPSILON with CHOUIA due to nag compiler problems
26 !
27 !history J-M HERVOUET (EDF R&D LNHE)
28 !+ 07/01/2013
29 !+ V6P3
30 !+ X and Y are now given per element.
31 !
32 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
33 !| F |-->| FUNCTION USED IN THE VECTOR FORMULA
34 !| G |-->| FUNCTION USED IN THE VECTOR FORMULA
35 !| H |-->| FUNCTION USED IN THE VECTOR FORMULA
36 !| IKLE1 |-->| FIRST POINT OF PRISMS
37 !| IKLE2 |-->| SECOND POINT OF PRISMS
38 !| IKLE3 |-->| THIRD POINT OF PRISMS
39 !| IKLE4 |-->| FOURTH POINT OF PRISMS
40 !| IKLE5 |-->| FIFTH POINT OF PRISMS
41 !| IKLE6 |-->| SIXTH POINT OF PRISMS
42 !| NELEM |-->| NUMBER OF ELEMENTS
43 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
44 !| SF |-->| BIEF_OBJ STRUCTURE OF F
45 !| SG |-->| BIEF_OBJ STRUCTURE OF G
46 !| SH |-->| BIEF_OBJ STRUCTURE OF H
47 !| SURFAC |-->| AREA OF TRIANGLES
48 !| W1 |<--| RESULT IN NON ASSEMBLED FORM
49 !| W2 |<--| RESULT IN NON ASSEMBLED FORM
50 !| W3 |<--| RESULT IN NON ASSEMBLED FORM
51 !| W4 |<--| RESULT IN NON ASSEMBLED FORM
52 !| W5 |<--| RESULT IN NON ASSEMBLED FORM
53 !| W6 |<--| RESULT IN NON ASSEMBLED FORM
54 !| X |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
55 !| Y |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
56 !| XMUL |-->| MULTIPLICATION COEFFICIENT
57 !| Z |-->| ELEVATIONS OF POINTS, PER POINTS
58 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
59 !
60  USE bief ! , EX_VC04PP => VC04PP
61 !
62  USE declarations_special
63  IMPLICIT NONE
64 !
65 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
66 !
67  INTEGER, INTENT(IN) :: NELEM,NELMAX
68  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax),IKLE3(nelmax)
69  INTEGER, INTENT(IN) :: IKLE4(nelmax),IKLE5(nelmax),IKLE6(nelmax)
70 !
71  DOUBLE PRECISION, INTENT(IN) ::X(nelmax,6),Y(nelmax,6),Z(*)
72  DOUBLE PRECISION, INTENT(IN) ::SURFAC(*)
73  DOUBLE PRECISION, INTENT(INOUT)::W1(nelmax),W2(nelmax),W3(nelmax)
74  DOUBLE PRECISION, INTENT(INOUT)::W4(nelmax),W5(nelmax),W6(nelmax)
75  DOUBLE PRECISION, INTENT(IN) ::XMUL
76 !
77  CHARACTER(LEN=16), INTENT(IN) :: FORMUL
78 !
79 ! STRUCTURES DE U,V ET DONNEES REELLES
80 !
81  TYPE(bief_obj), INTENT(IN) :: SU,SF,SG,SH
82  DOUBLE PRECISION, INTENT(IN) :: U(*),F(*),G(*),H(*)
83 !
84 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
85 !
86  DOUBLE PRECISION NF1,NF2,NF3,NG1,NG2,NG3,NH1,NH2,NH3,XS03
87  DOUBLE PRECISION XS06,X2,X3,Y2,Y3,H1,H2,H3,SOMVX,SOMVY,AUX
88  DOUBLE PRECISION X2X3,Y2Y3,X2AUX,X3AUX,Y2AUX,Y3AUX,GX(3),GY(3)
89  DOUBLE PRECISION NPX,NPY,NUXMOY,NUYMOY,XS24,NPXL,NPXU,NPYL,NPYU
90  DOUBLE PRECISION XM01,XM02,XM03,XM04,XM05,XM06,XM07,XM08,XM09,XM10
91  DOUBLE PRECISION XM11,XM12,XM13,XM14,XM15,XM16,XM17,XM18,XM19,XM20
92  DOUBLE PRECISION XM21,XM22,XM23,XM24,XM25,XM26,XM27,XM28,XM29,XM30
93  DOUBLE PRECISION NPX2,NPY2
94  INTEGER I1,I2,I3,I4,I5,I6,IELEM
95  LOGICAL INCHYD
96 !
97  DOUBLE PRECISION :: CHOUIA = 1.d-4
98 !
99 !***********************************************************************
100 !
101  inchyd=.false.
102  IF(formul(7:7).EQ.'2') inchyd=.true.
103  xs03 = xmul / 3.d0
104  xs06 = xmul / 6.d0
105  xs24 = xmul /24.d0
106 !
107 !-----------------------------------------------------------------------
108 !
109  IF(sf%ELM.NE.41) THEN
110  WRITE(lu,1001) sf%ELM
111 1001 FORMAT(1x,'VC02PP_STAR (BIEF) : TYPE OF F NOT IMPLEMENTED: ',i6)
112  CALL plante(1)
113  stop
114  ENDIF
115  IF(sg%ELM.NE.41) THEN
116  WRITE(lu,2001) sg%ELM
117 2001 FORMAT(1x,'VC02PP_STAR (BIEF) : TYPE OF G NOT IMPLEMENTED: ',i6)
118  CALL plante(1)
119  stop
120  ENDIF
121  IF(sh%ELM.NE.41) THEN
122  WRITE(lu,3001) sh%ELM
123 3001 FORMAT(1x,'VC02PP_STAR (BIEF) : TYPE OF H NOT IMPLEMENTED: ',i6)
124  CALL plante(1)
125  stop
126  ENDIF
127  IF(su%ELM.NE.41) THEN
128  WRITE(lu,4001) su%ELM
129 4001 FORMAT(1x,'VC02PP_STAR (BIEF) : TYPE OF U NOT IMPLEMENTED: ',i6)
130  CALL plante(1)
131  stop
132  ENDIF
133 !
134 !-----------------------------------------------------------------------
135 !
136 ! TERMES HORIZONTAUX
137 !
138  IF(formul(14:16).EQ.'HOR') THEN
139 !
140 !-----------------------------------------------------------------------
141 !
142  DO ielem = 1 , nelem
143 !
144  i1 = ikle1(ielem)
145  i2 = ikle2(ielem)
146  i3 = ikle3(ielem)
147  i4 = ikle4(ielem)
148  i5 = ikle5(ielem)
149  i6 = ikle6(ielem)
150 !
151  h1 = z(i4) - z(i1)
152  h2 = z(i5) - z(i2)
153  h3 = z(i6) - z(i3)
154 !
155 ! X2 = X(I2)-X(I1)
156 ! X3 = X(I3)-X(I1)
157 ! Y2 = Y(I2)-Y(I1)
158 ! Y3 = Y(I3)-Y(I1)
159 !
160  x2=x(ielem,2)
161  x3=x(ielem,3)
162  y2=y(ielem,2)
163  y3=y(ielem,3)
164 !
165  IF((inchyd.AND.max(z(i1),z(i2),z(i3)).GT.
166  & min(z(i4),z(i5),z(i6))) .OR.
167  & h1.LT.chouia.OR.h2.LT.chouia.OR.h3.LT.chouia ) THEN
168  aux = 0.d0
169  nuxmoy=0.d0
170  nuymoy=0.d0
171  ELSE
172  aux = xs24 / surfac(ielem)
173 ! TAKING INTO ACCOUNT AN AVERAGE DIFFUSION ON THE PRISM
174  nuxmoy=(f(i1)+f(i2)+f(i3)+f(i4)+f(i5)+f(i6))/6.d0
175  nuymoy=(g(i1)+g(i2)+g(i3)+g(i4)+g(i5)+g(i6))/6.d0
176  ENDIF
177 !
178  x2x3 = x2*x3
179  y2y3 = y2*y3
180  x2aux = x2x3-x2**2
181  x3aux = x3**2-x2x3
182  y2aux = y2y3-y2**2
183  y3aux = y3**2-y2y3
184 !
185 ! 2D DIFFUSION, LOWER LEVEL (SEE MT02AA)
186 ! LIKE IN 2D BUT MULTIPLIED BY DELTA(Z)/2
187  somvx = ( f(i1)*h1+f(i2)*h2+f(i3)*h3 ) * aux
188  somvy = ( g(i1)*h1+g(i2)*h2+g(i3)*h3 ) * aux
189 !
190 ! OFF-DIAGONAL TERMS FOR POINTS OF LOWER LEVEL
191 ! WITH MONOTONICITY ENSURED
192 !
193 ! TERM 1
194 !
195  xm01 = - somvy * x3aux - somvx * y3aux
196  xm02 = somvy * x2aux + somvx * y2aux
197  xm06 = - somvy * x2x3 - somvx * y2y3
198  xm16 = xm01
199  xm17 = xm02
200  xm21 = xm06
201 !
202 ! 2D DIFFUSION, UPPER LEVEL
203  somvx = ( f(i4)*h1+f(i5)*h2+f(i6)*h3 ) * aux
204  somvy = ( g(i4)*h1+g(i5)*h2+g(i6)*h3 ) * aux
205 !
206 ! OFF-DIAGONAL TERMS FOR POINTS OF UPPER LEVEL
207 ! WITH MONOTONICITY ENSURED
208 !
209 ! TERM 1
210 !
211  xm13 = - somvy * x3aux - somvx * y3aux
212  xm14 = somvy * x2aux + somvx * y2aux
213  xm15 = - somvy * x2x3 - somvx * y2y3
214  xm28 = xm13
215  xm29 = xm14
216  xm30 = xm15
217 !
218 ! AVERAGE OF NORMAL VECTOR TO PLANES (NOT NORMED)
219 ! ONE CAN CHECK THAT WE GET -1 0 OR 0 -1 WITH Z=X OR Z=Y
220 ! NPX=-DZ/DX NPY=-DZ/DY
221 !
222  npxl=-0.5d0*(y2*(z(i1)-z(i3))+y3*(z(i2)-z(i1)))
223  npxu=-0.5d0*(y2*(z(i4)-z(i6))+y3*(z(i5)-z(i4)))
224  npyl=-0.5d0*(x2*(z(i3)-z(i1))+x3*(z(i1)-z(i2)))
225  npyu=-0.5d0*(x2*(z(i6)-z(i4))+x3*(z(i4)-z(i5)))
226  npx=0.5d0*(npxl+npxu)/surfac(ielem)
227  npy=0.5d0*(npyl+npyu)/surfac(ielem)
228 !
229  npx=npx*nuxmoy
230  npy=npy*nuymoy
231 !
232 ! 2D GRADIENT MATRIX GX(3,3) AND GY(3,3)
233 ! BUT GX(1,J)=GX(2,J)=GX(3,J)
234 ! AND GY(1,J)=GY(2,J)=GY(3,J), HENCE ONLY GX(3) AND GY(3)
235 !
236  gx(2) = y3*xs06
237  gx(3) = - y2*xs06
238  gx(1) = - gx(2) - gx(3)
239 !
240  gy(2) = - x3 * xs06
241  gy(3) = x2 * xs06
242  gy(1) = - gy(2) - gy(3)
243 !
244 ! TERM 3
245 !
246 ! TERM 1-2 (01)
247  xm01=xm01 +0.5d0*( - ( npx*gx(1)+npy*gy(1)) )
248 ! TERM 1-3 (02)
249  xm02=xm02 +0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
250 ! TERM 1-4 (03)
251  xm03= +0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
252 ! TERM 1-5 (04)
253  xm04= +0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
254 ! TERM 1-6 (05)
255  xm05= +0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
256 ! TERM 2-3 (06)
257  xm06=xm06 +0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
258 ! TERM 2-4 (07)
259  xm07= +0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
260 ! TERM 2-5 (08)
261  xm08= +0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
262 ! TERM 2-6 (09)
263  xm09= +0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
264 ! TERM 3-4 (10)
265  xm10= +0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
266 ! TERM 3-5 (11)
267  xm11= +0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
268 ! TERM 3-6 (12)
269  xm12= +0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
270 ! TERME 4-5 (13)
271  xm13=xm13 +0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
272 ! TERM 4-6 (14)
273  xm14=xm14 +0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
274 ! TERM 5-6 (15)
275  xm15=xm15 +0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
276 !
277 ! TERM 2-1 (16)
278  xm16=xm16 +0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
279 ! TERM 3-1 (17)
280  xm17=xm17 +0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
281 ! TERM 4-1 (18)
282  xm18= 0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
283 ! TERM 5-1 (19)
284  xm19= 0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
285 ! TERM 6-1 (20)
286  xm20= 0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
287 ! TERM 3-2 (21)
288  xm21=xm21 +0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
289 ! TERM 4-2 (22)
290  xm22= 0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
291 ! TERM 5-2 (23)
292  xm23= 0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
293 ! TERM 6-2 (24)
294  xm24= 0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
295 ! TERM 4-3 (25)
296  xm25= 0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
297 ! TERM 5-3 (26)
298  xm26= 0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
299 ! TERM 5-4 (28)
300  xm28=xm28 +0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
301 ! TERM 6-3 (27)
302  xm27= 0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
303 ! TERM 6-4 (29)
304  xm29=xm29 +0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
305 ! TERM 6-5 (30)
306  xm30=xm30 +0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
307 !
308  w1(ielem) = + xm01 * ( u(i2)-u(i1) )
309  & + xm02 * ( u(i3)-u(i1) )
310  & + xm03 * ( u(i4)-u(i1) )
311  & + xm04 * ( u(i5)-u(i1) )
312  & + xm05 * ( u(i6)-u(i1) )
313  w2(ielem) = + xm16 * ( u(i1)-u(i2) )
314  & + xm06 * ( u(i3)-u(i2) )
315  & + xm07 * ( u(i4)-u(i2) )
316  & + xm08 * ( u(i5)-u(i2) )
317  & + xm09 * ( u(i6)-u(i2) )
318  w3(ielem) = + xm17 * ( u(i1)-u(i3) )
319  & + xm21 * ( u(i2)-u(i3) )
320  & + xm10 * ( u(i4)-u(i3) )
321  & + xm11 * ( u(i5)-u(i3) )
322  & + xm12 * ( u(i6)-u(i3) )
323  w4(ielem) = + xm18 * ( u(i1)-u(i4) )
324  & + xm22 * ( u(i2)-u(i4) )
325  & + xm25 * ( u(i3)-u(i4) )
326  & + xm13 * ( u(i5)-u(i4) )
327  & + xm14 * ( u(i6)-u(i4) )
328  w5(ielem) = + xm19 * ( u(i1)-u(i5) )
329  & + xm23 * ( u(i2)-u(i5) )
330  & + xm26 * ( u(i3)-u(i5) )
331  & + xm28 * ( u(i4)-u(i5) )
332  & + xm15 * ( u(i6)-u(i5) )
333  w6(ielem) = + xm20 * ( u(i1)-u(i6) )
334  & + xm24 * ( u(i2)-u(i6) )
335  & + xm27 * ( u(i3)-u(i6) )
336  & + xm29 * ( u(i4)-u(i6) )
337  & + xm30 * ( u(i5)-u(i6) )
338 !
339  ENDDO
340 !
341 !-----------------------------------------------------------------------
342 !
343  ELSEIF(formul(14:16).EQ.'VER') THEN
344 !
345 !-----------------------------------------------------------------------
346 !
347  DO ielem = 1 , nelem
348 !
349  i1 = ikle1(ielem)
350  i2 = ikle2(ielem)
351  i3 = ikle3(ielem)
352  i4 = ikle4(ielem)
353  i5 = ikle5(ielem)
354  i6 = ikle6(ielem)
355 !
356  h1 = z(i4) - z(i1)
357  h2 = z(i5) - z(i2)
358  h3 = z(i6) - z(i3)
359 !
360 ! X2 = X(I2)-X(I1)
361 ! X3 = X(I3)-X(I1)
362 ! Y2 = Y(I2)-Y(I1)
363 ! Y3 = Y(I3)-Y(I1)
364 !
365  x2=x(ielem,2)
366  x3=x(ielem,3)
367  y2=y(ielem,2)
368  y3=y(ielem,3)
369 !
370  IF((inchyd.AND.max(z(i1),z(i2),z(i3)).GT.
371  & min(z(i4),z(i5),z(i6))) .OR.
372  & h1.LT.chouia.OR.h2.LT.chouia.OR.h3.LT.chouia ) THEN
373  nf1=0.d0
374  nf2=0.d0
375  nf3=0.d0
376  ng1=0.d0
377  ng2=0.d0
378  ng3=0.d0
379  nh1=0.d0
380  nh2=0.d0
381  nh3=0.d0
382  nuxmoy=0.d0
383  nuymoy=0.d0
384  ELSE
385 ! SUIVANT LES CAS (II4=I1 OU I4, ETC.)
386 ! VARIANTE AVEC VISCOSITE VERTICALE LINEAIRE (II4=I4)
387 ! VARIANTE AVEC VISCOSITE VERTICALE P0 SUR LA VERTICALE (II4=I1)
388  nf1=0.5d0*(f(i1)+f(i4))/h1
389  nf2=0.5d0*(f(i2)+f(i5))/h2
390  nf3=0.5d0*(f(i3)+f(i6))/h3
391  ng1=0.5d0*(g(i1)+g(i4))/h1
392  ng2=0.5d0*(g(i2)+g(i5))/h2
393  ng3=0.5d0*(g(i3)+g(i6))/h3
394  nh1=0.5d0*(h(i1)+h(i4))/h1
395  nh2=0.5d0*(h(i2)+h(i5))/h2
396  nh3=0.5d0*(h(i3)+h(i6))/h3
397 ! TAKING INTO ACCOUNT AN AVERAGE DIFFUSION ON THE PRISM
398  nuxmoy=(f(i1)+f(i2)+f(i3)+f(i4)+f(i5)+f(i6))/6.d0
399  nuymoy=(g(i1)+g(i2)+g(i3)+g(i4)+g(i5)+g(i6))/6.d0
400  ENDIF
401 !
402 ! AVERAGE OF NORMAL VECTOR TO PLANES (NOT NORMED)
403 ! ONE CAN CHECK THAT WE GET -1 0 OR 0 -1 WITH Z=X OR Z=Y
404 ! NPX=-DZ/DX NPY=-DZ/DY
405 !
406  npxl=-0.5d0*(y2*(z(i1)-z(i3))+y3*(z(i2)-z(i1)))
407  npxu=-0.5d0*(y2*(z(i4)-z(i6))+y3*(z(i5)-z(i4)))
408  npyl=-0.5d0*(x2*(z(i3)-z(i1))+x3*(z(i1)-z(i2)))
409  npyu=-0.5d0*(x2*(z(i6)-z(i4))+x3*(z(i4)-z(i5)))
410  npx=0.5d0*(npxl+npxu)/surfac(ielem)
411  npy=0.5d0*(npyl+npyu)/surfac(ielem)
412  npx2=npx**2
413  npy2=npy**2
414 !
415  npx=npx*nuxmoy
416  npy=npy*nuymoy
417 !
418 ! 2D GRADIENT MATRIX GX(3,3) AND GY(3,3)
419 ! BUT GX(1,J)=GX(2,J)=GX(3,J)
420 ! AND GY(1,J)=GY(2,J)=GY(3,J), HENCE ONLY GX(3) AND GY(3)
421 !
422  gx(2) = y3*xs06
423  gx(3) = - y2*xs06
424  gx(1) = - gx(2) - gx(3)
425 !
426  gy(2) = - x3 * xs06
427  gy(3) = x2 * xs06
428  gy(1) = - gy(2) - gy(3)
429 !
430 !
431 ! TERMS 2
432 !
433 ! TERM 1-2 (01)
434  xm01=0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
435 ! TERM 1-3 (02)
436  xm02=0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
437 ! TERM 1-4 (03)
438  xm03=0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
439 ! TERM 1-5 (04)
440  xm04=0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
441 ! TERM 1-6 (05)
442  xm05=0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
443 ! TERM 2-3 (06)
444  xm06=0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
445 ! TERM 2-4 (07)
446  xm07=0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
447 ! TERM 2-5 (08)
448  xm08=0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
449 ! TERM 2-6 (09)
450  xm09=0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
451 ! TERM 3-4 (10)
452  xm10=0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
453 ! TERM 3-5 (11)
454  xm11=0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
455 ! TERM 3-6 (12)
456  xm12=0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
457 ! TERM 4-5 (13)
458  xm13=0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
459 ! TERM 4-6 (14)
460  xm14=0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
461 ! TERM 5-6 (15)
462  xm15=0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
463 !
464 ! TERM 2-1 (16)
465  xm16=0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
466 ! TERM 3-1 (17)
467  xm17=0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
468 ! TERM 4-1 (18)
469  xm18=0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
470 ! TERM 5-1 (19)
471  xm19=0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
472 ! TERM 6-1 (20)
473  xm20=0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
474 ! TERM 3-2 (21)
475  xm21=0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
476 ! TERM 4-2 (22)
477  xm22=0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
478 ! TERM 5-2 (23)
479  xm23=0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
480 ! TERM 6-2 (24)
481  xm24=0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
482 ! TERM 4-3 (25)
483  xm25=0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
484 ! TERM 5-3 (26)
485  xm26=0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
486 ! TERM 6-3 (27)
487  xm27=0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
488 ! TERM 5-4 (28)
489  xm28=0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
490 ! TERM 6-4 (29)
491  xm29=0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
492 ! TERM 6-5 (30)
493  xm30=0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
494 !
495 ! TERM 4
496 !
497  xm03=xm03-(npx2*nf1+npy2*ng1+nh1)*surfac(ielem)*xs03
498  xm08=xm08-(npx2*nf2+npy2*ng2+nh2)*surfac(ielem)*xs03
499  xm12=xm12-(npx2*nf3+npy2*ng3+nh3)*surfac(ielem)*xs03
500  xm18=xm18-(npx2*nf1+npy2*ng1+nh1)*surfac(ielem)*xs03
501  xm23=xm23-(npx2*nf2+npy2*ng2+nh2)*surfac(ielem)*xs03
502  xm27=xm27-(npx2*nf3+npy2*ng3+nh3)*surfac(ielem)*xs03
503 !
504  w1(ielem) = + xm01 * ( u(i2)-u(i1) )
505  & + xm02 * ( u(i3)-u(i1) )
506  & + xm03 * ( u(i4)-u(i1) )
507  & + xm04 * ( u(i5)-u(i1) )
508  & + xm05 * ( u(i6)-u(i1) )
509  w2(ielem) = + xm16 * ( u(i1)-u(i2) )
510  & + xm06 * ( u(i3)-u(i2) )
511  & + xm07 * ( u(i4)-u(i2) )
512  & + xm08 * ( u(i5)-u(i2) )
513  & + xm09 * ( u(i6)-u(i2) )
514  w3(ielem) = + xm17 * ( u(i1)-u(i3) )
515  & + xm21 * ( u(i2)-u(i3) )
516  & + xm10 * ( u(i4)-u(i3) )
517  & + xm11 * ( u(i5)-u(i3) )
518  & + xm12 * ( u(i6)-u(i3) )
519  w4(ielem) = + xm18 * ( u(i1)-u(i4) )
520  & + xm22 * ( u(i2)-u(i4) )
521  & + xm25 * ( u(i3)-u(i4) )
522  & + xm13 * ( u(i5)-u(i4) )
523  & + xm14 * ( u(i6)-u(i4) )
524  w5(ielem) = + xm19 * ( u(i1)-u(i5) )
525  & + xm23 * ( u(i2)-u(i5) )
526  & + xm26 * ( u(i3)-u(i5) )
527  & + xm28 * ( u(i4)-u(i5) )
528  & + xm15 * ( u(i6)-u(i5) )
529  w6(ielem) = + xm20 * ( u(i1)-u(i6) )
530  & + xm24 * ( u(i2)-u(i6) )
531  & + xm27 * ( u(i3)-u(i6) )
532  & + xm29 * ( u(i4)-u(i6) )
533  & + xm30 * ( u(i5)-u(i6) )
534 !
535  ENDDO
536 !
537 !-----------------------------------------------------------------------
538 !
539  ELSE
540 !
541 !-----------------------------------------------------------------------
542 !
543  WRITE(lu,202) formul
544 202 FORMAT(1x,'VC02PP_STAR (BIEF) :',/,
545  & 1x,'HOR OR VER LACKING AT THE END OF THE FORMULA : ',a16)
546  CALL plante(1)
547  stop
548 !
549 !-----------------------------------------------------------------------
550 !
551  ENDIF
552 !
553 !-----------------------------------------------------------------------
554 !
555  RETURN
556  END
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
vc02pp_star
subroutine vc02pp_star(XMUL, SF, SG, SH, SU, F, G, H, U, X, Y, Z, SURFAC, IKLE1, IKLE2, IKLE3, IKLE4, IKLE5, IKLE6, NELEM, NELMAX, W1, W2, W3, W4, W5, W6, FORMUL)
Definition: vc02pp_star.f:9
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:50 for The TELEMAC-MASCARET system by doxygen 1.8.13