Télécharger cmodcl.eso

Retour à la liste

cmodcl
C CMODCL    SOURCE    PV        17/12/08    21:16:11     9660           
      SUBROUTINE CMODCL(WRK52,NMATT1,EPSP1,EPSP2,YLIM1,YLIM2,Z1,Z2,Y1,
     &                  Y2,JSIGN,SIG1,SIG2,CONV)
C MODCLB    SOURCE    AM        98/12/23    21:40:04     3409
C
C
C     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C     %                                                               %
C     %  INTEGRATION DU MODCLB UNILATERAL AVEC FERMETURE LINEAIRE     %
C     %          DE FISSURES EN CONTRAINTES PLANES                    %
C     %                (MODCLB NON ASSOCIE)                           %
C     %      A L"AIDE DE LA METHODE DE NEWTON                         %
C     %                                                               %
C     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
*      SUBROUTINE MODCLB(WRK0,NMATT,EPSP1,EPSP2,YLIM1,YLIM2,Z1,Z2,Y1,Y2,
*     1                  JSIGN,SIG1,SIG2,CONV)
C
C
C        VARIABLES D"ENTREES
C
      IMPLICIT INTEGER(I-N)
-INC DECHE
C
      REAL*8 EPSP1,EPSP2
C
C        VARIABLES DE SORTIE
C
C
      REAL*8 SIG1,SIG2,Z1,Z2,Y1,Y2
      INTEGER JSIGN(3)
      LOGICAL CONV
C
C
C        PARAMETRES D"EVOLUTION
C
C
      REAL*8 E,NU,A1,A2,B1,B2,Y01,Y02,SIMF
      REAL*8 BETA1,BETA2,YLIM1,YLIM2
C
C
C        VARIABLES INTERNES
C
C
      REAL*8 DELSIG
      INTEGER NDIM
      PARAMETER (NDIM=6)
      REAL*8 JACOB(NDIM,NDIM),J15,J16,J25,J26
      REAL*8 DELAFO(NDIM)
      REAL*8 X(NDIM),DX(NDIM),X0(NDIM),XX,XIG1,XIG2
      REAL*8 ERREUR,TEST
      REAL*8 H1,H2,HT,HT2,trf
      REAL*8 ALFA1,ALFA2
      REAL*8 SIGY
      REAL*8 UN,DEUX
      REAL*8 ZMAX
      INTEGER ITER ,I,DIMPB
      LOGICAL LDOM1,LDOM2
 
      UN=1.D0
      DEUX=2.D0
      ZMAX=1.e+04
      E  = XMAT(1)
      NU   = XMAT(2)
      Y01   = XMAT(5)
      Y02   = XMAT(6)
      A1    = XMAT(7)
      A2    = XMAT(8)
      B1    = XMAT(9)
      B2    = XMAT(10)
      BETA1 = XMAT(11)
      BETA2 = XMAT(12)
      SIMF  = XMAT(13)
      SIGY  =SQRT(2.D0*E*Y01)
      DELSIG=SIGY*1.D-2
*      DELSIG=0.d0
*      print*,'E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2'
*      print*,E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2
      TEST=1.D-12
      ERREUR=1.D0
      ITER=0
      X(1)=SIG1
      X(2)=SIG2
      X(3)=Y1
      X(4)=Y2
      X(5)=Z1
      X(6)=Z2
      DO I=1, NDIM
         X0(I)=X(I)
      END DO
C
C
C
      JACOB(1,1) = -UN
      JACOB(1,2) = 0.D0
      JACOB(1,3) = 0.D0
      JACOB(1,4) = 0.D0
      JACOB(2,1) = 0.D0
      JACOB(2,2) = -UN
      JACOB(2,3) = 0.D0
      JACOB(2,4) = 0.D0
      JACOB(3,3) = -UN
      JACOB(3,4) = 0.D0
      JACOB(3,6) = 0.D0
      JACOB(4,3) = 0.D0
      JACOB(4,4) = -UN
      JACOB(4,5) = 0.D0
      JACOB(5,1) = 0.D0
      JACOB(5,2) = 0.D0
      JACOB(5,4) = 0.D0
      JACOB(5,5) = -UN
      JACOB(5,6) = 0.D0
      JACOB(6,1) = 0.D0
      JACOB(6,2) = 0.D0
      JACOB(6,3) = 0.D0
      JACOB(6,5) = 0.D0
      JACOB(6,6) = -UN
      alfa1=UN
      alfa2=UN
 
      DO WHILE (ERREUR .GT. TEST .AND. ITER .LE. 50)
         ITER=ITER+1
         IF (X(1) .GE. DELSIG) THEN
            H1=UN
         ELSE IF (X(1) .LE. -DELSIG) THEN
            H1=0.D0
         ELSE
            H1=(X(1)+DELSIG)/(DEUX*DELSIG)
         END IF
         IF (X(2) .GE. DELSIG) THEN
            H2=UN
         ELSE IF (X(2) .LE. -DELSIG) THEN
            H2=0.D0
         ELSE
            H2=(X(2)+DELSIG)/(DEUX*DELSIG)
         END IF
