The TELEMAC-MASCARET system  trunk
vc08cc.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE vc08cc
3 ! *****************
4 !
5  &(xmul,sf,su,sv,f,u,v,xel,yel,
6  & ikle1,ikle2,ikle3,ikle4,ikle5,ikle6,nelem,nelmax,
7  & w1,w2,w3,w4,w5,w6,formul)
8 !
9 !***********************************************************************
10 ! BIEF V6P1 21/08/2010
11 !***********************************************************************
12 !
13 !brief COMPUTES THE FOLLOWING VECTOR IN FINITE ELEMENTS:
14 !code
15 !+ / DF DF
16 !+ V = XMUL / PSII * ( U -- + V -- ) D(OMEGA)
17 !+ I /OMEGA DX DY
18 !+
19 !+ PSI(I) IS A BASE OF TYPE P2
20 !
21 !warning THE JACOBIAN MUST BE POSITIVE
22 !warning THE RESULT IS IN W IN NOT ASSEMBLED FORM
23 !
24 !history A FROEHLY
25 !+ 01/07/08
26 !+ V5P9
27 !+
28 !
29 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
30 !+ 13/07/2010
31 !+ V6P0
32 !+ Translation of French comments within the FORTRAN sources into
33 !+ English comments
34 !
35 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
36 !+ 21/08/2010
37 !+ V6P0
38 !+ Creation of DOXYGEN tags for automated documentation and
39 !+ cross-referencing of the FORTRAN sources
40 !
41 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
42 !| F |-->| FUNCTION USED IN THE VECTOR FORMULA
43 !| FORMUL |-->| STRING WITH FORMULA OF VECTOR
44 !| IKLE1 |-->| FIRST POINT OF TRIANGLES
45 !| IKLE2 |-->| SECOND POINT OF TRIANGLES
46 !| IKLE3 |-->| THIRD POINT OF TRIANGLES
47 !| IKLE4 |-->| FOURTH POINT OF TRIANGLES (QUADRATIC)
48 !| IKLE5 |-->| FIFTH POINT OF TRIANGLES (QUADRATIC)
49 !| IKLE6 |-->| SIXTH POINT OF TRIANGLES (QUADRATIC)
50 !| NELEM |-->| NUMBER OF ELEMENTS
51 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
52 !| SF |-->| BIEF_OBJ STRUCTURE OF F
53 !| SU |-->| BIEF_OBJ STRUCTURE OF U
54 !| SV |-->| BIEF_OBJ STRUCTURE OF V
55 !| U |-->| FUNCTION USED IN THE VECTOR FORMULA
56 !| V |-->| FUNCTION USED IN THE VECTOR FORMULA
57 !| W1 |<--| RESULT IN NON ASSEMBLED FORM
58 !| W2 |<--| RESULT IN NON ASSEMBLED FORM
59 !| W3 |<--| RESULT IN NON ASSEMBLED FORM
60 !| W4 |<--| RESULT IN NON ASSEMBLED FORM
61 !| W5 |<--| RESULT IN NON ASSEMBLED FORM
62 !| W6 |<--| RESULT IN NON ASSEMBLED FORM
63 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
64 !| XMUL |-->| MULTIPLICATION COEFFICIENT
65 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
66 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
67 !
68  USE bief, ex_vc08cc => vc08cc
69 !
70  USE declarations_special
71  IMPLICIT NONE
72 !
73 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
74 !
75  INTEGER, INTENT(IN) :: NELEM,NELMAX
76  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax*3),YEL(nelmax*3),XMUL
77 ! W1 IS ALSO USED AS 1-DIMENSIONAL FOR ALL W
78  DOUBLE PRECISION, INTENT(INOUT) :: W1(6*nelmax),W2(nelmax)
79  DOUBLE PRECISION, INTENT(INOUT) :: W3(nelmax),W4(nelmax)
80  DOUBLE PRECISION, INTENT(INOUT) :: W5(nelmax),W6(nelmax)
81 ! IKLE1 IS ALSO USED AS A 1-DIMENSIONAL IKLE
82  INTEGER, INTENT(IN) :: IKLE1(6*nelmax),IKLE2(nelmax)
83  INTEGER, INTENT(IN) :: IKLE3(nelmax),IKLE4(nelmax)
84  INTEGER, INTENT(IN) :: IKLE5(nelmax),IKLE6(nelmax)
85  CHARACTER(LEN=16), INTENT(IN) :: FORMUL
86 !
87 ! STRUCTURES OF F, G, H, U, V, W AND REAL DATA
88 !
89  TYPE(bief_obj), INTENT(IN) :: SF,SU,SV
90  DOUBLE PRECISION, INTENT(IN) :: F(*),U(*),V(*)
91 !
92 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
93 !
94  INTEGER IELEM,IELMF,IELMU,IELMV
95 !
96  DOUBLE PRECISION K1,K2,K3
97 !
98 !-----------------------------------------------------------------------
99 !
100  DOUBLE PRECISION X2,Y2,X3,Y3,F1,F2,F3,F4,F5,F6
101  DOUBLE PRECISION U1,U2,U3,U4,U5,U6,V1,V2,V3,V4,V5,V6
102  DOUBLE PRECISION ANS1,SUR6
103  DOUBLE PRECISION PHIT,USUR2,VSUR2,XSU90,XSU360,XSU630,XSU2520
104  DOUBLE PRECISION L12,L13,L21,L23,L31,L32,BETAN1,BETAN2,BETAN3
105  INTRINSIC max,min
106 !
107 !-----------------------------------------------------------------------
108 !
109  sur6 = 1.d0 / 6.d0
110  xsu90 = xmul/90.d0
111  xsu360 = xmul/360.d0
112  xsu630 = xmul/630.d0
113  xsu2520= xmul/2520.d0
114 !
115  ielmf=sf%ELM
116  ielmu=su%ELM
117  ielmv=sv%ELM
118 !
119 !-----------------------------------------------------------------------
120 !
121 ! FUNCTION F AND VECTOR U ARE P2
122 !
123  IF(ielmf.EQ.13.AND.ielmu.EQ.13.AND.ielmv.EQ.13) THEN
124 !
125  IF(formul(14:16).EQ.'PSI') THEN
126 !
127 ! PSI SCHEME P1 AND LINEAR INTERPOLATION
128 !
129  DO ielem = 1 , nelem
130 !
131  x2 = xel(ielem+nelmax)
132  x3 = xel(ielem+2*nelmax)
133  y2 = yel(ielem+nelmax)
134  y3 = yel(ielem+2*nelmax)
135 !
136  f1 = f(ikle1(ielem))
137  f2 = f(ikle2(ielem))
138  f3 = f(ikle3(ielem))
139 !
140  u1 = u(ikle1(ielem))
141  u2 = u(ikle2(ielem))
142  u3 = u(ikle3(ielem))
143  v1 = v(ikle1(ielem))
144  v2 = v(ikle2(ielem))
145  v3 = v(ikle3(ielem))
146 !
147  usur2 = (u1+u2+u3)*sur6
148  vsur2 = (v1+v2+v3)*sur6
149 !
150  k1 = usur2 * (y2-y3) - vsur2 * (x2-x3)
151  k2 = usur2 * (y3 ) - vsur2 * (x3 )
152  k3 = usur2 * ( -y2) - vsur2 * ( -x2)
153 !
154  l12 = max( min(k1,-k2) , 0.d0 )
155  l13 = max( min(k1,-k3) , 0.d0 )
156  l21 = max( min(k2,-k1) , 0.d0 )
157  l23 = max( min(k2,-k3) , 0.d0 )
158  l31 = max( min(k3,-k1) , 0.d0 )
159  l32 = max( min(k3,-k2) , 0.d0 )
160 !
161  betan1 = l12*(f1-f2) + l13*(f1-f3)
162  betan2 = l21*(f2-f1) + l23*(f2-f3)
163  betan3 = l31*(f3-f1) + l32*(f3-f2)
164 !
165  phit = betan1 + betan2 + betan3
166 !
167  IF(phit.GT.0.d0) THEN
168  w1(ielem) = xmul * max( min( betan1, phit),0.d0 )
169  w2(ielem) = xmul * max( min( betan2, phit),0.d0 )
170  w3(ielem) = xmul * max( min( betan3, phit),0.d0 )
171  ELSE
172  w1(ielem) = - xmul * max( min(-betan1,-phit),0.d0 )
173  w2(ielem) = - xmul * max( min(-betan2,-phit),0.d0 )
174  w3(ielem) = - xmul * max( min(-betan3,-phit),0.d0 )
175  ENDIF
176  w4(ielem) = (w1(ielem)+ w2(ielem))/2.d0
177  w5(ielem) = (w2(ielem)+ w3(ielem))/2.d0
178  w6(ielem) = (w3(ielem)+ w1(ielem))/2.d0
179 !
180  ENDDO ! IELEM
181 !
182  ELSE
183 !
184 ! CLASSICAL COMPUTATION
185 !
186  DO ielem = 1 , nelem
187 !
188  x2 = xel(ielem+nelmax)
189  x3 = xel(ielem+2*nelmax)
190  y2 = yel(ielem+nelmax)
191  y3 = yel(ielem+2*nelmax)
192 !
193  u1 = u(ikle1(ielem))
194  u2 = u(ikle2(ielem))
195  u3 = u(ikle3(ielem))
196  u4 = u(ikle4(ielem))
197  u5 = u(ikle5(ielem))
198  u6 = u(ikle6(ielem))
199 !
200  v1 = v(ikle1(ielem))
201  v2 = v(ikle2(ielem))
202  v3 = v(ikle3(ielem))
203  v4 = v(ikle4(ielem))
204  v5 = v(ikle5(ielem))
205  v6 = v(ikle6(ielem))
206 !
207  f1 = f(ikle1(ielem))
208  f2 = f(ikle2(ielem)) - f1
209  f3 = f(ikle3(ielem)) - f1
210  f4 = f(ikle4(ielem)) - f1
211  f5 = f(ikle5(ielem)) - f1
212  f6 = f(ikle6(ielem)) - f1
213 !
214  ans1 =-20.d0*y2*f3*u5-64.d0*x2*f6*v4+16.d0*y2*f4*u6+
215  & 16.d0*x3*f6*v4+4.d0*y2*f5*u3-32.d0*y3*f4*u5-16.d0*x2*f4*v6-
216  & 16.d0*y2*f3*u4+32.d0*x3*f4*v5-32.d0*x2*f4*v4-
217  & 24.d0*y2*f4*u1+24.d0*y2*f5*u1-80.d0*y3*f4*u4-
218  & 24.d0*x2*f5*v1+16.d0*x2*f3*v4-9.d0*x3*f2*v2+
219  & 80.d0*x3*f4*v4-16.d0*y3*f6*u4+9.d0*x2*f3*v3+
220  & 4.d0*y3*f6*u2-48.d0*x2*f4*v5-32.d0*y3*f6*u6-
221  & 32.d0*y2*f5*u4+32.d0*y3*f2*u4-16.d0*y2*f5*u6+
222  & 11.d0*x3*f2*v3+24.d0*y3*f6*u1+16.d0*x2*f5*v6-
223  & 64.d0*y3*f4*u6-48.d0*y3*f6*u5+4.d0*x2*f4*v3-
224  & 48.d0*x3*f5*v5-18.d0*x3*f2*v1+16.d0*y3*f2*u6+
225  & 32.d0*y2*f4*u4+20.d0*y3*f4*u3+11.d0*y2*f3*u2+
226  & 20.d0*y3*f2*u5-4.d0*x3*f6*v2+32.d0*x2*f5*v4-11.d0*x2*f3*v2-
227  & 32.d0*x3*f2*v4-9.d0*y2*f3*u3+96.d0*y2*f6*u1+9.d0*y3*f2*u2-
228  & 24.d0*x3*f6*v1+32.d0*y3*f5*u6-24.d0*y3*f5*u1+
229  & 48.d0*y2*f4*u5+96.d0*x3*f4*v1-48.d0*y2*f5*u5+
230  & 48.d0*x2*f5*v5-32.d0*x3*f5*v6-20.d0*x3*f2*v5+4.d0*x3*f5*v2+
231  & 64.d0*x3*f4*v6-16.d0*x2*f4*v2+18.d0*x2*f3*v1+
232  & 24.d0*x3*f5*v1-96.d0*x2*f6*v1+32.d0*x3*f6*v6-
233  & 11.d0*y3*f2*u3-4.d0*x2*f5*v3+20.d0*x2*f3*v5+16.d0*y3*f5*u3
234  w1(ielem) =( -16.d0*y3*f6*u3-20.d0*x3*f4*v3-16.d0*x3*f5*v3-
235  & 16.d0*y2*f5*u2+16.d0*y2*f4*u2+16.d0*x2*f5*v2-
236  & 20.d0*y2*f6*u2+32.d0*y2*f6*u5+20.d0*x2*f6*v2+
237  & 80.d0*y2*f6*u6-96.d0*y3*f4*u1-80.d0*x2*f6*v6-
238  & 32.d0*x2*f6*v5-18.d0*y2*f3*u1-4.d0*y2*f4*u3+
239  & 32.d0*x2*f3*v6+16.d0*x3*f6*v3-32.d0*y2*f3*u6-
240  & 16.d0*x3*f2*v6-4.d0*y3*f5*u2+24.d0*x2*f4*v1+
241  & 64.d0*y2*f6*u4+16.d0*y3*f5*u4+48.d0*y3*f5*u5+
242  & 48.d0*x3*f6*v5+18.d0*y3*f2*u1-16.d0*x3*f5*v4+
243  & ans1) * (-xsu2520)
244 !
245  ans1 = 32.d0*y2*f3*u5-16.d0*x2*f6*v4-16.d0*y2*f4*u6-16.d0*x3*f6*v4
246  &+16.d0*y2*f5*u3-64.d0*y3*f4*u5+16.d0*x2*f4*v6+16.d0*y2*f3*u4
247  &+64.d0*x3*f4*v5-48.d0*x2*f4*v4-16.d0*y2*f4*u1+16.d0*y2*f5*u1
248  &-80.d0*y3*f4*u4-96.d0*y3*f4*u2-16.d0*x2*f5*v1-16.d0*x2*f3*v4
249  &-78.d0*x3*f2*v2+80.d0*x3*f4*v4+16.d0*y3*f6*u4-9.d0*x2*f3*v3
250  &-24.d0*y3*f6*u2-48.d0*x2*f4*v5+48.d0*y3*f6*u6-48.d0*y2*f5*u4
251  &+48.d0*y3*f2*u4+16.d0*y2*f5*u6+9.d0*x3*f2*v3-4.d0*y3*f6*u1
252  &-16.d0*x2*f5*v6-32.d0*y3*f4*u6+32.d0*y3*f6*u5+16.d0*x2*f4*v3
253  &+32.d0*x3*f5*v5+9.d0*x3*f2*v1+12.d0*y3*f2*u6-20.d0*y2*f6*u3
254  &+48.d0*y2*f4*u4+20.d0*y3*f4*u3+18.d0*y2*f3*u2+48.d0*y3*f2*u5
255  &+24.d0*x3*f6*v2+48.d0*x2*f5*v4-18.d0*x2*f3*v2-48.d0*x3*f2*v4
256  &+96.d0*x3*f4*v2+9.d0*y2*f3*u3+20.d0*y2*f6*u1+78.d0*y3*f2*u2
257  &+4.d0*x3*f6*v1-48.d0*y3*f5*u6+4.d0*y3*f5*u1+48.d0*y2*f4*u5
258  &-48.d0*y2*f5*u5+48.d0*x2*f5*v5+48.d0*x3*f5*v6-48.d0*x3*f2*v5
259  &-24.d0*x3*f5*v2+32.d0*x3*f4*v6-120.d0*x2*f4*v2+11.d0*x2*f3*v1
260  &-4.d0*x3*f5*v1-20.d0*x2*f6*v1-48.d0*x3*f6*v6-9.d0*y3*f2*u3
261  &-16.d0*x2*f5*v3-32.d0*x2*f3*v5-16.d0*y3*f5*u3+16.d0*y3*f6*u3
262  &-20.d0*x3*f4*v3+16.d0*x3*f5*v3-120.d0*y2*f5*u2+120.d0*y2*f4*u2
263  &+120.d0*x2*f5*v2-16.d0*y2*f6*u5+16.d0*x2*f6*v5
264  w2(ielem) = (-11.d0*y2*f3*u1-16.d0*y2*f4*u3-20.d0*x2*f3*v6-
265  & 16.d0*x3*f6*v3+20.d0*y2*f3*u6-12.d0*x3*f2*v6+
266  & 24.d0*y3*f5*u2+16.d0*x2*f4*v1+16.d0*y2*f6*u4-
267  & 16.d0*y3*f5*u4-32.d0*y3*f5*u5-32.d0*x3*f6*v5+
268  & 20.d0*x2*f6*v3-9.d0*y3*f2*u1+16.d0*x3*f5*v4+
269  & ans1) * xsu2520
270 !
271  ans1 = -48.d0*y2*f3*u5-32.d0*x2*f6*v4-16.d0*y2*f4*u6-
272  & 16.d0*x3*f6*v4-24.d0*y2*f5*u3+16.d0*y3*f4*u5+
273  & 16.d0*x2*f4*v6-12.d0*y2*f3*u4-16.d0*x3*f4*v5+
274  & 48.d0*x2*f4*v4+4.d0*y2*f4*u1-4.d0*y2*f5*u1+20.d0*y3*f4*u2+
275  & 4.d0*x2*f5*v1+12.d0*x2*f3*v4+9.d0*x3*f2*v2+16.d0*y3*f6*u4+
276  & 78.d0*x2*f3*v3+16.d0*y3*f6*u2+32.d0*x2*f4*v5-
277  & 48.d0*y3*f6*u6+48.d0*y2*f5*u4-20.d0*y3*f2*u4+
278  & 16.d0*y2*f5*u6+18.d0*x3*f2*v3+16.d0*y3*f6*u1-
279  & 16.d0*x2*f5*v6-16.d0*y3*f4*u6-48.d0*y3*f6*u5-
280  & 24.d0*x2*f4*v3-48.d0*x3*f5*v5-11.d0*x3*f2*v1-
281  & 16.d0*y3*f2*u6+96.d0*y2*f6*u3-48.d0*y2*f4*u4+
282  & 9.d0*y2*f3*u2-32.d0*y3*f2*u5-16.d0*x3*f6*v2-
283  & 48.d0*x2*f5*v4-9.d0*x2*f3*v2+20.d0*x3*f2*v4-
284  & 20.d0*x3*f4*v2-78.d0*y2*f3*u3-9.d0*y3*f2*u2-
285  & 16.d0*x3*f6*v1+48.d0*y3*f5*u6-16.d0*y3*f5*u1-
286  & 32.d0*y2*f4*u5+20.d0*x3*f4*v1+32.d0*y2*f5*u5-
287  & 32.d0*x2*f5*v5-48.d0*x3*f5*v6+32.d0*x3*f2*v5+
288  & 16.d0*x3*f5*v2+16.d0*x3*f4*v6+16.d0*x2*f4*v2-
289  & 9.d0*x2*f3*v1+16.d0*x3*f5*v1+48.d0*x3*f6*v6-
290  & 18.d0*y3*f2*u3+24.d0*x2*f5*v3+48.d0*x2*f3*v5
291  w3(ielem) =(120.d0*y3*f5*u3-120.d0*y3*f6*u3-120.d0*x3*f5*v3+
292  & 16.d0*y2*f5*u2-16.d0*y2*f4*u2-16.d0*x2*f5*v2-
293  & 20.d0*y2*f6*u2+64.d0*y2*f6*u5+20.d0*x2*f6*v2+
294  & 80.d0*y2*f6*u6-20.d0*y3*f4*u1-80.d0*x2*f6*v6-
295  & 64.d0*x2*f6*v5+9.d0*y2*f3*u1+24.d0*y2*f4*u3+
296  & 48.d0*x2*f3*v6+120.d0*x3*f6*v3-48.d0*y2*f3*u6+
297  & 16.d0*x3*f2*v6-16.d0*y3*f5*u2-4.d0*x2*f4*v1+
298  & 32.d0*y2*f6*u4-16.d0*y3*f5*u4+48.d0*y3*f5*u5+
299  & 48.d0*x3*f6*v5-96.d0*x2*f6*v3+11.d0*y3*f2*u1+
300  & 16.d0*x3*f5*v4+ans1) * xsu2520
301 !
302  ans1 = 4.d0*y2*f3*u5-64.d0*x2*f6*v4-32.d0*y2*f4*u6-
303  & 32.d0*x3*f6*v4-12.d0*y2*f5*u3+16.d0*y3*f4*u5+
304  & 32.d0*x2*f4*v6-24.d0*y2*f3*u4-16.d0*x3*f4*v5+
305  & 96.d0*x2*f4*v4+8.d0*y2*f4*u1-8.d0*y2*f5*u1+
306  & 20.d0*y3*f4*u2+8.d0*x2*f5*v1+24.d0*x2*f3*v4+
307  & 12.d0*x3*f2*v2+32.d0*y3*f6*u4-3.d0*x2*f3*v3-
308  & 4.d0*y3*f6*u2+48.d0*x2*f4*v5+32.d0*y3*f6*u6+
309  & 96.d0*y2*f5*u4-40.d0*y3*f2*u4+32.d0*y2*f5*u6-
310  & 5.d0*x3*f2*v3-4.d0*y3*f6*u1-32.d0*x2*f5*v6-16.d0*y3*f4*u6+
311  & 32.d0*y3*f6*u5-12.d0*x2*f4*v3+32.d0*x3*f5*v5-
312  & 8.d0*x3*f2*v1-4.d0*y3*f2*u6-8.d0*y2*f6*u3-96.d0*y2*f4*u4-
313  & 4.d0*y2*f3*u2-20.d0*y3*f2*u5+4.d0*x3*f6*v2-96.d0*x2*f5*v4+
314  & 4.d0*x2*f3*v2+40.d0*x3*f2*v4-20.d0*x3*f4*v2+3.d0*y2*f3*u3+
315  & 16.d0*y2*f6*u1-12.d0*y3*f2*u2+4.d0*x3*f6*v1-32.d0*y3*f5*u6+
316  & 4.d0*y3*f5*u1-48.d0*y2*f4*u5+20.d0*x3*f4*v1+48.d0*y2*f5*u5-
317  & 48.d0*x2*f5*v5+32.d0*x3*f5*v6+20.d0*x3*f2*v5-4.d0*x3*f5*v2+
318  & 16.d0*x3*f4*v6+12.d0*x2*f4*v2+4.d0*x2*f3*v1-4.d0*x3*f5*v1-
319  & 16.d0*x2*f6*v1-32.d0*x3*f6*v6+5.d0*y3*f2*u3+12.d0*x2*f5*v3-
320  & 4.d0*x2*f3*v5+4.d0*y3*f5*u3-4.d0*y3*f6*u3-4.d0*x3*f5*v3+
321  & 12.d0*y2*f5*u2-12.d0*y2*f4*u2-12.d0*x2*f5*v2-4.d0*y2*f6*u2
322  w4(ielem) = (4.d0*x2*f6*v2+16.d0*y2*f6*u6-20.d0*y3*f4*u1-
323  & 16.d0*x2*f6*v6-4.d0*y2*f3*u1+12.d0*y2*f4*u3-
324  & 4.d0*x2*f3*v6+4.d0*x3*f6*v3+4.d0*y2*f3*u6+
325  & 4.d0*x3*f2*v6+4.d0*y3*f5*u2-8.d0*x2*f4*v1+
326  & 64.d0*y2*f6*u4-32.d0*y3*f5*u4-32.d0*y3*f5*u5-
327  & 32.d0*x3*f6*v5+8.d0*x2*f6*v3+8.d0*y3*f2*u1+
328  & 32.d0*x3*f5*v4+ ans1) * (-xsu630)
329 !
330  ans1 = -40.d0*y2*f3*u5+32.d0*y2*f4*u6+32.d0*x3*f6*v4+
331  & 8.d0*y2*f5*u3-64.d0*y3*f4*u5-32.d0*x2*f4*v6-4.d0*y2*f3*u4+
332  & 64.d0*x3*f4*v5-48.d0*x2*f4*v4-12.d0*y2*f4*u1+
333  & 12.d0*y2*f5*u1-16.d0*y3*f4*u4-16.d0*y3*f4*u2-
334  & 12.d0*x2*f5*v1+4.d0*x2*f3*v4-12.d0*x3*f2*v2+
335  & 16.d0*x3*f4*v4-32.d0*y3*f6*u4+12.d0*x2*f3*v3+
336  & 8.d0*y3*f6*u2-96.d0*x2*f4*v5-48.d0*y3*f6*u6-
337  & 48.d0*y2*f5*u4+20.d0*y3*f2*u4-32.d0*y2*f5*u6+
338  & 8.d0*x3*f2*v3+12.d0*y3*f6*u1+32.d0*x2*f5*v6-
339  & 96.d0*y3*f6*u5+8.d0*x2*f4*v3-96.d0*x3*f5*v5+
340  & 5.d0*x3*f2*v1+4.d0*y3*f2*u6+16.d0*y2*f6*u3+48.d0*y2*f4*u4+
341  & 4.d0*y3*f4*u3+8.d0*y2*f3*u2+40.d0*y3*f2*u5-
342  & 8.d0*x3*f6*v2+48.d0*x2*f5*v4-8.d0*x2*f3*v2-
343  & 20.d0*x3*f2*v4+16.d0*x3*f4*v2-12.d0*y2*f3*u3-
344  & 8.d0*y2*f6*u1+12.d0*y3*f2*u2-12.d0*x3*f6*v1+
345  & 48.d0*y3*f5*u6-12.d0*y3*f5*u1+96.d0*y2*f4*u5-
346  & 8.d0*x3*f4*v1-96.d0*y2*f5*u5+96.d0*x2*f5*v5-
347  & 48.d0*x3*f5*v6-40.d0*x3*f2*v5+8.d0*x3*f5*v2-
348  & 12.d0*x2*f4*v2-5.d0*x2*f3*v1+12.d0*x3*f5*v1+
349  & 8.d0*x2*f6*v1+48.d0*x3*f6*v6-8.d0*y3*f2*u3-8.d0*x2*f5*v3
350  w5(ielem) = (40.d0*x2*f3*v5+12.d0*y3*f5*u3-12.d0*y3*f6*u3-
351  & 4.d0*x3*f4*v3-12.d0*x3*f5*v3-12.d0*y2*f5*u2+
352  & 12.d0*y2*f4*u2+12.d0*x2*f5*v2-4.d0*y2*f6*u2+
353  & 64.d0*y2*f6*u5+4.d0*x2*f6*v2+16.d0*y2*f6*u6+
354  & 8.d0*y3*f4*u1-16.d0*x2*f6*v6-64.d0*x2*f6*v5+
355  & 5.d0*y2*f3*u1-8.d0*y2*f4*u3+20.d0*x2*f3*v6+
356  & 12.d0*x3*f6*v3-20.d0*y2*f3*u6-4.d0*x3*f2*v6-
357  & 8.d0*y3*f5*u2+12.d0*x2*f4*v1+32.d0*y3*f5*u4+
358  & 96.d0*y3*f5*u5+96.d0*x3*f6*v5-16.d0*x2*f6*v3-
359  & 5.d0*y3*f2*u1-32.d0*x3*f5*v4+ ans1) * xsu630
360 !
361  ans1 = 20.d0*y2*f3*u5-16.d0*x2*f6*v4-32.d0*y2*f4*u6-
362  & 32.d0*x3*f6*v4-4.d0*y2*f5*u3+32.d0*x2*f4*v6+4.d0*y2*f3*u4+
363  & 32.d0*x2*f4*v4+4.d0*y2*f4*u1-4.d0*y2*f5*u1-16.d0*y3*f4*u4+
364  & 8.d0*y3*f4*u2+4.d0*x2*f5*v1-4.d0*x2*f3*v4+3.d0*x3*f2*v2+
365  & 16.d0*x3*f4*v4+32.d0*y3*f6*u4-12.d0*x2*f3*v3-
366  & 12.d0*y3*f6*u2+32.d0*x2*f4*v5+96.d0*y3*f6*u6+
367  & 32.d0*y2*f5*u4-4.d0*y3*f2*u4+32.d0*y2*f5*u6-
368  & 4.d0*x3*f2*v3-8.d0*y3*f6*u1-32.d0*x2*f5*v6-64.d0*y3*f4*u6+
369  & 48.d0*y3*f6*u5-4.d0*x2*f4*v3+48.d0*x3*f5*v5-4.d0*x3*f2*v1+
370  & 24.d0*y3*f2*u6-20.d0*y2*f6*u3-32.d0*y2*f4*u4+4.d0*y3*f4*u3-
371  & 5.d0*y2*f3*u2-4.d0*y3*f2*u5+12.d0*x3*f6*v2-32.d0*x2*f5*v4+
372  & 5.d0*x2*f3*v2+4.d0*x3*f2*v4-8.d0*x3*f4*v2+12.d0*y2*f3*u3+
373  & 20.d0*y2*f6*u1-3.d0*y3*f2*u2+8.d0*x3*f6*v1-96.d0*y3*f5*u6+
374  & 8.d0*y3*f5*u1-32.d0*y2*f4*u5+16.d0*x3*f4*v1+32.d0*y2*f5*u5-
375  & 32.d0*x2*f5*v5+96.d0*x3*f5*v6+4.d0*x3*f2*v5-12.d0*x3*f5*v2+
376  & 64.d0*x3*f4*v6-4.d0*x2*f4*v2+8.d0*x2*f3*v1-8.d0*x3*f5*v1-
377  & 20.d0*x2*f6*v1-96.d0*x3*f6*v6+4.d0*y3*f2*u3+4.d0*x2*f5*v3-
378  & 20.d0*x2*f3*v5-12.d0*y3*f5*u3+12.d0*y3*f6*u3-4.d0*x3*f4*v3+
379  & 12.d0*x3*f5*v3-4.d0*y2*f5*u2+4.d0*y2*f4*u2+4.d0*x2*f5*v2-
380  & 16.d0*y2*f6*u5-16.d0*y3*f4*u1+16.d0*x2*f6*v5-8.d0*y2*f3*u1
381  w6(ielem) = (4.d0*y2*f4*u3-40.d0*x2*f3*v6-12.d0*x3*f6*v3+
382  & 40.d0*y2*f3*u6-24.d0*x3*f2*v6+12.d0*y3*f5*u2-
383  & 4.d0*x2*f4*v1+16.d0*y2*f6*u4-32.d0*y3*f5*u4-
384  & 48.d0*y3*f5*u5-48.d0*x3*f6*v5+20.d0*x2*f6*v3+
385  & 4.d0*y3*f2*u1+32.d0*x3*f5*v4+ans1) * (-xsu630)
386 !
387  ENDDO ! IELEM
388 !
389  ENDIF
390 !
391 ! FUNCTION F IS P2 AND VECTOR U LINEAR
392 !
393  ELSEIF(ielmf.EQ.13.AND.ielmu.EQ.11.AND.ielmv.EQ.11) THEN
394 !
395  IF(formul(14:16).EQ.'PSI') THEN
396 !
397 ! PSI SCHEME P1 AND LINEAR INTERPOLATION
398 !
399  DO ielem = 1 , nelem
400 !
401  x2 = xel(ielem+nelmax)
402  x3 = xel(ielem+2*nelmax)
403  y2 = yel(ielem+nelmax)
404  y3 = yel(ielem+2*nelmax)
405 !
406  f1 = f(ikle1(ielem))
407  f2 = f(ikle2(ielem))
408  f3 = f(ikle3(ielem))
409 !
410  u1 = u(ikle1(ielem))
411  u2 = u(ikle2(ielem))
412  u3 = u(ikle3(ielem))
413  v1 = v(ikle1(ielem))
414  v2 = v(ikle2(ielem))
415  v3 = v(ikle3(ielem))
416 !
417  usur2 = (u1+u2+u3)*sur6
418  vsur2 = (v1+v2+v3)*sur6
419 !
420  k1 = usur2 * (y2-y3) - vsur2 * (x2-x3)
421  k2 = usur2 * (y3 ) - vsur2 * (x3 )
422  k3 = usur2 * ( -y2) - vsur2 * ( -x2)
423 !
424  l12 = max( min(k1,-k2) , 0.d0 )
425  l13 = max( min(k1,-k3) , 0.d0 )
426  l21 = max( min(k2,-k1) , 0.d0 )
427  l23 = max( min(k2,-k3) , 0.d0 )
428  l31 = max( min(k3,-k1) , 0.d0 )
429  l32 = max( min(k3,-k2) , 0.d0 )
430 !
431  betan1 = l12*(f1-f2) + l13*(f1-f3)
432  betan2 = l21*(f2-f1) + l23*(f2-f3)
433  betan3 = l31*(f3-f1) + l32*(f3-f2)
434 !
435  phit = betan1 + betan2 + betan3
436 !
437  IF(phit.GT.0.d0) THEN
438  w1(ielem) = xmul * max( min( betan1, phit),0.d0 )
439  w2(ielem) = xmul * max( min( betan2, phit),0.d0 )
440  w3(ielem) = xmul * max( min( betan3, phit),0.d0 )
441  ELSE
442  w1(ielem) = - xmul * max( min(-betan1,-phit),0.d0 )
443  w2(ielem) = - xmul * max( min(-betan2,-phit),0.d0 )
444  w3(ielem) = - xmul * max( min(-betan3,-phit),0.d0 )
445  ENDIF
446  w4(ielem) = (w1(ielem)+ w2(ielem))/2.d0
447  w5(ielem) = (w2(ielem)+ w3(ielem))/2.d0
448  w6(ielem) = (w3(ielem)+ w1(ielem))/2.d0
449 !
450  ENDDO ! IELEM
451 !
452  ELSE
453 !
454 ! CLASSICAL COMPUTATION
455 !
456  DO ielem = 1 , nelem
457 !
458  x2 = xel(ielem+nelmax)
459  x3 = xel(ielem+2*nelmax)
460  y2 = yel(ielem+nelmax)
461  y3 = yel(ielem+2*nelmax)
462 !
463  u1 = u(ikle1(ielem))
464  u2 = u(ikle2(ielem))
465  u3 = u(ikle3(ielem))
466  v1 = v(ikle1(ielem))
467  v2 = v(ikle2(ielem))
468  v3 = v(ikle3(ielem))
469 !
470  f1 = f(ikle1(ielem))
471  f2 = f(ikle2(ielem)) - f1
472  f3 = f(ikle3(ielem)) - f1
473  f4 = f(ikle4(ielem)) - f1
474  f5 = f(ikle5(ielem)) - f1
475  f6 = f(ikle6(ielem)) - f1
476 !
477  w1(ielem) = (-4.d0*x3*f6*v2-4.d0*y3*f5*u2-4.d0*y2*f6*u2-
478  & 8.d0*y2*f4*u2-4.d0*x2*f5*v3-4.d0*y2*f4*u3+
479  & 8.d0*y3*f6*u3-x2*f3*v2+4.d0*x3*f5*v2+4.d0*y2*f5*u3+
480  & 4.d0*y3*f6*u2+8.d0*x3*f5*v3+4.d0*x2*f4*v3-
481  & 6.d0*y3*f2*u1+24.d0*x2*f6*v1+y2*f3*u2+
482  & 6.d0*y2*f3*u1+8.d0*x2*f4*v2-8.d0*x2*f5*v2+
483  & 8.d0*y3*f4*u2-24.d0*x3*f4*v1+8.d0*y2*f5*u2-
484  & 5.d0*y3*f2*u2-8.d0*x3*f4*v2-24.d0*y2*f6*u1+
485  & 6.d0*x3*f2*v1-6.d0*x2*f3*v1-y3*f2*u3+5.d0*y2*f3*u3-
486  & 8.d0*y2*f6*u3+8.d0*x2*f6*v3+4.d0*x2*f6*v2-
487  & 4.d0*x3*f4*v3-8.d0*x3*f6*v3-5.d0*x2*f3*v3+
488  & 24.d0*y3*f4*u1+4.d0*y3*f4*u3+5.d0*x3*f2*v2-
489  & 8.d0*y3*f5*u3+x3*f2*v3) * xsu360
490 !
491  w2(ielem)= (-24.d0*y2*f4*u2-8.d0*y3*f6*u3+6.d0*x2*f3*v2-
492  & 4.d0*y3*f6*u1-8.d0*x3*f5*v3+4.d0*x3*f6*v1-
493  & 3.d0*y3*f2*u1+4.d0*x2*f6*v1-6.d0*y2*f3*u2-y2*f3*u1+
494  & 4.d0*y3*f5*u1+24.d0*x2*f4*v2-24.d0*x2*f5*v2+
495  & 24.d0*y3*f4*u2-8.d0*x3*f4*v1+24.d0*y2*f5*u2-
496  & 18.d0*y3*f2*u2-24.d0*x3*f4*v2-4.d0*y2*f6*u1+
497  & 3.d0*x3*f2*v1+x2*f3*v1-3.d0*y3*f2*u3-5.d0*y2*f3*u3-
498  & 4.d0*x3*f5*v1+4.d0*y2*f6*u3-4.d0*x2*f6*v3-
499  & 4.d0*x3*f4*v3+8.d0*x3*f6*v3+5.d0*x2*f3*v3+
500  & 8.d0*y3*f4*u1+4.d0*y3*f4*u3+18.d0*x3*f2*v2+
501  & 8.d0*y3*f5*u3+3.d0*x3*f2*v3) * (-xsu360)
502 !
503  w3(ielem)= (4.d0*x2*f5*v1-4.d0*y2*f6*u2-4.d0*y2*f5*u1+
504  & 8.d0*y2*f4*u2+24.d0*y3*f6*u3-3.d0*x2*f3*v2+
505  & 24.d0*x3*f5*v3+y3*f2*u1+8.d0*x2*f6*v1+4.d0*y2*f4*u1+
506  & 3.d0*y2*f3*u2+3.d0*y2*f3*u1-4.d0*x2*f4*v1-
507  & 8.d0*x2*f4*v2+8.d0*x2*f5*v2-4.d0*y3*f4*u2-
508  & 4.d0*x3*f4*v1-8.d0*y2*f5*u2+5.d0*y3*f2*u2+
509  & 4.d0*x3*f4*v2-8.d0*y2*f6*u1-x3*f2*v1-3.d0*x2*f3*v1+
510  & 6.d0*y3*f2*u3+18.d0*y2*f3*u3-24.d0*y2*f6*u3+
511  & 24.d0*x2*f6*v3+4.d0*x2*f6*v2-24.d0*x3*f6*v3-
512  & 18.d0*x2*f3*v3+4.d0*y3*f4*u1-5.d0*x3*f2*v2-
513  & 24.d0*y3*f5*u3-6.d0*x3*f2*v3)*(-xsu360)
514 !
515  w4(ielem) = (4.d0*x3*f6*v2+8.d0*x2*f5*v1+4.d0*y3*f5*u2-
516  & 4.d0*y2*f6*u2-8.d0*y2*f5*u1+12.d0*y2*f4*u2+
517  & 4.d0*x2*f5*v3+4.d0*y2*f4*u3-4.d0*y3*f6*u3-
518  & 2.d0*x2*f3*v2-4.d0*x3*f5*v2-4.d0*y3*f6*u1-
519  & 4.d0*y2*f5*u3-4.d0*y3*f6*u2-4.d0*x3*f5*v3+
520  & 4.d0*x3*f6*v1-4.d0*x2*f4*v3+2.d0*y3*f2*u1+
521  & 8.d0*x2*f6*v1+8.d0*y2*f4*u1+2.d0*y2*f3*u2+
522  & 2.d0*y2*f3*u1+4.d0*y3*f5*u1-8.d0*x2*f4*v1-
523  & 12.d0*x2*f4*v2+12.d0*x2*f5*v2-4.d0*y3*f4*u2-
524  & 4.d0*x3*f4*v1-12.d0*y2*f5*u2+6.d0*y3*f2*u2+
525  & 4.d0*x3*f4*v2-8.d0*y2*f6*u1-2.d0*x3*f2*v1-
526  & 2.d0*x2*f3*v1+y3*f2*u3-y2*f3*u3-4.d0*x3*f5*v1+
527  & 4.d0*x2*f6*v2+4.d0*x3*f6*v3+x2*f3*v3+4.d0*y3*f4*u1-
528  & 6.d0*x3*f2*v2+4.d0*y3*f5*u3-x3*f2*v3)*xsu90
529 !
530  w5(ielem) = (8.d0*x3*f6*v2+4.d0*x2*f5*v1+8.d0*y3*f5*u2+
531  & 4.d0*y2*f6*u2-4.d0*y2*f5*u1+12.d0*y2*f4*u2+
532  & 8.d0*x2*f5*v3+8.d0*y2*f4*u3-12.d0*y3*f6*u3+
533  & 2.d0*x2*f3*v2-8.d0*x3*f5*v2-4.d0*y3*f6*u1-
534  & 8.d0*y2*f5*u3-8.d0*y3*f6*u2-12.d0*x3*f5*v3+
535  & 4.d0*x3*f6*v1-8.d0*x2*f4*v3+y3*f2*u1+4.d0*y2*f4*u1-
536  & 2.d0*y2*f3*u2-y2*f3*u1+4.d0*y3*f5*u1-4.d0*x2*f4*v1-
537  & 12.d0*x2*f4*v2+12.d0*x2*f5*v2-8.d0*y3*f4*u2-
538  & 12.d0*y2*f5*u2+6.d0*y3*f2*u2+8.d0*x3*f4*v2-x3*f2*v1+
539  & x2*f3*v1+2.d0*y3*f2*u3-6.d0*y2*f3*u3-4.d0*x3*f5*v1+
540  & 8.d0*y2*f6*u3-8.d0*x2*f6*v3-4.d0*x2*f6*v2+
541  & 4.d0*x3*f4*v3+12.d0*x3*f6*v3+6.d0*x2*f3*v3-
542  & 4.d0*y3*f4*u3-6.d0*x3*f2*v2+12.d0*y3*f5*u3-
543  & 2.d0*x3*f2*v3) * xsu90
544  w6(ielem) = (4.d0*x3*f6*v2+4.d0*x2*f5*v1+4.d0*y3*f5*u2-
545  & 4.d0*y2*f5*u1+4.d0*y2*f4*u2+4.d0*x2*f5*v3+
546  & 4.d0*y2*f4*u3-12.d0*y3*f6*u3+x2*f3*v2-
547  & 4.d0*x3*f5*v2-8.d0*y3*f6*u1-4.d0*y2*f5*u3-
548  & 4.d0*y3*f6*u2-12.d0*x3*f5*v3+8.d0*x3*f6*v1-
549  & 4.d0*x2*f4*v3-2.d0*y3*f2*u1+4.d0*x2*f6*v1+
550  & 4.d0*y2*f4*u1-y2*f3*u2-2.d0*y2*f3*u1+8.d0*y3*f5*u1-
551  & 4.d0*x2*f4*v1-4.d0*x2*f4*v2+4.d0*x2*f5*v2-
552  & 8.d0*x3*f4*v1-4.d0*y2*f5*u2+y3*f2*u2-
553  & 4.d0*y2*f6*u1+2.d0*x3*f2*v1+2.d0*x2*f3*v1-
554  & 2.d0*y3*f2*u3-6.d0*y2*f3*u3-8.d0*x3*f5*v1+
555  & 4.d0*y2*f6*u3-4.d0*x2*f6*v3-4.d0*x3*f4*v3+
556  & 12.d0*x3*f6*v3+6.d0*x2*f3*v3+8.d0*y3*f4*u1+
557  & 4.d0*y3*f4*u3-x3*f2*v2+12.d0*y3*f5*u3+2.d0*x3*f2*v3)
558  & * xsu90
559 !
560  ENDDO
561 !
562  ENDIF
563 !
564 !-----------------------------------------------------------------------
565 !
566  ELSE
567 !
568 !-----------------------------------------------------------------------
569 !
570  WRITE(lu,101) ielmf,sf%NAME
571  WRITE(lu,201) ielmu,su%NAME
572  WRITE(lu,301)
573 101 FORMAT(1x,'VC08CC (BIEF) :',/,
574  & 1x,'DISCRETIZATION OF F:',1i6,
575  & 1x,'REAL NAME: ',a6)
576 201 FORMAT(1x,'DISCRETIZATION OF U:',1i6,
577  & 1x,'REAL NAME: ',a6)
578 301 FORMAT(1x,'CASE NOT IMPLEMENTED')
579  CALL plante(1)
580  stop
581 !
582  ENDIF
583 !
584 !-----------------------------------------------------------------------
585 !
586  RETURN
587  END
declarations_special
Definition: declarations_special.F:3
vc08cc
subroutine vc08cc(XMUL, SF, SU, SV, F, U, V, XEL, YEL, IKLE1, IKLE2, IKLE3, IKLE4, IKLE5, IKLE6, NELEM, NELMAX, W1, W2, W3, W4, W5, W6, FORMUL)
Definition: vc08cc.f:9
declarations_special::lu
integer lu
Definition: declarations_special.F:45
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:52 for The TELEMAC-MASCARET system by doxygen 1.8.13