1 : C EPPLUM SOURCE CHAT 05/01/12 23:41:13 5004 2 : SUBROUTINE EPPLUM(TENS,PPLUS,IG,VAL1,VP1,QPLUS,Q,VP, 3 : . P,S,TAMP) 4 : C====================================================================== 5 : C 6 : C SOUS PROGRAMME DE CALCUL 7 : C POUR 8 : C TENSEUR DE DEFORMATION 9 : C 10 : C VERSION 1.0 11 : C ----------- 12 : C 13 : C 14 : C CALCUL DE : 15 : C 16 : C 1- Valeurs et vecteurs propres 17 : C 2- Tenseur Q 18 : C 3- Tenseur Q+ 19 : C 4- Tenseur d ordre 4 P+ 20 : C 21 : C====================================================================== 22 : C 23 : C CREATION : F.CORMERY 24 : C E.N.S.M.A - LMPM 25 : C DEC 1992 26 : C 27 : C====================================================================== 28 : IMPLICIT INTEGER(I-N) 29 : IMPLICIT REAL*8 (A-H,O-Z) 30 : C********************************************************************** 31 : C DIMENSIONS ET DATA 32 : C********************************************************************** 33 : C N36 N72 N75 34 : CC DIMENSION PPLUS(3,3,3,3),TENS(6),QPLUS(6,6),Q(6,6),VP(3) 35 : C N84 N90 36 : CC * ,P(3,3),S(6), 37 : C N99 38 : CC * VP1(3),VAL1(3,3),TAMP(3,3) 39 : C 40 : DIMENSION PPLUS(3,3,3,*),TENS(*),QPLUS(6,*),Q(6,*),VP(*) 41 : * ,P(3,*),S(*), 42 : * VP1(*),VAL1(3,*),TAMP(3,*) 43 : DATA ZERO/0.D0/,UN/1.D0/, 44 : * PRECIS/1.D-08/,DPRECS/1.D-08/ 45 : C---------------------------------------------------------------------- 46 : AMAX1(X,Y,Z,U,V,W)=DMAX1(X,Y,Z,U,V,W) 47 : C------------------ 48 : MT=10 49 : C********************************************************************** 50 : C INITIALISATION 51 : C********************************************************************** 52 : DO 5 J=1,6 53 : DO 6 K=1,6 54 : * P(J,K)=ZERO 55 : Q(J,K)=ZERO 56 : QPLUS(J,K)=ZERO 57 : 6 CONTINUE 58 : 5 CONTINUE 59 : DO 55 J=1,3 60 : DO 55 K=1,3 61 : 55 P(J,K)=ZERO 62 : C---------------------------------------------------------------------- 63 : TENS(4)=TENS(4)/2 64 : TENS(5)=TENS(5)/2 65 : TENS(6)=TENS(6)/2 66 : C********************************************************************** 67 : C NORMALISATION DU TENSEUR A 68 : C********************************************************************** 69 : C 70 : C---------------------------------------------------------------------- 71 : C Trouver la valeur max de TENS(I) 72 : C---------------------------------------------------------------------- 73 : DO 3 I=1,6 74 : S(I)=ABS(TENS(I)) 75 : 3 CONTINUE 76 : C-------------- 77 : TMAX=AMAX1(S(1),S(2),S(3),S(4),S(5),S(6)) 78 : IF(TMAX.EQ.0.D0)TMAX=UN 79 : C---------------------------------------------------------------------- 80 : C Normaliser a un la composante de TENS(I) la plus grande 81 : C---------------------------------------------------------------------- 82 : DO 4 I=1,6 83 : TENS(I)=TENS(I)/TMAX 84 : IF(ABS(TENS(I)).LE.1E-15) TENS(I)=0.D0 85 : 4 CONTINUE 86 : C------------------------ cas axes principaux ------------------------- 87 : NN=0 88 : DO 234 IV=4,6 89 : IF(ABS(TENS(IV)).LE.1E-15) NN=NN+1 90 : 234 CONTINUE 91 : IF(NN.EQ.3)THEN 92 : VP(1)=TENS(1) 93 : VP(2)=TENS(2) 94 : VP(3)=TENS(3) 95 : DO 235 I=1,3 96 : P(I,1)=0 97 : P(I,2)=0 98 : P(I,3)=0 99 : P(I,I)=1 100 : 235 CONTINUE 101 : goto 98 102 : ENDIF 103 : C*********************************************************************** 104 : C CALCUL DES VALEURS PROPRES 105 : C*********************************************************************** 105 : CALL VALPRP(TENS(1),TENS(2),TENS(3),TENS(6),TENS(4),TENS(5), 106 : * VP(1),VP(2),VP(3)) 107 : C*********************************************************************** 108 : C CALCUL DES VECTEURS PROPRES 109 : C*********************************************************************** 110 : IM=2 111 : IF(ABS(VP(1)-VP(2)).LT.1E-08)THEN 112 : VP(1)=VP(3) 113 : VP(2)=VP(2) 114 : VP(3)=VP(2) 115 : IM=2 116 : ENDIF 117 : IF(ABS(VP(1)-VP(3)).LT.1E-08)THEN 118 : VP(1)=VP(2) 119 : VP(2)=VP(3) 120 : IM=2 121 : ENDIF 122 : IF(ABS(VP(2)-VP(3)).LT.1E-08)THEN 123 : IM=2 124 : ENDIF 125 : C---------------------------------------------------------------------- 126 : IMM=0 127 : C---------------------------------------------------------------------- 128 : DO 10 I=1,IM 129 : C------------------- 130 : SDET1=(TENS(2)-VP(I))*(TENS(3)-VP(I))-TENS(4)**2 131 : SDET2=(TENS(1)-VP(I))*(TENS(3)-VP(I))-TENS(5)**2 132 : SDET3=(TENS(1)-VP(I))*(TENS(2)-VP(I))-TENS(6)**2 133 : SDET4=TENS(6)*(TENS(3)-VP(I))-TENS(4)*TENS(5) 134 : SDET5=-TENS(5)*(TENS(2)-VP(I))+TENS(6)*TENS(4) 135 : SDET6=TENS(4)*(TENS(1)-VP(I))-TENS(6)*TENS(5) 136 : C---------------------------------------------------------------------- 137 : C WRITE(10,*)'MINEURS :' 138 : C WRITE(10,*)SDET1,SDET2,SDET3,SDET4,SDET5,SDET6 139 : C---------------------------------------------------------------------- 140 : 141 : IF (ABS(SDET1).GT.DPRECS) THEN 142 : P(I,1)=UN 143 : P(I,2)=((-TENS(6)*(TENS(3)-VP(I)))+TENS(4)*TENS(5))/SDET1 144 : P(I,3)=((-TENS(5)*(TENS(2)-VP(I)))+TENS(4)*TENS(6))/SDET1 145 : GOTO 96 146 : C------------------- 147 : ENDIF 148 : IF (ABS(SDET2).GT.DPRECS) THEN 149 : P(I,2)=UN 150 : P(I,1)=((-TENS(6)*(TENS(3)-VP(I)))+TENS(4)*TENS(5))/SDET2 151 : P(I,3)=((-TENS(4)*(TENS(1)-VP(I)))+TENS(5)*TENS(6))/SDET2 152 : GOTO 96 153 : C------------------- 154 : ENDIF 155 : IF (ABS(SDET3).GT.DPRECS) THEN 156 : P(I,3)=UN 157 : P(I,1)=((-TENS(5)*(TENS(2)-VP(I)))+TENS(4)*TENS(6))/SDET3 158 : P(I,2)=((-TENS(4)*(TENS(1)-VP(I)))+TENS(5)*TENS(6))/SDET3 159 : GOTO 96 160 : C-------------------- 161 : ENDIF 162 : IF (ABS(SDET4).GT.DPRECS) THEN 163 : P(I,1)=UN 164 : P(I,2)=((-(TENS(3)-vp(i))*(TENS(1)-VP(I)))+TENS(5)**2)/SDET4 165 : P(I,3)=((TENS(4)*(TENS(1)-VP(I)))-TENS(5)*TENS(6))/SDET4 166 : GOTO 96 167 : C-------------------- 168 : ENDIF 169 : IF (ABS(SDET5).GT.DPRECS) THEN 170 : P(I,1)=(-(tens(4)**2)+(tens(2)-vp(i)))/sdet5 171 : P(I,2)=(-tens(6)*(tens(3)-vp(i))+tens(5)*tens(4))/sdet5 172 : P(I,3)=1 173 : GOTO 96 174 : C-------------------- 175 : ENDIF 176 : IF (ABS(SDET6).GT.DPRECS) THEN 177 : P(I,3)=UN 178 : P(I,1)=((-TENS(5)*TENS(4))+(TENS(3)-vp(i))*TENS(6))/SDET6 179 : P(I,2)=((-(TENS(3)-vp(i))*(TENS(1)-VP(I)))+TENS(5)**2)/SDET6 180 : ENDIF 181 : C----------------------------------------------------------------------- 182 : SSDET1=TENS(1)-VP(I) 183 : SSDET2=TENS(2)-VP(I) 184 : SSDET3=TENS(3)-VP(I) 185 : C-------------------CAS PARTICULIERS------------------------------------ 186 : IF (ABS(SSDET1).LE.PRECIS) THEN 187 : P(I,1)=1 188 : P(I,2)=0 189 : P(I,3)=0 190 : GOTO 96 191 : ENDIF 192 : IF (ABS(SSDET2).LE.PRECIS) THEN 193 : P(I,1)=0 194 : P(I,2)=1 195 : P(I,3)=0 196 : GOTO 96 197 : ENDIF 198 : IF (ABS(SSDET3).LE.PRECIS) THEN 199 : P(I,1)=0 200 : P(I,2)=0 201 : P(I,3)=1 202 : GOTO 96 203 : ENDIF 204 : IF (ABS(SSDET1).GT.PRECIS) THEN 205 : P(I,1)=-(TENS(6)+TENS(5))/SSDET1 206 : P(I,2)=1 207 : P(I,3)=1 208 : C------------------- 209 : GOTO 96 210 : ENDIF 211 : IF (ABS(SSDET2).GT.PRECIS) THEN 212 : P(I,1)=1 213 : P(I,2)=-(TENS(6)+TENS(4))/SSDET2 214 : P(I,3)=1 215 : GOTO 96 216 : C------------------- 217 : ENDIF 218 : IF (ABS(SSDET3).GT.PRECIS) THEN 219 : P(I,1)=1 220 : P(I,2)=1 221 : P(I,3)=-(TENS(5)+TENS(4))/SSDET3 222 : GOTO 96 223 : ENDIF 224 : C----------------------------------------------------------------------- 225 : WRITE(MT,*)'ERREUR DANS VPROP.FOR' 226 : WRITE(MT,1010) 227 : 1010 FORMAT(1X,'TENSEUR A SYM. D ORDRE 2 :', 228 : * /1X,'--------------------------'/) 229 : WRITE(MT,1001)TENS(1),TENS(6),TENS(5) 230 : 1001 FORMAT(15X,'* ',3e20.7,' *') 231 : WRITE(MT,1002)TENS(6),TENS(2),TENS(4) 232 : 1002 FORMAT(15X,'* ',3e20.7,' *') 233 : WRITE(MT,1003)TENS(5),TENS(4),TENS(3) 234 : 1003 FORMAT(15X,'* ',3e20.7,' *'/) 235 : STOP 236 : C----------------------------------------------------------------------- 237 : 96 CONTINUE 238 : C----------------------------------------------------------------------- 239 : 10 CONTINUE 240 : 98 CONTINUE 241 : IF(IM.EQ.2)THEN 242 : P(3,1)=P(1,2)*P(2,3)-P(1,3)*P(2,2) 243 : P(3,2)=P(1,3)*P(2,1)-P(1,1)*P(2,3) 244 : P(3,3)=P(1,1)*P(2,2)-P(1,2)*P(2,1) 245 : ENDIF 246 : C------------------------------------------------------------------------ 247 : C On verifie que la base formee est bien directe 248 : C------------------------------------------------------------------------ 249 : DIR=P(3,1)*(P(1,2)*P(2,3)-P(1,3)*P(2,2))+ 250 : * P(3,2)*(P(1,3)*P(2,1)-P(1,1)*P(2,3))+ 251 : * P(3,3)*(P(1,1)*P(2,2)-P(1,2)*P(2,1)) 252 : IF(DIR.LT.ZERO)THEN 253 : DO 12 J=1,3 254 : TAMP1=VP(2) 255 : VP(2)=VP(1) 256 : VP(1)=TAMP1 257 : TAMP(1,J)=P(1,J) 258 : P(1,J)=P(2,J) 259 : P(2,J)=TAMP(1,J) 260 : 12 CONTINUE 261 : ENDIF 262 : C------------------------------------------------------------------------ 263 : C Normalisation des vecteurs propres 264 : C------------------------------------------------------------------------ 265 : DO 11 J=1,5 266 : DO 11 I=1,3 267 : IF(ABS(P(I,1)).LT.1.E-15)P(I,1)=0.D0 268 : IF(ABS(P(I,2)).LT.1.E-15)P(I,2)=0.D0 269 : IF(ABS(P(I,3)).LT.1.E-15)P(I,3)=0.D0 270 : RAC=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2) 271 : P(I,1)=P(I,1)/RAC 272 : P(I,2)=P(I,2)/RAC 273 : P(I,3)=P(I,3)/RAC 274 : 11 CONTINUE 275 : C------------------------------------------------------------------------ 276 : C Factorisation par le facteur de normalisation 277 : C------------------------------------------------------------------------ 278 : C VP1(I)=VP(I)*TMAX 279 : C************************************************************************ 280 : C Calcul de Q et Q+ 281 : C************************************************************************ 282 : DO 223 I=1,3 283 : DO 50 J=1,3 284 : DO 40 K=1,3 285 : Q(J,K)=Q(J,K)+P(I,J)*P(I,K) 286 : IF (VP(I).GT.PRECIS) THEN 287 : QPLUS(J,K)=QPLUS(J,K)+P(I,J)*P(I,K) 288 : ENDIF 289 : 40 CONTINUE 290 : 50 CONTINUE 291 : C------------------------------------------------------------------------- 292 : C Fin de la boucle sur les valeur propres 293 : C------------------------------------------------------------------------- 294 : 223 CONTINUE 295 : C----------------------------------------------------------------------- 296 : C Verification si la base orthogonale est bien construite 297 : C----------------------------------------------------------------------- 298 : DO 51 J=1,3 299 : DO 51 K=1,3 300 : IF (ABS(Q(J,K)).LE.1E-15)Q(J,K)=ZERO 301 : IF (ABS(QPLUS(J,K)).LE.1E-15)QPLUS(J,K)=ZERO 302 : C----------------------- 303 : IF (J.NE.K)THEN 304 : IF(ABS(Q(J,K)).GT.1E-3)THEN 305 : WRITE(MT,*)'ERREUR DANS VPROP.FOR:' 306 : WRITE(MT,*)'- LA BASE N EST PAS ORTHOGONALE.' 307 : C---------------------- 308 : DO 118 I=1,6 309 : TENS(I)=TENS(I) 310 : 118 CONTINUE 311 : C----------------------- 312 : WRITE(MT,1110) 313 : 1110 FORMAT(/1X,'TENSEUR A SYM. D ORDRE 2 :', 314 : * /1X,'--------------------------'/) 315 : WRITE(MT,1001)TENS(1),TENS(6),TENS(5) 316 : 1101 FORMAT(15X,'* ',3e20.7,' *') 317 : WRITE(MT,1002)TENS(6),TENS(2),TENS(4) 318 : 1102 FORMAT(15X,'* ',3e20.7,' *') 319 : WRITE(MT,1003)TENS(5),TENS(4),TENS(3) 320 : 1103 FORMAT(15X,'* ',3e20.7,' *') 321 : DO 31 I=1,3 322 : SDET1=(TENS(2)-VP(I))*(TENS(3)-VP(I))-TENS(4)**2 323 : SDET2=(TENS(1)-VP(I))*(TENS(3)-VP(I))-TENS(5)**2 324 : SDET3=(TENS(1)-VP(I))*(TENS(2)-VP(I))-TENS(6)**2 325 : SDET4=TENS(6)*(TENS(3)-VP(I))-TENS(4)*TENS(5) 326 : SDET5=-TENS(5)*(TENS(2)-VP(I))+TENS(6)*TENS(4) 327 : SDET6=TENS(4)*(TENS(1)-VP(I))-TENS(6)*TENS(5) 328 : write(mt,*)sdet1,sdet2,sdet3,sdet4,sdet5,sdet6 329 : WRITE(MT,5000)VP(I),(P(I,L),L=1,3) 330 : 5000 FORMAT(3X,'VALEUR PRO. :',E12.5,' VECTEUR PROP.:',3E12.5) 331 : 31 CONTINUE 332 : C---------------- 333 : WRITE(MT,*)'********** EPPLUSM********************' 334 : WRITE(MT,*)' TENSEURS Q :' 335 : WRITE(MT,*) 336 : C---------------- 337 : DO 131 I=1,3 338 : WRITE(MT,7000)(Q(I,L),L=1,3) 339 : 7000 FORMAT(5X,'* ',3E12.5,' *') 340 : 131 CONTINUE 341 : C---------------- 342 : WRITE(MT,*) 343 : WRITE(MT,*)' TENSEURS Q+ :' 344 : WRITE(MT,*) 345 : C---------------- 346 : DO 132 I=1,3 347 : WRITE(MT,7001)(QPLUS(I,L),L=1,3) 348 : 7001 FORMAT(5X,'* ',3E12.5,' *') 349 : 132 CONTINUE 350 : C---------------- 351 : STOP 352 : C---------------- 353 : ENDIF 354 : ENDIF 355 : C----------------------- 356 : IF (J.EQ.K)THEN 357 : IF(ABS(Q(J,K)-UN).GT.1E-3)THEN 358 : WRITE(MT,*)'ERREUR DANS VPROP.FOR:' 359 : WRITE(MT,*)'- LA BASE N EST PAS ORTHOGONALE.' 360 : ENDIF 361 : ENDIF 362 : 51 CONTINUE 363 : 56 CONTINUE 364 : C********************************************************************** 365 : C CALCUL DE L OPERATEUR P+ 366 : C********************************************************************** 367 : DO 100 I=1,3 368 : DO 100 J=1,3 369 : DO 100 K=1,3 370 : DO 100 L=1,3 371 : PPLUS(I,J,K,L)=ZERO 372 : do 114 M=1,3 373 : do 114 N=1,3 374 : PPLUS(I,J,K,L)=PPLUS(I,J,K,L)+QPLUS(I,M) 375 : * *QPLUS(J,N)*(Q(K,M)*Q(L,N)+Q(L,M)*Q(K,N))/4.D0+ 376 : * QPLUS(J,M)*QPLUS(I,N)*(Q(K,M)*Q(L,N)+Q(L,M)*Q(K,N)) 377 : * /4.D0 378 : 114 continue 379 : IF(ABS(PPLUS(I,J,K,L)).LE.1E-15)PPLUS(I,J,L,K)=0.D0 380 : 100 CONTINUE 381 : C----------------------- 382 : c if (ig.eq.1)then 383 : c WRITE(MT,1117) 384 : c1117 FORMAT(1X,'TENSEUR A SYM. D ORDRE 2 :') 385 : c WRITE(MT,1001)TENS(1),TENS(6),TENS(5) 386 : c WRITE(MT,1002)TENS(6),TENS(2),TENS(4) 387 : c WRITE(MT,1003)TENS(5),TENS(4),TENS(3) 388 : c DO 234 I=1,3 389 : c WRITE(MT,5000)VP(I),(P(I,L),L=1,3) 390 : c234 CONTINUE 391 : C---------------- 392 : c WRITE(MT,*)' TENSEURS Q :' 393 : C---------------- 394 : c DO 231 I=1,3 395 : c WRITE(MT,7000)(Q(I,L),L=1,3) 396 : c231 CONTINUE 397 : c endif 398 : C----------------------------------------------------------------------- 399 : C Multiplication par le facteur de normalisation 400 : C----------------------------------------------------------------------- 401 : DO 110 I=1,6 402 : TENS(I)=TENS(I)*TMAX 403 : 110 CONTINUE 404 : C----------------------------------------------------------------------- 405 : TENS(4)=2*TENS(4) 406 : TENS(5)=2*TENS(5) 407 : TENS(6)=2*TENS(6) 408 : C---------------------------------------------------------------------- 409 : DO 341 I=1,3 410 : DO 341 J=1,3 411 : VAL1(I,J)=P(I,J) 412 : 341 CONTINUE 413 : DO 342 I=1,3 414 : VP1(I)=VP(I)*tmax 415 : 342 CONTINUE 416 : C------------------------------------------------------------------------- 417 : RETURN 418 : END 419 : C----------------------------------------------------------------------- 420 : 421 :
© Cast3M 2003 - All rights reserved.
Disclaimer