c
c     smooth transition on ht ht2
c
         trf = x(1) + x(2)
         IF(trf.GE.delsig)THEN
            HT=UN
            HT2=0.D0
         else if(trf.ge.(-delsig)) then
            ht = (trf+delsig)/DEUX/delsig
            ht2 = UN-ht
         else if(trf.ge.(-simf+delsig)) then
            ht = 0.D+00
            ht2 = UN
         else if(trf.ge.(-simf-delsig)) then
            ht = 0.D+00
c           ht2 =UN- (trf-(simf-delsig))/(DEUX*simf)
            ht2 = (trf+simf+delsig)/DEUX/delsig
         else
            ht = 0.D+00
            ht2 = 0.D+00
         endif
         LDOM1= (X(3) .GT. YLIM1 .AND. X(3).GT. Y01
     1                           .AND. X(5) .LT. ZMAX)
         LDOM2= (X(4) .GT. YLIM2 .AND. X(4).GT. Y02
     1                           .AND. X(6) .LT. ZMAX)
*         print*,'y1',x(3),'ylim1',ylim1
*         print*,'ldom1',ldom1,'ldom2',ldom2
C        =====> PROBLEME ELASTIQUE
C
C
C        CALCUL DU JACOBIEN
C
C
         JACOB(3,1) = H1*X(1)*X(5)*X(5)/E
C
         JACOB(3,2) = H2*X(2)*X(5)*X(5)/E
C
         JACOB(4,1) = (UN-H1)*X(1)*X(6)*X(6)/E
C
         JACOB(4,2) = (UN-H2)*X(2)*X(6)*X(6)/E
C
C        CALCUL DU SECOND MEMBRE
C
      XIG1 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN)
     1  *BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)
     2  *HT2+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+
     3  BETA2*(X(6)*(X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1*
     4  (X(6)*(H2-UN)-X(5)*H2))*SIMF+(BETA1*E*(X(5)*EPSP2-EPSP2)+
     5  (UN-X(5))*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-UN)-
     6  X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIMF+
     7  (DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*
     8  (-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))
     9  *BETA1*H1)*HT2)
C
      XIG2 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN)
     1  *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)*
     2  X(6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+
     3  (X(5)-UN)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+
     4  (X(5)+UN)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))*
     5  BETA2)*SIMF+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1)
     6  *HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*
     7  (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIMF+(DEUX-DEUX*X(5))*
     8  BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*
     9  (H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)
 
         DELAFO(1) =X(1)-XIG1
C
         DELAFO(2) =X(2)-XIG2
C
         XX= (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)*X(5)/(DEUX*E)
         DELAFO(3) = X(3)-XX
C
         XX=((UN-H1)*X(1)*X(1)+(UN-H2)*X(2)*X(2) )
     1      *X(6)*X(6)/(DEUX*E)
         DELAFO(4) = X(4)-XX
C
         IF (LDOM1) THEN
C
C        =====> ENDOMMAGEMENT EN TRACTION
C
      J15 = (((BETA1*HT2+BETA1*HT)*NU+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2
     &   ))*HT2+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2))*HT+BETA2*(X(6)*H2-H
     &   2)-E*EPSP1*H2)*SIMF+(BETA1*E*EPSP2-BETA1*E*EPSP1)*HT2)/((NU**2+
     &   H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6
     &   )*H2+X(6)**2*(H2-UN))*SIMF+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1
     &   *(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN
     &   )*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX*H2-UN)-DEU
     &   X*X(5)*H2)-X(6)*H2)*SIMF-DEUX*BETA1*HT2*NU+(BETA1*(-DEUX*X(5)*H
     &   2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*HT2)*((((X(5
     &   )-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN)*BETA2)*NU+
     &   BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT2+BETA1*
     &   (X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+BETA2*(X(6)*(
     &   X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1*(X(6)*(H2-UN)-X
     &   (5)*H2))*SIMF+(BETA1*E*(X(5)*EPSP2-EPSP2)+(UN-X(5))*BETA1*E*EPS
     &   P1)*HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*
     &   (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIMF+(DEUX-DEUX*X(5))*BE
     &   TA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)
     &   +DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
C
C
      J25 = (((BETA1*HT2+BETA1*HT)*NU+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6
     &   )*BETA1)*HT2+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6)*BETA1)*HT+((X(
     &   6)-UN)*BETA2-E*EPSP2)*H1)*SIMF+(BETA1*E*EPSP1-BETA1*E*EPSP2)*HT
     &   2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2
     &   ))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIMF+(DEUX-DEUX*X(5))*BETA1*HT
     &   2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX)
     &   )+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX
     &   *H2-UN)-DEUX*X(5)*H2)-X(6)*H2)*SIMF-DEUX*BETA1*HT2*NU+(BETA1*(-
     &   DEUX*X(5)*H2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*H
     &   T2)*((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN)
     &   *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)*X(
     &   6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)
     &   *X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+(X(5)+UN)*
     &   X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))*BETA2)*SIMF+(B
     &   ETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1)*HT2)/((NU**2
     &   +H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(
     &   6)*H2+X(6)**2*(H2-UN))*SIMF+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA
     &   1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-U
     &   N)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
