1 : C DERIV1 SOURCE AM 09/12/15 21:15:15 6587 2 : SUBROUTINE DERIV1(TAU,SIG,EPSV,VAR,SIGD,EPSVD,VARD,EPSTHD, 3 : &DSPT,XMAT,NSTRS,NVARI,NCOMAT,DD,DDV,DDINV,DDINVp,YOUNG,NYOUNG, 4 : &XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC,TI,TPOINT,ITHHER,VARF, 5 : & IB,IGAU,KERRE) 6 : 7 : C----------------------------------------------------------------------- 8 : C Objet: Cette subroutine calcule les derivees des variables internes 9 : C pour un materiau viscoplastique a endommagement et ecrouissage 10 : C isotrope en regime anisotherme 11 : C----------------------------------------------------------------------- 12 : C 13 : C----------------------------------------------------------------------- 14 : C Entree: TAU pas de temps 15 : C TI temperature au debut du pas 16 : C TPOINT derivee de la temperature 17 : C SIG(NSTRS) tenseur des contraintes dans le repère général 18 : C EPSV(NSTRS) tenseur des deformations dans le repère général 19 : C VAR(NVARI) tableau contenant les variables internes 20 : C pilotant les equations 21 : C VAR(1)=p ; VAR(2)=D11 ; VAR(3)=D22; VAR(4)=D33; 22 : C VAR(5)=D12; VAR(6)=D13; VAR(7)=D23; VAR(8)=DMAX 23 : C EPSTHD vitesse de dilatation thermique au debut du pas 24 : C DSPT(NSTRS,NSTRS) vitesse de deformation totale 25 : C XMAT(NCOMAT) tableau des parametres scalaires du materiau 26 : C a une temperature TI donnee 27 : C XMAT(1)=YOUNG ;XMAT(2)=XNU ;XMAT(3)=N 28 : C XMAT(4)=M ;XMAT(5)=KK; XMAT(6)=ALF2; 29 : C XMAT(7)=BET2 ;XMAT(8)=r; XMAT(9)=A ; 30 : C XMAT(10)=q ;XMAT(11)=ALPHAT 31 : C XMAT(12)=RHO; XMAT(13)=SIGY 32 : C MFR indice de la formulation mecanique(seulement massif 33 : C ou coque pour les materiaux endommageables 34 : C ICARA nombre de caracteristiques geometriques des elements 35 : C finis 36 : C XCAR(ICARA) tableau des caracteristiques geometriques des 37 : C elements finis 38 : C IFOURB= -2 EN CONTR. PLANES 39 : C -1 EN DEFORM. PLANES 40 : C 0 EN AXISYMETRIE 41 : C 1 EN SERIE DE FOURIER 42 : C 2 EN TRIDIM 43 : C INDIC=0, 1 OU -1 pour plasticite avec endommagement 44 : C =2 OU -2 pour viscoplasticite avec endommagement 45 : C ITHHER = 0 pas de chargement thermique et materiau constant 46 : C = 1 chargement thermique et materiau constant 47 : C = 2 chargement thermique et materiau(T) 48 : C------------------------------------------------------------------------ 49 : C 50 : C------------------------------------------------------------------------ 51 : C Sortie: EPSVD(NSTRS) derivee du tenseur des deformations 52 : C VARD(NVARI) tableau contenant les derivees des variables 53 : C internes 54 : C VARF(NVARI) état des variables internes à la fin du sous programme 55 : C SIGD(NSTRS) derivee des contraintes 56 : C DD(NSTRS,NSTRS) matrice de Hooke au debut du pas 57 : C DDV(NSTRS,NSTRS) derivee de DD 58 : C DDINV(NSTRS,NSTRS) inverse de DD 59 : C------------------------------------------------------------------------ 60 : C 61 : C 62 : IMPLICIT INTEGER(I-N) 63 : IMPLICIT REAL*8(A-H,O-Z) 64 : DIMENSION SIG(*),EPSV(*),VAR(*),XCAR(*),YOUNG(*),XNU(*) 65 : DIMENSION SIGD(*),EPSVD(*),VARD(*),DDINV(NSTRS,NSTRS) 66 : DIMENSION XMAT(*),DSPT(*),EPSTHD(*),VARF(NVARI) 67 : DIMENSION DDV(NSTRS,NSTRS),DDINVp(NSTRS,NSTRS),DD(NSTRS,NSTRS) 68 : DIMENSION EPSED(6),XX(6),Di(3),VaP(3),AA(3),BB(6) 69 : DIMENSION CO(3,3),COM(3,3),DCOT(3,3),COE(3,3),DCOE(3,3) 70 : DIMENSION AAA(3,3),BBB(3,3),CCC(3,3),UN2(3,3) 71 : DIMENSION FTH(3,3),ROT(3,3),ROTT(3,3),VcP(3,3),S(3,3) 72 : DIMENSION D1(3,3),D2(3,3),DV1(3,3),DI1(3,3),DI2(3,3) 73 : PARAMETER (UN=1.D0, UNM=-1.D0) 74 : LOGICAL endonul 75 : C 76 : C 77 : endonul = .FALSE. 78 : IF ((VAR(2).EQ.0).AND.(VAR(3).EQ.0).AND.(VAR(4).EQ.0).AND. 79 : & (VAR(5).EQ.0).AND.(VAR(6).EQ.0).AND.(VAR(7).EQ.0).AND. 80 : & (VAR(8).EQ.0)) THEN 81 : endonul = .TRUE. 82 : ENDIF 83 : 84 : C 85 : C----------------------------------------------------------------- 86 : C 87 : C EXPRESSIONS DES TENSEURS UTILISES DANS LE MODELE 88 : C 89 : C----------------------------------------------------------------- 90 : C 91 : C---------------------------------------------- 92 : C TENSEUR DES CONTRAINTES (dans repere general) 93 : C---------------------------------------------- 94 : C 94 : CALL ZERO(CO,3,3) 95 : IF (mfr.eq.1) THEN 96 : DO I=1,3 97 : CO(I,I) = SIG(I) 98 : ENDDO 99 : CO(1,2) = SIG(4) 100 : IF(NSTRS.NE.6) GOTO 50 101 : CO(1,3) = SIG(5) 102 : CO(2,3) = SIG(6) 103 : 50 CONTINUE 104 : CO(2,1) = CO(1,2) 105 : CO(3,1) = CO(1,3) 106 : CO(3,2) = CO(2,3) 107 : ENDIF 108 : C 109 : C------------------------------------ 110 : C TENSEUR D ENDOMMAGEMENT ET 111 : C TENSEUR D ENDOMMAGEMENT DIAGONALISE 112 : C------------------------------------ 113 : C 114 : XD11=VAR(2) 115 : XD22=VAR(3) 116 : XD33=VAR(4) 117 : XD12=VAR(5) 118 : XD13=VAR(6) 119 : XD23=VAR(7) 119 : CALL ZERO(D1,3,3) 120 : D1(1,1) = XD11 121 : D1(2,2) = XD22 122 : D1(3,3) = XD33 123 : D1(1,2) = XD12 124 : D1(1,3) = XD13 125 : D1(2,3) = XD23 126 : D1(2,1) = D1(1,2) 127 : D1(3,1) = D1(1,3) 128 : D1(3,2) = D1(2,3) 129 : 130 : * 131 : * AM 14/12/09 TEST SELON NSTRS 132 : * 133 : IF(IFOURB.LE.0) THEN 134 : JDIM = 2 134 : CALL JACOB4(D1,JDIM,Di,ROT) 135 : ELSE 136 : JDIM = 3 136 : CALL JACOB4(D1,JDIM,Di,ROT) 137 : ENDIF 138 : * 138 : CALL ZERO(D2,3,3) 139 : D2(1,1) = Di(1) 140 : D2(2,2) = Di(2) 141 : D2(3,3) = Di(3) 142 : w1 = 1.D0 - Di(1) 143 : w2 = 1.D0 - Di(2) 144 : w3 = 1.D0 - Di(3) 145 : IF (D2(1,1).GE.1.D0) w1=1.D-6 146 : IF (D2(2,2).GE.1.D0) w2=1.D-6 147 : IF (D2(3,3).GE.1.D0) w3=1.D-6 148 : 149 : ADMAX = DMAX1(Di(1),Di(2),Di(3)) 150 : C 151 : C------------------------ 152 : C TENSEUR UNITE D ORDRE 2 153 : C------------------------ 154 : C 154 : CALL ZERO(UN2,3,3) 155 : DO I=1,3 156 : UN2(I,I) = 1.D0 157 : ENDDO 158 : C 159 : C------------------------------------------------------- 160 : C TENSEUR DES CONTRAINTES MOTRICES (dans repere general) 161 : C------------------------------------------------------- 162 : C 163 : EQ=XMAT(10) 164 : DO I=1,3 165 : DO J=1,3 166 : ROTT(I,J)=ROT(I,J) 167 : ENDDO 168 : ENDDO 168 : CALL INVALM (ROTT, 3 , 3 , KER, 1.D-12) 169 : C 170 : IF (endonul) THEN 171 : DO I=1,3 172 : DO J=1,3 173 : COM(I,J) = CO(I,J) 174 : ENDDO 175 : ENDDO 176 : ELSE 177 : C --- Matrice de [(I-D)/det(I-D)]**(q/2) dans le repère principal de D --- 177 : CALL ZERO(DI1,3,3) 178 : DO I=1,3 179 : XNB=(1.D0-Di(I))/(w1*w2*w3) 180 : IF (XNB.GE.0.D0) THEN 181 : DI1(I,I)=(XNB)**(EQ/2.D0) 182 : ELSE 183 : DI1(I,I)=0.D0 184 : WRITE(6,*) 'erreur1, IB = ', IB, ' Igau=',igau 185 : KERRE = 99 186 : GO TO 100 187 : c Do j = 1,7 188 : c write(6,*) 'var(',j,')=', var(j) 189 : c ENDDO 190 : ENDIF 191 : ENDDO 192 : C --- passage dans le repère de général --- 192 : CALL MULMAT(AAA,ROT,DI1,3,3,3) 192 : CALL MULMAT(DI2,AAA,ROTT,3,3,3) 193 : 193 : CALL MULMAT(AAA,DI2,CO,3,3,3) 193 : CALL MULMAT(COM,AAA,DI2,3,3,3) 194 : ENDIF 195 : C 196 : C----------------------------------- 197 : C TENSEUR DES CONTRAINTES EFFECTIVES 198 : C----------------------------------- 199 : C 200 : IF (endonul) THEN 201 : DO I=1,3 202 : DO J=1,3 203 : COE(I,J) = CO(I,J) 204 : ENDDO 205 : ENDDO 206 : ELSE 207 : C --- Matrice de [I-D]**(-1/2) dans le repère principal de D --- 207 : CALL ZERO(DI1,3,3) 207 : CALL ZERO(DI2,3,3) 208 : DO I=1,3 209 : XNB=1.D0-Di(I) 210 : IF (XNB.GT.0.D0) THEN 211 : DI1(I,I)=SQRT(1.D0/(XNB)) 212 : ELSE 213 : DI1(I,I)=0.D0 214 : WRITE(6,*) 'erreur2 IB=', IB, ' IGAU= ', IGAU 215 : KERRE=99 216 : GO TO 100 217 : ENDIF 218 : ENDDO 219 : C --- passage dans le repère de général --- 219 : CALL MULMAT(AAA,ROT,DI1,3,3,3) 219 : CALL MULMAT(DI2,AAA,ROTT,3,3,3) 220 : c DO i=1,3 221 : c DO j=1,3 222 : c WRITE(6,*) 'DI2(',i,',',j,')=',DI2(i,j) 223 : c enddo 224 : c enddo 225 : 225 : CALL MULMAT(AAA,DI2,CO,3,3,3) 225 : CALL MULMAT(COE,AAA,DI2,3,3,3) 226 : ENDIF 227 : C 228 : C---------------------------------- 229 : C DEVIATEUR DE CONTRAINTE EFFECTIVE 230 : C---------------------------------- 231 : C 231 : CALL ZERO(BBB,3,3) 232 : DO I=1,3 233 : BBB(I,I)=-(COE(1,1)+COE(2,2)+COE(3,3))/(3.D0) 234 : ENDDO 234 : CALL PLMATR(COE,UN,BBB,UN,3,3,DCOE) 235 : C --- Calcul de (I-D)**(-1/2).DCOE.(I-D)**(-1/2) --- 236 : IF (endonul) THEN 237 : DO I=1,3 238 : DO J=1,3 239 : DCOT(I,J) = DCOE(I,J) 240 : ENDDO 241 : ENDDO 242 : ELSE 242 : CALL MULMAT(AAA,DI2,DCOE,3,3,3) 242 : CALL MULMAT(DCOT,AAA,DI2,3,3,3) 243 : ENDIF 244 : C 245 : C---------------------------------------------------------------------- 246 : C 247 : C CALCUL DE J1,J2 248 : C J1: CONTRAINTE MOYENNE 249 : C J2: CONTRAINTE DE VON MISES 250 : C 251 : C---------------------------------------------------------------------- 252 : C 253 : AJ1 = (SIG(1)+SIG(2)+SIG(3))/(3.D0) 254 : DO I=1,3 255 : XX(I) = SIG(I)-AJ1 256 : ENDDO 257 : DO I=4,NSTRS 258 : XX(I) = SIG(I) 259 : ENDDO 260 : AJ2 = PROCON(XX,XX,NSTRS) 261 : AJ2 = SQRT(1.5D0*AJ2) 262 : C 263 : C---------------------- 264 : C FORCE THERMODYNAMIQUE 265 : C---------------------- 266 : C 267 : C --- Etat de triaxialité de la contrainte --- 268 : IF ((SIG(1).EQ.0).AND.(SIG(2).EQ.0).AND.(SIG(3).EQ.0)) THEN 269 : tri=1.D0 270 : ELSE 271 : tri = (AJ2)/(SIG(1)+SIG(2)+SIG(3)) 272 : ENDIF 273 : c write(6,*) 'AJ2= ',AJ2 274 : c write(6,*) 'SIG(1)= ',SIG(1) 275 : c write(6,*) 'SIG(2)= ',SIG(2) 276 : c write(6,*) 'SIG(3)= ',SIG(3) 277 : 278 : c write(6,*) 'tri= ',tri 279 : 280 : C --- Valeurs propres et vecteurs propres du tenseur des contraintes motrices --- 280 : CALL JACOB3(COM,3,VaP,VcP) 281 : ER=XMAT(8) 282 : EA=XMAT(9) 283 : DO I=1,3 284 : AA(I) = ((tri*(VaP(I))/(EA))+ABS(tri*(VaP(I))/(EA)))/(2.D0) 285 : AA(I) = (AA(I))**(ER) 286 : ENDDO 287 : c DO I=1,3 288 : c DO J=1,3 289 : c write(6,*) 'COM(',i,',',j,')= ',COM(I,J) 290 : c ENDDO 291 : c ENDDO 292 : C DO J=1,3 293 : C write(6,*) 'VaP(',j,')= ',VaP(J) 294 : C ENDDO 295 : c DO J=1,3 296 : c write(6,*) 'AA(',j,')= ',AA(J) 297 : c ENDDO 298 : c DO I=1,3 299 : c DO J=1,3 300 : c write(6,*) 'VcP(',i,',',j,')= ',VcP(I,J) 301 : c ENDDO 302 : c ENDDO 303 : C 304 : C -- Dans le repère général --- 304 : CALL ZERO(FTH,3,3) 305 : FTH(1,1) = AA(1)*(VcP(1,1))**2+AA(2)*(VcP(1,2))**2 306 : FTH(1,1) = FTH(1,1)+AA(3)*(VcP(1,3))**2 307 : FTH(2,2) = AA(1)*(VcP(2,1))**2+AA(2)*(VcP(2,2))**2 308 : FTH(2,2) = FTH(2,2)+AA(3)*(VcP(2,3))**2 309 : FTH(3,3) = AA(1)*(VcP(3,1))**2+AA(2)*(VcP(3,2))**2 310 : FTH(3,3) = FTH(3,3)+AA(3)*(VcP(3,3))**2 311 : FTH(1,2) = AA(1)*VcP(1,1)*VcP(2,1)+AA(2)*VcP(1,2)*VcP(2,2) 312 : FTH(1,2) = FTH(1,2)+AA(3)*VcP(1,3)*VcP(2,3) 313 : FTH(1,3) = AA(1)*VcP(1,1)*VcP(3,1)+AA(2)*VcP(1,2)*VcP(3,2) 314 : FTH(1,3) = FTH(1,3)+AA(3)*VcP(1,3)*VcP(3,3) 315 : FTH(2,3) = AA(1)*VcP(2,1)*VcP(3,1)+AA(2)*VcP(2,2)*VcP(3,2) 316 : FTH(2,3) = FTH(1,3)+AA(3)*VcP(2,3)*VcP(3,3) 317 : FTH(3,1) = FTH(1,3) 318 : FTH(2,1) = FTH(1,2) 319 : FTH(3,2) = FTH(2,3) 320 : c DO I=1,3 321 : c DO J=1,3 322 : c write(6,*) 'FTH(',i,',',i,')= ',FTH(I,i) 323 : c ENDDO 324 : c ENDDO 325 : C 326 : C--------------------------------------------------------------------- 327 : C 328 : C CALCUL DE LA VITESSE DE DEFORMATION PLASTIQUE CUMULEE VARD(1) 329 : C 330 : C--------------------------------------------------------------------- 331 : C 332 : C --- Paramètres matériau --- 333 : EM=XMAT(4) 334 : EN=XMAT(3) 335 : EK=XMAT(5) 336 : ALF2=XMAT(6) 337 : C --- Calcul de la contrainte équivalente effective --- 338 : AJ3 = (COE(1,1)+COE(2,2)+COE(3,3))/(3.D0) 339 : DO I=1,3 340 : XX(I) = COE(I,I)-AJ3 341 : ENDDO 342 : XX(4) = COE(1,2) 343 : XX(5) = COE(1,3) 344 : XX(6) = COE(2,3) 345 : AJ4 = PROCON(XX,XX,NSTRS) 346 : AJ4 = SQRT(1.5D0*AJ4) 347 : C 348 : C --- Calcul de la norme euclidienne --- 349 : DO I=1,3 350 : DO J=1,3 351 : AAA(I,J)=(DCOT(I,J)*3.D0)/(2.D0*AJ4) 352 : ENDDO 353 : ENDDO 353 : CALL PLMATR(UN2,UN,D1,UNM,3,3,BBB) 353 : CALL INVALM (BBB, 3 , 3 , KER, 1.D-12) 354 : DO I=1,3 355 : DO J=1,3 356 : BBB(I,J)=BBB(I,J)*ALF2/3.D0 357 : ENDDO 358 : ENDDO 358 : CALL ZERO(CCC,3,3) 358 : CALL PLMATR(AAA,UN,BBB,UN,3,3,CCC) 359 : c DO I=1,3 360 : c DO j=1,3 361 : c write(6,*) 'CCC(',i,',',j,')= ', CCC(i,j) 362 : c ENDDO 363 : c ENDDO 364 : 365 : val = 0.D0 366 : DO I=1,3 367 : DO J=1,3 368 : val = val+(CCC(I,J))**(2.D0) 369 : ENDDO 370 : ENDDO 371 : xnor = SQRT(val) 372 : C WRITE(6,*) 'xnor=',xnor 373 : C --- Calcul du premier terme du produit --- 374 : IF ( VAR(1).EQ.0.D0 ) VAR(1)=VAR(1)+1.D-12 375 : ter = (AJ4+(ALF2*AJ3))/(EK*(VAR(1)**(1.D0/EM))) 376 : ter = ter + ABS((AJ4+(ALF2*AJ3))/(EK*(VAR(1)**(1.D0/EM)))) 377 : ter = ter/(2.D0) 378 : ter = ter**(EN) 379 : ter = ter*SQRT(2.D0/3.D0) 380 : c WRITE(6,*) 'ter=',ter 381 : C --- Calcul de l evolution de p --- 382 : IF ((SIG(1).EQ.0).AND.(SIG(2).EQ.0).AND.(SIG(3).EQ.0)) THEN 383 : VARD(1) = 0.D0 384 : ELSE 385 : VARD(1) = ter*xnor 386 : ENDIF 387 : C write(6,*) 'vard(1)=', vard(1) 388 : C 389 : C--------------------------------------------------------------------- 390 : C 391 : C CALCUL DE LA VARIATION DE L ENDOMMAGEMENT VARD(I) (I=2..7) 392 : C 393 : C--------------------------------------------------------------------- 394 : 395 : BET2=XMAT(7) 396 : C --- Dans le repère principal de D --- 396 : CALL ZERO(DV1,3,3) 397 : DV1(1,1) = (BET2*FTH(1,1))+FTH(2,2)+FTH(3,3) 398 : DV1(2,2) = (BET2*FTH(2,2))+FTH(1,1)+FTH(3,3) 399 : DV1(3,3) = (BET2*FTH(3,3))+FTH(2,2)+FTH(1,1) 400 : DV1(1,2) = (BET2-1.D0)*FTH(1,2) 401 : DV1(1,3) = (BET2-1.D0)*FTH(1,3) 402 : DV1(2,3) = (BET2-1.D0)*FTH(2,3) 403 : DV1(2,1) = DV1(1,2) 404 : DV1(3,1) = DV1(1,3) 405 : DV1(3,2) = DV1(2,3) 406 : C 407 : C --- Expressions des VARD(I) (I=2..7) --- 408 : DO I=2,4 409 : VARD(I) = DV1(I-1,I-1) 410 : ENDDO 411 : VARD(5) = DV1(1,2) 412 : VARD(6) = DV1(1,3) 413 : VARD(7) = DV1(2,3) 414 : c DO I=2,4 415 : c write(6,*) 'vard(2)= ', vard(2) 416 : c ENDDO 417 : C 418 : C--------------------------------------------------------------------- 419 : C 420 : C CALCUL DE LA VARIATION DE DEFORMATION PLASTIQUE 421 : C 422 : C--------------------------------------------------------------------- 423 : 424 : IF ((SIG(1).EQ.0).AND.(SIG(2).EQ.0).AND.(SIG(3).EQ.0)) THEN 425 : DO I=1,NSTRS 426 : EPSVD(I) = 0.D0 427 : ENDDO 428 : ELSE 429 : DO I=1,3 430 : EPSVD(I) = CCC(I,I)*VARD(1)*SQRT(3.D0/2.D0)/xnor 431 : ENDDO 432 : EPSVD(4) = CCC(1,2)*VARD(1)*SQRT(3.D0/2.D0)/xnor 433 : * 434 : * AM CORRECTION LE 9/12/09 435 : * 436 : **** IF (IFOUR.GT.0) THEN 437 : IF (IFOURB.GT.0) THEN 438 : EPSVD(5) = CCC(1,3)*VARD(1)*SQRT(3.D0/2.D0)/xnor 439 : EPSVD(6) = CCC(2,3)*VARD(1)*SQRT(3.D0/2.D0)/xnor 440 : ENDIF 441 : ENDIF 442 : 443 : 444 : c DO I=1,3 445 : c write(6,*) 'EPSVD(',i,')= ', EPSVD(i) 446 : c ENDDO 447 : C 448 : C---------------------------------------------------------------------- 449 : C 450 : C CALCUL DE LA VARIATION DES CONTRAINTES 451 : C 452 : C---------------------------------------------------------------------- 453 : 453 : CALL DERTRA(NYOUNG,YOUNG,TI,YUNG,YUNGV,TINF,TSUP) 453 : CALL DERTRA(NXNU,XNU,TI,ENU,ENUV,TINF,TSUP) 454 : C 455 : C --- Calcul de l inverse de la matrice de Hooke au debut du pas --- 455 : CALL ELAST4(2,IFOURB,VAR,NVARI,XMAT,NCOMAT,YUNGV,ENUV, 456 : &XCAR,ICARA,MFR,NSTRS,DDINV,DDV,KERRE,INDIC,ITHHER) 457 : IF(KERRE.NE.0) GO TO 100 458 : C 459 : C --- Calcul de la matrice de Hooke et de sa derivee au debut du pas --- 459 : CALL ELAST4(1,IFOURB,VAR,NVARI,XMAT,NCOMAT,YUNGV,ENUV, 460 : &XCAR,ICARA,MFR,NSTRS,DD,DDV,KERRE,INDIC,ITHHER) 461 : IF(KERRE.NE.0) GO TO 100 462 : C 463 : C --- Détermination de la variation de la matrive de Hooke inversée --- 463 : CALL HOOKP(VAR,VARD,NVARI,XMAT,NCOMAT,MFR,NSTRS,DDINVp) 464 : C 465 : C --- Calcul de la variation de contrainte --- 465 : CALL ZDANUL(SIGD,NSTRS) 466 : DO I=1,3 467 : XX(I) = CO(I,I) 468 : ENDDO 469 : XX(4) = CO(1,2) 470 : XX(5) = CO(1,3) 471 : XX(6) = CO(2,3) 472 : 472 : CALL ZDANUL(EPSED,NSTRS) 473 : DO I=1,NSTRS 474 : EPSED(I)=DSPT(I)-EPSVD(I)-EPSTHD(I) 475 : C write(6,*) 'DSPT(',I,') =', DSPT(I) 476 : c write(6,*) 'EPSVD(',I,') =', EPSVD(I) 477 : c write(6,*) 'EPSTHD(',I,') =', EPSTHD(I) 478 : C write(6,*) 'EPSED(',I,') =', EPSED(I) 479 : ENDDO 480 : 480 : CALL MULMAT(BB,DDINVp,XX,NSTRS,1,NSTRS) 481 : 482 : DO I=1,NSTRS 483 : XX(I)=EPSED(I)-BB(I) 484 : ENDDO 485 : 485 : CALL MULMAT(SIGD,DD,XX,NSTRS,1,NSTRS) 486 : c DO I=1,3 487 : c IF (iGAU.EQ.1) THEN 488 : c write(6,*) 'SIGD(1) =', SIGD(1) 489 : c ENDIF 490 : c ENDDO 491 : C 492 : C---------------------------------------------------------------------- 493 : C 494 : C MISE A JOUR DES VARIABLES INTERNES 495 : C 496 : C---------------------------------------------------------------------- 497 : 498 : VARF(8)=ADMAX 499 : 500 : 100 CONTINUE 501 : RETURN 502 : END 503 : 504 : 505 : 506 :
© Cast3M 2003 - All rights reserved.
Disclaimer