1 : C DAMTG3D SOURCE AM 12/02/27 23:15:30 7289 2 : SUBROUTINE DAMTG3D(parahot3,idimpara3,H66t,rt1,rc1, 3 : . esigma6v,lg,ntot,lerror,i1,i2,i3,i4,i5,i6,H66, 4 : . lfulldamage,lcp) 5 : 6 : c SUBROUTINE added by THG 7 : 8 : c parahot3 is the array that contains the material parameters and state variables 9 : c idimpara3 is the number of parameters for parahot3 10 : c H66t is the tangent matrix 11 : c rt1 and rc1 are the vectors built in DPRAN3D and used in DAMTG3D for calculation of tangent matrix 12 : c esigma6v is the effective stress vector 13 : c lg and lerror are logical 14 : c ntot is the integration point (global number for the whole structure): counter 15 : 16 : c ========================================================================================== 17 : c ! DAMAGE BOX of the 3D plastic-damage concrete model by T. Gernay ! 18 : c ! This subroutine calculates ! 19 : c ! the nominal stress at the end of the time step ! 20 : c ! the nominal consistent algorithmic tangent modulus Ct ! 21 : c ! ! 22 : c ! The effective stress esigma6v and the eff. tangent modulus H66t given by plastic box ! 23 : C ! ! 24 : c ========================================================================================== 25 : 26 : IMPLICIT REAL*8 (A-B,D-H,O-Z) 27 : implicit integer (I-K,M,N) 28 : implicit logical (L) 29 : implicit character*10 (C) 30 : 31 : dimension H66t(6,6) 32 : c Tangent matrix (output) 33 : dimension parahot3(idimpara3) 34 : c parahot3 : Material prop. and various info. at elevated temp. (input and output) 35 : c parahot3(idimpara3-34): equivalent plastic strain in compression at the end of the current step 36 : c parahot3(idimpara3-35): equivalent plastic strain in tension at the end of the current step 37 : c parahot3(idimpara3-6:idimpara3-1) = effective stress obtained by dpran3D 38 : dimension esigma6v(6) 39 : c Effective stress (input) 40 : 41 : dimension rc1(1,6),rt1(1,6) 42 : dimension v1(3),v2(3),v3(3),rcossigmapr(3,3) 43 : dimension p1(3,3),p2(3,3),p3(3,3) 44 : dimension pt1(3,3,3,3),pt2(3,3,3,3),pt3(3,3,3,3) 45 : dimension projt(3,3,3,3),w4(3,3,3,3) 46 : dimension Qtens(3,3,3,3),Qcomp(3,3,3,3) 47 : dimension esigmat(6),esigmac(6) 48 : dimension sigma(6),projtm(6,6),projcm(6,6) 49 : dimension p12(3,3),p23(3,3),p13(3,3) 50 : dimension pp12(3,3,3,3),pp23(3,3,3,3),pp13(3,3,3,3) 51 : dimension W4m(6,6),H66temp(6,6) 52 : dimension sigrt(6,6),sigrc(6,6) 53 : dimension p121(3,3),p231(3,3),p131(3,3) 54 : dimension p122(3,3),p232(3,3),p132(3,3) 55 : dimension damtensor(6,6), damtensor2(6,6),rt1mat(3,3) 56 : dimension rc1mat(3,3),esigmatmat(3,3),esigmacmat(3,3) 57 : dimension sigrttens(3,3,3,3),sigrctens(3,3,3,3),H66(6,6) 58 : dimension Qtensm(6,6),Qcompm(6,6),rc1t(6) 59 : dimension U(6,6) 60 : * dimension id4(3,3,3,3) 61 : 62 : c Matrix of unity I, or U 63 : c ----------------------- 64 : 65 : do jloc=1,6 66 : do iloc=1,6 67 : if (iloc.eq.jloc) then 68 : U(iloc,jloc) = 1.0d0 69 : else 70 : U(iloc,jloc) = 0.0d0 71 : endif 72 : enddo 73 : enddo 74 : 75 : c Fourth identity tensor id4 76 : C do iloc=1,3 77 : C do jloc=1,3 78 : C do kloc=1,3 79 : C do mloc=1,3 80 : C if ((iloc.eq.jloc).and.(kloc.eq.mloc)) then 81 : C id4(iloc,jloc,kloc,mloc) = 1.0d0 82 : C else 83 : C id4(iloc,jloc,kloc,mloc) = 0.0d0 84 : C endif 85 : C end do 86 : C end do 87 : C end do 88 : C end do 89 : 90 : i1 = 1 91 : i2 = 2 92 : i3 = 3 93 : i4 = 4 94 : i5 = 5 95 : i6 = 6 96 : i7 = 7 97 : i8 = 8 98 : i9 = 9 99 : i10=10 100 : i11=11 101 : i12=12 102 : i13=13 103 : i14=14 104 : i15=15 105 : i16=16 106 : i17=17 107 : i18=18 108 : i19=19 109 : i20=20 110 : i21=21 111 : i22=22 112 : i23=23 113 : i24=24 114 : i25=25 115 : i26=26 116 : i27=27 117 : i28=28 118 : i29=29 119 : i30=30 120 : i31=31 121 : i32=32 122 : i33=33 123 : i34=34 124 : i35=35 125 : i36=36 126 : i37=37 127 : i38=38 128 : i39=39 129 : 130 : r0 = 0. 131 : r1 = 1. 132 : 133 : c ================================================ 134 : c ! ! 135 : c ! EIGEN VALUES OF THE EFFECTIVE STRESSES ! 136 : c ! ====================================== ! 137 : c ================================================ 138 : 139 : c Effective stresses 140 : do iloc=i1,i6 141 : ESIGMA6V(iloc) = parahot3(idimpara3-i7+iloc) 142 : end do 143 : 144 : 145 : c Eigen values and eigen vector 145 : call sigmapr3ETC(esigma6v,S1,S2,S3,rcossigmapr,lcp) 146 : do iloc=i1,i3 147 : v1(iloc) = rcossigmapr(iloc,i1) 148 : v2(iloc) = rcossigmapr(iloc,i2) 149 : v3(iloc) = rcossigmapr(iloc,i3) 150 : end do 151 : 152 : c =============================================== 153 : c ! ! 154 : c ! UPDATE DAMAGE VARIABLES ! 155 : c ! ======================= ! 156 : c =============================================== 157 : 158 : rkappac = parahot3(idimpara3-34) 159 : rkappat = parahot3(idimpara3-35) 160 : 161 : c Tensile damage 162 : at = parahot3(16) 163 : dt = 1.0d0 - ((0.5d0)*exp(-at*rkappat)+(0.5d0)* 164 : . exp(-6.0d0*at*rkappat)) 165 : rdti = parahot3(idimpara3-41) 166 : if (dt.lt.rdti) then 167 : dt = rdti 168 : endif 169 : parahot3(idimpara3-33) = dt 170 : 171 : 172 : c Compressive damage 173 : rkappaic = parahot3(idimpara3-36) 174 : rdeltakc = rkappac - rkappaic 175 : rdci = parahot3(idimpara3-40) 176 : rdamprev = dlog(1.0d0-rdci) 177 : dc = 1.0d0 - exp((rdamprev)-parahot3(i14)*rdeltakc) 178 : if (dc.lt.parahot3(idimpara3-40)) dc = parahot3(idimpara3-40) 179 : parahot3(idimpara3-32) = dc 180 : 181 : c =================================================================== 182 : c ! ! 183 : c ! PROJECTION OPERATORS OF THE EFFECTIVE STRESSES P+ AND P- ! 184 : c ! ======================================================== ! 185 : c =================================================================== 186 : 187 : c p1(i3,i3) = v1 X v1 (tensorial product) 187 : call mulABT(v1,v1,p1,i3,i1,i3,i1) 188 : c [3;3]=[3;1]x[1;3] 188 : call mulABT(v2,v2,p2,i3,i1,i3,i1) 189 : c [3;3]=[3;1]x[1;3] 189 : call mulABT(v3,v3,p3,i3,i1,i3,i1) 190 : c [3;3]=[3;1]x[1;3] 191 : 192 : c P+ calculated as the sum = H(sig,i) p,i X p,i 193 : c First we calculate the PT,i = p,i X p,i 194 : do mloc=i1,i3 195 : do kloc=i1,i3 196 : do jloc=i1,i3 197 : do iloc=i1,i3 198 : pt1(iloc,jloc,kloc,mloc) = p1(iloc,jloc)*p1(kloc,mloc) 199 : c [3;3;3;3] = [3;3] x [3;3] tensorial product 200 : pt2(iloc,jloc,kloc,mloc) = p2(iloc,jloc)*p2(kloc,mloc) 201 : pt3(iloc,jloc,kloc,mloc) = p3(iloc,jloc)*p3(kloc,mloc) 202 : end do 203 : end do 204 : end do 205 : end do 206 : 207 : if (S1.gt.r0) then 208 : h1 = r1 209 : else 210 : h1 = r0 211 : endif 212 : if (S2.gt.r0) then 213 : h2 = r1 214 : else 215 : h2 = r0 216 : endif 217 : if (S3.gt.r0) then 218 : h3 = r1 219 : else 220 : h3 = r0 221 : endif 222 : 223 : c Projection operators of the effective stresses P+ [3;3;3;3] 224 : do mloc=i1,i3 225 : do kloc=i1,i3 226 : do jloc=i1,i3 227 : do iloc=i1,i3 228 : projt(iloc,jloc,kloc,mloc) = h1*pt1(iloc,jloc,kloc,mloc)+ 229 : . h2*pt2(iloc,jloc,kloc,mloc) + h3*pt3(iloc,jloc,kloc,mloc) 230 : end do 231 : end do 232 : end do 233 : end do 234 : 235 : c Transform the [3;3;3;3] tensor into [6;6] matrix (Voigt notation) 236 : c projtm is the tensile projection matrix [6;6] of the effective stress 237 : do jloc=i1,i3 238 : do iloc=i1,i3 239 : projtm(iloc,jloc) = projt(iloc,iloc,jloc,jloc) 240 : end do 241 : end do 242 : do iloc=i1,i3 243 : projtm(iloc,i4) = 0.5d0*(projt(iloc,iloc,i1,i2) + 244 : . projt(iloc,iloc,i2,i1)) 245 : projtm(iloc,i5) = 0.5d0*(projt(iloc,iloc,i1,i3) + 246 : . projt(iloc,iloc,i3,i1)) 247 : c 5th component is eps,xz 248 : projtm(iloc,i6) = 0.5d0*(projt(iloc,iloc,i2,i3) + 249 : . projt(iloc,iloc,i3,i2)) 250 : c 6th component is eps,yz 251 : end do 252 : do jloc=i1,i3 253 : projtm(i4,jloc) = 0.5d0*(projt(i1,i2,jloc,jloc) + 254 : . projt(i2,i1,jloc,jloc)) 255 : projtm(i5,jloc) = 0.5d0*(projt(i1,i3,jloc,jloc) + 256 : . projt(i3,i1,jloc,jloc)) 257 : projtm(i6,jloc) = 0.5d0*(projt(i2,i3,jloc,jloc) + 258 : . projt(i3,i2,jloc,jloc)) 259 : 260 : end do 261 : projtm(i4,i4) = 0.5d0*(projt(i1,i2,i1,i2) + 262 : . projt(i1,i2,i2,i1)) 263 : projtm(i4,i5) = 0.5d0*(projt(i1,i2,i1,i3) + 264 : . projt(i1,i2,i3,i1)) 265 : projtm(i4,i6) = 0.5d0*(projt(i1,i2,i2,i3) + 266 : . projt(i1,i2,i3,i2)) 267 : projtm(i5,i4) = 0.5d0*(projt(i1,i3,i1,i2) + 268 : . projt(i1,i3,i2,i1)) 269 : projtm(i5,i5) = 0.5d0*(projt(i1,i3,i1,i3) + 270 : . projt(i1,i3,i3,i1)) 271 : projtm(i6,i5) = 0.5d0*(projt(i2,i3,i1,i3) + 272 : . projt(i2,i3,i3,i1)) 273 : projtm(i6,i4) = 0.5d0*(projt(i2,i3,i1,i2) + 274 : . projt(i2,i3,i2,i1)) 275 : projtm(i5,i6) = 0.5d0*(projt(i1,i3,i2,i3) + 276 : . projt(i1,i3,i3,i2)) 277 : projtm(i6,i6) = 0.5d0*(projt(i2,i3,i2,i3) + 278 : . projt(i2,i3,i3,i2)) 279 : 280 : c Identity matrix [6,6] U used to calculate P- and Q- 281 : c for P-: not really necessary as sig- = sig - sig+ 282 : c for Q-: necessary for calculation of consistent algorit. tangent 283 : c Projection operators of the effective stresses P- 284 : do jloc=i1,i6 285 : do iloc=i1,i6 286 : projcm(iloc,jloc) = U(iloc,jloc) - projtm(iloc,jloc) 287 : end do 288 : end do 289 : 290 : c =================================================================== 291 : c ! ! 292 : c ! POSITIVE AND NEGATIVE PARTS OF THE EFFECTIVE STRESS TENSOR ! 293 : c ! ========================================================== ! 294 : c =================================================================== 295 : 296 : c projtm has to be transposed 296 : call mulATB(projtm,esigma6v,esigmat,i6,i6,i6,i1) 297 : c [6;1] = [6;6] x [6x1] 298 : do iloc=i1,i6 299 : esigmac(iloc) = esigma6v(iloc) - esigmat(iloc) 300 : end do 301 : 302 : c Save the negative part of the effective stress tensor to be used in 303 : c the transient creep strain calculation 304 : do iloc=i1,i6 305 : parahot3(23+iloc) = esigmac(iloc) 306 : end do 307 : 308 : c ======================================== 309 : c ! ! 310 : c ! NOMINAL STRESSES ! 311 : c ! ================ ! 312 : c ======================================== 313 : 314 : c SIGMA = (1-dt) * SIG,eff,tension + (1-dc) * SIG,eff,comp 315 : 316 : c concrete completely damaged 317 : if ((lfulldamage).or.(dc.gt.0.98).or.(dt.gt.0.98)) then 318 : do iloc=1,6 319 : SIGMA(iloc) = 0.d0 320 : end do 321 : do iloc=i1,i6 322 : parahot3(idimpara3-i13+iloc) = sigma(iloc) 323 : end do 324 : do jloc=1,6 325 : do iloc=1,6 326 : H66t(iloc,jloc) = 0.d0 327 : end do 328 : end do 329 : return 330 : endif 331 : 332 : dt = parahot3(idimpara3-33) 333 : do iloc = i1,i6 334 : SIGMA(iloc) = (1.0d0-dt)*esigmat(iloc) + 335 : . (1.0d0-dc)*esigmac(iloc) 336 : end do 337 : 338 : do iloc=i1,i6 339 : parahot3(idimpara3-i13+iloc) = sigma(iloc) 340 : end do 341 : 342 : c *************************************************************** 343 : c ======================================== 344 : c ! ! 345 : c ! TANGENT MODULUS ! 346 : c ! =============== ! 347 : c ======================================== 348 : 349 : c Calculation of (nominal) consistent algorithmic tangent modulus Ct. 350 : c In DPRan3D we have calculated the effective elastoplastic tangent 351 : c modulus Dt. Now, we have to include the damage contribution. 352 : 353 : c APEX of Drucker-Prager or equivalent stress of Rankine = 0 354 : c First, verify that we are not at the apex of the surfaces 355 : c If we are, Dt = 0 and in this case we should have Ct = 0 too 356 : rloc = 0.d0 357 : do jloc=1,6 358 : do iloc=1,6 359 : rloc = rloc+dabs(H66t(iloc,jloc)) 360 : enddo 361 : enddo 362 : if (rloc.lt.1.e+4) then 363 : do jloc=1,6 364 : do iloc=1,6 365 : H66t(iloc,jloc) = 0.d0 366 : end do 367 : end do 368 : return 369 : endif 370 : c End of test APEX 371 : 372 : do iloc=i1,6 373 : rt1(i1,iloc) = ( (at/2.)*dexp(-at*rkappat)+(6.*at/2.)* 374 : . dexp(-6.*at*rkappat) ) * rt1(i1,iloc) 375 : rc1(i1,iloc) = ( parahot3(14)*dexp(-parahot3(14)*rkappac) ) 376 : . * rc1(i1,iloc) 377 : end do 378 : 379 : c Calculation of W4 380 : c Calculation of p12, p23, p13 381 : c p12(i3,i3) = (1/2) (v1 X v2 + v2 X v1) (tensorial product) 381 : call mulABT(v1,v2,p121,i3,i1,i3,i1) 382 : c [3;3]=[3;1]x[1;3] 382 : call mulABT(v2,v1,p122,i3,i1,i3,i1) 382 : call mulABT(v2,v3,p231,i3,i1,i3,i1) 382 : call mulABT(v3,v2,p232,i3,i1,i3,i1) 382 : call mulABT(v1,v3,p131,i3,i1,i3,i1) 382 : call mulABT(v3,v1,p132,i3,i1,i3,i1) 383 : do jloc=i1,i3 384 : do iloc=i1,i3 385 : p12(iloc,jloc) = (1.0d0/2.0d0) * (p121(iloc,jloc) + 386 : . p122(iloc,jloc)) 387 : p23(iloc,jloc) = (1.0d0/2.0d0) * (p231(iloc,jloc) + 388 : . p232(iloc,jloc)) 389 : p13(iloc,jloc) = (1.0d0/2.0d0) * (p131(iloc,jloc) + 390 : . p132(iloc,jloc)) 391 : end do 392 : end do 393 : 394 : c Calculation of PPijkl = pij X pij 395 : c [3;3;3;3]=[3;3] X [3;3] 396 : do mloc=i1,i3 397 : do kloc=i1,i3 398 : do jloc=i1,i3 399 : do iloc=i1,i3 400 : pp12(iloc,jloc,kloc,mloc)=p12(iloc,jloc)*p12(kloc,mloc) 401 : pp23(iloc,jloc,kloc,mloc)=p23(iloc,jloc)*p23(kloc,mloc) 402 : pp13(iloc,jloc,kloc,mloc)=p13(iloc,jloc)*p13(kloc,mloc) 403 : end do 404 : end do 405 : end do 406 : end do 407 : 408 : c Calculation of Q+ ! [3;3;3;3] 409 : do mloc=i1,i3 410 : do kloc=i1,i3 411 : do jloc=i1,i3 412 : do iloc=i1,i3 413 : Qtens(iloc,jloc,kloc,mloc) = r0 414 : end do 415 : end do 416 : end do 417 : end do 418 : if ((S1-S2).ne.r0) then 419 : rloc = (h1*S1) - (h2*S2) 420 : do mloc=i1,i3 421 : do kloc=i1,i3 422 : do jloc=i1,i3 423 : do iloc=i1,i3 424 : Qtens(iloc,jloc,kloc,mloc)=Qtens(iloc,jloc,kloc,mloc) 425 : . + 2.0d0 * (rloc/(S1-S2)) 426 : . * pp12(iloc,jloc,kloc,mloc) 427 : end do 428 : end do 429 : end do 430 : end do 431 : else 432 : continue 433 : endif 434 : if ((S2-S3).ne.r0) then 435 : rloc = (h2*S2) - (h3*S3) 436 : do mloc=i1,i3 437 : do kloc=i1,i3 438 : do jloc=i1,i3 439 : do iloc=i1,i3 440 : Qtens(iloc,jloc,kloc,mloc)=Qtens(iloc,jloc,kloc,mloc) 441 : . + 2.0d0 * (rloc/(S2-S3)) 442 : . * pp23(iloc,jloc,kloc,mloc) 443 : end do 444 : end do 445 : end do 446 : end do 447 : else 448 : continue 449 : endif 450 : if ((S1-S3).ne.r0) then 451 : rloc = (h1*S1) - (h3*S3) 452 : do mloc=i1,i3 453 : do kloc=i1,i3 454 : do jloc=i1,i3 455 : do iloc=i1,i3 456 : Qtens(iloc,jloc,kloc,mloc)=Qtens(iloc,jloc,kloc,mloc) 457 : . + 2.0d0 * (rloc/(S1-S3)) 458 : . * pp13(iloc,jloc,kloc,mloc) 459 : end do 460 : end do 461 : end do 462 : end do 463 : else 464 : continue 465 : endif 466 : do mloc=i1,i3 467 : do kloc=i1,i3 468 : do jloc=i1,i3 469 : do iloc=i1,i3 470 : Qtens(iloc,jloc,kloc,mloc) = projt(iloc,jloc,kloc,mloc) 471 : . + Qtens(iloc,jloc,kloc,mloc) 472 : c [3;3;3;3] 473 : end do 474 : end do 475 : end do 476 : end do 477 : 478 : c Calculation of the matrix [6;6] Qtens 479 : c Transform the [3;3;3;3] Qtens tensor into [6;6] matrix (Voigt notation) 480 : do jloc=1,3 481 : do iloc=1,3 482 : Qtensm(iloc,jloc) = Qtens(iloc,iloc,jloc,jloc) 483 : end do 484 : end do 485 : do iloc=1,3 486 : Qtensm(iloc,4) = 0.5d0*(Qtens(iloc,iloc,1,2) 487 : . + Qtens(iloc,iloc,2,1)) 488 : Qtensm(iloc,5) = 0.5d0*(Qtens(iloc,iloc,1,3) 489 : . + Qtens(iloc,iloc,3,1)) 490 : Qtensm(iloc,6) = 0.5d0*(Qtens(iloc,iloc,2,3) 491 : . + Qtens(iloc,iloc,3,2)) 492 : end do 493 : do jloc=1,3 494 : Qtensm(4,jloc) = 0.5d0*(Qtens(1,2,jloc,jloc) 495 : . + Qtens(2,1,jloc,jloc)) 496 : Qtensm(5,jloc) = 0.5d0*(Qtens(1,3,jloc,jloc) 497 : . + Qtens(3,1,jloc,jloc)) 498 : Qtensm(6,jloc) = 0.5d0*(Qtens(2,3,jloc,jloc) 499 : . + Qtens(3,2,jloc,jloc)) 500 : end do 501 : Qtensm(4,4) = 0.5d0*(Qtens(1,2,1,2) + Qtens(1,2,2,1)) 502 : Qtensm(4,5) = 0.5d0*(Qtens(1,2,1,3) + Qtens(1,2,3,1)) 503 : Qtensm(4,6) = 0.5d0*(Qtens(1,2,2,3) + Qtens(1,2,3,2)) 504 : Qtensm(5,4) = 0.5d0*(Qtens(1,3,1,2) + Qtens(1,3,2,1)) 505 : Qtensm(5,5) = 0.5d0*(Qtens(1,3,1,3) + Qtens(1,3,3,1)) 506 : Qtensm(5,6) = 0.5d0*(Qtens(1,3,2,3) + Qtens(1,3,3,2)) 507 : Qtensm(6,4) = 0.5d0*(Qtens(2,3,1,2) + Qtens(2,3,2,1)) 508 : Qtensm(6,5) = 0.5d0*(Qtens(2,3,1,3) + Qtens(2,3,3,1)) 509 : Qtensm(6,6) = 0.5d0*(Qtens(2,3,2,3) + Qtens(2,3,3,2)) 510 : 511 : do jloc=1,6 512 : do iloc=1,6 513 : Qcompm(iloc,jloc) = U(iloc,jloc) - Qtensm(iloc,jloc) 514 : enddo 515 : enddo 516 : do iloc=1,3 517 : Qcompm(iloc,4) = 0.5d0*U(iloc,4) - Qtensm(iloc,4) 518 : Qcompm(iloc,5) = 0.5d0*U(iloc,5) - Qtensm(iloc,5) 519 : Qcompm(iloc,6) = 0.5d0*U(iloc,6) - Qtensm(iloc,6) 520 : enddo 521 : do jloc=1,3 522 : Qcompm(4,jloc) = 0.5d0*U(4,jloc) - Qtensm(4,jloc) 523 : Qcompm(5,jloc) = 0.5d0*U(5,jloc) - Qtensm(5,jloc) 524 : Qcompm(6,jloc) = 0.5d0*U(6,jloc) - Qtensm(6,jloc) 525 : enddo 526 : do jloc=4,6 527 : do iloc=4,6 528 : Qcompm(iloc,jloc) = 0.5d0*U(iloc,jloc) - Qtensm(iloc,jloc) 529 : enddo 530 : enddo 531 : 532 : do jloc=i1,6 533 : do iloc=i1,6 534 : w4m(iloc,jloc) = dt*Qtensm(iloc,jloc) + dc*Qcompm(iloc,jloc) 535 : end do 536 : end do 537 : 538 : c Calculation of (nominal) consistent algorithmic tangent modulus Ct. 539 : c [6;6] 540 : c Ct = (U - w4m) * Dt - sig,eff,+ * rt1 - sig,eff,- * rc1 540 : c In this subroutine Dt is called H66t (comes from DPRan3d) 541 : c and at the end, Ct has to be stored into H66t 542 : 543 : c Calculation of the product (U - w4m) * Dt 544 : do jloc=i1,i6 545 : do iloc=i1,i6 546 : w4m(iloc,jloc) = U(iloc,jloc) - w4m(iloc,jloc) 547 : end do 548 : end do 548 : call mulAB(w4m,H66t,H66temp,i6,i6,i6,i6) 549 : c [6;6] 550 : do j=1,6 551 : do i=1,6 552 : H66t(i,j) = H66temp(i,j) 553 : enddo 554 : enddo 555 : 556 : c esigmat and esigmac have been calculated above 556 : call mulAB(esigmat,rt1,sigrt,i6,i1,i1,i6) 557 : c [6;6] = [6;1] * [1;6] 557 : call mulAB(esigmac,rc1,sigrc,i6,i1,i1,i6) 558 : 559 : do jloc=i1,i6 560 : do iloc=i1,i6 561 : H66t(iloc,jloc) = H66t(iloc,jloc) 562 : . - sigrt(iloc,jloc) 563 : . - sigrc(iloc,jloc) 564 : end do 565 : end do 566 : 567 : c END OF TANGENT MODULUS 568 : 569 : RETURN 570 : END 571 : 572 : 573 : 574 : 575 : 576 : 577 : 578 : 579 : 580 : 581 : 582 : 583 :
© Cast3M 2003 - All rights reserved.
Disclaimer