Sources Esope | Cast3M

ecubic
C ECUBIC    SOURCE    CHAT      05/01/12    23:28:10     5004
      SUBROUTINE ECUBIC (SOGMA,DSOGMA,COHE1,COHE2,CO11,SI11,
     1 CO22,SI22,CO21,SI21,TRA1,TRA2,YOUNG,XNU,DEFPL,TSUG,IDAM,
     2 DLAM1,DLAM2,PSI1,PHI1,PSI2,PHI2,NCAS,ICRI,HACHE1,HACHE2,KERRE)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION DSOGMA(4),SOGMA(4),TSOG(4),TSOG1(4),TSOG2(4),
     1 SIP1(4),SAP1(4),SAP(4),SIP(4),SIP2(4),SAP2(4),TSUG(4),
     2 DSUG(4),SOG1(4),SOG2(4),DSOG1(4),DSOG2(4),A(3,3),B(3),
     3 C(3),D(3),DEF(4),DEFPL(4),SAG1(4),SAG2(4),SEG(4)
      DO 100 ITYP=1,4
      TSOG(ITYP)=DSOGMA(ITYP)+SOGMA(ITYP)
      SEG(ITYP)=TSOG(ITYP)
  100 CONTINUE
      ITER=1
      ICHRI=0
      ICHRO=0
      ZER=0.D0
      IDAM=12
      ISOM1=0
      ISOM2=0
      DLA1=0.D0
      DLA2=0.D0
C------------------------------------------------------------------
C ON EFFECTUE LE COUPLAGE
C------------------------------------------------------------------
C
      CALL CHREPE(CO11,SI11,TSOG,TSOG1)
      CALL CHREPE(CO22,SI22,TSOG,TSOG2)
      CALL CHREPE(CO11,SI11,SOGMA,SOG1)
      CALL CHREPE(CO22,SI22,SOGMA,SOG2)
      CALL CHREPE(CO11,SI11,DSOGMA,DSOG1)
      CALL CHREPE(CO22,SI22,DSOGMA,DSOG2)
C---------------------------------------------------------------
C     SI ON EST A LA POINTE CAS PRATICULIER
C     IL FAUT REGARDER SI ON SORT DANS LE CONE DES NORMALES
C---------------------------------------------------------------
C
      IF(SOG1(2).GT.(1.D0-1.D-6)*TRA1) ISOM1=1
      IF(SOG2(2).GT.(1.D0-1.D-6)*TRA2) ISOM2=1
      IF((ISOM1+ISOM2).EQ.0) GO TO 90000
C
      CIS1=SIGN(1.D0,SOG1(4))
      CIS2=SIGN(1.D0,SOG2(4))
      IF(SOG1(2).GT.(1.D0-1.D-6)*TRA1) CIS1=SIGN(1.D0,TSOG1(4))
      IF(SOG2(2).GT.(1.D0-1.D-6)*TRA2) CIS2=SIGN(1.D0,TSOG2(4))
      VCRA=DSOG1(2)-PSI1*ABS(DSOG1(4))
      VCRB=DSOG2(2)-PSI2*ABS(DSOG2(4))
      IF(ISOM1.EQ.1.AND.VCRA.GT.0.D0) ICHRO=1
      IF(ISOM2.EQ.1.AND.VCRB.GT.0.D0) ICHRI=1
      IF(ISOM1.EQ.1.AND.VCRA.GT.0.D0) GO TO 20000
      IF(ISOM2.EQ.1.AND.VCRB.GT.0.D0) GO TO 20000
      GOTO 80000
C
90000 CIS1=SIGN(1.D0,SOG1(4))
      CIS3=SIGN(1.D0,TSOG1(4))
      CIS2=SIGN(1.D0,SOG2(4))
      CIS4=SIGN(1.D0,TSOG2(4))
