The TELEMAC-MASCARET system  trunk
vc21pp.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE vc21pp
3 ! *****************
4 !
5  &( xmul,sf,f,x,y,z,surfac,
6  & ikle1,ikle2,ikle3,ikle4,ikle5,ikle6,nelem,nelmax,
7  & w1,w2,w3,w4,w5,w6)
8 !
9 !***********************************************************************
10 ! BIEF V8P0 21/06/2018
11 !***********************************************************************
12 !
13 !brief COMPUTES THE FOLLOWING VECTOR IN FINITE ELEMENTS:
14 !code
15 !+ /
16 !+ VEC(I) = XMUL / PSI(I) * F D(OMEGA)
17 !+ /OMEGA
18 !+
19 !+ PSI(I) IS A BASE OF TYPE P1 PRISM
20 !+
21 !+ F IS A VECTOR OF TYPE IELMF
22 !
23 !warning THE JACOBIAN MUST BE POSITIVE
24 !warning THE RESULT IS IN W IN NOT ASSEMBLED FORM - REAL MESH
25 !
26 !history J-M HERVOUET (LNH) ; F LEPEINTRE (LNH)
27 !+ 09/12/94
28 !+ V5P1
29 !+
30 !
31 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
32 !+ 13/07/2010
33 !+ V6P0
34 !+ Translation of French comments within the FORTRAN sources into
35 !+ English comments
36 !
37 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
38 !+ 21/08/2010
39 !+ V6P0
40 !+ Creation of DOXYGEN tags for automated documentation and
41 !+ cross-referencing of the FORTRAN sources
42 !
43 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
44 !| F |-->| FUNCTION USED IN THE VECTOR FORMULA
45 !| IKLE1 |-->| FIRST POINT OF PRISMS
46 !| IKLE2 |-->| SECOND POINT OF PRISMS
47 !| IKLE3 |-->| THIRD POINT OF PRISMS
48 !| IKLE4 |-->| FOURTH POINT OF PRISMS
49 !| IKLE5 |-->| FIFTH POINT OF PRISMS
50 !| IKLE6 |-->| SIXTH POINT OF PRISMS
51 !| NELEM |-->| NUMBER OF ELEMENTS
52 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
53 !| SF |-->| BIEF_OBJ STRUCTURE OF F
54 !| SURFAC |-->| AREA OF TRIANGLES
55 !| W1 |<--| RESULT IN NON ASSEMBLED FORM
56 !| W2 |<--| RESULT IN NON ASSEMBLED FORM
57 !| W3 |<--| RESULT IN NON ASSEMBLED FORM
58 !| W4 |<--| RESULT IN NON ASSEMBLED FORM
59 !| W5 |<--| RESULT IN NON ASSEMBLED FORM
60 !| W6 |<--| RESULT IN NON ASSEMBLED FORM
61 !| XMUL |-->| MULTIPLICATION COEFFICIENT
62 !| Z |-->| ELEVATIONS OF POINTS
63 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
64 !
65  USE bief, ex_vc21pp => vc21pp
66 !
67  USE declarations_special
68  IMPLICIT NONE
69 !
70 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
71 !
72  INTEGER, INTENT(IN) :: NELEM,NELMAX
73  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax),IKLE3(nelmax)
74  INTEGER, INTENT(IN) :: IKLE4(nelmax),IKLE5(nelmax),IKLE6(nelmax)
75 !
76  DOUBLE PRECISION, INTENT(IN) :: X(nelmax,6),Y(nelmax,6),Z(*)
77  DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
78  DOUBLE PRECISION,INTENT(INOUT)::W1(nelmax),W2(nelmax),W3(nelmax)
79  DOUBLE PRECISION,INTENT(INOUT)::W4(nelmax),W5(nelmax),W6(nelmax)
80  DOUBLE PRECISION, INTENT(IN) :: XMUL
81 !
82 ! STRUCTURE OF F AND REAL DATA
83 !
84  TYPE(bief_obj), INTENT(IN) :: SF
85  DOUBLE PRECISION, INTENT(IN) :: F(*)
86 !
87 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
88 !
89  INTEGER IELEM,IELMF,I1,I2,I3,I4,I5,I6
90  DOUBLE PRECISION COEF,H1,H2,H3,SHT,DIFF,X2,Y2,X3,Y3,XS24,XS144
91  DOUBLE PRECISION Z2,Z3,Z4,Z5,Z6
92  DOUBLE PRECISION F1,F2,F3,F4,F5,F6
93 !
94 !***********************************************************************
95 !
96  ielmf=sf%ELM
97 !
98 !-----------------------------------------------------------------------
99 !
100 ! F IS LINEAR
101 
102  xs24 = xmul/24.d0
103  xs144 = xmul/144.d0
104 !
105  IF(ielmf.EQ.41) THEN
106  DO ielem = 1 , nelem
107 !
108  i1 = ikle1(ielem)
109  i2 = ikle2(ielem)
110  i3 = ikle3(ielem)
111  i4 = ikle4(ielem)
112  i5 = ikle5(ielem)
113  i6 = ikle6(ielem)
114 !
115  f1 = f(i1)
116  f2 = f(i2)
117  f3 = f(i3)
118  f4 = f(i4)
119  f5 = f(i5)
120  f6 = f(i6)
121 !
122 ! REAL COORDINATES OF THE POINTS OF THE ELEMENT (ORIGIN IN 1)
123 !
124 ! Y2 = Y(I2) - Y(I1)
125 ! Y3 = Y(I3) - Y(I1)
126  x2 = x(ielem,2)
127  x3 = x(ielem,3)
128  y2 = y(ielem,2)
129  y3 = y(ielem,3)
130 !
131  z2 = z(i2) - z(i1)
132  z3 = z(i3) - z(i1)
133  z4 = z(i4) - z(i1)
134  z5 = z(i5) - z(i1)
135  z6 = z(i6) - z(i1)
136 
137 
138  h1 = z(ikle4(ielem)) - z(ikle1(ielem))
139  h2 = z(ikle5(ielem)) - z(ikle2(ielem))
140  h3 = z(ikle6(ielem)) - z(ikle3(ielem))
141  sht = h1 + h2 + h3
142 !!
143  w1(ielem)=(((2*f1-f6)*y2*( z5+3*z4-3*z3 -z2)
144  & +(2*f1-f5)*y3*(-z6-3*z4 +z3+3*z2)
145  & +f2*y3*(2*z6+3*z5+3*z4-2*z3)
146  & +f3*y2*(-3*z6-2*z5-3*z4+2*z2)
147  & +(f3-f6)*y3*(z5-z4+2*z2)
148  & +f4*y2*(3*z6+z5+3*z3-z2)
149  & +f4*y3*(-z6-3*z5+z3-3*z2)
150  & +(f5-f2)*y2*(z6-z4+2*z3))*xs144)**2
151 
152  w2(ielem)=((f1*y2*(z6+4*z5+3*z4-4*z3-4*z2)
153  & +f1*y3*(-2*z6-3*z5-3*z4+2*z3+6*z2)
154  & +(2*f2-f4)*y3*(z6+3*z5-z3)
155  & +f3*y2*(-3*z6-4*z5-z4+4*z2)
156  & +f4*y2*(2*z6+2*z5+z3-2*z2)
157  & +2*(f5-f2)*y2*(z6-z4+2*z3)
158  & +f5*y3*(z6+3*z4-z3-6*z2)
159  & +f6*y2*(-2*z5-2*z4+3*z3+2*z2)
160  & +(f6-f3)*y3*(-z5+z4-2*z2) )*xs144)**2
161 
162  w3(ielem)=((f1*y2*(3*z6+2*z5+3*z4-6*z3-2*z2)
163  & +f1*y3*(-4*z6-z5-3*z4+4*z3+4*z2)
164  & +f2*y3*(4*z6+3*z5+z4-4*z3)
165  & +(2*f3-f4)*y2*(-3*z6-z5+z2)
166  & +f4*y3*(-2*z6-2*z5+2*z3-z2)
167  & +(f5-f2)*y2*(z6-z4+2*z3)
168  & +f5*y3*(2*z6+2*z4-2*z3-3*z2)
169  & +f6*y2*(-z5-3*z4+6*z3+z2)
170  & +2*(f6-f3)*y3*(-z5+z4-2*z2) )*xs144)**2
171 
172  w4(ielem)=((f1*y2*(-3*z6+z5+6*z4-3*z3-z2)
173  & +f1*y3*(-z6+3*z5-6*z4+z3+3*z2)
174  & +f2*y3*(z6+3*z5-z3)
175  & +(2*f4-f3)*y2*(3*z6+z5-z2)
176  & +2*f4*y3*(-z6-3*z5+z3)
177  & +(f5-f2)*y2*(2*z6-2*z4+z3)
178  & +f5*y3*(2*z6+6*z4-2*z3-3*z2)
179  & +f6*y2*(-2*z5-6*z4+3*z3+2*z2)
180  & +(f6-f3)*y3*(-2*z5+2*z4-z2) )*xs144)**2
181 
182  w5(ielem)=((f1*y2*(-z6+2*z5+3*z4-2*z3-2*z2)
183  & +f2*y3*(z6+6*z5-3*z4-z3)
184  & +f3*y2*(-3*z6-2*z5+z4+2*z2)
185  & +f4*y2*(4*z6+4*z5-z3-4*z2)
186  & +f4*y3*(-2*z6-6*z5+2*z3+3*z2)
187  & +2*(f5-f2)*y2*(2*z6-2*z4+z3)
188  & +(2*f5-f1)*y3*(z6+3*z4-z3-3*z2)
189  & +f6*y2*(-4*z5-4*z4+3*z3+4*z2)
190  & +(f6-f3)*y3*(-2*z5+2*z4-z2) )*xs144)**2
191 
192  w6(ielem)=((f1*y3*(-2*z6+z5-3*z4+2*z3+2*z2)
193  & +f2*y3*(2*z6+3*z5-z4-2*z3)
194  & +f3*y2*(-6*z6-z5+3*z4+z2)
195  & +f4*y2*(6*z6+2*z5-3*z3-2*z2)
196  & +f4*y3*(-4*z6-4*z5+4*z3+z2)
197  & +(f5-f2)*y2*(2*z6-2*z4+z3)
198  & +f5*y3*(4*z6+4*z4-4*z3-3*z2)
199  & +(2*f6-f1)*y2*(-z5-3*z4+3*z3+z2)
200  & +2*(f6-f3)*y3*(-2*z5+2*z4-z2))*xs144)**2
201 
202 
203 
204  w1(ielem)=w1(ielem)+(((2*f1-f6)*x2*(-z5-3*z4+3*z3+z2)
205  & +2*f1*x3*(z6+3*z4-z3-3*z2)
206  & +f2*x3*(-2*z6-3*z5-3*z4+2*z3)
207  & +f3*x2*(3*z6+2*z5+3*z4-2*z2)
208  & +f4*x2*(-3*z6-z5-3*z3+z2)
209  & +f4*x3*(z6+3*z5-z3+3*z2)
210  & +(f5-f2)*x2*(-z6+z4-2*z3)
211  & +f5*x3*(-z6-3*z4+z3+3*z2)
212  & +(f6-f3)*x3*(z5-z4+2*z2))*xs144)**2
213 
214  w2(ielem)=w2(ielem)+((f1*x2*(-z6-4*z5-3*z4+4*z3+4*z2)
215  & +f1*x3*(2*z6+3*z5+3*z4-2*z3-6*z2)
216  & +(2*f2-f4)*x3*(-z6-3*z5+z3)
217  & +f3*x2*(3*z6+4*z5+z4-4*z2)
218  & +f4*x2*(-2*z6-2*z5-z3+2*z2)
219  & +2*(f5-f2)*x2*(-z6+z4-2*z3)
220  & +f5*x3*(-z6-3*z4+z3+6*z2)
221  & +f6*x2*(2*z5+2*z4-3*z3-2*z2)
222  & +(f6-f3)*x3*(z5-z4+2*z2))*xs144)**2
223 
224  w3(ielem)=w3(ielem)+((f1*x2*(-3*z6-2*z5-3*z4+6*z3+2*z2)
225  & +f1*x3*(4*z6+z5+3*z4-4*z3-4*z2)
226  & +f2*x3*(-4*z6-3*z5-z4+4*z3)
227  & +(2*f3-f4)*x2*(3*z6+z5-z2)
228  & +f4*x3*(2*z6+2*z5-2*z3+z2)
229  & +(f5-f2)*x2*(-z6+z4-2*z3)
230  & +f5*x3*(-2*z6-2*z4+2*z3+3*z2)
231  & +f6*x2*(z5+3*z4-6*z3-z2)
232  & +2*(f6-f3)*x3*(z5-z4+2*z2))*xs144)**2
233 
234  w4(ielem)=w4(ielem)+((f1*x2*(3*z6-z5-6*z4+3*z3+z2)
235  & +f1*x3*(z6-3*z5+6*z4-z3-3*z2)
236  & +(2*f4-f3)*x2*(-3*z6-z5+z2)
237  & +(2*f4-f2)*x3*(z6+3*z5-z3)
238  & +(f5-f2)*x2*(-2*z6+2*z4-z3)
239  & +f5*x3*(-2*z6-6*z4+2*z3+3*z2)
240  & +f6*x2*(2*z5+6*z4-3*z3-2*z2)
241  & +(f6-f3)*x3*(2*z5-2*z4+z2))*xs144)**2
242 
243  w5(ielem)=w5(ielem)+((f1*x2*(z6-2*z5-3*z4+2*z3+2*z2)
244  & +f2*x3*(-z6-6*z5+3*z4+z3)
245  & +f3*x2*(3*z6+2*z5-z4-2*z2)
246  & +f4*x2*(-4*z6-4*z5+z3+4*z2)
247  & +f4*x3*(2*z6+6*z5-2*z3-3*z2)
248  & +2*(f5-f2)*x2*(-2*z6+2*z4-z3)
249  & +(2*f5-f1)*x3*(-z6-3*z4+z3+3*z2)
250  & +f6*x2*(4*z5+4*z4-3*z3-4*z2)
251  & +(f6-f3)*x3*(2*z5-2*z4+z2))*xs144)**2
252 
253  w6(ielem)=w6(ielem)+((f1*x3*(2*z6-z5+3*z4-2*z3-2*z2)
254  & +f2*x3*(-2*z6-3*z5+z4+2*z3)
255  & +f3*x2*(6*z6+z5-3*z4-z2)
256  & +f4*x2*(-6*z6-2*z5+3*z3+2*z2)
257  & +f4*x3*(4*z6+4*z5-4*z3-z2)
258  & +(f5-f2)*x2*(-2*z6+2*z4-z3)
259  & +f5*x3*(-4*z6-4*z4+4*z3+3*z2)
260  & +(2*f6-f1)*x2*(z5+3*z4-3*z3-z2)
261  & +2*(f6-f3)*x3*(2*z5-2*z4+z2))*xs144)**2
262 
263  diff = (f4+f5+f6)-(f1+f2+f3)
264  coef=xs24*surfac(ielem)
265 
266  w1(ielem)=w1(ielem)+((f4-f1+diff)*coef)**2
267  w2(ielem)=w2(ielem)+((f5-f2+diff)*coef)**2
268  w3(ielem)=w3(ielem)+((f6-f3+diff)*coef)**2
269  w4(ielem)=w4(ielem)+((f4-f1+diff)*coef)**2
270  w5(ielem)=w5(ielem)+((f5-f2+diff)*coef)**2
271  w6(ielem)=w6(ielem)+((f6-f3+diff)*coef)**2
272 
273  w1(ielem)=w1(ielem)/(coef*(sht+h1))
274  w2(ielem)=w2(ielem)/(coef*(sht+h2))
275  w3(ielem)=w3(ielem)/(coef*(sht+h3))
276  w4(ielem)=w4(ielem)/(coef*(sht+h1))
277  w5(ielem)=w5(ielem)/(coef*(sht+h2))
278  w6(ielem)=w6(ielem)/(coef*(sht+h3))
279 !
280  ENDDO ! IELEM
281 !
282 !-----------------------------------------------------------------------
283 !
284  ELSE
285 !
286  WRITE(lu,102) ielmf,sf%NAME
287 102 FORMAT(1x,'VC01PP (BIEF) :',/,
288  & 1x,'DISCRETISATION OF F : ',1i6,' NOT IMPLEMENTED',/,
289  & 1x,'REAL NAME OF F: ',a6)
290  CALL plante(1)
291  stop
292 !
293  ENDIF
294 !
295 !-----------------------------------------------------------------------
296 !
297  RETURN
298  END
declarations_special
Definition: declarations_special.F:3
vc21pp
subroutine vc21pp(XMUL, SF, F, X, Y, Z, SURFAC, IKLE1, IKLE2, IKLE3, IKLE4, IKLE5, IKLE6, NELEM, NELMAX, W1, W2, W3, W4, W5, W6)
Definition: vc21pp.f:9
declarations_special::lu
integer lu
Definition: declarations_special.F:45
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:58 for The TELEMAC-MASCARET system by doxygen 1.8.13