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 temploadfem(xforcold,xforc,xforcact,iamforc,nforc, 20 & xloadold,xload,xloadact,iamload,nload,ibody,xbody,nbody, 21 & xbodyold,xbodyact,t1old,t1,t1act,iamt1,nk, 22 & amta,namta,nam,ampli,time,reltime,ttime,dtime,ithermal,nmethod, 23 & xbounold,xboun,xbounact,iamboun,nboun, 24 & nodeboun,ndirboun,nodeforc,ndirforc,istep,iinc, 25 & co,vold,itg,ntg,amname,ikboun,ilboun,nelemload,sideload,mi, 26 & ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset, 27 & iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc, 28 & set,nset) 29! 30! calculates the loading at a given time 31! 32 implicit none 33! 34 logical gasnode 35! 36 character*1 entity 37 character*20 sideload(*) 38 character*80 amname(*) 39 character*81 tieset(3,*),set(*) 40! 41 integer iamforc(*),iamload(2,*),iamt1(*),nelemload(2,*), 42 & nam,i,istart,iend,id,nforc,nload,nk,namta(3,*),ithermal(*), 43 & nmethod,iamt1i,iamboun(*),nboun,iamforci,iambouni,nset, 44 & iamloadi1,iamloadi2,ibody(3,*),itg(*),ntg,idof,one, 45 & nbody,iambodyi,nodeboun(*),ndirboun(*),nodeforc(2,*), 46 & ndirforc(*),istep,iinc,msecpt,node,j,ikboun(*),ilboun(*), 47 & ipresboun,mi(*),ntrans,inotr(2,*),idummy,integerglob(*), 48 & istartset(*),iendset(*),ialset(*),ntie,iselect(1), 49 & nmpc,ikmpc(*),ilmpc(*),nodempc(3,*),k,ist,index,ipompc(*) 50! 51 real*8 xforc(*),xforcact(*),xload(2,*),xloadact(2,*), 52 & t1(*),t1act(*),amta(2,*),ampli(*),time,fixed_temp, 53 & xforcold(*),xloadold(2,*),t1old(*),reltime,coefmpc(*), 54 & xbounold(*),xboun(*),xbounact(*),ttime,dtime,reftime, 55 & xbody(7,*),xbodyold(7,*),xbodyact(7,*),co(3,*), 56 & vold(0:mi(2),*),abqtime(2),coords(3),trab(7,*), 57 & veold(0:mi(2),*),ddummy,doubleglob(*) 58! 59 data msecpt /1/ 60 data one /1/ 61! 62! if an amplitude is active, the loading is scaled according to 63! the actual time. If no amplitude is active, then the load is 64! - scaled according to the relative time for a static step 65! - applied as a step loading for a dynamic step 66! 67! calculating all amplitude values for the current time 68! 69 do i=1,nam 70 if(namta(3,i).lt.0) then 71 reftime=ttime+dtime 72 else 73 reftime=time 74 endif 75 if(abs(namta(3,i)).ne.i) then 76 reftime=reftime-amta(1,namta(1,i)) 77 istart=namta(1,abs(namta(3,i))) 78 iend=namta(2,abs(namta(3,i))) 79 if(istart.eq.0) then 80 call uamplitude(reftime,amname(namta(3,i)),ampli(i)) 81 cycle 82 endif 83 else 84 istart=namta(1,i) 85 iend=namta(2,i) 86 if(istart.eq.0) then 87 call uamplitude(reftime,amname(i),ampli(i)) 88 cycle 89 endif 90 endif 91 call identamta(amta,reftime,istart,iend,id) 92 if(id.lt.istart) then 93 ampli(i)=amta(2,istart) 94 elseif(id.eq.iend) then 95 ampli(i)=amta(2,iend) 96 else 97 ampli(i)=amta(2,id)+(amta(2,id+1)-amta(2,id)) 98 & *(reftime-amta(1,id))/(amta(1,id+1)-amta(1,id)) 99 endif 100 enddo 101! 102! scaling the boundary conditions 103! 104 do i=1,nboun 105 if((xboun(i).lt.1.2357111318d0).and. 106 & (xboun(i).gt.1.2357111316d0)) then 107! 108! user subroutine for boundary conditions 109! 110 node=nodeboun(i) 111! 112! check whether node is a gasnode 113! 114 gasnode=.false. 115 call nident(itg,node,ntg,id) 116 if(id.gt.0) then 117 if(itg(id).eq.node) then 118 gasnode=.true. 119 endif 120 endif 121! 122 abqtime(1)=time 123 abqtime(2)=ttime+dtime 124! 125! a gasnode cannot move (displacement DOFs are used 126! for other purposes, e.g. mass flow and pressure) 127! 128 if(gasnode) then 129 do j=1,3 130 coords(j)=co(j,node) 131 enddo 132 else 133 do j=1,3 134 coords(j)=co(j,node)+vold(j,node) 135 enddo 136 endif 137! 138 if(ndirboun(i).eq.0) then 139 call utemp(xbounact(i),msecpt,istep,iinc,abqtime,node, 140 & coords,vold,mi) 141 else 142 call uboun(xbounact(i),istep,iinc,abqtime,node, 143 & ndirboun(i),coords,vold,mi) 144 endif 145 cycle 146 endif 147 if((xboun(i).lt.1.9232931375d0).and. 148 & (xboun(i).gt.1.9232931373d0)) then 149! 150! boundary conditions for submodel 151! 152 node=nodeboun(i) 153! 154! check whether node is a gasnode 155! 156 gasnode=.false. 157 call nident(itg,node,ntg,id) 158 if(id.gt.0) then 159 if(itg(id).eq.node) then 160 gasnode=.true. 161 endif 162 endif 163! 164! for the interpolation of submodels the undeformed 165! geometry is taken 166! 167 do j=1,3 168 coords(j)=co(j,node) 169 enddo 170! 171 entity='N' 172 one=1 173 iselect(1)=ndirboun(i)+1 174 call interpolsubmodel(integerglob,doubleglob,xbounact(i), 175 & coords,iselect,one,node,tieset,istartset,iendset, 176 & ialset,ntie,entity,set,nset) 177c write(*,*) 'tempload ',node,ndirboun(i),xbounact(i) 178 cycle 179 endif 180! 181 if(nam.gt.0) then 182 iambouni=iamboun(i) 183 else 184 iambouni=0 185 endif 186 if(iambouni.gt.0) then 187 xbounact(i)=xboun(i)*ampli(iambouni) 188 elseif(nmethod.eq.1) then 189 xbounact(i)=xbounold(i)+ 190 & (xboun(i)-xbounold(i))*reltime 191 else 192 xbounact(i)=xboun(i) 193 endif 194 enddo 195! 196! scaling the loading 197! 198 do i=1,nforc 199 if(ndirforc(i).eq.0) then 200 if((xforc(i).lt.1.2357111318d0).and. 201 & (xforc(i).gt.1.2357111316d0)) then 202! 203! user subroutine for the concentrated heat flux 204! 205 node=nodeforc(1,i) 206! 207! check whether node is a gasnode 208! 209 gasnode=.false. 210 call nident(itg,node,ntg,id) 211 if(id.gt.0) then 212 if(itg(id).eq.node) then 213 gasnode=.true. 214 endif 215 endif 216! 217 abqtime(1)=time 218 abqtime(2)=ttime+dtime 219! 220! a gasnode cannot move (displacement DOFs are used 221! for other purposes, e.g. mass flow and pressure) 222! 223 if(gasnode) then 224 do j=1,3 225 coords(j)=co(j,node) 226 enddo 227 else 228 do j=1,3 229 coords(j)=co(j,node)+vold(j,node) 230 enddo 231 endif 232! 233 call cflux(xforcact(i),msecpt,istep,iinc,abqtime,node, 234 & coords,vold,mi) 235 cycle 236 endif 237 else 238 if((xforc(i).lt.1.2357111318d0).and. 239 & (xforc(i).gt.1.2357111316d0)) then 240! 241! user subroutine for the concentrated load 242! 243 node=nodeforc(1,i) 244! 245 abqtime(1)=time 246 abqtime(2)=ttime+dtime 247! 248 do j=1,3 249 coords(j)=co(j,node)+vold(j,node) 250 enddo 251! 252 call cload(xforcact(i),istep,iinc,abqtime,node, 253 & ndirforc(i),coords,vold,mi,ntrans,trab,inotr,veold, 254 & nmethod,idummy,ddummy,ddummy) 255 cycle 256 endif 257 endif 258 if(nam.gt.0) then 259 iamforci=iamforc(i) 260 else 261 iamforci=0 262 endif 263 if(iamforci.gt.0) then 264 xforcact(i)=xforc(i)*ampli(iamforci) 265 elseif(nmethod.eq.1) then 266 xforcact(i)=xforcold(i)+ 267 & (xforc(i)-xforcold(i))*reltime 268 else 269 xforcact(i)=xforc(i) 270 endif 271 enddo 272! 273 do i=1,nload 274 ipresboun=0 275! 276! check for pressure boundary conditions 277! 278 if(sideload(i)(3:4).eq.'NP') then 279 node=nelemload(2,i) 280 idof=8*(node-1)+2 281 call nident(ikboun,idof,nboun,id) 282 if(id.gt.0) then 283 if(ikboun(id).eq.idof) then 284 ipresboun=1 285 xloadact(1,i)=xbounact(ilboun(id)) 286 endif 287 endif 288 endif 289! 290 if(ipresboun.eq.0) then 291 if(nam.gt.0) then 292 iamloadi1=iamload(1,i) 293 iamloadi2=iamload(2,i) 294 else 295 iamloadi1=0 296 iamloadi2=0 297 endif 298 if(iamloadi1.gt.0) then 299 xloadact(1,i)=xload(1,i)*ampli(iamloadi1) 300 elseif(nmethod.eq.1) then 301 xloadact(1,i)=xloadold(1,i)+ 302 & (xload(1,i)-xloadold(1,i))*reltime 303 else 304 xloadact(1,i)=xload(1,i) 305 endif 306 if(iamloadi2.gt.0) then 307 xloadact(2,i)=xload(2,i)*ampli(iamloadi2) 308c elseif(nmethod.eq.1) then 309c xloadact(2,i)=xload(2,i) 310 else 311 xloadact(2,i)=xload(2,i) 312 endif 313 endif 314 enddo 315! 316 do i=1,nbody 317 if(nam.gt.0) then 318 iambodyi=ibody(2,i) 319 else 320 iambodyi=0 321 endif 322 if(iambodyi.gt.0) then 323 xbodyact(1,i)=xbody(1,i)*ampli(iambodyi) 324 elseif(nmethod.eq.1) then 325 xbodyact(1,i)=xbodyold(1,i)+ 326 & (xbody(1,i)-xbodyold(1,i))*reltime 327 else 328 xbodyact(1,i)=xbody(1,i) 329 endif 330 enddo 331! 332! scaling the temperatures 333! 334 if(ithermal(1).eq.1) then 335 do i=1,nk 336 if((t1(i).lt.1.2357111318d0).and. 337 & (t1(i).gt.1.2357111316d0)) then 338! 339 abqtime(1)=time 340 abqtime(2)=ttime+dtime 341! 342 do j=1,3 343 coords(j)=co(j,i)+vold(j,i) 344 enddo 345 call utemp(t1act(i),msecpt,istep,iinc,abqtime,i, 346 & coords,vold,mi) 347 cycle 348 endif 349! 350 if((t1(i).lt.1.9232931375d0).and. 351 & (t1(i).gt.1.9232931373d0)) then 352! 353! for the interpolation of submodels the undeformed 354! geometry is taken 355! 356 do j=1,3 357 coords(j)=co(j,i) 358 enddo 359! 360 entity='N' 361 one=1 362 iselect(1)=1 363 call interpolsubmodel(integerglob,doubleglob,t1act(i), 364 & coords,iselect,one,i,tieset,istartset,iendset, 365 & ialset,ntie,entity,set,nset) 366 cycle 367 endif 368! 369 if(nam.gt.0) then 370 iamt1i=iamt1(i) 371 else 372 iamt1i=0 373 endif 374 if(iamt1i.gt.0) then 375 t1act(i)=t1(i)*ampli(iamt1i) 376 elseif(nmethod.eq.1) then 377 t1act(i)=t1old(i)+(t1(i)-t1old(i))*reltime 378 else 379 t1act(i)=t1(i) 380 endif 381 enddo 382! 383! taking temperature MPC's into account 384! 385 do j=1,nmpc 386 k=mod(ikmpc(j),8) 387 if(k.ne.0) cycle 388 i=ilmpc(j) 389 ist=ipompc(i) 390 node=nodempc(1,ist) 391 index=nodempc(3,ist) 392 fixed_temp=0.d0 393 if(index.ne.0) then 394 do 395 fixed_temp=fixed_temp- 396 & coefmpc(index)*t1act(nodempc(1,index)) 397 index=nodempc(3,index) 398 if(index.eq.0) exit 399 enddo 400 endif 401 t1act(node)=fixed_temp/coefmpc(ist) 402 enddo 403 endif 404c write(*,*) 'nboun' 405c do i=1,nboun 406c write(*,'(i7,1x,e11.4,1x,e11.4)') i,xbounact(i),xboun(i) 407c enddo 408c write(*,*) 'nforc' 409c do i=1,nforc 410c write(*,'(i7,1x,e11.4,1x,e11.4)') i,xforcact(i),xforc(i) 411c enddo 412c write(*,*) 'nload' 413c do i=1,nload 414c write(*,'(i7,1x,e11.4,1x,e11.4)') i,xloadact(1,i),xload(1,i) 415c enddo 416! 417 return 418 end 419