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