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