Sources Esope | Cast3M

epsln2
C EPSLN2    SOURCE    CB215821  26/08/24    21:16:29     12622          
 
      SUBROUTINE EPSLN2(F,EPS,N)
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
 
-INC PPARAM
-INC CCOPTIO
 
      DIMENSION F(*),EPS(*)
      DIMENSION C(3,3),S(3,3),D(3)
*
* Affichage du gradient de la transformation
*
      IF (IIMPI.EQ.199) THEN
        WRITE(IOIMP,7771) N
 7771   FORMAT(2X,'EPSLN2 - N=',I3/)
        N2=N*N
        WRITE(IOIMP,7772) (F(I),I=1,N2)
 7772   FORMAT(2X,'F  '/(3(1X,1PE12.5)))
      ENDIF
*
*   REMPLISSAGE DE C = Ftrans.F
*
      CALL ZERO(C,3,3)
      DO 7773 I=1,N
        DO 1 J=1,N
          r_z=0.D0
          DO K=0,N-1
            KN=K*N
            r_z=r_z+F(KN+I)*F(KN+J)
          ENDDO
          C(I,J)=r_z
 1    CONTINUE
 7773 CONTINUE
*
      IF (IIMPI.EQ.199) THEN
        WRITE(IOIMP,7702) ((C(I,J),J=1,N),I=1,N)
      ENDIF
*
*    PETITE VERIFICATION DE LA SYMETRIE DE C
*
      TOL=1.D-10
      DO 7774 I=2,N
        DO 11 J=1,I-1
          IF(ABS(C(I,J)-C(J,I)).GE.TOL) THEN
            CALL ERREUR(26)
            RETURN
          ENDIF
 11   CONTINUE
 7774 CONTINUE
 
C*AV      DO I=1,N
C*AV        C(I,I)=C(I,I) - 1.D0
C*AV      ENDDO
*
*  Calcul des valeurs propres de C et des vecteurs propres
*
      NLOC = N
C* En modes 2D PLAN ou 2D AXI, on travaille sur une matrice 2x2 mais
C* on verifie la nullite des termes en C(3,i), i = 1 a 2
      IF (IFOUR.LE.0) THEN
        IF (N.EQ.3) THEN
          NLOC = 2
          IF (ABS(C(3,1)+C(1,3)).GE.TOL) THEN
            CALL ERREUR(26)
            RETURN
          ELSE IF (ABS(C(3,2)+C(2,3)).GE.TOL) THEN
            CALL ERREUR(26)
            RETURN
          ENDIF
        ENDIF
      ENDIF
 
      CALL JACOB3(C,NLOC,D,S)
 
      IF (IIMPI.EQ.199) THEN
        WRITE(IOIMP,7701) (D(K),K=1,N)
      ENDIF
*
*  Calcul de ln(U) = 1/2 ln(C) (valeurs propres)
*
      DO 2 I=1,N
        D(I) = 0.5D0*LOG(D(I))
C*AV        D(I) = 0.5D0*LOG(1.D0+D(I))
2     CONTINUE
*
      DO 7775 I=1,N
        DO 3 J=1,N
          r_z=0.D0
          DO 31 K=1,N
            r_z = r_z + S(I,K)*D(K)*S(J,K)
 31       CONTINUE
          C(I,J)=r_z
 3    CONTINUE
 7775 CONTINUE
*
      IF(IIMPI.EQ.199) THEN
        WRITE(IOIMP,7701) (D(K),K=1,N)
 7701   FORMAT(2X,'  D     '/(6(1X,1PE12.5)))
        WRITE(IOIMP,7702) ((C(I,J),J=1,N),I=1,N)
 7702   FORMAT(2X,'  C     '/(3(1X,1PE12.5)))
      ENDIF
*
*    RANGEMENT DANS EPS
*
      IF (N.EQ.2) THEN
        EPS(1)=C(1,1)
        EPS(2)=C(2,2)
        EPS(3)=C(2,1)+C(1,2)
        EPS(4)=0.D0
        EPS(5)=0.D0
        EPS(6)=0.D0
*
      ELSE IF (N.EQ.3) THEN
*
        IF (IFOUR.EQ.1.OR.IFOUR.EQ.2) THEN
          EPS(1)=C(1,1)
          EPS(2)=C(2,2)
          EPS(3)=C(3,3)
          EPS(4)=C(2,1)+C(1,2)
          EPS(5)=C(3,1)+C(1,3)
          EPS(6)=C(3,2)+C(2,3)
        ELSE IF (IFOUR.EQ.0.OR.IFOUR.EQ.-2.OR.IFOUR.EQ.-3) THEN
          EPS(1)=C(1,1)
          EPS(2)=C(2,2)
          EPS(3)=C(3,3)
          EPS(4)=C(2,1)+C(1,2)
          EPS(5)=0.D0
          EPS(6)=0.D0
        ELSE IF (IFOUR.EQ.-1) THEN
          IF (ABS(C(3,3)).GE.TOL) THEN
            CALL ERREUR(26)
            RETURN
          ENDIF
          EPS(1)=C(1,1)
          EPS(2)=C(2,2)
          EPS(3)=0.D0
          EPS(4)=C(2,1)+C(1,2)
          EPS(5)=0.D0
          EPS(6)=0.D0
* Modes de calcul 1D
        ELSE IF (IFOUR.GE.3.AND.IFOUR.LE.15) THEN
          EPS(1)=C(1,1)
          IF (IFOUR.EQ.3) THEN
            IF (ABS(C(2,2)).GE.TOL.OR.ABS(C(3,3)).GE.TOL) THEN
              CALL ERREUR(26)
              RETURN
            ENDIF
            EPS(2)=0.D0
            EPS(3)=0.D0
          ELSE IF (IFOUR.EQ.5.OR.IFOUR.EQ.7) THEN
            IF (ABS(C(3,3)).GE.TOL) THEN
              CALL ERREUR(26)
              RETURN
            ENDIF
            EPS(2)=C(2,2)
            EPS(3)=0.D0
          ELSE IF (IFOUR.EQ.4.OR.IFOUR.EQ.9.OR.IFOUR.EQ.12) THEN
            IF (ABS(C(2,2)).GE.TOL) THEN
              CALL ERREUR(26)
              RETURN
            ENDIF
            EPS(2)=0.D0
            EPS(3)=C(3,3)
          ELSE
            EPS(2)=C(2,2)
            EPS(3)=C(3,3)
          ENDIF
          EPS(4)=0.D0
          EPS(5)=0.D0
          EPS(6)=0.D0
        ELSE
          CALL ERREUR(19)
        ENDIF
*
      ELSE
        CALL ERREUR(19)
      ENDIF
*
      IF (IIMPI.EQ.199) THEN
        WRITE(IOIMP,7730) (EPS(K),K=1,6)
 7730   FORMAT(2X,'EPSLN2- EPS    '/(3(1X,1PE12.5)))
      ENDIF
*
      RETURN
      END