The TELEMAC-MASCARET system  trunk
vc08aa.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE vc08aa
3 ! *****************
4 !
5  &(xmul,sf,su,sv,f,u,v,xel,yel,ikle,
6  & nelem,nelmax,w1,w2,w3,formul)
7 !
8 !***********************************************************************
9 ! BIEF V6P1 21/08/2010
10 !***********************************************************************
11 !
12 !brief COMPUTES THE FOLLOWING VECTOR IN FINITE ELEMENTS:
13 !code
14 !+ / DF DF
15 !+ V = XMUL / PSII * ( U -- + V -- ) D(OMEGA)
16 !+ I /OMEGA DX DY
17 !+
18 !+ PSI(I) IS A BASE OF TYPE P1 TRIANGLE
19 !
20 !warning THE JACOBIAN MUST BE POSITIVE
21 !warning THE RESULT IS IN W IN NOT ASSEMBLED FORM
22 !
23 !history A FROEHLY (MATMECA)
24 !+ 01/07/08
25 !+ V5P9
26 !+
27 !
28 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
29 !+ 13/07/2010
30 !+ V6P0
31 !+ Translation of French comments within the FORTRAN sources into
32 !+ English comments
33 !
34 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
35 !+ 21/08/2010
36 !+ V6P0
37 !+ Creation of DOXYGEN tags for automated documentation and
38 !+ cross-referencing of the FORTRAN sources
39 !
40 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
41 !| F |-->| FUNCTION USED IN THE VECTOR FORMULA
42 !| FORMUL |-->| STRING WITH FORMULA OF VECTOR
43 !| IKLE |-->| CONNECTIVITY TABLE
44 !| NELEM |-->| NUMBER OF ELEMENTS
45 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
46 !| SF |-->| BIEF_OBJ STRUCTURE OF F
47 !| SU |-->| BIEF_OBJ STRUCTURE OF U
48 !| SV |-->| BIEF_OBJ STRUCTURE OF V
49 !| U |-->| FUNCTION USED IN THE VECTOR FORMULA
50 !| V |-->| FUNCTION USED IN THE VECTOR FORMULA
51 !| W1 |<--| RESULT IN NON ASSEMBLED FORM
52 !| W2 |<--| RESULT IN NON ASSEMBLED FORM
53 !| W3 |<--| RESULT IN NON ASSEMBLED FORM
54 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
55 !| XMUL |-->| MULTIPLICATION COEFFICIENT
56 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
57 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
58 !
59  USE bief, ex_vc08aa => vc08aa
60 !
61  USE declarations_special
62  IMPLICIT NONE
63 !
64 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
65 !
66  INTEGER, INTENT(IN) :: NELEM,NELMAX
67  INTEGER, INTENT(IN) :: IKLE(nelmax,*)
68 !
69  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,*),YEL(nelmax,*)
70  DOUBLE PRECISION, INTENT(INOUT) ::W1(nelmax),W2(nelmax),W3(nelmax)
71  DOUBLE PRECISION, INTENT(IN) :: XMUL
72 !
73 ! STRUCTURES OF F, U, V AND REAL DATA
74 !
75  TYPE(bief_obj), INTENT(IN) :: SF,SU,SV
76  DOUBLE PRECISION, INTENT(IN) :: F(*),U(*),V(*)
77 !
78  CHARACTER(LEN=16), INTENT(IN) :: FORMUL
79 !
80 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
81 !
82  INTEGER IELEM,IELMF,IELMU,IELMV
83 !
84  INTRINSIC min,max
85 !
86 !-----------------------------------------------------------------------
87 !
88  DOUBLE PRECISION X2,Y2,X3,Y3,F1,F2,F3
89  DOUBLE PRECISION U1,U2,U3,U4,U5,U6,V1,V2,V3,V4,V5,V6
90  DOUBLE PRECISION SUR6,XSUR24,XSUR120,XSU216,F1MF3,F2MF1
91  DOUBLE PRECISION K1,K2,K3,USUR2,VSUR2,PHIT
92  DOUBLE PRECISION L12,L13,L21,L23,L31,L32,BETAN1,BETAN2,BETAN3
93 !
94 !-----------------------------------------------------------------------
95 !
96  sur6 = 1.d0 / 6.d0
97  xsur24 = xmul/24.d0
98  xsur120= xmul/120.d0
99  xsu216 = xmul/216.d0
100 !
101  ielmf=sf%ELM
102  ielmu=su%ELM
103  ielmv=sv%ELM
104 !
105 !-----------------------------------------------------------------------
106 !
107 ! FUNCTION F AND VECTOR U ARE LINEAR
108 !
109  IF(ielmf.EQ.11.AND.ielmu.EQ.11.AND.ielmv.EQ.11) THEN
110 !
111  IF(formul(14:16).EQ.'PSI') THEN
112 !
113 ! PSI SCHEME
114 !
115  DO ielem = 1 , nelem
116 !
117  x2 = xel(ielem,2)
118  x3 = xel(ielem,3)
119  y2 = yel(ielem,2)
120  y3 = yel(ielem,3)
121 !
122  f1 = f(ikle(ielem,1))
123  f2 = f(ikle(ielem,2))
124  f3 = f(ikle(ielem,3))
125 !
126  u1 = u(ikle(ielem,1))
127  u2 = u(ikle(ielem,2))
128  u3 = u(ikle(ielem,3))
129  v1 = v(ikle(ielem,1))
130  v2 = v(ikle(ielem,2))
131  v3 = v(ikle(ielem,3))
132 !
133  usur2 = (u1+u2+u3)*sur6
134  vsur2 = (v1+v2+v3)*sur6
135 !
136  k1 = usur2 * (y2-y3) - vsur2 * (x2-x3)
137  k2 = usur2 * (y3 ) - vsur2 * (x3 )
138  k3 = usur2 * ( -y2) - vsur2 * ( -x2)
139 !
140  l12 = max( min(k1,-k2) , 0.d0 )
141  l13 = max( min(k1,-k3) , 0.d0 )
142  l21 = max( min(k2,-k1) , 0.d0 )
143  l23 = max( min(k2,-k3) , 0.d0 )
144  l31 = max( min(k3,-k1) , 0.d0 )
145  l32 = max( min(k3,-k2) , 0.d0 )
146 !
147  betan1 = l12*(f1-f2) + l13*(f1-f3)
148  betan2 = l21*(f2-f1) + l23*(f2-f3)
149  betan3 = l31*(f3-f1) + l32*(f3-f2)
150 !
151  phit = betan1 + betan2 + betan3
152 !
153  IF(phit.GT.0.d0) THEN
154  w1(ielem) = xmul * max( min( betan1, phit),0.d0 )
155  w2(ielem) = xmul * max( min( betan2, phit),0.d0 )
156  w3(ielem) = xmul * max( min( betan3, phit),0.d0 )
157  ELSE
158  w1(ielem) = - xmul * max( min(-betan1,-phit),0.d0 )
159  w2(ielem) = - xmul * max( min(-betan2,-phit),0.d0 )
160  w3(ielem) = - xmul * max( min(-betan3,-phit),0.d0 )
161  ENDIF
162 !
163  ENDDO ! IELEM
164 !
165  ELSE
166 !
167 ! NORMAL CENTERED SCHEME
168 !
169  DO ielem = 1 , nelem
170 !
171  x2 = xel(ielem,2)
172  x3 = xel(ielem,3)
173  y2 = yel(ielem,2)
174  y3 = yel(ielem,3)
175 !
176  f1mf3 = f(ikle(ielem,1)) - f(ikle(ielem,3))
177  f2mf1 = f(ikle(ielem,2)) - f(ikle(ielem,1))
178 !
179  u1 = u(ikle(ielem,1))
180  u2 = u(ikle(ielem,2))
181  u3 = u(ikle(ielem,3))
182  v1 = v(ikle(ielem,1))
183  v2 = v(ikle(ielem,2))
184  v3 = v(ikle(ielem,3))
185 !
186  w1(ielem)=( ( y2*f1mf3 + y3*f2mf1 ) * (u1+u1+u2+u3)
187  & - ( x2*f1mf3 + x3*f2mf1 ) * (v1+v1+v2+v3) )*xsur24
188 !
189  w2(ielem)=( ( y2*f1mf3 + y3*f2mf1 ) * (u1+u2+u2+u3)
190  & - ( x2*f1mf3 + x3*f2mf1 ) * (v1+v2+v2+v3) )*xsur24
191 !
192  w3(ielem)=( ( y2*f1mf3 + y3*f2mf1 ) * (u1+u2+u3+u3)
193  & - ( x2*f1mf3 + x3*f2mf1 ) * (v1+v2+v3+v3) )*xsur24
194 !
195  ENDDO ! IELEM
196 !
197  ENDIF
198 !
199 !-----------------------------------------------------------------------
200 !
201 ! FUNCTION F IS LINEAR AND VECTOR U QUASI-BUBBLE
202 !
203  ELSEIF(ielmf.EQ.11.AND.ielmu.EQ.12.AND.ielmv.EQ.12) THEN
204 !
205  IF(formul(14:16).EQ.'PSI') THEN
206 !
207 ! PSI SCHEME (U4 DISCARDED HERE)
208 ! (AS IF U WAS P1 )
209 !
210  DO ielem = 1 , nelem
211 !
212  x2 = xel(ielem,2)
213  x3 = xel(ielem,3)
214  y2 = yel(ielem,2)
215  y3 = yel(ielem,3)
216 !
217  f1 = f(ikle(ielem,1))
218  f2 = f(ikle(ielem,2))
219  f3 = f(ikle(ielem,3))
220 !
221  u1 = u(ikle(ielem,1))
222  u2 = u(ikle(ielem,2))
223  u3 = u(ikle(ielem,3))
224  v1 = v(ikle(ielem,1))
225  v2 = v(ikle(ielem,2))
226  v3 = v(ikle(ielem,3))
227 !
228  usur2 = (u1+u2+u3)*sur6
229  vsur2 = (v1+v2+v3)*sur6
230 !
231  k1 = usur2 * (y2-y3) - vsur2 * (x2-x3)
232  k2 = usur2 * (y3 ) - vsur2 * (x3 )
233  k3 = usur2 * ( -y2) - vsur2 * ( -x2)
234 !
235  l12 = max( min(k1,-k2) , 0.d0 )
236  l13 = max( min(k1,-k3) , 0.d0 )
237  l21 = max( min(k2,-k1) , 0.d0 )
238  l23 = max( min(k2,-k3) , 0.d0 )
239  l31 = max( min(k3,-k1) , 0.d0 )
240  l32 = max( min(k3,-k2) , 0.d0 )
241 !
242  betan1 = l12*(f1-f2) + l13*(f1-f3)
243  betan2 = l21*(f2-f1) + l23*(f2-f3)
244  betan3 = l31*(f3-f1) + l32*(f3-f2)
245 !
246  phit = betan1 + betan2 + betan3
247 !
248  IF(phit.GT.0.d0) THEN
249  w1(ielem) = xmul * max( min( betan1, phit),0.d0 )
250  w2(ielem) = xmul * max( min( betan2, phit),0.d0 )
251  w3(ielem) = xmul * max( min( betan3, phit),0.d0 )
252  ELSE
253  w1(ielem) = - xmul * max( min(-betan1,-phit),0.d0 )
254  w2(ielem) = - xmul * max( min(-betan2,-phit),0.d0 )
255  w3(ielem) = - xmul * max( min(-betan3,-phit),0.d0 )
256  ENDIF
257 !
258  ENDDO ! IELEM
259 !
260  ELSE
261 !
262 ! NORMAL CENTERED SCHEME
263 !
264  DO ielem = 1 , nelem
265 !
266  x2 = xel(ielem,2)
267  x3 = xel(ielem,3)
268  y2 = yel(ielem,2)
269  y3 = yel(ielem,3)
270 !
271  f1 = f(ikle(ielem,1))
272  f2 = f(ikle(ielem,2)) - f1
273  f3 = f(ikle(ielem,3)) - f1
274 !
275  u1 = u(ikle(ielem,1))
276  u2 = u(ikle(ielem,2))
277  u3 = u(ikle(ielem,3))
278  u4 = u(ikle(ielem,4))
279  v1 = v(ikle(ielem,1))
280  v2 = v(ikle(ielem,2))
281  v3 = v(ikle(ielem,3))
282  v4 = v(ikle(ielem,4))
283 !
284  w1(ielem)=(5*x2*f3*v3+12*x2*f3*v4+5*x2*f3*v2+14*x2*f3*v1-5
285  & *x3*f2*v3-12*x3*f2*v4-5*x3*f2*v2-14*x3*f2*v1-5*f3*u3*
286  & y2-12*f3*u4*y2-5*f3*u2*y2-14*f3*u1*y2+5*f2*u3*y3+12*
287  & f2*u4*y3+5*f2*u2*y3+14*f2*u1*y3)*xsu216
288  w2(ielem)=(5*x2*f3*v3+12*x2*f3*v4+14*x2*f3*v2+5*x2*f3*v1-5
289  & *x3*f2*v3-12*x3*f2*v4-14*x3*f2*v2-5*x3*f2*v1-5*f3*u3*
290  & y2-12*f3*u4*y2-14*f3*u2*y2-5*f3*u1*y2+5*f2*u3*y3+12*
291  & f2*u4*y3+14*f2*u2*y3+5*f2*u1*y3)*xsu216
292  w3(ielem)=(14*x2*f3*v3+12*x2*f3*v4+5*x2*f3*v2+5*x2*f3*v1-
293  & 14*x3*f2*v3-12*x3*f2*v4-5*x3*f2*v2-5*x3*f2*v1-14*f3*
294  & u3*y2-12*f3*u4*y2-5*f3*u2*y2-5*f3*u1*y2+14*f2*u3*y3+
295  & 12*f2*u4*y3+5*f2*u2*y3+5*f2*u1*y3)*xsu216
296 !
297  ENDDO ! IELEM
298 !
299  ENDIF
300 !
301 !-----------------------------------------------------------------------
302 !
303 ! FUNCTION F IS LINEAR AND VECTOR U P2
304 !
305  ELSEIF(ielmf.EQ.11.AND.ielmu.EQ.13.AND.ielmv.EQ.13) THEN
306 !
307  IF(formul(14:16).EQ.'PSI') THEN
308 !
309 ! PSI SCHEME (U4,U5 AND U6 DISCARDED HERE)
310 ! (AS IF U WAS P1 )
311 !
312  DO ielem = 1 , nelem
313 !
314  x2 = xel(ielem,2)
315  x3 = xel(ielem,3)
316  y2 = yel(ielem,2)
317  y3 = yel(ielem,3)
318 !
319  f1 = f(ikle(ielem,1))
320  f2 = f(ikle(ielem,2))
321  f3 = f(ikle(ielem,3))
322 !
323  u1 = u(ikle(ielem,1))
324  u2 = u(ikle(ielem,2))
325  u3 = u(ikle(ielem,3))
326  v1 = v(ikle(ielem,1))
327  v2 = v(ikle(ielem,2))
328  v3 = v(ikle(ielem,3))
329 !
330  usur2 = (u1+u2+u3)*sur6
331  vsur2 = (v1+v2+v3)*sur6
332 !
333  k1 = usur2 * (y2-y3) - vsur2 * (x2-x3)
334  k2 = usur2 * (y3 ) - vsur2 * (x3 )
335  k3 = usur2 * ( -y2) - vsur2 * ( -x2)
336 !
337  l12 = max( min(k1,-k2) , 0.d0 )
338  l13 = max( min(k1,-k3) , 0.d0 )
339  l21 = max( min(k2,-k1) , 0.d0 )
340  l23 = max( min(k2,-k3) , 0.d0 )
341  l31 = max( min(k3,-k1) , 0.d0 )
342  l32 = max( min(k3,-k2) , 0.d0 )
343 !
344  betan1 = l12*(f1-f2) + l13*(f1-f3)
345  betan2 = l21*(f2-f1) + l23*(f2-f3)
346  betan3 = l31*(f3-f1) + l32*(f3-f2)
347 !
348  phit = betan1 + betan2 + betan3
349 !
350  IF(phit.GT.0.d0) THEN
351  w1(ielem) = xmul * max( min( betan1, phit),0.d0 )
352  w2(ielem) = xmul * max( min( betan2, phit),0.d0 )
353  w3(ielem) = xmul * max( min( betan3, phit),0.d0 )
354  ELSE
355  w1(ielem) = - xmul * max( min(-betan1,-phit),0.d0 )
356  w2(ielem) = - xmul * max( min(-betan2,-phit),0.d0 )
357  w3(ielem) = - xmul * max( min(-betan3,-phit),0.d0 )
358  ENDIF
359 !
360  ENDDO ! IELEM
361 !
362  ELSE
363 !
364 ! NORMAL CENTERED SCHEME
365 !
366  DO ielem = 1 , nelem
367 !
368  x2 = xel(ielem,2)
369  x3 = xel(ielem,3)
370  y2 = yel(ielem,2)
371  y3 = yel(ielem,3)
372 !
373  f1 = f(ikle(ielem,1))
374  f2 = f(ikle(ielem,2)) - f1
375  f3 = f(ikle(ielem,3)) - f1
376 !
377  u1 = u(ikle(ielem,1))
378  u2 = u(ikle(ielem,2))
379  u3 = u(ikle(ielem,3))
380  u4 = u(ikle(ielem,4))
381  u5 = u(ikle(ielem,5))
382  u6 = u(ikle(ielem,6))
383  v1 = v(ikle(ielem,1))
384  v2 = v(ikle(ielem,2))
385  v3 = v(ikle(ielem,3))
386  v4 = v(ikle(ielem,4))
387  v5 = v(ikle(ielem,5))
388  v6 = v(ikle(ielem,6))
389 !
390  w1(ielem)=
391  & (2.d0*u1*y3*f2-2.d0*u1*y2*f3-v2*x2*f3-8.d0*v4*x3*f2+
392  & 8.d0*u4*y3*f2+v3*x3*f2-u2*y3*f2-v3*x2*f3+4.d0*u5*y3*f2-
393  & 4.d0*u5*y2*f3+u3*y2*f3-4.d0*v5*x3*f2+v2*x3*f2-
394  & 8.d0*u6*y2*f3-8.d0*u4*y2*f3+4.d0*v5*x2*f3+
395  & 2.d0*v1*x2*f3-u3*y3*f2+8.d0*v6*x2*f3-8.d0*v6*x3*f2+
396  & 8.d0*v4*x2*f3-2.d0*v1*x3*f2+8.d0*u6*y3*f2+
397  & u2*y2*f3)*xsur120
398 !
399  w2(ielem)=
400  & -(-8.d0*v5*x2*f3-4.d0*u6*y3*f2-v1*x3*f2+4.d0*v6*x3*f2+
401  & 8.d0*u4*y2*f3-4.d0*v6*x2*f3+2.d0*u2*y2*f3+
402  & 8.d0*u5*y2*f3+v3*x2*f3+8.d0*v5*x3*f2+2.d0*v2*x3*f2-
403  & 2.d0*u2*y3*f2-8.d0*u5*y3*f2+v1*x2*f3-2.d0*v2*x2*f3-
404  & v3*x3*f2-8.d0*v4*x2*f3+u3*y3*f2+u1*y3*f2-u1*y2*f3+
405  & 4.d0*u6*y2*f3-8.d0*u4*y3*f2+8.d0*v4*x3*f2-u3*y2*f3)
406  & *xsur120
407 !
408  w3(ielem) =
409  & (-v5*x3*f2*8.d0+v2*x3*f2-u6*y2*f3*8.d0-u4*y2*f3*4.d0+
410  & v5*x2*f3*8.d0-v1*x2*f3+u3*y3*f2*2.d0+v6*x2*f3*8.d0-
411  & v6*x3*f2*8.d0+v4*x2*f3*4.d0+v1*x3*f2+u6*y3*f2*8.d0+
412  & u2*y2*f3-u1*y3*f2+u1*y2*f3-v2*x2*f3-v4*x3*f2*4.d0+
413  & u4*y3*f2*4.d0-v3*x3*f2*2.d0-u2*y3*f2+v3*x2*f3*2.d0+
414  & u5*y3*f2*8.d0-u5*y2*f3*8.d0-u3*y2*f3*2.d0)*xsur120
415 !
416  ENDDO ! IELEM
417 !
418  ENDIF
419 !
420 !-----------------------------------------------------------------------
421 !
422  ELSE
423 !
424 !-----------------------------------------------------------------------
425 !
426  WRITE(lu,101) ielmf,sf%NAME
427  WRITE(lu,201) ielmu,su%NAME
428  WRITE(lu,301)
429 101 FORMAT(1x,'VC08AA (BIEF) :',/,
430  & 1x,'DISCRETIZATION OF F:',1i6,
431  & 1x,'REAL NAME: ',a6)
432 201 FORMAT(1x,'DISCRETIZATION OF U:',1i6,
433  & 1x,'REAL NAME: ',a6)
434 301 FORMAT(1x,'CASE NOT IMPLEMENTED')
435  CALL plante(1)
436  stop
437 !
438  ENDIF
439 !
440 !-----------------------------------------------------------------------
441 !
442  RETURN
443  END
declarations_special
Definition: declarations_special.F:3
vc08aa
subroutine vc08aa(XMUL, SF, SU, SV, F, U, V, XEL, YEL, IKLE, NELEM, NELMAX, W1, W2, W3, FORMUL)
Definition: vc08aa.f:8
declarations_special::lu
integer lu
Definition: declarations_special.F:45
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:51 for The TELEMAC-MASCARET system by doxygen 1.8.13