The TELEMAC-MASCARET system  trunk
vc19aa.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE vc19aa
3 ! *****************
4 !
5  &(xmul,sf,sg,sh,su,sv,f,g,h,u,v,
6  & xel,yel,ikle1,ikle2,ikle3,nelem,nelmax,w1,w2,w3,formul)
7 !
8 !***********************************************************************
9 ! BIEF V7P0 21/08/2010
10 !***********************************************************************
11 !
12 !brief COMPUTES THE FOLLOWING VECTOR IN FINITE ELEMENTS:
13 !code
14 !+ / DPSII DPSII
15 !+ V = XMUL / F * ( U ----- + V ----- ) D(OMEGA)
16 !+ I /OMEGA DX DY
17 !+
18 !+
19 !+ PSI(I) IS A BASE OF TYPE P1 TRIANGLE
20 !+
21 !+ F, U AND V ARE VECTORS
22 !+ ->
23 !+ BEWARE: IF FORMUL='HUGRADP ' THE VELOCITY IS: U
24 !+ -> --->
25 !+ IF FORMUL='HUGRADP2' THE VELOCITY IS: U + G * GRAD(H)
26 !
27 !warning THE JACOBIAN MUST BE POSITIVE
28 !warning THE RESULT IS IN W IN NOT ASSEMBLED FORM
29 !warning U AND V ARE QUASI-BUBBLE; TREATED AS IF LINEAR!!!!!!!!!!!!
30 !
31 !history C-T PHAM (LNHE)
32 !+ 09/01/08
33 !+ V5P8
34 !+
35 !
36 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
37 !+ 13/07/2010
38 !+ V6P0
39 !+ Translation of French comments within the FORTRAN sources into
40 !+ English comments
41 !
42 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
43 !+ 21/08/2010
44 !+ V6P0
45 !+ Creation of DOXYGEN tags for automated documentation and
46 !+ cross-referencing of the FORTRAN sources
47 !
48 !history J-M HERVOUET (EDF LAB, LNHE)
49 !+ 13/05/2014
50 !+ V7P0
51 !+ Discontinuous elements better treated: new types 15, 16 and 17 for
52 !+ discontinuous linear, quasi-bubble, and quadratic, rather than
53 !+ using component DIMDISC=11, 12 or 13.
54 !
55 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
56 !| F |-->| FUNCTION USED IN THE VECTOR FORMULA
57 !| FORMUL |-->| STRING WITH THE FORMULA DESCRIBING THE VECTOR
58 !| G |-->| FUNCTION USED IN THE VECTOR FORMULA
59 !| H |-->| FUNCTION USED IN THE VECTOR FORMULA
60 !| IKLE1 |-->| FIRST POINT OF TRIANGLES
61 !| IKLE2 |-->| SECOND POINT OF TRIANGLES
62 !| IKLE3 |-->| THIRD POINT OF TRIANGLES
63 !| NELEM |-->| NUMBER OF ELEMENTS
64 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
65 !| SF |-->| BIEF_OBJ STRUCTURE OF F
66 !| SG |-->| BIEF_OBJ STRUCTURE OF G
67 !| SH |-->| BIEF_OBJ STRUCTURE OF H
68 !| SU |-->| BIEF_OBJ STRUCTURE OF U
69 !| SV |-->| BIEF_OBJ STRUCTURE OF V
70 !| U |-->| FUNCTION USED IN THE VECTOR FORMULA
71 !| V |-->| FUNCTION USED IN THE VECTOR FORMULA
72 !| W1 |<--| RESULT IN NON ASSEMBLED FORM
73 !| W2 |<--| RESULT IN NON ASSEMBLED FORM
74 !| W3 |<--| RESULT IN NON ASSEMBLED FORM
75 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
76 !| XMUL |-->| MULTIPLICATION COEFFICIENT
77 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
78 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
79 !
80  USE bief, ex_vc19aa => vc19aa
81 !
82  USE declarations_special
83  IMPLICIT NONE
84 !
85 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
86 !
87  INTEGER, INTENT(IN) :: NELEM,NELMAX
88  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax),IKLE3(nelmax)
89  DOUBLE PRECISION, INTENT(IN) :: XMUL
90  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,*),YEL(nelmax,*)
91  DOUBLE PRECISION, INTENT(INOUT) :: W1(*),W2(*),W3(*)
92 !
93 ! STRUCTURES OF H, U, V AND REAL DATA
94 !
95  TYPE(bief_obj), INTENT(IN) :: SF,SH,SG,SU,SV
96  DOUBLE PRECISION, INTENT(IN) :: F(*),G(*),H(*),U(*),V(*)
97  CHARACTER(LEN=16), INTENT(IN) :: FORMUL
98 !
99 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
100 !
101  INTEGER IELEM,IELMU,IELMV,IELMF,IELMH,IELMG,I1,I2,I3
102 !
103  DOUBLE PRECISION X2,Y2,X3,Y3,Z1,Z2,Z3,ZX,ZY,DET
104  DOUBLE PRECISION H1,H2,H3,U1,U2,U3,V1,V2,V3
105  DOUBLE PRECISION H123,U123,V123,HU123,HV123
106  DOUBLE PRECISION XSUR24
107 !
108 !-----------------------------------------------------------------------
109 !
110  xsur24 = xmul / 24.d0
111 !
112 !-----------------------------------------------------------------------
113 !
114  ielmu=su%ELM
115  ielmv=sv%ELM
116  ielmf=sf%ELM
117  ielmg=sg%ELM
118  ielmh=sh%ELM
119 !
120 !-----------------------------------------------------------------------
121 !
122  IF(formul(1:8).EQ.'HUGRADP ') THEN
123 !
124 ! F LINEAR, U AND V ARE LINEAR OR QUASI-BUBBLE
125 ! BUT THE FOURTH POINT DISCARDED
126 !
127  IF( ielmf.EQ.11
128  & .AND.(ielmu.EQ.11.OR.ielmu.EQ.12)
129  & .AND.(ielmv.EQ.11.OR.ielmv.EQ.12) ) THEN
130 !
131  DO ielem = 1 , nelem
132 !
133  x2 = xel(ielem,2)
134  x3 = xel(ielem,3)
135  y2 = yel(ielem,2)
136  y3 = yel(ielem,3)
137 !
138  i1 = ikle1(ielem)
139  i2 = ikle2(ielem)
140  i3 = ikle3(ielem)
141  h1 = f(i1)
142  h2 = f(i2)
143  h3 = f(i3)
144  u1 = u(i1)
145  u2 = u(i2)
146  u3 = u(i3)
147  v1 = v(i1)
148  v2 = v(i2)
149  v3 = v(i3)
150 !
151  h123 = h1+h2+h3
152  u123 = u1+u2+u3
153  v123 = v1+v2+v3
154 !
155  hu123 = h1*u1+h2*u2+h3*u3
156  hv123 = h1*v1+h2*v2+h3*v3
157 !
158  w1(ielem) = ( (y2-y3)*(h123*u123+hu123)
159  & +(x3-x2)*(h123*v123+hv123) )*xsur24
160  w2(ielem) = ( y3 *(h123*u123+hu123)
161  & -x3 *(h123*v123+hv123) )*xsur24
162  w3(ielem) = ( -y2 *(h123*u123+hu123)
163  & +x2 *(h123*v123+hv123) )*xsur24
164 !
165  ENDDO
166 !
167 !-----------------------------------------------------------------------
168 !
169  ELSE
170 !
171 !-----------------------------------------------------------------------
172 !
173  WRITE(lu,101) ielmf,sf%NAME
174  WRITE(lu,401) ielmu,su%NAME
175  WRITE(lu,501) ielmv,sv%NAME
176  WRITE(lu,601)
177 101 FORMAT(1x,'VC19AA (BIEF) :',/,
178  & 1x,'DISCRETIZATION OF F:',1i6,
179  & 1x,'REAL NAME: ',a6)
180 201 FORMAT(1x,'DISCRETIZATION OF G:',1i6,
181  & 1x,'REAL NAME: ',a6)
182 301 FORMAT(1x,'DISCRETIZATION OF H:',1i6,
183  & 1x,'REAL NAME: ',a6)
184 401 FORMAT(1x,'DISCRETIZATION OF U:',1i6,
185  & 1x,'REAL NAME: ',a6)
186 501 FORMAT(1x,'DISCRETIZATION OF V:',1i6,
187  & 1x,'REAL NAME: ',a6)
188 601 FORMAT(1x,'CASE NOT IMPLEMENTED')
189  CALL plante(1)
190  stop
191 !
192  ENDIF
193 !
194 !-----------------------------------------------------------------------
195 !
196  ELSEIF(formul(1:8).EQ.'HUGRADP2') THEN
197 !
198 !-----------------------------------------------------------------------
199 !
200 !
201 ! F AND G ARE LINEAR, U, V ARE LINEAR OR QUASI-BUBBLE
202 ! H IS PIECEWISE LINEAR
203 !
204  IF( ielmf.EQ.11
205  & .AND.ielmg.EQ.11
206  & .AND.ielmh.EQ.15
207  & .AND.(ielmu.EQ.11.OR.ielmu.EQ.12)
208  & .AND.(ielmv.EQ.11.OR.ielmv.EQ.12) ) THEN
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  i1 = ikle1(ielem)
218  i2 = ikle2(ielem)
219  i3 = ikle3(ielem)
220  h1 = f(i1)
221  h2 = f(i2)
222  h3 = f(i3)
223  det=x2*y3-x3*y2
224  z1=h(ielem)
225  z2=h(ielem+ nelmax)-z1
226  z3=h(ielem+2*nelmax)-z1
227  zx=(z2*y3-z3*y2)/det
228  zy=(x2*z3-x3*z2)/det
229  u1 = u(i1) + g(i1)*zx
230  u2 = u(i2) + g(i2)*zx
231  u3 = u(i3) + g(i3)*zx
232  v1 = v(i1) + g(i1)*zy
233  v2 = v(i2) + g(i2)*zy
234  v3 = v(i3) + g(i3)*zy
235 !
236  h123 = h1+h2+h3
237  u123 = u1+u2+u3
238  v123 = v1+v2+v3
239 !
240  hu123 = h1*u1+h2*u2+h3*u3
241  hv123 = h1*v1+h2*v2+h3*v3
242 !
243  w1(ielem) = ( (y2-y3)*(h123*u123+hu123)
244  & +(x3-x2)*(h123*v123+hv123) )*xsur24
245  w2(ielem) = ( y3 *(h123*u123+hu123)
246  & -x3 *(h123*v123+hv123) )*xsur24
247  w3(ielem) = ( -y2 *(h123*u123+hu123)
248  & +x2 *(h123*v123+hv123) )*xsur24
249 !
250  ENDDO
251 !
252 ! F, G AND H ARE LINEAR, U AND V ARE LINEAR OR QUASI-BUBBLE
253 !
254  ELSEIF( ielmf.EQ.11
255  & .AND.ielmg.EQ.11
256  & .AND.ielmh.EQ.11
257  & .AND.(ielmu.EQ.11.OR.ielmu.EQ.12)
258  & .AND.(ielmv.EQ.11.OR.ielmv.EQ.12) ) THEN
259 !
260  DO ielem = 1 , nelem
261 !
262  x2 = xel(ielem,2)
263  x3 = xel(ielem,3)
264  y2 = yel(ielem,2)
265  y3 = yel(ielem,3)
266 !
267  i1 = ikle1(ielem)
268  i2 = ikle2(ielem)
269  i3 = ikle3(ielem)
270  h1 = f(i1)
271  h2 = f(i2)
272  h3 = f(i3)
273  det=x2*y3-x3*y2
274  z1=h(i1)
275  z2=h(i2)-z1
276  z3=h(i3)-z1
277  zx=(z2*y3-z3*y2)/det
278  zy=(x2*z3-x3*z2)/det
279  u1 = u(i1) + g(i1)*zx
280  u2 = u(i2) + g(i2)*zx
281  u3 = u(i3) + g(i3)*zx
282  v1 = v(i1) + g(i1)*zy
283  v2 = v(i2) + g(i2)*zy
284  v3 = v(i3) + g(i3)*zy
285 !
286  h123 = h1+h2+h3
287  u123 = u1+u2+u3
288  v123 = v1+v2+v3
289 !
290  hu123 = h1*u1+h2*u2+h3*u3
291  hv123 = h1*v1+h2*v2+h3*v3
292 !
293  w1(ielem) = ( (y2-y3)*(h123*u123+hu123)
294  & +(x3-x2)*(h123*v123+hv123) )*xsur24
295  w2(ielem) = ( y3 *(h123*u123+hu123)
296  & -x3 *(h123*v123+hv123) )*xsur24
297  w3(ielem) = ( -y2 *(h123*u123+hu123)
298  & +x2 *(h123*v123+hv123) )*xsur24
299 !
300  ENDDO
301 !
302 !-----------------------------------------------------------------------
303 !
304  ELSE
305 !
306 !-----------------------------------------------------------------------
307 !
308  WRITE(lu,101) ielmf,sf%NAME
309  WRITE(lu,201) ielmg,sg%NAME
310  WRITE(lu,301) ielmh,sh%NAME
311  WRITE(lu,401) ielmu,su%NAME
312  WRITE(lu,401) ielmv,sv%NAME
313  WRITE(lu,601)
314  CALL plante(1)
315  stop
316 !
317  ENDIF
318 !
319 !-----------------------------------------------------------------------
320 !
321  ELSEIF(formul(1:8).EQ.'HUGRADP3') THEN
322 !
323 !-----------------------------------------------------------------------
324 !
325 !
326 ! F AND G ARE LINEAR, H IS PIECEWISE LINEAR
327 !
328  IF(ielmf.EQ.11.AND.ielmg.EQ.11.AND.ielmh.EQ.15) THEN
329 !
330  DO ielem = 1 , nelem
331 !
332  x2 = xel(ielem,2)
333  x3 = xel(ielem,3)
334  y2 = yel(ielem,2)
335  y3 = yel(ielem,3)
336 !
337  i1 = ikle1(ielem)
338  i2 = ikle2(ielem)
339  i3 = ikle3(ielem)
340  h1 = f(i1)
341  h2 = f(i2)
342  h3 = f(i3)
343  det=x2*y3-x3*y2
344  z1=h(ielem)
345  z2=h(ielem+ nelmax)-z1
346  z3=h(ielem+2*nelmax)-z1
347  zx=(z2*y3-z3*y2)/det
348  zy=(x2*z3-x3*z2)/det
349  u1 = g(i1)*zx
350  u2 = g(i2)*zx
351  u3 = g(i3)*zx
352  v1 = g(i1)*zy
353  v2 = g(i2)*zy
354  v3 = g(i3)*zy
355 !
356  h123 = h1+h2+h3
357  u123 = u1+u2+u3
358  v123 = v1+v2+v3
359 !
360  hu123 = h1*u1+h2*u2+h3*u3
361  hv123 = h1*v1+h2*v2+h3*v3
362 !
363  w1(ielem) = ( (y2-y3)*(h123*u123+hu123)
364  & +(x3-x2)*(h123*v123+hv123) )*xsur24
365  w2(ielem) = ( y3 *(h123*u123+hu123)
366  & -x3 *(h123*v123+hv123) )*xsur24
367  w3(ielem) = ( -y2 *(h123*u123+hu123)
368  & +x2 *(h123*v123+hv123) )*xsur24
369 !
370  ENDDO
371 !
372  ELSEIF( ielmf.EQ.11
373  & .AND.ielmg.EQ.11
374  & .AND.ielmh.EQ.11 ) THEN
375 !
376 ! F, G AND H ARE LINEAR
377 !
378  DO ielem = 1 , nelem
379 !
380  x2 = xel(ielem,2)
381  x3 = xel(ielem,3)
382  y2 = yel(ielem,2)
383  y3 = yel(ielem,3)
384 !
385  i1 = ikle1(ielem)
386  i2 = ikle2(ielem)
387  i3 = ikle3(ielem)
388  h1 = f(i1)
389  h2 = f(i2)
390  h3 = f(i3)
391  det=x2*y3-x3*y2
392  z1=h(i1)
393  z2=h(i2)-z1
394  z3=h(i3)-z1
395  zx=(z2*y3-z3*y2)/det
396  zy=(x2*z3-x3*z2)/det
397  u1 = g(i1)*zx
398  u2 = g(i2)*zx
399  u3 = g(i3)*zx
400  v1 = g(i1)*zy
401  v2 = g(i2)*zy
402  v3 = g(i3)*zy
403 !
404  h123 = h1+h2+h3
405  u123 = u1+u2+u3
406  v123 = v1+v2+v3
407 !
408  hu123 = h1*u1+h2*u2+h3*u3
409  hv123 = h1*v1+h2*v2+h3*v3
410 !
411  w1(ielem) = ( (y2-y3)*(h123*u123+hu123)
412  & +(x3-x2)*(h123*v123+hv123) )*xsur24
413  w2(ielem) = ( y3 *(h123*u123+hu123)
414  & -x3 *(h123*v123+hv123) )*xsur24
415  w3(ielem) = ( -y2 *(h123*u123+hu123)
416  & +x2 *(h123*v123+hv123) )*xsur24
417 !
418  ENDDO
419 !
420  ELSE
421  WRITE(lu,101) ielmf,sf%NAME
422  WRITE(lu,201) ielmg,sg%NAME
423  WRITE(lu,301) ielmh,sh%NAME
424  WRITE(lu,601)
425  CALL plante(1)
426  stop
427  ENDIF
428 !
429 !-----------------------------------------------------------------------
430 !
431  ELSE
432 !
433  WRITE(lu,2000) formul
434 2000 FORMAT(1x,'VC19AA (BIEF):',/,
435  & 1x,'FORMULA: ',a16,' UNEXPECTED')
436  CALL plante(1)
437  stop
438 !
439  ENDIF
440 !
441 !-----------------------------------------------------------------------
442 !
443  RETURN
444  END
declarations_special
Definition: declarations_special.F:3
vc19aa
subroutine vc19aa(XMUL, SF, SG, SH, SU, SV, F, G, H, U, V, XEL, YEL, IKLE1, IKLE2, IKLE3, NELEM, NELMAX, W1, W2, W3, FORMUL)
Definition: vc19aa.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:57 for The TELEMAC-MASCARET system by doxygen 1.8.13