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