1 : C CHABOK SOURCE CHAT 05/01/12 21:53:52 5004 2 : SUBROUTINE CHABOK(YUNG,XNU,IA,EI, 3 : 1 XMAT,ALPHA1,JA,IBOU,SI,DEPS,EPST,EPSTAR,AMTRI,KERRE, 4 : 2 SN,ALPHA2,NUMCHA,ecou,necou) 5 : C*********************************************************************** 6 : C INTEGRATION MODELE DE CHABOCHE 7 : C*********************************************************************** 8 : IMPLICIT INTEGER(I-N) 9 : IMPLICIT REAL *8(A-H,O-Z) 10 : DIMENSION XMAT(*),EI(*),ALPHA1(*),ALPHA2(*) 11 : DIMENSION AMTRI(*),YK5(18),DYK1(18),DYK2(18),DYK3(18), 12 : . DYK4(18),EINT(6),W1INT(6),W2INT(6) 13 : C 14 : * Commun ECOU: sert de fourre-tout pour les tableaux 15 : * 16 : SEGMENT ECOU 17 : *** COMMON/ECOU/TEST,ALFAH, 18 : REAL*8 TEST, ALFAH, 19 : 1 HPAS,TEMPS,COVNMS(6),VECPRO(9),VALPRO(6), 20 : 2 CVNMSD(12),STOT(6),SIGEL(6),DSIGP(6),SIGT(6),W1(6),W2(6), 21 : 1 DALPHA(6),DSIGO(6),E(12),XINV(3), 22 : 2 SIPLAD(6),DSIGP0(6),TET,TETI 23 : ENDSEGMENT 24 : 25 : C COMMON/ECOU/TEST,ALFAH, 26 : C 1 HPAS,TEMPS,COVNMS(6),VECPRO(9),VALPRO(6), 27 : C 2 CVNMSD(12),STOT(6),SIGEL(6),DSIGP(6),SIGT(6),W1(6),W2(6), 28 : C 1 DALPHA(6),DSIGO(6),E(12),XINV(3), 29 : C 2 SIPLAD(6),DSIGP0(6),TET,TETI 30 : C 31 : * 32 : * Commun NECOU utilisé dans ECOINC 33 : * 34 : SEGMENT NECOU 35 : * COMMON/NECOU/NCOURB,IPLAST,IT,IMAPLA,ISOTRO, 36 : INTEGER NCOURB,IPLAST,IT,IMAPLA,ISOTRO, 37 : . ITYP,IFOUR,IFLUAG, 38 : . ICINE,ITHER,IFLUPL,ICYCL,IBI, 39 : . JFLUAG,KFLUAG,LFLUAG, 40 : . IRELAX,JNTRIN,MFLUAG,JSOUFL,JGRDEF 41 : ENDSEGMENT 42 : C COMMON/NECOU/NCOURB,IPLAST,IT,IMAPLA,ISOTRO, 43 : C . ITYP,IFOUR,IFLUAG, 44 : C . ICINE,ITHER,IFLUPL,ICYCL,IBI, 45 : C . JFLUAG,KFLUAG,LFLUAG, 46 : C . IRELAX,JNTRIN,MFLUAG,JSOUFL,JGRDEF 47 : C 48 : C EN SORTIE 49 : C KERRE = 2 SI ON ATTEINT LE NOMBRE MAX D'ITERATIONS INTERNES 50 : C 51 : DATA ITMAX/15/ 52 : DATA IDECOU/20/ 53 : DATA RAC2/0.707106781D00/ 54 : PREDEF=1.D-6 55 : PRESIG=YUNG*PREDEF 56 : PREPH=1.D-3 57 : PRERUN=TEST 58 : PREMIN=3.D-1*PRERUN 59 : C 60 : C INITIALISATION A 0. DES VECTEURS POUR TREANOR 61 : C 62 : NVY=18 63 : DO 3341 IB=1,NVY 64 : YK5(IB)=0. 65 : DYK1(IB)=0. 66 : DYK2(IB)=0. 67 : DYK3(IB)=0. 68 : DYK4(IB)=0. 69 : 3341 CONTINUE 70 : A2=0. 71 : C2=0. 72 : B=0. 73 : RM=0. 74 : PHI=1.D00 75 : ICOD=0 75 : CALL CHALIM(EPSTAR,R,XMAT,TET,ICOD,A1,C1,A2,C2,R0,RM,B, 76 : . PHI,PSI,OME,ICENT2,IDIAM,NUMCHA) 77 : ELT=YUNG/(1.+XNU) 78 : IF(ITHER.NE.0) ELT=ELT*EI(IA)/YUNG 79 : G=ELT*0.5 80 : SAC1=A1*C1 81 : AC1=SAC1*0.66666667 82 : SAC2=A2*C2 83 : AC2=SAC2*0.66666667 84 : SAC12=SAC1+SAC2 85 : AC12=AC1+AC2 86 : GO TO (101,102,102,103,105,102,103,103,109,999, 87 : . 999,999,102,101),ITYP 88 : 101 E(7)=ELT*(1.-XNU)/(1.-2.*XNU) 89 : E(8)=ELT*XNU/(1.-2.*XNU) 90 : E(9)=AC12 91 : E(10)=AC1 92 : E(11)=AC2 93 : GO TO 200 94 : 102 E(7)=ELT/(1.-XNU) 95 : E(8)=E(7)*XNU 96 : E(9)=SAC12 97 : E(10)=SAC1 98 : E(11)=SAC2 99 : GO TO 200 100 : 103 E(7)=ELT*(1.+XNU) 101 : E(8)=0. 102 : E(9)=SAC12 103 : E(10)=SAC1 104 : E(11)=SAC2 105 : GO TO 200 106 : 105 CONTINUE 107 : 109 CONTINUE 108 : C 109 : C INITIALISATIONS 110 : C 111 : 200 ITER=0 111 : CALL CHAREM(IDIAM,R0,RM,B,R,EPSTAR,EPST,DEPS,STOT, 112 : .W1,W2,E,IBOU,ICENT2,JA,ALPHA1,ALPHA2,PHI,PSI,OME) 113 : SI=VONMIS(E,ITYP,ALFAH,COVNMS) 114 : C 115 : C PREMIER DEPS 116 : C 117 : ICAS=1 117 : CALL CHAFON(EPST,AMTRI,AMTRI(7),AMTRI(13),DYK1, 118 : . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT, 119 : . DEPS,R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou) 120 : 670 DEPS0=DEPS 121 : DEPSI=DEPS 122 : IF(ITER.EQ.0) GO TO 666 122 : CALL CHAREM(IDIAM,R0,RM,B,R,EPSTAR,EPST,DBID,STOT, 123 : . W1,W2,E,IBOU,ICENT2,JA,ALPHA1,ALPHA2,PHI,PSI,OME) 124 : SI=VONMIS(E,ITYP,ALFAH,COVNMS) 125 : 666 CONTINUE 126 : C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 127 : C ITERATIONS INTERNES 128 : C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 129 : 555 ITER=ITER+1 130 : C 131 : C CALCUL PAR LA METHODE DE TREANOR 132 : C 133 : DO 413 IB=1,IBOU 134 : EINT(IB)=E(IB) 135 : W1INT(IB)=W1(IB) 136 : IF(ICENT2.EQ.0) GO TO 413 137 : W2INT(IB)=W2(IB) 138 : 413 CONTINUE 139 : SINT=SI 140 : DEDEP=0. 141 : HDEP=DEPSI 142 : IF(ITER.EQ.1) HDEP0=HDEP 143 : EPSINT=EPST 144 : C----------------------------------------------------------------------- 145 : C BOUCLE SUR LES SOUS-PAS 146 : C----------------------------------------------------------------------- 147 : DO 416 IDECO=1,IDECOU 148 : JESSAY=0 149 : DEPRES=DEPSI-DEDEP 150 : IF(ABS(HDEP).GT.ABS(DEPRES)) HDEP=DEPRES 151 : IF(IDECO.EQ.1) GO TO 417 152 : C 153 : 424 CONTINUE 154 : JESSAY=JESSAY+1 155 : DO 418 IB=1,IBOU 156 : E(IB)=EINT(IB) 157 : W1(IB)=W1INT(IB) 158 : IF(ICENT2.EQ.0) GO TO 418 159 : W2(IB)=W2INT(IB) 160 : 418 CONTINUE 161 : SI=SINT 162 : 417 ICAS=2 163 : WEP=0.5D00 163 : CALL CHAFON( EPSINT ,AMTRI,AMTRI(7),AMTRI(13),DYK1, 164 : . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,HDEP, 165 : . R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou) 166 : C 167 : EPST=EPSINT+0.5D00*HDEP 167 : CALL CHAFON( EPST ,AMTRI,AMTRI(7),AMTRI(13),DYK2, 168 : . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,HDEP, 169 : . R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou) 170 : C 171 : WEP=1.D00 171 : CALL CHAFON( EPST ,AMTRI,AMTRI(7),AMTRI(13),DYK3, 172 : . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,HDEP, 173 : . R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou) 174 : C 175 : PHMAX=0. 176 : DO 4101 IB=1,NVY 177 : IF(ABS(DYK2(IB)-DYK1(IB)).LT.PRESIG) GO TO 4101 178 : PH=2.*(DYK2(IB)-DYK3(IB))/(DYK2(IB)-DYK1(IB)) 179 : IF(PHMAX.LT.PH) PHMAX=PH 180 : 4101 CONTINUE 181 : IF(PHMAX.LE.0.) GO TO 4102 182 : IF(PHMAX.GE.PREPH) GO TO 4123 183 : XL1=1.-PHMAX*0.5 184 : XL2=0.5-PHMAX*0.166666666666666667 185 : XL3=0.166666666666666667-PHMAX*0.04166666666666666667 186 : XL4=-PHMAX*0.166666666666666667 187 : GO TO 4103 188 : 4123 CONTINUE 189 : XL1=(1.-EXP(-PHMAX))/PHMAX 190 : XL2=(1.-XL1)/PHMAX 191 : XL3=(0.5-XL2)/PHMAX 192 : XL4=XL1-2.*XL2 193 : GO TO 4103 194 : 4102 PHMAX=0. 195 : XL1=1. 196 : XL2=0.5 197 : XL3=0.166666666666666667 198 : XL4=0. 199 : 4103 CONTINUE 200 : XM1=-XL2+4.*XL3 201 : XM2=2.*(XL2-2.*XL3) 202 : XM3=4.*XL3-3.*XL2 203 : C 204 : C NOUVELLE ESTIMATION DE LA SOLUTION 205 : C 206 : DO 4104 IB=1,IBOU 207 : E(IB)=EINT(IB)+2.*XL2*DYK3(IB)+XL2*PHMAX*DYK2(IB)+XL4*DYK1(IB) 208 : YK5(IB)=E(IB) 209 : W1(IB)=W1INT(IB)+2.*XL2*DYK3(6+IB) 210 : . +XL2*PHMAX*DYK2(6+IB)+XL4*DYK1(6+IB) 211 : YK5(6+IB)=W1(IB) 212 : IF(ICENT2.EQ.0) GO TO 4104 213 : W2(IB)=W2INT(IB)+2.*XL2*DYK3(12+IB) 214 : . +XL2*PHMAX*DYK2(12+IB)+XL4*DYK1(12+IB) 215 : YK5(12+IB)=W2(IB) 216 : 4104 CONTINUE 217 : SI=VONMIS(E,ITYP,ALFAH,COVNMS) 218 : C 219 : EPST=EPSINT+HDEP 220 : WEP=1.D00 220 : CALL CHAFON( EPST ,AMTRI,AMTRI(7),AMTRI(13),DYK4, 221 : . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,HDEP, 222 : . R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou) 223 : C 224 : C ESTIMATION FINALE DE LA SOLUTION 225 : C 226 : FAC1=1.+PHMAX*XM3+2.*PHMAX*XM2 227 : FAC2=XL1+XM3+0.5*XM2*PHMAX 228 : FAC3=XM2*(1.+0.5*PHMAX) 229 : DO 4105 IB=1,IBOU 230 : E(IB)=EINT(IB)*FAC1+DYK1(IB)*FAC2+DYK2(IB)*FAC3 231 : . +DYK3(IB)*XM2+DYK4(IB)*XM1+YK5(IB)*XM1*PHMAX 232 : SIGT(IB)=E(IB) 233 : W1(IB)=W1INT(IB)*FAC1+DYK1(6+IB)*FAC2 234 : . +DYK2(6+IB)*FAC3 +DYK3(6+IB)*XM2 +DYK4(6+IB)*XM1 235 : . +YK5(6+IB)*XM1*PHMAX 236 : IF(ICENT2.EQ.0) GO TO 4105 237 : W2(IB)= W2INT(IB)*FAC1 +DYK1(12+IB)*FAC2 238 : . +DYK2(12+IB)*FAC3 +DYK3(12+IB)*XM2 +DYK4(12+IB)*XM1 239 : . +YK5(12+IB)*XM1*PHMAX 240 : 4105 CONTINUE 241 : SI=VONMIS(SIGT,ITYP,ALFAH,COVNMS) 242 : C 243 : C CALCUL DE L'ERREUR 244 : C 245 : ERMAX=0. 246 : DO 4106 IB=1,IBOU 247 : ERM=ABS(SIGT(IB)-YK5(IB))/MAX(PRESIG, 248 : . ABS(SIGT(IB)-EINT(IB))) 249 : IF(ERMAX.LT.ERM) ERMAX=ERM 250 : ERM=ABS(W1(IB)-YK5(6+IB))/MAX(PRESIG, 251 : . ABS(W1(IB)-W1INT(IB))) 252 : IF(ERMAX.LT.ERM) ERMAX=ERM 253 : IF(ICENT2.EQ.0) GO TO 4106 254 : ERM=ABS(W2(IB)-YK5(12+IB))/MAX(PRESIG, 255 : . ABS(W2(IB)- W2INT(IB))) 256 : IF(ERMAX.LT.ERM) ERMAX=ERM 257 : 4106 CONTINUE 258 : C 259 : IF(ERMAX.LE.PRERUN) GO TO 1240 260 : C----------------------------------------------------------------------- 261 : C POUR ITER >1 ET IDECO=1 , ON REPREND LE HDEP0 SI ON N'A PAS REUSSI 262 : C A INTEGRER TOUT LE DEPSI EN UN SEUL SOUS-PAS 263 : C----------------------------------------------------------------------- 264 : HDEP=HDEP*0.5D00 265 : IF(JESSAY.GT.1) GO TO 424 266 : IF(IDECO.EQ.1.AND.ITER.GT.1.AND.ABS(HDEP).GT.ABS(HDEP0)) 267 : . HDEP=HDEP0 268 : GO TO 424 269 : 1240 DEDEP=DEDEP+HDEP 270 : C----------------------------------------------------------------------- 271 : C SI ON A FINI , ON SORT EN 514 . SINON ON CONTINUE 272 : C----------------------------------------------------------------------- 273 : IF(ABS(DEDEP).GE.ABS(DEPSI)) GO TO 514 274 : IF(ERMAX.LE.PREMIN) HDEP=2.D00*HDEP 275 : EPSINT=EPST 276 : SINT=SI 277 : DO 502 IB=1,IBOU 278 : EINT(IB)=SIGT(IB) 279 : W1INT(IB)=W1(IB) 280 : IF(ICENT2.EQ.0) GO TO 502 281 : W2INT(IB)=W2(IB) 282 : 502 CONTINUE 283 : C....................................................................... 284 : C FIN DE LA BOUCLE IDECO 285 : C....................................................................... 286 : 416 CONTINUE 287 : C 288 : C ON CONTINUE AVEC UN DEPS DIMINUE 289 : C 290 : DEPS=DEPS-DEPSI+DEDEP 291 : C----------------------------------------------------------------------- 292 : C FIN DE TREANOR 293 : C----------------------------------------------------------------------- 294 : 514 CONTINUE 295 : HDEP0=HDEP 296 : IF(ABS(SI-R)/R.LT.TEST) GO TO 57 297 : IF(ITER.GE.ITMAX) GO TO 56 298 : C 299 : C CALCUL DE LA PENTE ET DU DEPSI 300 : C 301 : ICAS=3 302 : WEP=0.0D00 302 : CALL CHAFON(EPST,AMTRI,AMTRI(7),AMTRI(13),DYK4, 303 : . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,DEPS, 304 : . R0,RM,B,ICAS,WEP,E,W1,W2,PSI,OME,ecou) 305 : DEPSI=(SI-R)/WEP 306 : DEPS=DEPS+DEPSI 307 : IF(DEPS.LT.0.) GO TO 54 308 : GO TO 555 309 : C 310 : 54 DEPS=0.5D00*DEPS0 311 : GO TO 670 312 : C 313 : 56 KERRE=2 314 : 57 JX=JA 315 : DO 301 I=1,IBOU 316 : SIGEL(I)=SIGT(I) 317 : DALPHA(I)=W1(I)-ALPHA1(JX) 318 : IF(ICENT2.EQ.0) GO TO 301 319 : DALPHA(I)=DALPHA(I)+W2(I) 320 : 301 JX=JX+1 321 : C 322 : C LES DEUX CENTRES SONT CUMULES DANS ALPHA1 323 : C 324 : SN=R 325 : C 326 : C MISE A JOUR DES CENTRES DES SPHERES 327 : C 328 : JX=JA 329 : DO 303 I=1,IBOU 330 : ALPHA1(JX)=W1(I) 331 : IF(ICENT2.EQ.0) GO TO 303 332 : ALPHA1(JX)=ALPHA1(JX)+W2(I) 333 : ALPHA2(JX)=W2(I) 334 : 303 JX=JX+1 335 : RETURN 336 : 999 WRITE(6,7999) 337 : 7999 FORMAT('0 CHABOK - CAS NON IMPLEMENTE '/) 338 : RETURN 339 : END 340 : 341 : 342 : 343 : 344 : 345 : 346 :
© Cast3M 2003 - All rights reserved.
Disclaimer