The TELEMAC-MASCARET system  trunk
vc04aa.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE vc04aa
3 ! *****************
4 !
5  &(xmul,su,sv,u,v,xel,yel,
6  & ikle1,ikle2,ikle3,nelem,nelmax,w1,w2,w3,specad)
7 !
8 !***********************************************************************
9 ! BIEF V7P0 21/08/2010
10 !***********************************************************************
11 !
12 !brief COMPUTES THE FOLLOWING VECTOR IN FINITE ELEMENTS:
13 !code
14 !+ / DPSI DPSI
15 !+ V = XMUL / ( U -- + V -- ) D(OMEGA)
16 !+ I /OMEGA DX DY
17 !+
18 !+ PSI(I) IS A BASE OF TYPE P1 TRIANGLE
19 !
20 !note SEE TESTS ON FORMUL TO UNDERSTAND HOW U AND V ARE DEALT WITH:
21 !+
22 !+ IF FORMUL(7:7).EQ.' ' : VELOCITY WITH COMPONENTS U AND V LINEAR.
23 !+
24 !+ IF FORMUL(7:7).EQ.'2' : VELOCITY EQUALS U*GRAD(V) U AND V LINEAR
25 !+ BUT V MAY BE PIECE-WISE LINEAR SO V%DIMDISC
26 !+ IS LOOKED AT IN THIS CASE.
27 !
28 !warning THE JACOBIAN MUST BE POSITIVE
29 !warning THE RESULT IS IN W IN NOT ASSEMBLED FORM
30 !
31 !history J-M HERVOUET (LNH)
32 !+ 01/06/06
33 !+ V5P7
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 !+ 12/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 !| IKLE1 |-->| FIRST POINT OF TRIANGLES
57 !| IKLE2 |-->| SECOND POINT OF TRIANGLES
58 !| IKLE3 |-->| THIRD POINT OF TRIANGLES
59 !| NELEM |-->| NUMBER OF ELEMENTS
60 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
61 !| SPECAD |-->| IF YES, SPECIAL ADVECTION FIELD, SEE ABOVE
62 !| SU |-->| BIEF_OBJ STRUCTURE OF U
63 !| SV |-->| BIEF_OBJ STRUCTURE OF V
64 !| SURFAC |-->| AREA OF TRIANGLES
65 !| U |-->| FUNCTION USED IN THE VECTOR FORMULA
66 !| V |-->| FUNCTION USED IN THE VECTOR FORMULA
67 !| W1 |<--| RESULT IN NON ASSEMBLED FORM
68 !| W2 |<--| RESULT IN NON ASSEMBLED FORM
69 !| W3 |<--| RESULT IN NON ASSEMBLED FORM
70 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
71 !| XMUL |-->| MULTIPLICATION COEFFICIENT
72 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
73 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
74 !
75  USE bief, ex_vc04aa => vc04aa
76 !
77  USE declarations_special
78  IMPLICIT NONE
79 !
80 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
81 !
82  INTEGER, INTENT(IN) :: NELEM,NELMAX
83  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax),IKLE3(nelmax)
84 !
85  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,*),YEL(nelmax,*)
86  DOUBLE PRECISION, INTENT(INOUT):: W1(nelmax),W2(nelmax),W3(nelmax)
87  DOUBLE PRECISION, INTENT(IN) :: XMUL
88 !
89 ! STRUCTURES OF F, U, V AND REAL DATA
90 !
91  TYPE(bief_obj), INTENT(IN) :: SU,SV
92  DOUBLE PRECISION, INTENT(IN) :: U(*),V(*)
93 !
94  LOGICAL, INTENT(IN) :: SPECAD
95 !
96 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
97 !
98  INTEGER IELEM,IELMU,IELMV
99  DOUBLE PRECISION XSUR6,X2,Y2,X3,Y3,U1,U2,U3,V1,V2,V3
100  DOUBLE PRECISION GRADVX,GRADVY,DET
101 !
102 !-----------------------------------------------------------------------
103 !
104  ielmu=su%ELM
105  ielmv=sv%ELM
106 !
107  xsur6 = xmul/6.d0
108 !
109 !-----------------------------------------------------------------------
110 !
111  IF(.NOT.specad) THEN
112 !
113 ! VELOCITY WITH LINEAR COMPONENTS U AND V
114 !
115  IF(ielmu.NE.11.OR.ielmv.NE.11) THEN
116  WRITE(lu,201) ielmu,su%NAME
117 201 FORMAT(1x,'VC04AA: DISCRETIZATION OF U:',1i6,
118  & 1x,'REAL NAME: ',a6)
119 301 FORMAT(1x,'CASE NOT IMPLEMENTED')
120  CALL plante(1)
121  stop
122  ENDIF
123 !
124  DO ielem = 1 , nelem
125  x2 = xel(ielem,2)
126  x3 = xel(ielem,3)
127  y2 = yel(ielem,2)
128  y3 = yel(ielem,3)
129  u1 = u(ikle1(ielem))
130  u2 = u(ikle2(ielem))
131  u3 = u(ikle3(ielem))
132  v1 = v(ikle1(ielem))
133  v2 = v(ikle2(ielem))
134  v3 = v(ikle3(ielem))
135  w1(ielem) = (u3*y2-u3*y3-v3*x2+v3*x3+u2*y2-u2*y3
136  & -v2*x2+v2*x3+u1*y2-u1*y3-v1*x2+v1*x3)*xsur6
137  w2(ielem) = (u1*y3+u2*y3-v1*x3-v2*x3+u3*y3-v3*x3)*xsur6
138  w3(ielem) = (-u1*y2-u2*y2+v1*x2+v2*x2-u3*y2+v3*x2)*xsur6
139  ENDDO
140 !
141 !-----------------------------------------------------------------------
142 !
143  ELSE
144 !
145 ! VELOCITY EQUALS U * GRAD(V)
146 ! WITH VARIOUS DISCRETISATIONS OF U AND V
147 !
148  IF(ielmv.EQ.11) THEN
149 !
150  IF(ielmu.EQ.11) THEN
151 ! V LINEAR, U LINEAR
152  DO ielem = 1 , nelem
153  x2 = xel(ielem,2)
154  x3 = xel(ielem,3)
155  y2 = yel(ielem,2)
156  y3 = yel(ielem,3)
157  det=x2*y3-x3*y2
158  v1 = v(ikle1(ielem))
159  v2 = v(ikle2(ielem))
160  v3 = v(ikle3(ielem))
161  gradvx=((v2-v1)*y3-(v3-v1)*y2)/det
162  gradvy=(x2*(v3-v1)-x3*(v2-v1))/det
163  u1 = u(ikle1(ielem))*gradvx
164  u2 = u(ikle2(ielem))*gradvx
165  u3 = u(ikle3(ielem))*gradvx
166  v1 = u(ikle1(ielem))*gradvy
167  v2 = u(ikle2(ielem))*gradvy
168  v3 = u(ikle3(ielem))*gradvy
169  w1(ielem) = (u3*y2-u3*y3-v3*x2+v3*x3+u2*y2-u2*y3
170  & -v2*x2+v2*x3+u1*y2-u1*y3-v1*x2+v1*x3)*xsur6
171  w2(ielem) = (u1*y3+u2*y3-v1*x3-v2*x3+u3*y3-v3*x3)*xsur6
172  w3(ielem) = (-u1*y2-u2*y2+v1*x2+v2*x2-u3*y2+v3*x2)*xsur6
173  ENDDO
174  ELSEIF(ielmu.EQ.10) THEN
175 ! V LINEAR, U PIECE-WISE CONSTANT
176  DO ielem = 1 , nelem
177  x2 = xel(ielem,2)
178  x3 = xel(ielem,3)
179  y2 = yel(ielem,2)
180  y3 = yel(ielem,3)
181  det=x2*y3-x3*y2
182  v1 = v(ikle1(ielem))
183  v2 = v(ikle2(ielem))
184  v3 = v(ikle3(ielem))
185  gradvx=((v2-v1)*y3-(v3-v1)*y2)/det
186  gradvy=(x2*(v3-v1)-x3*(v2-v1))/det
187  u1 = u(ielem)*gradvx
188  u2 = u(ielem)*gradvx
189  u3 = u(ielem)*gradvx
190  v1 = u(ielem)*gradvy
191  v2 = u(ielem)*gradvy
192  v3 = u(ielem)*gradvy
193  w1(ielem) = (u3*y2-u3*y3-v3*x2+v3*x3+u2*y2-u2*y3
194  & -v2*x2+v2*x3+u1*y2-u1*y3-v1*x2+v1*x3)*xsur6
195  w2(ielem) = (u1*y3+u2*y3-v1*x3-v2*x3+u3*y3-v3*x3)*xsur6
196  w3(ielem) = (-u1*y2-u2*y2+v1*x2+v2*x2-u3*y2+v3*x2)*xsur6
197  ENDDO
198  ELSE
199  WRITE(lu,*) 'WRONG DISCRETISATION OF U IN VC04AA'
200  CALL plante(1)
201  stop
202  ENDIF
203 !
204  ELSEIF(ielmv.EQ.15) THEN
205 !
206  IF(ielmu.EQ.11) THEN
207 ! V PIECE-WISE LINEAR, U LINEAR
208  DO ielem=1,nelem
209  x2 = xel(ielem,2)
210  x3 = xel(ielem,3)
211  y2 = yel(ielem,2)
212  y3 = yel(ielem,3)
213  det=x2*y3-x3*y2
214  v1 = v(ielem )
215  v2 = v(ielem+ nelmax)
216  v3 = v(ielem+2*nelmax)
217  gradvx=((v2-v1)*y3-(v3-v1)*y2)/det
218  gradvy=(x2*(v3-v1)-x3*(v2-v1))/det
219  u1 = u(ikle1(ielem))*gradvx
220  u2 = u(ikle2(ielem))*gradvx
221  u3 = u(ikle3(ielem))*gradvx
222  v1 = u(ikle1(ielem))*gradvy
223  v2 = u(ikle2(ielem))*gradvy
224  v3 = u(ikle3(ielem))*gradvy
225  w1(ielem) = (u3*y2-u3*y3-v3*x2+v3*x3+u2*y2-u2*y3
226  & -v2*x2+v2*x3+u1*y2-u1*y3-v1*x2+v1*x3)*xsur6
227  w2(ielem) = (u1*y3+u2*y3-v1*x3-v2*x3+u3*y3-v3*x3)*xsur6
228  w3(ielem) = (-u1*y2-u2*y2+v1*x2+v2*x2-u3*y2+v3*x2)*xsur6
229  ENDDO
230  ELSEIF(ielmu.EQ.10) THEN
231 ! V PIECE-WISE LINEAR, U LINEAR
232  DO ielem=1,nelem
233  x2 = xel(ielem,2)
234  x3 = xel(ielem,3)
235  y2 = yel(ielem,2)
236  y3 = yel(ielem,3)
237  det=x2*y3-x3*y2
238  v1 = v(ielem )
239  v2 = v(ielem+ nelmax)
240  v3 = v(ielem+2*nelmax)
241  gradvx=((v2-v1)*y3-(v3-v1)*y2)/det
242  gradvy=(x2*(v3-v1)-x3*(v2-v1))/det
243  u1 = u(ielem)*gradvx
244  u2 = u(ielem)*gradvx
245  u3 = u(ielem)*gradvx
246  v1 = u(ielem)*gradvy
247  v2 = u(ielem)*gradvy
248  v3 = u(ielem)*gradvy
249  w1(ielem) = (u3*y2-u3*y3-v3*x2+v3*x3+u2*y2-u2*y3
250  & -v2*x2+v2*x3+u1*y2-u1*y3-v1*x2+v1*x3)*xsur6
251  w2(ielem) = (u1*y3+u2*y3-v1*x3-v2*x3+u3*y3-v3*x3)*xsur6
252  w3(ielem) = (-u1*y2-u2*y2+v1*x2+v2*x2-u3*y2+v3*x2)*xsur6
253  ENDDO
254  ELSE
255  WRITE(lu,*) 'WRONG DISCRETISATION OF U IN VC04AA'
256  CALL plante(1)
257  stop
258  ENDIF
259 !
260  ELSE
261  WRITE(lu,501) sv%DIMDISC
262  WRITE(lu,301)
263 501 FORMAT(1x,'VC04AA: V%DIMDISC=',1i6)
264  CALL plante(1)
265  stop
266  ENDIF
267 !
268  ENDIF
269 !
270 !-----------------------------------------------------------------------
271 !
272  RETURN
273  END
declarations_special
Definition: declarations_special.F:3
vc04aa
subroutine vc04aa(XMUL, SU, SV, U, V, XEL, YEL, IKLE1, IKLE2, IKLE3, NELEM, NELMAX, W1, W2, W3, SPECAD)
Definition: vc04aa.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:50 for The TELEMAC-MASCARET system by doxygen 1.8.13