1! CalculiX - A 3-dimensional finite element program 2! Copyright (C) 1998-2021 Guido Dhondt 3! 4! This program is free software; you can redistribute it and/or 5! modify it under the terms of the GNU General Public License as 6! published by the Free Software Foundation(version 2); 7! 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12! GNU General Public License for more details. 13! 14! You should have received a copy of the GNU General Public License 15! along with this program; if not, write to the Free Software 16! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 17! 18 subroutine us3_materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon, 19 & nalcon,imat,amat,iorien,pgauss,orab,ntmat_,elas,rho,iel, 20 & ithermal,alzero,mattyp,t0l,t1l,ihyper,istiff,elconloc,eth, 21 & kode,plicon,nplicon,plkcon,nplkcon,npmat_,plconloc,mi,dtime, 22 & iint,xstiff,ncmat_) 23! 24 implicit none 25! 26! determines the material data for element iel 27! 28! istiff=0: only interpolation of material data 29! istiff=1: copy the consistent tangent matrix from the field 30! xstiff and check for zero entries 31! 32 character*80 amat 33! 34 integer nelcon(2,*),nrhcon(*),nalcon(2,*), 35 & imat,iorien,ithermal(*),j,k,mattyp,kal(2,6),j1,j2,j3,j4, 36 & jj,ntmat_,istiff,nelconst,ihyper,kode,itemp,kin,nelas, 37 & iel,iint,mi(*),ncmat_,id,two,seven, 38 & nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_ 39! 40 real*8 elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*), 41 & alcon(0:6,ntmat_,*),eth(6),xstiff(27,mi(1),*), 42 & orab(7,*),elas(21),alph(6),alzero(*),rho,t0l,t1l, 43 & skl(3,3),xa(3,3),elconloc(*),emax,pgauss(3), 44 & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*), 45 & plconloc(802),dtime 46! 47! 48! 49 kal=reshape((/1,1,2,2,3,3,1,2,1,3,2,3/),(/2,6/)) 50! 51 two=2 52 seven=7 53! 54! nelconst: # constants read from file 55! nelas: # constants in the local tangent stiffness matrix 56! 57 if(istiff.eq.1) then 58 59 nelas=nelcon(1,imat) 60 if((nelas.lt.0).or.((nelas.ne.2).and.(iorien.ne.0))) nelas=21 61! 62! calculating the density (needed for the mass matrix and 63! gravity or centrifugal loading) 64! 65 if(ithermal(1).eq.0) then 66 rho=rhcon(1,1,imat) 67 else 68 call ident2(rhcon(0,1,imat),t1l,nrhcon(imat),two,id) 69 if(nrhcon(imat).eq.0) then 70 continue 71 elseif(nrhcon(imat).eq.1) then 72 rho=rhcon(1,1,imat) 73 elseif(id.eq.0) then 74 rho=rhcon(1,1,imat) 75 elseif(id.eq.nrhcon(imat)) then 76 rho=rhcon(1,id,imat) 77 else 78 rho=rhcon(1,id,imat)+ 79 & (rhcon(1,id+1,imat)-rhcon(1,id,imat))* 80 & (t1l-rhcon(0,id,imat))/ 81 & (rhcon(0,id+1,imat)-rhcon(0,id,imat)) 82 endif 83 endif 84! 85! for nonlinear behavior (nonlinear geometric or 86! nonlinear material behavior): copy the stiffness matrix 87! from the last stress calculation 88! 89 do j=1,21 90 elas(j)=xstiff(j,iint,iel) 91 enddo 92! 93! check whether the fully anisotropic case can be 94! considered as orthotropic 95! 96 if(nelas.eq.21) then 97 emax=0.d0 98 do j=1,9 99 emax=max(emax,dabs(elas(j))) 100 enddo 101 do j=10,21 102 if(dabs(elas(j)).gt.emax*1.d-10) then 103 emax=-1.d0 104 exit 105 endif 106 enddo 107 if(emax.gt.0.d0) nelas=9 108 endif 109! 110! determining the type: isotropic, orthotropic or anisotropic 111! 112 if(nelas.le.2) then 113 mattyp=1 114 elseif(nelas.le.9) then 115 mattyp=2 116 else 117 mattyp=3 118 endif 119! 120 else 121! 122 nelconst=nelcon(1,imat) 123! 124 if(nelconst.lt.0) then 125! 126! inelastic material or user material 127! 128 if(nelconst.eq.-1) then 129 nelconst=3 130 elseif(nelconst.eq.-2) then 131 nelconst=3 132 elseif(nelconst.eq.-3) then 133 nelconst=2 134 elseif(nelconst.eq.-4) then 135 nelconst=3 136 elseif(nelconst.eq.-5) then 137 nelconst=6 138 elseif(nelconst.eq.-6) then 139 nelconst=9 140 elseif(nelconst.eq.-7) then 141 nelconst=3 142 elseif(nelconst.eq.-8) then 143 nelconst=7 144 elseif(nelconst.eq.-9) then 145 nelconst=12 146 elseif(nelconst.eq.-10) then 147 nelconst=2 148 elseif(nelconst.eq.-11) then 149 nelconst=4 150 elseif(nelconst.eq.-12) then 151 nelconst=6 152 elseif(nelconst.eq.-13) then 153 nelconst=5 154 elseif(nelconst.eq.-14) then 155 nelconst=6 156 elseif(nelconst.eq.-15) then 157 nelconst=3 158 elseif(nelconst.eq.-16) then 159 nelconst=6 160 elseif(nelconst.eq.-17) then 161 nelconst=9 162 elseif(nelconst.eq.-50) then 163 nelconst=5 164 elseif(nelconst.eq.-51) then 165 nelconst=2 166 elseif(nelconst.eq.-52) then 167 nelconst=5 168 elseif(nelconst.le.-100) then 169 nelconst=-nelconst-100 170 endif 171! 172 endif 173! 174! in case no initial temperatures are defined, the calculation 175! is assumed athermal, and the first available set material 176! constants are used 177! 178 if(ithermal(1).eq.0) then 179 if(ihyper.ne.1) then 180 do k=1,nelconst 181 elconloc(k)=elcon(k,1,imat) 182 enddo 183 else 184 do k=1,nelconst 185 elconloc(k)=elcon(k,1,imat) 186 enddo 187! 188 itemp=1 189! 190 if((kode.lt.-50).and.(kode.gt.-100)) then 191 plconloc(1)=0.d0 192 plconloc(2)=0.d0 193 plconloc(3)=0.d0 194 plconloc(801)=nplicon(1,imat)+0.5d0 195 plconloc(802)=nplkcon(1,imat)+0.5d0 196! 197! isotropic hardening 198! 199 if(nplicon(1,imat).ne.0) then 200 kin=0 201 call plcopy(plicon,nplicon,plconloc,npmat_,ntmat_, 202 & imat,itemp,iel,kin) 203 endif 204! 205! kinematic hardening 206! 207 if(nplkcon(1,imat).ne.0) then 208 kin=1 209 call plcopy(plkcon,nplkcon,plconloc,npmat_,ntmat_, 210 & imat,itemp,iel,kin) 211 endif 212! 213 endif 214! 215 endif 216 else 217! 218! calculating the expansion coefficients 219! 220! 221 222 call ident2(alcon(0,1,imat),t1l,nalcon(2,imat),seven,id) 223 if(nalcon(2,imat).eq.0) then 224 do k=1,6 225 alph(k)=0.d0 226 enddo 227 continue 228 elseif(nalcon(2,imat).eq.1) then 229 do k=1,nalcon(1,imat) 230 alph(k)=alcon(k,1,imat)!*(t1l-alzero(imat)) 231 enddo 232 elseif(id.eq.0) then 233 do k=1,nalcon(1,imat) 234 alph(k)=alcon(k,1,imat)!*(t1l-alzero(imat)) 235 enddo 236 elseif(id.eq.nalcon(2,imat)) then 237 do k=1,nalcon(1,imat) 238 alph(k)=alcon(k,id,imat)!*(t1l-alzero(imat)) 239 enddo 240 else 241 do k=1,nalcon(1,imat) 242 alph(k)=(alcon(k,id,imat)+ 243 & (alcon(k,id+1,imat)-alcon(k,id,imat))* 244 & (t1l-alcon(0,id,imat))/ 245 & (alcon(0,id+1,imat)-alcon(0,id,imat))) 246! & *(t1l-alzero(imat)) 247 enddo 248 endif 249! 250! subtracting the initial temperature influence 251! 252! call ident2(alcon(0,1,imat),t0l,nalcon(2,imat),seven,id) 253! if(nalcon(2,imat).eq.0) then 254! continue 255! elseif(nalcon(2,imat).eq.1) then 256! do k=1,nalcon(1,imat) 257! alph(k)=alph(k)-alcon(k,1,imat)*(t0l-alzero(imat)) 258! enddo 259! elseif(id.eq.0) then 260! do k=1,nalcon(1,imat) 261! alph(k)=alph(k)-alcon(k,1,imat)*(t0l-alzero(imat)) 262! enddo 263! elseif(id.eq.nalcon(2,imat)) then 264! do k=1,nalcon(1,imat) 265! alph(k)=alph(k)-alcon(k,id,imat)*(t0l-alzero(imat)) 266! enddo 267! else 268! do k=1,nalcon(1,imat) 269! alph(k)=alph(k)-(alcon(k,id,imat)+ 270! & (alcon(k,id+1,imat)-alcon(k,id,imat))* 271! & (t0l-alcon(0,id,imat))/ 272! & (alcon(0,id+1,imat)-alcon(0,id,imat))) 273! & *(t0l-alzero(imat)) 274! enddo 275! endif 276! 277! storing the thermal strains 278! 279 if(nalcon(1,imat).eq.1) then 280 do k=1,3 281 eth(k)=alph(1) 282 enddo 283 do k=4,6 284 eth(k)=0.d0 285 enddo 286 elseif(nalcon(1,imat).eq.3) then 287 do k=1,3 288 eth(k)=alph(k) 289 enddo 290 do k=4,6 291 eth(k)=0.d0 292 enddo 293 else 294 do k=1,6 295 eth(k)=alph(k) 296 enddo 297 endif 298! 299! calculating the hardening coefficients 300! 301! for the calculation of the stresses, the whole curve 302! has to be stored: 303! plconloc(2*k-1), k=1...20: equivalent plastic strain values (iso) 304! plconloc(2*k),k=1...20: corresponding stresses (iso) 305! plconloc(39+2*k),k=1...20: equivalent plastic strain values (kin) 306! plconloc(40+2*k),k=1...20: corresponding stresses (kin) 307! 308! initialization 309! 310 if((kode.lt.-50).and.(kode.gt.-100)) then 311 if(npmat_.eq.0) then 312 plconloc(801)=0.5d0 313 plconloc(802)=0.5d0 314 else 315 plconloc(1)=0.d0 316 plconloc(2)=0.d0 317 plconloc(3)=0.d0 318 plconloc(801)=nplicon(1,imat)+0.5d0 319 plconloc(802)=nplkcon(1,imat)+0.5d0 320! 321! isotropic hardening 322! 323 if(nplicon(1,imat).ne.0) then 324! 325 if(nplicon(0,imat).eq.1) then 326 id=-1 327 else 328 call ident2(plicon(0,1,imat),t1l, 329 & nplicon(0,imat),2*npmat_+1,id) 330 endif 331! 332 if(nplicon(0,imat).eq.0) then 333 continue 334 elseif((nplicon(0,imat).eq.1).or.(id.eq.0).or. 335 & (id.eq.nplicon(0,imat))) then 336 if(id.le.0) then 337 itemp=1 338 else 339 itemp=id 340 endif 341 kin=0 342 call plcopy(plicon,nplicon,plconloc,npmat_, 343 & ntmat_,imat,itemp,iel,kin) 344 if((id.eq.0).or.(id.eq.nplicon(0,imat))) then 345 endif 346 else 347 kin=0 348 call plmix(plicon,nplicon,plconloc,npmat_, 349 & ntmat_,imat,id+1,t1l,iel,kin) 350 endif 351 endif 352! 353! kinematic hardening 354! 355 if(nplkcon(1,imat).ne.0) then 356! 357 if(nplkcon(0,imat).eq.1) then 358 id=-1 359 else 360 call ident2(plkcon(0,1,imat),t1l, 361 & nplkcon(0,imat),2*npmat_+1,id) 362 endif 363! 364 if(nplkcon(0,imat).eq.0) then 365 continue 366 elseif((nplkcon(0,imat).eq.1).or.(id.eq.0).or. 367 & (id.eq.nplkcon(0,imat))) then 368 if(id.le.0)then 369 itemp=1 370 else 371 itemp=id 372 endif 373 kin=1 374 call plcopy(plkcon,nplkcon,plconloc,npmat_, 375 & ntmat_,imat,itemp,iel,kin) 376 if((id.eq.0).or.(id.eq.nplkcon(0,imat))) then 377 endif 378 else 379 kin=1 380 call plmix(plkcon,nplkcon,plconloc,npmat_, 381 & ntmat_,imat,id+1,t1l,iel,kin) 382 endif 383 endif 384 endif 385 endif 386! 387! calculating the elastic constants 388! 389 call ident2(elcon(0,1,imat),t1l,nelcon(2,imat),ncmat_+1,id) 390 if(nelcon(2,imat).eq.0) then 391 continue 392 elseif(nelcon(2,imat).eq.1) then 393 do k=1,nelconst 394 elconloc(k)=elcon(k,1,imat) 395 enddo 396 elseif(id.eq.0) then 397 do k=1,nelconst 398 elconloc(k)=elcon(k,1,imat) 399 enddo 400 elseif(id.eq.nelcon(2,imat)) then 401 do k=1,nelconst 402 elconloc(k)=elcon(k,id,imat) 403 enddo 404 else 405 do k=1,nelconst 406 elconloc(k)=elcon(k,id,imat)+ 407 & (elcon(k,id+1,imat)-elcon(k,id,imat))* 408 & (t1l-elcon(0,id,imat))/ 409 & (elcon(0,id+1,imat)-elcon(0,id,imat)) 410 enddo 411 endif 412! 413! modifying the thermal constants if anisotropic and 414! a transformation was defined 415! 416 if((iorien.ne.0).and.(nalcon(1,imat).gt.1)) then 417! 418! calculating the transformation matrix 419! 420 call transformatrix(orab(1,iorien),pgauss,skl) 421! 422! transforming the thermal strain 423! 424 xa(1,1)=eth(1) 425 xa(1,2)=eth(4) 426 xa(1,3)=eth(5) 427 xa(2,1)=eth(4) 428 xa(2,2)=eth(2) 429 xa(2,3)=eth(6) 430 xa(3,1)=eth(5) 431 xa(3,2)=eth(6) 432 xa(3,3)=eth(3) 433! 434 do jj=1,6 435 eth(jj)=0.d0 436 j1=kal(1,jj) 437 j2=kal(2,jj) 438 do j3=1,3 439 do j4=1,3 440 eth(jj)=eth(jj)+ 441 & xa(j3,j4)*skl(j1,j3)*skl(j2,j4) 442 enddo 443 enddo 444 enddo 445 endif 446 endif 447 endif 448! 449 return 450 end 451 452 ! 453 SUBROUTINE us3_csys_cr(xg,tm,tmg) 454 IMPLICIT NONE 455 REAL*8, INTENT(IN) :: xg(3,3) ! element coordinates in global csys 456 REAL*8, INTENT(OUT) :: tm(3,3),tmg(18,18) ! transformation matrix (e1,e2,e3)T 457 REAL*8 :: e1(3),e2(3),e3(3), dl 458 ! 459 INTEGER :: j,i 460 ! 461 tm(:,:) = 0.d0 462 ! element frame (e1,e2,e3) 463 ! e1 = 1 -> 2 464 do j = 1,3 465 e1(j) = xg(2,j)-xg(1,j) 466 enddo 467 ! norm it 468 dl = dsqrt(e1(1)*e1(1)+e1(2)*e1(2)+e1(3)*e1(3)) 469 do j = 1,3 470 e1(j) = e1(j)/dl 471 enddo 472 ! 473 ! e2 = 2 -> 3 (temp) 474 ! 475 do j = 1,3 476 e2(j) = xg(3,j)-xg(1,j) 477 enddo 478 ! 479 ! e3 = e1 x e2 480 ! 481 e3(1) = e1(2)*e2(3) - e1(3)*e2(2) 482 e3(2) = e1(3)*e2(1) - e1(1)*e2(3) 483 e3(3) = e1(1)*e2(2) - e1(2)*e2(1) 484 ! norm it 485 dl = dsqrt(e3(1)*e3(1) + e3(2)*e3(2) + e3(3)*e3(3)) 486 do j = 1,3 487 e3(j) = e3(j)/dl 488 enddo 489 ! 490 ! e2 = e3 x e1 491 ! 492 e2(:) = 0.d0 493 e2(1) = e3(2)*e1(3) - e3(3)*e1(2) 494 e2(2) = e3(3)*e1(1) - e3(1)*e1(3) 495 e2(3) = e3(1)*e1(2) - e3(2)*e1(1) 496 ! norm it 497 dl = dsqrt(e2(1)*e2(1) + e2(2)*e2(2) + e2(3)*e2(3)) 498 do j = 1,3 499 e2(j) = e2(j)/dl 500 enddo 501 ! transformation matrix from the global into the local system 502 do j = 1,3 503 tm(1,j) = e1(j) 504 tm(2,j) = e2(j) 505 tm(3,j) = e3(j) 506 enddo 507 ! 508 tmg(:,:) = 0.d0 509 do i = 1,3 510 do j = 1,3 511 tmg(i,j) = tm(i,j) ! f1 512 tmg(i+3,j+3) = tm(i,j) ! r1 513 tmg(i+6,j+6) = tm(i,j) ! f2 514 tmg(i+9,j+9) = tm(i,j) ! r2 515 tmg(i+12,j+12) = tm(i,j) ! f3 516 tmg(i+15,j+15) = tm(i,j) ! r3 517 enddo 518 enddo 519 ! 520 RETURN 521 END 522 ! 523 SUBROUTINE us3_csys(xg,tm,tmg) 524 IMPLICIT NONE 525 REAL*8, INTENT(IN) :: xg(3,3) ! element coordinates in global csys 526 REAL*8, INTENT(OUT) :: tm(3,3),tmg(18,18) ! transformation matrix (e1,e2,e3)T 527 REAL*8 :: e1(3),e2(3),e3(3),dl,dd,xno(3) 528 REAL*8 :: xi,et,xl(3,8),xs(3,7),p(3),shp(7,8) 529 REAL*8 :: a(3,3) 530 ! 531 INTEGER :: iflag,j,i,l,k 532 ! 533 tm(:,:) = 0.d0 534 ! 535 do i=1,3 536 do j=1,3 537 xl(i,j) = xg(j,i) 538 enddo 539 enddo 540 xi = 0.d0 541 et = 0.d0 542 iflag = 2 543 call shape3tri(xi,et,xl,xno,xs,shp,iflag) 544 dd = dsqrt(xno(1)*xno(1)+xno(2)*xno(2)+xno(3)*xno(3)) 545 do l = 1,3 546 xno(l)=xno(l)/dd 547 enddo 548 ! coordinates of the point at (xi,et)=(0.,0.) 549 do l = 1,3 550 p(l) = 0.d0 551 do k = 1,3 552 p(l) = p(l) + shp(4,k)*xl(l,k) 553 enddo 554 enddo 555 ! unit matrix 556 do k = 1,3 557 do l = 1,3 558 a(k,l) = 0.d0 559 enddo 560 a(k,k) = 1.d0 561 enddo 562 dd = a(1,1)*xno(1)+a(2,1)*xno(2)+a(3,1)*xno(3) 563 if(dabs(dd).gt.0.999999999536d0) then 564 ! project the z-axis 565 dd = a(1,3)*xno(1)+a(2,3)*xno(2)+a(3,3)*xno(3) 566 do l = 1,3 567 e1(l) = a(l,3)-dd*xno(l) 568 enddo 569 else 570 ! project the x-axis 571 do l =1 ,3 572 e1(l) = a(l,1)-dd*xno(l) 573 enddo 574 endif 575 ! 576 dd=dsqrt(e1(1)*e1(1)+e1(2)*e1(2)+e1(3)*e1(3)) 577 do l = 1,3 578 e1(l) = e1(l)/dd 579 enddo 580 ! e2 = n x e1 581 e2(1) = xno(2)*e1(3) - e1(2)*xno(3) 582 e2(2) = xno(3)*e1(1) - e1(3)*xno(1) 583 e2(3) = xno(1)*e1(2) - e1(1)*xno(2) 584 do l = 1,3 585 e3(l) = xno(l) 586 enddo 587 ! put in tm and out 588 do j = 1,3 589 tm(1,j) = e1(j) 590 tm(2,j) = e2(j) 591 tm(3,j) = e3(j) 592 enddo 593 tmg(:,:) = 0.d0 594 do i = 1,3 595 do j = 1,3 596 tmg(i,j) = tm(i,j) ! f1 597 tmg(i+3,j+3) = tm(i,j) ! r1 598 tmg(i+6,j+6) = tm(i,j) ! f2 599 tmg(i+9,j+9) = tm(i,j) ! r2 600 tmg(i+12,j+12) = tm(i,j) ! f3 601 tmg(i+15,j+15) = tm(i,j) ! r3 602 enddo 603 enddo 604 ! 605 RETURN 606 END 607 ! 608 SUBROUTINE us3_Bs(X,Bs) 609 IMPLICIT NONE 610 REAL*8, INTENT(IN) :: X(3,3) 611 REAL*8, INTENT(OUT) :: Bs(2,18) 612 ! 613 REAL*8 :: x1,x2,x3,y1,y2,y3,x31,y31,y12,x21 614 REAL*8 :: Ae,x0,y0,X_1(3),Y_1(3),X_2(3),Y_2(3),Ai 615 REAL*8 :: X_3(3),Y_3(3),bs1(2,6),bs2(2,6),bs3(2,6) 616 REAL*8 :: B1(2,18),B2(2,18),B3(2,18),a3 617 ! 618 Bs1(:,:) = 0.d0 619 Bs2(:,:) = 0.d0 620 Bs3(:,:) = 0.d0 621 Bs(:,:) = 0.d0 622 ! 623 a3 = 1.d0/3.d0 624 ! 625 x1 = X(1,1) 626 x2 = X(2,1) 627 x3 = X(3,1) 628 y1 = X(1,2) 629 y2 = X(2,2) 630 y3 = X(3,2) 631 x31 = x3-x1 632 y31 = y3-y1 633 y12 = y1-y2 634 x21 = x2-x1 635 ! 636 Ae = 0.5d0*(x21*y31-x31*(-y12)) ! area 637 ! 638 x0 = (x1+x2+x3)/3.d0 639 y0 = (y1+y2+y3)/3.d0 640 X_1 = (/X0,x1,x2/) 641 Y_1 = (/Y0,y1,y2/) ! coords tri 1 642 X_2 = (/X0,x2,x3/) 643 Y_2 = (/Y0,y2,y3/) ! coords tri 2 644 X_3 = (/X0,x3,x1/) 645 Y_3 = (/Y0,y3,y1/) ! coords tri 3 646 ! 647 ! cell smoothing 648 ! 649 call us3_CS(X_1,Y_1,bs1,bs2,bs3,Ai) 650 B1(:,1:6) = a3*bs1(:,:) + bs2(:,:) 651 B1(:,7:12) = a3*bs1(:,:) + bs3(:,:) 652 B1(:,13:18) = a3*bs1(:,:) 653 B1 = B1*Ai 654 ! 655 call us3_CS(X_2,Y_2,bs1,bs2,bs3,Ai) 656 B2(:,1:6) = a3*bs1(:,:) 657 B2(:,7:12) = a3*bs1(:,:) + bs2(:,:) 658 B2(:,13:18) = a3*bs1(:,:) + bs3(:,:) 659 B2 = B2*Ai 660 ! 661 call us3_CS(X_3,Y_3,bs1,bs2,bs3,Ai) 662 B3(:,1:6) = a3*bs1(:,:) + bs3(:,:) 663 B3(:,7:12) = a3*bs1(:,:) 664 B3(:,13:18) = a3*bs1(:,:) + bs2(:,:) 665 B3 = B3*Ai 666 ! 667 Bs = (1.d0/Ae)*(B1(:,:)+B2(:,:)+B3(:,:)) 668 ! 669 RETURN 670 END 671 ! 672 SUBROUTINE us3_Kp(X,Db,Ds,Kp) 673 IMPLICIT NONE 674 REAL*8, INTENT(IN) :: X(3,3) 675 REAL*8, INTENT(IN) :: Db(3,3) 676 REAL*8, INTENT(IN) :: Ds(2,2) 677 REAL*8, INTENT(OUT) :: Kp(18,18) 678 ! 679 REAL*8 :: Bs(2,18),Bb(3,18),Ae 680 ! 681 Ae = 0.5d0*((X(2,1)-X(1,1))*(X(3,2)-X(1,2)) 682 & -(X(3,1)-X(1,1))*(-(X(1,2)-X(2,2)))) ! area 683 ! 684 call us3_Bs(X,Bs) 685 call us3_Bb(X(:,1),X(:,2),Bb) 686 ! 687 Kp = (matmul(matmul(transpose(Bs),Ds),Bs) + 688 & matmul(matmul(transpose(Bb),Db),Bb))*Ae 689 ! 690 RETURN 691 END 692 ! 693 SUBROUTINE us3_CS(X,Y,bs1,bs2,bs3,Ae) 694 IMPLICIT NONE 695 REAL*8, INTENT(IN) :: X(3),Y(3) 696 REAL*8, INTENT(OUT) :: bs1(2,6),bs2(2,6),bs3(2,6),Ae 697 ! 698 REAL*8 :: x13,x32,x21,y31,y12,y23 699 REAL*8 :: a1,a2,a3,a4 700 ! 701 bs1(:,:) = 0.d0 702 bs2(:,:) = 0.d0 703 bs3(:,:) = 0.d0 704 ! 705 x21 = X(2)-X(1) 706 x13 = X(1)-X(3) 707 y31 = Y(3)-Y(1) 708 y12 = Y(1)-Y(2) 709 x32 = X(3)-X(2) 710 y23 = Y(2)-Y(3) 711 ! triangle area - XY-tri 712 Ae = 0.5d0*(x21*y31-x13*y12) 713 a1 = 0.5d0*y12*x13 714 a2 = 0.5d0*y31*x21 715 a3 = 0.5d0*x21*x13 716 a4 = 0.5d0*y12*y31 717 bs1(1,3) = +0.5d0*x32/Ae 718 bs1(1,4) = -0.5d0 719 bs1(2,3) = +0.5d0*y23/Ae 720 bs1(2,5) = +0.5d0 721 ! 722 bs2(1,3) = +0.5d0*x13/Ae 723 bs2(1,4) = +0.5d0*a1/Ae 724 bs2(1,5) = +0.5d0*a3/Ae 725 bs2(2,3) = +0.5d0*y31/Ae 726 bs2(2,4) = +0.5d0*a4/Ae 727 bs2(2,5) = +0.5d0*a2/Ae 728 ! 729 bs3(1,3) = +0.5d0*x21/Ae 730 bs3(1,4) = -0.5d0*a2/Ae 731 bs3(1,5) = -0.5d0*a3/Ae 732 bs3(2,3) = +0.5d0*y12/Ae 733 bs3(2,4) = -0.5d0*a4/Ae 734 bs3(2,5) = -0.5d0*a1/Ae 735 ! 736 RETURN 737 END 738 ! 739 SUBROUTINE us3_Bb(X,Y,bb) 740 IMPLICIT NONE 741 REAL*8, INTENT(IN) :: X(3),Y(3) 742 REAL*8, INTENT(OUT) :: bb(3,18) 743 ! 744 REAL*8 :: x13,x31,x21,y31,y12 745 REAL*8 :: y23,x32,Ae 746 REAL*8 :: dNdx1,dNdx2,dNdx3 747 REAL*8 :: dNdy1,dNdy2,dNdy3 748 ! 749 bb(:,:) = 0.d0 750 x21 = X(2)-X(1) 751 x31 = X(3)-X(1) 752 x32 = X(3)-X(2) 753 y23 = Y(2)-Y(3) 754 y31 = Y(3)-Y(1) 755 y12 = Y(1)-Y(2) 756 ! 757 Ae = 0.5d0*(x21*y31-x31*(-y12)) 758 ! 759 dNdx1 = y23/(2.d0*Ae) 760 dNdy1 = x32/(2.d0*Ae) 761 dNdx2 = y31/(2.d0*Ae) 762 dNdy2 = -x31/(2.d0*Ae) 763 dNdx3 = y12/(2.d0*Ae) 764 dNdy3 = x21/(2.d0*Ae) 765 ! 766 bb(1,5) = +dNdx1 767 bb(1,11) = +dNdx2 768 bb(1,17) = +dNdx3 769 bb(2,4) = -dNdy1 770 bb(2,10) = -dNdy2 771 bb(2,16) = -dNdy3 772 bb(3,4) = -dNdx1 773 bb(3,5) = +dNdy1 774 bb(3,10) = -dNdx2 775 bb(3,11) = +dNdy2 776 bb(3,16) = -dNdx3 777 bb(3,17) = +dNdy3 778 ! 779 RETURN 780 END 781 ! 782 ! 783 SUBROUTINE us3_Km(X,K,Qin,h) 784 IMPLICIT NONE 785 REAL*8, INTENT(IN) :: X(3,3),Qin(3,3),h 786 REAL*8, INTENT(OUT) :: K(18,18) 787 ! 788 REAL*8 :: alpha,ab,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9 789 REAL*8 :: x12,x23,x31,y12,y23,y31,A14 790 REAL*8 :: x21,x32,x13,y21,y32,y13,T0(3,9),Te(3,3) 791 REAL*8 :: Ae,A2,A4,LL21,LL32,LL13,h2,L(9,3),V 792 REAL*8 :: Q1(3,3),Q2(3,3),Q3(3,3),Q4(3,3),Q5(3,3),Q6(3,3) 793 REAL*8 :: Enat(3,3),KO(3,3),Kb(9,9),Kh(9,9),Km(9,9) 794 ! 795 INTEGER :: i 796 ! 797 K(:,:) = 0.d0 798 Km(:,:) = 0.d0 799 Kh(:,:) = 0.d0 800 Kb(:,:) = 0.d0 801 L(:,:) = 0.d0 802 T0(:,:) = 0.d0 803 Te(:,:) = 0.d0 804 Q1(:,:) = 0.d0 805 Q2(:,:) = 0.d0 806 Q3(:,:) = 0.d0 807 Q4(:,:) = 0.d0 808 Q5(:,:) = 0.d0 809 Q6(:,:) = 0.d0 810 KO(:,:) = 0.d0 811 ! 812 alpha = 1.d0/8.d0 813 ab = alpha/6.d0 814 b0 = (alpha**2)/4.d0 ! (4.0*(1.0+factor**2)) 815 b1 = 1.d0 816 b2 = 2.d0 817 b3 = 1.d0 818 b4 = 0.d0 819 b5 = 1.d0 820 b6 = -1.d0 821 b7 = -1.d0 822 b8 = -1.d0 823 b9 = -2. 824 ! 825 x12 = X(1,1) - X(2,1) 826 x23 = X(2,1) - X(3,1) 827 x31 = X(3,1) - X(1,1) 828 y12 = X(1,2) - X(2,2) 829 y23 = X(2,2) - X(3,2) 830 y31 = X(3,2) - X(1,2) 831 x21 = -x12 832 x32 = -x23 833 x13 = -x31 834 y21 = -y12 835 y32 = -y23 836 y13 = -y31 837 ! 838 Ae = 0.5d0*(x21*y31-x31*(-y12)) 839 A2 = 2.d0*Ae 840 A4 = 4.d0*Ae 841 h2 = 0.5d0*h 842 V = Ae*h 843 ! 844 LL21 = (x21**2) + (y21**2) 845 LL32 = (x32**2) + (y32**2) 846 LL13 = (x13**2) + (y13**2) 847 ! 848 ! lumping matrix 849 ! 850 L(1,1) = h2*y23 851 L(1,3) = h2*x32 852 L(2,2) = h2*x32 853 L(2,3) = h2*y23 854 L(3,1) = h2*y23*(y13-y21)*ab 855 L(3,2) = h2*x32*(x31-x12)*ab 856 L(3,3) = h2*(x31*y13-x12*y21)*2.d0*ab 857 ! 858 L(4,1) = h2*y31 859 L(4,3) = h2*x13 860 L(5,2) = h2*x13 861 L(5,3) = h2*y31 862 L(6,1) = h2*y31*(y21-y32)*ab 863 L(6,2) = h2*x13*(x12-x23)*ab 864 L(6,3) = h2*(x12*y21-x23*y32)*2.d0*ab 865 ! 866 L(7,1) = h2*y12 867 L(7,3) = h2*x21 868 L(8,2) = h2*x21 869 L(8,3) = h2*y12 870 L(9,1) = h2*y12*(y32-y13)*ab 871 L(9,2) = h2*x21*(x23-x31)*ab 872 L(9,3) = h2*(x23*y32-x31*y13)*2.d0*ab 873 ! 874 ! basic stiffness 875 ! 876 Kb = matmul(matmul(L,Qin),transpose(L))/V 877 ! 878 ! trasformation hierachical rotations 879 ! 880 T0(1,1) = x32/A4 881 T0(1,2) = y32/A4 882 T0(1,3) = 1.d0 883 T0(1,4) = x13/A4 884 T0(1,5) = y13/A4 885 T0(1,7) = x21/A4 886 T0(1,8) = y21/A4 887 ! 888 T0(2,1) = x32/A4 889 T0(2,2) = y32/A4 890 T0(2,4) = x13/A4 891 T0(2,5) = y13/A4 892 T0(2,6) = 1.d0 893 T0(2,7) = x21/A4 894 T0(2,8) = y21/A4 895 ! 896 T0(3,1) = x32/A4 897 T0(3,2) = y32/A4 898 T0(3,4) = x13/A4 899 T0(3,5) = y13/A4 900 T0(3,7) = x21/A4 901 T0(3,8) = y21/A4 902 T0(3,9) = 1.d0 903 ! 904 ! transformation natural pattern 905 ! 906 A14 = (1.d0/(Ae*A4)) 907 Te(1,1) = A14*y23*y13*LL21 908 Te(1,2) = A14*y31*y21*LL32 909 Te(1,3) = A14*y12*y32*LL13 910 Te(2,1) = A14*x23*x13*LL21 911 Te(2,2) = A14*x31*x21*LL32 912 Te(2,3) = A14*x12*x32*LL13 913 Te(3,1) = A14*(y23*x31+x32*y13)*LL21 914 Te(3,2) = A14*(y31*x12+x13*y21)*LL32 915 Te(3,3) = A14*(y12*x23+x21*y32)*LL13 916 ! 917 A14 = (A2/3.d0) 918 ! 919 ! nodal strain-displ.-matrix 920 ! 921 Q1(1,1) = A14*b1/LL21 922 Q1(1,2) = A14*b2/LL21 923 Q1(1,3) = A14*b3/LL21 924 Q1(2,1) = A14*b4/LL32 925 Q1(2,2) = A14*b5/LL32 926 Q1(2,3) = A14*b6/LL32 927 Q1(3,1) = A14*b7/LL13 928 Q1(3,2) = A14*b8/LL13 929 Q1(3,3) = A14*b9/LL13 930 ! 931 Q2(1,1) = A14*b9/LL21 932 Q2(1,2) = A14*b7/LL21 933 Q2(1,3) = A14*b8/LL21 934 Q2(2,1) = A14*b3/LL32 935 Q2(2,2) = A14*b1/LL32 936 Q2(2,3) = A14*b2/LL32 937 Q2(3,1) = A14*b6/LL13 938 Q2(3,2) = A14*b4/LL13 939 Q2(3,3) = A14*b5/LL13 940 ! 941 Q3(1,1) = A14*b5/LL21 942 Q3(1,2) = A14*b6/LL21 943 Q3(1,3) = A14*b4/LL21 944 Q3(2,1) = A14*b8/LL32 945 Q3(2,2) = A14*b9/LL32 946 Q3(2,3) = A14*b7/LL32 947 Q3(3,1) = A14*b2/LL13 948 Q3(3,2) = A14*b3/LL13 949 Q3(3,3) = A14*b1/LL13 950 ! 951 Q4 = (Q1 + Q2)*0.5d0 952 Q5 = (Q2 + Q3)*0.5d0 953 Q6 = (Q3 + Q1)*0.5d0 954 ! 955 Enat = matmul(matmul(transpose(Te),Qin),Te) 956 ! 957 ! higher stiffness with respect to hier...rots 958 ! 959 KO = (3.d0/4.d0)*b0*V*(matmul(matmul(transpose(Q4),Enat),Q4) 960 & + matmul(matmul(transpose(Q5),Enat),Q5) 961 & + matmul(matmul(transpose(Q6),Enat),Q6)) 962 ! higher stiffness [18x18) 963 Kh = matmul(matmul(transpose(T0),KO),T0) 964 Km = Kb + Kh 965 ! 966 do i=1,3 967 K(1+(i-1)*6,1:2) = Km(1+(i-1)*3,1:2) 968 K(1+(i-1)*6,6) = Km(1+(i-1)*3,3) 969 K(1+(i-1)*6,7:8) = Km(1+(i-1)*3,4:5) 970 K(1+(i-1)*6,12) = Km(1+(i-1)*3,6) 971 K(1+(i-1)*6,13:14) = Km(1+(i-1)*3,7:8) 972 K(1+(i-1)*6,18) = Km(1+(i-1)*3,9) 973 !zeile 2,8,14: 974 K(2+(i-1)*6,1:2) = Km(2+(i-1)*3,1:2) 975 K(2+(i-1)*6,6) = Km(2+(i-1)*3,3) 976 K(2+(i-1)*6,7:8) = Km(2+(i-1)*3,4:5) 977 K(2+(i-1)*6,12) = Km(2+(i-1)*3,6) 978 K(2+(i-1)*6,13:14) = Km(2+(i-1)*3,7:8) 979 K(2+(i-1)*6,18) = Km(2+(i-1)*3,9) 980 ! 981 K(6+(i-1)*6,1:2) = Km(3+(i-1)*3,1:2) 982 K(6+(i-1)*6,6) = Km(3+(i-1)*3,3) 983 K(6+(i-1)*6,7:8) = Km(3+(i-1)*3,4:5) 984 K(6+(i-1)*6,12) = Km(3+(i-1)*3,6) 985 K(6+(i-1)*6,13:14) = Km(3+(i-1)*3,7:8) 986 K(6+(i-1)*6,18) = Km(3+(i-1)*3,9) 987 enddo 988 ! 989 RETURN 990 END 991 ! 992 SUBROUTINE bm_ANDES(X,bm,Qin,h,r,s) 993 IMPLICIT NONE 994 REAL*8, INTENT(IN) :: X(3,3),Qin(3,3),h,r,s 995 REAL*8, INTENT(OUT) :: bm(3,18) 996 ! 997 REAL*8 :: alpha,ab,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9 998 REAL*8 :: x12,x23,x31,y12,y23,y31,A14 999 REAL*8 :: x21,x32,x13,y21,y32,y13,T0(3,9),Te(3,3) 1000 REAL*8 :: Ae,A2,A4,LL21,LL32,LL13,h2,L(9,3),V 1001 REAL*8 :: Q1(3,3),Q2(3,3),Q3(3,3),Q4(3,3),Q5(3,3),Q6(3,3) 1002 REAL*8 :: Enat(3,3),matQ(3,3),bm2(3,9) 1003 ! 1004 INTEGER :: i,i3,i6 1005 ! 1006 L(:,:) = 0.d0 1007 T0(:,:) = 0.d0 1008 Te(:,:) = 0.d0 1009 Q1(:,:) = 0.d0 1010 Q2(:,:) = 0.d0 1011 Q3(:,:) = 0.d0 1012 Q4(:,:) = 0.d0 1013 Q5(:,:) = 0.d0 1014 Q6(:,:) = 0.d0 1015 ! 1016 alpha = 1.d0/8.d0 1017 ab = alpha/6.d0 1018 b0 = (alpha**2)/4.d0 ! (4.0*(1.0+factor**2)) 1019 b1 = 1.d0 1020 b2 = 2.d0 1021 b3 = 1.d0 1022 b4 = 0.d0 1023 b5 = 1.d0 1024 b6 = -1.d0 1025 b7 = -1.d0 1026 b8 = -1.d0 1027 b9 = -2.d0 1028 ! 1029 x12 = X(1,1) - X(2,1) 1030 x23 = X(2,1) - X(3,1) 1031 x31 = X(3,1) - X(1,1) 1032 y12 = X(1,2) - X(2,2) 1033 y23 = X(2,2) - X(3,2) 1034 y31 = X(3,2) - X(1,2) 1035 x21 = -x12 1036 x32 = -x23 1037 x13 = -x31 1038 y21 = -y12 1039 y32 = -y23 1040 y13 = -y31 1041 ! 1042 Ae = 0.5d0*(x21*y31-x31*(-y12)) 1043 A2 = 2.d0*Ae 1044 A4 = 4.d0*Ae 1045 h2 = 0.5d0*h 1046 V = Ae*h 1047 ! 1048 LL21 = (x21**2) + (y21**2) 1049 LL32 = (x32**2) + (y32**2) 1050 LL13 = (x13**2) + (y13**2) 1051 ! 1052 ! lumping matrix 1053 ! 1054 L(1,1) = h2*y23 1055 L(1,3) = h2*x32 1056 L(2,2) = h2*x32 1057 L(2,3) = h2*y23 1058 L(3,1) = h2*y23*(y13-y21)*ab 1059 L(3,2) = h2*x32*(x31-x12)*ab 1060 L(3,3) = h2*(x31*y13-x12*y21)*2.d0*ab 1061 ! 1062 L(4,1) = h2*y31 1063 L(4,3) = h2*x13 1064 L(5,2) = h2*x13 1065 L(5,3) = h2*y31 1066 L(6,1) = h2*y31*(y21-y32)*ab 1067 L(6,2) = h2*x13*(x12-x23)*ab 1068 L(6,3) = h2*(x12*y21-x23*y32)*2.d0*ab 1069 ! 1070 L(7,1) = h2*y12 1071 L(7,3) = h2*x21 1072 L(8,2) = h2*x21 1073 L(8,3) = h2*y12 1074 L(9,1) = h2*y12*(y32-y13)*ab 1075 L(9,2) = h2*x21*(x23-x31)*ab 1076 L(9,3) = h2*(x23*y32-x31*y13)*2.d0*ab 1077 ! 1078 ! 1079 ! trasformation hierachical rotations 1080 ! 1081 T0(1,1) = x32/A4 1082 T0(1,2) = y32/A4 1083 T0(1,3) = 1.d0 1084 T0(1,4) = x13/A4 1085 T0(1,5) = y13/A4 1086 T0(1,7) = x21/A4 1087 T0(1,8) = y21/A4 1088 ! 1089 T0(2,1) = x32/A4 1090 T0(2,2) = y32/A4 1091 T0(2,4) = x13/A4 1092 T0(2,5) = y13/A4 1093 T0(2,6) = 1.d0 1094 T0(2,7) = x21/A4 1095 T0(2,8) = y21/A4 1096 ! 1097 T0(3,1) = x32/A4 1098 T0(3,2) = y32/A4 1099 T0(3,4) = x13/A4 1100 T0(3,5) = y13/A4 1101 T0(3,7) = x21/A4 1102 T0(3,8) = y21/A4 1103 T0(3,9) = 1.d0 1104 ! 1105 ! transformation natural pattern 1106 ! 1107 A14 = (1.d0/(Ae*A4)) 1108 Te(1,1) = A14*y23*y13*LL21 1109 Te(1,2) = A14*y31*y21*LL32 1110 Te(1,3) = A14*y12*y32*LL13 1111 Te(2,1) = A14*x23*x13*LL21 1112 Te(2,2) = A14*x31*x21*LL32 1113 Te(2,3) = A14*x12*x32*LL13 1114 Te(3,1) = A14*(y23*x31+x32*y13)*LL21 1115 Te(3,2) = A14*(y31*x12+x13*y21)*LL32 1116 Te(3,3) = A14*(y12*x23+x21*y32)*LL13 1117 ! 1118 A14 = (A2/3.d0) 1119 ! 1120 ! nodal strain-displ.-matrix 1121 ! 1122 1123 Q1(1,1) = A14*b1/LL21 1124 Q1(1,2) = A14*b2/LL21 1125 Q1(1,3) = A14*b3/LL21 1126 Q1(2,1) = A14*b4/LL32 1127 Q1(2,2) = A14*b5/LL32 1128 Q1(2,3) = A14*b6/LL32 1129 Q1(3,1) = A14*b7/LL13 1130 Q1(3,2) = A14*b8/LL13 1131 Q1(3,3) = A14*b9/LL13 1132 ! 1133 Q2(1,1) = A14*b9/LL21 1134 Q2(1,2) = A14*b7/LL21 1135 Q2(1,3) = A14*b8/LL21 1136 Q2(2,1) = A14*b3/LL32 1137 Q2(2,2) = A14*b1/LL32 1138 Q2(2,3) = A14*b2/LL32 1139 Q2(3,1) = A14*b6/LL13 1140 Q2(3,2) = A14*b4/LL13 1141 Q2(3,3) = A14*b5/LL13 1142 ! 1143 Q3(1,1) = A14*b5/LL21 1144 Q3(1,2) = A14*b6/LL21 1145 Q3(1,3) = A14*b4/LL21 1146 Q3(2,1) = A14*b8/LL32 1147 Q3(2,2) = A14*b9/LL32 1148 Q3(2,3) = A14*b7/LL32 1149 Q3(3,1) = A14*b2/LL13 1150 Q3(3,2) = A14*b3/LL13 1151 Q3(3,3) = A14*b1/LL13 1152 ! 1153 Q4 = (Q1 + Q2)*0.5d0 1154 Q5 = (Q2 + Q3)*0.5d0 1155 Q6 = (Q3 + Q1)*0.5d0 1156 ! 1157 Enat = matmul(matmul(transpose(Te),Qin),Te) 1158 ! 1159 matQ = (1.d0-r-s)*Q1 + r*Q2 + s*Q3 1160 bm2 = (3.d0/2.d0)*(b0**0.5)*matmul(matmul(Te,matQ),T0) 1161 bm2 = bm2 + transpose(L)/V 1162 ! 1163 bm(:,:) = 0.d0 1164 ! 1165 do i=1,3 1166 ! 1167 i3 = (i-1)*3 1168 i6 = (i-1)*6 1169 ! 1170 bm(1:3,1+i6:2+i6) = bm2(1:3,1+i3:2+i3) 1171 bm(1:3,6+i6) = bm2(1:3,3+i3) 1172 ! 1173 enddo 1174 ! 1175 RETURN 1176 END 1177 ! 1178 SUBROUTINE us3_matma(e,un,h,Dm,Db,Ds,Qin,Qs,x,iflag) 1179 IMPLICIT NONE 1180 REAL*8, INTENT(IN) :: e,un,h,x(3,3) 1181 INTEGER, INTENT(IN) :: iflag 1182 REAL*8, INTENT(OUT) :: Dm(3,3),Db(3,3),Ds(2,2),Qin(3,3),Qs(2,2) 1183 REAL*8 :: q1,kap,hei(3),he 1184 INTEGER :: k 1185 ! 1186 ! Reduced material matrix (S33=0)+ integration in thichkness dirc. 1187 ! pre- in thickness dirc. integrated: 1188 if(iflag.eq.1) then 1189 Qin(:,:) = 0.d0 1190 q1 = e/(1.d0-un**2) 1191 Qin(1,1) = q1 1192 Qin(1,2) = q1*un 1193 Qin(2,1) = q1*un 1194 Qin(2,2) = q1 1195 Qin(3,3) = q1*(1.d0-un)/2.d0 1196 ! membrane: 1197 Dm = Qin*h 1198 ! bending: 1199 Db = Qin*h**3/12.d0 1200 ! shear 1201 1202 ! 1203 ! shear correction 1204 ! 1205 ! l1 = 1->2 1206! hei(1) = dsqrt((x(2,1)-x(1,1))**2+(x(2,2)-x(1,2))**2) 1207! ! l1 = 2->3 1208! hei(2) = dsqrt((x(3,1)-x(2,1))**2+(x(3,2)-x(2,2))**2) 1209! ! l1 = 3->1 1210! hei(3) = dsqrt((x(3,1)-x(1,1))**2+(x(3,2)-x(1,2))**2) 1211! ! 1212! he = 0.d0 1213! do k = 1,3 1214! if(he.LT.abs(hei(k))) then 1215! he = abs(hei(k)) 1216! endif 1217! enddo 1218 ! 1219! kap = 5.d0/6.d0 1220! q1 = (h**3*kap/(h**2+0.1d0*he**2)*e/2.0/(1+un)) 1221! Ds(:,:) = 0.0 1222! Ds(1,1) = q1 1223! Ds(2,2) = q1 1224! Qs(:,:) = 0.d0 1225! q1 = (h**2*kap/(h**2+0.1d0*he**2)*e/2.0/(1+un)) 1226! Qs(1,1) = q1 1227! Qs(2,2) = q1 1228 1229 Qs(:,:) = 0.d0 1230 kap = 5.d0/6.d0 1231 q1 = e/(2.d0*(1.d0+un)) 1232 Qs(1,1) = q1*kap 1233 Qs(2,2) = q1*kap 1234 Ds = Qs*h 1235 else 1236 Qin(:,:) = 0.d0 1237 q1 = e/(1.d0-un**2) 1238 Qin(1,1) = q1 1239 Qin(1,2) = q1*un 1240 Qin(2,1) = q1*un 1241 Qin(2,2) = q1 1242 Qin(3,3) = q1*(1.d0-un)/2.d0 1243 Qs(:,:) = 0.d0 1244 kap = 5.d0/6.d0 1245 q1 = e/(2.d0*(1.d0+un)) 1246 Qs(1,1) = q1*kap 1247 Qs(2,2) = q1*kap 1248 endif 1249 RETURN 1250 END 1251 ! 1252 SUBROUTINE us3_M(X,h,rho,M) 1253 IMPLICIT NONE 1254 REAL*8, INTENT(IN) :: X(3,3),rho,h 1255 REAL*8, INTENT(OUT) :: M(18,18) 1256 REAL*8 :: MLST(12,12),Ae,Ai,TCH(12,9),alpha 1257 REAL*8 :: Mm(9,9),points3(3,2),w3,r,s,Nrs(3) 1258 REAL*8 :: N_u(6,18),q1,m_3t(6,6) 1259 REAL*8 :: x12,x23,x31,y12,y23,y31 1260 REAL*8 :: x21,x32,x13,y21,y32,y13 1261 INTEGER :: i,j 1262 ! 1263 alpha = 1.d0/8.d0 1264 ! 1265 points3(1,1) = 0.1666666666667 1266 points3(1,2) = 0.1666666666667 1267 points3(2,1) = 0.6666666666667 1268 points3(2,2) = 0.1666666666667 1269 points3(3,1) = 0.1666666666667 1270 points3(3,2) = 0.6666666666667 1271 w3 = 0.333333333333333 1272 ! 1273 x12 = X(1,1) - X(2,1) 1274 x23 = X(2,1) - X(3,1) 1275 x31 = X(3,1) - X(1,1) 1276 y12 = X(1,2) - X(2,2) 1277 y23 = X(2,2) - X(3,2) 1278 y31 = X(3,2) - X(1,2) 1279 x21 = -x12 1280 x32 = -x23 1281 x13 = -x31 1282 y21 = -y12 1283 y32 = -y23 1284 y13 = -y31 1285! ! 1286 Ae = 0.5d0*(x21*y31-x31*(-y12)) 1287! Ai = Ae*h*rho/180.d0 1288! ! 1289! MLST(:,:) = 0.d0 1290! MLST(1,1)=+Ai*6.d0 1291! MLST(1,3)=-Ai*1.d0 1292! MLST(1,5)=-Ai*1.d0 1293! MLST(1,9)=-Ai*4.d0 1294! MLST(2,2)=+Ai*6.d0 1295! MLST(2,4)=-Ai*1.d0 1296! MLST(2,6)=-Ai*1.d0 1297! MLST(2,10)=-Ai*4.d0 1298! MLST(3,5)=-Ai*4.d0 1299! MLST(3,7)=+Ai*32.d0 1300! MLST(3,9)=+Ai*16.d0 1301! MLST(3,11)=+Ai*16.d0 1302! MLST(4,6)=-Ai*4.d0 1303! MLST(4,8)=+Ai*32.d0 1304! MLST(4,10)=+Ai*16.d0 1305! MLST(4,12)=+Ai*16.d0 1306! MLST(5,1)=-Ai*1.d0 1307! MLST(5,3) =+Ai*6.d0 1308! MLST(5,5)=-Ai*1.d0 1309! MLST(5,11)=-Ai*4.d0 1310! MLST(6,2)=-Ai*1.d0 1311! MLST(6,4) =+Ai*6.d0 1312! MLST(6,6)=-Ai*1.d0 1313! MLST(6,12)=-Ai*4.d0 1314! MLST(7,1)=-Ai*4.d0 1315! MLST(7,7) =+Ai*16.d0 1316! MLST(7,9)=+Ai*32.d0 1317! MLST(7,11)=+Ai*16.d0 1318! MLST(8,2)=-Ai*4.d0 1319! MLST(8,8) =+Ai*16.d0 1320! MLST(8,10)=+Ai*32.d0 1321! MLST(8,12)=+Ai*16.d0 1322! MLST(9,1)=-Ai*1.d0 1323! MLST(9,3) =-Ai*1.d0 1324! MLST(9,5)=+Ai*6.d0 1325! MLST(9,7) =-Ai*4.d0 1326! MLST(10,2)=-Ai*1.d0 1327! MLST(10,4) =-Ai*1.d0 1328! MLST(10,6)=+Ai*6.d0 1329! MLST(10,8) =-Ai*4.d0 1330! MLST(11,3)=-Ai*4.d0 1331! MLST(11,7) =+Ai*16.d0 1332! MLST(11,9)=+Ai*16.d0 1333! MLST(11,11) =+Ai*32.d0 1334! MLST(12,4)=-Ai*4.d0 1335! MLST(12,8) =+Ai*16.d0 1336! MLST(12,10)=+Ai*16.d0 1337! MLST(12,12) =+Ai*32.d0 1338! ! 1339! TCH(:,:) = 0.d0 1340! TCH(1,1) = 1.d0 1341! TCH(2,2) = 1.d0 1342! TCH(3,1) = 0.5d0 1343! TCH(3,3) = 0.125d0*alpha*y12 1344! TCH(3,6) = 0.125d0*alpha*y21 1345! TCH(3,7) = 0.5d0 1346! TCH(4,2) = 0.5d0 1347! TCH(4,3) = 0.125d0*alpha*x21 1348! TCH(4,6) = 0.125d0*alpha*x12 1349! TCH(4,8) = 0.5d0 1350! TCH(5,4) = 1.d0 1351! TCH(6,5) = 1.d0 1352! TCH(7,4) = 0.5d0 1353! TCH(7,6) = 0.125d0*alpha*y23 1354! TCH(7,9) = 0.125d0*alpha*y32 1355! TCH(7,7) = 0.5d0 1356! TCH(8,5) = 0.5d0 1357! TCH(8,6) = 0.125d0*alpha*x32 1358! TCH(8,9) = 0.125d0*alpha*x23 1359! TCH(8,8) = 0.5d0 1360! TCH(9,7) = 1.d0 1361! TCH(10,8) = 1.d0 1362! TCH(11,1) = 0.5d0 1363! TCH(11,3) = 0.125d0*alpha*y13 1364! TCH(11,9) = 0.125d0*alpha*y31 1365! TCH(11,7) = 0.5d0 1366! TCH(12,2) = 0.5d0 1367! TCH(12,3) = 0.125d0*alpha*x31 1368! TCH(12,9) = 0.125d0*alpha*x13 1369! TCH(12,8) = 0.5d0 1370! ! 1371! Mm = matmul(matmul(transpose(TCH),MLST),TCH) 1372! ! 1373! do i=1,3 1374! M(1+(i-1)*6,1:2) = Mm(1+(i-1)*3,1:2) 1375! M(1+(i-1)*6,6) = Mm(1+(i-1)*3,3) 1376! M(1+(i-1)*6,7:8) = Mm(1+(i-1)*3,4:5) 1377! M(1+(i-1)*6,12) = Mm(1+(i-1)*3,6) 1378! M(1+(i-1)*6,13:14) = Mm(1+(i-1)*3,7:8) 1379! M(1+(i-1)*6,18) = Mm(1+(i-1)*3,9) 1380! !zeile 2,8,14: 1381! M(2+(i-1)*6,1:2) = Mm(2+(i-1)*3,1:2) 1382! M(2+(i-1)*6,6) = Mm(2+(i-1)*3,3) 1383! M(2+(i-1)*6,7:8) = Mm(2+(i-1)*3,4:5) 1384! M(2+(i-1)*6,12) = Mm(2+(i-1)*3,6) 1385! M(2+(i-1)*6,13:14) = Mm(2+(i-1)*3,7:8) 1386! M(2+(i-1)*6,18) = Mm(2+(i-1)*3,9) 1387! ! 1388! M(6+(i-1)*6,1:2) = Mm(3+(i-1)*3,1:2) 1389! M(6+(i-1)*6,6) = Mm(3+(i-1)*3,3) 1390! M(6+(i-1)*6,7:8) = Mm(3+(i-1)*3,4:5) 1391! M(6+(i-1)*6,12) = Mm(3+(i-1)*3,6) 1392! M(6+(i-1)*6,13:14) = Mm(3+(i-1)*3,7:8) 1393! M(6+(i-1)*6,18) = Mm(3+(i-1)*3,9) 1394! enddo 1395! ! 1396 M(:,:) = 0.d0 1397 ! 1398 N_u(:,:) = 0.d0 1399 m_3t(:,:) = 0.d0 1400 q1 = rho*h 1401 m_3t(1,1) = q1 1402 m_3t(2,2) = q1 1403 m_3t(3,3) = q1 1404 q1 = (rho*h**3)/12.d0 1405 m_3t(4,4) = q1 1406 m_3t(5,5) = q1 1407 ! 1408 do i=1,3 1409 r = points3(i,1) 1410 s = points3(i,2) 1411 Nrs(1) = 1.d0-r-s 1412 Nrs(2) = r 1413 Nrs(3) = s 1414 ! 1415 do j = 1,3 1416 N_u(1,1+(j-1)*6) = Nrs(j) 1417 N_u(2,2+(j-1)*6) = Nrs(j) 1418 N_u(3,3+(j-1)*6) = Nrs(j) 1419 N_u(4,4+(j-1)*6) = Nrs(j) 1420 N_u(5,5+(j-1)*6) = Nrs(j) 1421 enddo 1422 ! 1423 M = M + matmul(matmul(transpose(N_u),m_3t),N_u)*Ae*w3 1424 enddo 1425 ! 1426 RETURN 1427 END 1428 ! 1429 ! 1430 SUBROUTINE us3_ae(X,Ae) 1431 IMPLICIT NONE 1432 REAL*8, INTENT(IN) :: X(3,3) 1433 REAL*8, INTENT(OUT) :: Ae 1434 ! 1435 Ae = 0.5d0*((X(2,1)-X(1,1))*(X(3,2)-X(1,2)) 1436 & -(X(3,1)-X(1,1))*(-(X(1,2)-X(2,2)))) ! area 1437 ! 1438 RETURN 1439 END 1440 ! 1441 SUBROUTINE us3_xu(x,ushell,ueg,xg,tm,vl) 1442 IMPLICIT NONE 1443 REAL*8, INTENT(IN) :: xg(3,3),tm(3,3),vl(6,3) 1444 REAL*8, INTENT(OUT) :: x(3,3),ushell(18),ueg(18) 1445 ! 1446 x(:,:) = 0.d0 1447 x(1,:) = matmul(tm,xg(1,:)) 1448 x(2,:) = matmul(tm,xg(2,:)) 1449 x(3,:) = matmul(tm,xg(3,:)) 1450 ! 1451 ushell(1:3) = matmul(tm,vl(1:3,1)) 1452 ushell(4:6) = matmul(tm,vl(4:6,1)) 1453 ushell(7:9) = matmul(tm,vl(1:3,2)) 1454 ushell(10:12) = matmul(tm,vl(4:6,2)) 1455 ushell(13:15) = matmul(tm,vl(1:3,3)) 1456 ushell(16:18) = matmul(tm,vl(4:6,3)) 1457 ! 1458 ueg(:) = 0.d0 1459 ueg(1:3) = vl(1:3,1) 1460 ueg(4:6) = vl(4:6,1) 1461 ueg(7:9) = vl(1:3,2) 1462 ueg(10:12) = vl(4:6,2) 1463 ueg(13:15) = vl(1:3,3) 1464 ueg(16:18) = vl(4:6,3) 1465 ! 1466 RETURN 1467 END 1468 ! 1469 SUBROUTINE us3_linel_Qi(e,un,Qin,Qs) 1470 IMPLICIT NONE 1471 REAL*8, INTENT(IN) :: e,un 1472 REAL*8, INTENT(OUT) :: Qin(3,3),Qs(2,2) 1473 REAL*8 :: q1,kap 1474 INTEGER :: k 1475 ! 1476 ! Reduced material matrix (S33=0) 1477 ! 1478 Qin(:,:) = 0.d0 1479 q1 = e/(1.d0-un**2) 1480 Qin(1,1) = q1 1481 Qin(1,2) = q1*un 1482 Qin(2,1) = q1*un 1483 Qin(2,2) = q1 1484 Qin(3,3) = q1*(1.d0-un)/2.d0 1485 ! 1486 Qs(:,:) = 0.d0 1487 kap = 5.d0/6.d0 1488 q1 = e/(2.d0*(1.d0+un)) 1489 Qs(1,1) = q1*kap 1490 Qs(2,2) = q1*kap 1491 RETURN 1492 END 1493 ! 1494