80000 IF(NCAS.EQ.1) CALL CPHOOK(ZER,PSI1,ZER,CIS1,SIP1(1),SIP1(2)
     1 ,SIP1(3),SIP1(4),XNU,YOUNG)
      ICHRO=0
      ICHRI=0
      IF(NCAS.EQ.1) CALL CPHOOK(ZER,PSI2,ZER,CIS2,SAP2(1),SAP2(2)
     1 ,SAP2(3),SAP2(4),XNU,YOUNG)
      IF(NCAS.NE.1) CALL DPHOOK(ZER,PSI1,ZER,CIS1,SIP1(1),SIP1(2)
     1 ,SIP1(3),SIP1(4),XNU,YOUNG)
      IF(NCAS.NE.1) CALL DPHOOK(ZER,PSI2,ZER,CIS2,SAP2(1),SAP2(2)
     1 ,SAP2(3),SAP2(4),XNU,YOUNG)
      CALL CHREPE(CO21,SI21,SIP1,SIP2)
      CALL CHREPE(CO21,-SI21,SAP2,SAP1)
      CALL CHREPE(CO11,-SI11,SIP1,SIP)
      CALL CHREPE(CO22,-SI22,SAP2,SAP)
      ALPH1=PHI1*SIP1(2)+CIS1*SIP1(4)
      ALPH2=PHI1*SAP1(2)+CIS1*SAP1(4)
      ALPH3=-COHE1+PHI1*TSOG1(2)+CIS1*TSOG1(4)
      ALPH4=PHI2*SIP2(2)+CIS2*SIP2(4)
      ALPH5=PHI2*SAP2(2)+CIS2*SAP2(4)
      ALPH6=-COHE2+PHI2*TSOG2(2)+CIS2*TSOG2(4)
      DET=ALPH1*ALPH5-ALPH4*ALPH2
      REF=YOUNG*YOUNG*1.D-7
      IF(DET.LE.REF) GO TO 30000
C------------------------------------------------------------------
C SI LE DETERMINANT EST NUL ON NE PEUT EFFECTUER LE COUPLAGE
C ON ECOULE DONC SUIVANT UN SEUL CRITERE  (20000)
C------------------------------------------------------------------
C
      DLA1=(ALPH3*ALPH5-ALPH6*ALPH2)/DET
      DLA2=(ALPH1*ALPH6-ALPH3*ALPH4)/DET
      DO 110 ITYP=1,4
      DSUG(ITYP)=DLA1*SIP(ITYP)+DLA2*SAP(ITYP)
      TSUG(ITYP)=TSOG(ITYP)-DSUG(ITYP)
      DSUG(ITYP)=SEG(ITYP)-TSUG(ITYP)
  110 CONTINUE
      CALL CHREPE(CO11,SI11,TSUG,SAG1)
      CALL CHREPE(CO22,SI22,TSUG,SAG2)
C-----------------------------------------------------
C   DANS LE CAS OU ON DEPASSE LA LIMITE EN TRACTION
C   POUR LE CRITERE 1 OU 2 ON CALCULE DES NOUVEAUX
C   DLAMDA
C-----------------------------------------------------
      IF(SAG1(2).GT.TRA1*(1.D0+1.D-6).OR.SAG2(2).GT.TRA2*(1.D0+1.D-6))
     1 GO TO 500
      DLAM1=DLA1+DLAM1
      DLAM2=DLA2+DLAM2
      IF(NCAS.EQ.1) CALL CPCONT(DSUG,DEFPL,XNU,YOUNG)
      IF(NCAS.NE.1) CALL DPCONT(DSUG,DEFPL,XNU,YOUNG)
      RETURN
  500 C(3)=100.D0
      D(3)=100.D0
      IF(SAG1(2).LE.(TRA1*(1.D0+1.D-6))) GO TO 600
      B(1)=-COHE1+PHI1*SOG1(2)+CIS1*SOG1(4)
      B(2)=-COHE2+PHI2*SOG2(2)+CIS2*SOG2(4)
      B(3)=SOG1(2)-TRA1
      A(1,3)=-(PHI1*DSOG1(2)+CIS1*DSOG1(4))
      A(2,3)=-(PHI2*DSOG2(2)+CIS2*DSOG2(4))
      A(3,3)=-DSOG1(2)
      A(1,1)=ALPH1
      A(1,2)=ALPH2
      A(2,1)=ALPH4
      A(2,2)=ALPH5
      A(3,1)=SIP1(2)
      A(3,2)=SAP1(2)
      CALL RESNEQ(A,B,C,3,KERRE)
      IF(SAG2(2).LE.(TRA2*(1.D0-1.D-6))) GO TO 800
 600  B(1)=COHE1-PHI1*SOG1(2)-CIS1*SOG1(4)
      B(2)=COHE2-PHI2*SOG2(2)-CIS2*SOG2(4)
      B(3)=SOG2(2)-TRA2
      A(1,3)=-(PHI1*DSOG1(2)+CIS1*DSOG1(4))
      A(2,3)=-(PHI2*DSOG2(2)+CIS2*DSOG2(4))
      A(3,3)=-DSOG2(2)
      A(1,1)=ALPH1
      A(1,2)=ALPH2
      A(2,1)=ALPH4
      A(2,2)=ALPH5
      A(3,1)=SIP2(2)
      A(3,2)=SAP2(2)
      CALL RESNEQ(A,B,D,3,KERRE)
 800  IF(D(3).LE.C(3)) DLA1=D(1)
      IF(D(3).GT.C(3)) ICHRO=1
      IF(D(3).LE.C(3)) ICHRI=1
      IF(ABS(D(3)-C(3)).LE.(ABS(C(3))*1.D-6)) ICHRI=1
      IF(ABS(D(3)-C(3)).LE.(ABS(C(3))*1.D-6)) ICHRO=1
      IF(D(3).LE.C(3)) DLA2=D(2)
      IF(D(3).LE.C(3)) XAT=D(3)
      IF(D(3).GT.C(3)) DLA1=C(1)
      IF(D(3).GT.C(3)) DLA2=C(2)
      IF(D(3).GT.C(3)) XAT=C(3)
      DO 111 ITYP=1,4
      DSUG(ITYP)=DLA1*SIP(ITYP)+DLA2*SAP(ITYP)
      SOGMA(ITYP)=SOGMA(ITYP)-DSUG(ITYP)+XAT*DSOGMA(ITYP)
      DSOGMA(ITYP)=(1.D0-XAT)*DSOGMA(ITYP)
      TSOG(ITYP)=DSOGMA(ITYP)+SOGMA(ITYP)
  111 CONTINUE
      DLAM1=DLAM1+DLA1
      DLAM2=DLAM2+DLA2
      CALL CHREPE(CO11,SI11,TSOG,TSOG1)
      CALL CHREPE(CO22,SI22,TSOG,TSOG2)
      CALL CHREPE(CO11,SI11,SOGMA,SOG1)
      VCRIT1=VCRITE(SOG1(2),SOG1(4),PHI1,COHE1)
      CALL CHREPE(CO22,SI22,SOGMA,SOG2)
      VCRIT2=VCRITE(SOG2(2),SOG2(4),PHI2,COHE2)
      CALL CHREPE(CO22,SI22,DSOGMA,DSOG2)
      CALL CHREPE(CO11,SI11,DSOGMA,DSOG1)
      VCRA=DSOG1(2)-PSI1*ABS(DSOG1(4))
      VCRB=DSOG2(2)-PSI2*ABS(DSOG2(4))
      ITER=ITER+1
      IF(ICHRO.EQ.1.AND.VCRA.GT.0.D0) GO TO 20000
      IF(ICHRI.EQ.1.AND.VCRB.GT.0.D0) GO TO 20000
      IF(ITER.EQ.7) GO TO 20000
C
      IF(ICHRO.EQ.1) CIS1=SIGN(1.D0,TSOG1(4))
      IF(ICHRO.NE.1) CIS1=SIGN(1.D0,SOG1(4))
      IF(ICHRI.EQ.1) CIS2=SIGN(1.D0,TSOG2(4))
      IF(ICHRI.NE.1) CIS2=SIGN(1.D0,SOG2(4))
      GO TO 80000
C
C---------------------------------------------------------
C    CAS OU LES DEUX CRITERES SONT CONFONDUS
C---------------------------------------------------------
C
30000 IF(ICRI.EQ.1) CALL CHREPE(CO11,SI11,SOGMA,SOG1)
      IF(ICRI.EQ.1) CALL CHREPE(CO11,SI11,DSOGMA,DSOG1)
      IF(ICRI.EQ.2) CALL CHREPE(CO22,SI22,DSOGMA,DSOG2)
      IF(ICRI.EQ.2) CALL CHREPE(CO22,SI22,SOGMA,SOG2)
      GAMMA1=0.D0
      GAMMA2=0.D0
      IF(ICRI.EQ.1) CALL ECUBI1(SOG1,DSOG1,GAMMA1,PHI1,PHI1,NCAS,TRA1,
     1 XNU,YOUNG,COHE1,TSUG,DEFPL,DLA1,CO11,SI11,HACHE1,IDAM,KERRE)
      IF(ICRI.EQ.1) DLAM1=DLA1+DLAM1
      IF(ICRI.EQ.2) CALL ECUBI1(SOG2,DSOG2,GAMMA2,PHI2,PHI2,NCAS,TRA2,
     1 XNU,YOUNG,COHE2,TSUG,DEFPL,DLA2,CO22,SI22,HACHE2,IDAM,KERRE)
      IF(ICRI.EQ.2) DLAM2=DLA2+DLAM2
      DO 126 ITYP=1,4
      DSUG(ITYP)=SEG(ITYP)-TSUG(ITYP)
  126 CONTINUE
      IF(NCAS.EQ.1) CALL CPCONT(DSUG,DEFPL,XNU,YOUNG)
      IF(NCAS.NE.1) CALL DPCONT(DSUG,DEFPL,XNU,YOUNG)
      RETURN
