The TELEMAC-MASCARET system  trunk
vc08bb.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE vc08bb
3 ! *****************
4 !
5  &(xmul,sf,su,sv,f,u,v,xel,yel,
6  & ikle1,ikle2,ikle3,ikle4,nelem,nelmax,w1,w2,w3,w4,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 QUASI-BUBBLE TYPE
19 !
20 !warning THE JACOBIAN MUST BE POSITIVE
21 !warning THE RESULT IS IN W IN NOT ASSEMBLED FORM
22 !
23 !history J-M HERVOUET (LNH) ; F LEPEINTRE (LNH)
24 !+ 29/12/05
25 !+ V5P6
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 !history S.E.BOURBAN (HRW)
41 !+ 21/03/2017
42 !+ V7P3
43 !+ Replacement of the DATA declarations by the PARAMETER associates
44 !
45 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
46 !| F |-->| FUNCTION USED IN THE VECTOR FORMULA
47 !| FORMUL |-->| STRING WITH FORMULA OF VECTOR
48 !| IKLE1 |-->| FIRST POINT OF TRIANGLES
49 !| IKLE2 |-->| SECOND POINT OF TRIANGLES
50 !| IKLE3 |-->| THIRD POINT OF TRIANGLES
51 !| IKLE4 |-->| QUASI-BUBBLE POINT OF TRIANGLES
52 !| NELEM |-->| NUMBER OF ELEMENTS
53 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
54 !| SF |-->| BIEF_OBJ STRUCTURE OF F
55 !| SU |-->| BIEF_OBJ STRUCTURE OF U
56 !| SV |-->| BIEF_OBJ STRUCTURE OF V
57 !| U |-->| FUNCTION USED IN THE VECTOR FORMULA
58 !| V |-->| FUNCTION USED IN THE VECTOR FORMULA
59 !| W1 |<--| RESULT IN NON ASSEMBLED FORM
60 !| W2 |<--| RESULT IN NON ASSEMBLED FORM
61 !| W3 |<--| RESULT IN NON ASSEMBLED FORM
62 !| W4 |<--| 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_vc08bb => vc08bb
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(4*nelmax),W2(nelmax)
79  DOUBLE PRECISION, INTENT(INOUT) :: W3(nelmax),W4(nelmax)
80 ! IKLE1 IS ALSO USED AS A 1-DIMENSIONAL IKLE
81  INTEGER, INTENT(IN) :: IKLE1(4*nelmax)
82  INTEGER, INTENT(IN) :: IKLE2(nelmax),IKLE3(nelmax),IKLE4(nelmax)
83  CHARACTER(LEN=16), INTENT(IN) :: FORMUL
84 !
85 ! STRUCTURES OF F, G, H, U, V, W AND REAL DATA
86 !
87  TYPE(bief_obj), INTENT(IN) :: SF,SU,SV
88  DOUBLE PRECISION, INTENT(IN) :: F(*),U(*),V(*)
89 !
90 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
91 !
92  INTEGER IELEM,IELMF,IELMU,IELMV
93  INTEGER IG1,IG2,IG3,IT,IAD1,IAD2,IAD3
94 !
95  DOUBLE PRECISION K1,K2,K3
96 !
97 !-----------------------------------------------------------------------
98 !
99  DOUBLE PRECISION X2,Y2,X3,Y3,F1,F2,F3,F4,U1,U2,U3,U4,V1,V2,V3,V4
100  DOUBLE PRECISION XSUR72,XSU216,X1,Y1,TIERS
101  DOUBLE PRECISION PHIT,SUR6,USUR2,VSUR2
102  DOUBLE PRECISION L12,L13,L21,L23,L31,L32,SIG,BETAN1,BETAN2,BETAN3
103 !
104  INTRINSIC max,min,sign
105 !
106 !-----------------------------------------------------------------------
107 !
108 ! FOR A QUASI-BUBBLE TRIANGLE : NODE NUMBERS OF THE SUB-TRIANGLES
109 ! IN THE INITIAL TRIANGLE
110 ! IL(NUMBER OF THE SUB-TRIANGLE,LOCAL NUMBER IN THE SUB-TRIANGLE)
111 !
112  INTEGER :: IL(3,3)
113  parameter( il = reshape( (/
114  & 1,2,3,2,3,1,4,4,4 /), shape=(/ 3,3 /) ) )
115 !
116 !-----------------------------------------------------------------------
117 !
118  xsur72 = xmul/72.d0
119  xsu216 = xmul/216.d0
120  tiers = 1.d0 / 3.d0
121  sur6 = 1.d0 / 6.d0
122 !
123  ielmf=sf%ELM
124  ielmu=su%ELM
125  ielmv=sv%ELM
126 !
127 !-----------------------------------------------------------------------
128 !
129 ! FUNCTION F AND VECTOR U ARE QUASI-BUBBLE
130 !
131  IF(ielmf.EQ.12.AND.ielmu.EQ.12.AND.ielmv.EQ.12) THEN
132 !
133  IF(formul(14:16).EQ.'PSI') THEN
134 !
135 ! INITIALISES W
136 !
137  DO ielem = 1 , nelem
138  w1(ielem) = 0.d0
139  w2(ielem) = 0.d0
140  w3(ielem) = 0.d0
141  w4(ielem) = 0.d0
142  ENDDO ! IELEM
143 !
144 ! PSI SCHEME, LOOP ON THE 3 SUB-TRIANGLES AND
145 ! PREASSEMBLY
146 !
147  DO it=1,3
148  DO ielem = 1 , nelem
149 !
150 ! ADDRESSES IN ARRAY (NELMAX,*)
151  iad1= ielem + (il(it,1)-1)*nelmax
152  iad2= ielem + (il(it,2)-1)*nelmax
153  iad3= ielem + (il(it,3)-1)*nelmax
154 ! GLOBAL NUMBERS IN THE INITIAL TRIANGLE
155  ig1 = ikle1(iad1)
156  ig2 = ikle1(iad2)
157  ig3 = ikle1(iad3)
158 ! COORDINATES OF THE SUB-TRIANGLE NODES
159  x1 = xel(iad1)
160  x2 = xel(iad2) - x1
161 ! POINT 3 IS ALWAYS THE CENTRE OF THE INITIAL TRIANGLE
162  x3=tiers*(xel(ielem)+xel(ielem+nelmax)+xel(ielem+2*nelmax))-x1
163  y1 = yel(iad1)
164  y2 = yel(iad2) - y1
165 ! POINT 3 IS ALWAYS THE CENTRE OF THE INITIAL TRIANGLE
166  y3=tiers*(yel(ielem)+yel(ielem+nelmax)+yel(ielem+2*nelmax))-y1
167 ! F VALUES IN THE SUB-TRIANGLE
168  f1 = f(ig1)
169  f2 = f(ig2)
170  f3 = f(ig3)
171 !
172  usur2 = (u(ig1)+u(ig2)+u(ig3))*sur6
173  vsur2 = (v(ig1)+v(ig2)+v(ig3))*sur6
174 !
175  k1 = usur2 * (y2-y3) - vsur2 * (x2-x3)
176  k2 = usur2 * (y3 ) - vsur2 * (x3 )
177  k3 = usur2 * ( -y2) - vsur2 * ( -x2)
178 !
179  l12 = max( min(k1,-k2) , 0.d0 )
180  l13 = max( min(k1,-k3) , 0.d0 )
181  l21 = max( min(k2,-k1) , 0.d0 )
182  l23 = max( min(k2,-k3) , 0.d0 )
183  l31 = max( min(k3,-k1) , 0.d0 )
184  l32 = max( min(k3,-k2) , 0.d0 )
185 !
186  betan1 = l12*(f1-f2) + l13*(f1-f3)
187  betan2 = l21*(f2-f1) + l23*(f2-f3)
188  betan3 = l31*(f3-f1) + l32*(f3-f2)
189 !
190  phit = betan1 + betan2 + betan3
191  sig = sign(1.d0,phit)
192 !
193  w1(iad1)=w1(iad1)+ xmul * sig * max(min(sig*betan1,sig*phit),0.d0)
194  w1(iad2)=w1(iad2)+ xmul * sig * max(min(sig*betan2,sig*phit),0.d0)
195  w1(iad3)=w1(iad3)+ xmul * sig * max(min(sig*betan3,sig*phit),0.d0)
196 !
197  ENDDO ! IELEM
198  ENDDO ! IT
199 !
200  ELSE
201 !
202 ! CLASSICAL COMPUTATION
203 !
204  DO ielem = 1 , nelem
205 !
206  x2 = xel(ielem+nelmax)
207  x3 = xel(ielem+2*nelmax)
208  y2 = yel(ielem+nelmax)
209  y3 = yel(ielem+2*nelmax)
210 !
211  u1 = u(ikle1(ielem))
212  u2 = u(ikle2(ielem))
213  u3 = u(ikle3(ielem))
214  u4 = u(ikle4(ielem))
215  v1 = v(ikle1(ielem))
216  v2 = v(ikle2(ielem))
217  v3 = v(ikle3(ielem))
218  v4 = v(ikle4(ielem))
219 !
220  f1 = f(ikle1(ielem))
221  f2 = f(ikle2(ielem)) - f1
222  f3 = f(ikle3(ielem)) - f1
223  f4 = f(ikle4(ielem)) - f1
224 !
225  w1(ielem)= (((v3+v4+2*v1)*x2+(v3+v4+2*v1)*x3-(y3+y2)*u3-(y3+
226  & y2)*u4-2*(y3+y2)*u1)*f3-3*((v3+v4+2*v1)*x3-(v4+v2+2*
227  & v1)*x2-(y3-y2)*u4-2*(y3-y2)*u1-u3*y3+u2*y2)*f4-((v4+v2+
228  & 2*v1)*x2+(v4+v2+2*v1)*x3-(y3+y2)*u4-(y3+y2)*u2-2*(y3+y2
229  & )*u1)*f2)*xsur72
230 !
231  w2(ielem)= (-(((2*v3+3*v4+6*v2+v1)*x3-(v3-v1)*x2-(2*y3-y2)
232  & *u3-(y3+y2)*u1-3*u4*y3-6*u2*y3)*f2-(2*(v3+v4+2*v2)*x2
233  & -(v3+v4+2*v2)*x3+(y3-2*y2)*u3+(y3-2*y2)*u4+2*(y3-2*
234  & y2)*u2)*f3-3*((v3+v4+2*v2)*x3-(v3-v1)*x2-(y3-y2)*u3-u4*
235  & y3-2*u2*y3-u1*y2)*f4))*xsur72
236 !
237  w3(ielem)= (((6*v3+3*v4+2*v2+v1)*x2-(v2-v1)*x3-(y3+y2)*u1+(
238  & y3-2*y2)*u2-6*u3*y2-3*u4*y2)*f3+((2*v3+v4+v2)*x2-2*(
239  & 2*v3+v4+v2)*x3+2*(2*y3-y2)*u3+(2*y3-y2)*u4+(2*y3-y2)*
240  & u2)*f2-3*((2*v3+v4+v2)*x2-(v2-v1)*x3+(y3-y2)*u2-2*u3*
241  & y2-u4*y2-u1*y3)*f4)*xsur72
242 !
243  w4(ielem)= (((3*v3+6*v4+2*v2+v1)*x2-(v2-v1)*x3-(y3+y2)*u1+(
244  & y3-2*y2)*u2-3*u3*y2-6*u4*y2)*f3-((2*v3+6*v4+3*v2+v1
245  & )*x3-(v3-v1)*x2-(2*y3-y2)*u3-(y3+y2)*u1-6*u4*y3-3*u2*
246  & y3)*f2-3*((v3-v1)*x2-(v2-v1)*x3-(y3-y2)*u1-u3*y2+u2*y3)*
247  & f4)*xsur72
248 !
249  ENDDO ! IELEM
250 !
251  ENDIF
252 !
253 ! FUNCTION F IS QUASI-BUBBLE AND VECTOR U LINEAR
254 !
255  ELSEIF(ielmf.EQ.12.AND.ielmu.EQ.11.AND.ielmv.EQ.11) THEN
256 !
257  IF(formul(14:16).EQ.'PSI') THEN
258 !
259  WRITE(lu,401)
260 401 FORMAT(1x,'VC08BB (BIEF) : PSI NOT IMPLEMENTED FOR QUASI-BUBBLE')
261  CALL plante(1)
262  stop
263 !
264  ELSE
265 !
266 ! CLASSICAL COMPUTATION
267 !
268  DO ielem = 1 , nelem
269 !
270  x2 = xel(ielem+nelmax)
271  x3 = xel(ielem+2*nelmax)
272  y2 = yel(ielem+nelmax)
273  y3 = yel(ielem+2*nelmax)
274 !
275  u1 = u(ikle1(ielem))
276  u2 = u(ikle2(ielem))
277  u3 = u(ikle3(ielem))
278  v1 = v(ikle1(ielem))
279  v2 = v(ikle2(ielem))
280  v3 = v(ikle3(ielem))
281 !
282  f1 = f(ikle1(ielem))
283  f2 = f(ikle2(ielem)) - f1
284  f3 = f(ikle3(ielem)) - f1
285  f4 = f(ikle4(ielem)) - f1
286 !
287  w1(ielem)=(((4*v3+v2+7*v1)*x2+(4*v3+v2+7*v1)*x3-4*(y3+y2
288  & )*u3-(y3+y2)*u2-7*(y3+y2)*u1)*f3-3*((4*v3+v2+7*v1)*x3
289  & -(v3+4*v2+7*v1)*x2-(4*y3-y2)*u3-7*(y3-y2)*u1-(y3-4*
290  & y2)*u2)*f4-((v3+4*v2+7*v1)*x2+(v3+4*v2+7*v1)*x3-(y3+
291  & y2)*u3-4*(y3+y2)*u2-7*(y3+y2)*u1)*f2)*xsu216
292  w2(ielem)=((2*(4*v3+7*v2+v1)*x2-(4*v3+7*v2+v1)*x3+4*(y3
293  & -2*y2)*u3+7*(y3-2*y2)*u2+(y3-2*y2)*u1)*f3+3*((4*v3+
294  & 7*v2+v1)*x3-3*(v3-v1)*x2-(4*y3-3*y2)*u3-(y3+3*y2)*u1-
295  & 7*u2*y3)*f4-3*((3*v3+7*v2+2*v1)*x3-(v3-v1)*x2-(3*y3-
296  & y2)*u3-(2*y3+y2)*u1-7*u2*y3)*f2)*xsu216
297  w3(ielem)=(((7*v3+4*v2+v1)*x2-2*(7*v3+4*v2+v1)*x3+7*(2
298  & *y3-y2)*u3+4*(2*y3-y2)*u2+(2*y3-y2)*u1)*f2-3*((7*v3+
299  & 4*v2+v1)*x2-3*(v2-v1)*x3-(3*y3+y2)*u1+(3*y3-4*y2)*u2-
300  & 7*u3*y2)*f4+3*((7*v3+3*v2+2*v1)*x2-(v2-v1)*x3-(y3+2*
301  & y2)*u1+(y3-3*y2)*u2-7*u3*y2)*f3)*xsu216
302  w4(ielem)=(((5*v3+4*v2+3*v1)*x2-(v2-v1)*x3-(y3+3*y2)*u1+(
303  & y3-4*y2)*u2-5*u3*y2)*f3-((4*v3+5*v2+3*v1)*x3-(v3-v1)
304  & *x2-(4*y3-y2)*u3-(3*y3+y2)*u1-5*u2*y3)*f2-3*((v3-v1)*
305  & x2-(v2-v1)*x3-(y3-y2)*u1-u3*y2+u2*y3)*f4)*xsur72
306 !
307  ENDDO ! IELEM
308 !
309  ENDIF
310 !
311 !-----------------------------------------------------------------------
312 !
313  ELSE
314 !
315 !-----------------------------------------------------------------------
316 !
317  WRITE(lu,101) ielmf,sf%NAME
318  WRITE(lu,201) ielmu,su%NAME
319  WRITE(lu,301)
320 101 FORMAT(1x,'VC08BB (BIEF) :',/,
321  & 1x,'DISCRETIZATION OF F:',1i6,
322  & 1x,'REAL NAME: ',a6)
323 201 FORMAT(1x,'DISCRETIZATION OF U:',1i6,
324  & 1x,'REAL NAME: ',a6)
325 301 FORMAT(1x,'CASE NOT IMPLEMENTED')
326  CALL plante(1)
327  stop
328 !
329  ENDIF
330 !
331 !-----------------------------------------------------------------------
332 !
333  RETURN
334  END
declarations_special
Definition: declarations_special.F:3
vc08bb
subroutine vc08bb(XMUL, SF, SU, SV, F, U, V, XEL, YEL, IKLE1, IKLE2, IKLE3, IKLE4, NELEM, NELMAX, W1, W2, W3, W4, FORMUL)
Definition: vc08bb.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