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 resultsmech_se(co,kon,ipkon,lakon,ne,v, 20 & stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero, 21 & ielmat,ielorien,norien,orab,ntmat_,t0,t1,ithermal,prestr, 22 & iprestr,eme,iperturb,fn,iout,vold,nmethod, 23 & veold,dtime,time,ttime,plicon,nplicon,plkcon,nplkcon, 24 & xstateini,xstiff,xstate,npmat_,matname,mi,ielas,icmd, 25 & ncmat_,nstate_,stiini,vini,ener,eei,enerini,istep,iinc, 26 & springarea,reltime,calcul_fn,calcul_cauchy,nener, 27 & ikin,ne0,thicke,emeini,pslavsurf, 28 & pmastsurf,mortar,clearini,nea,neb,ielprop,prop,dfn, 29 & idesvar,nodedesi,fn0,sti,icoordinate,dxstiff, 30 & ialdesi,xdesi) 31! 32! calculates stresses and the material tangent at the integration 33! points and the internal forces at the nodes 34! 35 implicit none 36! 37 integer cauchy 38! 39 character*8 lakon(*),lakonl 40 character*80 amat,matname(*) 41! 42 integer kon(*),konl(26),nea,neb,mi(*),mint2d,nopes, 43 & nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*), 44 & ielorien(mi(3),*),ntmat_,ipkon(*),ne0,iflag,null, 45 & istep,iinc,mt,ne,mattyp,ithermal(*),iprestr,i,j,k,m1,m2,jj, 46 & i1,m3,m4,kk,nener,indexe,nope,norien,iperturb(*),iout, 47 & icmd,ihyper,nmethod,kode,imat,mint3d,iorien,ielas, 48 & istiff,ncmat_,nstate_,ikin,ilayer,nlayer,ki,kl,ielprop(*), 49 & nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,calcul_fn, 50 & calcul_cauchy,nopered,mortar,jfaces,igauss, 51 & idesvar,node,nodedesi(*),kscale,idir,nlgeom_undo, 52 & iactive,icoordinate,ialdesi(*),ii, 53 & node1,node2,ifaceqexp(2,20),ifacewexp(2,15) 54! 55 real*8 co(3,*),v(0:mi(2),*),shp(4,26),stiini(6,mi(1),*), 56 & stx(6,mi(1),*),xl(3,26),vl(0:mi(2),26),stre(6),prop(*), 57 & elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),xs2(3,7), 58 & alcon(0:6,ntmat_,*),vini(0:mi(2),*),thickness, 59 & alzero(*),orab(7,*),elas(21),rho,fn(0:mi(2),*), 60 & fnl(3,10),skl(3,3),beta(6),xl2(3,8),qa(4), 61 & vkl(0:3,3),t0(*),t1(*),prestr(6,mi(1),*),eme(6,mi(1),*), 62 & ckl(3,3),vold(0:mi(2),*),eloc(9),veold(0:mi(2),*), 63 & springarea(2,*),elconloc(21),eth(6),xkl(3,3),voldl(0:mi(2),26), 64 & xikl(3,3),ener(mi(1),*),emec(6),eei(6,mi(1),*),enerini(mi(1),*), 65 & emec0(6),veoldl(0:mi(2),26),xsj2(3),shp2(7,8), 66 & e,un,al,um,am1,xi,et,ze,tt,exx,eyy,ezz,exy,exz,eyz, 67 & xsj,vj,t0l,t1l,dtime,weight,pgauss(3),vij,time,ttime, 68 & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*), 69 & xstiff(27,mi(1),*),xstate(nstate_,mi(1),*),plconloc(802), 70 & vokl(3,3),xstateini(nstate_,mi(1),*),vikl(3,3), 71 & gs(8,4),a,reltime,tlayer(4),dlayer(4),xlayer(mi(3),4), 72 & thicke(mi(3),*),emeini(6,mi(1),*),clearini(3,9,*), 73 & pslavsurf(3,*),pmastsurf(6,*),sti(6,mi(1),*), 74 & fn0(0:mi(2),*),dfn(0:mi(2),*),hglf(3,4),ahr, 75 & dxstiff(27,mi(1),ne,*),xdesi(3,*) 76! 77! 78! 79 include "gauss.f" 80! 81 iflag=3 82 null=0 83 qa(3)=-1.d0 84 qa(4)=0.d0 85! 86 mt=mi(2)+1 87! 88! nodes in expansion direction of hex element 89! 90 data ifaceqexp /5,17, 91 & 6,18, 92 & 7,19, 93 & 8,20, 94 & 1,17, 95 & 2,18, 96 & 3,19, 97 & 4,20, 98 & 13,0, 99 & 14,0, 100 & 15,0, 101 & 16,0, 102 & 9,0, 103 & 10,0, 104 & 11,0, 105 & 12,0, 106 & 0,0, 107 & 0,0, 108 & 0,0, 109 & 0,0/ 110! 111! nodes in expansion direction of wedge element 112! 113 data ifacewexp /4,13, 114 & 5,14, 115 & 6,16, 116 & 1,13, 117 & 2,14, 118 & 3,15, 119 & 10,0, 120 & 11,0, 121 & 12,0, 122 & 7,0, 123 & 8,0, 124 & 9,0, 125 & 0,0, 126 & 0,0, 127 & 0,0/ 128! 129! ------------------------------------------------------------- 130! Initialisation of the loop for the variation of 131! the internal forces 132! ------------------------------------------------------------- 133! 134! 135! Loop over all elements in thread 136 do ii=nea,neb 137 if(idesvar.gt.0) then 138 i=ialdesi(ii) 139 else 140 i=ii 141 endif 142! 143 lakonl=lakon(i) 144! 145 if(ipkon(i).lt.0) cycle 146! 147! no 3D-fluid elements 148! 149 if(lakonl(1:1).eq.'F') cycle 150 if(lakonl(1:7).eq.'DCOUP3D') cycle 151! 152 if(lakonl(7:8).ne.'LC') then 153! 154 imat=ielmat(1,i) 155 amat=matname(imat) 156 if(norien.gt.0) then 157 iorien=max(0,ielorien(1,i)) 158 else 159 iorien=0 160 endif 161! 162 if(nelcon(1,imat).lt.0) then 163 ihyper=1 164 else 165 ihyper=0 166 endif 167 elseif(lakonl(4:5).eq.'20') then 168! 169! composite materials 170! 171 mint2d=4 172 nopes=8 173! determining the number of layers 174! 175 nlayer=0 176 do k=1,mi(3) 177 if(ielmat(k,i).ne.0) then 178 nlayer=nlayer+1 179 endif 180 enddo 181! 182! determining the layer thickness and global thickness 183! at the shell integration points 184! 185 iflag=1 186 indexe=ipkon(i) 187 do kk=1,mint2d 188 xi=gauss3d2(1,kk) 189 et=gauss3d2(2,kk) 190 call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag) 191 tlayer(kk)=0.d0 192 do k=1,nlayer 193 thickness=0.d0 194 do j=1,nopes 195 thickness=thickness+thicke(k,indexe+j)*shp2(4,j) 196 enddo 197 tlayer(kk)=tlayer(kk)+thickness 198 xlayer(k,kk)=thickness 199 enddo 200 enddo 201 iflag=3 202! 203 ilayer=0 204 do k=1,4 205 dlayer(k)=0.d0 206 enddo 207 elseif(lakonl(4:5).eq.'15') then 208! 209! composite materials 210! 211 mint2d=3 212 nopes=6 213! determining the number of layers 214! 215 nlayer=0 216 do k=1,mi(3) 217 if(ielmat(k,i).ne.0) then 218 nlayer=nlayer+1 219 endif 220 enddo 221! 222! determining the layer thickness and global thickness 223! at the shell integration points 224! 225 iflag=1 226 indexe=ipkon(i) 227 do kk=1,mint2d 228 xi=gauss3d10(1,kk) 229 et=gauss3d10(2,kk) 230 call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 231 tlayer(kk)=0.d0 232 do k=1,nlayer 233 thickness=0.d0 234 do j=1,nopes 235 thickness=thickness+thicke(k,indexe+j)*shp2(4,j) 236 enddo 237 tlayer(kk)=tlayer(kk)+thickness 238 xlayer(k,kk)=thickness 239 enddo 240 enddo 241 iflag=3 242! 243 ilayer=0 244 do k=1,3 245 dlayer(k)=0.d0 246 enddo 247! 248 endif 249! 250 indexe=ipkon(i) 251c Bernhardi start 252 if(lakonl(1:5).eq.'C3D8I') then 253 nope=11 254 elseif(lakonl(4:5).eq.'20') then 255c Bernhardi end 256 nope=20 257 elseif(lakonl(4:4).eq.'8') then 258 nope=8 259 elseif(lakonl(4:5).eq.'10') then 260 nope=10 261 elseif(lakonl(4:4).eq.'4') then 262 nope=4 263 elseif(lakonl(4:5).eq.'15') then 264 nope=15 265 elseif(lakonl(4:4).eq.'6') then 266 nope=6 267 elseif((lakonl(1:1).eq.'E').and.(lakonl(7:7).ne.'F')) then 268! 269! spring elements, contact spring elements and 270! dashpot elements 271! 272 if(lakonl(7:7).eq.'C') then 273! 274! contact spring elements 275! 276 if(mortar.eq.1) then 277! 278! face-to-face penalty 279! 280 nope=kon(ipkon(i)) 281 elseif(mortar.eq.0) then 282! 283! node-to-face penalty 284! 285 nope=ichar(lakonl(8:8))-47 286 konl(nope+1)=kon(indexe+nope+1) 287 endif 288 else 289! 290! genuine spring elements and dashpot elements 291! 292 nope=ichar(lakonl(8:8))-47 293 endif 294 else 295 cycle 296 endif 297! 298! computation of the coordinates of the local nodes w/ variation 299! if the designnode belongs to the considered element 300! 301 if(idesvar.gt.0) then 302! 303 if(icoordinate.eq.1) then 304! 305! the coordinates of the nodes are the design variables 306! 307 do j=1,nope 308 node=kon(indexe+j) 309 if(node.eq.nodedesi(idesvar)) then 310 iactive=j 311 exit 312 endif 313 enddo 314 endif 315! 316 endif 317! 318 if(lakonl(4:5).eq.'8R') then 319 mint3d=1 320 elseif(lakonl(4:7).eq.'20RB') then 321 if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then 322 mint3d=50 323 else 324 call beamintscheme(lakonl,mint3d,ielprop(i),prop, 325 & null,xi,et,ze,weight) 326 endif 327 elseif((lakonl(4:4).eq.'8').or. 328 & (lakonl(4:6).eq.'20R')) then 329 if(lakonl(7:8).eq.'LC') then 330 mint3d=8*nlayer 331 else 332 mint3d=8 333 endif 334 elseif(lakonl(4:4).eq.'2') then 335 mint3d=27 336 elseif(lakonl(4:5).eq.'10') then 337 mint3d=4 338 elseif(lakonl(4:4).eq.'4') then 339 mint3d=1 340 elseif(lakonl(4:5).eq.'15') then 341 if(lakonl(7:8).eq.'LC') then 342 mint3d=6*nlayer 343 else 344 mint3d=9 345 endif 346 elseif(lakonl(4:4).eq.'6') then 347 mint3d=2 348 elseif(lakonl(1:1).eq.'E') then 349 mint3d=0 350 endif 351 352 do j=1,nope 353 konl(j)=kon(indexe+j) 354 do k=1,3 355 xl(k,j)=co(k,konl(j)) 356 vl(k,j)=v(k,konl(j)) 357 voldl(k,j)=vold(k,konl(j)) 358 enddo 359 enddo 360! 361! computation of the coordinates of the local nodes w/ variation 362! if the designnode belongs to the considered element 363! 364 if((idesvar.gt.0).and.(icoordinate.eq.1)) then 365 do j=1,3 366 xl(j,iactive)=xl(j,iactive)+xdesi(j,idesvar) 367 enddo 368 if(lakonl(1:5).eq.'C3D20') then 369 if((lakonl(7:7).eq.'A').or. 370 & (lakonl(7:7).eq.'L').or. 371 & (lakonl(7:7).eq.'S').or. 372 & (lakonl(7:7).eq.'E')) then 373 node1=ifaceqexp(1,iactive) 374 node2=ifaceqexp(2,iactive) 375 do j=1,3 376 if(iactive.le.8) then 377 xl(j,node1)=xl(j,node1)+xdesi(j,idesvar) 378 xl(j,node2)=xl(j,node2)+xdesi(j,idesvar) 379 else 380 xl(j,node1)=xl(j,node1)+xdesi(j,idesvar) 381 endif 382 enddo 383 endif 384 elseif(lakonl(1:5).eq.'C3D15') then 385 if((lakonl(7:7).eq.'A').or. 386 & (lakonl(7:7).eq.'L').or. 387 & (lakonl(7:7).eq.'S').or. 388 & (lakonl(7:7).eq.'E')) then 389 node1=ifacewexp(1,iactive) 390 node2=ifacewexp(2,iactive) 391 do j=1,3 392 if(iactive.le.6) then 393 xl(j,node1)=xl(j,node1)+xdesi(j,idesvar) 394 xl(j,node2)=xl(j,node2)+xdesi(j,idesvar) 395 else 396 xl(j,node1)=xl(j,node1)+xdesi(j,idesvar) 397 endif 398 enddo 399 endif 400 elseif(lakonl(1:4).eq.'C3D8') then 401 if((lakonl(7:7).eq.'A').or. 402 & (lakonl(7:7).eq.'L').or. 403 & (lakonl(7:7).eq.'S').or. 404 & (lakonl(7:7).eq.'E')) then 405 node1=ifaceqexp(1,iactive) 406 do j=1,3 407 xl(j,node1)=xl(j,node1)+xdesi(j,idesvar) 408 enddo 409 endif 410 elseif(lakonl(1:4).eq.'C3D6') then 411 if((lakonl(7:7).eq.'A').or. 412 & (lakonl(7:7).eq.'L').or. 413 & (lakonl(7:7).eq.'S').or. 414 & (lakonl(7:7).eq.'E')) then 415 node1=ifaceqexp(1,iactive) 416 do j=1,3 417 xl(j,node1)=xl(j,node1)+xdesi(j,idesvar) 418 enddo 419 endif 420 endif 421 endif 422! 423! calculating the forces for the contact elements 424! 425 if(mint3d.eq.0) then 426! 427! "normal" spring and dashpot elements 428! 429 kode=nelcon(1,imat) 430 if(lakonl(7:7).eq.'A') then 431 t0l=0.d0 432 t1l=0.d0 433 if(ithermal(1).eq.1) then 434 t0l=(t0(konl(1))+t0(konl(2)))/2.d0 435 t1l=(t1(konl(1))+t1(konl(2)))/2.d0 436 elseif(ithermal(1).ge.2) then 437 t0l=(t0(konl(1))+t0(konl(2)))/2.d0 438 t1l=(vold(0,konl(1))+vold(0,konl(2)))/2.d0 439 endif 440 endif 441! 442! spring elements (including contact springs) 443! 444 if(lakonl(2:2).eq.'S') then 445 if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'1').or. 446 & (lakonl(7:7).eq.'2').or.((mortar.eq.0).and. 447 & ((nmethod.ne.1).or.(iperturb(1).ge.2).or.(iout.ne.-1)))) 448 & then 449 call springforc_n2f(xl,konl,vl,imat,elcon,nelcon,elas, 450 & fnl,ncmat_,ntmat_,nope,lakonl,t1l,kode,elconloc, 451 & plicon,nplicon,npmat_,ener(1,i),nener, 452 & stx(1,1,i),mi,springarea(1,konl(nope+1)),nmethod, 453 & ne0,nstate_,xstateini,xstate,reltime, 454 & ielas,ener(1,i+ne),ielorien,orab,norien,i) 455 elseif((mortar.eq.1).and. 456 & ((nmethod.ne.1).or.(iperturb(1).ge.2).or.(iout.ne.-1))) 457 & then 458 jfaces=kon(indexe+nope+2) 459 igauss=kon(indexe+nope+1) 460 call springforc_f2f(xl,vl,imat,elcon,nelcon,elas, 461 & fnl,ncmat_,ntmat_,nope,lakonl,t1l,kode,elconloc, 462 & plicon,nplicon,npmat_,ener(1,i),nener, 463 & stx(1,1,i),mi,springarea(1,igauss),nmethod, 464 & ne0,nstate_,xstateini,xstate,reltime, 465 & ielas,jfaces,igauss,pslavsurf,pmastsurf, 466 & clearini,ener(1,i+ne),kscale,konl,iout,i) 467 endif 468! 469! next lines are not executed in linstatic.c before the 470! setup of the stiffness matrix (i.e. nmethod=1 and 471! iperturb(1)<1 and iout=-1). 472! 473 if((lakonl(7:7).eq.'A').or. 474 & ((nmethod.ne.1).or.(iperturb(1).ge.2).or.(iout.ne.-1))) 475 & then 476 if(idesvar.eq.0) then 477 do j=1,nope 478 do k=1,3 479 fn0(k,indexe+j)=fn0(k,indexe+j)+fnl(k,j) 480 enddo 481 enddo 482 else 483 do j=1,nope 484 do k=1,3 485 dfn(k,konl(j))=dfn(k,konl(j))+fnl(k,j) 486 enddo 487 enddo 488 endif 489 endif 490! 491! dashpot elements (including contact dashpots) 492! 493 elseif((nmethod.eq.4).or.(nmethod.eq.5).or. 494 & ((abs(nmethod).eq.1).and.(iperturb(1).ge.2))) then 495 endif 496 elseif(ikin.eq.1) then 497 do j=1,nope 498 do k=1,3 499 veoldl(k,j)=veold(k,konl(j)) 500 enddo 501 enddo 502 endif 503! 504 do jj=1,mint3d 505 if(lakonl(4:5).eq.'8R') then 506 xi=gauss3d1(1,jj) 507 et=gauss3d1(2,jj) 508 ze=gauss3d1(3,jj) 509 weight=weight3d1(jj) 510 elseif(lakonl(4:7).eq.'20RB') then 511 if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then 512 xi=gauss3d13(1,jj) 513 et=gauss3d13(2,jj) 514 ze=gauss3d13(3,jj) 515 weight=weight3d13(jj) 516 else 517 call beamintscheme(lakonl,mint3d,ielprop(i),prop, 518 & jj,xi,et,ze,weight) 519 endif 520 elseif((lakonl(4:4).eq.'8').or. 521 & (lakonl(4:6).eq.'20R')) 522 & then 523 if(lakonl(7:8).ne.'LC') then 524 xi=gauss3d2(1,jj) 525 et=gauss3d2(2,jj) 526 ze=gauss3d2(3,jj) 527 weight=weight3d2(jj) 528 else 529 kl=mod(jj,8) 530 if(kl.eq.0) kl=8 531! 532 xi=gauss3d2(1,kl) 533 et=gauss3d2(2,kl) 534 ze=gauss3d2(3,kl) 535 weight=weight3d2(kl) 536! 537 ki=mod(jj,4) 538 if(ki.eq.0) ki=4 539! 540 if(kl.eq.1) then 541 ilayer=ilayer+1 542 if(ilayer.gt.1) then 543 do k=1,4 544 dlayer(k)=dlayer(k)+xlayer(ilayer-1,k) 545 enddo 546 endif 547 endif 548 ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/ 549 & tlayer(ki)-1.d0 550 weight=weight*xlayer(ilayer,ki)/tlayer(ki) 551! 552! material and orientation 553! 554 imat=ielmat(ilayer,i) 555 amat=matname(imat) 556 if(norien.gt.0) then 557 iorien=max(0,ielorien(ilayer,i)) 558 else 559 iorien=0 560 endif 561! 562 if(nelcon(1,imat).lt.0) then 563 ihyper=1 564 else 565 ihyper=0 566 endif 567 endif 568 elseif(lakonl(4:4).eq.'2') then 569 xi=gauss3d3(1,jj) 570 et=gauss3d3(2,jj) 571 ze=gauss3d3(3,jj) 572 weight=weight3d3(jj) 573 elseif(lakonl(4:5).eq.'10') then 574 xi=gauss3d5(1,jj) 575 et=gauss3d5(2,jj) 576 ze=gauss3d5(3,jj) 577 weight=weight3d5(jj) 578 elseif(lakonl(4:4).eq.'4') then 579 xi=gauss3d4(1,jj) 580 et=gauss3d4(2,jj) 581 ze=gauss3d4(3,jj) 582 weight=weight3d4(jj) 583 elseif(lakonl(4:5).eq.'15') then 584 if(lakonl(7:8).ne.'LC') then 585 xi=gauss3d8(1,jj) 586 et=gauss3d8(2,jj) 587 ze=gauss3d8(3,jj) 588 weight=weight3d8(jj) 589 else 590 kl=mod(jj,6) 591 if(kl.eq.0) kl=6 592! 593 xi=gauss3d10(1,kl) 594 et=gauss3d10(2,kl) 595 ze=gauss3d10(3,kl) 596 weight=weight3d10(kl) 597! 598 ki=mod(jj,3) 599 if(ki.eq.0) ki=3 600! 601 if(kl.eq.1) then 602 ilayer=ilayer+1 603 if(ilayer.gt.1) then 604 do k=1,3 605 dlayer(k)=dlayer(k)+xlayer(ilayer-1,k) 606 enddo 607 endif 608 endif 609 ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/ 610 & tlayer(ki)-1.d0 611 weight=weight*xlayer(ilayer,ki)/tlayer(ki) 612! 613! material and orientation 614! 615 imat=ielmat(ilayer,i) 616 amat=matname(imat) 617 if(norien.gt.0) then 618 iorien=max(0,ielorien(ilayer,i)) 619 else 620 iorien=0 621 endif 622! 623 if(nelcon(1,imat).lt.0) then 624 ihyper=1 625 else 626 ihyper=0 627 endif 628 endif 629 elseif(lakonl(4:4).eq.'6') then 630 xi=gauss3d7(1,jj) 631 et=gauss3d7(2,jj) 632 ze=gauss3d7(3,jj) 633 weight=weight3d7(jj) 634 endif 635! 636c Bernhardi start 637 if(lakonl(1:5).eq.'C3D8R') then 638 call shape8hr(xl,xsj,shp,gs,a) 639 elseif(lakonl(1:5).eq.'C3D8I') then 640 call shape8hu(xi,et,ze,xl,xsj,shp,iflag) 641 elseif(nope.eq.20) then 642c Bernhardi end 643 if(lakonl(7:7).eq.'A') then 644 call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag) 645 elseif((lakonl(7:7).eq.'E').or. 646 & (lakonl(7:7).eq.'S')) then 647 call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag) 648 else 649 call shape20h(xi,et,ze,xl,xsj,shp,iflag) 650 endif 651 elseif(nope.eq.8) then 652 call shape8h(xi,et,ze,xl,xsj,shp,iflag) 653 elseif(nope.eq.10) then 654 call shape10tet(xi,et,ze,xl,xsj,shp,iflag) 655 elseif(nope.eq.4) then 656 call shape4tet(xi,et,ze,xl,xsj,shp,iflag) 657 elseif(nope.eq.15) then 658 call shape15w(xi,et,ze,xl,xsj,shp,iflag) 659 else 660 call shape6w(xi,et,ze,xl,xsj,shp,iflag) 661 endif 662! 663! vkl(m2,m3) contains the derivative of the m2- 664! component of the displacement with respect to 665! direction m3 666! 667 do m2=1,3 668 do m3=1,3 669 vkl(m2,m3)=0.d0 670 enddo 671 enddo 672! 673 do m1=1,nope 674 do m2=1,3 675 do m3=1,3 676 vkl(m2,m3)=vkl(m2,m3)+shp(m3,m1)*vl(m2,m1) 677 enddo 678c write(*,*) 'vnoeie',i,konl(m1),(vkl(m2,k),k=1,3) 679 enddo 680 enddo 681! 682! for frequency analysis or buckling with preload the 683! strains are calculated with respect to the deformed 684! configuration 685! for a linear iteration within a nonlinear increment: 686! the tangent matrix is calculated at strain at the end 687! of the previous increment 688! 689 if((iperturb(1).eq.1).or.(iperturb(1).eq.-1))then 690 do m2=1,3 691 do m3=1,3 692 vokl(m2,m3)=0.d0 693 enddo 694 enddo 695! 696 do m1=1,nope 697 do m2=1,3 698 do m3=1,3 699 vokl(m2,m3)=vokl(m2,m3)+ 700 & shp(m3,m1)*voldl(m2,m1) 701 enddo 702 enddo 703 enddo 704 endif 705! 706 kode=nelcon(1,imat) 707! 708! calculating the strain 709! 710! attention! exy,exz and eyz are engineering strains! 711! 712 exx=vkl(1,1) 713 eyy=vkl(2,2) 714 ezz=vkl(3,3) 715 exy=vkl(1,2)+vkl(2,1) 716 exz=vkl(1,3)+vkl(3,1) 717 eyz=vkl(2,3)+vkl(3,2) 718! 719 if(iperturb(2).eq.1) then 720! 721! Lagrangian strain 722! 723 exx=exx+(vkl(1,1)**2+vkl(2,1)**2+vkl(3,1)**2)/2.d0 724 eyy=eyy+(vkl(1,2)**2+vkl(2,2)**2+vkl(3,2)**2)/2.d0 725 ezz=ezz+(vkl(1,3)**2+vkl(2,3)**2+vkl(3,3)**2)/2.d0 726 exy=exy+vkl(1,1)*vkl(1,2)+vkl(2,1)*vkl(2,2)+ 727 & vkl(3,1)*vkl(3,2) 728 exz=exz+vkl(1,1)*vkl(1,3)+vkl(2,1)*vkl(2,3)+ 729 & vkl(3,1)*vkl(3,3) 730 eyz=eyz+vkl(1,2)*vkl(1,3)+vkl(2,2)*vkl(2,3)+ 731 & vkl(3,2)*vkl(3,3) 732! 733! for frequency analysis or buckling with preload the 734! strains are calculated with respect to the deformed 735! configuration 736! 737 elseif(iperturb(1).eq.1) then 738 exx=exx+vokl(1,1)*vkl(1,1)+vokl(2,1)*vkl(2,1)+ 739 & vokl(3,1)*vkl(3,1) 740 eyy=eyy+vokl(1,2)*vkl(1,2)+vokl(2,2)*vkl(2,2)+ 741 & vokl(3,2)*vkl(3,2) 742 ezz=ezz+vokl(1,3)*vkl(1,3)+vokl(2,3)*vkl(2,3)+ 743 & vokl(3,3)*vkl(3,3) 744 exy=exy+vokl(1,1)*vkl(1,2)+vokl(1,2)*vkl(1,1)+ 745 & vokl(2,1)*vkl(2,2)+vokl(2,2)*vkl(2,1)+ 746 & vokl(3,1)*vkl(3,2)+vokl(3,2)*vkl(3,1) 747 exz=exz+vokl(1,1)*vkl(1,3)+vokl(1,3)*vkl(1,1)+ 748 & vokl(2,1)*vkl(2,3)+vokl(2,3)*vkl(2,1)+ 749 & vokl(3,1)*vkl(3,3)+vokl(3,3)*vkl(3,1) 750 eyz=eyz+vokl(1,2)*vkl(1,3)+vokl(1,3)*vkl(1,2)+ 751 & vokl(2,2)*vkl(2,3)+vokl(2,3)*vkl(2,2)+ 752 & vokl(3,2)*vkl(3,3)+vokl(3,3)*vkl(3,2) 753 endif 754! 755! storing the local strains 756! 757 if(iperturb(1).ne.-1) then 758 eloc(1)=exx 759 eloc(2)=eyy 760 eloc(3)=ezz 761 eloc(4)=exy/2.d0 762 eloc(5)=exz/2.d0 763 eloc(6)=eyz/2.d0 764 else 765! 766! linear iteration within a nonlinear increment: 767! 768 eloc(1)=vokl(1,1)+ 769 & (vokl(1,1)**2+vokl(2,1)**2+vokl(3,1)**2)/2.d0 770 eloc(2)=vokl(2,2)+ 771 & (vokl(1,2)**2+vokl(2,2)**2+vokl(3,2)**2)/2.d0 772 eloc(3)=vokl(3,3)+ 773 & (vokl(1,3)**2+vokl(2,3)**2+vokl(3,3)**2)/2.d0 774 eloc(4)=(vokl(1,2)+vokl(2,1)+vokl(1,1)*vokl(1,2)+ 775 & vokl(2,1)*vokl(2,2)+vokl(3,1)*vokl(3,2))/2.d0 776 eloc(5)=(vokl(1,3)+vokl(3,1)+vokl(1,1)*vokl(1,3)+ 777 & vokl(2,1)*vokl(2,3)+vokl(3,1)*vokl(3,3))/2.d0 778 eloc(6)=(vokl(2,3)+vokl(3,2)+vokl(1,2)*vokl(1,3)+ 779 & vokl(2,2)*vokl(2,3)+vokl(3,2)*vokl(3,3))/2.d0 780 endif 781! 782! calculating the deformation gradient (needed to 783! convert the element stiffness matrix from spatial 784! coordinates to material coordinates 785! deformation plasticity) 786! 787 if((kode.eq.-50).or.(kode.le.-100)) then 788! 789! calculating the deformation gradient 790! 791c Bernhardi start 792 xkl(1,1)=vkl(1,1)+1.0d0 793 xkl(2,2)=vkl(2,2)+1.0d0 794 xkl(3,3)=vkl(3,3)+1.0d0 795c Bernhardi end 796 xkl(1,2)=vkl(1,2) 797 xkl(1,3)=vkl(1,3) 798 xkl(2,3)=vkl(2,3) 799 xkl(2,1)=vkl(2,1) 800 xkl(3,1)=vkl(3,1) 801 xkl(3,2)=vkl(3,2) 802! 803! calculating the Jacobian 804! 805 vj=xkl(1,1)*(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2)) 806 & -xkl(1,2)*(xkl(2,1)*xkl(3,3)-xkl(2,3)*xkl(3,1)) 807 & +xkl(1,3)*(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1)) 808! 809! inversion of the deformation gradient (only for 810! deformation plasticity) 811! 812 if(kode.eq.-50) then 813! 814 ckl(1,1)=(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2))/vj 815 ckl(2,2)=(xkl(1,1)*xkl(3,3)-xkl(1,3)*xkl(3,1))/vj 816 ckl(3,3)=(xkl(1,1)*xkl(2,2)-xkl(1,2)*xkl(2,1))/vj 817 ckl(1,2)=(xkl(1,3)*xkl(3,2)-xkl(1,2)*xkl(3,3))/vj 818 ckl(1,3)=(xkl(1,2)*xkl(2,3)-xkl(2,2)*xkl(1,3))/vj 819 ckl(2,3)=(xkl(2,1)*xkl(1,3)-xkl(1,1)*xkl(2,3))/vj 820 ckl(2,1)=(xkl(3,1)*xkl(2,3)-xkl(2,1)*xkl(3,3))/vj 821 ckl(3,1)=(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1))/vj 822 ckl(3,2)=(xkl(3,1)*xkl(1,2)-xkl(1,1)*xkl(3,2))/vj 823! 824! converting the Lagrangian strain into Eulerian 825! strain (only for deformation plasticity) 826! 827 cauchy=0 828 call str2mat(eloc,ckl,vj,cauchy) 829 endif 830! 831 endif 832! 833! calculating fields for incremental plasticity 834! 835 if(kode.le.-100) then 836! 837! calculating the deformation gradient at the 838! start of the increment 839! 840! calculating the displacement gradient at the 841! start of the increment 842! 843 do m2=1,3 844 do m3=1,3 845 vikl(m2,m3)=0.d0 846 enddo 847 enddo 848! 849 do m1=1,nope 850 do m2=1,3 851 do m3=1,3 852 vikl(m2,m3)=vikl(m2,m3) 853 & +shp(m3,m1)*vini(m2,konl(m1)) 854 enddo 855 enddo 856 enddo 857! 858! calculating the deformation gradient of the old 859! fields 860! 861 xikl(1,1)=vikl(1,1)+1.d0 862 xikl(2,2)=vikl(2,2)+1.d0 863 xikl(3,3)=vikl(3,3)+1.d0 864 xikl(1,2)=vikl(1,2) 865 xikl(1,3)=vikl(1,3) 866 xikl(2,3)=vikl(2,3) 867 xikl(2,1)=vikl(2,1) 868 xikl(3,1)=vikl(3,1) 869 xikl(3,2)=vikl(3,2) 870! 871! calculating the Jacobian 872! 873 vij=xikl(1,1)*(xikl(2,2)*xikl(3,3) 874 & -xikl(2,3)*xikl(3,2)) 875 & -xikl(1,2)*(xikl(2,1)*xikl(3,3) 876 & -xikl(2,3)*xikl(3,1)) 877 & +xikl(1,3)*(xikl(2,1)*xikl(3,2) 878 & -xikl(2,2)*xikl(3,1)) 879! 880! stresses at the start of the increment 881! 882 do m1=1,6 883 stre(m1)=sti(m1,jj,i) 884 enddo 885! 886 endif 887! 888! prestress values 889! 890 if(iprestr.eq.1) then 891 do kk=1,6 892 beta(kk)=-prestr(kk,jj,i) 893 enddo 894 else 895 do kk=1,6 896 beta(kk)=0.d0 897 enddo 898 endif 899! 900 if(ithermal(1).ge.1) then 901! 902! calculating the temperature difference in 903! the integration point 904! 905 t0l=0.d0 906 t1l=0.d0 907 if(ithermal(1).eq.1) then 908 if((lakonl(4:5).eq.'8 ').or. 909 & (lakonl(4:5).eq.'8I')) then 910 do i1=1,8 911 t0l=t0l+t0(konl(i1))/8.d0 912 t1l=t1l+t1(konl(i1))/8.d0 913 enddo 914 elseif(lakonl(4:6).eq.'20 ') then 915 nopered=20 916 call lintemp(t0,konl,nopered,jj,t0l) 917 call lintemp(t1,konl,nopered,jj,t1l) 918 elseif(lakonl(4:6).eq.'10T') then 919 call linscal10(t0,konl,t0l,null,shp) 920 call linscal10(t1,konl,t1l,null,shp) 921 else 922 do i1=1,nope 923 t0l=t0l+shp(4,i1)*t0(konl(i1)) 924 t1l=t1l+shp(4,i1)*t1(konl(i1)) 925 enddo 926 endif 927 elseif(ithermal(1).ge.2) then 928 if((lakonl(4:5).eq.'8 ').or. 929 & (lakonl(4:5).eq.'8I')) then 930 do i1=1,8 931 t0l=t0l+t0(konl(i1))/8.d0 932 t1l=t1l+vold(0,konl(i1))/8.d0 933 enddo 934 elseif(lakonl(4:6).eq.'20 ') then 935 nopered=20 936 call lintemp_th0(t0,konl,nopered,jj,t0l,mi) 937 call lintemp_th1(vold,konl,nopered,jj,t1l,mi) 938 elseif(lakonl(4:6).eq.'10T') then 939 call linscal10(t0,konl,t0l,null,shp) 940 call linscal10(vold,konl,t1l,mi(2),shp) 941 else 942 do i1=1,nope 943 t0l=t0l+shp(4,i1)*t0(konl(i1)) 944 t1l=t1l+shp(4,i1)*vold(0,konl(i1)) 945 enddo 946 endif 947 endif 948 tt=t1l-t0l 949 endif 950! 951! calculating the coordinates of the integration point 952! for material orientation purposes (for cylindrical 953! coordinate systems) 954! 955 if((iorien.gt.0).or.(kode.le.-100)) then 956 do j=1,3 957 pgauss(j)=0.d0 958 do i1=1,nope 959 pgauss(j)=pgauss(j)+shp(4,i1)*co(j,konl(i1)) 960 enddo 961 enddo 962 endif 963! 964! material data; for linear elastic materials 965! this includes the calculation of the stiffness 966! matrix 967! 968 istiff=0 969! 970 call materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon, 971 & nalcon,imat,amat,iorien,pgauss,orab,ntmat_, 972 & elas,rho,i,ithermal,alzero,mattyp,t0l,t1l,ihyper, 973 & istiff,elconloc,eth,kode,plicon,nplicon, 974 & plkcon,nplkcon,npmat_,plconloc,mi(1),dtime,jj, 975 & xstiff,ncmat_) 976! 977! determining the mechanical strain 978! 979 if(ithermal(1).ne.0) then 980 do m1=1,6 981 emec(m1)=eloc(m1)-eth(m1) 982c emec0(m1)=emeini(m1,jj,i) 983 enddo 984 else 985 do m1=1,6 986 emec(m1)=eloc(m1) 987c emec0(m1)=emeini(m1,jj,i) 988 enddo 989 endif 990 if(kode.le.-100) then 991 do m1=1,6 992 emec0(m1)=emeini(m1,jj,i) 993 enddo 994 endif 995! 996! subtracting the plastic initial strains 997! 998 if(iprestr.eq.2) then 999 do m1=1,6 1000 emec(m1)=emec(m1)-prestr(m1,jj,i) 1001 enddo 1002 endif 1003! 1004! calculating the local stiffness and stress 1005! 1006 nlgeom_undo=0 1007 call mechmodel(elconloc,elas,emec,kode,emec0,ithermal, 1008 & icmd,beta,stre,xkl,ckl,vj,xikl,vij, 1009 & plconloc,xstate,xstateini,ielas, 1010 & amat,t1l,dtime,time,ttime,i,jj,nstate_,mi(1), 1011 & iorien,pgauss,orab,eloc,mattyp,qa(3),istep,iinc, 1012 & ipkon,nmethod,iperturb,qa(4),nlgeom_undo) 1013! 1014 if(((nmethod.ne.4).or.(iperturb(1).ne.0)).and. 1015 & (nmethod.ne.5)) then 1016 if(idesvar.eq.0) then 1017 do m1=1,21 1018 xstiff(m1,jj,i)=elas(m1) 1019 enddo 1020 elseif(icoordinate.ne.1) then 1021! 1022! for orientation design variables: store the modified 1023! stiffness matrix for use in mafillsmse 1024! 1025 idir=idesvar-3*((idesvar-1)/3) 1026 do m1=1,21 1027 dxstiff(m1,jj,i,idir)=elas(m1) 1028 enddo 1029 endif 1030 endif 1031! 1032 if((iperturb(1).eq.-1).and.(nlgeom_undo.eq.0)) then 1033! 1034! if the forced displacements were changed at 1035! the start of a nonlinear step, the nodal 1036! forces due do this displacements are 1037! calculated in a purely linear way, and 1038! the first iteration is purely linear in order 1039! to allow the displacements to redistribute 1040! in a quasi-static way (only applies to 1041! quasi-static analyses (*STATIC)) 1042! 1043 eloc(1)=exx-vokl(1,1) 1044 eloc(2)=eyy-vokl(2,2) 1045 eloc(3)=ezz-vokl(3,3) 1046 eloc(4)=exy-(vokl(1,2)+vokl(2,1)) 1047 eloc(5)=exz-(vokl(1,3)+vokl(3,1)) 1048 eloc(6)=eyz-(vokl(2,3)+vokl(3,2)) 1049! 1050 if(mattyp.eq.1) then 1051 e=elas(1) 1052 un=elas(2) 1053 um=e/(1.d0+un) 1054 al=un*um/(1.d0-2.d0*un) 1055 um=um/2.d0 1056 am1=al*(eloc(1)+eloc(2)+eloc(3)) 1057 stre(1)=am1+2.d0*um*eloc(1) 1058 stre(2)=am1+2.d0*um*eloc(2) 1059 stre(3)=am1+2.d0*um*eloc(3) 1060 stre(4)=um*eloc(4) 1061 stre(5)=um*eloc(5) 1062 stre(6)=um*eloc(6) 1063 elseif(mattyp.eq.2) then 1064 stre(1)=eloc(1)*elas(1)+eloc(2)*elas(2) 1065 & +eloc(3)*elas(4) 1066 stre(2)=eloc(1)*elas(2)+eloc(2)*elas(3) 1067 & +eloc(3)*elas(5) 1068 stre(3)=eloc(1)*elas(4)+eloc(2)*elas(5) 1069 & +eloc(3)*elas(6) 1070 stre(4)=eloc(4)*elas(7) 1071 stre(5)=eloc(5)*elas(8) 1072 stre(6)=eloc(6)*elas(9) 1073 elseif(mattyp.eq.3) then 1074 stre(1)=eloc(1)*elas(1)+eloc(2)*elas(2)+ 1075 & eloc(3)*elas(4)+eloc(4)*elas(7)+ 1076 & eloc(5)*elas(11)+eloc(6)*elas(16) 1077 stre(2)=eloc(1)*elas(2)+eloc(2)*elas(3)+ 1078 & eloc(3)*elas(5)+eloc(4)*elas(8)+ 1079 & eloc(5)*elas(12)+eloc(6)*elas(17) 1080 stre(3)=eloc(1)*elas(4)+eloc(2)*elas(5)+ 1081 & eloc(3)*elas(6)+eloc(4)*elas(9)+ 1082 & eloc(5)*elas(13)+eloc(6)*elas(18) 1083 stre(4)=eloc(1)*elas(7)+eloc(2)*elas(8)+ 1084 & eloc(3)*elas(9)+eloc(4)*elas(10)+ 1085 & eloc(5)*elas(14)+eloc(6)*elas(19) 1086 stre(5)=eloc(1)*elas(11)+eloc(2)*elas(12)+ 1087 & eloc(3)*elas(13)+eloc(4)*elas(14)+ 1088 & eloc(5)*elas(15)+eloc(6)*elas(20) 1089 stre(6)=eloc(1)*elas(16)+eloc(2)*elas(17)+ 1090 & eloc(3)*elas(18)+eloc(4)*elas(19)+ 1091 & eloc(5)*elas(20)+eloc(6)*elas(21) 1092 endif 1093 endif 1094! 1095 skl(1,1)=stre(1) 1096 skl(2,2)=stre(2) 1097 skl(3,3)=stre(3) 1098 skl(2,1)=stre(4) 1099 skl(3,1)=stre(5) 1100 skl(3,2)=stre(6) 1101! 1102 skl(1,2)=skl(2,1) 1103 skl(1,3)=skl(3,1) 1104 skl(2,3)=skl(3,2) 1105! 1106! calculation of the nodal forces 1107! 1108 if(calcul_fn.eq.1)then 1109 if(idesvar.eq.0) then 1110! 1111! unperturbed state 1112! 1113 do m1=1,nope 1114 do m2=1,3 1115! 1116! linear elastic part 1117! 1118 do m3=1,3 1119 fn0(m2,indexe+m1)=fn0(m2,indexe+m1)+ 1120 & xsj*skl(m2,m3)*shp(m3,m1)*weight 1121 enddo 1122! 1123! nonlinear geometric part 1124! 1125 if((iperturb(2).eq.1).and.(nlgeom_undo.eq.0)) 1126 & then 1127 do m3=1,3 1128 do m4=1,3 1129 fn0(m2,indexe+m1)=fn0(m2,indexe+m1)+ 1130 & xsj*skl(m4,m3)*weight* 1131 & (vkl(m2,m4)*shp(m3,m1)+ 1132 & vkl(m2,m3)*shp(m4,m1))/2.d0 1133 enddo 1134 enddo 1135 endif 1136! 1137 enddo 1138 enddo 1139c Bernhardi start 1140 if(lakonl(1:5).eq.'C3D8R') then 1141 ahr=elas(1)*a 1142! 1143 do i1=1,3 1144 do k=1,4 1145 hglf(i1,k)=0.0d0 1146 do j=1,8 1147 hglf(i1,k)=hglf(i1,k)+gs(j,k)*vl(i1,j) 1148 enddo 1149 hglf(i1,k)=hglf(i1,k)*ahr 1150 enddo 1151 enddo 1152 do i1=1,3 1153 do j=1,8 1154 do k=1,4 1155 fn0(i1,indexe+j)=fn0(i1,indexe+j)+ 1156 & hglf(i1,k)*gs(j,k) 1157 enddo 1158 enddo 1159 enddo 1160 endif 1161c Bernhardi end 1162 else 1163! 1164! perturbed state 1165! 1166 do m1=1,nope 1167 do m2=1,3 1168! 1169! linear elastic part 1170! 1171 do m3=1,3 1172 dfn(m2,konl(m1))=dfn(m2,konl(m1))+ 1173 & xsj*skl(m2,m3)*shp(m3,m1)*weight 1174 enddo 1175! 1176! nonlinear geometric part 1177! 1178 if((iperturb(2).eq.1).and.(nlgeom_undo.eq.0)) 1179 & then 1180 do m3=1,3 1181 do m4=1,3 1182 dfn(m2,konl(m1))=dfn(m2,konl(m1))+ 1183 & xsj*skl(m4,m3)*weight* 1184 & (vkl(m2,m4)*shp(m3,m1)+ 1185 & vkl(m2,m3)*shp(m4,m1))/2.d0 1186 enddo 1187 enddo 1188 endif 1189! 1190 enddo 1191 enddo 1192c Bernhardi start 1193 if(lakonl(1:5).eq.'C3D8R') then 1194! 1195 ahr=elas(1)*a 1196! 1197 do i1=1,3 1198 do k=1,4 1199 hglf(i1,k)=0.0d0 1200 do j=1,8 1201 hglf(i1,k)=hglf(i1,k)+gs(j,k)*vl(i1,j) 1202 enddo 1203 hglf(i1,k)=hglf(i1,k)*ahr 1204 enddo 1205 enddo 1206 do i1=1,3 1207 do j=1,8 1208 do k=1,4 1209 dfn(i1,konl(j))=dfn(i1,konl(j)) 1210 & +hglf(i1,k)*gs(j,k) 1211 enddo 1212 enddo 1213 enddo 1214 endif 1215c Bernhardi end 1216 endif 1217 endif 1218! 1219 enddo ! enddo loop over the integration points 1220! 1221! subtracting the unperturbed state of the element at stake 1222! 1223 if(idesvar.gt.0) then 1224 do m1=1,nope 1225 do m2=1,3 1226 dfn(m2,konl(m1))=dfn(m2,konl(m1)) 1227 & -fn0(m2,indexe+m1) 1228 enddo 1229 enddo 1230 endif 1231! 1232! end of loop over all elements in thread 1233 enddo 1234! 1235 return 1236 end 1237