1! 2! CalculiX - A 3-dimensional finite element program 3! Copyright (C) 1998-2021 Guido Dhondt 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU General Public License as 7! published by the Free Software Foundation(version 2); 8! 9! 10! This program is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU General Public License for more details. 14! 15! You should have received a copy of the GNU General Public License 16! along with this program; if not, write to the Free Software 17! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18! 19 subroutine umat_lin_el_corot( 20 & amat,iel,iint,kode,elconloc,emec,emec0, 21 & beta,xokl,voj,xkl,vj,ithermal,t1l,dtime,time,ttime, 22 & icmd,ielas,mi,nstate_,xstateini,xstate,stre,stiff, 23 & iorien,pgauss,orab,eloc,nlgeom_undo) 24! 25! calculates stiffness and stresses for a user defined material 26! law 27! 28! icmd=3: calculates stress at mechanical strain 29! else: calculates stress at mechanical strain and the stiffness 30! matrix 31! 32! INPUT: 33! 34! amat material name 35! iel element number 36! iint integration point number 37! 38! kode material type (-100-#of constants entered 39! under *USER MATERIAL): can be used for materials 40! with varying number of constants 41! 42! elconloc(21) user defined constants defined by the keyword 43! card *USER MATERIAL (max. 21, actual # = 44! -kode-100), interpolated for the 45! actual temperature t1l 46! 47! emec(6) Lagrange mechanical strain tensor (component order: 48! 11,22,33,12,13,23) at the end of the increment 49! (thermal strains are subtracted) 50! emec0(6) Lagrange mechanical strain tensor at the start of the 51! increment (thermal strains are subtracted) 52! beta(6) residual stress tensor (the stress entered under 53! the keyword *INITIAL CONDITIONS,TYPE=STRESS) 54! 55! xokl(3,3) deformation gradient at the start of the increment 56! voj Jacobian at the start of the increment 57! xkl(3,3) deformation gradient at the end of the increment 58! vj Jacobian at the end of the increment 59! 60! ithermal 0: no thermal effects are taken into account 61! >0: thermal effects are taken into account (triggered 62! by the keyword *INITIAL CONDITIONS,TYPE=TEMPERATURE) 63! t1l temperature at the end of the increment 64! dtime time length of the increment 65! time step time at the end of the current increment 66! ttime total time at the start of the current step 67! 68! icmd not equal to 3: calculate stress and stiffness 69! 3: calculate only stress 70! ielas 0: no elastic iteration: irreversible effects 71! are allowed 72! 1: elastic iteration, i.e. no irreversible 73! deformation allowed 74! 75! mi(1) max. # of integration points per element in the 76! model 77! nstate_ max. # of state variables in the model 78! 79! xstateini(nstate_,mi(1),# of elements) 80! state variables at the start of the increment 81! xstate(nstate_,mi(1),# of elements) 82! state variables at the end of the increment 83! 84! stre(6) Piola-Kirchhoff stress of the second kind 85! at the start of the increment 86! 87! iorien number of the local coordinate axis system 88! in the integration point at stake (takes the value 89! 0 if no local system applies) 90! pgauss(3) global coordinates of the integration point 91! orab(7,*) description of all local coordinate systems. 92! If a local coordinate system applies the global 93! tensors can be obtained by premultiplying the local 94! tensors with skl(3,3). skl is determined by calling 95! the subroutine transformatrix: 96! call transformatrix(orab(1,iorien),pgauss,skl) 97! eloc(6) Lagrange total strain tensor (component order: 98! 11,22,33,12,13,23) at the end of the increment 99! (thermal strains are subtracted) 100! 101! 102! OUTPUT: 103! 104! xstate(nstate_,mi(1),# of elements) 105! updated state variables at the end of the increment 106! stre(6) Piola-Kirchhoff stress of the second kind at the 107! end of the increment 108! stiff(21): consistent tangent stiffness matrix in the material 109! frame of reference at the end of the increment. In 110! other words: the derivative of the PK2 stress with 111! respect to the Lagrangian strain tensor. The matrix 112! is supposed to be symmetric, only the upper half is 113! to be given in the same order as for a fully 114! anisotropic elastic material (*ELASTIC,TYPE=ANISO). 115! Notice that the matrix is an integral part of the 116! fourth order material tensor, i.e. the Voigt notation 117! is not used. 118! eloc(6) linear total strain tensor for large 119! rotations (component order: 120! 11,22,33,12,13,23) at the end of the increment 121! (thermal strains are subtracted) 122! nlgeom_undo 0: Lagrange strain goes out 123! 1: linear strain for large rotations goes out 124! 125 implicit none 126! 127 character*80 amat 128! 129 integer ithermal(*),icmd,kode,ielas,iel,iint,nstate_,mi(*),iorien, 130 & kel(4,21),n,matz,ier,i,j,kal(2,6),j1,j2,j3,j4,nlgeom_undo, 131 & nconstants 132! 133 real*8 elconloc(*),stiff(21),emec(6),emec0(6),beta(6),stre(6), 134 & vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*), 135 & time,ttime,pnewdt,xstate(nstate_,mi(1),*),e(3,3),we(3), 136 & xstateini(nstate_,mi(1),*),wc(3),z(3,3),fv1(3),fv2(3),eps, 137 & pi,young,fla(3),xm1(3,3),xm2(3,3),xm3(3,3),dfla(3),a(21),b(21), 138 & d(3,3),c(3,3),xmm1(21),xmm2(21),xmm3(21),ca(3),cb(3),u(6), 139 & dude(21),um,un,al,am1,am2,ee,eloc(6),eth(6),elin(6) 140! 141 kal=reshape((/1,1,2,2,3,3,1,2,1,3,2,3/),(/2,6/)) 142! 143 kel=reshape((/1,1,1,1,1,1,2,2,2,2,2,2,1,1,3,3,2,2,3,3,3,3,3,3, 144 & 1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3, 145 & 3,3,1,3,1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3, 146 & 1,2,2,3,1,3,2,3,2,3,2,3/),(/4,21/)) 147! 148 d=reshape((/1.d0,0.d0,0.d0,0.d0,1.d0,0.d0,0.d0,0.d0,1.d0/), 149 & (/3,3/)) 150! 151 nlgeom_undo=1 152! 153! corotational linear elastic material. The constants are the 154! the same as for linear elastic material (2 for isotropic, 155! 9 for orthotropic and 21 for completely anisotropic material) 156! 157! PK2=(Hooke matrix)*(U-I) 158! 159! U is the right stretch tensor = sum sqrt(Lambda_i)*M_i 160! = square root of C 161! 162 pi=4.d0*datan(1.d0) 163! 164c if(elconloc(1).lt.1.d-30) then 165c write(*,*) '*ERROR in umat_lin_el_corot: the Young modulus' 166c write(*,*) ' is too small' 167c call exit(201) 168c endif 169c! 170c if(elconloc(2).lt.1.d-30) then 171c write(*,*) '*ERROR in umat_lin_el_corot: maximum pressure' 172c write(*,*) ' value is too small' 173c call exit(201) 174c endif 175c young=elconloc(1) 176c eps=elconloc(2)*pi/young 177! 178! calculation of the eigenvalues and eigenvectors of the 179! Lagrange strain 180! 181 e(1,1)=emec(1) 182 e(2,2)=emec(2) 183 e(3,3)=emec(3) 184 e(1,2)=emec(4) 185 e(1,3)=emec(5) 186 e(2,3)=emec(6) 187 e(2,1)=emec(4) 188 e(3,1)=emec(5) 189 e(3,2)=emec(6) 190c if(iint.eq.1) write(*,*) 'emec' 191c if(iint.eq.1) write(*,100) (emec(i),i=1,6) 192! 193 n=3 194 matz=1 195 ier=0 196! 197 call rs(n,n,e,we,matz,z,fv1,fv2,ier) 198! 199 if(ier.ne.0) then 200 write(*,*) ' 201 & *ERROR calculating the eigenvalues/vectors in umat_tension' 202 call exit(201) 203 endif 204! 205! calculating the eigenvalues of the Cauchy tensor 206! at the end of the increment (eigenvalues of the Cauchy tensor) 207! 208 do i=1,3 209 wc(i)=2.d0*we(i)+1.d0 210c if(iint.eq.1) 211c & write(*,*) 'eigenvalues',i,wc(i),z(1,i),z(2,i),z(3,i) 212 enddo 213! 214! check for 3 equal eigenvalues 215! 216 if((dabs(we(3)-we(2)).lt.1.d-10).and. 217 & (dabs(we(2)-we(1)).lt.1.d-10)) then 218 fla(1)=dsqrt(wc(1)) 219 do i=1,3 220 u(i)=fla(1) 221 enddo 222 do i=4,6 223 u(i)=0.d0 224 enddo 225c if(iint.eq.1) write(*,*) '1 diff u' 226c if(iint.eq.1) write(*,100) (u(i),i=1,6) 227! 228 if(icmd.ne.3) then 229 dfla(1)=1.d0/(2.d0*fla(1)) 230! 231 do j=1,21 232 j1=kel(1,j) 233 j2=kel(2,j) 234 j3=kel(3,j) 235 j4=kel(4,j) 236 dude(j)=dfla(1)*(d(j3,j1)*d(j4,j2)+d(j3,j2)*d(j4,j1)) 237 enddo 238c if(iint.eq.1) write(*,*) '1 diff dude' 239c if(iint.eq.1) write(*,100) (dude(i),i=1,21) 240 endif 241! 242! 2 equal eigenvalues: w2=w3 243! 244 elseif(dabs(we(3)-we(2)).lt.1.d-10) then 245! 246! calculating the structural tensor for the unequal 247! eigenvalue of the Cauchy and the Lagrange tensor 248! 249 do i=1,3 250 do j=1,3 251 xm1(i,j)=z(i,1)*z(j,1) 252 enddo 253 enddo 254! 255! calculating the modified Young's modulus (due to the 256! tension-only modification) 257! 258 do i=1,2 259 fla(i)=dsqrt(wc(i)) 260 enddo 261! 262! calculating the right stretch tensor 263! 264 do j=1,6 265 j1=kal(1,j) 266 j2=kal(2,j) 267 u(j)=fla(1)*xm1(j1,j2)+fla(2)*(d(j1,j2)-xm1(j1,j2)) 268 enddo 269c if(iint.eq.1) write(*,*) '2 diff u (w2=w3)' 270c if(iint.eq.1) write(*,100) (u(i),i=1,6) 271! 272 if(icmd.ne.3) then 273! 274! calculating the matrix dU/dE 275! 276! derivative of the modified young's modulus w.r.t. 277! the Cauchy eigenvalues 278! 279 do i=1,2 280 dfla(i)=1.d0/(2.d0*fla(i)) 281 enddo 282! 283! matrix dU/dE 284! 285 do j=1,21 286 j1=kel(1,j) 287 j2=kel(2,j) 288 j3=kel(3,j) 289 j4=kel(4,j) 290! 291! calculating the auxiliary field xmm1 292! 293 xmm1(j)=xm1(j1,j2)*xm1(j3,j4) 294! 295! calculating the auxiliary fields a and b 296! 297 a(j)=(d(j3,j1)*d(j4,j2)+d(j3,j2)*d(j4,j1))/2.d0 298 b(j)=(d(j3,j1)*xm1(j4,j2)+d(j4,j1)*xm1(j3,j2)+ 299 & d(j4,j2)*xm1(j1,j3)+d(j3,j2)*xm1(j1,j4))/2.d0 300! 301! calculating the matrix dU/dE 302! 303 dude(j)=2.d0*(dfla(1)*xmm1(j) 304 & +dfla(2)*(a(j)+xmm1(j)-b(j)) 305 & +(fla(1)-fla(2))*(b(j)-2.d0*xmm1(j)) 306 & /(2.d0*(we(1)-we(2)))) 307 enddo 308c if(iint.eq.1) write(*,*) '2 diff dude (w2=w3)' 309c if(iint.eq.1) write(*,100) (dude(i),i=1,21) 310 endif 311! 312! 2 equal eigenvalues: w1=w2 313! 314 elseif(dabs(we(2)-we(1)).lt.1.d-10) then 315! 316! calculating the structural tensor for the unequal 317! eigenvalue of the Cauchy and the Lagrange tensor 318! 319 do i=1,3 320 do j=1,3 321 xm3(i,j)=z(i,3)*z(j,3) 322 enddo 323 enddo 324! 325! calculating the modified Young's modulus (due to the 326! tension-only modification) 327! 328 do i=2,3 329 fla(i)=dsqrt(wc(i)) 330 enddo 331! 332! calculating the right stretch tensor 333! 334 do j=1,6 335 j1=kal(1,j) 336 j2=kal(2,j) 337 u(j)=fla(3)*xm3(j1,j2)+fla(2)*(d(j1,j2)-xm3(j1,j2)) 338 enddo 339c if(iint.eq.1) write(*,*) '2 diff u (w1=w2)' 340c if(iint.eq.1) write(*,100) (u(i),i=1,6) 341! 342 if(icmd.ne.3) then 343! 344! calculating the matrix dU/dE 345! 346! derivative of the modified young's modulus w.r.t. 347! the Cauchy eigenvalues 348! 349 do i=2,3 350 dfla(i)=1.d0/(2.d0*fla(i)) 351 enddo 352! 353! matrix dU/dE 354! 355 do j=1,21 356 j1=kel(1,j) 357 j2=kel(2,j) 358 j3=kel(3,j) 359 j4=kel(4,j) 360! 361! calculating the auxiliary field xmm3 362! 363 xmm3(j)=xm3(j1,j2)*xm3(j3,j4) 364! 365! calculating the auxiliary fields a and b 366! 367 a(j)=(d(j3,j1)*d(j4,j2)+d(j3,j2)*d(j4,j1))/2.d0 368 b(j)=(d(j3,j1)*xm3(j4,j2)+d(j4,j1)*xm3(j3,j2)+ 369 & d(j4,j2)*xm3(j1,j3)+d(j3,j2)*xm3(j1,j4))/2.d0 370! 371! calculating the matrix dU/dE 372! 373 dude(j)=2.d0*(dfla(3)*xmm3(j) 374 & +dfla(2)*(a(j)+xmm3(j)-b(j)) 375 & +(fla(2)-fla(3))*(b(j)-2.d0*xmm3(j)) 376 & /(2.d0*(we(2)-we(3)))) 377 enddo 378c if(iint.eq.1) write(*,*) '2 diff dude (w1=w2)' 379c if(iint.eq.1) write(*,100) (dude(i),i=1,21) 380 endif 381! 382! 2 equal eigenvalues: w1=w3 383! 384 elseif(dabs(we(3)-we(1)).lt.1.d-10) then 385! 386! calculating the structural tensor for the unequal 387! eigenvalue of the Cauchy and the Lagrange tensor 388! 389 do i=1,3 390 do j=1,3 391 xm3(i,j)=z(i,3)*z(j,3) 392 enddo 393 enddo 394! 395! calculating the modified Young's modulus (due to the 396! tension-only modification) 397! 398 do i=2,3 399 fla(i)=dsqrt(wc(i)) 400 enddo 401! 402! calculating the right stretch tensor 403! 404 do j=1,6 405 j1=kal(1,j) 406 j2=kal(2,j) 407 u(j)=fla(2)*xm3(j1,j2)+fla(3)*(d(j1,j2)-xm3(j1,j2)) 408 enddo 409c if(iint.eq.1) write(*,*) '2 diff u (w1=w2)' 410c if(iint.eq.1) write(*,100) (u(i),i=1,6) 411! 412 if(icmd.ne.3) then 413! 414! calculating the matrix dU/dE 415! 416! derivative of the modified young's modulus w.r.t. 417! the Cauchy eigenvalues 418! 419 do i=2,3 420 dfla(i)=1.d0/(2.d0*fla(i)) 421 enddo 422! 423! matrix dU/dE 424! 425 do j=1,21 426 j1=kel(1,j) 427 j2=kel(2,j) 428 j3=kel(3,j) 429 j4=kel(4,j) 430! 431! calculating the auxiliary field xmm3 432! 433 xmm3(j)=xm3(j1,j2)*xm3(j3,j4) 434! 435! calculating the auxiliary fields a and b 436! 437 a(j)=(d(j3,j1)*d(j4,j2)+d(j3,j2)*d(j4,j1))/2.d0 438 b(j)=(d(j3,j1)*xm3(j4,j2)+d(j4,j1)*xm3(j3,j2)+ 439 & d(j4,j2)*xm3(j1,j3)+d(j3,j2)*xm3(j1,j4))/2.d0 440! 441! calculating the matrix dU/dE 442! 443 dude(j)=2.d0*(dfla(2)*xmm3(j) 444 & +dfla(3)*(a(j)+xmm3(j)-b(j)) 445 & +(fla(3)-fla(2))*(b(j)-2.d0*xmm3(j)) 446 & /(2.d0*(we(3)-we(2)))) 447 enddo 448c if(iint.eq.1) write(*,*) '2 diff dude (w1=w2)' 449c if(iint.eq.1) write(*,100) (dude(i),i=1,21) 450 endif 451! 452! 3 different eigenvalues 453! 454 else 455! 456! calculating the structural tensors of the Cauchy and the 457! Lagrange tensor 458! 459 do i=1,3 460 do j=1,3 461 xm1(i,j)=z(i,1)*z(j,1) 462 xm2(i,j)=z(i,2)*z(j,2) 463 xm3(i,j)=z(i,3)*z(j,3) 464 enddo 465 enddo 466! 467! calculating the modified Young's modulus (due to the 468! tension-only modification) 469! 470 do i=1,3 471 fla(i)=dsqrt(wc(i)) 472 enddo 473! 474! calculating the right stretch tensor 475! 476 do j=1,6 477 j1=kal(1,j) 478 j2=kal(2,j) 479 u(j)=fla(1)*xm1(j1,j2)+fla(2)*xm2(j1,j2) 480 & +fla(3)*xm3(j1,j2) 481 enddo 482c if(iint.eq.1) write(*,*) '3 diff u' 483c if(iint.eq.1) write(*,100) (u(i),i=1,6) 484! 485 if(icmd.ne.3) then 486! 487! calculating the matrix dU/dE 488! 489! derivative of the modified young's modulus w.r.t. 490! the Cauchy eigenvalues 491! 492 do i=1,3 493 dfla(i)=1.d0/(2.d0*fla(i)) 494 enddo 495! 496 cb(1)=1.d0/(4.d0*(we(1)-we(2))*(we(1)-we(3))) 497 ca(1)=-(wc(2)+wc(3))*cb(1) 498 cb(2)=1.d0/(4.d0*(we(2)-we(3))*(we(2)-we(1))) 499 ca(2)=-(wc(3)+wc(1))*cb(2) 500 cb(3)=1.d0/(4.d0*(we(3)-we(1))*(we(3)-we(2))) 501 ca(3)=-(wc(1)+wc(2))*cb(3) 502! 503! Cauchy tensor 504! 505 do i=1,3 506 do j=1,3 507 c(i,j)=2.d0*e(i,j) 508 enddo 509 c(i,i)=c(i,i)+1.d0 510 enddo 511! 512! matrix dU/dE 513! 514 do j=1,21 515 j1=kel(1,j) 516 j2=kel(2,j) 517 j3=kel(3,j) 518 j4=kel(4,j) 519! 520! calculating the auxiliary fields xmm1,xmm2,xmm3 521! 522 xmm1(j)=xm1(j1,j2)*xm1(j3,j4) 523 xmm2(j)=xm2(j1,j2)*xm2(j3,j4) 524 xmm3(j)=xm3(j1,j2)*xm3(j3,j4) 525! 526! calculating the auxiliary fields a and b 527! (Dhondt, G., The Finite Element Method for Three- 528! Dimensional Thermomechanical Applications, Wiley 2004, 529! formulae (4.203) and (4.204)) 530! 531 a(j)=(d(j3,j1)*d(j4,j2)+d(j3,j2)*d(j4,j1))/2.d0 532 & -xmm1(j)-xmm2(j)-xmm3(j) 533 b(j)=(d(j3,j1)*c(j4,j2)+d(j4,j1)*c(j3,j2)+ 534 & d(j4,j2)*c(j1,j3)+d(j3,j2)*c(j1,j4))/2.d0 535 & -2.d0*(wc(1)*xmm1(j)+wc(2)*xmm2(j)+wc(3)*xmm3(j)) 536! 537! calculating the matrix dU/dE 538! 539 dude(j)=2.d0*(dfla(1)*xmm1(j)+dfla(2)*xmm2(j)+ 540 & dfla(3)*xmm3(j) 541 & +fla(1)*(cb(1)*b(j)+ca(1)*a(j)) 542 & +fla(2)*(cb(2)*b(j)+ca(2)*a(j)) 543 & +fla(3)*(cb(3)*b(j)+ca(3)*a(j))) 544 enddo 545c if(iint.eq.1) write(*,*) '3 diff dude' 546c if(iint.eq.1) write(*,100) (dude(i),i=1,21) 547 endif 548 endif 549! 550 nconstants=-kode-100 551! 552 if(nconstants.eq.2) then 553! 554! isotropic 555! 556! calculating the elastic constants 557! 558 ee=elconloc(1) 559 un=elconloc(2) 560 um=ee/(1.d0+un) 561 al=un*um/(1.d0-2.d0*un) 562 um=um/2.d0 563 am1=al+2.d0*um 564 am2=2.d0*um 565! 566! determining the thermal strain 567! 568c do i=1,6 569c eth(i)=eloc(i)-emec(i) 570c enddo 571! 572! calculating the strain tensor 573! 574 do i=1,3 575 elin(i)=u(i)-1.d0 576 enddo 577 do i=4,6 578 elin(i)=u(i) 579 enddo 580! 581! determining the new total strain 582! 583c do i=1,6 584c eloc(i)=elin(i)+eth(i) 585c enddo 586! 587! calculating the stresses 588! 589 stre(1)=am1*elin(1)+al*(elin(2)+elin(3))-beta(1) 590 stre(2)=am1*elin(2)+al*(elin(1)+elin(3))-beta(2) 591 stre(3)=am1*elin(3)+al*(elin(1)+elin(2))-beta(3) 592 stre(4)=am2*elin(4)-beta(4) 593 stre(5)=am2*elin(5)-beta(5) 594 stre(6)=am2*elin(6)-beta(6) 595! 596 if(icmd.ne.3) then 597c if(iint.eq.1) then 598cc write(*,100) (stre(i),i=1,6) 599c write(*,*) 'dude' 600c write(*,100) (dude(i),i=1,6) 601c write(*,100) (dude(i),i=7,10) 602c write(*,100) (dude(i),i=11,15) 603c write(*,100) (dude(i),i=16,21) 604c endif 605! 606! calculating the stiffness matrix 607! 608! stiff_klpq=lambda*(delta_kl*dude_mmpq+delta_pq*dude_mmkl)/2 609! +2*mu*dude_klpq 610! 611 do i=1,21 612 stiff(i)=2.d0*um*dude(i) 613 enddo 614! 615 stiff( 1)=stiff( 1)+al*( 616 & dude( 1)+ 617 & dude( 2)+ 618 & dude( 4)+ 619 & dude( 1)+ 620 & dude( 2)+ 621 & dude( 4) 622 & )/2.d0 623 stiff( 2)=stiff( 2)+al*( 624 & dude( 2)+ 625 & dude( 3)+ 626 & dude( 5)+ 627 & dude( 1)+ 628 & dude( 2)+ 629 & dude( 4) 630 & )/2.d0 631 stiff( 3)=stiff( 3)+al*( 632 & dude( 2)+ 633 & dude( 3)+ 634 & dude( 5)+ 635 & dude( 2)+ 636 & dude( 3)+ 637 & dude( 5) 638 & )/2.d0 639 stiff( 4)=stiff( 4)+al*( 640 & dude( 4)+ 641 & dude( 5)+ 642 & dude( 6)+ 643 & dude( 1)+ 644 & dude( 2)+ 645 & dude( 4) 646 & )/2.d0 647 stiff( 5)=stiff( 5)+al*( 648 & dude( 4)+ 649 & dude( 5)+ 650 & dude( 6)+ 651 & dude( 2)+ 652 & dude( 3)+ 653 & dude( 5) 654 & )/2.d0 655 stiff( 6)=stiff( 6)+al*( 656 & dude( 4)+ 657 & dude( 5)+ 658 & dude( 6)+ 659 & dude( 4)+ 660 & dude( 5)+ 661 & dude( 6) 662 & )/2.d0 663 stiff( 7)=stiff( 7)+al*( 664 & dude( 7)+ 665 & dude( 8)+ 666 & dude( 9) 667 & )/2.d0 668 stiff( 8)=stiff( 8)+al*( 669 & dude( 7)+ 670 & dude( 8)+ 671 & dude( 9) 672 & )/2.d0 673 stiff( 9)=stiff( 9)+al*( 674 & dude( 7)+ 675 & dude( 8)+ 676 & dude( 9) 677 & )/2.d0 678 stiff(11)=stiff(11)+al*( 679 & dude(11)+ 680 & dude(12)+ 681 & dude(13) 682 & )/2.d0 683 stiff(12)=stiff(12)+al*( 684 & dude(11)+ 685 & dude(12)+ 686 & dude(13) 687 & )/2.d0 688 stiff(13)=stiff(13)+al*( 689 & dude(11)+ 690 & dude(12)+ 691 & dude(13) 692 & )/2.d0 693 stiff(16)=stiff(16)+al*( 694 & dude(16)+ 695 & dude(17)+ 696 & dude(18) 697 & )/2.d0 698 stiff(17)=stiff(17)+al*( 699 & dude(16)+ 700 & dude(17)+ 701 & dude(18) 702 & )/2.d0 703 stiff(18)=stiff(18)+al*( 704 & dude(16)+ 705 & dude(17)+ 706 & dude(18) 707 & )/2.d0 708c stiff(1)=stiff(1)+al*(dude(1)+dude(2)+dude(4)) 709c stiff(2)=stiff(2)+al*(dude(2)+dude(3)+dude(5)) 710c stiff(3)=stiff(3)+al*(dude(4)+dude(5)+dude(6)) 711c stiff(4)=stiff(4)+al*(dude(2)+dude(3)+dude(5)) 712c stiff(5)=stiff(5)+al*(dude(4)+dude(5)+dude(6)) 713c stiff(6)=stiff(6)+al*(dude(4)+dude(5)+dude(6)) 714c if(iint.eq.1) then 715cc write(*,100) (stre(i),i=1,6) 716c write(*,*) 'stiff' 717c write(*,100) (stiff(i),i=1,6) 718c write(*,100) (stiff(i),i=7,10) 719c write(*,100) (stiff(i),i=11,15) 720c write(*,100) (stiff(i),i=16,21) 721c endif 722 endif 723 endif 724 100 format(6(1x,e15.8)) 725! 726 return 727 end 728 729 730