The TELEMAC-MASCARET system  trunk
vc13pp.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE vc13pp
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,icoord,formul)
8 !
9 !***********************************************************************
10 ! BIEF V6P3 21/08/2010
11 !***********************************************************************
12 !
13 !brief COMPUTES THE FOLLOWING VECTOR IN FINITE ELEMENTS:
14 !code
15 !+ (EXAMPLE OF THE X COMPONENT, WHICH CORRESPONDS TO ICOORD=1)
16 !+
17 !+ / DF
18 !+ VEC(I) = XMUL / ( P *( -- )) D(OMEGA)
19 !+ /OMEGA I DX
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 !+ HERE, IF F IS P0, IT REALLY MEANS THAT F IS
28 !+ P1, BUT GIVEN BY ELEMENTS.
29 !+ THE SIZE OF F SHOULD THEN BE : F(NELMAX,3).
30 !
31 !warning THE JACOBIAN MUST BE POSITIVE
32 !warning THE RESULT IS IN W IN NOT ASSEMBLED FORM - REAL MESH
33 !
34 !history J-M HERVOUET (LNH)
35 !+ 19/05/2010
36 !+ V6P0
37 !+ NEW FILTER FOR PARTLY CRUSHED ELEMENTS
38 !
39 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
40 !+ 13/07/2010
41 !+ V6P0
42 !+ Translation of French comments within the FORTRAN sources into
43 !+ English comments
44 !
45 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
46 !+ 21/08/2010
47 !+ V6P0
48 !+ Creation of DOXYGEN tags for automated documentation and
49 !+ cross-referencing of the FORTRAN sources
50 !
51 !history J-M HERVOUET (EDF R&D LNHE)
52 !+ 07/01/2013
53 !+ V6P3
54 !+ X and Y are now given per element.
55 !
56 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
57 !| F |-->| FUNCTION USED IN THE VECTOR FORMULA
58 !| FORMUL |-->| SEE AT THE END OF THE SUBROUTINE
59 !| ICOORD |-->| 1: DERIVATIVE ALONG X, 2: ALONG Y
60 !| IKLE1 |-->| FIRST POINT OF PRISMS
61 !| IKLE2 |-->| SECOND POINT OF PRISMS
62 !| IKLE3 |-->| THIRD POINT OF PRISMS
63 !| IKLE4 |-->| FOURTH POINT OF PRISMS
64 !| IKLE5 |-->| FIFTH POINT OF PRISMS
65 !| IKLE6 |-->| SIXTH POINT OF PRISMS
66 !| NELEM |-->| NUMBER OF ELEMENTS
67 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
68 !| SF |-->| BIEF_OBJ STRUCTURE OF F
69 !| SURFAC |-->| AREA OF TRIANGLES
70 !| W1 |<--| RESULT IN NON ASSEMBLED FORM
71 !| W2 |<--| RESULT IN NON ASSEMBLED FORM
72 !| W3 |<--| RESULT IN NON ASSEMBLED FORM
73 !| W4 |<--| RESULT IN NON ASSEMBLED FORM
74 !| W5 |<--| RESULT IN NON ASSEMBLED FORM
75 !| W6 |<--| RESULT IN NON ASSEMBLED FORM
76 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
77 !| XMUL |-->| MULTIPLICATION COEFFICIENT
78 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
79 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
80 !
81  USE bief, ex_vc13pp => vc13pp
82 !
84  IMPLICIT NONE
85 !
86 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
87 !
88  INTEGER, INTENT(IN) :: NELEM,NELMAX,ICOORD
89  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax),IKLE3(nelmax)
90  INTEGER, INTENT(IN) :: IKLE4(nelmax),IKLE5(nelmax),IKLE6(nelmax)
91 ! NPOIN
92  DOUBLE PRECISION, INTENT(IN) :: X(nelmax,6),Y(nelmax,6),Z(*)
93  DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
94  DOUBLE PRECISION, INTENT(INOUT) ::W1(nelmax),W2(nelmax),W3(nelmax)
95  DOUBLE PRECISION, INTENT(INOUT) ::W4(nelmax),W5(nelmax),W6(nelmax)
96  DOUBLE PRECISION, INTENT(IN) :: XMUL
97 !
98 ! STRUCTURE OF F AND REAL DATA
99 !
100  TYPE(bief_obj), INTENT(IN) :: SF
101  DOUBLE PRECISION, INTENT(IN) :: F(*)
102  CHARACTER(LEN=16), INTENT(IN) :: FORMUL
103 !
104 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
105 !
106  DOUBLE PRECISION XS24,XS144,F1,F2,F3,F4,F5,F6,XMU
107  DOUBLE PRECISION X2,X3,Y2,Y3,Z1,Z2,Z3,Z4,Z5,Z6
108  INTEGER I1,I2,I3,I4,I5,I6,IELEM,IELMF
109 !
110  INTRINSIC max,min
111 !
112 !-----------------------------------------------------------------------
113 !
114  xs24 = xmul/24.d0
115  xs144 = xmul/144.d0
116 !
117 !-----------------------------------------------------------------------
118 !
119  ielmf=sf%ELM
120 !
121 !=======================================================================
122 !
123 ! F IS LINEAR
124 !
125  IF(ielmf.EQ.41) THEN
126 !
127  IF(icoord.EQ.1) THEN
128 !
129 !-----------------------------------------------------------------------
130 !
131 ! DERIVATIVE WRT X
132 !
133  DO ielem = 1 , nelem
134 !
135  i1 = ikle1(ielem)
136  i2 = ikle2(ielem)
137  i3 = ikle3(ielem)
138  i4 = ikle4(ielem)
139  i5 = ikle5(ielem)
140  i6 = ikle6(ielem)
141 !
142  f1 = f(i1)
143  f2 = f(i2)
144  f3 = f(i3)
145  f4 = f(i4)
146  f5 = f(i5)
147  f6 = f(i6)
148 !
149 ! REAL COORDINATES OF THE POINTS OF THE ELEMENT (ORIGIN IN 1)
150 !
151 ! Y2 = Y(I2) - Y(I1)
152 ! Y3 = Y(I3) - Y(I1)
153  y2 = y(ielem,2)
154  y3 = y(ielem,3)
155 !
156  z2 = z(i2) - z(i1)
157  z3 = z(i3) - z(i1)
158  z4 = z(i4) - z(i1)
159  z5 = z(i5) - z(i1)
160  z6 = z(i6) - z(i1)
161 !
162  w1(ielem)=( (2*f1-f6)*y2*( z5+3*z4-3*z3 -z2)
163  & +(2*f1-f5)*y3*(-z6-3*z4 +z3+3*z2)
164  & +f2*y3*(2*z6+3*z5+3*z4-2*z3)
165  & +f3*y2*(-3*z6-2*z5-3*z4+2*z2)
166  & +(f3-f6)*y3*(z5-z4+2*z2)
167  & +f4*y2*(3*z6+z5+3*z3-z2)
168  & +f4*y3*(-z6-3*z5+z3-3*z2)
169  & +(f5-f2)*y2*(z6-z4+2*z3) )*xs144
170  w2(ielem)=( f1*y2*(z6+4*z5+3*z4-4*z3-4*z2)
171  & +f1*y3*(-2*z6-3*z5-3*z4+2*z3+6*z2)
172  & +(2*f2-f4)*y3*(z6+3*z5-z3)
173  & +f3*y2*(-3*z6-4*z5-z4+4*z2)
174  & +f4*y2*(2*z6+2*z5+z3-2*z2)
175  & +2*(f5-f2)*y2*(z6-z4+2*z3)
176  & +f5*y3*(z6+3*z4-z3-6*z2)
177  & +f6*y2*(-2*z5-2*z4+3*z3+2*z2)
178  & +(f6-f3)*y3*(-z5+z4-2*z2) )*xs144
179  w3(ielem)=( f1*y2*(3*z6+2*z5+3*z4-6*z3-2*z2)
180  & +f1*y3*(-4*z6-z5-3*z4+4*z3+4*z2)
181  & +f2*y3*(4*z6+3*z5+z4-4*z3)
182  & +(2*f3-f4)*y2*(-3*z6-z5+z2)
183  & +f4*y3*(-2*z6-2*z5+2*z3-z2)
184  & +(f5-f2)*y2*(z6-z4+2*z3)
185  & +f5*y3*(2*z6+2*z4-2*z3-3*z2)
186  & +f6*y2*(-z5-3*z4+6*z3+z2)
187  & +2*(f6-f3)*y3*(-z5+z4-2*z2) )*xs144
188  w4(ielem)=( f1*y2*(-3*z6+z5+6*z4-3*z3-z2)
189  & +f1*y3*(-z6+3*z5-6*z4+z3+3*z2)
190  & +f2*y3*(z6+3*z5-z3)
191  & +(2*f4-f3)*y2*(3*z6+z5-z2)
192  & +2*f4*y3*(-z6-3*z5+z3)
193  & +(f5-f2)*y2*(2*z6-2*z4+z3)
194  & +f5*y3*(2*z6+6*z4-2*z3-3*z2)
195  & +f6*y2*(-2*z5-6*z4+3*z3+2*z2)
196  & +(f6-f3)*y3*(-2*z5+2*z4-z2) )*xs144
197  w5(ielem)=( f1*y2*(-z6+2*z5+3*z4-2*z3-2*z2)
198  & +f2*y3*(z6+6*z5-3*z4-z3)
199  & +f3*y2*(-3*z6-2*z5+z4+2*z2)
200  & +f4*y2*(4*z6+4*z5-z3-4*z2)
201  & +f4*y3*(-2*z6-6*z5+2*z3+3*z2)
202  & +2*(f5-f2)*y2*(2*z6-2*z4+z3)
203  & +(2*f5-f1)*y3*(z6+3*z4-z3-3*z2)
204  & +f6*y2*(-4*z5-4*z4+3*z3+4*z2)
205  & +(f6-f3)*y3*(-2*z5+2*z4-z2) )*xs144
206  w6(ielem)=( +f1*y3*(-2*z6+z5-3*z4+2*z3+2*z2)
207  & +f2*y3*(2*z6+3*z5-z4-2*z3)
208  & +f3*y2*(-6*z6-z5+3*z4+z2)
209  & +f4*y2*(6*z6+2*z5-3*z3-2*z2)
210  & +f4*y3*(-4*z6-4*z5+4*z3+z2)
211  & +(f5-f2)*y2*(2*z6-2*z4+z3)
212  & +f5*y3*(4*z6+4*z4-4*z3-3*z2)
213  & +(2*f6-f1)*y2*(-z5-3*z4+3*z3+z2)
214  & +2*(f6-f3)*y3*(-2*z5+2*z4-z2) )*xs144
215 !
216  ENDDO ! IELEM
217 !
218  ELSEIF(icoord.EQ.2) THEN
219 !
220 !-----------------------------------------------------------------------
221 !
222 ! DERIVATIVE WRT Y
223 !
224  DO ielem = 1 , nelem
225 !
226  i1 = ikle1(ielem)
227  i2 = ikle2(ielem)
228  i3 = ikle3(ielem)
229  i4 = ikle4(ielem)
230  i5 = ikle5(ielem)
231  i6 = ikle6(ielem)
232 !
233  f1 = f(i1)
234  f2 = f(i2)
235  f3 = f(i3)
236  f4 = f(i4)
237  f5 = f(i5)
238  f6 = f(i6)
239 !
240 ! REAL COORDINATES OF THE POINTS OF THE ELEMENT (ORIGIN IN 1)
241 !
242 ! X2 = X(I2) - X(I1)
243 ! X3 = X(I3) - X(I1)
244  x2 = x(ielem,2)
245  x3 = x(ielem,3)
246 !
247  z2 = z(i2) - z(i1)
248  z3 = z(i3) - z(i1)
249  z4 = z(i4) - z(i1)
250  z5 = z(i5) - z(i1)
251  z6 = z(i6) - z(i1)
252 !
253  w1(ielem)=( (2*f1-f6)*x2*(-z5-3*z4+3*z3+z2)
254  & +2*f1*x3*(z6+3*z4-z3-3*z2)
255  & +f2*x3*(-2*z6-3*z5-3*z4+2*z3)
256  & +f3*x2*(3*z6+2*z5+3*z4-2*z2)
257  & +f4*x2*(-3*z6-z5-3*z3+z2)
258  & +f4*x3*(z6+3*z5-z3+3*z2)
259  & +(f5-f2)*x2*(-z6+z4-2*z3)
260  & +f5*x3*(-z6-3*z4+z3+3*z2)
261  & +(f6-f3)*x3*(z5-z4+2*z2) )*xs144
262  w2(ielem)=( f1*x2*(-z6-4*z5-3*z4+4*z3+4*z2)
263  & +f1*x3*(2*z6+3*z5+3*z4-2*z3-6*z2)
264  & +(2*f2-f4)*x3*(-z6-3*z5+z3)
265  & +f3*x2*(3*z6+4*z5+z4-4*z2)
266  & +f4*x2*(-2*z6-2*z5-z3+2*z2)
267  & +2*(f5-f2)*x2*(-z6+z4-2*z3)
268  & +f5*x3*(-z6-3*z4+z3+6*z2)
269  & +f6*x2*(2*z5+2*z4-3*z3-2*z2)
270  & +(f6-f3)*x3*(z5-z4+2*z2) )*xs144
271  w3(ielem)=( f1*x2*(-3*z6-2*z5-3*z4+6*z3+2*z2)
272  & +f1*x3*(4*z6+z5+3*z4-4*z3-4*z2)
273  & +f2*x3*(-4*z6-3*z5-z4+4*z3)
274  & +(2*f3-f4)*x2*(3*z6+z5-z2)
275  & +f4*x3*(2*z6+2*z5-2*z3+z2)
276  & +(f5-f2)*x2*(-z6+z4-2*z3)
277  & +f5*x3*(-2*z6-2*z4+2*z3+3*z2)
278  & +f6*x2*(z5+3*z4-6*z3-z2)
279  & +2*(f6-f3)*x3*(z5-z4+2*z2) )*xs144
280  w4(ielem)=( f1*x2*(3*z6-z5-6*z4+3*z3+z2)
281  & +f1*x3*(z6-3*z5+6*z4-z3-3*z2)
282  & +(2*f4-f3)*x2*(-3*z6-z5+z2)
283  & +(2*f4-f2)*x3*(z6+3*z5-z3)
284  & +(f5-f2)*x2*(-2*z6+2*z4-z3)
285  & +f5*x3*(-2*z6-6*z4+2*z3+3*z2)
286  & +f6*x2*(2*z5+6*z4-3*z3-2*z2)
287  & +(f6-f3)*x3*(2*z5-2*z4+z2) )*xs144
288  w5(ielem)=( f1*x2*(z6-2*z5-3*z4+2*z3+2*z2)
289  & +f2*x3*(-z6-6*z5+3*z4+z3)
290  & +f3*x2*(3*z6+2*z5-z4-2*z2)
291  & +f4*x2*(-4*z6-4*z5+z3+4*z2)
292  & +f4*x3*(2*z6+6*z5-2*z3-3*z2)
293  & +2*(f5-f2)*x2*(-2*z6+2*z4-z3)
294  & +(2*f5-f1)*x3*(-z6-3*z4+z3+3*z2)
295  & +f6*x2*(4*z5+4*z4-3*z3-4*z2)
296  & +(f6-f3)*x3*(2*z5-2*z4+z2) )*xs144
297  w6(ielem)=(+f1*x3*(2*z6-z5+3*z4-2*z3-2*z2)
298  & +f2*x3*(-2*z6-3*z5+z4+2*z3)
299  & +f3*x2*(6*z6+z5-3*z4-z2)
300  & +f4*x2*(-6*z6-2*z5+3*z3+2*z2)
301  & +f4*x3*(4*z6+4*z5-4*z3-z2)
302  & +(f5-f2)*x2*(-2*z6+2*z4-z3)
303  & +f5*x3*(-4*z6-4*z4+4*z3+3*z2)
304  & +(2*f6-f1)*x2*(z5+3*z4-3*z3-z2)
305  & +2*(f6-f3)*x3*(2*z5-2*z4+z2) )*xs144
306 !
307 !
308  ENDDO ! IELEM
309 !
310  ELSEIF(icoord.EQ.3) THEN
311 !
312 !-----------------------------------------------------------------------
313 !
314 ! DERIVATIVE WRT Z
315 !
316  DO ielem = 1 , nelem
317 !
318  i1 = ikle1(ielem)
319  i2 = ikle2(ielem)
320  i3 = ikle3(ielem)
321  i4 = ikle4(ielem)
322  i5 = ikle5(ielem)
323  i6 = ikle6(ielem)
324 !
325  f1 = f(i1)
326  f2 = f(i2)
327  f3 = f(i3)
328  f4 = f(i4)
329  f5 = f(i5)
330  f6 = f(i6)
331 !
332 ! REAL COORDINATES OF THE POINTS OF THE ELEMENT (ORIGIN IN 1)
333 !
334 ! X2 = X(I2) - X(I1)
335 ! X3 = X(I3) - X(I1)
336 ! Y2 = Y(I2) - Y(I1)
337 ! Y3 = Y(I3) - Y(I1)
338 !
339 ! XMU = XS48*(X2*Y3-X3*Y2)
340  xmu = xs24*surfac(ielem)
341 !
342 ! NOT LUMPED VERSION
343 ! DIFF = (F4+F5+F6) - (F1+F2+F3)
344 ! W1(IELEM)=(F4-F1+DIFF)*XMU
345 ! W2(IELEM)=(F5-F2+DIFF)*XMU
346 ! W3(IELEM)=(F6-F3+DIFF)*XMU
347 ! LUMPED VERSION (LIKE THE DIFFUSION MATRIX)
348 ! SEE W COMPUTATION IN TELEMAC-3D IN PROVEL
349  w1(ielem)=4*(f4-f1)*xmu
350  w2(ielem)=4*(f5-f2)*xmu
351  w3(ielem)=4*(f6-f3)*xmu
352 !
353  w4(ielem)=w1(ielem)
354  w5(ielem)=w2(ielem)
355  w6(ielem)=w3(ielem)
356 !
357  ENDDO ! IELEM
358 !
359  ELSE
360 !
361 !-----------------------------------------------------------------------
362 !
363  WRITE(lu,201) icoord
364 201 FORMAT(1x,'VC13PP (BIEF) : IMPOSSIBLE COMPONENT ',
365  & 1i6,' CHECK ICOORD')
366  CALL plante(1)
367  stop
368 !
369  ENDIF
370 !
371 !=======================================================================
372 !
373  ELSE
374 !
375 !=======================================================================
376 !
377  WRITE(lu,102) ielmf,sf%NAME
378 102 FORMAT(1x,'VC13PP (BIEF) :',/,
379  & 1x,'DISCRETISATION OF F : ',1i6,' NOT IMPLEMENTED',/,
380  & 1x,'REAL NAME OF F: ',a6)
381  CALL plante(1)
382  stop
383 !
384  ENDIF
385 !
386 !=======================================================================
387 !
388 ! HYDROSTATIC INCONSISTENCIES
389 !
390 ! COMMON TREATMENT FOR FILTERS 2,3 AND 4
391 !
392  IF(formul(6:6).EQ.'2'.OR.formul(6:6).EQ.'3'.OR.
393  & formul(6:6).EQ.'4' ) THEN
394 !
395  DO ielem = 1 , nelem
396 !
397  i1 = ikle1(ielem)
398  i2 = ikle2(ielem)
399  i3 = ikle3(ielem)
400  i4 = ikle4(ielem)
401  i5 = ikle5(ielem)
402  i6 = ikle6(ielem)
403 !
404  IF(max(z(i1),z(i2),z(i3)).GT.min(z(i4),z(i5),z(i6))) THEN
405  w1(ielem)=0.d0
406  w2(ielem)=0.d0
407  w3(ielem)=0.d0
408  w4(ielem)=0.d0
409  w5(ielem)=0.d0
410  w6(ielem)=0.d0
411  ENDIF
412 !
413  ENDDO
414 !
415  ENDIF
416 !
417 ! FILTER 3
418 !
419  IF(formul(6:6).EQ.'3'.AND.(icoord.EQ.1.OR.icoord.EQ.2)) THEN
420 !
421  DO ielem = 1 , nelem
422 !
423  i1 = ikle1(ielem)
424  i2 = ikle2(ielem)
425  i3 = ikle3(ielem)
426  i4 = ikle4(ielem)
427  i5 = ikle5(ielem)
428  i6 = ikle6(ielem)
429 !
430  f1 = f(i1)
431  f2 = f(i2)
432  f3 = f(i3)
433  f4 = f(i4)
434  f5 = f(i5)
435  f6 = f(i6)
436 !
437 ! IF THERE IS A POSSIBILITY OF STRATIFICATION
438 ! GRADIENTS CANCELLED
439 !
440  IF( min(max(f1,f4),max(f2,f5),max(f3,f6)).GE.
441  & max(min(f1,f4),min(f2,f5),min(f3,f6)) ) THEN
442  w1(ielem)=0.d0
443  w2(ielem)=0.d0
444  w3(ielem)=0.d0
445  w4(ielem)=0.d0
446  w5(ielem)=0.d0
447  w6(ielem)=0.d0
448  ENDIF
449 !
450  ENDDO
451 !
452  ENDIF
453 !
454 ! FILTER 4
455 !
456  IF(formul(6:6).EQ.'4'.AND.(icoord.EQ.1.OR.icoord.EQ.2)) THEN
457 !
458  DO ielem = 1 , nelem
459 !
460  i1 = ikle1(ielem)
461  i2 = ikle2(ielem)
462  i3 = ikle3(ielem)
463  i4 = ikle4(ielem)
464  i5 = ikle5(ielem)
465  i6 = ikle6(ielem)
466 !
467  f1 = f(i1)
468  f2 = f(i2)
469  f3 = f(i3)
470  f4 = f(i4)
471  f5 = f(i5)
472  f6 = f(i6)
473 !
474  z1 = z(i1)
475  z2 = z(i2)
476  z3 = z(i3)
477  z4 = z(i4)
478  z5 = z(i5)
479  z6 = z(i6)
480 !
481 ! CHECKS IF A STRATIFICATION IS IMPOSSIBLE
482 ! IN THIS CASE (GO TO 1000) GRADIENTS ARE KEPT
483 !
484 ! 1 IN BETWEEN 3 AND 6
485  IF(z1.GE.z3.AND.z1.LE.z6) THEN
486  IF(f1.LT.min(f3,f6).OR.f1.GT.max(f3,f6)) GO TO 1000
487  ENDIF
488 ! 1 IN BETWEEN 2 AND 5
489  IF(z1.GE.z2.AND.z1.LE.z5) THEN
490  IF(f1.LT.min(f2,f5).OR.f1.GT.max(f2,f5)) GO TO 1000
491  ENDIF
492 ! 2 IN BETWEEN 1 AND 4
493  IF(z2.GE.z1.AND.z2.LE.z4) THEN
494  IF(f2.LT.min(f1,f4).OR.f2.GT.max(f1,f4)) GO TO 1000
495  ENDIF
496 ! 2 IN BETWEEN 3 AND 6
497  IF(z2.GE.z3.AND.z2.LE.z6) THEN
498  IF(f2.LT.min(f3,f6).OR.f2.GT.max(f3,f6)) GO TO 1000
499  ENDIF
500 ! 3 IN BETWEEN 1 AND 4
501  IF(z3.GE.z1.AND.z3.LE.z4) THEN
502  IF(f3.LT.min(f1,f4).OR.f3.GT.max(f1,f4)) GO TO 1000
503  ENDIF
504 ! 3 IN BETWEEN 2 AND 5
505  IF(z3.GE.z2.AND.z3.LE.z5) THEN
506  IF(f3.LT.min(f2,f5).OR.f3.GT.max(f2,f5)) GO TO 1000
507  ENDIF
508 ! 4 IN BETWEEN 2 AND 5
509  IF(z4.GE.z2.AND.z4.LE.z5) THEN
510  IF(f4.LT.min(f2,f5).OR.f4.GT.max(f2,f5)) GO TO 1000
511  ENDIF
512 ! 4 IN BETWEEN 3 AND 6
513  IF(z4.GE.z3.AND.z4.LE.z6) THEN
514  IF(f4.LT.min(f3,f6).OR.f4.GT.max(f3,f6)) GO TO 1000
515  ENDIF
516 ! 5 IN BETWEEN 1 AND 4
517  IF(z5.GE.z1.AND.z5.LE.z4) THEN
518  IF(f5.LT.min(f1,f4).OR.f5.GT.max(f1,f4)) GO TO 1000
519  ENDIF
520 ! 5 IN BETWEEN 3 AND 6
521  IF(z5.GE.z3.AND.z5.LE.z6) THEN
522  IF(f5.LT.min(f3,f6).OR.f5.GT.max(f3,f6)) GO TO 1000
523  ENDIF
524 ! 6 IN BETWEEN 1 AND 4
525  IF(z6.GE.z1.AND.z6.LE.z4) THEN
526  IF(f6.LT.min(f1,f4).OR.f6.GT.max(f1,f4)) GO TO 1000
527  ENDIF
528 ! 6 IN BETWEEN 2 AND 5
529  IF(z6.GE.z2.AND.z6.LE.z5) THEN
530  IF(f6.LT.min(f2,f5).OR.f6.GT.max(f2,f5)) GO TO 1000
531  ENDIF
532 !
533 ! SO THERE IS A POSSIBILITY OF STRATIFICATION
534 ! GRADIENTS CANCELLED
535 !
536  w1(ielem)=0.d0
537  w2(ielem)=0.d0
538  w3(ielem)=0.d0
539  w4(ielem)=0.d0
540  w5(ielem)=0.d0
541  w6(ielem)=0.d0
542 !
543 1000 CONTINUE
544 !
545  ENDDO
546 !
547  ENDIF
548 !
549 ! FILTER FOR PARTLY CRUSHED PRISMS
550 !
551  IF(formul(7:7).EQ.'2') THEN
552 !
553  DO ielem = 1 , nelem
554  i1 = ikle1(ielem)
555  i2 = ikle2(ielem)
556  i3 = ikle3(ielem)
557  i4 = ikle4(ielem)
558  i5 = ikle5(ielem)
559  i6 = ikle6(ielem)
560  IF(z(i4)-z(i1).LT.1.d-3.OR.
561  & z(i5)-z(i2).LT.1.d-3.OR.
562  & z(i6)-z(i3).LT.1.d-3 ) THEN
563  w1(ielem)=0.d0
564  w2(ielem)=0.d0
565  w3(ielem)=0.d0
566  w4(ielem)=0.d0
567  w5(ielem)=0.d0
568  w6(ielem)=0.d0
569  ENDIF
570  ENDDO
571 !
572  ENDIF
573 !
574 !=======================================================================
575 !
576  RETURN
577  END
subroutine vc13pp(XMUL, SF, F, X, Y, Z, SURFAC, IKLE1, IKLE2, IKLE3, IKLE4, IKLE5, IKLE6, NELEM, NELMAX, W1, W2, W3, W4, W5, W6, ICOORD, FORMUL)
Definition: vc13pp.f:9
Definition: bief.f:3