Esope source files

   1 : C EPSLN2    SOURCE    FANDEUR   12/03/15    21:24:18     7312           
   2 : 
   3 :       SUBROUTINE EPSLN2(F,EPS,N)
   4 : 
   5 :       IMPLICIT INTEGER(I-N)
   6 :       IMPLICIT REAL*8(A-H,O-Z)
   7 : 
   8 : -INC CCOPTIO
   9 : 
  10 :       DIMENSION F(*),EPS(*)
  11 :       DIMENSION C(3,3),S(3,3),D(3)
  12 : *
  13 : * Affichage du gradient de la transformation
  14 : *
  15 :       IF (IIMPI.EQ.199) THEN
  16 :         WRITE(IOIMP,7771) N
  17 :  7771   FORMAT(2X,'EPSLN2 - N=',I3/)
  18 :         N2=N*N
  19 :         WRITE(IOIMP,7772) (F(I),I=1,N2)
  20 :  7772   FORMAT(2X,'F  '/(3(1X,1PE12.5)))
  21 :       ENDIF
  22 : *
  23 : *   REMPLISSAGE DE C = Ftrans.F
  24 : *
  24 :       CALL ZERO(C,3,3)
  25 :       DO 1 I=1,N
  26 :         DO 1 J=1,N
  27 :           r_z=0.D0
  28 :           DO K=0,N-1
  29 :             KN=K*N
  30 :             r_z=r_z+F(KN+I)*F(KN+J)
  31 :           ENDDO
  32 :           C(I,J)=r_z
  33 :  1    CONTINUE
  34 : *
  35 :       IF (IIMPI.EQ.199) THEN
  36 :         WRITE(IOIMP,7702) ((C(I,J),J=1,N),I=1,N)
  37 :       ENDIF
  38 : *
  39 : *    PETITE VERIFICATION DE LA SYMETRIE DE C
  40 : *
  41 :       TOL=1.D-10
  42 :       DO 11 I=2,N
  43 :         DO 11 J=1,I-1
  44 :           IF(ABS(C(I,J)-C(J,I)).GE.TOL) THEN
  44 :             CALL ERREUR(26)
  45 :             RETURN
  46 :           ENDIF
  47 :  11   CONTINUE
  48 : 
  49 : C*AV      DO I=1,N
  50 : C*AV        C(I,I)=C(I,I) - 1.D0
  51 : C*AV      ENDDO
  52 : *
  53 : *  Calcul des valeurs propres de C et des vecteurs propres
  54 : *
  55 :       NLOC = N
  56 : C* En modes 2D PLAN ou 2D AXI, on travaille sur une matrice 2x2 mais
  57 : C* on verifie la nullite des termes en C(3,i), i = 1 a 2
  58 :       IF (IFOUR.LE.0) THEN
  59 :         IF (N.EQ.3) THEN
  60 :           NLOC = 2
  61 :           IF (ABS(C(3,1)+C(1,3)).GE.TOL) THEN
  61 :             CALL ERREUR(26)
  62 :             RETURN
  63 :           ELSE IF (ABS(C(3,2)+C(2,3)).GE.TOL) THEN
  63 :             CALL ERREUR(26)
  64 :             RETURN
  65 :           ENDIF
  66 :         ENDIF
  67 :       ENDIF
  68 : 
  68 :       CALL JACOB3(C,NLOC,D,S)
  69 : 
  70 :       IF (IIMPI.EQ.199) THEN
  71 :         WRITE(IOIMP,7701) (D(K),K=1,N)
  72 :       ENDIF
  73 : *
  74 : *  Calcul de ln(U) = 1/2 ln(C) (valeurs propres)
  75 : *
  76 :       DO 2 I=1,N
  77 :         D(I) = 0.5D0*LOG(D(I))
  78 : C*AV        D(I) = 0.5D0*LOG(1.D0+D(I))
  79 : 2     CONTINUE
  80 : *
  81 :       DO 3 I=1,N
  82 :         DO 3 J=1,N
  83 :           r_z=0.D0
  84 :           DO 31 K=1,N
  85 :             r_z = r_z + S(I,K)*D(K)*S(J,K)
  86 :  31       CONTINUE
  87 :           C(I,J)=r_z
  88 :  3    CONTINUE
  89 : *
  90 :       IF(IIMPI.EQ.199) THEN
  91 :         WRITE(IOIMP,7701) (D(K),K=1,N)
  92 :  7701   FORMAT(2X,'  D     '/(6(1X,1PE12.5)))
  93 :         WRITE(IOIMP,7702) ((C(I,J),J=1,N),I=1,N)
  94 :  7702   FORMAT(2X,'  C     '/(3(1X,1PE12.5)))
  95 :       ENDIF
  96 : *
  97 : *    RANGEMENT DANS EPS
  98 : *
  99 :       IF (N.EQ.2) THEN
 100 :         EPS(1)=C(1,1)
 101 :         EPS(2)=C(2,2)
 102 :         EPS(3)=C(2,1)+C(1,2)
 103 :         EPS(4)=0.D0
 104 :         EPS(5)=0.D0
 105 :         EPS(6)=0.D0
 106 : *
 107 :       ELSE IF (N.EQ.3) THEN
 108 : *
 109 :         IF (IFOUR.EQ.1.OR.IFOUR.EQ.2) THEN
 110 :           EPS(1)=C(1,1)
 111 :           EPS(2)=C(2,2)
 112 :           EPS(3)=C(3,3)
 113 :           EPS(4)=C(2,1)+C(1,2)
 114 :           EPS(5)=C(3,1)+C(1,3)
 115 :           EPS(6)=C(3,2)+C(2,3)
 116 :         ELSE IF (IFOUR.EQ.0.OR.IFOUR.EQ.-2.OR.IFOUR.EQ.-3) THEN
 117 :           EPS(1)=C(1,1)
 118 :           EPS(2)=C(2,2)
 119 :           EPS(3)=C(3,3)
 120 :           EPS(4)=C(2,1)+C(1,2)
 121 :           EPS(5)=0.D0
 122 :           EPS(6)=0.D0
 123 :         ELSE IF (IFOUR.EQ.-1) THEN
 124 :           IF (ABS(C(3,3)).GE.TOL) THEN
 124 :             CALL ERREUR(26)
 125 :             RETURN
 126 :           ENDIF
 127 :           EPS(1)=C(1,1)
 128 :           EPS(2)=C(2,2)
 129 :           EPS(3)=0.D0
 130 :           EPS(4)=C(2,1)+C(1,2)
 131 :           EPS(5)=0.D0
 132 :           EPS(6)=0.D0
 133 : * Modes de calcul 1D
 134 :         ELSE IF (IFOUR.GE.3.AND.IFOUR.LE.15) THEN
 135 :           EPS(1)=C(1,1)
 136 :           IF (IFOUR.EQ.3) THEN
 137 :             IF (ABS(C(2,2)).GE.TOL.OR.ABS(C(3,3)).GE.TOL) THEN
 137 :               CALL ERREUR(26)
 138 :               RETURN
 139 :             ENDIF
 140 :             EPS(2)=0.D0
 141 :             EPS(3)=0.D0
 142 :           ELSE IF (IFOUR.EQ.5.OR.IFOUR.EQ.7) THEN
 143 :             IF (ABS(C(3,3)).GE.TOL) THEN
 143 :               CALL ERREUR(26)
 144 :               RETURN
 145 :             ENDIF
 146 :             EPS(2)=C(2,2)
 147 :             EPS(3)=0.D0
 148 :           ELSE IF (IFOUR.EQ.4.OR.IFOUR.EQ.9.OR.IFOUR.EQ.12) THEN
 149 :             IF (ABS(C(2,2)).GE.TOL) THEN
 149 :               CALL ERREUR(26)
 150 :               RETURN
 151 :             ENDIF
 152 :             EPS(2)=0.D0
 153 :             EPS(3)=C(3,3)
 154 :           ELSE
 155 :             EPS(2)=C(2,2)
 156 :             EPS(3)=C(3,3)
 157 :           ENDIF
 158 :           EPS(4)=0.D0
 159 :           EPS(5)=0.D0
 160 :           EPS(6)=0.D0
 161 :         ELSE
 161 :           CALL ERREUR(19)
 162 :         ENDIF
 163 : *
 164 :       ELSE
 164 :         CALL ERREUR(19)
 165 :       ENDIF
 166 : *
 167 :       IF (IIMPI.EQ.199) THEN
 168 :         WRITE(IOIMP,7730) (EPS(K),K=1,6)
 169 :  7730   FORMAT(2X,'EPSLN2- EPS    '/(3(1X,1PE12.5)))
 170 :       ENDIF
 171 : *
 172 :       RETURN
 173 :       END
 174 : 
 175 :