The TELEMAC-MASCARET system  trunk
vc11tt.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE vc11tt
3 ! *****************
4 !
5  &( xmul,sf,sg,f,g,x,y,z,ikle1,ikle2,ikle3,ikle4,nelem,nelmax,
6  & w1,w2,w3,w4,icoord )
7 !
8 !***********************************************************************
9 ! BIEF V6P1 21/08/2010
10 !***********************************************************************
11 !
12 !brief COMPUTES THE FOLLOWING VECTOR IN FINITE ELEMENTS:
13 !code
14 !+ (EXAMPLE OF THE X COMPONENT, WHICH CORRESPONDS TO ICOORD=1)
15 !+
16 !+ / DF
17 !+ VEC(I) = XMUL / ( G P *( -- )) D(OMEGA)
18 !+ /OMEGA I DX
19 !+
20 !+
21 !+ P IS A LINEAR BASE
22 !+ I
23 !+
24 !+ F IS A VECTOR OF TYPE P1 OR OTHER
25 !
26 !note IMPORTANT : IF F IS OF TYPE P0, THE RESULT IS 0.
27 !+
28 !+ HERE, IF F IS P0, IT REALLY MEANS THAT F IS
29 !+ P1, BUT GIVEN BY ELEMENTS.
30 !+
31 !+ THE SIZE OF F SHOULD THEN BE : F(NELMAX,4).
32 !
33 !warning THE JACOBIAN MUST BE POSITIVE
34 !warning THE RESULT IS IN W IN NOT ASSEMBLED FORM - REAL MESH
35 !
36 !history ARNAUD DESITTER - UNIVERSITY OF BRISTOL
37 !+ **/04/98
38 !+
39 !+
40 !
41 !history J-M HERVOUET (LNH)
42 !+ 25/03/02
43 !+ V5P3
44 !+
45 !
46 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
47 !+ 13/07/2010
48 !+ V6P0
49 !+ Translation of French comments within the FORTRAN sources into
50 !+ English comments
51 !
52 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
53 !+ 21/08/2010
54 !+ V6P0
55 !+ Creation of DOXYGEN tags for automated documentation and
56 !+ cross-referencing of the FORTRAN sources
57 !
58 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
59 !| F |-->| FUNCTION USED IN THE VECTOR FORMULA
60 !| G |-->| FUNCTION USED IN THE VECTOR FORMULA
61 !| ICOORD |-->| 1: DERIVATIVE ALONG X, 2: ALONG Y
62 !| IKLE1 |-->| FIRST POINT OF TETRAHEDRA
63 !| IKLE2 |-->| SECOND POINT OF TETRAHEDRA
64 !| IKLE3 |-->| THIRD POINT OF TETRAHEDRA
65 !| IKLE4 |-->| FOURTH POINT OF TETRAHEDRA
66 !| NELEM |-->| NUMBER OF ELEMENTS
67 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
68 !| SF |-->| BIEF_OBJ STRUCTURE OF F
69 !| SG |-->| BIEF_OBJ STRUCTURE OF G
70 !| SURFAC |-->| AREA OF TRIANGLES
71 !| W1 |<--| RESULT IN NON ASSEMBLED FORM
72 !| W2 |<--| RESULT IN NON ASSEMBLED FORM
73 !| W3 |<--| RESULT IN NON ASSEMBLED FORM
74 !| W4 |<--| 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
81 !
83  IMPLICIT NONE
84 !
85 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
86 !
87  INTEGER, INTENT(IN) :: NELEM,NELMAX,ICOORD
88  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax)
89  INTEGER, INTENT(IN) :: IKLE3(nelmax),IKLE4(nelmax)
90 !
91  DOUBLE PRECISION, INTENT(IN) :: X(*),Y(*),Z(*),XMUL
92  DOUBLE PRECISION, INTENT(INOUT) :: W1(nelmax),W2(nelmax)
93  DOUBLE PRECISION, INTENT(INOUT) :: W3(nelmax),W4(nelmax)
94 !
95 ! STRUCTURES OF F, G, H, U, V, W AND REAL DATA
96 !
97  TYPE(bief_obj), INTENT(IN) :: SF,SG
98  DOUBLE PRECISION, INTENT(IN) :: F(*),G(*)
99 !
100 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
101 !
102 ! LOCAL VARIABLES
103 !
104  INTEGER IELEM,IELMF,IELMG
105  DOUBLE PRECISION F1,F2,F3,F4
106  DOUBLE PRECISION G1,G2,G3,G4,X2,X3,X4,Y2,Y3,Y4,Z2,Z3,Z4
107  INTEGER I1,I2,I3,I4
108 !
109  DOUBLE PRECISION XSUR120,F2MF1,F3MF1,F4MF1,G2MG1,G3MG1,G4MG1
110 !
111 !-----------------------------------------------------------------------
112 ! INITIALISES
113 !
114  xsur120 = xmul/120.d0
115 !
116  ielmf = sf%ELM
117  ielmg = sg%ELM
118 !
119 !-----------------------------------------------------------------------
120 ! F AND G ARE LINEAR
121 !
122  IF ((ielmf.EQ.31.AND.((ielmg.EQ.31).OR.(ielmg.EQ.30))).OR.
123  & (ielmf.EQ.51.AND.ielmg.EQ.51) ) THEN
124 !
125  IF (icoord.EQ.1) THEN
126 !
127 !-----------------------------------------------------------------------
128 ! DERIVATIVE WRT X
129 !
130  DO ielem = 1 , nelem
131 !
132  i1 = ikle1(ielem)
133  i2 = ikle2(ielem)
134  i3 = ikle3(ielem)
135  i4 = ikle4(ielem)
136 !
137  f1 = f(i1)
138  f2 = f(i2)
139  f3 = f(i3)
140  f4 = f(i4)
141 !
142  IF (ielmg.EQ.31) THEN
143  g1 = g(i1)
144  g2 = g(i2)
145  g3 = g(i3)
146  g4 = g(i4)
147  ELSE
148  g1 = g(ielem)
149  g2 = g1
150  g3 = g1
151  g4 = g1
152  ENDIF
153 
154 !
155  f2mf1 = f2-f1
156  f3mf1 = f3-f1
157  f4mf1 = f4-f1
158  g2mg1 = g2-g1
159  g3mg1 = g3-g1
160  g4mg1 = g4-g1
161 !
162 ! REAL COORDINATES OF THE POINTS OF THE ELEMENT (ORIGIN IN 1)
163 !
164  y2 = y(i2) - y(i1)
165  y3 = y(i3) - y(i1)
166  y4 = y(i4) - y(i1)
167  z2 = z(i2) - z(i1)
168  z3 = z(i3) - z(i1)
169  z4 = z(i4) - z(i1)
170 !
171  w1(ielem) = (
172  & (5*f2mf1*g1+f2mf1*g2mg1+f2mf1*g3mg1+f2mf1*g4mg1)*(y3*z4-y4*z3)
173  &+(5*f3mf1*g1+f3mf1*g2mg1+f3mf1*g3mg1+f3mf1*g4mg1)*(z2*y4-y2*z4)
174  &+(5*f4mf1*g1+f4mf1*g2mg1+f4mf1*g3mg1+f4mf1*g4mg1)*(y2*z3-z2*y3)
175  & ) * xsur120
176 !
177  w2(ielem) = (
178  &-f4mf1*z2*y3*g4mg1+f4mf1*y2*z3*g4mg1+f3mf1*z2*y4*g4mg1
179  &-f3mf1*y2*z4*g4mg1+f2mf1*y3*z4*g4mg1-f2mf1*y4*z3*g4mg1
180  &+5*f2mf1*y3*z4*g1+2*f2mf1*y3*z4*g2mg1+f2mf1*y3*z4*g3mg1
181  &-5*f2mf1*y4*z3*g1-2*f2mf1*y4*z3*g2mg1-f2mf1*y4*z3*g3mg1
182  &-5*f3mf1*y2*z4*g1-2*f3mf1*y2*z4*g2mg1-f3mf1*y2*z4*g3mg1
183  &+5*f3mf1*z2*y4*g1+2*f3mf1*z2*y4*g2mg1+f3mf1*z2*y4*g3mg1
184  &+5*f4mf1*y2*z3*g1+2*f4mf1*y2*z3*g2mg1+f4mf1*y2*z3*g3mg1
185  &-5*f4mf1*z2*y3*g1-2*f4mf1*z2*y3*g2mg1-f4mf1*z2*y3*g3mg1
186  & ) * xsur120
187  w3(ielem) = (
188  &-(-f2mf1*y3*z4+f2mf1*y4*z3+f3mf1*y2*z4
189  & -f3mf1*z2*y4-f4mf1*y2*z3+f4mf1*z2*y3)
190  & *(2*g3mg1+g2mg1+g4mg1+5*g1)
191  & ) * xsur120
192  w4(ielem) = (
193  &-2*f4mf1*z2*y3*g4mg1
194  &+2*f4mf1*y2*z3*g4mg1
195  &+2*f3mf1*z2*y4*g4mg1
196  &-2*f3mf1*y2*z4*g4mg1
197  &+2*f2mf1*y3*z4*g4mg1
198  &-2*f2mf1*y4*z3*g4mg1
199  &+5*f2mf1*y3*z4*g1
200  &+f2mf1*y3*z4*g2mg1
201  &+f2mf1*y3*z4*g3mg1-5*f2mf1*y4*z3*g1-f2mf1*y4*z3*g2mg1
202  &-f2mf1*y4*z3*g3mg1-5*f3mf1*y2*z4*g1-f3mf1*y2*z4*g2mg1
203  &-f3mf1*y2*z4*g3mg1+5*f3mf1*z2*y4*g1+f3mf1*z2*y4*g2mg1
204  &+f3mf1*z2*y4*g3mg1+5*f4mf1*y2*z3*g1+f4mf1*y2*z3*g2mg1
205  &+f4mf1*y2*z3*g3mg1-5*f4mf1*z2*y3*g1-f4mf1*z2*y3*g2mg1
206  &-f4mf1*z2*y3*g3mg1
207  & ) * xsur120
208 !
209  ENDDO
210 !
211  ELSE IF (icoord.EQ.2) THEN
212 !
213 !-----------------------------------------------------------------------
214 ! DERIVATIVE WRT Y
215 !
216  DO ielem = 1 , nelem
217 !
218  i1 = ikle1(ielem)
219  i2 = ikle2(ielem)
220  i3 = ikle3(ielem)
221  i4 = ikle4(ielem)
222 !
223  f1 = f(i1)
224  f2 = f(i2)
225  f3 = f(i3)
226  f4 = f(i4)
227 !
228  IF (ielmg.EQ.31) THEN
229  g1 = g(i1)
230  g2 = g(i2)
231  g3 = g(i3)
232  g4 = g(i4)
233  ELSE
234  g1 = g(ielem)
235  g2 = g1
236  g3 = g1
237  g4 = g1
238  ENDIF
239 !
240  f2mf1 = f2-f1
241  f3mf1 = f3-f1
242  f4mf1 = f4-f1
243  g2mg1 = g2-g1
244  g3mg1 = g3-g1
245  g4mg1 = g4-g1
246 !
247 ! REAL COORDINATES OF THE POINTS OF THE ELEMENT (ORIGIN IN 1)
248 !
249  x2 = x(i2) - x(i1)
250  x3 = x(i3) - x(i1)
251  x4 = x(i4) - x(i1)
252  z2 = z(i2) - z(i1)
253  z3 = z(i3) - z(i1)
254  z4 = z(i4) - z(i1)
255 !
256  w1(ielem) = (
257  &-f2mf1*x3*z4*g2mg1+f3mf1*x2*z4*g3mg1+5*f3mf1*x2*z4*g1
258  &-f2mf1*x3*z4*g3mg1-5*f2mf1*x3*z4*g1
259  &+f2mf1*x4*z3*g2mg1+f2mf1*x4*z3*g3mg1+5*f2mf1*x4*z3*g1
260  &+f3mf1*x2*z4*g2mg1+f4mf1*z2*x3*g2mg1+f4mf1*z2*x3*g3mg1
261  &+5*f4mf1*z2*x3*g1-f4mf1*x2*z3*g2mg1-f4mf1*x2*z3*g3mg1
262  &-5*f4mf1*x2*z3*g1-f3mf1*z2*x4*g2mg1-f3mf1*z2*x4*g3mg1
263  &-5*f3mf1*z2*x4*g1+f2mf1*x4*z3*g4mg1+f3mf1*x2*z4*g4mg1
264  &-f4mf1*x2*z3*g4mg1-f3mf1*z2*x4*g4mg1-f2mf1*x3*z4*g4mg1
265  &+f4mf1*z2*x3*g4mg1 ) * xsur120
266  w2(ielem) = (
267  & -2*f2mf1*x3*z4*g2mg1+f3mf1*x2*z4*g3mg1+5*f3mf1*x2*z4*g1-f
268  &2mf1*x3*z4*g3mg1-5*f2mf1*x3*z4*g1+2*f2mf1*x4*z3*g2mg1+f2mf1*x4*z3*
269  &g3mg1+5*f2mf1*x4*z3*g1+2*f3mf1*x2*z4*g2mg1+2*f4mf1*z2*x3*g2mg1+f4m
270  &f1*z2*x3*g3mg1+5*f4mf1*z2*x3*g1-2*f4mf1*x2*z3*g2mg1-f4mf1*x2*z3*g3
271  &mg1-5*f4mf1*x2*z3*g1-2*f3mf1*z2*x4*g2mg1-f3mf1*z2*x4*g3mg1-5*f3mf1
272  &*z2*x4*g1+f2mf1*x4*z3*g4mg1+f3mf1*x2*z4*g4mg1-f4mf1*x2*z3*g4mg1-f3
273  &mf1*z2*x4*g4mg1-f2mf1*x3*z4*g4mg1+f4mf1*z2*x3*g4mg1 ) * xsur120
274  w3(ielem) = (
275  & -(f2mf1*x3*z4-f2mf1*x4*z3-f3mf1*x2*z4+f3mf1*z2*x4+f4mf1*x
276  &2*z3-f4mf1*z2*x3)*(2*g3mg1+g2mg1+g4mg1+5*g1) ) * xsur120
277  w4(ielem) = (
278  & -f2mf1*x3*z4*g2mg1+f3mf1*x2*z4*g3mg1+5*f3mf1*x2*z4*g1-f2m
279  &f1*x3*z4*g3mg1-5*f2mf1*x3*z4*g1+f2mf1*x4*z3*g2mg1+f2mf1*x4*z3*g3mg
280  &1+5*f2mf1*x4*z3*g1+f3mf1*x2*z4*g2mg1+f4mf1*z2*x3*g2mg1+f4mf1*z2*x3
281  &*g3mg1+5*f4mf1*z2*x3*g1-f4mf1*x2*z3*g2mg1-f4mf1*x2*z3*g3mg1-5*f4mf
282  &1*x2*z3*g1-f3mf1*z2*x4*g2mg1-f3mf1*z2*x4*g3mg1-5*f3mf1*z2*x4*g1+2*
283  &f2mf1*x4*z3*g4mg1+2*f3mf1*x2*z4*g4mg1-2*f4mf1*x2*z3*g4mg1-2*f3mf1*
284  &z2*x4*g4mg1-2*f2mf1*x3*z4*g4mg1+2*f4mf1*z2*x3*g4mg1 ) * xsur120
285 !
286  ENDDO
287 !
288  ELSE IF (icoord.EQ.3) THEN
289 !-----------------------------------------------------------------------
290 ! DERIVATIVE WRT Z
291 !
292  DO ielem = 1 , nelem
293 !
294  i1 = ikle1(ielem)
295  i2 = ikle2(ielem)
296  i3 = ikle3(ielem)
297  i4 = ikle4(ielem)
298 !
299  f1 = f(i1)
300  f2 = f(i2)
301  f3 = f(i3)
302  f4 = f(i4)
303 !
304  IF (ielmg.EQ.31) THEN
305  g1 = g(i1)
306  g2 = g(i2)
307  g3 = g(i3)
308  g4 = g(i4)
309  ELSE
310  g1 = g(ielem)
311  g2 = g1
312  g3 = g1
313  g4 = g1
314  ENDIF
315 !
316  f2mf1 = f2-f1
317  f3mf1 = f3-f1
318  f4mf1 = f4-f1
319  g2mg1 = g2-g1
320  g3mg1 = g3-g1
321  g4mg1 = g4-g1
322 !
323 ! REAL COORDINATES OF THE POINTS OF THE ELEMENT
324 !
325  x2 = x(i2) - x(i1)
326  x3 = x(i3) - x(i1)
327  x4 = x(i4) - x(i1)
328  y2 = y(i2) - y(i1)
329  y3 = y(i3) - y(i1)
330  y4 = y(i4) - y(i1)
331 !
332  w1(ielem) = (
333  & 5*f2mf1*x3*y4*g1+f2mf1*x3*y4*g2mg1+f2mf1*x3*y4*g3mg1-5*f2
334  &mf1*x4*y3*g1-f2mf1*x4*y3*g2mg1-f2mf1*x4*y3*g3mg1-5*f3mf1*x2*y4*g1-
335  &f3mf1*x2*y4*g2mg1-f3mf1*x2*y4*g3mg1+5*f3mf1*y2*x4*g1+f3mf1*y2*x4*g
336  &2mg1+f3mf1*y2*x4*g3mg1+5*f4mf1*x2*y3*g1+f4mf1*x2*y3*g2mg1-5*f4mf1*
337  &y2*x3*g1-f4mf1*y2*x3*g2mg1-f4mf1*y2*x3*g3mg1+f4mf1*x2*y3*g3mg1-f4m
338  &f1*y2*x3*g4mg1-f3mf1*x2*y4*g4mg1+f4mf1*x2*y3*g4mg1+f2mf1*x3*y4*g4m
339  &g1+f3mf1*y2*x4*g4mg1-f2mf1*x4*y3*g4mg1 ) * xsur120
340  w2(ielem) = (
341  & 5*f2mf1*x3*y4*g1+2*f2mf1*x3*y4*g2mg1+f2mf1*x3*y4*g3mg1-5*
342  &f2mf1*x4*y3*g1-2*f2mf1*x4*y3*g2mg1-f2mf1*x4*y3*g3mg1-5*f3mf1*x2*y4
343  &*g1-2*f3mf1*x2*y4*g2mg1-f3mf1*x2*y4*g3mg1+5*f3mf1*y2*x4*g1+2*f3mf1
344  &*y2*x4*g2mg1+f3mf1*y2*x4*g3mg1+5*f4mf1*x2*y3*g1+2*f4mf1*x2*y3*g2mg
345  &1-5*f4mf1*y2*x3*g1-2*f4mf1*y2*x3*g2mg1-f4mf1*y2*x3*g3mg1+f4mf1*x2*
346  &y3*g3mg1-f4mf1*y2*x3*g4mg1-f3mf1*x2*y4*g4mg1+f4mf1*x2*y3*g4mg1+f2m
347  &f1*x3*y4*g4mg1+f3mf1*y2*x4*g4mg1-f2mf1*x4*y3*g4mg1 ) * xsur120
348  w3(ielem) = (
349  & -(-f2mf1*x3*y4+f2mf1*x4*y3+f3mf1*x2*y4-f3mf1*y2*x4-f4mf1*
350  &x2*y3+f4mf1*y2*x3)*(2*g3mg1+g2mg1+g4mg1+5*g1) ) * xsur120
351  w4(ielem) = (
352  & 5*f2mf1*x3*y4*g1+f2mf1*x3*y4*g2mg1+f2mf1*x3*y4*g3mg1-5*f2
353  &mf1*x4*y3*g1-f2mf1*x4*y3*g2mg1-f2mf1*x4*y3*g3mg1-5*f3mf1*x2*y4*g1-
354  &f3mf1*x2*y4*g2mg1-f3mf1*x2*y4*g3mg1+5*f3mf1*y2*x4*g1+f3mf1*y2*x4*g
355  &2mg1+f3mf1*y2*x4*g3mg1+5*f4mf1*x2*y3*g1+f4mf1*x2*y3*g2mg1-5*f4mf1*
356  &y2*x3*g1-f4mf1*y2*x3*g2mg1-f4mf1*y2*x3*g3mg1+f4mf1*x2*y3*g3mg1-2*f
357  &4mf1*y2*x3*g4mg1-2*f3mf1*x2*y4*g4mg1+2*f4mf1*x2*y3*g4mg1+2*f2mf1*x
358  &3*y4*g4mg1+2*f3mf1*y2*x4*g4mg1-2*f2mf1*x4*y3*g4mg1 ) * xsur120
359 !
360  ENDDO
361 !
362  ELSE
363 !
364 !-----------------------------------------------------------------------
365 !
366  WRITE(lu,201) icoord
367  201 FORMAT(1x,'VC11TT (BIEF) : IMPOSSIBLE COMPONENT ',
368  & 1i6,' CHECK ICOORD')
369  CALL plante(1)
370  stop
371 !
372  ENDIF
373 !
374 !-----------------------------------------------------------------------
375 ! ERROR
376 !
377  ELSE
378 !
379 !-----------------------------------------------------------------------
380 !
381  WRITE(lu,1101) ielmf,sf%NAME
382  WRITE(lu,1201) ielmg,sg%NAME
383  WRITE(lu,1301)
384  CALL plante(1)
385  stop
386  1101 FORMAT(1x,'VC11TT (BIEF) :',/,
387  & 1x,'DISCRETIZATION OF F:',1i6,
388  & 1x,'REAL NAME: ',a6)
389  1201 FORMAT(1x,'DISCRETIZATION OF G:',1i6,
390  & 1x,'REAL NAME: ',a6)
391  1301 FORMAT(1x,'CASE NOT IMPLEMENTED')
392 !
393  ENDIF
394 !
395 !-----------------------------------------------------------------------
396 !
397  RETURN
398  END
subroutine vc11tt(XMUL, SF, SG, F, G, X, Y, Z, IKLE1, IKLE2, IKLE3, IKLE4, NELEM, NELMAX, W1, W2, W3, W4, ICOORD)
Definition: vc11tt.f:8
Definition: bief.f:3