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 e_c3d_rhs(co,nk,konl,lakonl,p1,p2,omx,bodyfx,nbody, 20 & ff,nelem,nmethod,rhcon,ielmat,ntmat_,vold,iperturb,nelemload, 21 & sideload,xload,nload,idist,ttime,time,istep,iinc,dtime, 22 & xloadold,reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc, 23 & veold,matname,mi,ielprop,prop) 24! 25! computation of the rhs for the element with 26! the topology in konl: only for nonlinear calculations (i.e. 27! the field ff contains no temperature and eigenstress 28! contributions) 29! 30! RESTRICTION: the only material parameter occurring in the 31! loading is the density. In the present subroutine the density 32! is assumed to be temperature INDEPENDENT. 33! 34 implicit none 35! 36 logical ivolumeforce 37! 38 character*8 lakonl 39 character*20 sideload(*) 40 character*80 matname(*),amat 41! 42 integer mi(*),nk,konl(20),ielmat(mi(3),*),nbody,ifaceq(8,6), 43 & nelemload(2,*),ielprop(*),null, 44 & nelem,nmethod,iperturb(*),nload,idist,i,i1,j1,jj,jj1,id,kk, 45 & ipointer,nope,nopes,j,k,ntmat_,i2,imat,ii,ig,mint2d,mint3d, 46 & ifacet(6,4),ifacew(8,5),istep,iinc,layer,kspt,jltyp,iflag, 47 & ipompc(*),nodempc(3,*),nmpc,ikmpc(*),ilmpc(*),iscale 48! 49 real*8 co(3,*),p1(3,2),p2(3,2),omx(2),bodyfx(3),veold(0:mi(2),*), 50 & rhcon(0:1,ntmat_,*),xs2(3,7),xloadold(2,*),reltime, 51 & bodyf(3),om(2),rho,bf(3),q(3),shpj(4,20),xl(3,20), 52 & shp(4,20),voldl(3,20),xl2(3,8),xsj2(3),shp2(7,8), 53 & vold(0:mi(2),*),prop(*), 54 & xload(2,*),xi,et,ze,const,xsj,ff(60),weight,ttime,time,tvar(2), 55 & coords(3),dtime,coefmpc(*) 56! 57 include "gauss.f" 58! 59 data ifaceq /4,3,2,1,11,10,9,12, 60 & 5,6,7,8,13,14,15,16, 61 & 1,2,6,5,9,18,13,17, 62 & 2,3,7,6,10,19,14,18, 63 & 3,4,8,7,11,20,15,19, 64 & 4,1,5,8,12,17,16,20/ 65 data ifacet /1,3,2,7,6,5, 66 & 1,2,4,5,9,8, 67 & 2,3,4,6,10,9, 68 & 1,4,3,8,10,7/ 69 data ifacew /1,3,2,9,8,7,0,0, 70 & 4,5,6,10,11,12,0,0, 71 & 1,2,5,4,7,14,10,13, 72 & 2,3,6,5,8,15,11,14, 73 & 4,6,3,1,12,15,9,13/ 74 data iflag /3/ 75 data null /0/ 76! 77 tvar(1)=time 78 tvar(2)=ttime+time 79! 80 if(lakonl(4:4).eq.'2') then 81 nope=20 82 nopes=8 83 elseif(lakonl(4:4).eq.'8') then 84 nope=8 85 nopes=4 86 elseif(lakonl(4:5).eq.'10') then 87 nope=10 88 nopes=6 89 elseif(lakonl(4:4).eq.'4') then 90 nope=4 91 nopes=3 92 elseif(lakonl(4:5).eq.'15') then 93 nope=15 94 else 95 nope=6 96 endif 97! 98 if(lakonl(4:5).eq.'8R') then 99 mint2d=1 100 mint3d=1 101 elseif(lakonl(4:7).eq.'20RB') then 102 if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then 103 mint2d=4 104 mint3d=50 105 else 106 mint2d=4 107 call beamintscheme(lakonl,mint3d,ielprop(nelem),prop, 108 & null,xi,et,ze,weight) 109 endif 110 elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) then 111 mint2d=4 112 mint3d=8 113 elseif(lakonl(4:4).eq.'2') then 114 mint2d=9 115 mint3d=27 116 elseif(lakonl(4:5).eq.'10') then 117 mint2d=3 118 mint3d=4 119 elseif(lakonl(4:4).eq.'4') then 120 mint2d=1 121 mint3d=1 122 elseif(lakonl(4:5).eq.'15') then 123 mint3d=9 124 else 125 mint3d=2 126 endif 127! 128! computation of the coordinates of the local nodes 129! 130 do i=1,nope 131 do j=1,3 132 xl(j,i)=co(j,konl(i)) 133 enddo 134 enddo 135! 136! initialisation for distributed forces 137! 138c if(idist.ne.0) then 139 do i=1,3*nope 140 ff(i)=0.d0 141 enddo 142c endif 143! 144! displacements for 2nd order static and modal theory 145! 146c if(iperturb.ne.0) then 147 if((iperturb(1).eq.1).or.(iperturb(2).eq.1)) then 148 do i1=1,nope 149 do i2=1,3 150 voldl(i2,i1)=vold(i2,konl(i1)) 151 enddo 152 enddo 153 endif 154! 155! calculating the density: no temperature dependence assumed! 156! otherwise cfr. e_c3d.f 157! 158 imat=ielmat(1,nelem) 159 amat=matname(imat) 160 rho=rhcon(1,1,imat) 161! 162! computation of the body forces 163! 164 ivolumeforce=.false. 165 if(nbody.ne.0) then 166 ivolumeforce=.true. 167 do ii=1,2 168 om(ii)=omx(ii)*rho 169 enddo 170 do ii=1,3 171 bodyf(ii)=bodyfx(ii)*rho 172 enddo 173 elseif(nload.gt.0) then 174 call nident2(nelemload,nelem,nload,id) 175 do 176 if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit 177 if((sideload(id)(1:2).ne.'BX').and. 178 & (sideload(id)(1:2).ne.'BY').and. 179 & (sideload(id)(1:2).ne.'BZ')) then 180 id=id-1 181 cycle 182 else 183 ivolumeforce=.true. 184 exit 185 endif 186 enddo 187 endif 188! 189! computation of the rhs: loop over the Gauss points 190! 191 if(ivolumeforce) then 192 do kk=1,mint3d 193 if(lakonl(4:5).eq.'8R') then 194 xi=gauss3d1(1,kk) 195 et=gauss3d1(2,kk) 196 ze=gauss3d1(3,kk) 197 weight=weight3d1(kk) 198 elseif(lakonl(4:7).eq.'20RB') then 199 if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then 200 xi=gauss3d13(1,kk) 201 et=gauss3d13(2,kk) 202 ze=gauss3d13(3,kk) 203 weight=weight3d13(kk) 204 else 205 call beamintscheme(lakonl,mint3d,ielprop(nelem),prop, 206 & kk,xi,et,ze,weight) 207 endif 208 elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) 209 & then 210 xi=gauss3d2(1,kk) 211 et=gauss3d2(2,kk) 212 ze=gauss3d2(3,kk) 213 weight=weight3d2(kk) 214 elseif(lakonl(4:4).eq.'2') then 215 xi=gauss3d3(1,kk) 216 et=gauss3d3(2,kk) 217 ze=gauss3d3(3,kk) 218 weight=weight3d3(kk) 219 elseif(lakonl(4:5).eq.'10') then 220 xi=gauss3d5(1,kk) 221 et=gauss3d5(2,kk) 222 ze=gauss3d5(3,kk) 223 weight=weight3d5(kk) 224 elseif(lakonl(4:4).eq.'4') then 225 xi=gauss3d4(1,kk) 226 et=gauss3d4(2,kk) 227 ze=gauss3d4(3,kk) 228 weight=weight3d4(kk) 229 elseif(lakonl(4:5).eq.'15') then 230 xi=gauss3d8(1,kk) 231 et=gauss3d8(2,kk) 232 ze=gauss3d8(3,kk) 233 weight=weight3d8(kk) 234 else 235 xi=gauss3d7(1,kk) 236 et=gauss3d7(2,kk) 237 ze=gauss3d7(3,kk) 238 weight=weight3d7(kk) 239 endif 240! 241! calculation of the shape functions and their derivatives 242! in the gauss point 243! 244 if(nope.eq.20) then 245 if(lakonl(7:7).eq.'A') then 246 call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag) 247 elseif((lakonl(7:7).eq.'E').or.(lakonl(7:7).eq.'S')) then 248 call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag) 249 else 250 call shape20h(xi,et,ze,xl,xsj,shp,iflag) 251 endif 252 elseif(nope.eq.8) then 253 call shape8h(xi,et,ze,xl,xsj,shp,iflag) 254 elseif(nope.eq.10) then 255 call shape10tet(xi,et,ze,xl,xsj,shp,iflag) 256 elseif(nope.eq.4) then 257 call shape4tet(xi,et,ze,xl,xsj,shp,iflag) 258 elseif(nope.eq.15) then 259 call shape15w(xi,et,ze,xl,xsj,shp,iflag) 260 else 261 call shape6w(xi,et,ze,xl,xsj,shp,iflag) 262 endif 263! 264! check the jacobian determinant 265! 266 if(xsj.lt.1.d-20) then 267 write(*,*) '*ERROR in e_c3d_rhs: nonpositive jacobian' 268 write(*,*) ' determinant in element',nelem 269 write(*,*) 270 xsj=dabs(xsj) 271 nmethod=0 272 endif 273! 274! computation of the right hand side 275! 276! incorporating the jacobian determinant in the shape 277! functions 278! 279 if((nload.gt.0).or.(nbody.ne.0)) then 280 do i1=1,nope 281 shpj(4,i1)=shp(4,i1)*xsj 282 enddo 283 endif 284! 285! distributed body flux 286! 287 if(nload.gt.0) then 288 call nident2(nelemload,nelem,nload,id) 289 do 290 if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit 291 if((sideload(id)(1:2).ne.'BX').and. 292 & (sideload(id)(1:2).ne.'BY').and. 293 & (sideload(id)(1:2).ne.'BZ')) then 294 id=id-1 295 cycle 296 endif 297 if(sideload(id)(3:4).eq.'NU') then 298 do j=1,3 299 coords(j)=0.d0 300 do i1=1,nope 301 coords(j)=coords(j)+ 302 & shp(4,i1)*xl(j,i1) 303 enddo 304 enddo 305 if(sideload(id)(1:2).eq.'BX') then 306 jltyp=1 307 elseif(sideload(id)(1:2).eq.'BY') then 308 jltyp=2 309 elseif(sideload(id)(1:2).eq.'BZ') then 310 jltyp=3 311 endif 312 iscale=1 313 call dload(xload(1,id),istep,iinc,tvar,nelem,i, 314 & layer,kspt,coords,jltyp,sideload(id),vold,co, 315 & lakonl,konl,ipompc,nodempc,coefmpc,nmpc,ikmpc, 316 & ilmpc,iscale,veold,rho,amat,mi) 317 if((nmethod.eq.1).and.(iscale.ne.0)) 318 & xload(1,id)=xloadold(1,id)+ 319 & (xload(1,id)-xloadold(1,id))*reltime 320 endif 321 jj1=1 322 do jj=1,nope 323 if(sideload(id)(1:2).eq.'BX') 324 & ff(jj1)=ff(jj1)+xload(1,id)*shpj(4,jj)*weight 325 if(sideload(id)(1:2).eq.'BY') 326 & ff(jj1+1)=ff(jj1+1)+xload(1,id)*shpj(4,jj) 327 & *weight 328 if(sideload(id)(1:2).eq.'BZ') 329 & ff(jj1+2)=ff(jj1+2)+xload(1,id)*shpj(4,jj) 330 & *weight 331 jj1=jj1+3 332 enddo 333 id=id-1 334 enddo 335 endif 336! 337! body forces 338! 339 if(nbody.ne.0) then 340! 341 do i1=1,3 342 bf(i1)=0.d0 343 enddo 344 do jj1=1,2 345 if(dabs(om(jj1)).lt.1.d-20) cycle 346 do i1=1,3 347! 348! computation of the global coordinates of the gauss 349! point 350! 351 q(i1)=0.d0 352c if(iperturb.eq.0) then 353 if((iperturb(1).ne.1).and.(iperturb(2).ne.1)) 354 & then 355 do j1=1,nope 356 q(i1)=q(i1)+shp(4,j1)*xl(i1,j1) 357 enddo 358 else 359 do j1=1,nope 360 q(i1)=q(i1)+shp(4,j1)* 361 & (xl(i1,j1)+voldl(i1,j1)) 362 enddo 363 endif 364! 365 q(i1)=q(i1)-p1(i1,jj1) 366 enddo 367 const=q(1)*p2(1,jj1)+q(2)*p2(2,jj1)+q(3)*p2(3,jj1) 368! 369! inclusion of the centrifugal force into the body force 370! 371 do i1=1,3 372 bf(i1)=bf(i1)+(q(i1)-const*p2(i1,jj1))*om(jj1) 373 enddo 374 enddo 375! 376 do i1=1,3 377 bf(i1)=bf(i1)+bodyf(i1) 378 enddo 379! 380 jj1=1 381 do jj=1,nope 382 ff(jj1)=ff(jj1)+bf(1)*shpj(4,jj)*weight 383 ff(jj1+1)=ff(jj1+1)+bf(2)*shpj(4,jj)*weight 384 ff(jj1+2)=ff(jj1+2)+bf(3)*shpj(4,jj)*weight 385 jj1=jj1+3 386 enddo 387 endif 388! 389 enddo 390 endif 391! 392! distributed loads 393! 394 if(nload.eq.0) then 395 return 396 endif 397 call nident2(nelemload,nelem,nload,id) 398 do 399 if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit 400 if(sideload(id)(1:1).ne.'P') then 401 id=id-1 402 cycle 403 endif 404 read(sideload(id)(2:2),'(i1)') ig 405! 406! treatment of wedge faces 407! 408 if(lakonl(4:4).eq.'6') then 409 mint2d=1 410 if(ig.le.2) then 411 nopes=3 412 else 413 nopes=4 414 endif 415 endif 416 if(lakonl(4:5).eq.'15') then 417 if(ig.le.2) then 418 mint2d=3 419 nopes=6 420 else 421 mint2d=4 422 nopes=8 423 endif 424 endif 425! 426 if((nope.eq.20).or.(nope.eq.8)) then 427c if(iperturb.eq.0) then 428 if((iperturb(1).ne.1).and.(iperturb(2).ne.1)) then 429 do i=1,nopes 430 do j=1,3 431 xl2(j,i)=co(j,konl(ifaceq(i,ig))) 432 enddo 433 enddo 434 else 435 do i=1,nopes 436 do j=1,3 437 xl2(j,i)=co(j,konl(ifaceq(i,ig)))+ 438 & vold(j,konl(ifaceq(i,ig))) 439 enddo 440 enddo 441 endif 442 elseif((nope.eq.10).or.(nope.eq.4)) then 443c if(iperturb.eq.0) then 444 if((iperturb(1).ne.1).and.(iperturb(2).ne.1)) then 445 do i=1,nopes 446 do j=1,3 447 xl2(j,i)=co(j,konl(ifacet(i,ig))) 448 enddo 449 enddo 450 else 451 do i=1,nopes 452 do j=1,3 453 xl2(j,i)=co(j,konl(ifacet(i,ig)))+ 454 & vold(j,konl(ifacet(i,ig))) 455 enddo 456 enddo 457 endif 458 else 459c if(iperturb.eq.0) then 460 if((iperturb(1).ne.1).and.(iperturb(2).ne.1)) then 461 do i=1,nopes 462 do j=1,3 463 xl2(j,i)=co(j,konl(ifacew(i,ig))) 464 enddo 465 enddo 466 else 467 do i=1,nopes 468 do j=1,3 469 xl2(j,i)=co(j,konl(ifacew(i,ig)))+ 470 & vold(j,konl(ifacew(i,ig))) 471 enddo 472 enddo 473 endif 474 endif 475! 476 do i=1,mint2d 477 if((lakonl(4:5).eq.'8R').or. 478 & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then 479 xi=gauss2d1(1,i) 480 et=gauss2d1(2,i) 481 weight=weight2d1(i) 482 elseif((lakonl(4:4).eq.'8').or. 483 & (lakonl(4:6).eq.'20R').or. 484 & ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then 485 xi=gauss2d2(1,i) 486 et=gauss2d2(2,i) 487 weight=weight2d2(i) 488 elseif(lakonl(4:4).eq.'2') then 489 xi=gauss2d3(1,i) 490 et=gauss2d3(2,i) 491 weight=weight2d3(i) 492 elseif((lakonl(4:5).eq.'10').or. 493 & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then 494 xi=gauss2d5(1,i) 495 et=gauss2d5(2,i) 496 weight=weight2d5(i) 497 elseif((lakonl(4:4).eq.'4').or. 498 & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then 499 xi=gauss2d4(1,i) 500 et=gauss2d4(2,i) 501 weight=weight2d4(i) 502 endif 503! 504 if(nopes.eq.8) then 505 call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag) 506 elseif(nopes.eq.4) then 507 call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag) 508 elseif(nopes.eq.6) then 509 call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 510 else 511 call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 512 endif 513! 514! for nonuniform load: determine the coordinates of the 515! point (transferred into the user subroutine) 516! 517 if(sideload(id)(3:4).eq.'NU') then 518 do k=1,3 519 coords(k)=0.d0 520 do j=1,nopes 521 coords(k)=coords(k)+xl2(k,j)*shp2(4,j) 522 enddo 523 enddo 524 read(sideload(id)(2:2),'(i1)') jltyp 525 jltyp=jltyp+20 526 iscale=1 527 call dload(xload(1,id),istep,iinc,tvar,nelem,i,layer, 528 & kspt,coords,jltyp,sideload(id),vold,co,lakonl, 529 & konl,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc, 530 & iscale,veold,rho,amat,mi) 531 if((nmethod.eq.1).and.(iscale.ne.0)) 532 & xload(1,id)=xloadold(1,id)+ 533 & (xload(1,id)-xloadold(1,id))*reltime 534 endif 535! 536 do k=1,nopes 537 if((nope.eq.20).or.(nope.eq.8)) then 538 ipointer=(ifaceq(k,ig)-1)*3 539 elseif((nope.eq.10).or.(nope.eq.4)) then 540 ipointer=(ifacet(k,ig)-1)*3 541 else 542 ipointer=(ifacew(k,ig)-1)*3 543 endif 544 ff(ipointer+1)=ff(ipointer+1)-shp2(4,k)*xload(1,id) 545 & *xsj2(1)*weight 546 ff(ipointer+2)=ff(ipointer+2)-shp2(4,k)*xload(1,id) 547 & *xsj2(2)*weight 548 ff(ipointer+3)=ff(ipointer+3)-shp2(4,k)*xload(1,id) 549 & *xsj2(3)*weight 550 enddo 551 enddo 552! 553 id=id-1 554 enddo 555! 556 return 557 end 558 559