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

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