Sources Esope | Cast3M

cfltai
C CFLTAI    SOURCE    CB215821  17/11/30    21:15:27     9639           
      SUBROUTINE CFLTAI(ITAILL,IPMAIL,MELE)
*---------------------------------------------------------------------*
*
*          calcul du pas de temps CFL
*          CFLTAI calcule le paramètre de taille
*
* entree :
*         ipmail maillage actif en entree et en sortie
*         mele   numero de l'element fini dans nomtp
* sortie
*         itaill pointeur sur segment melval : actif et vide en entree
*                                                   et actif en sortie
*
*---------------------------------------------------------------------*
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
*
 
-INC PPARAM
-INC CCOPTIO
-INC SMELEME
-INC SMCOORD
-INC SMCHAML
-INC CCHAMP
*
      DIMENSION X(3),X1(2),X2(2)
*
      DIMENSION V24(2),V34(2),V41(2)
      DIMENSION V35(2),V57(2),V71(2)
      DIMENSION V13(3),V12(3),V23(3), DELTA(3)
*
*
      MELEME = IPMAIL
      MELVAL = ITAILL
*
*
*    branchement en fonction de l'élément fini
*
*           1            5              0              5             0
      GOTO (99,99,99, 4,99, 6,99, 4,99,10,99,99,99,99,99,99,99,99,99,99,
     2      99,99,99,99,99,99,27,27,46,99,99,99,99,99,99,99,99,99,99,99,
     4      99,46,99,44,99,46,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
     6      99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
     8      99,99,99,46,99,99,99,99,99,99,99,99,27,99,99,99,99,99,99,99,
     1      99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
     2      99,99,99,99,99,99,99),MELE
