The TELEMAC-MASCARET system  trunk
vc13bb.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE vc13bb
3 ! *****************
4 !
5  &(xmul,sf,f,xel,yel,
6  & ikle1,ikle2,ikle3,ikle4,nelem,nelmax,w1,w2,w3,w4,icoord)
7 !
8 !***********************************************************************
9 ! BIEF V7P0 21/08/2010
10 !***********************************************************************
11 !
12 !brief COMPUTES THE FOLLOWING VECTOR IN FINITE ELEMENTS:
13 !code
14 !+ / --->
15 !+ VEC(I) = XMUL / PSI(I) * GRAD(F) D(OMEGA)
16 !+ /OMEGA
17 !+
18 !+ PSI(I) IS A BASE OF QUASI-BUBBLE TYPE
19 !
20 !note IMPORTANT : IF F IS OF TYPE P0, THE RESULT IS 0.
21 !+ HERE, IF F IS P0, IT REALLY MEANS THAT F IS
22 !+ P1, BUT GIVEN BY ELEMENTS.
23 !+ THE SIZE OF F SHOULD THEN BE : F(NELMAX,4).
24 !
25 !warning THE JACOBIAN MUST BE POSITIVE
26 !warning THE RESULT IS IN W IN NOT ASSEMBLED FORM
27 !
28 !history J-M HERVOUET (LNH) ; C MOULIN (LNH)
29 !+ 10/01/95
30 !+ V5P1
31 !+
32 !
33 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
34 !+ 13/07/2010
35 !+ V6P0
36 !+ Translation of French comments within the FORTRAN sources into
37 !+ English comments
38 !
39 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
40 !+ 21/08/2010
41 !+ V6P0
42 !+ Creation of DOXYGEN tags for automated documentation and
43 !+ cross-referencing of the FORTRAN sources
44 !
45 !history J-M HERVOUET (EDF LAB, LNHE)
46 !+ 12/05/2014
47 !+ V7P0
48 !+ Discontinuous elements better treated: new types 15, 16 and 17 for
49 !+ discontinuous linear, quasi-bubble, and quadratic, rather than
50 !+ using component DIMDISC=11, 12 or 13.
51 !
52 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
53 !| F |-->| FUNCTION USED IN THE VECTOR FORMULA
54 !| ICOORD |-->| 1: DERIVATIVE ALONG X, 2: ALONG Y
55 !| IKLE1 |-->| FIRST POINT OF TRIANGLES
56 !| IKLE2 |-->| SECOND POINT OF TRIANGLES
57 !| IKLE3 |-->| THIRD POINT OF TRIANGLES
58 !| IKLE4 |-->| QUASI-BUBBLE POINT OF TRIANGLES
59 !| NELEM |-->| NUMBER OF ELEMENTS
60 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
61 !| SF |-->| BIEF_OBJ STRUCTURE OF F
62 !| SURFAC |-->| AREA OF TRIANGLES
63 !| W1 |<--| RESULT IN NON ASSEMBLED FORM
64 !| W2 |<--| RESULT IN NON ASSEMBLED FORM
65 !| W3 |<--| RESULT IN NON ASSEMBLED FORM
66 !| W4 |<--| RESULT IN NON ASSEMBLED FORM
67 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
68 !| XMUL |-->| MULTIPLICATION COEFFICIENT
69 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
70 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
71 !
72  USE bief, ex_vc13bb => vc13bb
73 !
74  USE declarations_special
75  IMPLICIT NONE
76 !
77 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
78 !
79  INTEGER, INTENT(IN) :: NELEM,NELMAX,ICOORD
80  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax)
81  INTEGER, INTENT(IN) :: IKLE3(nelmax),IKLE4(nelmax)
82 !
83  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,*),YEL(nelmax,*)
84  DOUBLE PRECISION, INTENT(INOUT) :: W1(nelmax),W2(nelmax)
85  DOUBLE PRECISION, INTENT(INOUT) :: W3(nelmax),W4(nelmax)
86  DOUBLE PRECISION, INTENT(IN) :: XMUL
87 !
88 ! STRUCTURE OF F AND REAL DATA
89 !
90  TYPE(bief_obj), INTENT(IN) :: SF
91  DOUBLE PRECISION, INTENT(IN) :: F(*)
92 !
93 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
94 !
95  INTEGER IELEM,IELMF
96  DOUBLE PRECISION F1,F2,F3,F4,X2,X3,Y2,Y3,XSUR9,XSUR6,XSUR18
97 !
98 !-----------------------------------------------------------------------
99 !
100  ielmf=sf%ELM
101  xsur9 = xmul / 9.d0
102  xsur6 = xmul / 6.d0
103  xsur18= xmul /18.d0
104 !
105 !-----------------------------------------------------------------------
106 ! F OF TYPE P1
107 !-----------------------------------------------------------------------
108 !
109  IF(ielmf.EQ.11) THEN
110 !
111 !================================
112 ! DERIVATIVE WRT X =
113 !================================
114 !
115  IF(icoord.EQ.1) THEN
116 !
117  DO ielem = 1 , nelem
118 !
119  f1 = f(ikle1(ielem))
120  f2 = f(ikle2(ielem))
121  f3 = f(ikle3(ielem))
122 !
123  y2 = yel(ielem,2)
124  y3 = yel(ielem,3)
125 !
126  w1(ielem) = ( - (f3-f1)*y2 + (f2-f1)*y3 ) * xsur9
127  w2(ielem) = w1(ielem)
128  w3(ielem) = w1(ielem)
129  w4(ielem) = 1.5d0*w1(ielem)
130 !
131  ENDDO ! IELEM
132 !
133 !================================
134 ! DERIVATIVE WRT Y =
135 !================================
136 !
137  ELSEIF(icoord.EQ.2) THEN
138 !
139  DO ielem = 1 , nelem
140 !
141  f1 = f(ikle1(ielem))
142  f2 = f(ikle2(ielem))
143  f3 = f(ikle3(ielem))
144 !
145  x2 = xel(ielem,2)
146  x3 = xel(ielem,3)
147 !
148  w1(ielem) = ( x2*(f3-f1) - x3*(f2-f1) ) * xsur9
149  w2(ielem) = w1(ielem)
150  w3(ielem) = w1(ielem)
151  w4(ielem) = 1.5d0*w1(ielem)
152 !
153  ENDDO
154 !
155  ELSE
156 !
157  WRITE(lu,201) icoord
158  CALL plante(1)
159  stop
160 !
161  ENDIF
162 !
163 !-----------------------------------------------------------------------
164 !
165 ! BEWARE: HERE F IS LINEAR BUT DISCONTINUOUS BETWEEN THE ELEMENTS
166 !
167  ELSEIF(ielmf.EQ.15) THEN
168 !
169 ! X COORDINATE
170 !
171  IF(icoord.EQ.1) THEN
172 !
173  DO ielem = 1 , nelem
174 !
175  f1 = f(ielem)
176  f2 = f(ielem+nelmax)
177  f3 = f(ielem+2*nelmax)
178 !
179  y2 = yel(ielem,2)
180  y3 = yel(ielem,3)
181 !
182  w1(ielem) = ( - (f3-f1)*y2 + (f2-f1)*y3 ) * xsur9
183  w2(ielem) = w1(ielem)
184  w3(ielem) = w1(ielem)
185  w4(ielem) = 1.5d0*w1(ielem)
186 !
187  ENDDO
188 !
189  ELSEIF(icoord.EQ.2) THEN
190 !
191 ! Y COORDINATE
192 !
193  DO ielem = 1 , nelem
194 !
195  f1 = f(ielem)
196  f2 = f(ielem+nelmax)
197  f3 = f(ielem+2*nelmax)
198 !
199  x2 = xel(ielem,2)
200  x3 = xel(ielem,3)
201 !
202  w1(ielem) = ( x2*(f3-f1) - x3*(f2-f1) ) * xsur9
203  w2(ielem) = w1(ielem)
204  w3(ielem) = w1(ielem)
205  w4(ielem) = 1.5d0*w1(ielem)
206 !
207  ENDDO
208 !
209  ELSE
210 !
211  WRITE(lu,201) icoord
212  CALL plante(1)
213  stop
214 !
215  ENDIF
216 !
217 !-----------------------------------------------------------------------
218 ! F IS QUASI-BUBBLE
219 !-----------------------------------------------------------------------
220 !
221  ELSEIF(ielmf.EQ.12) THEN
222 !
223 !================================
224 ! DERIVATIVE WRT X =
225 !================================
226 !
227  IF(icoord.EQ.1) THEN
228 !
229  DO ielem = 1 , nelem
230 !
231  f1 = f(ikle1(ielem))
232  f2 = f(ikle2(ielem))
233  f3 = f(ikle3(ielem))
234  f4 = f(ikle4(ielem))
235 !
236  y2 = yel(ielem,2)
237  y3 = yel(ielem,3)
238 !
239  w1(ielem) = -(( y3 +y2)*(f3-f2)-3*(y3-y2)*(f4-f1))*xsur18
240  w2(ielem) = (( y3-2*y2)*(f3-f1)-3*(f4-f2)*y3 )*xsur18
241  w3(ielem) = ((2*y3 -y2)*(f2-f1)+3*(f4-f3)*y2)*xsur18
242  w4(ielem) = ( -(y3 -y2)* f1-f3*y2+f2*y3)*xsur6
243 !
244  ENDDO ! IELEM
245 !
246 !================================
247 ! DERIVATIVE WRT Y =
248 !================================
249 !
250  ELSEIF(icoord.EQ.2) THEN
251 !
252  DO ielem = 1 , nelem
253 !
254  f1 = f(ikle1(ielem))
255  f2 = f(ikle2(ielem))
256  f3 = f(ikle3(ielem))
257  f4 = f(ikle4(ielem))
258 !
259  x2 = xel(ielem,2)
260  x3 = xel(ielem,3)
261 !
262  w1(ielem)=((x2+x3)*(f3-f2)+3*(x2-x3)*(f4-f1))*xsur18
263  w2(ielem)=((2*x2-x3)*(f3-f1)+3*x3*(f4-f2))*xsur18
264  w3(ielem)=((x2-2*x3)*(f2-f1)-3*x2*(f4-f3))*xsur18
265  w4(ielem)=(x2*(f3-f1)-x3*(f2-f1))*xsur6
266 !
267  ENDDO
268 !
269  ELSE
270 !
271  WRITE(lu,201) icoord
272  CALL plante(1)
273  stop
274 !
275  ENDIF
276 !
277 !-----------------------------------------------------------------------
278 ! F IS QUASI-BUBBLE BUT DISCONTINUOUS BETWEEN THE ELEMENTS
279 !-----------------------------------------------------------------------
280 !
281  ELSEIF(ielmf.EQ.16) THEN
282 !
283 !================================
284 ! DERIVATIVE WRT X =
285 !================================
286 !
287  IF(icoord.EQ.1) THEN
288 !
289  DO ielem = 1 , nelem
290 !
291  f1 = f(ielem )
292  f2 = f(ielem+ nelmax)
293  f3 = f(ielem+2*nelmax)
294  f4 = f(ielem+3*nelmax)
295 !
296  y2 = yel(ielem,2)
297  y3 = yel(ielem,3)
298 !
299  w1(ielem)=-((y3+y2)*(f3-f2)-3*(y3-y2)*(f4-f1))*xsur18
300  w2(ielem)=((y3-2*y2)*(f3-f1)-3*(f4-f2)*y3)*xsur18
301  w3(ielem)=((2*y3-y2)*(f2-f1)+3*(f4-f3)*y2)*xsur18
302  w4(ielem)=(-((y3-y2)*f1+f3*y2-f2*y3))*xsur6
303 !
304  ENDDO
305 !
306 !================================
307 ! DERIVATIVE WRT Y =
308 !================================
309 !
310  ELSEIF(icoord.EQ.2) THEN
311 !
312  DO ielem = 1 , nelem
313 !
314  f1 = f(ielem )
315  f2 = f(ielem+ nelmax)
316  f3 = f(ielem+2*nelmax)
317  f4 = f(ielem+3*nelmax)
318 !
319  x2 = xel(ielem,2)
320  x3 = xel(ielem,3)
321 !
322  w1(ielem)=((x2+x3)*(f3-f2)+3*(x2-x3)*(f4-f1))*xsur18
323  w2(ielem)=((2*x2-x3)*(f3-f1)+3*x3*(f4-f2))*xsur18
324  w3(ielem)=((x2-2*x3)*(f2-f1)-3*x2*(f4-f3))*xsur18
325  w4(ielem)=(x2*(f3-f1)-x3*(f2-f1))*xsur6
326 !
327  ENDDO
328 !
329  ELSE
330 !
331  WRITE(lu,201) icoord
332  CALL plante(1)
333  stop
334 !
335  ENDIF
336 !
337 !-----------------------------------------------------------------------
338 ! OTHERS
339 ! ELSEIF
340 !-----------------------------------------------------------------------
341 !
342  ELSE
343 !
344  WRITE(lu,101) ielmf,sf%NAME
345 101 FORMAT(1x,'VC13BB (BIEF) :',/,
346  & 1x,'DISCRETIZATION OF F NOT AVAILABLE:',1i6,
347  & 1x,'REAL NAME: ',a6)
348  CALL plante(1)
349  stop
350 !
351  ENDIF
352 !
353 201 FORMAT(1x,'VC13BB (BIEF) : IMPOSSIBLE COMPONENT ',
354  & 1i6,' CHECK ICOORD')
355 !
356 !-----------------------------------------------------------------------
357 !
358  RETURN
359  END
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
vc13bb
subroutine vc13bb(XMUL, SF, F, XEL, YEL, IKLE1, IKLE2, IKLE3, IKLE4, NELEM, NELMAX, W1, W2, W3, W4, ICOORD)
Definition: vc13bb.f:8
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:55 for The TELEMAC-MASCARET system by doxygen 1.8.13