mt12ac.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\mt12ac.f
00002 !
00073                      SUBROUTINE MT12AC
00074 !                    *****************
00075 !
00076      &(  A11 , A12 , A13 , A14 , A15, A16,
00077      &   A21 , A22 , A23 , A24 , A25, A26,
00078      &   A31 , A32 , A33 , A34 , A35, A36,
00079      &   XMUL,SF,SU,SV,F,U,V,
00080      &   XEL,YEL,SURFAC,
00081      &   IKLE1,IKLE2,IKLE3,IKLE4,IKLE5,IKLE6,
00082      &   NELEM,NELMAX,ICOORD)
00083 !
00084 !***********************************************************************
00085 ! BIEF   V6P1                                   21/08/2010
00086 !***********************************************************************
00087 !
00088 !
00089 !
00090 !
00091 !
00092 !
00093 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00094 !| A11            |<--| ELEMENTS OF MATRIX
00095 !| ...            |<--| ELEMENTS OF MATRIX
00096 !| A36            |<--| ELEMENTS OF MATRIX
00097 !| F              |-->| FUNCTION USED IN THE FORMULA
00098 !| ICOORD         |-->| 1: DERIVATIVE ALONG X, 2: ALONG Y
00099 !| IKLE1          |-->| FIRST POINTS OF TRIANGLES
00100 !| IKLE2          |-->| SECOND POINTS OF TRIANGLES
00101 !| IKLE3          |-->| THIRD POINTS OF TRIANGLES
00102 !| IKLE4          |-->| QUASI-BUBBLE POINT
00103 !| NELEM          |-->| NUMBER OF ELEMENTS
00104 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
00105 !| SF             |-->| STRUCTURE OF FUNCTIONS F
00106 !| SU             |-->| BIEF_OBJ STRUCTURE OF U
00107 !| SURFAC         |-->| AREA OF TRIANGLES
00108 !| SV             |-->| BIEF_OBJ STRUCTURE OF V
00109 !| U              |-->| FUNCTION U USED IN THE FORMULA
00110 !| V              |-->| FUNCTION V USED IN THE FORMULA
00111 !| XEL            |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
00112 !| YEL            |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
00113 !| XMUL           |-->| MULTIPLICATION FACTOR
00114 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00115 !
00116       USE BIEF !, EX_MT12AC => MT12AC
00117 !
00118       IMPLICIT NONE
00119       INTEGER LNG,LU
00120       COMMON/INFO/LNG,LU
00121 !
00122 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00123 !
00124       INTEGER, INTENT(IN) :: NELEM,NELMAX,ICOORD
00125       INTEGER, INTENT(IN) :: IKLE1(NELMAX),IKLE2(NELMAX)
00126       INTEGER, INTENT(IN) :: IKLE3(NELMAX),IKLE4(NELMAX)
00127       INTEGER, INTENT(IN) :: IKLE5(NELMAX),IKLE6(NELMAX)
00128 !
00129       DOUBLE PRECISION, INTENT(INOUT) :: A11(*),A12(*),A13(*)
00130       DOUBLE PRECISION, INTENT(INOUT) :: A14(*),A15(*),A16(*)
00131       DOUBLE PRECISION, INTENT(INOUT) :: A21(*),A22(*),A23(*)
00132       DOUBLE PRECISION, INTENT(INOUT) :: A24(*),A25(*),A26(*)
00133       DOUBLE PRECISION, INTENT(INOUT) :: A31(*),A32(*),A33(*)
00134       DOUBLE PRECISION, INTENT(INOUT) :: A34(*),A35(*),A36(*)
00135       DOUBLE PRECISION, INTENT(IN) :: XMUL
00136       DOUBLE PRECISION, INTENT(IN) :: F(*),U(*),V(*)
00137 !
00138 !     STRUCTURES OF F, U, V
00139       TYPE(BIEF_OBJ), INTENT(IN) :: SF,SU,SV
00140 !
00141       DOUBLE PRECISION, INTENT(IN) :: XEL(NELMAX,3),YEL(NELMAX,3)
00142       DOUBLE PRECISION, INTENT(IN) :: SURFAC(NELMAX)
00143 !
00144 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00145 !
00146       INTEGER IELEM,IELMF,IELMU,IELMV
00147       DOUBLE PRECISION X2,X3,Y2,Y3,F1,F2,F3
00148       DOUBLE PRECISION U1,U2,U3,U4,U5,U6,V1,V2,V3,V4,V5,V6,UX,UY
00149       DOUBLE PRECISION XSUR12 ,XSUR180,XSUR720
00150       DOUBLE PRECISION AUX12  ,AUX180 ,AUX720 , UNSURF
00151 !
00152 !-----------------------------------------------------------------------
00153 !
00154       IELMF=SF%ELM
00155       IELMU=SU%ELM
00156       IELMV=SV%ELM
00157 !
00158       XSUR12  = XMUL /  12.D0
00159       XSUR720 = XMUL / 720.D0
00160       XSUR180 = XMUL / 180.D0
00161 !
00162 !-----------------------------------------------------------------------
00163 !  CASE WHERE F IS OF TYPE P1
00164 !-----------------------------------------------------------------------
00165 !
00166       IF(IELMF.EQ.11) THEN
00167 !
00168       IF(IELMU.EQ.10.AND.IELMV.EQ.10) THEN
00169 !
00170 !================================
00171 !  DERIVATIVE WRT X  =
00172 !================================
00173 !
00174         IF(ICOORD.EQ.1) THEN
00175 !
00176 !   LOOP ON THE ELEMENTS
00177 !
00178         DO IELEM = 1 , NELEM
00179 !
00180 !   INITIALISES THE GEOMETRICAL VARIABLES
00181 !
00182         X2 = XEL(IELEM,2)
00183         X3 = XEL(IELEM,3)
00184         Y2 = YEL(IELEM,2)
00185         Y3 = YEL(IELEM,3)
00186 !
00187         F1  =  F(IKLE1(IELEM))
00188         F2  =  F(IKLE2(IELEM)) - F1
00189         F3  =  F(IKLE3(IELEM)) - F1
00190 !
00191         UX  =  U(IELEM)
00192         UY  =  V(IELEM)
00193 !
00194         UNSURF= 1.D0 / SURFAC(IELEM)
00195         AUX12 = XSUR12 * UNSURF
00196 !
00197 !   EXTRADIAGONAL TERMS
00198 !
00199         A12(IELEM)= 0.D0
00200         A13(IELEM)= 0.D0
00201         A21(IELEM)= 0.D0
00202         A23(IELEM)= 0.D0
00203         A24(IELEM)= (Y3*UX-UY*X3) * (Y3*F2-Y2*F3) * AUX12
00204         A25(IELEM)= A24(IELEM)
00205         A26(IELEM)= A24(IELEM)
00206         A31(IELEM)= 0.D0
00207         A32(IELEM)= 0.D0
00208         A34(IELEM)= (X2*UY-UX*Y2) * (Y3*F2-Y2*F3) * AUX12
00209         A35(IELEM)= A34(IELEM)
00210         A36(IELEM)= A34(IELEM)
00211 !
00212 !   THE SUM OF EACH COLUMN IS 0
00213 !
00214         A14(IELEM)= - A24(IELEM) - A34(IELEM)
00215         A15(IELEM)= A14(IELEM)
00216         A16(IELEM)= A14(IELEM)
00217 !
00218 !C DIAGONAL ELEMENTS
00219 !
00220         A11(IELEM) = 0.D0
00221         A22(IELEM) = 0.D0
00222         A33(IELEM) = 0.D0
00223 !
00224       ENDDO ! IELEM
00225 !
00226         ELSEIF(ICOORD.EQ.2) THEN
00227 !
00228 !================================
00229 !  DERIVATIVE WRT Y  =
00230 !================================
00231 !
00232         DO IELEM = 1 , NELEM
00233 !
00234 !   INITIALISES THE GEOMETRICAL VARIABLES
00235 !
00236         X2  =  XEL(IELEM,2)
00237         X3  =  XEL(IELEM,3)
00238         Y2  =  YEL(IELEM,2)
00239         Y3  =  YEL(IELEM,3)
00240 !
00241         F1  =  F(IKLE1(IELEM))
00242         F2  =  F(IKLE2(IELEM)) - F1
00243         F3  =  F(IKLE3(IELEM)) - F1
00244 !
00245         UX  =  U(IELEM)
00246         UY  =  V(IELEM)
00247 !
00248         UNSURF = 1.D0 / SURFAC(IELEM)
00249         AUX12 = XSUR12 * UNSURF
00250 !
00251 !   EXTRADIAGONAL TERMS
00252 !
00253         A12(IELEM)= 0.D0
00254         A13(IELEM)= 0.D0
00255         A21(IELEM)= 0.D0
00256         A23(IELEM)= 0.D0
00257         A24(IELEM)= (X3*F2-X2*F3) * (UY*X3-UX*Y3) * AUX12
00258         A25(IELEM)=  A24(IELEM)
00259         A26(IELEM)=  A24(IELEM)
00260         A31(IELEM)= 0.D0
00261         A32(IELEM)= 0.D0
00262         A34(IELEM)= (X3*F2-X2*F3) * (UX*Y2-UY*X2) * AUX12
00263         A35(IELEM)=  A34(IELEM)
00264         A36(IELEM)=  A34(IELEM)
00265 !
00266 !   THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
00267 !
00268         A14(IELEM)= - A24(IELEM) - A34(IELEM)
00269         A15(IELEM)=   A14(IELEM)
00270         A16(IELEM)=   A14(IELEM)
00271 !
00272 !C   DIAGONAL TERMS
00273 !
00274         A11(IELEM) = 0.D0
00275         A22(IELEM) = 0.D0
00276         A33(IELEM) = 0.D0
00277 !
00278         ENDDO ! IELEM
00279 !
00280         ELSE
00281 !
00282           IF (LNG.EQ.1) WRITE(LU,200) ICOORD
00283           IF (LNG.EQ.2) WRITE(LU,201) ICOORD
00284           CALL PLANTE(1)
00285           STOP
00286         ENDIF
00287 !
00288       ELSEIF(IELMU.EQ.13) THEN
00289 !
00290 !================================
00291 !  DERIVATIVE WRT X  =
00292 !================================
00293 !
00294         IF(ICOORD.EQ.1) THEN
00295 !
00296 !   LOOP ON THE ELEMENTS
00297 !
00298         DO IELEM = 1 , NELEM
00299 !
00300 !   INITIALISES THE GEOMETRICAL VARIABLES
00301 !
00302         X2 = XEL(IELEM,2)
00303         X3 = XEL(IELEM,3)
00304         Y2 = YEL(IELEM,2)
00305         Y3 = YEL(IELEM,3)
00306 !
00307         F1  =  F(IKLE1(IELEM))
00308         F2  =  F(IKLE2(IELEM)) - F1
00309         F3  =  F(IKLE3(IELEM)) - F1
00310 !
00311         U1  =  U(IKLE1(IELEM))
00312         U2  =  U(IKLE2(IELEM))
00313         U3  =  U(IKLE3(IELEM))
00314         U4  =  U(IKLE4(IELEM))
00315         U5  =  U(IKLE5(IELEM))
00316         U6  =  U(IKLE6(IELEM))
00317         V1  =  V(IKLE1(IELEM))
00318         V2  =  V(IKLE2(IELEM))
00319         V3  =  V(IKLE3(IELEM))
00320         V4  =  V(IKLE4(IELEM))
00321         V5  =  V(IKLE5(IELEM))
00322         V6  =  V(IKLE6(IELEM))
00323 !
00324         UNSURF = 1.D0 / SURFAC(IELEM)
00325         AUX720 = XSUR720 * UNSURF
00326         AUX180 = XSUR180 * UNSURF
00327 !
00328 !   EXTRADIAGONAL TERMS
00329 !
00330         A12(IELEM)= (Y3*F2-Y2*F3) * ((U3+U1-6.D0*U2+4.D0*U6) * (Y3-Y2)
00331      &            + (V3+V1-6.D0*V2+4.D0*V6) * (X2-X3)) * AUX720
00332 !
00333         A13(IELEM)= (Y3*F2-Y2*F3) * ((4.D0*V4-6.D0*V3+V2+V1) * (X2-X3)
00334      &            + (6.D0*U3-4.D0*U4-U2-U1) * (Y2-Y3)) * AUX720
00335 !
00336         A21(IELEM)= (Y3*F2-Y2*F3) * ((4.D0*V5+V3+V2-6.D0*V1) * X3
00337      &            - (U3-6.D0*U1+4.D0*U5+U2) * Y3     ) * AUX720
00338 !
00339         A23(IELEM)= (Y3*F2-Y2*F3) * ((V1-6.D0*V3+4.D0*V4+V2) * X3
00340      &            + (6.D0*U3-4.D0*U4-U2-U1) * Y3     ) * AUX720
00341 !
00342         A24(IELEM)= (Y3*F2-Y2*F3) * ((V3-4.D0*V6-8.D0*V4-4.D0*V5) * X3
00343      &            - (U3-8.D0*U4-4.D0*U6-4.D0*U5) * Y3) * AUX180
00344 !
00345         A25(IELEM)= (Y3*F2-Y2*F3) * ((V1-4.D0*V6-8.D0*V5-4.D0*V4) * X3
00346      &            + (4.D0*U6-U1+4.D0*U4+8.D0*U5) * Y3) * AUX180
00347 !
00348         A26(IELEM)= (Y3*F2-Y2*F3) * ((V2-4.D0*V4-4.D0*V5-8.D0*V6) * X3
00349      &            + (8.D0*U6+4.D0*U4+4.D0*U5-U2) * Y3) * AUX180
00350 !
00351         A31(IELEM)= (Y3*F2-Y2*F3) * ((6.D0*V1-4.D0*V5-V3-V2) * X2
00352      &            + (U3-6.D0*U1+4.D0*U5+U2) * Y2         ) * AUX720
00353 !
00354         A32(IELEM)= (Y3*F2-Y2*F3) * ((-4.D0*V6+6.D0*V2-V1-V3) * X2
00355      &            + (U3+U1+4.D0*U6-6.D0*U2) * Y2        ) * AUX720
00356 !
00357         A34(IELEM)= (Y3*F2-Y2*F3) * ((4.D0*V6+8.D0*V4+4.D0*V5-V3) * X2
00358      &            + (U3-8.D0*U4-4.D0*U6-4.D0*U5) * Y2    ) * AUX180
00359 !
00360         A35(IELEM)= (Y2*F3-Y3*F2) * ((V1-8.D0*V5-4.D0*V4-4.D0*V6) * X2
00361      &            + (4.D0*U6-U1+4.D0*U4+8.D0*U5) * Y2   ) * AUX180
00362 !
00363         A36(IELEM)= (Y2*F3-Y3*F2) * ((V2-4.D0*V4-4.D0*V5-8.D0*V6) * X2
00364      &            + (8.D0*U6+4.D0*U4+4.D0*U5-U2) * Y2   ) * AUX180
00365 !
00366 !   THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
00367 !
00368         A14(IELEM) = - A24(IELEM) - A34(IELEM)
00369         A15(IELEM) = - A25(IELEM) - A35(IELEM)
00370         A16(IELEM) = - A26(IELEM) - A36(IELEM)
00371 !
00372 !   DIAGONAL TERMS
00373 !   THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
00374 !
00375         A11(IELEM) = - A21(IELEM) - A31(IELEM)
00376         A22(IELEM) = - A12(IELEM) - A32(IELEM)
00377         A33(IELEM) = - A13(IELEM) - A23(IELEM)
00378 !
00379       ENDDO ! IELEM
00380 !
00381         ELSEIF(ICOORD.EQ.2) THEN
00382 !
00383 !================================
00384 !  DERIVATIVE WRT Y  =
00385 !================================
00386 !
00387         DO IELEM = 1 , NELEM
00388 !
00389 !   INITIALISES THE GEOMETRICAL VARIABLES
00390 !
00391         X2  =  XEL(IELEM,2)
00392         X3  =  XEL(IELEM,3)
00393         Y2  =  YEL(IELEM,2)
00394         Y3  =  YEL(IELEM,3)
00395 !
00396         F1  =  F(IKLE1(IELEM))
00397         F2  =  F(IKLE2(IELEM)) - F1
00398         F3  =  F(IKLE3(IELEM)) - F1
00399 !
00400         U1  =  U(IKLE1(IELEM))
00401         U2  =  U(IKLE2(IELEM))
00402         U3  =  U(IKLE3(IELEM))
00403         U4  =  U(IKLE4(IELEM))
00404         U5  =  U(IKLE5(IELEM))
00405         U6  =  U(IKLE6(IELEM))
00406         V1  =  V(IKLE1(IELEM))
00407         V2  =  V(IKLE2(IELEM))
00408         V3  =  V(IKLE3(IELEM))
00409         V4  =  V(IKLE4(IELEM))
00410         V5  =  V(IKLE5(IELEM))
00411         V6  =  V(IKLE6(IELEM))
00412 !
00413         UNSURF= 1.D0 / SURFAC(IELEM)
00414         AUX180= XSUR180 * UNSURF
00415         AUX720= XSUR720 * UNSURF
00416 !
00417 !   EXTRADIAGONAL TERMS
00418 !
00419         A12(IELEM)= (X2*F3-X3*F2) * ((4.D0*V6-6.D0*V2+V3+V1) * (X2-X3)
00420      &            + (U1+U3+4.D0*U6-6.D0*U2) * (Y3-Y2)) * AUX720
00421 !
00422         A13(IELEM)= (X2*F3-X3*F2) * ((V2-6.D0*V3+4.D0*V4+V1) * (X2-X3)
00423      &            + (6.D0*U3-U1-U2-4.D0*U4) * (Y2-Y3)) * AUX720
00424 !
00425         A21(IELEM)= (X3*F2-X2*F3) * ((6.D0*V1-V3-V2-4.D0*V5) * X3
00426      &            + (U3-6.D0*U1+U2+4.D0*U5) * Y3      ) * AUX720
00427 !
00428         A23(IELEM)= (X2*F3-X3*F2) * ((V2-6.D0*V3+4.D0*V4+V1) * X3
00429      &            + (6.D0*U3-U1-U2-4.D0*U4) * Y3      ) * AUX720
00430 !
00431         A24(IELEM)= (X2*F3-X3*F2) * ((V3-4.D0*V6-4.D0*V5-8.D0*V4) * X3
00432      &            + (4.D0*U5+8.D0*U4-U3+4.D0*U6) * Y3 ) * AUX180
00433 !
00434         A25(IELEM)= (X2*F3-X3*F2) * ((V1-8.D0*V5-4.D0*V4-4.D0*V6) * X3
00435      &            + (4.D0*U4-U1+8.D0*U5+4.D0*U6) * Y3 ) * AUX180
00436 !
00437         A26(IELEM)= (X2*F3-X3*F2) * ((V2-4.D0*V5-4.D0*V4-8.D0*V6) * X3
00438      &            + (4.D0*U5-U2+8.D0*U6+4.D0*U4) * Y3 ) * AUX180
00439 !
00440         A31(IELEM)= (X2*F3-X3*F2) * ((6.D0*V1-V3-V2-4.D0*V5) * X2
00441      &            + (U3-6.D0*U1+U2+4.D0*U5) * Y2      ) * AUX720
00442 !
00443         A32(IELEM)= (X3*F2-X2*F3) * ((4.D0*V6-6.D0*V2+V3+V1) * X2
00444      &            + (6.D0*U2-U1-U3-4.D0*U6) * Y2      ) * AUX720
00445 !
00446         A34(IELEM)= (X3*F2-X2*F3) * ((V3-4.D0*V6-4.D0*V5-8.D0*V4) * X2
00447      &            + (4.D0*U5+8.D0*U4-U3+4.D0*U6) * Y2 ) * AUX180
00448 !
00449         A35(IELEM)= (X3*F2-X2*F3) * ((V1-8.D0*V5-4.D0*V4-4.D0*V6) * X2
00450      &            + (4.D0*U4-U1+8.D0*U5+4.D0*U6) * Y2 ) * AUX180
00451 !
00452         A36(IELEM)= (X3*F2-X2*F3) * ((V2-4.D0*V5-4.D0*V4-8.D0*V6) * X2
00453      &            + (4.D0*U5-U2+8.D0*U6+4.D0*U4) * Y2 ) * AUX180
00454 !
00455 !   THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
00456 !
00457         A14(IELEM) = - A24(IELEM) - A34(IELEM)
00458         A15(IELEM) = - A25(IELEM) - A35(IELEM)
00459         A16(IELEM) = - A26(IELEM) - A36(IELEM)
00460 !
00461 !   DIAGONAL TERMS
00462 !   THE SUM OF THE MATRIX ROWS IS 0 (VECTOR)
00463 !
00464         A11(IELEM) = - A21(IELEM) - A31(IELEM)
00465         A22(IELEM) = - A12(IELEM) - A32(IELEM)
00466         A33(IELEM) = - A13(IELEM) - A23(IELEM)
00467 !
00468         ENDDO ! IELEM
00469 !
00470         ELSE
00471           IF (LNG.EQ.1) WRITE(LU,200) ICOORD
00472           IF (LNG.EQ.2) WRITE(LU,201) ICOORD
00473           CALL PLANTE(1)
00474           STOP
00475         ENDIF
00476 !
00477       ELSE
00478 !
00479         IF (LNG.EQ.1) WRITE(LU,300) IELMU
00480         IF (LNG.EQ.2) WRITE(LU,301) IELMU
00481         CALL PLANTE(1)
00482         STOP
00483 !
00484       ENDIF
00485 !
00486 !-----------------------------------------------------------------------
00487 !
00488       ELSE
00489         IF (LNG.EQ.1) WRITE(LU,100) IELMF
00490         IF (LNG.EQ.2) WRITE(LU,101) IELMF
00491 100     FORMAT(1X,'MT12AC (BIEF) :',/,
00492      &         1X,'DISCRETISATION DE F : ',1I6,' NON PREVUE')
00493 101     FORMAT(1X,'MT12AC (BIEF) :',/,
00494      &         1X,'DISCRETIZATION OF F : ',1I6,' NOT AVAILABLE')
00495         CALL PLANTE(1)
00496         STOP
00497       ENDIF
00498 !
00499 200   FORMAT(1X,'MT12AC (BIEF) : COMPOSANTE IMPOSSIBLE ',
00500      &          1I6,' VERIFIER ICOORD')
00501 201   FORMAT(1X,'MT12AC (BIEF) : IMPOSSIBLE COMPONENT ',
00502      &          1I6,' CHECK ICOORD')
00503 !
00504 300   FORMAT(1X,'MT12AC (BIEF) :',/,
00505      &       1X,'DISCRETISATION DE U : ',1I6,' NON PREVUE')
00506 301   FORMAT(1X,'MT12AC (BIEF) :',/,
00507      &       1X,'DISCRETIZATION OF U : ',1I6,' NOT AVAILABLE')
00508 !
00509 !-----------------------------------------------------------------------
00510 !
00511       RETURN
00512       END

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