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_th(co,nk,konl,lakonl, 20 & ff,nelem,nmethod,t0,t1,vold,nelemload, 21 & sideload,xload,nload,idist,dtime, 22 & ttime,time,istep,iinc,xloadold,reltime, 23 & ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,mi, 24 & ielprop,prop,sti,xstateini,xstate,nstate_) 25! 26! computation of the rhs for the element with 27! the topology in konl 28! 29! ff: rhs without temperature and eigenstress contribution 30! 31 implicit none 32! 33 logical ivolumeforce 34! 35 character*8 lakonl 36 character*20 sideload(*) 37! 38 integer konl(20),ifaceq(8,6),nelemload(2,*),nk,nelem,nmethod, 39 & nload,idist,i,j,k,i1,iflag,ipompc(*),nodempc(3,*),nmpc, 40 & jj,id,ipointer,ig,kk,nope,nopes,mint2d,ikmpc(*),ilmpc(*), 41 & mint3d,ifacet(6,4),nopev,ifacew(8,5),iinc,istep,jltyp, 42 & iscale,mi(*),ielprop(*),null,nstate_ 43! 44 real*8 co(3,*),xl(3,20),shp(4,20),xs2(3,7),xloadold(2,*), 45 & ff(60),shpj(4,20),dxsj2,temp,press,t0(*),t1(*),coords(3), 46 & xl2(3,8),xsj2(3),shp2(7,8),vold(0:mi(2),*),xload(2,*), 47 & xi,et,ze,xsj,xsjj,t1l,ttime,time,weight,pgauss(3),tvar(2), 48 & reltime,areaj,coefmpc(*),tl2(8),prop(*),sti(6,mi(1),*), 49 & xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*) 50! 51 real*8 dtime 52! 53 include "gauss.f" 54! 55 data ifaceq /4,3,2,1,11,10,9,12, 56 & 5,6,7,8,13,14,15,16, 57 & 1,2,6,5,9,18,13,17, 58 & 2,3,7,6,10,19,14,18, 59 & 3,4,8,7,11,20,15,19, 60 & 4,1,5,8,12,17,16,20/ 61 data ifacet /1,3,2,7,6,5, 62 & 1,2,4,5,9,8, 63 & 2,3,4,6,10,9, 64 & 1,4,3,8,10,7/ 65 data ifacew /1,3,2,9,8,7,0,0, 66 & 4,5,6,10,11,12,0,0, 67 & 1,2,5,4,7,14,10,13, 68 & 2,3,6,5,8,15,11,14, 69 & 4,6,3,1,12,15,9,13/ 70 data iflag /3/ 71 data null /0/ 72! 73 tvar(1)=time 74 tvar(2)=ttime+time 75! 76 if(lakonl(4:4).eq.'2') then 77 nope=20 78 nopev=8 79 nopes=8 80 elseif(lakonl(4:4).eq.'8') then 81 nope=8 82 nopev=8 83 nopes=4 84 elseif(lakonl(4:5).eq.'10') then 85 nope=10 86 nopev=4 87 nopes=6 88 elseif(lakonl(4:4).eq.'4') then 89 nope=4 90 nopev=4 91 nopes=3 92 elseif(lakonl(4:5).eq.'15') then 93 nope=15 94 nopev=6 95 else 96 nope=6 97 nopev=6 98 endif 99! 100 if(lakonl(4:5).eq.'8R') then 101 mint2d=1 102 mint3d=1 103 elseif(lakonl(4:7).eq.'20RB') then 104 if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then 105 mint2d=4 106 mint3d=50 107 else 108 mint2d=4 109 call beamintscheme(lakonl,mint3d,ielprop(nelem),prop, 110 & null,xi,et,ze,weight) 111 endif 112 elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) then 113 mint2d=4 114 mint3d=8 115 elseif(lakonl(4:4).eq.'2') then 116 mint2d=9 117 mint3d=27 118 elseif(lakonl(4:5).eq.'10') then 119 mint2d=3 120 mint3d=4 121 elseif(lakonl(4:4).eq.'4') then 122 mint2d=1 123 mint3d=1 124 elseif(lakonl(4:5).eq.'15') then 125 mint3d=9 126 else 127 mint3d=2 128 endif 129! 130! computation of the coordinates of the local nodes 131! 132 do i=1,nope 133 do j=1,3 134 xl(j,i)=co(j,konl(i)) 135 enddo 136 enddo 137! 138! initialisation for distributed forces 139! 140c if(idist.ne.0) then 141 do i=1,nope 142 ff(i)=0.d0 143 enddo 144c endif 145! 146! computation of the body forces 147! 148 ivolumeforce=.false. 149c if(nload.gt.0) then 150 call nident2(nelemload,nelem,nload,id) 151 do 152 if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit 153 if(sideload(id)(1:2).ne.'BF') then 154 id=id-1 155 cycle 156 else 157 ivolumeforce=.true. 158 exit 159 endif 160 enddo 161c endif 162! 163! computation of the matrix: loop over the Gauss points 164! 165 if(ivolumeforce) then 166 do kk=1,mint3d 167 if(lakonl(4:5).eq.'8R') then 168 xi=gauss3d1(1,kk) 169 et=gauss3d1(2,kk) 170 ze=gauss3d1(3,kk) 171 weight=weight3d1(kk) 172 elseif(lakonl(4:7).eq.'20RB') then 173 if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then 174 xi=gauss3d13(1,kk) 175 et=gauss3d13(2,kk) 176 ze=gauss3d13(3,kk) 177 weight=weight3d13(kk) 178 else 179 call beamintscheme(lakonl,mint3d,ielprop(nelem),prop, 180 & kk,xi,et,ze,weight) 181 endif 182 elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) 183 & then 184 xi=gauss3d2(1,kk) 185 et=gauss3d2(2,kk) 186 ze=gauss3d2(3,kk) 187 weight=weight3d2(kk) 188 elseif(lakonl(4:4).eq.'2') then 189 xi=gauss3d3(1,kk) 190 et=gauss3d3(2,kk) 191 ze=gauss3d3(3,kk) 192 weight=weight3d3(kk) 193 elseif(lakonl(4:5).eq.'10') then 194 xi=gauss3d5(1,kk) 195 et=gauss3d5(2,kk) 196 ze=gauss3d5(3,kk) 197 weight=weight3d5(kk) 198 elseif(lakonl(4:4).eq.'4') then 199 xi=gauss3d4(1,kk) 200 et=gauss3d4(2,kk) 201 ze=gauss3d4(3,kk) 202 weight=weight3d4(kk) 203 elseif(lakonl(4:5).eq.'15') then 204 xi=gauss3d8(1,kk) 205 et=gauss3d8(2,kk) 206 ze=gauss3d8(3,kk) 207 weight=weight3d8(kk) 208 else 209 xi=gauss3d7(1,kk) 210 et=gauss3d7(2,kk) 211 ze=gauss3d7(3,kk) 212 weight=weight3d7(kk) 213 endif 214! 215! calculation of the shape functions and their derivatives 216! in the gauss point 217! 218 if(nope.eq.20) then 219 if(lakonl(7:7).eq.'A') then 220 call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag) 221 elseif((lakonl(7:7).eq.'E').or.(lakonl(7:7).eq.'S')) then 222 call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag) 223 else 224 call shape20h(xi,et,ze,xl,xsj,shp,iflag) 225 endif 226 elseif(nope.eq.8) then 227 call shape8h(xi,et,ze,xl,xsj,shp,iflag) 228 elseif(nope.eq.10) then 229 call shape10tet(xi,et,ze,xl,xsj,shp,iflag) 230 elseif(nope.eq.4) then 231 call shape4tet(xi,et,ze,xl,xsj,shp,iflag) 232 elseif(nope.eq.15) then 233 call shape15w(xi,et,ze,xl,xsj,shp,iflag) 234 else 235 call shape6w(xi,et,ze,xl,xsj,shp,iflag) 236 endif 237! 238! check the jacobian determinant 239! 240 if(xsj.lt.1.d-20) then 241 write(*,*) '*ERROR in e_c3d_rhs_th: nonpositive jacobian' 242 write(*,*) ' determinant in element',nelem 243 write(*,*) 244 xsj=dabs(xsj) 245 nmethod=0 246 endif 247! 248! calculating the temperature in the integration 249! point 250! 251 t1l=0.d0 252! 253 do i1=1,nope 254 t1l=t1l+shp(4,i1)*vold(0,konl(i1)) 255 enddo 256! 257! incorporating the jacobian determinant in the shape 258! functions 259! 260 xsjj=dsqrt(xsj) 261 do i1=1,nope 262 shpj(1,i1)=shp(1,i1)*xsjj 263 shpj(2,i1)=shp(2,i1)*xsjj 264 shpj(3,i1)=shp(3,i1)*xsjj 265 shpj(4,i1)=shp(4,i1)*xsj 266 enddo 267! 268! computation of the right hand side 269! 270! distributed heat flux 271! 272c if(nload.gt.0) then 273 call nident2(nelemload,nelem,nload,id) 274 areaj=xsj*weight 275 do 276 if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit 277 if(sideload(id)(1:2).ne.'BF') then 278 id=id-1 279 cycle 280 endif 281 if(sideload(id)(3:4).eq.'NU') then 282 do j=1,3 283 pgauss(j)=0.d0 284 do i1=1,nope 285 pgauss(j)=pgauss(j)+ 286 & shp(4,i1)*co(j,konl(i1)) 287 enddo 288 enddo 289 jltyp=1 290 iscale=1 291 call dflux(xload(1,id),t1l,istep,iinc,tvar, 292 & nelem,kk,pgauss,jltyp,temp,press,sideload(id), 293 & areaj,vold,co,lakonl,konl,ipompc,nodempc,coefmpc, 294 & nmpc,ikmpc,ilmpc,iscale,mi) 295 if((nmethod.eq.1).and.(iscale.ne.0)) 296 & xload(1,id)=xloadold(1,id)+ 297 & (xload(1,id)-xloadold(1,id))*reltime 298 endif 299 do jj=1,nope 300 ff(jj)=ff(jj)+xload(1,id)*shpj(4,jj)*weight 301 enddo 302 exit 303 enddo 304c endif 305! 306 enddo 307 endif 308! 309! distributed loads 310! 311c if(nload.eq.0) then 312c return 313c endif 314 call nident2(nelemload,nelem,nload,id) 315 do 316 if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit 317 if((sideload(id)(1:1).ne.'F').and. 318 & (sideload(id)(1:1).ne.'R').and. 319 & (sideload(id)(1:1).ne.'S')) then 320 id=id-1 321 cycle 322 endif 323 read(sideload(id)(2:2),'(i1)') ig 324! 325! treatment of wedge faces 326! 327 if(lakonl(4:4).eq.'6') then 328 mint2d=1 329 if(ig.le.2) then 330 nopes=3 331 else 332 nopes=4 333 endif 334 endif 335 if(lakonl(4:5).eq.'15') then 336 if(ig.le.2) then 337 mint2d=3 338 nopes=6 339 else 340 mint2d=4 341 nopes=8 342 endif 343 endif 344! 345 if((nope.eq.20).or.(nope.eq.8)) then 346 do i=1,nopes 347 tl2(i)=vold(0,konl(ifaceq(i,ig))) 348 do j=1,3 349 xl2(j,i)=co(j,konl(ifaceq(i,ig)))+ 350 & vold(j,konl(ifaceq(i,ig))) 351 enddo 352 enddo 353 elseif((nope.eq.10).or.(nope.eq.4)) then 354 do i=1,nopes 355 tl2(i)=vold(0,konl(ifacet(i,ig))) 356 do j=1,3 357 xl2(j,i)=co(j,konl(ifacet(i,ig)))+ 358 & vold(j,konl(ifacet(i,ig))) 359 enddo 360 enddo 361 else 362 do i=1,nopes 363 tl2(i)=vold(0,konl(ifacew(i,ig))) 364 do j=1,3 365 xl2(j,i)=co(j,konl(ifacew(i,ig)))+ 366 & vold(j,konl(ifacew(i,ig))) 367 enddo 368 enddo 369 endif 370! 371 do i=1,mint2d 372 if((lakonl(4:5).eq.'8R').or. 373 & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then 374 xi=gauss2d1(1,i) 375 et=gauss2d1(2,i) 376 weight=weight2d1(i) 377 elseif((lakonl(4:4).eq.'8').or. 378 & (lakonl(4:6).eq.'20R').or. 379 & ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then 380 xi=gauss2d2(1,i) 381 et=gauss2d2(2,i) 382 weight=weight2d2(i) 383 elseif(lakonl(4:4).eq.'2') then 384 xi=gauss2d3(1,i) 385 et=gauss2d3(2,i) 386 weight=weight2d3(i) 387 elseif((lakonl(4:5).eq.'10').or. 388 & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then 389 xi=gauss2d5(1,i) 390 et=gauss2d5(2,i) 391 weight=weight2d5(i) 392 elseif((lakonl(4:4).eq.'4').or. 393 & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then 394 xi=gauss2d4(1,i) 395 et=gauss2d4(2,i) 396 weight=weight2d4(i) 397 endif 398! 399 if(nopes.eq.8) then 400 call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag) 401 elseif(nopes.eq.4) then 402 call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag) 403 elseif(nopes.eq.6) then 404 call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 405 else 406 call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 407 endif 408! 409 dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+ 410 & xsj2(3)*xsj2(3)) 411 areaj=dxsj2*weight 412! 413 temp=0.d0 414 do j=1,nopes 415 temp=temp+tl2(j)*shp2(4,j) 416 enddo 417! 418! for nonuniform load: determine the coordinates of the 419! point (transferred into the user subroutine) 420! 421 if(sideload(id)(3:4).eq.'NU') then 422 do k=1,3 423 coords(k)=0.d0 424 do j=1,nopes 425 coords(k)=coords(k)+xl2(k,j)*shp2(4,j) 426 enddo 427 enddo 428 read(sideload(id)(2:2),'(i1)') jltyp 429 jltyp=jltyp+10 430 if(sideload(id)(1:1).eq.'S') then 431 iscale=1 432 call dflux(xload(1,id),temp,istep,iinc,tvar, 433 & nelem,i,coords,jltyp,temp,press,sideload(id), 434 & areaj,vold,co,lakonl,konl,ipompc,nodempc, 435 & coefmpc,nmpc,ikmpc,ilmpc,iscale,mi) 436 if((nmethod.eq.1).and.(iscale.ne.0)) 437 & xload(1,id)=xloadold(1,id)+ 438 & (xload(1,id)-xloadold(1,id))*reltime 439 endif 440 endif 441! 442 do k=1,nopes 443 if((nope.eq.20).or.(nope.eq.8)) then 444 ipointer=ifaceq(k,ig) 445 elseif((nope.eq.10).or.(nope.eq.4)) then 446 ipointer=ifacet(k,ig) 447 else 448 ipointer=ifacew(k,ig) 449 endif 450 if(sideload(id)(1:1).eq.'S') then 451! 452! flux INTO the face is positive (input deck convention) 453! this is different from the convention in the theory 454! 455 ff(ipointer)=ff(ipointer)+shp2(4,k)*xload(1,id) 456 & *dxsj2*weight 457 elseif(sideload(id)(1:1).eq.'F') then 458 write(*,*) '*ERROR in e_c3d_rhs_th.f: no' 459 write(*,*) ' film conditions allowed' 460 write(*,*) ' in an modal dynamic calculation' 461 call exit(201) 462 elseif(sideload(id)(1:1).eq.'R') then 463 write(*,*) '*ERROR in e_c3d_rhs_th.f: no' 464 write(*,*) ' radiation conditions allowed' 465 write(*,*) ' in an modal dynamic calculation' 466 call exit(201) 467 endif 468 enddo 469! 470 enddo 471! 472 id=id-1 473 enddo 474! 475 return 476 end 477 478 479 480