1 : C VECT SOURCE CHAT 05/01/13 04:05:16 5004 2 : SUBROUTINE VECT(TENS,VAL,VP,P,S,TAMP) 3 : C====================================================================== 4 : C 5 : C SOUS PROGRAMME DE CALCUL 6 : C POUR 7 : C TENSEUR DE DEFORMATION 8 : C 9 : C VERSION 1.0 10 : C ----------- 11 : C 12 : C 13 : C CALCUL DE : 14 : C 15 : C 1- Valeurs et vecteurs propres 16 : C 2- Tenseur Q 17 : C 18 : C====================================================================== 19 : C 20 : C CREATION : F.CORMERY 21 : C E.N.S.M.A - LMPM 22 : C JUILLET 1993 23 : C 24 : C====================================================================== 25 : IMPLICIT INTEGER(I-N) 26 : IMPLICIT REAL*8 (A-H,O-Z) 27 : C********************************************************************** 28 : C DIMENSIONS ET DATA 29 : C********************************************************************** 30 : C 31 : CC DIMENSION TENS(6),VP(3) 32 : C N9 N15 33 : CC * ,P(3,3),S(6), 34 : C N24 35 : CC * VAL(3,3),TAMP(3,3) 36 : C 37 : DIMENSION TENS(*),VP(*) 38 : * ,P(3,*),S(*), 39 : * VAL(3,*),TAMP(3,*) 40 : DATA ZERO/0.D0/,UN/1.D0/, 41 : * PRECIS/1.D-08/,DPRECS/1.D-08/ 42 : C---------------------------------------------------------------------- 43 : AMAX1(X,Y,Z,U,V,W)=DMAX1(X,Y,Z,U,V,W) 44 : C------------------ 45 : MT=10 46 : C********************************************************************** 47 : C INITIALISATION 48 : C********************************************************************** 49 : * DO 5 J=1,6 50 : * DO 6 K=1,6 51 : * P(J,K)=ZERO 52 : * Q(J,K)=ZERO 53 : * QPLUS(J,K)=ZERO 54 : *6 CONTINUE 55 : 5 CONTINUE 56 : DO 55 J=1,3 57 : DO 55 K=1,3 58 : 55 P(J,K)=ZERO 59 : C---------------------------------------------------------------------- 60 : C********************************************************************** 61 : C NORMALISATION DU TENSEUR A 62 : C********************************************************************** 63 : C 64 : C---------------------------------------------------------------------- 65 : C Trouver la valeur max de TENS(I) 66 : C---------------------------------------------------------------------- 67 : DO 3 I=1,6 68 : S(I)=ABS(TENS(I)) 69 : 3 CONTINUE 70 : C-------------- 71 : TMAX=AMAX1(S(1),S(2),S(3),S(4),S(5),S(6)) 72 : IF(TMAX.EQ.0.D0)TMAX=UN 73 : C---------------------------------------------------------------------- 74 : C Normaliser a un la composante de TENS(I) la plus grande 75 : C---------------------------------------------------------------------- 76 : DO 4 I=1,6 77 : TENS(I)=TENS(I)/TMAX 78 : IF(ABS(TENS(I)).LE.1E-15) TENS(I)=0.D0 79 : 4 CONTINUE 80 : C------------------------ cas axes principaux ------------------------- 81 : NN=0 82 : DO 234 IV=4,6 83 : IF(ABS(TENS(IV)).LE.1E-15) NN=NN+1 84 : 234 CONTINUE 85 : IF(NN.EQ.3)THEN 86 : VP(1)=TENS(1) 87 : VP(2)=TENS(2) 88 : VP(3)=TENS(3) 89 : DO 235 I=1,3 90 : P(I,1)=0 91 : P(I,2)=0 92 : P(I,3)=0 93 : P(I,I)=1 94 : 235 CONTINUE 95 : goto 98 96 : ENDIF 97 : C*********************************************************************** 98 : C CALCUL DES VALEURS PROPRES 99 : C*********************************************************************** 99 : CALL VALPRP(TENS(1),TENS(2),TENS(3),TENS(6),TENS(4),TENS(5), 100 : * VP(1),VP(2),VP(3)) 101 : C*********************************************************************** 102 : C CALCUL DES VECTEURS PROPRES 103 : C*********************************************************************** 104 : IM=2 105 : IF(ABS(VP(1)-VP(2)).LT.1E-8)THEN 106 : VP(1)=VP(3) 107 : VP(2)=VP(2) 108 : VP(3)=VP(2) 109 : IM=2 110 : ENDIF 111 : IF(ABS(VP(1)-VP(3)).LT.1E-8)THEN 112 : VP(1)=VP(2) 113 : VP(2)=VP(3) 114 : IM=2 115 : ENDIF 116 : IF(ABS(VP(2)-VP(3)).LT.1E-8)THEN 117 : IM=2 118 : ENDIF 119 : C---------------------------------------------------------------------- 120 : IMM=0 121 : C---------------------------------------------------------------------- 122 : DO 10 I=1,IM 123 : C------------------- 124 : SDET1=(TENS(2)-VP(I))*(TENS(3)-VP(I))-TENS(4)**2 125 : SDET2=(TENS(1)-VP(I))*(TENS(3)-VP(I))-TENS(5)**2 126 : SDET3=(TENS(1)-VP(I))*(TENS(2)-VP(I))-TENS(6)**2 127 : SDET4=TENS(6)*(TENS(3)-VP(I))-TENS(4)*TENS(5) 128 : SDET5=-TENS(5)*(TENS(2)-VP(I))+TENS(6)*TENS(4) 129 : SDET6=TENS(4)*(TENS(1)-VP(I))-TENS(6)*TENS(5) 130 : C---------------------------------------------------------------------- 131 : C WRITE(10,*)'MINEURS :' 132 : C WRITE(10,*)SDET1,SDET2,SDET3,SDET4,SDET5,SDET6 133 : C---------------------------------------------------------------------- 134 : 135 : IF (ABS(SDET1).GT.DPRECS) THEN 136 : P(I,1)=UN 137 : P(I,2)=((-TENS(6)*(TENS(3)-VP(I)))+TENS(4)*TENS(5))/SDET1 138 : P(I,3)=((-TENS(5)*(TENS(2)-VP(I)))+TENS(4)*TENS(6))/SDET1 139 : GOTO 96 140 : C------------------- 141 : ENDIF 142 : IF (ABS(SDET2).GT.DPRECS) THEN 143 : P(I,2)=UN 144 : P(I,1)=((-TENS(6)*(TENS(3)-VP(I)))+TENS(4)*TENS(5))/SDET2 145 : P(I,3)=((-TENS(4)*(TENS(1)-VP(I)))+TENS(5)*TENS(6))/SDET2 146 : GOTO 96 147 : C------------------- 148 : ENDIF 149 : IF (ABS(SDET3).GT.DPRECS) THEN 150 : P(I,3)=UN 151 : P(I,1)=((-TENS(5)*(TENS(2)-VP(I)))+TENS(4)*TENS(6))/SDET3 152 : P(I,2)=((-TENS(4)*(TENS(1)-VP(I)))+TENS(5)*TENS(6))/SDET3 153 : GOTO 96 154 : C-------------------- 155 : ENDIF 156 : IF (ABS(SDET4).GT.DPRECS) THEN 157 : P(I,1)=UN 158 : P(I,2)=((-(TENS(3)-vp(i))*(TENS(1)-VP(I)))+TENS(5)**2)/SDET4 159 : P(I,3)=((TENS(4)*(TENS(1)-VP(I)))-TENS(5)*TENS(6))/SDET4 160 : GOTO 96 161 : C-------------------- 162 : ENDIF 163 : IF (ABS(SDET5).GT.DPRECS) THEN 164 : P(I,1)=(-(tens(4)**2)+(tens(2)-vp(i)))/sdet5 165 : P(I,2)=(-tens(6)*(tens(3)-vp(i))+tens(5)*tens(4))/sdet5 166 : P(I,3)=1 167 : GOTO 96 168 : C-------------------- 169 : ENDIF 170 : IF (ABS(SDET6).GT.DPRECS) THEN 171 : P(I,3)=UN 172 : P(I,1)=((-TENS(5)*TENS(4))+(TENS(3)-vp(i))*TENS(6))/SDET6 173 : P(I,2)=((-(TENS(3)-vp(i))*(TENS(1)-VP(I)))+TENS(5)**2)/SDET6 174 : ENDIF 175 : C----------------------------------------------------------------------- 176 : SSDET1=TENS(1)-VP(I) 177 : SSDET2=TENS(2)-VP(I) 178 : SSDET3=TENS(3)-VP(I) 179 : C-------------------CAS PARTICULIERS------------------------------------ 180 : IF (ABS(SSDET1).LE.PRECIS) THEN 181 : P(I,1)=1 182 : P(I,2)=0 183 : P(I,3)=0 184 : GOTO 96 185 : ENDIF 186 : IF (ABS(SSDET2).LE.PRECIS) THEN 187 : P(I,1)=0 188 : P(I,2)=1 189 : P(I,3)=0 190 : GOTO 96 191 : ENDIF 192 : IF (ABS(SSDET3).LE.PRECIS) THEN 193 : P(I,1)=0 194 : P(I,2)=0 195 : P(I,3)=1 196 : GOTO 96 197 : ENDIF 198 : IF (ABS(SSDET1).GT.PRECIS) THEN 199 : P(I,1)=-(TENS(6)+TENS(5))/SSDET1 200 : P(I,2)=1 201 : P(I,3)=1 202 : C------------------- 203 : GOTO 96 204 : ENDIF 205 : IF (ABS(SSDET2).GT.PRECIS) THEN 206 : P(I,1)=1 207 : P(I,2)=-(TENS(6)+TENS(4))/SSDET2 208 : P(I,3)=1 209 : GOTO 96 210 : C------------------- 211 : ENDIF 212 : IF (ABS(SSDET3).GT.PRECIS) THEN 213 : P(I,1)=1 214 : P(I,2)=1 215 : P(I,3)=-(TENS(5)+TENS(4))/SSDET3 216 : GOTO 96 217 : ENDIF 218 : C----------------------------------------------------------------------- 219 : WRITE(MT,*)'ERREUR DANS VPROP.FOR' 220 : WRITE(MT,1010) 221 : 1010 FORMAT(1X,'TENSEUR A SYM. D ORDRE 2 :', 222 : * /1X,'--------------------------'/) 223 : WRITE(MT,1001)TENS(1),TENS(6),TENS(5) 224 : 1001 FORMAT(15X,'* ',3e20.7,' *') 225 : WRITE(MT,1002)TENS(6),TENS(2),TENS(4) 226 : 1002 FORMAT(15X,'* ',3e20.7,' *') 227 : WRITE(MT,1003)TENS(5),TENS(4),TENS(3) 228 : 1003 FORMAT(15X,'* ',3e20.7,' *'/) 229 : STOP 230 : C----------------------------------------------------------------------- 231 : 96 CONTINUE 232 : C----------------------------------------------------------------------- 233 : 10 CONTINUE 234 : 98 CONTINUE 235 : IF(IM.EQ.2)THEN 236 : P(3,1)=P(1,2)*P(2,3)-P(1,3)*P(2,2) 237 : P(3,2)=P(1,3)*P(2,1)-P(1,1)*P(2,3) 238 : P(3,3)=P(1,1)*P(2,2)-P(1,2)*P(2,1) 239 : ENDIF 240 : C------------------------------------------------------------------------ 241 : C On verifie que la base formee est bien directe 242 : C------------------------------------------------------------------------ 243 : DIR=P(3,1)*(P(1,2)*P(2,3)-P(1,3)*P(2,2))+ 244 : * P(3,2)*(P(1,3)*P(2,1)-P(1,1)*P(2,3))+ 245 : * P(3,3)*(P(1,1)*P(2,2)-P(1,2)*P(2,1)) 246 : IF(DIR.LT.ZERO)THEN 247 : DO 12 J=1,3 248 : TAMP1=VP(2) 249 : VP(2)=VP(1) 250 : VP(1)=TAMP1 251 : TAMP(1,J)=P(1,J) 252 : P(1,J)=P(2,J) 253 : P(2,J)=TAMP(1,J) 254 : 12 CONTINUE 255 : ENDIF 256 : C------------------------------------------------------------------------ 257 : C Normalisation des vecteurs propres 258 : C------------------------------------------------------------------------ 259 : DO 11 J=1,5 260 : DO 11 I=1,3 261 : IF(ABS(P(I,1)).LT.1.E-15)P(I,1)=0.D0 262 : IF(ABS(P(I,2)).LT.1.E-15)P(I,2)=0.D0 263 : IF(ABS(P(I,3)).LT.1.E-15)P(I,3)=0.D0 264 : RAC=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2) 265 : P(I,1)=P(I,1)/RAC 266 : P(I,2)=P(I,2)/RAC 267 : P(I,3)=P(I,3)/RAC 268 : 11 CONTINUE 269 : C----------------------------------------------------------------------- 270 : C Multiplication par le facteur de normalisation 271 : C----------------------------------------------------------------------- 272 : DO 110 I=1,6 273 : TENS(I)=TENS(I)*TMAX 274 : 110 CONTINUE 275 : C----------------------------------------------------------------------- 276 : DO 341 I=1,3 277 : DO 341 J=1,3 278 : VAL(I,J)=P(I,J) 279 : 341 CONTINUE 280 : DO 342 I=1,3 281 : VP(I)=VP(I)*tmax 282 : 342 CONTINUE 283 : C------------------------------------------------------------------------- 284 : RETURN 285 : END 286 : C----------------------------------------------------------------------- 287 : 288 :
© Cast3M 2003 - All rights reserved.
Disclaimer