Télécharger la plaquette Cast3M


fr en

Télécharger chabok.eso

Retour à la liste

Numérotation des lignes : 

chabok

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales