Download ecubic.eso

Back to the list

   1 : C ECUBIC    SOURCE    CHAT      05/01/12    23:28:10     5004           
   2 :       SUBROUTINE ECUBIC (SOGMA,DSOGMA,COHE1,COHE2,CO11,SI11,                    
   3 :      1 CO22,SI22,CO21,SI21,TRA1,TRA2,YOUNG,XNU,DEFPL,TSUG,IDAM,                 
   4 :      2 DLAM1,DLAM2,PSI1,PHI1,PSI2,PHI2,NCAS,ICRI,HACHE1,HACHE2,KERRE)           
   5 :       IMPLICIT INTEGER(I-N)
   6 :       IMPLICIT REAL*8(A-H,O-Z)                                                  
   7 :       DIMENSION DSOGMA(4),SOGMA(4),TSOG(4),TSOG1(4),TSOG2(4),                   
   8 :      1 SIP1(4),SAP1(4),SAP(4),SIP(4),SIP2(4),SAP2(4),TSUG(4),                   
   9 :      2 DSUG(4),SOG1(4),SOG2(4),DSOG1(4),DSOG2(4),A(3,3),B(3),                   
  10 :      3 C(3),D(3),DEF(4),DEFPL(4),SAG1(4),SAG2(4),SEG(4)                         
  11 :       DO 100 ITYP=1,4                                                           
  12 :       TSOG(ITYP)=DSOGMA(ITYP)+SOGMA(ITYP)                                       
  13 :       SEG(ITYP)=TSOG(ITYP)                                                      
  14 :   100 CONTINUE                                                                  
  15 :       ITER=1                                                                    
  16 :       ICHRI=0                                                                   
  17 :       ICHRO=0                                                                   
  18 :       ZER=0.D0                                                                  
  19 :       IDAM=12                                                                   
  20 :       ISOM1=0                                                                   
  21 :       ISOM2=0                                                                   
  22 :       DLA1=0.D0                                                                 
  23 :       DLA2=0.D0                                                                 
  24 : C------------------------------------------------------------------             
  25 : C ON EFFECTUE LE COUPLAGE                                                       
  26 : C------------------------------------------------------------------             
  27 : C                                                                               
  27 :       CALL CHREPE(CO11,SI11,TSOG,TSOG1)                                         
  27 :       CALL CHREPE(CO22,SI22,TSOG,TSOG2)                                         
  27 :       CALL CHREPE(CO11,SI11,SOGMA,SOG1)                                         
  27 :       CALL CHREPE(CO22,SI22,SOGMA,SOG2)                                         
  27 :       CALL CHREPE(CO11,SI11,DSOGMA,DSOG1)                                       
  27 :       CALL CHREPE(CO22,SI22,DSOGMA,DSOG2)                                       
  28 : C---------------------------------------------------------------                
  29 : C     SI ON EST A LA POINTE CAS PRATICULIER                                     
  30 : C     IL FAUT REGARDER SI ON SORT DANS LE CONE DES NORMALES                     
  31 : C---------------------------------------------------------------                
  32 : C                                                                               
  33 :       IF(SOG1(2).GT.(1.D0-1.D-6)*TRA1) ISOM1=1                                  
  34 :       IF(SOG2(2).GT.(1.D0-1.D-6)*TRA2) ISOM2=1                                  
  35 :       IF((ISOM1+ISOM2).EQ.0) GO TO 90000                                        
  36 : C                                                                               
  37 :       CIS1=SIGN(1.D0,SOG1(4))                                                   
  38 :       CIS2=SIGN(1.D0,SOG2(4))                                                   
  39 :       IF(SOG1(2).GT.(1.D0-1.D-6)*TRA1) CIS1=SIGN(1.D0,TSOG1(4))                 
  40 :       IF(SOG2(2).GT.(1.D0-1.D-6)*TRA2) CIS2=SIGN(1.D0,TSOG2(4))                 
  41 :       VCRA=DSOG1(2)-PSI1*ABS(DSOG1(4))                                          
  42 :       VCRB=DSOG2(2)-PSI2*ABS(DSOG2(4))                                          
  43 :       IF(ISOM1.EQ.1.AND.VCRA.GT.0.D0) ICHRO=1                                   
  44 :       IF(ISOM2.EQ.1.AND.VCRB.GT.0.D0) ICHRI=1                                   
  45 :       IF(ISOM1.EQ.1.AND.VCRA.GT.0.D0) GO TO 20000                               
  46 :       IF(ISOM2.EQ.1.AND.VCRB.GT.0.D0) GO TO 20000                               
  47 :       GOTO 80000                                                                
  48 : C                                                                               
  49 : 90000 CIS1=SIGN(1.D0,SOG1(4))                                                   
  50 :       CIS3=SIGN(1.D0,TSOG1(4))                                                  
  51 :       CIS2=SIGN(1.D0,SOG2(4))                                                   
  52 :       CIS4=SIGN(1.D0,TSOG2(4))                                                  
  52 : 80000 IF(NCAS.EQ.1) CALL CPHOOK(ZER,PSI1,ZER,CIS1,SIP1(1),SIP1(2)               
  53 :      1 ,SIP1(3),SIP1(4),XNU,YOUNG)                                              
  54 :       ICHRO=0                                                                   
  55 :       ICHRI=0                                                                   
  55 :       IF(NCAS.EQ.1) CALL CPHOOK(ZER,PSI2,ZER,CIS2,SAP2(1),SAP2(2)               
  56 :      1 ,SAP2(3),SAP2(4),XNU,YOUNG)                                              
  56 :       IF(NCAS.NE.1) CALL DPHOOK(ZER,PSI1,ZER,CIS1,SIP1(1),SIP1(2)               
  57 :      1 ,SIP1(3),SIP1(4),XNU,YOUNG)                                              
  57 :       IF(NCAS.NE.1) CALL DPHOOK(ZER,PSI2,ZER,CIS2,SAP2(1),SAP2(2)               
  58 :      1 ,SAP2(3),SAP2(4),XNU,YOUNG)                                              
  58 :       CALL CHREPE(CO21,SI21,SIP1,SIP2)                                          
  58 :       CALL CHREPE(CO21,-SI21,SAP2,SAP1)                                         
  58 :       CALL CHREPE(CO11,-SI11,SIP1,SIP)                                          
  58 :       CALL CHREPE(CO22,-SI22,SAP2,SAP)                                          
  59 :       ALPH1=PHI1*SIP1(2)+CIS1*SIP1(4)                                           
  60 :       ALPH2=PHI1*SAP1(2)+CIS1*SAP1(4)                                           
  61 :       ALPH3=-COHE1+PHI1*TSOG1(2)+CIS1*TSOG1(4)                                  
  62 :       ALPH4=PHI2*SIP2(2)+CIS2*SIP2(4)                                           
  63 :       ALPH5=PHI2*SAP2(2)+CIS2*SAP2(4)                                           
  64 :       ALPH6=-COHE2+PHI2*TSOG2(2)+CIS2*TSOG2(4)                                  
  65 :       DET=ALPH1*ALPH5-ALPH4*ALPH2                                               
  66 :       REF=YOUNG*YOUNG*1.D-7                                                     
  67 :       IF(DET.LE.REF) GO TO 30000                                                
  68 : C------------------------------------------------------------------             
  69 : C SI LE DETERMINANT EST NUL ON NE PEUT EFFECTUER LE COUPLAGE                    
  70 : C ON ECOULE DONC SUIVANT UN SEUL CRITERE  (20000)                               
  71 : C------------------------------------------------------------------             
  72 : C                                                                               
  73 :       DLA1=(ALPH3*ALPH5-ALPH6*ALPH2)/DET                                        
  74 :       DLA2=(ALPH1*ALPH6-ALPH3*ALPH4)/DET                                        
  75 :       DO 110 ITYP=1,4                                                           
  76 :       DSUG(ITYP)=DLA1*SIP(ITYP)+DLA2*SAP(ITYP)                                  
  77 :       TSUG(ITYP)=TSOG(ITYP)-DSUG(ITYP)                                          
  78 :       DSUG(ITYP)=SEG(ITYP)-TSUG(ITYP)                                           
  79 :   110 CONTINUE                                                                  
  79 :       CALL CHREPE(CO11,SI11,TSUG,SAG1)                                          
  79 :       CALL CHREPE(CO22,SI22,TSUG,SAG2)                                          
  80 : C-----------------------------------------------------                          
  81 : C   DANS LE CAS OU ON DEPASSE LA LIMITE EN TRACTION                             
  82 : C   POUR LE CRITERE 1 OU 2 ON CALCULE DES NOUVEAUX                              
  83 : C   DLAMDA                                                                      
  84 : C-----------------------------------------------------                          
  85 :       IF(SAG1(2).GT.TRA1*(1.D0+1.D-6).OR.SAG2(2).GT.TRA2*(1.D0+1.D-6))          
  86 :      1 GO TO 500                                                                
  87 :       DLAM1=DLA1+DLAM1                                                          
  88 :       DLAM2=DLA2+DLAM2                                                          
  88 :       IF(NCAS.EQ.1) CALL CPCONT(DSUG,DEFPL,XNU,YOUNG)                           
  88 :       IF(NCAS.NE.1) CALL DPCONT(DSUG,DEFPL,XNU,YOUNG)                           
  89 :       RETURN                                                                    
  90 :   500 C(3)=100.D0                                                               
  91 :       D(3)=100.D0                                                               
  92 :       IF(SAG1(2).LE.(TRA1*(1.D0+1.D-6))) GO TO 600                              
  93 :       B(1)=-COHE1+PHI1*SOG1(2)+CIS1*SOG1(4)                                     
  94 :       B(2)=-COHE2+PHI2*SOG2(2)+CIS2*SOG2(4)                                     
  95 :       B(3)=SOG1(2)-TRA1                                                         
  96 :       A(1,3)=-(PHI1*DSOG1(2)+CIS1*DSOG1(4))                                     
  97 :       A(2,3)=-(PHI2*DSOG2(2)+CIS2*DSOG2(4))                                     
  98 :       A(3,3)=-DSOG1(2)                                                          
  99 :       A(1,1)=ALPH1                                                              
 100 :       A(1,2)=ALPH2                                                              
 101 :       A(2,1)=ALPH4                                                              
 102 :       A(2,2)=ALPH5                                                              
 103 :       A(3,1)=SIP1(2)                                                            
 104 :       A(3,2)=SAP1(2)                                                            
 104 :       CALL RESNEQ(A,B,C,3,KERRE)                                                
 105 :       IF(SAG2(2).LE.(TRA2*(1.D0-1.D-6))) GO TO 800                              
 106 :  600  B(1)=COHE1-PHI1*SOG1(2)-CIS1*SOG1(4)                                      
 107 :       B(2)=COHE2-PHI2*SOG2(2)-CIS2*SOG2(4)                                      
 108 :       B(3)=SOG2(2)-TRA2                                                         
 109 :       A(1,3)=-(PHI1*DSOG1(2)+CIS1*DSOG1(4))                                     
 110 :       A(2,3)=-(PHI2*DSOG2(2)+CIS2*DSOG2(4))                                     
 111 :       A(3,3)=-DSOG2(2)                                                          
 112 :       A(1,1)=ALPH1                                                              
 113 :       A(1,2)=ALPH2                                                              
 114 :       A(2,1)=ALPH4                                                              
 115 :       A(2,2)=ALPH5                                                              
 116 :       A(3,1)=SIP2(2)                                                            
 117 :       A(3,2)=SAP2(2)                                                            
 117 :       CALL RESNEQ(A,B,D,3,KERRE)                                                
 118 :  800  IF(D(3).LE.C(3)) DLA1=D(1)                                                
 119 :       IF(D(3).GT.C(3)) ICHRO=1                                                  
 120 :       IF(D(3).LE.C(3)) ICHRI=1                                                  
 121 :       IF(ABS(D(3)-C(3)).LE.(ABS(C(3))*1.D-6)) ICHRI=1                           
 122 :       IF(ABS(D(3)-C(3)).LE.(ABS(C(3))*1.D-6)) ICHRO=1                           
 123 :       IF(D(3).LE.C(3)) DLA2=D(2)                                                
 124 :       IF(D(3).LE.C(3)) XAT=D(3)                                                 
 125 :       IF(D(3).GT.C(3)) DLA1=C(1)                                                
 126 :       IF(D(3).GT.C(3)) DLA2=C(2)                                                
 127 :       IF(D(3).GT.C(3)) XAT=C(3)                                                 
 128 :       DO 111 ITYP=1,4                                                           
 129 :       DSUG(ITYP)=DLA1*SIP(ITYP)+DLA2*SAP(ITYP)                                  
 130 :       SOGMA(ITYP)=SOGMA(ITYP)-DSUG(ITYP)+XAT*DSOGMA(ITYP)                       
 131 :       DSOGMA(ITYP)=(1.D0-XAT)*DSOGMA(ITYP)                                      
 132 :       TSOG(ITYP)=DSOGMA(ITYP)+SOGMA(ITYP)                                       
 133 :   111 CONTINUE                                                                  
 134 :       DLAM1=DLAM1+DLA1                                                          
 135 :       DLAM2=DLAM2+DLA2                                                          
 135 :       CALL CHREPE(CO11,SI11,TSOG,TSOG1)                                         
 135 :       CALL CHREPE(CO22,SI22,TSOG,TSOG2)                                         
 135 :       CALL CHREPE(CO11,SI11,SOGMA,SOG1)                                         
 136 :       VCRIT1=VCRITE(SOG1(2),SOG1(4),PHI1,COHE1)                                 
 136 :       CALL CHREPE(CO22,SI22,SOGMA,SOG2)                                         
 137 :       VCRIT2=VCRITE(SOG2(2),SOG2(4),PHI2,COHE2)                                 
 137 :       CALL CHREPE(CO22,SI22,DSOGMA,DSOG2)                                       
 137 :       CALL CHREPE(CO11,SI11,DSOGMA,DSOG1)                                       
 138 :       VCRA=DSOG1(2)-PSI1*ABS(DSOG1(4))                                          
 139 :       VCRB=DSOG2(2)-PSI2*ABS(DSOG2(4))                                          
 140 :       ITER=ITER+1                                                               
 141 :       IF(ICHRO.EQ.1.AND.VCRA.GT.0.D0) GO TO 20000                               
 142 :       IF(ICHRI.EQ.1.AND.VCRB.GT.0.D0) GO TO 20000                               
 143 :       IF(ITER.EQ.7) GO TO 20000                                                 
 144 : C                                                                               
 145 :       IF(ICHRO.EQ.1) CIS1=SIGN(1.D0,TSOG1(4))                                   
 146 :       IF(ICHRO.NE.1) CIS1=SIGN(1.D0,SOG1(4))                                    
 147 :       IF(ICHRI.EQ.1) CIS2=SIGN(1.D0,TSOG2(4))                                   
 148 :       IF(ICHRI.NE.1) CIS2=SIGN(1.D0,SOG2(4))                                    
 149 :       GO TO 80000                                                               
 150 : C                                                                               
 151 : C---------------------------------------------------------                      
 152 : C    CAS OU LES DEUX CRITERES SONT CONFONDUS                                    
 153 : C---------------------------------------------------------                      
 154 : C                                                                               
 154 : 30000 IF(ICRI.EQ.1) CALL CHREPE(CO11,SI11,SOGMA,SOG1)                           
 154 :       IF(ICRI.EQ.1) CALL CHREPE(CO11,SI11,DSOGMA,DSOG1)                         
 154 :       IF(ICRI.EQ.2) CALL CHREPE(CO22,SI22,DSOGMA,DSOG2)                         
 154 :       IF(ICRI.EQ.2) CALL CHREPE(CO22,SI22,SOGMA,SOG2)                           
 155 :       GAMMA1=0.D0                                                               
 156 :       GAMMA2=0.D0                                                               
 156 :       IF(ICRI.EQ.1) CALL ECUBI1(SOG1,DSOG1,GAMMA1,PHI1,PHI1,NCAS,TRA1,          
 157 :      1 XNU,YOUNG,COHE1,TSUG,DEFPL,DLA1,CO11,SI11,HACHE1,IDAM,KERRE)             
 158 :       IF(ICRI.EQ.1) DLAM1=DLA1+DLAM1                                            
 158 :       IF(ICRI.EQ.2) CALL ECUBI1(SOG2,DSOG2,GAMMA2,PHI2,PHI2,NCAS,TRA2,          
 159 :      1 XNU,YOUNG,COHE2,TSUG,DEFPL,DLA2,CO22,SI22,HACHE2,IDAM,KERRE)             
 160 :       IF(ICRI.EQ.2) DLAM2=DLA2+DLAM2                                            
 161 :       DO 126 ITYP=1,4                                                           
 162 :       DSUG(ITYP)=SEG(ITYP)-TSUG(ITYP)                                           
 163 :   126 CONTINUE                                                                  
 163 :       IF(NCAS.EQ.1) CALL CPCONT(DSUG,DEFPL,XNU,YOUNG)                           
 163 :       IF(NCAS.NE.1) CALL DPCONT(DSUG,DEFPL,XNU,YOUNG)                           
 164 :       RETURN                                                                    
 165 : C                                                                               
 166 : C                                                                               
 167 : C---------------------------------------------                                  
 168 : C    CAS OU RAMENERAIT MAL A LA POINTE                                          
 169 : C     SIGMA0 EN DEHORS DU CRITERE AVANT                                         
 170 : C     DE RAMENER AU SOMMET                                                      
 171 : C---------------------------------------------                                  
 172 : C20000 IF(ICHRI.EQ.2) GOTO 621                                                  
 173 : C     COEF1=PHI2*SI21*SI21-CIS2*SI21*CO21                                       
 174 : C     IF(ABS(COEF1).LE.1.D-6) GOTO 622                                          
 175 : C     COEF2=COHE2-(PHI2*CO21*CO21-CIS2*SI21*CO21)*TRA1                          
 176 : C     SOG1(1)=COEF2/COEF1                                                       
 177 : C     SOG1(2)=TRA1                                                              
 178 : C     SOG1(4)=0.D0                                                              
 178 : C     CALL CHREPE(CO11,-SI11,SOG1,SOGMA)                                        
 179 : C     DO 623 ITYP=1,4                                                           
 180 : C     DSOGMA(ITYP)=SEG(ITYP)-SOGMA(ITYP)                                        
 181 : C 623 CONTINUE                                                                  
 182 : C     GOTO 622                                                                  
 183 : C 621 COEF1=PHI1*SI21*SI21+CIS1*SI21*CO21                                       
 184 : C     IF(ABS(COEF1).LE.1.D-6) GOTO 622                                          
 185 : C     COEF2=COHE1-(PHI1*CO21*CO21+CIS1*SI21*CO21)*TRA2                          
 186 : C     SOG2(1)=COEF2/COEF1                                                       
 187 : C     SOG2(2)=TRA2                                                              
 188 : C     SOG2(4)=0.D0                                                              
 188 : C     CALL CHREPE(CO22,-SI22,SOG2,SOGMA)                                        
 189 : C     DO 624 ITYP=1,4                                                           
 190 : C     DSOGMA(ITYP)=SEG(ITYP)-SOGMA(ITYP)                                        
 191 : C 624 CONTINUE                                                                  
 191 : C     IF(NCAS.NE.1) CALL DPCONT(DSOGMA,DEF,XNU,YOUNG)                           
 192 : C     DLA=SQRT(DEF(1)*DEF(1)+DEF(2)*DEF(2)+DEF(3)*DEF(3)                        
 193 : C    1 +DEF(4)*DEF(4))                                                          
 194 : C     DLAM1=DLA/2.D0+DLAM1                                                      
 195 : C     DLAM2=DLA/2.D0+DLAM2                                                      
 196 : C     DO 125 ITYP=1,4                                                           
 197 : C     DSUG(ITYP)=SEG(ITYP)-SOGMA(ITYP)                                          
 198 : C     TSUG(ITYP)=SOGMA(ITYP)                                                    
 199 : C 125 CONTINUE                                                                  
 199 : C     IF(NCAS.NE.1) CALL DPCONT(DSUG,DEFPL,XNU,YOUNG)                           
 200 : C     RETURN                                                                    
 201 : C     END                                                                       
 202 : C                                                                               
 203 : C-------------------------------------------------------------                  
 204 : C    CAS OU ON EST A LA LIMITE EN TRACTION D'UN CRITERE                         
 205 : C       ET A LA SURFACE DE L'AUTRE  (ON RAMENE A LA POINTE)                     
 206 : C-------------------------------------------------------------                  
 206 : 20000 IF(NCAS.EQ.1) CALL CPCONT(DSOGMA,DEF,XNU,YOUNG)                           
 206 :       IF(NCAS.NE.1) CALL DPCONT(DSOGMA,DEF,XNU,YOUNG)                           
 207 :       DLA=SQRT(DEF(1)*DEF(1)+DEF(2)*DEF(2)+DEF(3)*DEF(3)                        
 208 :      1 +DEF(4)*DEF(4))                                                          
 209 :       DLAM1=DLA/2.D0+DLAM1                                                      
 210 :       DLAM2=DLA/2.D0+DLAM2                                                      
 211 :       DO 125 ITYP=1,4                                                           
 212 :       DSUG(ITYP)=SEG(ITYP)-SOGMA(ITYP)                                          
 213 :       TSUG(ITYP)=SOGMA(ITYP)                                                    
 214 :   125 CONTINUE                                                                  
 214 :       IF(NCAS.EQ.1) CALL CPCONT(DSUG,DEFPL,XNU,YOUNG)                           
 214 :       IF(NCAS.NE.1) CALL DPCONT(DSUG,DEFPL,XNU,YOUNG)                           
 215 :       RETURN                                                                    
 216 :       END                                                                       
 217 :  

© Cast3M 2003 - All rights reserved.
Disclaimer