C EPPLUM SOURCE GF238795 18/02/05 21:15:14 9726 SUBROUTINE EPPLUM(TENS,PPLUS,IG,VAL1,VP1,QPLUS,Q,VP, . P,S,TAMP) C====================================================================== C C SOUS PROGRAMME DE CALCUL C POUR C TENSEUR DE DEFORMATION C C VERSION 1.0 C ----------- C C C CALCUL DE : C C 1- Valeurs et vecteurs propres C 2- Tenseur Q C 3- Tenseur Q+ C 4- Tenseur d ordre 4 P+ C C====================================================================== C C CREATION : F.CORMERY C E.N.S.M.A - LMPM C DEC 1992 C C====================================================================== IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C********************************************************************** C DIMENSIONS ET DATA C********************************************************************** C N36 N72 N75 CC DIMENSION PPLUS(3,3,3,3),TENS(6),QPLUS(6,6),Q(6,6),VP(3) C N84 N90 CC * ,P(3,3),S(6), C N99 CC * VP1(3),VAL1(3,3),TAMP(3,3) C DIMENSION PPLUS(3,3,3,*),TENS(*),QPLUS(6,*),Q(6,*),VP(*) * ,P(3,*),S(*), * VP1(*),VAL1(3,*),TAMP(3,*) DATA ZERO/0.D0/,UN/1.D0/, * PRECIS/1.D-08/,DPRECS/1.D-08/ INTEGER IM C---------------------------------------------------------------------- AMAX1(X,Y,Z,U,V,W)= MAX(X,Y,Z,U,V,W) C------------------ IM=0 MT=10 C********************************************************************** C INITIALISATION C********************************************************************** DO 5 J=1,6 DO 6 K=1,6 * P(J,K)=ZERO Q(J,K)=ZERO QPLUS(J,K)=ZERO 6 CONTINUE 5 CONTINUE DO 55 J=1,3 DO 55 K=1,3 55 P(J,K)=ZERO C---------------------------------------------------------------------- TENS(4)=TENS(4)/2 TENS(5)=TENS(5)/2 TENS(6)=TENS(6)/2 C********************************************************************** C NORMALISATION DU TENSEUR A C********************************************************************** C C---------------------------------------------------------------------- C Trouver la valeur max de TENS(I) C---------------------------------------------------------------------- DO 3 I=1,6 S(I)=ABS(TENS(I)) 3 CONTINUE C-------------- TMAX=AMAX1(S(1),S(2),S(3),S(4),S(5),S(6)) IF(TMAX.EQ.0.D0)TMAX=UN C---------------------------------------------------------------------- C Normaliser a un la composante de TENS(I) la plus grande C---------------------------------------------------------------------- DO 4 I=1,6 TENS(I)=TENS(I)/TMAX IF(ABS(TENS(I)).LE.1E-15) TENS(I)=0.D0 4 CONTINUE C------------------------ cas axes principaux ------------------------- NN=0 DO 234 IV=4,6 IF(ABS(TENS(IV)).LE.1E-15) NN=NN+1 234 CONTINUE IF(NN.EQ.3)THEN VP(1)=TENS(1) VP(2)=TENS(2) VP(3)=TENS(3) DO 235 I=1,3 P(I,1)=0 P(I,2)=0 P(I,3)=0 P(I,I)=1 235 CONTINUE goto 98 ENDIF C*********************************************************************** C CALCUL DES VALEURS PROPRES C*********************************************************************** CALL VALPRP(TENS(1),TENS(2),TENS(3),TENS(6),TENS(4),TENS(5), * VP(1),VP(2),VP(3)) C*********************************************************************** C CALCUL DES VECTEURS PROPRES C*********************************************************************** IM=2 IF(ABS(VP(1)-VP(2)).LT.1E-08)THEN VP(1)=VP(3) VP(2)=VP(2) VP(3)=VP(2) IM=2 ENDIF IF(ABS(VP(1)-VP(3)).LT.1E-08)THEN VP(1)=VP(2) VP(2)=VP(3) IM=2 ENDIF IF(ABS(VP(2)-VP(3)).LT.1E-08)THEN IM=2 ENDIF C---------------------------------------------------------------------- IMM=0 C---------------------------------------------------------------------- DO 10 I=1,IM C------------------- SDET1=(TENS(2)-VP(I))*(TENS(3)-VP(I))-TENS(4)**2 SDET2=(TENS(1)-VP(I))*(TENS(3)-VP(I))-TENS(5)**2 SDET3=(TENS(1)-VP(I))*(TENS(2)-VP(I))-TENS(6)**2 SDET4=TENS(6)*(TENS(3)-VP(I))-TENS(4)*TENS(5) SDET5=-TENS(5)*(TENS(2)-VP(I))+TENS(6)*TENS(4) SDET6=TENS(4)*(TENS(1)-VP(I))-TENS(6)*TENS(5) C---------------------------------------------------------------------- C WRITE(10,*)'MINEURS :' C WRITE(10,*)SDET1,SDET2,SDET3,SDET4,SDET5,SDET6 C---------------------------------------------------------------------- IF (ABS(SDET1).GT.DPRECS) THEN P(I,1)=UN P(I,2)=((-TENS(6)*(TENS(3)-VP(I)))+TENS(4)*TENS(5))/SDET1 P(I,3)=((-TENS(5)*(TENS(2)-VP(I)))+TENS(4)*TENS(6))/SDET1 GOTO 96 C------------------- ENDIF IF (ABS(SDET2).GT.DPRECS) THEN P(I,2)=UN P(I,1)=((-TENS(6)*(TENS(3)-VP(I)))+TENS(4)*TENS(5))/SDET2 P(I,3)=((-TENS(4)*(TENS(1)-VP(I)))+TENS(5)*TENS(6))/SDET2 GOTO 96 C------------------- ENDIF IF (ABS(SDET3).GT.DPRECS) THEN P(I,3)=UN P(I,1)=((-TENS(5)*(TENS(2)-VP(I)))+TENS(4)*TENS(6))/SDET3 P(I,2)=((-TENS(4)*(TENS(1)-VP(I)))+TENS(5)*TENS(6))/SDET3 GOTO 96 C-------------------- ENDIF IF (ABS(SDET4).GT.DPRECS) THEN P(I,1)=UN P(I,2)=((-(TENS(3)-vp(i))*(TENS(1)-VP(I)))+TENS(5)**2)/SDET4 P(I,3)=((TENS(4)*(TENS(1)-VP(I)))-TENS(5)*TENS(6))/SDET4 GOTO 96 C-------------------- ENDIF IF (ABS(SDET5).GT.DPRECS) THEN P(I,1)=(-(tens(4)**2)+(tens(2)-vp(i)))/sdet5 P(I,2)=(-tens(6)*(tens(3)-vp(i))+tens(5)*tens(4))/sdet5 P(I,3)=1 GOTO 96 C-------------------- ENDIF IF (ABS(SDET6).GT.DPRECS) THEN P(I,3)=UN P(I,1)=((-TENS(5)*TENS(4))+(TENS(3)-vp(i))*TENS(6))/SDET6 P(I,2)=((-(TENS(3)-vp(i))*(TENS(1)-VP(I)))+TENS(5)**2)/SDET6 ENDIF C----------------------------------------------------------------------- SSDET1=TENS(1)-VP(I) SSDET2=TENS(2)-VP(I) SSDET3=TENS(3)-VP(I) C-------------------CAS PARTICULIERS------------------------------------ IF (ABS(SSDET1).LE.PRECIS) THEN P(I,1)=1 P(I,2)=0 P(I,3)=0 GOTO 96 ENDIF IF (ABS(SSDET2).LE.PRECIS) THEN P(I,1)=0 P(I,2)=1 P(I,3)=0 GOTO 96 ENDIF IF (ABS(SSDET3).LE.PRECIS) THEN P(I,1)=0 P(I,2)=0 P(I,3)=1 GOTO 96 ENDIF IF (ABS(SSDET1).GT.PRECIS) THEN P(I,1)=-(TENS(6)+TENS(5))/SSDET1 P(I,2)=1 P(I,3)=1 C------------------- GOTO 96 ENDIF IF (ABS(SSDET2).GT.PRECIS) THEN P(I,1)=1 P(I,2)=-(TENS(6)+TENS(4))/SSDET2 P(I,3)=1 GOTO 96 C------------------- ENDIF IF (ABS(SSDET3).GT.PRECIS) THEN P(I,1)=1 P(I,2)=1 P(I,3)=-(TENS(5)+TENS(4))/SSDET3 GOTO 96 ENDIF C----------------------------------------------------------------------- WRITE(MT,*)'ERREUR DANS VPROP.FOR' WRITE(MT,1010) 1010 FORMAT(1X,'TENSEUR A SYM. D ORDRE 2 :', * /1X,'--------------------------'/) WRITE(MT,1001)TENS(1),TENS(6),TENS(5) 1001 FORMAT(15X,'* ',3e20.7,' *') WRITE(MT,1002)TENS(6),TENS(2),TENS(4) 1002 FORMAT(15X,'* ',3e20.7,' *') WRITE(MT,1003)TENS(5),TENS(4),TENS(3) 1003 FORMAT(15X,'* ',3e20.7,' *'/) STOP C----------------------------------------------------------------------- 96 CONTINUE C----------------------------------------------------------------------- 10 CONTINUE 98 CONTINUE IF(IM.EQ.2)THEN P(3,1)=P(1,2)*P(2,3)-P(1,3)*P(2,2) P(3,2)=P(1,3)*P(2,1)-P(1,1)*P(2,3) P(3,3)=P(1,1)*P(2,2)-P(1,2)*P(2,1) ENDIF C------------------------------------------------------------------------ C On verifie que la base formee est bien directe C------------------------------------------------------------------------ DIR=P(3,1)*(P(1,2)*P(2,3)-P(1,3)*P(2,2))+ * P(3,2)*(P(1,3)*P(2,1)-P(1,1)*P(2,3))+ * P(3,3)*(P(1,1)*P(2,2)-P(1,2)*P(2,1)) IF(DIR.LT.ZERO)THEN DO 12 J=1,3 TAMP1=VP(2) VP(2)=VP(1) VP(1)=TAMP1 TAMP(1,J)=P(1,J) P(1,J)=P(2,J) P(2,J)=TAMP(1,J) 12 CONTINUE ENDIF C------------------------------------------------------------------------ C Normalisation des vecteurs propres C------------------------------------------------------------------------ DO 11 J=1,5 DO 11 I=1,3 IF(ABS(P(I,1)).LT.1.E-15)P(I,1)=0.D0 IF(ABS(P(I,2)).LT.1.E-15)P(I,2)=0.D0 IF(ABS(P(I,3)).LT.1.E-15)P(I,3)=0.D0 RAC=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2) P(I,1)=P(I,1)/RAC P(I,2)=P(I,2)/RAC P(I,3)=P(I,3)/RAC 11 CONTINUE C------------------------------------------------------------------------ C Factorisation par le facteur de normalisation C------------------------------------------------------------------------ C VP1(I)=VP(I)*TMAX C************************************************************************ C Calcul de Q et Q+ C************************************************************************ DO 223 I=1,3 DO 50 J=1,3 DO 40 K=1,3 Q(J,K)=Q(J,K)+P(I,J)*P(I,K) IF (VP(I).GT.PRECIS) THEN QPLUS(J,K)=QPLUS(J,K)+P(I,J)*P(I,K) ENDIF 40 CONTINUE 50 CONTINUE C------------------------------------------------------------------------- C Fin de la boucle sur les valeur propres C------------------------------------------------------------------------- 223 CONTINUE C----------------------------------------------------------------------- C Verification si la base orthogonale est bien construite C----------------------------------------------------------------------- DO 51 J=1,3 DO 51 K=1,3 IF (ABS(Q(J,K)).LE.1E-15)Q(J,K)=ZERO IF (ABS(QPLUS(J,K)).LE.1E-15)QPLUS(J,K)=ZERO C----------------------- IF (J.NE.K)THEN IF(ABS(Q(J,K)).GT.1E-3)THEN WRITE(MT,*)'ERREUR DANS VPROP.FOR:' WRITE(MT,*)'- LA BASE N EST PAS ORTHOGONALE.' C---------------------- DO 118 I=1,6 TENS(I)=TENS(I) 118 CONTINUE C----------------------- WRITE(MT,1110) 1110 FORMAT(/1X,'TENSEUR A SYM. D ORDRE 2 :', * /1X,'--------------------------'/) WRITE(MT,1001)TENS(1),TENS(6),TENS(5) 1101 FORMAT(15X,'* ',3e20.7,' *') WRITE(MT,1002)TENS(6),TENS(2),TENS(4) 1102 FORMAT(15X,'* ',3e20.7,' *') WRITE(MT,1003)TENS(5),TENS(4),TENS(3) 1103 FORMAT(15X,'* ',3e20.7,' *') DO 31 I=1,3 SDET1=(TENS(2)-VP(I))*(TENS(3)-VP(I))-TENS(4)**2 SDET2=(TENS(1)-VP(I))*(TENS(3)-VP(I))-TENS(5)**2 SDET3=(TENS(1)-VP(I))*(TENS(2)-VP(I))-TENS(6)**2 SDET4=TENS(6)*(TENS(3)-VP(I))-TENS(4)*TENS(5) SDET5=-TENS(5)*(TENS(2)-VP(I))+TENS(6)*TENS(4) SDET6=TENS(4)*(TENS(1)-VP(I))-TENS(6)*TENS(5) write(mt,*)sdet1,sdet2,sdet3,sdet4,sdet5,sdet6 WRITE(MT,5000)VP(I),(P(I,L),L=1,3) 5000 FORMAT(3X,'VALEUR PRO. :',E12.5,' VECTEUR PROP.:',3E12.5) 31 CONTINUE C---------------- WRITE(MT,*)'********** EPPLUSM********************' WRITE(MT,*)' TENSEURS Q :' WRITE(MT,*) C---------------- DO 131 I=1,3 WRITE(MT,7000)(Q(I,L),L=1,3) 7000 FORMAT(5X,'* ',3E12.5,' *') 131 CONTINUE C---------------- WRITE(MT,*) WRITE(MT,*)' TENSEURS Q+ :' WRITE(MT,*) C---------------- DO 132 I=1,3 WRITE(MT,7001)(QPLUS(I,L),L=1,3) 7001 FORMAT(5X,'* ',3E12.5,' *') 132 CONTINUE C---------------- STOP C---------------- ENDIF ENDIF C----------------------- IF (J.EQ.K)THEN IF(ABS(Q(J,K)-UN).GT.1E-3)THEN WRITE(MT,*)'ERREUR DANS VPROP.FOR:' WRITE(MT,*)'- LA BASE N EST PAS ORTHOGONALE.' ENDIF ENDIF 51 CONTINUE 56 CONTINUE C********************************************************************** C CALCUL DE L OPERATEUR P+ C********************************************************************** DO 100 I=1,3 DO 100 J=1,3 DO 100 K=1,3 DO 100 L=1,3 PPLUS(I,J,K,L)=ZERO do 114 M=1,3 do 114 N=1,3 PPLUS(I,J,K,L)=PPLUS(I,J,K,L)+QPLUS(I,M) * *QPLUS(J,N)*(Q(K,M)*Q(L,N)+Q(L,M)*Q(K,N))/4.D0+ * QPLUS(J,M)*QPLUS(I,N)*(Q(K,M)*Q(L,N)+Q(L,M)*Q(K,N)) * /4.D0 114 continue IF(ABS(PPLUS(I,J,K,L)).LE.1E-15)PPLUS(I,J,L,K)=0.D0 100 CONTINUE C----------------------- c if (ig.eq.1)then c WRITE(MT,1117) c1117 FORMAT(1X,'TENSEUR A SYM. D ORDRE 2 :') c WRITE(MT,1001)TENS(1),TENS(6),TENS(5) c WRITE(MT,1002)TENS(6),TENS(2),TENS(4) c WRITE(MT,1003)TENS(5),TENS(4),TENS(3) c DO 234 I=1,3 c WRITE(MT,5000)VP(I),(P(I,L),L=1,3) c234 CONTINUE C---------------- c WRITE(MT,*)' TENSEURS Q :' C---------------- c DO 231 I=1,3 c WRITE(MT,7000)(Q(I,L),L=1,3) c231 CONTINUE c endif C----------------------------------------------------------------------- C Multiplication par le facteur de normalisation C----------------------------------------------------------------------- DO 110 I=1,6 TENS(I)=TENS(I)*TMAX 110 CONTINUE C----------------------------------------------------------------------- TENS(4)=2*TENS(4) TENS(5)=2*TENS(5) TENS(6)=2*TENS(6) C---------------------------------------------------------------------- DO 341 I=1,3 DO 341 J=1,3 VAL1(I,J)=P(I,J) 341 CONTINUE DO 342 I=1,3 VP1(I)=VP(I)*tmax 342 CONTINUE C------------------------------------------------------------------------- RETURN END C-----------------------------------------------------------------------