C
C
C---------------------------------------------
C    CAS OU RAMENERAIT MAL A LA POINTE
C     SIGMA0 EN DEHORS DU CRITERE AVANT
C     DE RAMENER AU SOMMET
C---------------------------------------------
C20000 IF(ICHRI.EQ.2) GOTO 621
C     COEF1=PHI2*SI21*SI21-CIS2*SI21*CO21
C     IF(ABS(COEF1).LE.1.D-6) GOTO 622
C     COEF2=COHE2-(PHI2*CO21*CO21-CIS2*SI21*CO21)*TRA1
C     SOG1(1)=COEF2/COEF1
C     SOG1(2)=TRA1
C     SOG1(4)=0.D0
C     CALL CHREPE(CO11,-SI11,SOG1,SOGMA)
C     DO 623 ITYP=1,4
C     DSOGMA(ITYP)=SEG(ITYP)-SOGMA(ITYP)
C 623 CONTINUE
C     GOTO 622
C 621 COEF1=PHI1*SI21*SI21+CIS1*SI21*CO21
C     IF(ABS(COEF1).LE.1.D-6) GOTO 622
C     COEF2=COHE1-(PHI1*CO21*CO21+CIS1*SI21*CO21)*TRA2
C     SOG2(1)=COEF2/COEF1
C     SOG2(2)=TRA2
C     SOG2(4)=0.D0
C     CALL CHREPE(CO22,-SI22,SOG2,SOGMA)
C     DO 624 ITYP=1,4
C     DSOGMA(ITYP)=SEG(ITYP)-SOGMA(ITYP)
C 624 CONTINUE
C     IF(NCAS.NE.1) CALL DPCONT(DSOGMA,DEF,XNU,YOUNG)
C     DLA=SQRT(DEF(1)*DEF(1)+DEF(2)*DEF(2)+DEF(3)*DEF(3)
C    1 +DEF(4)*DEF(4))
C     DLAM1=DLA/2.D0+DLAM1
C     DLAM2=DLA/2.D0+DLAM2
C     DO 125 ITYP=1,4
C     DSUG(ITYP)=SEG(ITYP)-SOGMA(ITYP)
C     TSUG(ITYP)=SOGMA(ITYP)
C 125 CONTINUE
C     IF(NCAS.NE.1) CALL DPCONT(DSUG,DEFPL,XNU,YOUNG)
C     RETURN
C     END
C
C-------------------------------------------------------------
C    CAS OU ON EST A LA LIMITE EN TRACTION D'UN CRITERE
C       ET A LA SURFACE DE L'AUTRE  (ON RAMENE A LA POINTE)
C-------------------------------------------------------------
20000 IF(NCAS.EQ.1) CALL CPCONT(DSOGMA,DEF,XNU,YOUNG)
      IF(NCAS.NE.1) CALL DPCONT(DSOGMA,DEF,XNU,YOUNG)
      DLA=SQRT(DEF(1)*DEF(1)+DEF(2)*DEF(2)+DEF(3)*DEF(3)
     1 +DEF(4)*DEF(4))
      DLAM1=DLA/2.D0+DLAM1
      DLAM2=DLA/2.D0+DLAM2
      DO 125 ITYP=1,4
      DSUG(ITYP)=SEG(ITYP)-SOGMA(ITYP)
      TSUG(ITYP)=SOGMA(ITYP)
  125 CONTINUE
      IF(NCAS.EQ.1) CALL CPCONT(DSUG,DEFPL,XNU,YOUNG)
      IF(NCAS.NE.1) CALL DPCONT(DSUG,DEFPL,XNU,YOUNG)
      RETURN
      END