vc08bb.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\vc08bb.f
00002 !
00068                      SUBROUTINE VC08BB
00069 !                    *****************
00070 !
00071      &(XMUL,SF,SU,SV,F,U,V,XEL,YEL,
00072      & IKLE1,IKLE2,IKLE3,IKLE4,NELEM,NELMAX,W1,W2,W3,W4,FORMUL)
00073 !
00074 !***********************************************************************
00075 ! BIEF   V6P1                                   21/08/2010
00076 !***********************************************************************
00077 !
00078 !
00079 !
00080 !
00081 !
00082 !
00083 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00084 !| F              |-->| FUNCTION USED IN THE VECTOR FORMULA
00085 !| FORMUL         |-->| STRING WITH FORMULA OF VECTOR
00086 !| IKLE1          |-->| FIRST POINT OF TRIANGLES
00087 !| IKLE2          |-->| SECOND POINT OF TRIANGLES
00088 !| IKLE3          |-->| THIRD POINT OF TRIANGLES
00089 !| IKLE4          |-->| QUASI-BUBBLE POINT OF TRIANGLES
00090 !| NELEM          |-->| NUMBER OF ELEMENTS
00091 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
00092 !| SF             |-->| BIEF_OBJ STRUCTURE OF F
00093 !| SU             |-->| BIEF_OBJ STRUCTURE OF U
00094 !| SV             |-->| BIEF_OBJ STRUCTURE OF V
00095 !| U              |-->| FUNCTION USED IN THE VECTOR FORMULA
00096 !| V              |-->| FUNCTION USED IN THE VECTOR FORMULA
00097 !| W1             |<--| RESULT IN NON ASSEMBLED FORM
00098 !| W2             |<--| RESULT IN NON ASSEMBLED FORM
00099 !| W3             |<--| RESULT IN NON ASSEMBLED FORM
00100 !| W4             |<--| RESULT IN NON ASSEMBLED FORM
00101 !| XEL            |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
00102 !| XMUL           |-->| MULTIPLICATION COEFFICIENT
00103 !| YEL            |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
00104 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00105 !
00106       USE BIEF, EX_VC08BB => VC08BB
00107 !
00108       IMPLICIT NONE
00109       INTEGER LNG,LU
00110       COMMON/INFO/LNG,LU
00111 !
00112 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00113 !
00114       INTEGER, INTENT(IN) :: NELEM,NELMAX
00115       DOUBLE PRECISION, INTENT(IN) :: XEL(NELMAX*3),YEL(NELMAX*3),XMUL
00116 !     W1 IS ALSO USED AS 1-DIMENSIONAL FOR ALL W
00117       DOUBLE PRECISION, INTENT(INOUT) :: W1(4*NELMAX),W2(NELMAX)
00118       DOUBLE PRECISION, INTENT(INOUT) :: W3(NELMAX),W4(NELMAX)
00119 !     IKLE1 IS ALSO USED AS A 1-DIMENSIONAL IKLE
00120       INTEGER, INTENT(IN) :: IKLE1(4*NELMAX)
00121       INTEGER, INTENT(IN) :: IKLE2(NELMAX),IKLE3(NELMAX),IKLE4(NELMAX)
00122       CHARACTER(LEN=16), INTENT(IN) :: FORMUL
00123 !
00124 !     STRUCTURES OF F, G, H, U, V, W AND REAL DATA
00125 !
00126       TYPE(BIEF_OBJ), INTENT(IN) :: SF,SU,SV
00127       DOUBLE PRECISION, INTENT(IN) :: F(*),U(*),V(*)
00128 !
00129 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00130 !
00131       INTEGER IELEM,IELMF,IELMU,IELMV,IL(3,3)
00132       INTEGER IG1,IG2,IG3,IT,IAD1,IAD2,IAD3
00133 !
00134       DOUBLE PRECISION K1,K2,K3
00135 !
00136 !-----------------------------------------------------------------------
00137 !
00138       DOUBLE PRECISION X2,Y2,X3,Y3,F1,F2,F3,F4,U1,U2,U3,U4,V1,V2,V3,V4
00139       DOUBLE PRECISION XSUR72,XSU216,X1,Y1,TIERS
00140       DOUBLE PRECISION PHIT,SUR6,USUR2,VSUR2
00141       DOUBLE PRECISION L12,L13,L21,L23,L31,L32,SIG,BETAN1,BETAN2,BETAN3
00142 !
00143       INTRINSIC MAX,MIN,SIGN
00144 !
00145 !-----------------------------------------------------------------------
00146 !
00147 !     FOR A QUASI-BUBBLE TRIANGLE : NODE NUMBERS OF THE SUB-TRIANGLES
00148 !     IN THE INITIAL TRIANGLE
00149 !     IL(NUMBER OF THE SUB-TRIANGLE,LOCAL NUMBER IN THE SUB-TRIANGLE)
00150 !
00151       DATA IL /1,2,3,2,3,1,4,4,4/
00152 !
00153 !-----------------------------------------------------------------------
00154 !
00155       XSUR72 = XMUL/72.D0
00156       XSU216 = XMUL/216.D0
00157       TIERS  = 1.D0 / 3.D0
00158       SUR6   = 1.D0 / 6.D0
00159 !
00160       IELMF=SF%ELM
00161       IELMU=SU%ELM
00162       IELMV=SV%ELM
00163 !
00164 !-----------------------------------------------------------------------
00165 !
00166 !     FUNCTION F AND VECTOR U ARE QUASI-BUBBLE
00167 !
00168       IF(IELMF.EQ.12.AND.IELMU.EQ.12.AND.IELMV.EQ.12) THEN
00169 !
00170       IF(FORMUL(14:16).EQ.'PSI') THEN
00171 !
00172 !     INITIALISES W
00173 !
00174       DO IELEM = 1 , NELEM
00175         W1(IELEM) = 0.D0
00176         W2(IELEM) = 0.D0
00177         W3(IELEM) = 0.D0
00178         W4(IELEM) = 0.D0
00179       ENDDO ! IELEM
00180 !
00181 !     PSI SCHEME, LOOP ON THE 3 SUB-TRIANGLES AND
00182 !     PREASSEMBLY
00183 !
00184       DO IT=1,3
00185       DO IELEM = 1 , NELEM
00186 !
00187 !     ADDRESSES IN ARRAY (NELMAX,*)
00188       IAD1= IELEM + (IL(IT,1)-1)*NELMAX
00189       IAD2= IELEM + (IL(IT,2)-1)*NELMAX
00190       IAD3= IELEM + (IL(IT,3)-1)*NELMAX
00191 !     GLOBAL NUMBERS IN THE INITIAL TRIANGLE
00192       IG1 = IKLE1(IAD1)
00193       IG2 = IKLE1(IAD2)
00194       IG3 = IKLE1(IAD3)
00195 !     COORDINATES OF THE SUB-TRIANGLE NODES
00196       X1 = XEL(IAD1)
00197       X2 = XEL(IAD2) - X1
00198 !     POINT 3 IS ALWAYS THE CENTRE OF THE INITIAL TRIANGLE
00199       X3=TIERS*(XEL(IELEM)+XEL(IELEM+NELMAX)+XEL(IELEM+2*NELMAX))-X1
00200       Y1 = YEL(IAD1)
00201       Y2 = YEL(IAD2) - Y1
00202 !     POINT 3 IS ALWAYS THE CENTRE OF THE INITIAL TRIANGLE
00203       Y3=TIERS*(YEL(IELEM)+YEL(IELEM+NELMAX)+YEL(IELEM+2*NELMAX))-Y1
00204 !     F VALUES IN THE SUB-TRIANGLE
00205       F1 = F(IG1)
00206       F2 = F(IG2)
00207       F3 = F(IG3)
00208 !
00209       USUR2 = (U(IG1)+U(IG2)+U(IG3))*SUR6
00210       VSUR2 = (V(IG1)+V(IG2)+V(IG3))*SUR6
00211 !
00212       K1 = USUR2 * (Y2-Y3) - VSUR2 * (X2-X3)
00213       K2 = USUR2 * (Y3   ) - VSUR2 * (X3   )
00214       K3 = USUR2 * (  -Y2) - VSUR2 * (  -X2)
00215 !
00216       L12 = MAX(  MIN(K1,-K2) , 0.D0 )
00217       L13 = MAX(  MIN(K1,-K3) , 0.D0 )
00218       L21 = MAX(  MIN(K2,-K1) , 0.D0 )
00219       L23 = MAX(  MIN(K2,-K3) , 0.D0 )
00220       L31 = MAX(  MIN(K3,-K1) , 0.D0 )
00221       L32 = MAX(  MIN(K3,-K2) , 0.D0 )
00222 !
00223       BETAN1 = L12*(F1-F2) + L13*(F1-F3)
00224       BETAN2 = L21*(F2-F1) + L23*(F2-F3)
00225       BETAN3 = L31*(F3-F1) + L32*(F3-F2)
00226 !
00227       PHIT = BETAN1 + BETAN2 + BETAN3
00228       SIG  = SIGN(1.D0,PHIT)
00229 !
00230       W1(IAD1)=W1(IAD1)+ XMUL * SIG * MAX(MIN(SIG*BETAN1,SIG*PHIT),0.D0)
00231       W1(IAD2)=W1(IAD2)+ XMUL * SIG * MAX(MIN(SIG*BETAN2,SIG*PHIT),0.D0)
00232       W1(IAD3)=W1(IAD3)+ XMUL * SIG * MAX(MIN(SIG*BETAN3,SIG*PHIT),0.D0)
00233 !
00234       ENDDO ! IELEM
00235       ENDDO ! IT
00236 !
00237       ELSE
00238 !
00239 !     CLASSICAL COMPUTATION
00240 !
00241       DO IELEM = 1 , NELEM
00242 !
00243         X2 = XEL(IELEM+NELMAX)
00244         X3 = XEL(IELEM+2*NELMAX)
00245         Y2 = YEL(IELEM+NELMAX)
00246         Y3 = YEL(IELEM+2*NELMAX)
00247 !
00248         U1 = U(IKLE1(IELEM))
00249         U2 = U(IKLE2(IELEM))
00250         U3 = U(IKLE3(IELEM))
00251         U4 = U(IKLE4(IELEM))
00252         V1 = V(IKLE1(IELEM))
00253         V2 = V(IKLE2(IELEM))
00254         V3 = V(IKLE3(IELEM))
00255         V4 = V(IKLE4(IELEM))
00256 !
00257         F1 = F(IKLE1(IELEM))
00258         F2 = F(IKLE2(IELEM)) - F1
00259         F3 = F(IKLE3(IELEM)) - F1
00260         F4 = F(IKLE4(IELEM)) - F1
00261 !
00262         W1(IELEM)= (((V3+V4+2*V1)*X2+(V3+V4+2*V1)*X3-(Y3+Y2)*U3-(Y3+
00263      &    Y2)*U4-2*(Y3+Y2)*U1)*F3-3*((V3+V4+2*V1)*X3-(V4+V2+2*
00264      &    V1)*X2-(Y3-Y2)*U4-2*(Y3-Y2)*U1-U3*Y3+U2*Y2)*F4-((V4+V2+
00265      &    2*V1)*X2+(V4+V2+2*V1)*X3-(Y3+Y2)*U4-(Y3+Y2)*U2-2*(Y3+Y2
00266      &    )*U1)*F2)*XSUR72
00267 !
00268         W2(IELEM)= (-(((2*V3+3*V4+6*V2+V1)*X3-(V3-V1)*X2-(2*Y3-Y2)
00269      &    *U3-(Y3+Y2)*U1-3*U4*Y3-6*U2*Y3)*F2-(2*(V3+V4+2*V2)*X2
00270      &    -(V3+V4+2*V2)*X3+(Y3-2*Y2)*U3+(Y3-2*Y2)*U4+2*(Y3-2*
00271      &    Y2)*U2)*F3-3*((V3+V4+2*V2)*X3-(V3-V1)*X2-(Y3-Y2)*U3-U4*
00272      &    Y3-2*U2*Y3-U1*Y2)*F4))*XSUR72
00273 !
00274         W3(IELEM)= (((6*V3+3*V4+2*V2+V1)*X2-(V2-V1)*X3-(Y3+Y2)*U1+(
00275      &    Y3-2*Y2)*U2-6*U3*Y2-3*U4*Y2)*F3+((2*V3+V4+V2)*X2-2*(
00276      &    2*V3+V4+V2)*X3+2*(2*Y3-Y2)*U3+(2*Y3-Y2)*U4+(2*Y3-Y2)*
00277      &    U2)*F2-3*((2*V3+V4+V2)*X2-(V2-V1)*X3+(Y3-Y2)*U2-2*U3*
00278      &    Y2-U4*Y2-U1*Y3)*F4)*XSUR72
00279 !
00280         W4(IELEM)= (((3*V3+6*V4+2*V2+V1)*X2-(V2-V1)*X3-(Y3+Y2)*U1+(
00281      &    Y3-2*Y2)*U2-3*U3*Y2-6*U4*Y2)*F3-((2*V3+6*V4+3*V2+V1
00282      &    )*X3-(V3-V1)*X2-(2*Y3-Y2)*U3-(Y3+Y2)*U1-6*U4*Y3-3*U2*
00283      &    Y3)*F2-3*((V3-V1)*X2-(V2-V1)*X3-(Y3-Y2)*U1-U3*Y2+U2*Y3)*
00284      &    F4)*XSUR72
00285 !
00286       ENDDO ! IELEM
00287 !
00288       ENDIF
00289 !
00290 !     FUNCTION F IS QUASI-BUBBLE AND VECTOR U LINEAR
00291 !
00292       ELSEIF(IELMF.EQ.12.AND.IELMU.EQ.11.AND.IELMV.EQ.11) THEN
00293 !
00294       IF(FORMUL(14:16).EQ.'PSI') THEN
00295 !
00296         IF (LNG.EQ.1) WRITE(LU,400)
00297         IF (LNG.EQ.2) WRITE(LU,401)
00298 400   FORMAT(1X,'VC08BB (BIEF) : PSI NON PROGRAMME EN QUASI-BULLE')
00299 401   FORMAT(1X,'VC08BB (BIEF) : PSI NOT IMPLEMENTED FOR QUASI-BUBBLE')
00300         CALL PLANTE(1)
00301         STOP
00302 !
00303       ELSE
00304 !
00305 !     CLASSICAL COMPUTATION
00306 !
00307       DO IELEM = 1 , NELEM
00308 !
00309         X2 = XEL(IELEM+NELMAX)
00310         X3 = XEL(IELEM+2*NELMAX)
00311         Y2 = YEL(IELEM+NELMAX)
00312         Y3 = YEL(IELEM+2*NELMAX)
00313 !
00314         U1 = U(IKLE1(IELEM))
00315         U2 = U(IKLE2(IELEM))
00316         U3 = U(IKLE3(IELEM))
00317         V1 = V(IKLE1(IELEM))
00318         V2 = V(IKLE2(IELEM))
00319         V3 = V(IKLE3(IELEM))
00320 !
00321         F1 = F(IKLE1(IELEM))
00322         F2 = F(IKLE2(IELEM)) - F1
00323         F3 = F(IKLE3(IELEM)) - F1
00324         F4 = F(IKLE4(IELEM)) - F1
00325 !
00326         W1(IELEM)=(((4*V3+V2+7*V1)*X2+(4*V3+V2+7*V1)*X3-4*(Y3+Y2
00327      &   )*U3-(Y3+Y2)*U2-7*(Y3+Y2)*U1)*F3-3*((4*V3+V2+7*V1)*X3
00328      &   -(V3+4*V2+7*V1)*X2-(4*Y3-Y2)*U3-7*(Y3-Y2)*U1-(Y3-4*
00329      &   Y2)*U2)*F4-((V3+4*V2+7*V1)*X2+(V3+4*V2+7*V1)*X3-(Y3+
00330      &   Y2)*U3-4*(Y3+Y2)*U2-7*(Y3+Y2)*U1)*F2)*XSU216
00331         W2(IELEM)=((2*(4*V3+7*V2+V1)*X2-(4*V3+7*V2+V1)*X3+4*(Y3
00332      &   -2*Y2)*U3+7*(Y3-2*Y2)*U2+(Y3-2*Y2)*U1)*F3+3*((4*V3+
00333      &   7*V2+V1)*X3-3*(V3-V1)*X2-(4*Y3-3*Y2)*U3-(Y3+3*Y2)*U1-
00334      &   7*U2*Y3)*F4-3*((3*V3+7*V2+2*V1)*X3-(V3-V1)*X2-(3*Y3-
00335      &   Y2)*U3-(2*Y3+Y2)*U1-7*U2*Y3)*F2)*XSU216
00336         W3(IELEM)=(((7*V3+4*V2+V1)*X2-2*(7*V3+4*V2+V1)*X3+7*(2
00337      &   *Y3-Y2)*U3+4*(2*Y3-Y2)*U2+(2*Y3-Y2)*U1)*F2-3*((7*V3+
00338      &   4*V2+V1)*X2-3*(V2-V1)*X3-(3*Y3+Y2)*U1+(3*Y3-4*Y2)*U2-
00339      &   7*U3*Y2)*F4+3*((7*V3+3*V2+2*V1)*X2-(V2-V1)*X3-(Y3+2*
00340      &   Y2)*U1+(Y3-3*Y2)*U2-7*U3*Y2)*F3)*XSU216
00341         W4(IELEM)=(((5*V3+4*V2+3*V1)*X2-(V2-V1)*X3-(Y3+3*Y2)*U1+(
00342      &   Y3-4*Y2)*U2-5*U3*Y2)*F3-((4*V3+5*V2+3*V1)*X3-(V3-V1)
00343      &   *X2-(4*Y3-Y2)*U3-(3*Y3+Y2)*U1-5*U2*Y3)*F2-3*((V3-V1)*
00344      &   X2-(V2-V1)*X3-(Y3-Y2)*U1-U3*Y2+U2*Y3)*F4)*XSUR72
00345 !
00346       ENDDO ! IELEM
00347 !
00348       ENDIF
00349 !
00350 !-----------------------------------------------------------------------
00351 !
00352       ELSE
00353 !
00354 !-----------------------------------------------------------------------
00355 !
00356         IF (LNG.EQ.1) WRITE(LU,100) IELMF,SF%NAME
00357         IF (LNG.EQ.1) WRITE(LU,200) IELMU,SU%NAME
00358         IF (LNG.EQ.1) WRITE(LU,300)
00359         IF (LNG.EQ.2) WRITE(LU,101) IELMF,SF%NAME
00360         IF (LNG.EQ.1) WRITE(LU,201) IELMU,SU%NAME
00361         IF (LNG.EQ.1) WRITE(LU,301)
00362 100     FORMAT(1X,'VC08BB (BIEF) :',/,
00363      &         1X,'DISCRETISATION DE F : ',1I6,
00364      &         1X,'NOM REEL : ',A6)
00365 200     FORMAT(1X,'DISCRETISATION DE U : ',1I6,
00366      &         1X,'NOM REEL : ',A6)
00367 300     FORMAT(1X,'CAS NON PREVU')
00368 101     FORMAT(1X,'VC08BB (BIEF) :',/,
00369      &         1X,'DISCRETIZATION OF F:',1I6,
00370      &         1X,'REAL NAME: ',A6)
00371 201     FORMAT(1X,'DISCRETIZATION OF U:',1I6,
00372      &         1X,'REAL NAME: ',A6)
00373 301     FORMAT(1X,'CASE NOT IMPLEMENTED')
00374         CALL PLANTE(1)
00375         STOP
00376 !
00377       ENDIF
00378 !
00379 !-----------------------------------------------------------------------
00380 !
00381       RETURN
00382       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0