C CHABOK SOURCE OF166741 25/11/04 21:15:20 12349 SUBROUTINE CHABOK(YUNG,XNU,IA,EI, 1 XMAT,ALPHA1,JA,IBOU,SI,DEPS,EPST,EPSTAR,AMTRI,KERRE, 2 SN,ALPHA2,NUMCHA,ecou,necou) C*********************************************************************** C INTEGRATION MODELE DE CHABOCHE C*********************************************************************** IMPLICIT INTEGER(I-N) IMPLICIT REAL *8(A-H,O-Z) DIMENSION XMAT(*),EI(*),ALPHA1(*),ALPHA2(*),AMTRI(*) DIMENSION YK5(18),DYK1(18),DYK2(18),DYK3(18), & DYK4(18),EINT(6),W1INT(6),W2INT(6) -INC TECOU C C EN SORTIE C KERRE = 2 SI ON ATTEINT LE NOMBRE MAX D'ITERATIONS INTERNES C DATA ITMAX/15/ DATA IDECOU/20/ DATA RAC2/0.707106781D00/ PREDEF=1.D-6 PRESIG=YUNG*PREDEF PREPH=1.D-3 PRERUN= ecou.ECTEST PREMIN=3.D-1*PRERUN C C INITIALISATION A 0. DES VECTEURS POUR TREANOR C NVY=18 DO IB=1,NVY YK5(IB)=0.D0 DYK1(IB)=0.D0 DYK2(IB)=0.D0 DYK3(IB)=0.D0 DYK4(IB)=0.D0 ENDDO A2=0. C2=0. B=0. RM=0. PHI=1.D00 ICOD=0 CALL CHALIM(EPSTAR,R,XMAT,TET,ICOD,A1,C1,A2,C2,R0,RM,B, . PHI,PSI,OME,ICENT2,IDIAM,NUMCHA) ELT=YUNG/(1.+XNU) IF(ITHER.NE.0) ELT=ELT*EI(IA)/YUNG G=ELT*0.5 SAC1=A1*C1 AC1=SAC1*0.66666667 SAC2=A2*C2 AC2=SAC2*0.66666667 SAC12=SAC1+SAC2 AC12=AC1+AC2 GO TO (101,102,102,103,105,102,103,103,109,999, . 999,999,102,101),ITYP 101 E(7)=ELT*(1.-XNU)/(1.-2.*XNU) E(8)=ELT*XNU/(1.-2.*XNU) E(9)=AC12 E(10)=AC1 E(11)=AC2 GO TO 200 102 E(7)=ELT/(1.-XNU) E(8)=E(7)*XNU E(9)=SAC12 E(10)=SAC1 E(11)=SAC2 GO TO 200 103 E(7)=ELT*(1.+XNU) E(8)=0. E(9)=SAC12 E(10)=SAC1 E(11)=SAC2 GO TO 200 105 CONTINUE 109 CONTINUE C C INITIALISATIONS C 200 CONTINUE ITER=0 CALL CHAREM(IDIAM,R0,RM,B,R,EPSTAR,EPST,DEPS,STOT, .W1,W2,E,IBOU,ICENT2,JA,ALPHA1,ALPHA2,PHI,PSI,OME) SI=VONMIS(E,ITYP,ALFAH,COVNMS) C C PREMIER DEPS C ICAS=1 CALL CHAFON(EPST,AMTRI,AMTRI(7),AMTRI(13),DYK1, . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT, . DEPS,R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou) 670 DEPS0=DEPS DEPSI=DEPS IF(ITER.EQ.0) GO TO 666 CALL CHAREM(IDIAM,R0,RM,B,R,EPSTAR,EPST,DBID,STOT, . W1,W2,E,IBOU,ICENT2,JA,ALPHA1,ALPHA2,PHI,PSI,OME) SI=VONMIS(E,ITYP,ALFAH,COVNMS) 666 CONTINUE C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C ITERATIONS INTERNES C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 555 CONTINUE ITER=ITER+1 C C CALCUL PAR LA METHODE DE TREANOR C DO IB=1,IBOU EINT(IB)=E(IB) W1INT(IB)=W1(IB) IF (ICENT2.NE.0) W2INT(IB)=W2(IB) ENDDO SINT=SI DEDEP=0. HDEP=DEPSI IF (ITER.EQ.1) HDEP0=HDEP EPSINT=EPST C----------------------------------------------------------------------- C BOUCLE SUR LES SOUS-PAS C----------------------------------------------------------------------- DO 416 IDECO=1,IDECOU JESSAY=0 DEPRES=DEPSI-DEDEP IF(ABS(HDEP).GT.ABS(DEPRES)) HDEP=DEPRES IF(IDECO.EQ.1) GO TO 417 C 424 CONTINUE JESSAY=JESSAY+1 DO 418 IB=1,IBOU E(IB)=EINT(IB) W1(IB)=W1INT(IB) IF(ICENT2.EQ.0) GO TO 418 W2(IB)=W2INT(IB) 418 CONTINUE SI=SINT 417 ICAS=2 WEP=0.5D00 CALL CHAFON( EPSINT ,AMTRI,AMTRI(7),AMTRI(13),DYK1, . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,HDEP, . R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou) C EPST=EPSINT+0.5D00*HDEP CALL CHAFON( EPST ,AMTRI,AMTRI(7),AMTRI(13),DYK2, . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,HDEP, . R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou) C WEP=1.D00 CALL CHAFON( EPST ,AMTRI,AMTRI(7),AMTRI(13),DYK3, . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,HDEP, . R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou) C PHMAX=0. DO 4101 IB=1,NVY IF(ABS(DYK2(IB)-DYK1(IB)).LT.PRESIG) GO TO 4101 PH=2.*(DYK2(IB)-DYK3(IB))/(DYK2(IB)-DYK1(IB)) IF(PHMAX.LT.PH) PHMAX=PH 4101 CONTINUE IF(PHMAX.LE.0.) GO TO 4102 IF(PHMAX.GE.PREPH) GO TO 4123 XL1=1.-PHMAX*0.5 XL2=0.5-PHMAX*0.166666666666666667 XL3=0.166666666666666667-PHMAX*0.04166666666666666667 XL4=-PHMAX*0.166666666666666667 GO TO 4103 4123 CONTINUE XL1=(1.-EXP(-PHMAX))/PHMAX XL2=(1.-XL1)/PHMAX XL3=(0.5-XL2)/PHMAX XL4=XL1-2.*XL2 GO TO 4103 4102 PHMAX=0. XL1=1. XL2=0.5 XL3=0.166666666666666667 XL4=0. 4103 CONTINUE XM1=-XL2+4.*XL3 XM2=2.*(XL2-2.*XL3) XM3=4.*XL3-3.*XL2 C C NOUVELLE ESTIMATION DE LA SOLUTION C DO 4104 IB=1,IBOU E(IB)=EINT(IB)+2.*XL2*DYK3(IB)+XL2*PHMAX*DYK2(IB)+XL4*DYK1(IB) YK5(IB)=E(IB) W1(IB)=W1INT(IB)+2.*XL2*DYK3(6+IB) . +XL2*PHMAX*DYK2(6+IB)+XL4*DYK1(6+IB) YK5(6+IB)=W1(IB) IF(ICENT2.EQ.0) GO TO 4104 W2(IB)=W2INT(IB)+2.*XL2*DYK3(12+IB) . +XL2*PHMAX*DYK2(12+IB)+XL4*DYK1(12+IB) YK5(12+IB)=W2(IB) 4104 CONTINUE SI=VONMIS(E,ITYP,ALFAH,COVNMS) C EPST=EPSINT+HDEP WEP=1.D00 CALL CHAFON( EPST ,AMTRI,AMTRI(7),AMTRI(13),DYK4, . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,HDEP, . R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou) C C ESTIMATION FINALE DE LA SOLUTION C FAC1=1.+PHMAX*XM3+2.*PHMAX*XM2 FAC2=XL1+XM3+0.5*XM2*PHMAX FAC3=XM2*(1.+0.5*PHMAX) DO 4105 IB=1,IBOU E(IB)=EINT(IB)*FAC1+DYK1(IB)*FAC2+DYK2(IB)*FAC3 . +DYK3(IB)*XM2+DYK4(IB)*XM1+YK5(IB)*XM1*PHMAX SIGT(IB)=E(IB) W1(IB)=W1INT(IB)*FAC1+DYK1(6+IB)*FAC2 . +DYK2(6+IB)*FAC3 +DYK3(6+IB)*XM2 +DYK4(6+IB)*XM1 . +YK5(6+IB)*XM1*PHMAX IF(ICENT2.EQ.0) GO TO 4105 W2(IB)= W2INT(IB)*FAC1 +DYK1(12+IB)*FAC2 . +DYK2(12+IB)*FAC3 +DYK3(12+IB)*XM2 +DYK4(12+IB)*XM1 . +YK5(12+IB)*XM1*PHMAX 4105 CONTINUE SI=VONMIS(SIGT,ITYP,ALFAH,COVNMS) C C CALCUL DE L'ERREUR C ERMAX=0. DO 4106 IB=1,IBOU ERM=ABS(SIGT(IB)-YK5(IB))/MAX(PRESIG, . ABS(SIGT(IB)-EINT(IB))) IF(ERMAX.LT.ERM) ERMAX=ERM ERM=ABS(W1(IB)-YK5(6+IB))/MAX(PRESIG, . ABS(W1(IB)-W1INT(IB))) IF(ERMAX.LT.ERM) ERMAX=ERM IF(ICENT2.EQ.0) GO TO 4106 ERM=ABS(W2(IB)-YK5(12+IB))/MAX(PRESIG, . ABS(W2(IB)- W2INT(IB))) IF(ERMAX.LT.ERM) ERMAX=ERM 4106 CONTINUE C IF(ERMAX.LE.PRERUN) GO TO 1240 C----------------------------------------------------------------------- C POUR ITER >1 ET IDECO=1 , ON REPREND LE HDEP0 SI ON N'A PAS REUSSI C A INTEGRER TOUT LE DEPSI EN UN SEUL SOUS-PAS C----------------------------------------------------------------------- HDEP=HDEP*0.5D00 IF(JESSAY.GT.1) GO TO 424 IF(IDECO.EQ.1.AND.ITER.GT.1.AND.ABS(HDEP).GT.ABS(HDEP0)) . HDEP=HDEP0 GO TO 424 1240 DEDEP=DEDEP+HDEP C----------------------------------------------------------------------- C SI ON A FINI , ON SORT EN 514 . SINON ON CONTINUE C----------------------------------------------------------------------- IF(ABS(DEDEP).GE.ABS(DEPSI)) GO TO 514 IF(ERMAX.LE.PREMIN) HDEP=2.D00*HDEP EPSINT=EPST SINT=SI DO 502 IB=1,IBOU EINT(IB)=SIGT(IB) W1INT(IB)=W1(IB) IF(ICENT2.EQ.0) GO TO 502 W2INT(IB)=W2(IB) 502 CONTINUE C....................................................................... C FIN DE LA BOUCLE IDECO C....................................................................... 416 CONTINUE C C ON CONTINUE AVEC UN DEPS DIMINUE C DEPS=DEPS-DEPSI+DEDEP C----------------------------------------------------------------------- C FIN DE TREANOR C----------------------------------------------------------------------- 514 CONTINUE HDEP0=HDEP IF(ABS(SI-R)/R.LT.ECTEST) GO TO 57 IF(ITER.GE.ITMAX) GO TO 56 C C CALCUL DE LA PENTE ET DU DEPSI C ICAS=3 WEP=0.0D00 CALL CHAFON(EPST,AMTRI,AMTRI(7),AMTRI(13),DYK4, . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,DEPS, . R0,RM,B,ICAS,WEP,E,W1,W2,PSI,OME,ecou) DEPSI=(SI-R)/WEP DEPS=DEPS+DEPSI IF(DEPS.LT.0.) GO TO 54 GO TO 555 C 54 DEPS=0.5D00*DEPS0 GO TO 670 C 56 KERRE=2 57 JX=JA DO I=1,IBOU SIGEL(I)=SIGT(I) DALPHA(I)=W1(I)-ALPHA1(JX) IF (ICENT2.NE.0) DALPHA(I)=DALPHA(I)+W2(I) JX=JX+1 ENDDO C C LES DEUX CENTRES SONT CUMULES DANS ALPHA1 C SN=R C C MISE A JOUR DES CENTRES DES SPHERES C JX=JA DO I=1,IBOU ALPHA1(JX)=W1(I) IF (ICENT2.NE.0) THEN ALPHA1(JX)=ALPHA1(JX)+W2(I) ALPHA2(JX)=W2(I) ENDIF JX=JX+1 ENDDO RETURN 999 WRITE(6,7999) 7999 FORMAT('0 CHABOK - CAS NON IMPLEMENTE '/) RETURN END