1 : C CFLTAI SOURCE PV 09/03/12 21:16:35 6325 2 : SUBROUTINE CFLTAI(ITAILL,IPMAIL,MELE) 3 : *---------------------------------------------------------------------* 4 : * 5 : * calcul du pas de temps CFL 6 : * CFLTAI calcule le paramètre de taille 7 : * 8 : * entree : 9 : * ipmail maillage actif en entree et en sortie 10 : * mele numero de l'element fini dans nomtp 11 : * sortie 12 : * itaill pointeur sur segment melval : actif et vide en entree 13 : * et actif en sortie 14 : * 15 : *---------------------------------------------------------------------* 16 : IMPLICIT INTEGER(I-N) 17 : IMPLICIT REAL*8(A-H,O-Z) 18 : * 19 : -INC CCOPTIO 20 : -INC SMELEME 21 : -INC SMCOORD 22 : -INC SMCHAML 23 : -INC CCHAMP 24 : * 25 : DIMENSION X(3),X1(2),X2(2) 26 : * 27 : DIMENSION V24(2),V34(2),V41(2) 28 : DIMENSION V35(2),V57(2),V71(2) 29 : DIMENSION V13(3),V12(3),V23(3), DELTA(3) 30 : * 31 : * 32 : MELEME = IPMAIL 33 : MELVAL = ITAILL 34 : * 35 : * 36 : * branchement en fonction de l'élément fini 37 : * 38 : * 1 5 0 5 0 39 : GOTO (99,99,99, 4,99, 6,99, 4,99,10,99,99,99,99,99,99,99,99,99,99, 40 : 2 99,99,99,99,99,99,27,27,46,99,99,99,99,99,99,99,99,99,99,99, 41 : 4 99,46,99,44,99,46,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 42 : 6 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 43 : 8 99,99,99,46,99,99,99,99,99,99,99,99,27,99,99,99,99,99,99,99, 44 : 1 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 45 : 2 99,99,99,99,99,99,99),MELE 46 : * 47 : * cas du triangle à trois noeuds 48 : * 49 : 4 CONTINUE 50 : NEL = NUM(/2) 51 : DO 402 I=1,NEL 52 : V13(1)=XCOOR((IDIM+1)*(NUM(3,I)-1)+1)- 53 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+1) 54 : V13(2)=XCOOR((IDIM+1)*(NUM(3,I)-1)+2)- 55 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+2) 56 : V12(1)=XCOOR((IDIM+1)*(NUM(2,I)-1)+1)- 57 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+1) 58 : V12(2)=XCOOR((IDIM+1)*(NUM(2,I)-1)+2)- 59 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+2) 60 : V23(1)=XCOOR((IDIM+1)*(NUM(3,I)-1)+1)- 61 : & XCOOR((IDIM+1)*(NUM(2,I)-1)+1) 62 : V23(2)=XCOOR((IDIM+1)*(NUM(3,I)-1)+2)- 63 : & XCOOR((IDIM+1)*(NUM(2,I)-1)+2) 64 : V12L=V12(1)*V12(1)+V12(2)*V12(2) 65 : V13L=V13(1)*V13(1)+V13(2)*V13(2) 66 : V23L=V23(1)*V23(1)+V23(2)*V23(2) 67 : XLON=MIN(V12L,V13L,V23L) 68 : VOL = ABS(V12(1)*V13(2) - V12(2)*V13(1)) 69 : XL2 = VOL/MAX(V12L,V13L,V23L) 70 : VELCHE(1,I)=MIN(SQRT(XLON),XL2) 71 : 402 CONTINUE 72 : RETURN 73 : * 74 : * cas du triangle à six noeuds 75 : * 76 : 6 CONTINUE 77 : NEL = NUM(/2) 78 : DO 602 I=1,NEL 79 : * on reprend la formulation du tri3 en considérant les noeuds sommets 80 : * du triangle 1 3 5 81 : V13(1)=XCOOR((IDIM+1)*(NUM(5,I)-1)+1)- 82 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+1) 83 : V13(2)=XCOOR((IDIM+1)*(NUM(5,I)-1)+2)- 84 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+2) 85 : V12(1)=XCOOR((IDIM+1)*(NUM(3,I)-1)+1)- 86 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+1) 87 : V12(2)=XCOOR((IDIM+1)*(NUM(3,I)-1)+2)- 88 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+2) 89 : V23(1)=XCOOR((IDIM+1)*(NUM(5,I)-1)+1)- 90 : & XCOOR((IDIM+1)*(NUM(3,I)-1)+1) 91 : V23(2)=XCOOR((IDIM+1)*(NUM(5,I)-1)+2)- 92 : & XCOOR((IDIM+1)*(NUM(3,I)-1)+2) 93 : V12L=V12(1)*V12(1)+V12(2)*V12(2) 94 : V13L=V13(1)*V13(1)+V13(2)*V13(2) 95 : V23L=V23(1)*V23(1)+V23(2)*V23(2) 96 : XLON=MIN(V12L,V13L,V23L) 97 : VOL = ABS(V12(1)*V13(2) - V12(2)*V13(1)) 98 : XL2 = VOL/MAX(V12L,V13L,V23L) 99 : * pour tenir compte des 3 autres noeuds on divise par 4. 100 : VELCHE(1,I)=MIN(SQRT(XLON),XL2)*0.5D0 101 : 602 CONTINUE 102 : RETURN 103 : * 104 : * cas du qua4 massif 105 : * 106 : 8 CONTINUE 107 : NEL = NUM(/2) 108 : DO 802 I=1,NEL 109 : V13(1)=XCOOR((IDIM+1)*(NUM(3,I)-1)+1)- 110 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+1) 111 : V13(2)=XCOOR((IDIM+1)*(NUM(3,I)-1)+2)- 112 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+2) 113 : V24(1)=XCOOR((IDIM+1)*(NUM(4,I)-1)+1)- 114 : & XCOOR((IDIM+1)*(NUM(2,I)-1)+1) 115 : V24(2)=XCOOR((IDIM+1)*(NUM(4,I)-1)+2)- 116 : & XCOOR((IDIM+1)*(NUM(2,I)-1)+2) 117 : V12(1)=XCOOR((IDIM+1)*(NUM(2,I)-1)+1)- 118 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+1) 119 : V12(2)=XCOOR((IDIM+1)*(NUM(2,I)-1)+2)- 120 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+2) 121 : V23(1)=XCOOR((IDIM+1)*(NUM(3,I)-1)+1)- 122 : & XCOOR((IDIM+1)*(NUM(2,I)-1)+1) 123 : V23(2)=XCOOR((IDIM+1)*(NUM(3,I)-1)+2)- 124 : & XCOOR((IDIM+1)*(NUM(2,I)-1)+2) 125 : V34(1)=XCOOR((IDIM+1)*(NUM(4,I)-1)+1)- 126 : & XCOOR((IDIM+1)*(NUM(3,I)-1)+1) 127 : V34(2)=XCOOR((IDIM+1)*(NUM(4,I)-1)+2)- 128 : & XCOOR((IDIM+1)*(NUM(3,I)-1)+2) 129 : V41(1)=XCOOR((IDIM+1)*(NUM(1,I)-1)+1)- 130 : & XCOOR((IDIM+1)*(NUM(4,I)-1)+1) 131 : V41(2)=XCOOR((IDIM+1)*(NUM(4,I)-1)+2)- 132 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+2) 133 : V12L=V12(1)*V12(1)+V12(2)*V12(2) 134 : V23L=V23(1)*V23(1)+V23(2)*V23(2) 135 : V34L=V34(1)*V34(1)+V34(2)*V34(2) 136 : V41L=V41(1)*V41(1)+V41(2)*V41(2) 137 : V13L=V13(1)*V13(1)+V13(2)*V13(2) 138 : V24L=V24(1)*V24(1)+V24(2)*V24(2) 139 : XLON = MIN(V12L,V23L,V34L,V41L,V13L,V24L) 140 : VELCHE(1,I) = SQRT(XLON) * 0.8D0 141 : 802 CONTINUE 142 : RETURN 143 : * 144 : * cas du qua8 145 : * 146 : 10 CONTINUE 147 : NEL = NUM(/2) 148 : DO 1002 I=1,NEL 149 : * on cherche le plus petit des quatres coins 150 : V13(1) = XCOOR((IDIM+1)*(NUM(3,I)-1)+1)- 151 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+1) 152 : V13(2) = XCOOR((IDIM+1)*(NUM(3,I)-1)+2)- 153 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+2) 154 : V35(1) = XCOOR((IDIM+1)*(NUM(5,I)-1)+1)- 155 : & XCOOR((IDIM+1)*(NUM(3,I)-1)+1) 156 : V35(2) = XCOOR((IDIM+1)*(NUM(5,I)-1)+2)- 157 : & XCOOR((IDIM+1)*(NUM(3,I)-1)+2) 158 : V57(1) = XCOOR((IDIM+1)*(NUM(7,I)-1)+1)- 159 : & XCOOR((IDIM+1)*(NUM(5,I)-1)+1) 160 : V57(2) = XCOOR((IDIM+1)*(NUM(7,I)-1)+2)- 161 : & XCOOR((IDIM+1)*(NUM(5,I)-1)+2) 162 : V71(1) = XCOOR((IDIM+1)*(NUM(1,I)-1)+1)- 163 : & XCOOR((IDIM+1)*(NUM(7,I)-1)+1) 164 : V71(2) = XCOOR((IDIM+1)*(NUM(1,I)-1)+2)- 165 : & XCOOR((IDIM+1)*(NUM(7,I)-1)+2) 166 : DV71 = SQRT(V71(1)*V71(1) + V71(2)*V71(2)) 167 : DV57 = SQRT(V57(1)*V57(1) + V57(2)*V57(2)) 168 : DV35 = SQRT(V35(1)*V35(1) + V35(2)*V35(2)) 169 : DV13 = SQRT(V13(1)*V13(1) + V13(2)*V13(2)) 170 : XLON = MIN(DV71,DV57,DV35,DV13) 171 : * de plus on regarde les médianes pour considérer le 'cisaillement' 172 : * eventuel de l'élément 173 : XM1=(XCOOR((IDIM+1)*(NUM(3,I)-1)+1)+ 174 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+1))*0.5D0 175 : YM1=(XCOOR((IDIM+1)*(NUM(3,I)-1)+2)+ 176 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+2))*0.5D0 177 : XM2=(XCOOR((IDIM+1)*(NUM(5,I)-1)+1)+ 178 : & XCOOR((IDIM+1)*(NUM(3,I)-1)+1))*0.5D0 179 : YM2=(XCOOR((IDIM+1)*(NUM(5,I)-1)+2)+ 180 : & XCOOR((IDIM+1)*(NUM(3,I)-1)+2))*0.5D0 181 : XM3=(XCOOR((IDIM+1)*(NUM(7,I)-1)+1)+ 182 : & XCOOR((IDIM+1)*(NUM(5,I)-1)+1))*0.5D0 183 : YM3=(XCOOR((IDIM+1)*(NUM(7,I)-1)+2)+ 184 : & XCOOR((IDIM+1)*(NUM(5,I)-1)+2))*0.5D0 185 : XM4=(XCOOR((IDIM+1)*(NUM(1,I)-1)+1)+ 186 : & XCOOR((IDIM+1)*(NUM(7,I)-1)+1))*0.5D0 187 : YM4=(XCOOR((IDIM+1)*(NUM(1,I)-1)+2)+ 188 : & XCOOR((IDIM+1)*(NUM(7,I)-1)+2))*0.5D0 189 : AX=XM2-XM4 190 : AY=YM2-YM4 191 : BX=XM3-XM1 192 : BY=YM3-YM1 193 : AL=DSQRT(AX*AX+AY*AY) 194 : BL=DSQRT(BX*BX+BY*BY) 195 : ANX=-AY/AL 196 : ANY= AX/AL 197 : BNX=-BY/BL 198 : BNY= BX/BL 199 : H1=DABS(BX*ANX+BY*ANY) 200 : H2=DABS(AX*BNX+AY*BNY) 201 : XLON=DMIN1(XLON,H1) 202 : XLON=DMIN1(XLON,H2) 203 : * on divise par 2 pour prendre en compte le fait qu'il y a deux noeuds par coté 204 : * on prend un facteur de securité de 0.5 205 : VELCHE(1,I) = XLON * 0.25D0 206 : 1002 CONTINUE 207 : RETURN 208 : * 209 : * cas de la barre basé sur un seg2 210 : * 211 : 46 CONTINUE 212 : NEL = NUM(/2) 213 : DO 4602 I=1,NEL 214 : DO 4601 J=1,IDIM 215 : X(J)=XCOOR((IDIM+1)*(NUM(2,I)-1)+J)- 216 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+J) 217 : VELCHE(1,I) = VELCHE(1,I) + X(J)*X(J) 218 : 4601 CONTINUE 219 : VELCHE(1,I) = SQRT(VELCHE(1,I)) 220 : 4602 CONTINUE 221 : RETURN 222 : * 223 : * cas du coq2 basé sur un seg2 (en 2d: axisymétrique) 224 : * la taille diffère du cas barre lorsqu'on est proche de l'axe 225 : 44 CONTINUE 226 : NEL = NUM(/2) 227 : DO 4402 I=1,NEL 228 : XTEST1 = 0.D0 229 : XTEST2 = 0.D0 230 : DO 4401 J=1,IDIM 231 : X2(J)=XCOOR((IDIM+1)*(NUM(2,I)-1)+J) 232 : X1(J)=XCOOR((IDIM+1)*(NUM(1,I)-1)+J) 233 : X(J) = X2(J) - X1(J) 234 : VELCHE(1,I) = VELCHE(1,I) + X(J)*X(J) 235 : XTEST1 = XTEST1 + X1(J)*X1(J) 236 : XTEST2 = XTEST2 + X2(J)*X2(J) 237 : 4401 CONTINUE 238 : XTEST1 = MIN(XTEST1,XTEST2) 239 : IF ((XTEST1.LT.VELCHE(1,I)).AND.(IDIM.EQ.2)) THEN 240 : VELCHE(1,I) = VELCHE(1,I) / 4 241 : ENDIF 242 : VELCHE(1,I) = SQRT(VELCHE(1,I)) 243 : 4402 CONTINUE 244 : RETURN 245 : * 246 : * cas du coq3 basé sur un seg3 247 : * 248 : 27 CONTINUE 249 : NEL = NUM(/2) 250 : DO 2705 I=1,NEL 251 : AIRE = 0.D0 252 : V12L = 0.D0 253 : V13L = 0.D0 254 : V23L = 0.D0 255 : DO 2701 J=1,IDIM 256 : V13(J)=XCOOR((IDIM+1)*(NUM(3,I)-1)+J)- 257 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+J) 258 : V12(J)=XCOOR((IDIM+1)*(NUM(2,I)-1)+J)- 259 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+J) 260 : V23(J)=XCOOR((IDIM+1)*(NUM(3,I)-1)+J)- 261 : & XCOOR((IDIM+1)*(NUM(2,I)-1)+J) 262 : V12L = V12L + V12(J)*V12(J) 263 : V13L = V13L + V13(J)*V13(J) 264 : V23L = V23L + V23(J)*V23(J) 265 : 2701 CONTINUE 266 : DO 2702 J=1,IDIM 267 : J0 = MOD(J,IDIM) 268 : J1 = MOD((J+1),IDIM) 269 : DELTA(J) = V12(J0+1)*V13(J1+1) - V12(J1+1)*V13(J0+1) 270 : IF (IDIM.EQ.3) THEN 271 : AIRE = AIRE + DELTA(J)*DELTA(J) 272 : ENDIF 273 : 2702 CONTINUE 274 : IF (IDIM.EQ.3) THEN 275 : AIRE = 0.5D0 * SQRT(AIRE) 276 : ELSE 277 : AIRE = 0.5D0 * DELTA(1) 278 : ENDIF 279 : VELCHE(1,I) = 2.D0 * AIRE / SQRT(MAX(V12L,V13L,V23L)) 280 : 2705 CONTINUE 281 : RETURN 282 : * 283 : 99 CONTINUE 284 : MOTERR(1:4)=NOMTP(MELE) 285 : MOTERR(9:12)='CFLT' 285 : CALL ERREUR(86) 286 : * 287 : RETURN 288 : END 289 : 290 : 291 : 292 : 293 : 294 : 295 : 296 : 297 :
© Cast3M 2003 - All rights reserved.
Disclaimer