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

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