C
            JACOB(1,5) = J15
            JACOB(2,5) = J25
            JACOB(3,5) = (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)/E
            JACOB(4,5) = 0.D0
            JACOB(5,5) = -1.D0
            JACOB(5,1) = 0.D0
            JACOB(5,2) = 0.D0
            JACOB(5,3) = A1*B1*(A1*(X(3)-Y01))**(B1-UN)
            JACOB(5,4) = 0.D0
            DELAFO(5) = X(5)-((A1*(X(3)-Y01))**B1+UN)
            IF (LDOM2) THEN
C
C        =====> ENDOMMAGEMENT EN COMPRESSION
C
      J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2)
     1   -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIMF/
     2   ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)*
     3   X(6)*H2+X(6)**2*(H2-1))*SIMF+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X
     4   (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5)
     5   **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5
     6   )*H2+2*X(6)*(H2-1))*SIMF+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE
     7   TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X
     8   (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)
     9   *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+
     :   BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6
     ;   )*(H2-1)-X(5)*H2))*SIMF+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE
     <   TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6
     =   )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIMF+(2-2*X(5))*BETA1
     >   *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+(
     ?   (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
C
      J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B
     1   ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1-
     2   E*EPSP2+(2*X(6)-1)*BETA2)*SIMF/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X
     3   (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIMF+(2-
     4   2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)
     5   *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X
     6   (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIMF+(BETA1*
     7   (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H
     8   T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6)
     9   +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6
     :   )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2
     ;   -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS
     <   P2+(X(6)**2-X(6))*BETA2)*SIMF+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)
     =   -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*
     >   H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIMF+(2-2*X(5))
     ?   *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)
     @   +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
C
C
               JACOB(1,6) = J16
               JACOB(2,6) = J26
               JACOB(3,6) = 0.D0
               JACOB(4,6) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2)
     &                      *X(2)*X(2))*X(6)/E
               JACOB(5,6) = 0.D0
               JACOB(6,6) = -1.D0
               JACOB(6,1) = 0.D0
               JACOB(6,2) = 0.D0
               JACOB(6,3) = 0.D0
*               JACOB(6,4) = (A2*B2*(X(4)-Y02)**(B2-1.0))
               JACOB(6,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0)
               JACOB(6,5) = 0.D0
*               DELAFO(6) = X(6)-(A2*(X(4)-Y02)**B2+1.0)
               DELAFO(6) = X(6)-((A2*(X(4)-Y02))**B2+1.0)
               DIMPB=6
            ELSE
               DIMPB=5
            END IF
         ELSE
C
C        =====> PAS D'ENDOMMAGEMENT EN TRACTION
C
            IF (LDOM2) THEN
C
C        =====> ENDOMMAGEMENT EN COMPRESSION
C
      J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2)
     1   -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIMF/
     2   ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)*
     3   X(6)*H2+X(6)**2*(H2-1))*SIMF+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X
     4   (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5)
     5   **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5
     6   )*H2+2*X(6)*(H2-1))*SIMF+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE
     7   TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X
     8   (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)
     9   *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+
     :   BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6
     ;   )*(H2-1)-X(5)*H2))*SIMF+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE
     <   TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6
     =   )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIMF+(2-2*X(5))*BETA1
     >   *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+(
     ?   (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
C
      J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B
     1   ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1-
     2   E*EPSP2+(2*X(6)-1)*BETA2)*SIMF/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X
     3   (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIMF+(2-
     4   2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)
     5   *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X
     6   (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIMF+(BETA1*
     7   (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H
     8   T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6)
     9   +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6
     :   )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2
     ;   -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS
     <   P2+(X(6)**2-X(6))*BETA2)*SIMF+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)
     =   -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*
     >   H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIMF+(2-2*X(5))
     ?   *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)
     @   +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
C
C
 
               JACOB(1,5) = J16
               JACOB(2,5) = J26
               JACOB(3,5) = 0.D0
               JACOB(4,5) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2)
     &                      *X(2)*X(2))*X(6)/E
               JACOB(5,5) = -1.D0
               JACOB(5,1) = 0.D0
               JACOB(5,2) = 0.D0
               JACOB(5,3) = 0.D0