*
* cas du triangle à trois noeuds
*
 4    CONTINUE
        NEL = NUM(/2)
        DO 402 I=1,NEL
          V13(1)=XCOOR((IDIM+1)*(NUM(3,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+1)
          V13(2)=XCOOR((IDIM+1)*(NUM(3,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+2)
          V12(1)=XCOOR((IDIM+1)*(NUM(2,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+1)
          V12(2)=XCOOR((IDIM+1)*(NUM(2,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+2)
          V23(1)=XCOOR((IDIM+1)*(NUM(3,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(2,I)-1)+1)
          V23(2)=XCOOR((IDIM+1)*(NUM(3,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(2,I)-1)+2)
         V12L=V12(1)*V12(1)+V12(2)*V12(2)
         V13L=V13(1)*V13(1)+V13(2)*V13(2)
         V23L=V23(1)*V23(1)+V23(2)*V23(2)
         XLON=MIN(V12L,V13L,V23L)
         VOL = ABS(V12(1)*V13(2) - V12(2)*V13(1))
         XL2 = VOL/MAX(V12L,V13L,V23L)
         VELCHE(1,I)=MIN(SQRT(XLON),XL2)
 402   CONTINUE
       RETURN
*
* cas du triangle à six noeuds
*
 6    CONTINUE
        NEL = NUM(/2)
        DO 602 I=1,NEL
* on reprend la formulation du tri3 en considérant les noeuds sommets
* du triangle 1 3 5
          V13(1)=XCOOR((IDIM+1)*(NUM(5,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+1)
          V13(2)=XCOOR((IDIM+1)*(NUM(5,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+2)
          V12(1)=XCOOR((IDIM+1)*(NUM(3,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+1)
          V12(2)=XCOOR((IDIM+1)*(NUM(3,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+2)
          V23(1)=XCOOR((IDIM+1)*(NUM(5,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(3,I)-1)+1)
          V23(2)=XCOOR((IDIM+1)*(NUM(5,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(3,I)-1)+2)
         V12L=V12(1)*V12(1)+V12(2)*V12(2)
         V13L=V13(1)*V13(1)+V13(2)*V13(2)
         V23L=V23(1)*V23(1)+V23(2)*V23(2)
         XLON=MIN(V12L,V13L,V23L)
         VOL = ABS(V12(1)*V13(2) - V12(2)*V13(1))
         XL2 = VOL/MAX(V12L,V13L,V23L)
* pour tenir compte des 3 autres noeuds on divise par 4.
         VELCHE(1,I)=MIN(SQRT(XLON),XL2)*0.5D0
 602   CONTINUE
       RETURN
*
* cas du qua4 massif
*
 8    CONTINUE
        NEL = NUM(/2)
        DO 802 I=1,NEL
          V13(1)=XCOOR((IDIM+1)*(NUM(3,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+1)
          V13(2)=XCOOR((IDIM+1)*(NUM(3,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+2)
          V24(1)=XCOOR((IDIM+1)*(NUM(4,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(2,I)-1)+1)
          V24(2)=XCOOR((IDIM+1)*(NUM(4,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(2,I)-1)+2)
          V12(1)=XCOOR((IDIM+1)*(NUM(2,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+1)
          V12(2)=XCOOR((IDIM+1)*(NUM(2,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+2)
          V23(1)=XCOOR((IDIM+1)*(NUM(3,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(2,I)-1)+1)
          V23(2)=XCOOR((IDIM+1)*(NUM(3,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(2,I)-1)+2)
          V34(1)=XCOOR((IDIM+1)*(NUM(4,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(3,I)-1)+1)
          V34(2)=XCOOR((IDIM+1)*(NUM(4,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(3,I)-1)+2)
          V41(1)=XCOOR((IDIM+1)*(NUM(1,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(4,I)-1)+1)
          V41(2)=XCOOR((IDIM+1)*(NUM(4,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+2)
          V12L=V12(1)*V12(1)+V12(2)*V12(2)
          V23L=V23(1)*V23(1)+V23(2)*V23(2)
          V34L=V34(1)*V34(1)+V34(2)*V34(2)
          V41L=V41(1)*V41(1)+V41(2)*V41(2)
          V13L=V13(1)*V13(1)+V13(2)*V13(2)
          V24L=V24(1)*V24(1)+V24(2)*V24(2)
          XLON = MIN(V12L,V23L,V34L,V41L,V13L,V24L)
          VELCHE(1,I) = SQRT(XLON) * 0.8D0
 802    CONTINUE
        RETURN
*
* cas du qua8
*
 10   CONTINUE
        NEL = NUM(/2)
        DO 1002 I=1,NEL
* on cherche le plus petit des quatres coins
          V13(1) = XCOOR((IDIM+1)*(NUM(3,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+1)
          V13(2) = XCOOR((IDIM+1)*(NUM(3,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+2)
          V35(1) = XCOOR((IDIM+1)*(NUM(5,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(3,I)-1)+1)
          V35(2) = XCOOR((IDIM+1)*(NUM(5,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(3,I)-1)+2)
          V57(1) = XCOOR((IDIM+1)*(NUM(7,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(5,I)-1)+1)
          V57(2) = XCOOR((IDIM+1)*(NUM(7,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(5,I)-1)+2)
          V71(1) = XCOOR((IDIM+1)*(NUM(1,I)-1)+1)-
     &                  XCOOR((IDIM+1)*(NUM(7,I)-1)+1)
          V71(2) = XCOOR((IDIM+1)*(NUM(1,I)-1)+2)-
     &                  XCOOR((IDIM+1)*(NUM(7,I)-1)+2)
          DV71 = SQRT(V71(1)*V71(1) + V71(2)*V71(2))
          DV57 = SQRT(V57(1)*V57(1) + V57(2)*V57(2))
          DV35 = SQRT(V35(1)*V35(1) + V35(2)*V35(2))
          DV13 = SQRT(V13(1)*V13(1) + V13(2)*V13(2))
          XLON = MIN(DV71,DV57,DV35,DV13)
* de plus on regarde les médianes pour considérer le 'cisaillement'
* eventuel de l'élément
          XM1=(XCOOR((IDIM+1)*(NUM(3,I)-1)+1)+
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+1))*0.5D0
          YM1=(XCOOR((IDIM+1)*(NUM(3,I)-1)+2)+
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+2))*0.5D0
          XM2=(XCOOR((IDIM+1)*(NUM(5,I)-1)+1)+
     &                  XCOOR((IDIM+1)*(NUM(3,I)-1)+1))*0.5D0
          YM2=(XCOOR((IDIM+1)*(NUM(5,I)-1)+2)+
     &                  XCOOR((IDIM+1)*(NUM(3,I)-1)+2))*0.5D0
          XM3=(XCOOR((IDIM+1)*(NUM(7,I)-1)+1)+
     &                  XCOOR((IDIM+1)*(NUM(5,I)-1)+1))*0.5D0
          YM3=(XCOOR((IDIM+1)*(NUM(7,I)-1)+2)+
     &                  XCOOR((IDIM+1)*(NUM(5,I)-1)+2))*0.5D0
          XM4=(XCOOR((IDIM+1)*(NUM(1,I)-1)+1)+
     &                  XCOOR((IDIM+1)*(NUM(7,I)-1)+1))*0.5D0
          YM4=(XCOOR((IDIM+1)*(NUM(1,I)-1)+2)+
     &                  XCOOR((IDIM+1)*(NUM(7,I)-1)+2))*0.5D0
          AX=XM2-XM4
          AY=YM2-YM4
          BX=XM3-XM1
          BY=YM3-YM1
          AL=SQRT(AX*AX+AY*AY)
          BL=SQRT(BX*BX+BY*BY)
          ANX=-AY/AL
          ANY= AX/AL
          BNX=-BY/BL
          BNY= BX/BL
          H1=ABS(BX*ANX+BY*ANY)
          H2=ABS(AX*BNX+AY*BNY)
          XLON=MIN(XLON,H1)
          XLON=MIN(XLON,H2)
* on divise par 2 pour prendre en compte le fait qu'il y a deux noeuds par coté
* on prend un facteur de securité de 0.5
          VELCHE(1,I) = XLON * 0.25D0
 1002   CONTINUE
        RETURN
*
* cas de la barre basé sur un seg2
*
 46   CONTINUE
        NEL = NUM(/2)
        DO 4602 I=1,NEL
            DO 4601 J=1,IDIM
              X(J)=XCOOR((IDIM+1)*(NUM(2,I)-1)+J)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+J)
              VELCHE(1,I) = VELCHE(1,I) + X(J)*X(J)
 4601      CONTINUE
           VELCHE(1,I) =  SQRT(VELCHE(1,I))
 4602   CONTINUE
        RETURN
*
* cas du coq2 basé sur un seg2 (en 2d: axisymétrique)
* la taille diffère du cas barre lorsqu'on est proche de l'axe
 44   CONTINUE
        NEL = NUM(/2)
        DO 4402 I=1,NEL
            XTEST1 = 0.D0
            XTEST2 = 0.D0
            DO 4401 J=1,IDIM
               X2(J)=XCOOR((IDIM+1)*(NUM(2,I)-1)+J)
               X1(J)=XCOOR((IDIM+1)*(NUM(1,I)-1)+J)
               X(J) = X2(J) - X1(J)
              VELCHE(1,I) = VELCHE(1,I) + X(J)*X(J)
              XTEST1 = XTEST1 + X1(J)*X1(J)
              XTEST2 = XTEST2 + X2(J)*X2(J)
 4401      CONTINUE
           XTEST1 = MIN(XTEST1,XTEST2)
           IF ((XTEST1.LT.VELCHE(1,I)).AND.(IDIM.EQ.2)) THEN
             VELCHE(1,I) = VELCHE(1,I) / 4
           ENDIF
           VELCHE(1,I) =  SQRT(VELCHE(1,I))
 4402   CONTINUE
        RETURN
*
* cas du coq3 basé sur un seg3
*
 27   CONTINUE
        NEL = NUM(/2)
        DO 2705 I=1,NEL
            AIRE = 0.D0
            V12L = 0.D0
            V13L = 0.D0
            V23L = 0.D0
            DO 2701 J=1,IDIM
               V13(J)=XCOOR((IDIM+1)*(NUM(3,I)-1)+J)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+J)
               V12(J)=XCOOR((IDIM+1)*(NUM(2,I)-1)+J)-
     &                  XCOOR((IDIM+1)*(NUM(1,I)-1)+J)
               V23(J)=XCOOR((IDIM+1)*(NUM(3,I)-1)+J)-
     &                  XCOOR((IDIM+1)*(NUM(2,I)-1)+J)
               V12L = V12L + V12(J)*V12(J)
               V13L = V13L + V13(J)*V13(J)
               V23L = V23L + V23(J)*V23(J)
 2701      CONTINUE
           DO 2702 J=1,IDIM
              J0 = MOD(J,IDIM)
              J1 = MOD((J+1),IDIM)
              DELTA(J) = V12(J0+1)*V13(J1+1) - V12(J1+1)*V13(J0+1)
              IF (IDIM.EQ.3) THEN
                AIRE = AIRE + DELTA(J)*DELTA(J)
              ENDIF
 2702      CONTINUE
           IF (IDIM.EQ.3) THEN
             AIRE = 0.5D0 * SQRT(AIRE)
           ELSE
             AIRE = 0.5D0 * DELTA(1)
           ENDIF
           VELCHE(1,I) = 2.D0 * AIRE / SQRT(MAX(V12L,V13L,V23L))
 2705   CONTINUE
        RETURN
*
   99 CONTINUE
      MOTERR(1:4)=NOMTP(MELE)
      MOTERR(9:12)='CFLT'
      CALL ERREUR(86)
*
      RETURN
      END