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 rhs(co,nk,kon,ipkon,lakon,ne, 20 & ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc, 21 & nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr, 22 & fext,nactdof,neq,nmethod, 23 & ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,alcon, 24 & nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,t0,t1, 25 & ithermal,iprestr,vold,iperturb,iexpl,plicon, 26 & nplicon,plkcon,nplkcon,npmat_,ttime,time,istep,iinc,dtime, 27 & physcon,ibody,xloadold,reltime,veold,matname,mi,ikactmech, 28 & nactmech,ielprop,prop,sti,xstateini,xstate,nstate_, 29 & ntrans,inotr,trab,nea,neb) 30! 31! filling the right hand side load vector b 32! 33! b contains the contributions due to mechanical forces only 34! 35 implicit none 36! 37 character*8 lakon(*) 38 character*20 sideload(*) 39 character*80 matname(*) 40! 41 integer kon(*),ipompc(*),nodempc(3,*),ipobody(2,*),nbody, 42 & nodeforc(2,*),ndirforc(*),nelemload(2,*),ikmpc(*),mi(*), 43 & ilmpc(*),nactdof(0:mi(2),*),konl(20),nelcon(2,*),ibody(3,*), 44 & nrhcon(*),nalcon(2,*),ielmat(mi(3),*),ielorien(mi(3),*), 45 & ipkon(*),ielprop(*),nstate_,idummy,rhsi,ntrans,inotr(2,*), 46 & nk,ne,nmpc,nforc,nload,neq,nmethod,nom,m,idm,idir,itr, 47 & ithermal(*),iprestr,iperturb(*),i,j,k,idist,jj,node, 48 & id,ist,index,jdof1,jdof,node1,ntmat_,indexe,nope,norien, 49 & iexpl,idof1,iinc,istep,icalccg,nplicon(0:ntmat_,*), 50 & nplkcon(0:ntmat_,*),npmat_,ikactmech(*),nactmech,nea,neb 51! 52 real*8 co(3,*),coefmpc(*),xforc(*),xload(2,*),p1(3,2), 53 & p2(3,2),fext(*),bodyf(3),elcon(0:21,ntmat_,*),a(3,3), 54 & rhcon(0:1,ntmat_,*),xloadold(2,*),reltime,prop(*), 55 & alcon(0:6,ntmat_,*),alzero(*),orab(7,*),xbody(7,*),cgr(4,*), 56 & t0(*),t1(*),vold(0:mi(2),*),ff(60),time,ttime,dtime, 57 & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*), 58 & om(2),physcon(*),veold(0:mi(2),*),sti(6,mi(1),*), 59 & xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*), 60 & fnext,trab(7,*) 61! 62 icalccg=0 63! 64c if((nmethod.ge.4).and.(iperturb(1).lt.2).and.(nactmech.lt.neq/2)) 65c & then 66c ! 67c ! modal dynamics and steady state dynamics: 68c ! only nonzeros are reset to zero 69c ! 70c do i=1,nactmech 71c fext(ikactmech(i)+1)=0.d0 72c enddo 73c else 74c do i=1,neq 75c fext(i)=0.d0 76c enddo 77c endif 78 nactmech=0 79! 80! distributed forces (body forces or thermal loads or 81! residual stresses or distributed face loads) 82! 83 if((nbody.ne.0).or.(nload.ne.0)) then 84 idist=1 85 else 86 idist=0 87 endif 88! 89 if(((ithermal(1).le.1).or.(ithermal(1).eq.3)).and.(idist.ne.0)) 90 & then 91! 92! mechanical analysis: loop over all elements 93! 94 do i=nea,neb 95! 96 if(ipkon(i).lt.0) cycle 97 indexe=ipkon(i) 98 if(lakon(i)(4:4).eq.'2') then 99 nope=20 100 elseif(lakon(i)(4:4).eq.'8') then 101 nope=8 102 elseif(lakon(i)(4:5).eq.'10') then 103 nope=10 104 elseif(lakon(i)(4:4).eq.'4') then 105 nope=4 106 elseif(lakon(i)(4:5).eq.'15') then 107 nope=15 108 elseif(lakon(i)(4:4).eq.'6') then 109 nope=6 110 else 111 cycle 112 endif 113! 114 do j=1,nope 115 konl(j)=kon(indexe+j) 116 enddo 117! 118! assigning centrifugal forces 119! 120 if(nbody.gt.0) then 121 nom=0 122 om(1)=0.d0 123 om(2)=0.d0 124 bodyf(1)=0.d0 125 bodyf(2)=0.d0 126 bodyf(3)=0.d0 127! 128 index=i 129 do 130 j=ipobody(1,index) 131 if(j.eq.0) exit 132 if(ibody(1,j).eq.1) then 133 nom=nom+1 134 if(nom.gt.2) then 135 write(*,*) 136 & '*ERROR in rhs: no more than two centri-' 137 write(*,*) 138 & ' fugal loading cards allowed' 139 call exit(201) 140 endif 141 om(nom)=xbody(1,j) 142 p1(1,nom)=xbody(2,j) 143 p1(2,nom)=xbody(3,j) 144 p1(3,nom)=xbody(4,j) 145 p2(1,nom)=xbody(5,j) 146 p2(2,nom)=xbody(6,j) 147 p2(3,nom)=xbody(7,j) 148! 149! assigning gravity forces 150! 151 elseif(ibody(1,j).eq.2) then 152 bodyf(1)=bodyf(1)+xbody(1,j)*xbody(2,j) 153 bodyf(2)=bodyf(2)+xbody(1,j)*xbody(3,j) 154 bodyf(3)=bodyf(3)+xbody(1,j)*xbody(4,j) 155! 156! assigning newton gravity forces 157! 158 elseif(ibody(1,j).eq.3) then 159 call newton(icalccg,ne,ipkon,lakon,kon,t0,co,rhcon, 160 & nrhcon,ntmat_,physcon,i,cgr,bodyf,ielmat, 161 & ithermal,vold) 162 endif 163 index=ipobody(2,index) 164 if(index.eq.0) exit 165 enddo 166 endif 167! 168 if(idist.ne.0) 169 & call e_c3d_rhs(co,nk,konl,lakon(i),p1,p2,om,bodyf,nbody, 170 & ff,i,nmethod,rhcon,ielmat,ntmat_,vold,iperturb, 171 & nelemload,sideload,xload,nload,idist,ttime,time,istep, 172 & iinc,dtime,xloadold,reltime,ipompc,nodempc,coefmpc,nmpc, 173 & ikmpc,ilmpc,veold,matname,mi,ielprop,prop) 174! 175! modal dynamics and steady state dynamics: 176! location of nonzeros is stored 177! 178 if((nmethod.ge.4).and.(iperturb(1).lt.2)) then 179 do jj=1,3*nope 180! 181 j=(jj-1)/3+1 182 k=jj-3*(j-1) 183! 184 node1=kon(indexe+j) 185 jdof1=nactdof(k,node1) 186! 187! distributed forces 188! 189 if(idist.ne.0) then 190 if(dabs(ff(jj)).lt.1.d-30) cycle 191 if(jdof1.le.0) then 192 if(nmpc.ne.0) then 193 idof1=(node1-1)*8+k 194 call nident(ikmpc,idof1,nmpc,id) 195 if((id.gt.0).and.(ikmpc(id).eq.idof1)) then 196 id=ilmpc(id) 197 ist=ipompc(id) 198 index=nodempc(3,ist) 199 do 200 jdof1=nactdof(nodempc(2,index), 201 & nodempc(1,index)) 202 if(jdof1.gt.0) then 203 fext(jdof1)=fext(jdof1) 204 & -coefmpc(index)*ff(jj)/ 205 & coefmpc(ist) 206 call nident(ikactmech,jdof1-1, 207 & nactmech,idm) 208 do 209 if(idm.gt.0) then 210 if(ikactmech(idm).eq.jdof1-1) 211 & exit 212 endif 213 nactmech=nactmech+1 214 do m=nactmech,idm+2,-1 215 ikactmech(m)=ikactmech(m-1) 216 enddo 217 ikactmech(idm+1)=jdof1-1 218 exit 219 enddo 220 endif 221 index=nodempc(3,index) 222 if(index.eq.0) exit 223 enddo 224 endif 225 endif 226 cycle 227 endif 228 fext(jdof1)=fext(jdof1)+ff(jj) 229 call nident(ikactmech,jdof1-1,nactmech, 230 & idm) 231 do 232 if(idm.gt.0) then 233 if(ikactmech(idm).eq.jdof1-1) exit 234 endif 235 nactmech=nactmech+1 236 do m=nactmech,idm+2,-1 237 ikactmech(m)=ikactmech(m-1) 238 enddo 239 ikactmech(idm+1)=jdof1-1 240 exit 241 enddo 242 endif 243! 244 enddo 245! 246! other procedures 247! 248 else 249 do jj=1,3*nope 250! 251 j=(jj-1)/3+1 252 k=jj-3*(j-1) 253! 254 node1=kon(indexe+j) 255 jdof1=nactdof(k,node1) 256! 257! distributed forces 258! 259 if(idist.ne.0) then 260 if(jdof1.le.0) then 261 if(nmpc.ne.0) then 262 idof1=(node1-1)*8+k 263 call nident(ikmpc,idof1,nmpc,id) 264 if((id.gt.0).and.(ikmpc(id).eq.idof1)) then 265 id=ilmpc(id) 266 ist=ipompc(id) 267 index=nodempc(3,ist) 268 do 269 jdof1=nactdof(nodempc(2,index), 270 & nodempc(1,index)) 271 if(jdof1.gt.0) then 272 fext(jdof1)=fext(jdof1) 273 & -coefmpc(index)*ff(jj)/ 274 & coefmpc(ist) 275 endif 276 index=nodempc(3,index) 277 if(index.eq.0) exit 278 enddo 279 endif 280 endif 281 cycle 282 endif 283 fext(jdof1)=fext(jdof1)+ff(jj) 284 endif 285! 286 enddo 287 endif 288 enddo 289! 290 elseif((ithermal(1).eq.2).and.(nload.gt.0)) then 291! 292! thermal analysis: loop over all elements 293! 294 do i=nea,neb 295! 296 if(ipkon(i).lt.0) cycle 297 indexe=ipkon(i) 298 if(lakon(i)(4:4).eq.'2') then 299 nope=20 300 elseif(lakon(i)(4:4).eq.'8') then 301 nope=8 302 elseif(lakon(i)(4:5).eq.'10') then 303 nope=10 304 elseif(lakon(i)(4:4).eq.'4') then 305 nope=4 306 elseif(lakon(i)(4:5).eq.'15') then 307 nope=15 308 elseif(lakon(i)(4:4).eq.'6') then 309 nope=6 310 else 311 cycle 312 endif 313! 314 do j=1,nope 315 konl(j)=kon(indexe+j) 316 enddo 317! 318 if(nload.gt.0) 319 & call e_c3d_rhs_th(co,nk,konl,lakon(i), 320 & ff,i,nmethod,t0,t1,vold,nelemload, 321 & sideload,xload,nload,idist,dtime, 322 & ttime,time,istep,iinc,xloadold,reltime, 323 & ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,mi, 324 & ielprop,prop,sti,xstateini,xstate,nstate_) 325! 326! modal dynamics and steady state dynamics: 327! location of nonzeros is stored 328! 329 if((nmethod.ge.4.and.(iperturb(1).lt.2))) then 330 do jj=1,nope 331! 332 j=jj 333! 334 node1=kon(indexe+j) 335 jdof1=nactdof(0,node1) 336! 337! distributed forces 338! 339 if(idist.ne.0) then 340 if(dabs(ff(jj)).lt.1.d-30) cycle 341 if(jdof1.le.0) then 342 if(nmpc.ne.0) then 343 idof1=(node1-1)*8 344 call nident(ikmpc,idof1,nmpc,id) 345 if((id.gt.0).and.(ikmpc(id).eq.idof1)) then 346 id=ilmpc(id) 347 ist=ipompc(id) 348 index=nodempc(3,ist) 349 do 350 jdof1=nactdof(nodempc(2,index), 351 & nodempc(1,index)) 352 if(jdof1.gt.0) then 353 fext(jdof1)=fext(jdof1) 354 & -coefmpc(index)*ff(jj)/ 355 & coefmpc(ist) 356 call nident(ikactmech,jdof1-1, 357 & nactmech,idm) 358 do 359 if(idm.gt.0) then 360 if(ikactmech(idm).eq.jdof1-1) 361 & exit 362 endif 363 nactmech=nactmech+1 364 do m=nactmech,idm+2,-1 365 ikactmech(m)=ikactmech(m-1) 366 enddo 367 ikactmech(idm+1)=jdof1-1 368 exit 369 enddo 370 endif 371 index=nodempc(3,index) 372 if(index.eq.0) exit 373 enddo 374 endif 375 endif 376 cycle 377 endif 378 fext(jdof1)=fext(jdof1)+ff(jj) 379 call nident(ikactmech,jdof1-1,nactmech, 380 & idm) 381 do 382 if(idm.gt.0) then 383 if(ikactmech(idm).eq.jdof1-1) exit 384 endif 385 nactmech=nactmech+1 386 do m=nactmech,idm+2,-1 387 ikactmech(m)=ikactmech(m-1) 388 enddo 389 ikactmech(idm+1)=jdof1-1 390 exit 391 enddo 392 endif 393! 394 enddo 395! 396! other procedures 397! 398 else 399 do jj=1,nope 400! 401 j=jj 402! 403 node1=kon(indexe+j) 404 jdof1=nactdof(0,node1) 405! 406! distributed forces 407! 408 if(idist.ne.0) then 409 if(jdof1.le.0) then 410 if(nmpc.ne.0) then 411 idof1=(node1-1)*8 412 call nident(ikmpc,idof1,nmpc,id) 413 if((id.gt.0).and.(ikmpc(id).eq.idof1)) then 414 id=ilmpc(id) 415 ist=ipompc(id) 416 index=nodempc(3,ist) 417 do 418 jdof1=nactdof(nodempc(2,index), 419 & nodempc(1,index)) 420 if(jdof1.gt.0) then 421 fext(jdof1)=fext(jdof1) 422 & -coefmpc(index)*ff(jj)/ 423 & coefmpc(ist) 424 endif 425 index=nodempc(3,index) 426 if(index.eq.0) exit 427 enddo 428 endif 429 endif 430 cycle 431 endif 432 fext(jdof1)=fext(jdof1)+ff(jj) 433 endif 434! 435 enddo 436 endif 437 enddo 438! 439 endif 440 441 return 442 end 443