*               JACOB(5,4) = (A2*B2*(X(4)-Y02)**(B2-1.0))
               JACOB(5,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0)
               JACOB(5,5) = -1.D0
*               DELAFO(5) = X(6)-(A2*(X(4)-Y02)**B2+1.0)
               DELAFO(5) = X(6)-((A2*(X(4)-Y02))**B2+1.0)
               DIMPB=5
            ELSE
               DIMPB=4
            END IF
         END IF
C
C        RESOLUTION
C
         CALL GAUSCL(JACOB,DX,DELAFO,NDIM,DIMPB,0)
C
C        CALCUL DE L"ERREUR ET DE LA NOUVELLE SOLUTION
C
         ERREUR=0.D0
         DO I=1,2
            ERREUR=ERREUR+DX(I)*DX(I)/(SIGY*SIGY)
            X(I)=X(I)+DX(I)
         END DO
         X(3)=X(3)+DX(3)
         ERREUR=ERREUR+DX(3)*DX(3)/(Y01*Y01)
         X(4)=X(4)+DX(4)
         ERREUR=ERREUR+DX(4)*DX(4)/(Y02*Y02)
         IF (LDOM1) THEN
           IF ((X(5)+DX(5)).LT.Z1) THEN
              DX(5)=Z1-X(5)
           END IF
           IF ( X(5) .GT. ZMAX)THEN
             DX(5)=ZMAX + 0.1D0 - X(5)
           END IF
           X(5)=X(5)+DX(5)
           ERREUR=ERREUR+DX(5)*DX(5)/(X(5)*X(5))
           IF (LDOM2) THEN
           IF ((X(6)+DX(6)).LT.Z2) THEN
              DX(6)=Z2-X(6)
           END IF
           IF ( X(6) .GT. ZMAX)THEN
             DX(6)=ZMAX + 0.1D0 - X(6)
           END IF
           X(6)=X(6)+DX(6)
           ERREUR=ERREUR+DX(6)*DX(6)/(X(6)*X(6))
           END IF
         ELSE IF (LDOM2) THEN
           IF ((X(6)+DX(5)).LT.Z2) THEN
              DX(5)=Z2-X(6)
           END IF
           IF ( X(6) .GT. ZMAX)THEN
             DX(5)=ZMAX + 0.1D0 - X(6)
           END IF
           X(6)=X(6)+DX(5)
           ERREUR=ERREUR+DX(5)*DX(5)/(X(6)*X(6))
         END IF
         IF (ERREUR .GT. 1.D+100)THEN
            CONV=.FALSE.
           print*,'erreur dans modele',erreur
           print*,'iteration num',iter
           print*,'delafo',delafo
            RETURN
         END IF
C
C        QUELQUES IMPRESSIONS EN CAS D'ERREUR
C
         if (.not. conv )then
            PRINT*,'ITERATION  ',ITER
            print*,'LDOM1  :',LDOM1,'   LDOM2  :',LDOM2
            PRINT*,'EPS1= ',EPSP1,'  EPS2= ',EPSP2
            PRINT*,'SIG1:',X(1) ,'   DSIG1:',DX(1),'ERR :',DELAFO(1)
            PRINT*,'SIG2:',X(2) ,'   DSIG2:',DX(2),'ERR :',DELAFO(2)
            PRINT*,'Y1  :',X(3) ,'   DY1  :',DX(3),'ERR :',DELAFO(3)
            PRINT*,'Y2  :',X(4) ,'   DY2  :',DX(4),'ERR :',DELAFO(4)
            PRINT*,'Z1  :',X(5) ,'   DZ1  :',DX(5),'ERR :',DELAFO(5)
            PRINT*,'Z2  :',X(6) ,'   DZ2  :',DX(6),'ERR :',DELAFO(6)
            PRINT*,'H1',H1,'H2',H2,ht,ht2
            PRINT*
            write(*,200)jacob
200         format(6(6(2x,e11.4),/))
         END IF
 
      END DO
      IF (ERREUR .GT. TEST ) THEN
         CONV=.FALSE.
      ELSE
         CONV=.TRUE.
         SIG1= X(1)
         SIG2= X(2)
         Y1    = X(3)
         Y2    = X(4)
         Z1    = X(5)
         Z2    = X(6)
         IF (SIG1 .GE. 0.D0) THEN
            JSIGN(1)=1
         ELSE
            JSIGN(1)=0
         END IF
         IF (SIG2.GE. 0.D0)THEN
            JSIGN(2)=1
         ELSE
            JSIGN(2)=0
         END IF
         JSIGN(3)=1
      END IF
      RETURN
      END
 
 
 
 
 
 
 
 
 
 
 
 
 

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