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 springstiff_n2f(xl,elas,konl,voldl,s,imat,elcon,nelcon, 20 & ncmat_,ntmat_,nope,lakonl,t1l,kode,elconloc,plicon, 21 & nplicon,npmat_,iperturb,springarea,nmethod,mi,ne0, 22 & nstate_,xstateini,xstate,reltime,nasym,ielorien,orab, 23 & norien,nelem) 24! 25! calculates the stiffness of a spring (node-to-face penalty) 26! (User's manual -> theory -> boundary conditions -> 27! node-to-face penalty contact) 28! 29 implicit none 30! 31 character*8 lakonl 32! 33 integer konl(20),i,j,imat,ncmat_,ntmat_,k,l,nope,nterms,iflag, 34 & i1,kode,niso,id,nplicon(0:ntmat_,*),npmat_,nelcon(2,*), 35 & iperturb(*),nmethod,mi(*),ne0,nstate_,nasym,ielorien(mi(3),*), 36 & norien,nelem,idof,idof1,idof2,iorien 37! 38 real*8 xl(3,10),elas(21),ratio(9),pproj(3),val,shp2(7,9), 39 & al(3),s(60,60),voldl(0:mi(2),10),pl(3,10),xn(3),dm, 40 & c1,c2,c3,c4,alpha,beta,elcon(0:ncmat_,ntmat_,*),xm(3), 41 & xmu(3,3,10),dxmu(3,10),dval(3,10),fpu(3,3,10),xi,et, 42 & xs2(3,7),t1l,elconloc(21),plconloc(802),xk,fk, 43 & xiso(200),yiso(200),dd0,plicon(0:2*npmat_,ntmat_,*), 44 & a11,a12,a22,b1(3,10),b2(3,10),dal(3,3,10),qxxy(3),fnl(3), 45 & qxyy(3),dxi(3,10),det(3,10),determinant,c11,c12,c22, 46 & qxyx(3),qyxy(3),springarea(2),dd,dist,t(3),tu(3,3,10), 47 & xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*), 48 & um,eps,pi,dftdt(3,3),tp(3),te(3),ftrial(3),clear, 49 & dftrial,dfnl,dfshear,dg,dte,alnew(3),dfn(3,10),reltime, 50 & overlap,pres,dpresdoverlap,overclosure,orab(7,*), 51 & xn1(3),xn2(3),a(3,3) 52! 53! 54! 55 iflag=4 56! 57! actual positions of the nodes belonging to the contact spring 58! (otherwise no contact force) 59! 60 do i=1,nope 61 do j=1,3 62 pl(j,i)=xl(j,i)+voldl(j,i) 63 enddo 64 enddo 65! 66 if(lakonl(7:7).eq.'A') then 67 if(lakonl(4:6).eq.'RNG') then 68! 69! SPRINGA-element 70! 71 dd0=dsqrt((xl(1,2)-xl(1,1))**2 72 & +(xl(2,2)-xl(2,1))**2 73 & +(xl(3,2)-xl(3,1))**2) 74 dd=dsqrt((pl(1,2)-pl(1,1))**2 75 & +(pl(2,2)-pl(2,1))**2 76 & +(pl(3,2)-pl(3,1))**2) 77 if(dd.le.0.d0) then 78 write(*,*) '*ERROR in springstiff_n2f: spring connects' 79 write(*,*) ' two nodes at the same location:' 80 write(*,*) ' spring length is zero' 81 call exit(201) 82 endif 83 do i=1,3 84 xn(i)=(pl(i,2)-pl(i,1))/dd 85 enddo 86 val=dd-dd0 87! 88! interpolating the material data 89! 90 call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l, 91 & elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_) 92! 93! calculating the spring force and the spring constant 94! 95 if(kode.eq.2)then 96 xk=elconloc(1) 97 fk=xk*val 98 else 99 niso=int(plconloc(801)) 100 do i=1,niso 101 xiso(i)=plconloc(2*i-1) 102 yiso(i)=plconloc(2*i) 103 enddo 104 call ident(xiso,val,niso,id) 105 if(id.eq.0) then 106 xk=0.d0 107 fk=yiso(1) 108 elseif(id.eq.niso) then 109 xk=0.d0 110 fk=yiso(niso) 111 else 112 xk=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id)) 113 fk=yiso(id)+xk*(val-xiso(id)) 114 endif 115 endif 116! 117 c1=fk/dd 118 c2=xk-c1 119 do i=1,3 120 do j=1,3 121 s(i,j)=c2*xn(i)*xn(j) 122 enddo 123 s(i,i)=s(i,i)+c1 124 enddo 125 else 126! 127! GAP-element 128! interpolating the material data 129! 130 call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l, 131 & elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_) 132! 133 dd=elconloc(1) 134 xn(1)=elconloc(2) 135 xn(2)=elconloc(3) 136 xn(3)=elconloc(4) 137 xk=elconloc(5) 138 pi=4.d0*datan(1.d0) 139 eps=pi*elconloc(6)/xk 140 overclosure=-dd-xn(1)*(voldl(1,2)-voldl(1,1)) 141 & -xn(2)*(voldl(2,2)-voldl(2,1)) 142 & -xn(3)*(voldl(3,2)-voldl(3,1)) 143 fk=-xk*overclosure*(0.5d0+datan(overclosure/eps)/pi) 144 c2=xk*((0.5d0+datan(overclosure/eps)/pi) 145 & +overclosure/(pi*eps*(1.d0+(overclosure/eps)**2))) 146 do i=1,3 147 do j=1,3 148 s(i,j)=c2*xn(i)*xn(j) 149 enddo 150 enddo 151 endif 152! 153 do i=1,3 154 do j=1,3 155 s(i+3,j)=-s(i,j) 156 s(i,j+3)=-s(i,j) 157 s(i+3,j+3)=s(i,j) 158 enddo 159 enddo 160 return 161 elseif(lakonl(7:7).eq.'1') then 162! 163! spring1-element 164! determine the direction of action of the spring 165! 166 idof=nint(elcon(3,1,imat)) 167! 168 if(norien.gt.0) then 169 iorien=max(0,ielorien(1,nelem)) 170 else 171 iorien=0 172 endif 173! 174 if(iorien.eq.0) then 175 do i=1,3 176 xn(i)=0.d0 177 enddo 178 xn(idof)=1.d0 179 else 180 call transformatrix(orab(1,iorien),xl(1,1),a) 181 do i=1,3 182 xn(i)=a(i,idof) 183 enddo 184 endif 185! 186! interpolating the material data 187! 188 call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l, 189 & elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_) 190! 191! calculating the spring constant 192! 193 if(kode.eq.2)then 194 xk=elconloc(1) 195 else 196 niso=int(plconloc(801)) 197 do i=1,niso 198 xiso(i)=plconloc(2*i-1) 199 yiso(i)=plconloc(2*i) 200 enddo 201 call ident(xiso,val,niso,id) 202 if(id.eq.0) then 203 xk=0.d0 204 elseif(id.eq.niso) then 205 xk=0.d0 206 else 207 xk=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id)) 208 endif 209 endif 210! 211! assembling the stiffness matrix 212! 213 do i=1,3 214 do j=1,3 215 s(i,j)=xk*xn(i)*xn(j) 216 enddo 217 enddo 218 return 219 elseif(lakonl(7:7).eq.'2') then 220! 221! spring2-element 222! determine the direction of action of the spring 223! 224 idof1=nint(elcon(3,1,imat)) 225 idof2=nint(elcon(4,1,imat)) 226! 227 if(norien.gt.0) then 228 iorien=max(0,ielorien(1,nelem)) 229 else 230 iorien=0 231 endif 232! 233 if(iorien.eq.0) then 234 do i=1,3 235 xn1(i)=0.d0 236 xn2(i)=0.d0 237 enddo 238 xn1(idof1)=1.d0 239 xn2(idof2)=1.d0 240 else 241 call transformatrix(orab(1,iorien),xl(1,1),a) 242 do i=1,3 243 xn1(i)=a(i,idof1) 244 enddo 245 call transformatrix(orab(1,iorien),xl(1,2),a) 246 do i=1,3 247 xn2(i)=a(i,idof2) 248 enddo 249 endif 250! 251! interpolating the material data 252! 253 call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l, 254 & elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_) 255! 256! calculating the spring constant 257! 258 if(kode.eq.2)then 259 xk=elconloc(1) 260 else 261 niso=int(plconloc(801)) 262 do i=1,niso 263 xiso(i)=plconloc(2*i-1) 264 yiso(i)=plconloc(2*i) 265 enddo 266 call ident(xiso,val,niso,id) 267 if(id.eq.0) then 268 xk=0.d0 269 elseif(id.eq.niso) then 270 xk=0.d0 271 else 272 xk=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id)) 273 endif 274 endif 275! 276! assembling the stiffness matrix 277! 278 do i=1,3 279 do j=1,3 280 s(i,j)=xk*xn1(i)*xn1(j) 281 s(i,j+3)=-xk*xn1(i)*xn2(j) 282 s(i+3,j)=-xk*xn2(i)*xn1(j) 283 s(i+3,j+3)=xk*xn2(i)*xn2(j) 284 enddo 285 enddo 286 return 287! 288 endif 289! 290! contact springs 291! 292 nterms=nope-1 293! 294! vector al connects the actual position of the slave node 295! with its projection on the master face = vec_r (User's 296! manual -> theory -> boundary conditions -> node-to-face 297! penalty contact) 298! 299 do i=1,3 300 pproj(i)=pl(i,nope) 301 enddo 302 call attach_2d(pl,pproj,nterms,ratio,dist,xi,et) 303 do i=1,3 304 al(i)=pl(i,nope)-pproj(i) 305 enddo 306! 307! determining the jacobian vector on the surface 308! 309 if(nterms.eq.8) then 310 call shape8q(xi,et,pl,xm,xs2,shp2,iflag) 311 elseif(nterms.eq.4) then 312 call shape4q(xi,et,pl,xm,xs2,shp2,iflag) 313 elseif(nterms.eq.6) then 314 call shape6tri(xi,et,pl,xm,xs2,shp2,iflag) 315 else 316 call shape3tri(xi,et,pl,xm,xs2,shp2,iflag) 317 endif 318! 319! d xi / d vec_u_j 320! d eta / d vec_u_j 321! 322! dxi and det are determined from the orthogonality 323! condition 324! 325 a11=-(xs2(1,1)*xs2(1,1)+xs2(2,1)*xs2(2,1)+xs2(3,1)*xs2(3,1)) 326 & +al(1)*xs2(1,5)+al(2)*xs2(2,5)+al(3)*xs2(3,5) 327 a12=-(xs2(1,1)*xs2(1,2)+xs2(2,1)*xs2(2,2)+xs2(3,1)*xs2(3,2)) 328 & +al(1)*xs2(1,6)+al(2)*xs2(2,6)+al(3)*xs2(3,6) 329 a22=-(xs2(1,2)*xs2(1,2)+xs2(2,2)*xs2(2,2)+xs2(3,2)*xs2(3,2)) 330 & +al(1)*xs2(1,7)+al(2)*xs2(2,7)+al(3)*xs2(3,7) 331! 332 do i=1,3 333 do j=1,nterms 334 b1(i,j)=shp2(4,j)*xs2(i,1)-shp2(1,j)*al(i) 335 b2(i,j)=shp2(4,j)*xs2(i,2)-shp2(2,j)*al(i) 336 enddo 337 b1(i,nope)=-xs2(i,1) 338 b2(i,nope)=-xs2(i,2) 339 enddo 340! 341 determinant=a11*a22-a12*a12 342 c11=a22/determinant 343 c12=-a12/determinant 344 c22=a11/determinant 345! 346 do i=1,3 347 do j=1,nope 348 dxi(i,j)=c11*b1(i,j)+c12*b2(i,j) 349 det(i,j)=c12*b1(i,j)+c22*b2(i,j) 350 enddo 351 enddo 352! 353! d vec_r / d vec_u_k 354! 355 do i=1,nope 356 do j=1,3 357 do k=1,3 358 dal(j,k,i)=-xs2(j,1)*dxi(k,i)-xs2(j,2)*det(k,i) 359 enddo 360 enddo 361 enddo 362 do i=1,nterms 363 do j=1,3 364 dal(j,j,i)=dal(j,j,i)-shp2(4,i) 365 enddo 366 enddo 367 do j=1,3 368 dal(j,j,nope)=dal(j,j,nope)+1.d0 369 enddo 370! 371! d2q/dxx x dq/dy 372! 373 qxxy(1)=xs2(2,5)*xs2(3,2)-xs2(3,5)*xs2(2,2) 374 qxxy(2)=xs2(3,5)*xs2(1,2)-xs2(1,5)*xs2(3,2) 375 qxxy(3)=xs2(1,5)*xs2(2,2)-xs2(2,5)*xs2(1,2) 376! 377! dq/dx x d2q/dyy 378! 379 qxyy(1)=xs2(2,1)*xs2(3,7)-xs2(3,1)*xs2(2,7) 380 qxyy(2)=xs2(3,1)*xs2(1,7)-xs2(1,1)*xs2(3,7) 381 qxyy(3)=xs2(1,1)*xs2(2,7)-xs2(2,1)*xs2(1,7) 382! 383! Modified by Stefan Sicklinger 384! 385! dq/dx x d2q/dxy 386! 387 qxyx(1)=xs2(2,1)*xs2(3,6)-xs2(3,1)*xs2(2,6) 388 qxyx(2)=xs2(3,1)*xs2(1,6)-xs2(1,1)*xs2(3,6) 389 qxyx(3)=xs2(1,1)*xs2(2,6)-xs2(2,1)*xs2(1,6) 390! 391! d2q/dxy x dq/dy 392! 393 qyxy(1)=xs2(2,6)*xs2(3,2)-xs2(3,6)*xs2(2,2) 394 qyxy(2)=xs2(3,6)*xs2(1,2)-xs2(1,6)*xs2(3,2) 395 qyxy(3)=xs2(1,6)*xs2(2,2)-xs2(2,6)*xs2(1,2) 396! 397! 398! End modifications 399! 400! normal on the surface 401! 402 dm=dsqrt(xm(1)*xm(1)+xm(2)*xm(2)+xm(3)*xm(3)) 403 do i=1,3 404 xn(i)=xm(i)/dm 405 enddo 406! 407! distance from surface along normal (= clearance) 408! 409 clear=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3) 410 if(nmethod.eq.1) then 411 clear=clear-springarea(2)*(1.d0-reltime) 412 endif 413! 414! representative area: usually the slave surface stored in 415! springarea; however, if no area was assigned because the 416! node does not belong to any element, the master surface 417! is used 418! 419 if(springarea(1).le.0.d0) then 420 if(nterms.eq.3) then 421 springarea(1)=dm/2.d0 422 else 423 springarea(1)=dm*4.d0 424 endif 425 endif 426! 427! alpha and beta, taking the representative area into account 428! (conversion of pressure into force) 429! 430 if(int(elcon(3,1,imat)).eq.1) then 431! 432! exponential overclosure 433! 434 if(dabs(elcon(2,1,imat)).lt.1.d-30) then 435 elas(1)=0.d0 436 elas(2)=0.d0 437 else 438 alpha=elcon(2,1,imat)*springarea(1) 439 beta=elcon(1,1,imat) 440 if(-beta*clear.gt.23.d0-dlog(alpha)) then 441 beta=(dlog(alpha)-23.d0)/clear 442 endif 443 elas(1)=dexp(-beta*clear+dlog(alpha)) 444 elas(2)=-beta*elas(1) 445 endif 446 elseif(int(elcon(3,1,imat)).eq.2) then 447! 448! linear overclosure 449! 450 pi=4.d0*datan(1.d0) 451 eps=elcon(1,1,imat)*pi/elcon(2,1,imat) 452 elas(1)=(-springarea(1)*elcon(2,1,imat)*clear* 453 & (0.5d0+datan(-clear/eps)/pi)) 454 elas(2)=-springarea(1)*elcon(2,1,imat)* 455 & ((0.5d0+datan(-clear/eps)/pi)- 456 & clear/(pi*eps*(1.d0+(clear/eps)**2))) 457 elseif(int(elcon(3,1,imat)).eq.3) then 458! 459! tabular overclosure 460! 461! interpolating the material data 462! 463 call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l, 464 & elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_) 465 overlap=-clear 466 niso=int(plconloc(801)) 467 do i=1,niso 468 xiso(i)=plconloc(2*i-1) 469 yiso(i)=plconloc(2*i) 470 enddo 471 call ident(xiso,overlap,niso,id) 472 if(id.eq.0) then 473 dpresdoverlap=0.d0 474 pres=yiso(1) 475 elseif(id.eq.niso) then 476 dpresdoverlap=0.d0 477 pres=yiso(niso) 478 else 479 dpresdoverlap=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id)) 480 pres=yiso(id)+dpresdoverlap*(overlap-xiso(id)) 481 endif 482 elas(1)=springarea(1)*pres 483 elas(2)=-springarea(1)*dpresdoverlap 484 endif 485! 486! contact force 487! 488 do i=1,3 489 fnl(i)=-elas(1)*xn(i) 490 enddo 491! 492! derivatives of the jacobian vector w.r.t. the displacement 493! vectors (d m / d u_k) 494! 495 do k=1,nterms 496 xmu(1,1,k)=0.d0 497 xmu(2,2,k)=0.d0 498 xmu(3,3,k)=0.d0 499 xmu(1,2,k)=shp2(1,k)*xs2(3,2)-shp2(2,k)*xs2(3,1) 500 xmu(2,3,k)=shp2(1,k)*xs2(1,2)-shp2(2,k)*xs2(1,1) 501 xmu(3,1,k)=shp2(1,k)*xs2(2,2)-shp2(2,k)*xs2(2,1) 502 xmu(1,3,k)=-xmu(3,1,k) 503 xmu(2,1,k)=-xmu(1,2,k) 504 xmu(3,2,k)=-xmu(2,3,k) 505 enddo 506 do i=1,3 507 do j=1,3 508 xmu(i,j,nope)=0.d0 509 enddo 510 enddo 511! 512! correction due to change of xi and eta 513! 514 do k=1,nope 515 do j=1,3 516 do i=1,3 517! 518! modified by Stefan Sicklinger 519! 520 xmu(i,j,k)=xmu(i,j,k)+(qxxy(i)+qxyx(i))*dxi(j,k) 521 & +(qxyy(i)+qyxy(i))*det(j,k) 522 enddo 523 enddo 524 enddo 525! 526! derivatives of the size of the jacobian vector w.r.t. the 527! displacement vectors (d ||m||/d u_k) 528! 529 do k=1,nope 530 do i=1,3 531 dxmu(i,k)=xn(1)*xmu(1,i,k)+xn(2)*xmu(2,i,k)+ 532 & xn(3)*xmu(3,i,k) 533 enddo 534! 535! auxiliary variable: (d clear d u_k)*||m|| 536! 537 do i=1,3 538 dval(i,k)=al(1)*xmu(1,i,k)+al(2)*xmu(2,i,k)+ 539 & al(3)*xmu(3,i,k)-clear*dxmu(i,k)+ 540 & xm(1)*dal(1,i,k)+xm(2)*dal(2,i,k)+xm(3)*dal(3,i,k) 541 enddo 542! 543 enddo 544! 545 c1=1.d0/dm 546 c2=c1*c1 547 c3=elas(2)*c2 548 c4=elas(1)*c1 549! 550! derivatives of the forces w.r.t. the displacement vectors 551! 552 do k=1,nope 553 do j=1,3 554 do i=1,3 555 fpu(i,j,k)=-c3*xm(i)*dval(j,k) 556 & +c4*(xn(i)*dxmu(j,k)-xmu(i,j,k)) 557 enddo 558 enddo 559 enddo 560! 561! Coulomb friction for static calculations 562! 563 if(ncmat_.ge.7) then 564 um=elcon(6,1,imat) 565 if(um.gt.0.d0) then 566! 567! stiffness of shear stress versus slip curve 568! 569 xk=elcon(7,1,imat)*springarea(1) 570! 571! calculating the relative displacement between the slave node 572! and its projection on the master surface 573! 574 do i=1,3 575 alnew(i)=voldl(i,nope) 576 do j=1,nterms 577 alnew(i)=alnew(i)-ratio(j)*voldl(i,j) 578 enddo 579 enddo 580! 581! calculating the difference in relative displacement since 582! the start of the increment = vec_s 583! 584 do i=1,3 585 al(i)=alnew(i)-xstateini(3+i,1,ne0+konl(nope+1)) 586 enddo 587! 588! s=||vec_s|| 589! 590 val=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3) 591! 592! update the relative tangential displacement vec_t 593! 594 do i=1,3 595 t(i)=xstateini(6+i,1,ne0+konl(nope+1))+al(i)-val*xn(i) 596 enddo 597! 598! store the actual relative displacement and 599! the actual relative tangential displacement 600! 601 do i=1,3 602 xstate(3+i,1,ne0+konl(nope+1))=alnew(i) 603 xstate(6+i,1,ne0+konl(nope+1))=t(i) 604 enddo 605! 606! d vec_s/ d vec_u_k 607! notice: xi & et are const. 608! 609 do k=1,nope 610 do i=1,3 611 do j=1,3 612 dal(i,j,k)=0.d0 613 enddo 614 enddo 615 enddo 616! 617 do i=1,nterms 618 do j=1,3 619 dal(j,j,i)=-shp2(4,i) 620 enddo 621 enddo 622! 623 do j=1,3 624 dal(j,j,nope)=1.d0 625 enddo 626! 627! d s/ dvec_u_k.||m|| 628! 629 do k=1,nope 630 do i=1,3 631 dval(i,k)=al(1)*xmu(1,i,k)+al(2)*xmu(2,i,k) 632 & +al(3)*xmu(3,i,k)-val*dxmu(i,k) 633 & +xm(1)*dal(1,i,k)+xm(2)*dal(2,i,k) 634 & +xm(3)*dal(3,i,k) 635 enddo 636 enddo 637! 638! d vec_t/d vec_u_k 639! 640 do k=1,nope 641 do j=1,3 642 do i=1,3 643 tu(i,j,k)=dal(i,j,k) 644 & -c1*(xn(i)*(dval(j,k)-val*dxmu(j,k)) 645 & +val*xmu(i,j,k)) 646 enddo 647 enddo 648 enddo 649! 650! size of normal force 651! 652 dfnl=dsqrt(fnl(1)**2+fnl(2)**2+fnl(3)**2) 653! 654! maximum size of shear force 655! 656 dfshear=um*dfnl 657! 658! plastic and elastic slip 659! 660 do i=1,3 661 tp(i)=xstateini(i,1,ne0+konl(nope+1)) 662 te(i)=t(i)-tp(i) 663 enddo 664! 665! the force due to normal contact is in -xn 666! direction (internal force) -> minus signs 667! 668 do k=1,nope 669 do i=1,3 670 dfn(i,k)=-xn(1)*fpu(1,i,k)-xn(2)*fpu(2,i,k)- 671 & xn(3)*fpu(3,i,k) 672 enddo 673 enddo 674! 675 dte=dsqrt(te(1)*te(1)+te(2)*te(2)+te(3)*te(3)) 676! 677! trial force 678! 679 do i=1,3 680 ftrial(i)=xk*te(i) 681 enddo 682 dftrial=dsqrt(ftrial(1)**2+ftrial(2)**2+ftrial(3)**2) 683! 684! check whether stick or slip 685! 686 if((dftrial.lt.dfshear) .or. (dftrial.le.0.d0)) then 687! 688! stick force 689! 690 do i=1,3 691 fnl(i)=fnl(i)+ftrial(i) 692 xstate(i,1,ne0+konl(nope+1))=tp(i) 693 enddo 694! 695! stick stiffness 696! 697 do k=1,nope 698 do j=1,3 699 do i=1,3 700 fpu(i,j,k)=fpu(i,j,k)+xk*tu(i,j,k) 701 enddo 702 enddo 703 enddo 704 else 705! 706! slip force 707! 708 dg=(dftrial-dfshear)/xk 709 do i=1,3 710 ftrial(i)=te(i)/dte 711 fnl(i)=fnl(i)+dfshear*ftrial(i) 712 xstate(i,1,ne0+konl(nope+1))=tp(i)+dg*ftrial(i) 713 enddo 714! 715! slip stiffness 716! 717 c1=xk*dfshear/dftrial 718 do i=1,3 719 do j=1,3 720 dftdt(i,j)=-c1*ftrial(i)*ftrial(j) 721 enddo 722 dftdt(i,i)=dftdt(i,i)+c1 723 enddo 724! 725 do k=1,nope 726 do j=1,3 727 do i=1,3 728 do l=1,3 729 fpu(i,j,k)=fpu(i,j,k)+dftdt(i,l)*tu(l,j,k) 730 enddo 731 if((nmethod.ne.4).or.(iperturb(1).gt.1)) then 732 fpu(i,j,k)=fpu(i,j,k)+um*ftrial(i)*dfn(j,k) 733 endif 734 enddo 735 enddo 736 enddo 737 endif 738 endif 739 endif 740! 741! determining the stiffness matrix contributions 742! 743! complete field shp2 744! 745 shp2(1,nope)=0.d0 746 shp2(2,nope)=0.d0 747 shp2(4,nope)=-1.d0 748! 749 do k=1,nope 750 do l=1,nope 751 do i=1,3 752 i1=i+(k-1)*3 753 do j=1,3 754 s(i1,j+(l-1)*3)=-shp2(4,k)*fpu(i,j,l) 755 & -shp2(1,k)*fnl(i)*dxi(j,l) 756 & -shp2(2,k)*fnl(i)*det(j,l) 757 enddo 758 enddo 759 enddo 760 enddo 761! 762! symmetrizing the matrix 763! this is done in the absence of friction or for modal dynamic 764! calculations 765! 766 if((nasym.eq.0).or.((nmethod.eq.4).and.(iperturb(1).le.1))) then 767 do j=1,3*nope 768 do i=1,j-1 769 s(i,j)=(s(i,j)+s(j,i))/2.d0 770 enddo 771 enddo 772 endif 773! 774 return 775 end 776 777