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 temploaddiff(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,
23     &     nmethod,xbounold,xboun,xbounact,iamboun,nboun,
24     &     nodeboun,ndirboun,nodeforc,ndirforc,istep,iinc,
25     &     co,vold,itg,ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
26     &     xforcdiff,xloaddiff,xbodydiff,t1diff,xboundiff,iabsload,
27     &     iprescribedboundary,ntrans,trab,inotr,veold,nactdof,bcont,fn,
28     &     ipobody,iponoel,inoel,ipkon,kon,lakon,ielprop,prop,ielmat,
29     &     shcon,nshcon,rhcon,nrhcon,ntmat_,cocon,ncocon)
30!
31!     calculates the loading at a given time and the difference with
32!     the last call of temploaddiff: is needed in the modal dynamic
33!     procedure (dyna.c, dynacont.c; speeds up the calculation)
34!
35      implicit none
36!
37      logical gasnode
38!
39      character*8 lakon(*)
40      character*20 sideload(*)
41      character*80 amname(*)
42!
43      integer iamforc(*),iamload(2,*),iamt1(*),nelemload(2,*),
44     &     nam,i,istart,iend,id,nforc,nload,nk,namta(3,*),ithermal(*),
45     &     nmethod,iamt1i,iamboun(*),nboun,iamforci,iambouni,
46     &     iamloadi1,iamloadi2,ibody(3,*),itg(*),ntg,idof,inotr(2,*),
47     &     nbody,iambodyi,nodeboun(*),ndirboun(*),nodeforc(2,*),
48     &     ndirforc(*),istep,iinc,msecpt,node,j,ikboun(*),ilboun(*),
49     &     ipresboun,mi(*),iabsload,iprescribedboundary,ntrans,
50     &     nactdof(0:mi(2),*),ipobody(2,*),iponoel(*),inoel(2,*),
51     &     ipkon(*),kon(*),ielprop(*),ielmat(mi(3),*),nshcon(*),
52     &     nrhcon(*),ncocon(2,*),ntmat_
53!
54      real*8 xforc(*),xforcact(*),xload(2,*),xloadact(2,*),
55     &     t1(*),t1act(*),amta(2,*),ampli(*),time,xforcdiff(*),
56     &     xforcold(*),xloadold(2,*),t1old(*),reltime,xloaddiff(2,*),
57     &     xbounold(*),xboun(*),xbounact(*),ttime,dtime,reftime,
58     &     xbody(7,*),xbodyold(7,*),xbodydiff(7,*),t1diff(*),
59     &     xbodyact(7,*),co(3,*),vold(0:mi(2),*),abqtime(2),coords(3),
60     &     xboundiff(*),trab(7,*), veold(0:mi(2),*),bcont(*),
61     &     prop(*),shcon(0:3,ntmat_,*),rhcon(0:1,ntmat_,*),
62     &     cocon(0:6,ntmat_,*),fn(0:mi(2),*)
63!
64!
65!
66      data msecpt /1/
67!
68!     if an amplitude is active, the loading is scaled according to
69!     the actual time. If no amplitude is active, then the load is
70!     - scaled according to the relative time for a static step
71!     - applied as a step loading for a dynamic step
72!
73!     calculating all amplitude values for the current time
74!
75      do i=1,nam
76        if(namta(3,i).lt.0) then
77          reftime=ttime+time
78        else
79          reftime=time
80        endif
81        if(abs(namta(3,i)).ne.i) then
82          reftime=reftime-amta(1,namta(1,i))
83          istart=namta(1,abs(namta(3,i)))
84          iend=namta(2,abs(namta(3,i)))
85          if(istart.eq.0) then
86            call uamplitude(reftime,amname(namta(3,i)),ampli(i))
87            cycle
88          endif
89        else
90          istart=namta(1,i)
91          iend=namta(2,i)
92          if(istart.eq.0) then
93            call uamplitude(reftime,amname(i),ampli(i))
94            cycle
95          endif
96        endif
97        call identamta(amta,reftime,istart,iend,id)
98        if(id.lt.istart) then
99          ampli(i)=amta(2,istart)
100        elseif(id.eq.iend) then
101          ampli(i)=amta(2,iend)
102        else
103          ampli(i)=amta(2,id)+(amta(2,id+1)-amta(2,id))
104     &         *(reftime-amta(1,id))/(amta(1,id+1)-amta(1,id))
105        endif
106      enddo
107!
108!     scaling the boundary conditions
109!
110      if(iprescribedboundary.eq.1) then
111        do i=1,nboun
112          if((xboun(i).lt.1.2357111318d0).and.
113     &         (xboun(i).gt.1.2357111316d0)) then
114!
115!     user subroutine for boundary conditions
116!
117            node=nodeboun(i)
118!
119!     check whether node is a gasnode
120!
121            gasnode=.false.
122            call nident(itg,node,ntg,id)
123            if(id.gt.0) then
124              if(itg(id).eq.node) then
125                gasnode=.true.
126              endif
127            endif
128!
129            abqtime(1)=time
130            abqtime(2)=ttime+time
131!
132!     a gasnode cannot move (displacement DOFs are used
133!     for other purposes, e.g. mass flow and pressure)
134!
135            if(gasnode) then
136              do j=1,3
137                coords(j)=co(j,node)
138              enddo
139            else
140              do j=1,3
141                coords(j)=co(j,node)+vold(j,node)
142              enddo
143            endif
144!
145            if(iabsload.eq.0) then
146              xboundiff(i)=xbounact(i)
147            else
148              xboundiff(i)=xbounact(i)-xboundiff(i)
149            endif
150            if(ndirboun(i).eq.0) then
151              call utemp(xbounact(i),msecpt,istep,iinc,abqtime,node,
152     &             coords,vold,mi,iponoel,inoel,
153     &             ipobody,xbodyact,ibody,ipkon,kon,
154     &             lakon,ielprop,prop,ielmat,
155     &             shcon,nshcon,rhcon,nrhcon,ntmat_,cocon,ncocon)
156            else
157              call uboun(xbounact(i),istep,iinc,abqtime,node,
158     &             ndirboun(i),coords,vold,mi,iponoel,inoel,
159     &             ipobody,xbodyact,ibody,ipkon,kon,
160     &             lakon,ielprop,prop,ielmat,
161     &             shcon,nshcon,rhcon,nrhcon,ntmat_,cocon,ncocon)
162            endif
163            xboundiff(i)=xbounact(i)-xboundiff(i)
164            cycle
165          endif
166!
167          if(nam.gt.0) then
168            iambouni=iamboun(i)
169          else
170            iambouni=0
171          endif
172!
173          if(iabsload.eq.0) then
174            xboundiff(i)=xbounact(i)
175          else
176            xboundiff(i)=xbounact(i)-xboundiff(i)
177          endif
178          if(iambouni.gt.0) then
179            if(amname(iambouni).eq.'RAMP12357111317') then
180              xbounact(i)=xbounold(i)+
181     &             (xboun(i)-xbounold(i))*reltime
182            else
183              xbounact(i)=xboun(i)*ampli(iambouni)
184            endif
185          elseif(nmethod.eq.1) then
186            xbounact(i)=xbounold(i)+
187     &           (xboun(i)-xbounold(i))*reltime
188          else
189            xbounact(i)=xboun(i)
190          endif
191          xboundiff(i)=xbounact(i)-xboundiff(i)
192        enddo
193      endif
194!
195!     scaling the loading
196!
197      do i=1,nforc
198        if(ndirforc(i).eq.0) then
199          if((xforc(i).lt.1.2357111318d0).and.
200     &         (xforc(i).gt.1.2357111316d0)) then
201            iabsload=2
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+time
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            if(iabsload.eq.0) then
234              xforcdiff(i)=xforcact(i)
235            else
236              xforcdiff(i)=xforcact(i)-xforcdiff(i)
237            endif
238            call cflux(xforcact(i),msecpt,istep,iinc,abqtime,node,
239     &           coords,vold,mi,ipkon,kon,lakon,iponoel,inoel,
240     &           ielprop,prop,ielmat,shcon,nshcon,rhcon,nrhcon,
241     &           ntmat_,cocon,ncocon)
242            xforcdiff(i)=xforcact(i)-xforcdiff(i)
243            cycle
244          endif
245        else
246          if((xforc(i).lt.1.2357111318d0).and.
247     &         (xforc(i).gt.1.2357111316d0)) then
248            iabsload=2
249!
250!     user subroutine for the concentrated force
251!
252            node=nodeforc(1,i)
253!
254            abqtime(1)=time
255            abqtime(2)=ttime+time
256!
257            do j=1,3
258              coords(j)=co(j,node)+vold(j,node)
259            enddo
260!
261            if(iabsload.eq.0) then
262              xforcdiff(i)=xforcact(i)
263            else
264              xforcdiff(i)=xforcact(i)-xforcdiff(i)
265            endif
266            call cload(xforcact(i),istep,iinc,abqtime,node,
267     &           ndirforc(i),coords,vold,mi,ntrans,trab,inotr,veold)
268            xforcdiff(i)=xforcact(i)-xforcdiff(i)
269            cycle
270          endif
271        endif
272        if(nam.gt.0) then
273          iamforci=iamforc(i)
274        else
275          iamforci=0
276        endif
277!
278        if(iabsload.eq.0) then
279          xforcdiff(i)=xforcact(i)
280        else
281          xforcdiff(i)=xforcact(i)-xforcdiff(i)
282        endif
283        if(iamforci.gt.0) then
284          if(amname(iamforci).eq.'RAMP12357111317') then
285            xforcact(i)=xforcold(i)+
286     &           (xforc(i)-xforcold(i))*reltime
287          else
288            xforcact(i)=xforc(i)*ampli(iamforci)
289          endif
290        elseif(nmethod.eq.1) then
291          xforcact(i)=xforcold(i)+
292     &         (xforc(i)-xforcold(i))*reltime
293        else
294          xforcact(i)=xforc(i)
295        endif
296        xforcdiff(i)=xforcact(i)-xforcdiff(i)
297      enddo
298!
299      do i=1,nload
300!
301!     check for dload subroutine
302!
303        if(sideload(i)(3:4).eq.'NU') then
304          iabsload=2
305          cycle
306        endif
307!
308        ipresboun=0
309!
310!     check for pressure boundary conditions
311!
312        if(sideload(i)(3:4).eq.'NP') then
313          node=nelemload(2,i)
314          idof=8*(node-1)+2
315          call nident(ikboun,idof,nboun,id)
316          if(id.gt.0) then
317            if(ikboun(id).eq.idof) then
318              ipresboun=1
319              if(iabsload.eq.0) then
320                xloaddiff(1,i)=xloadact(1,i)
321              else
322                xloaddiff(1,i)=xloadact(1,i)-xloaddiff(1,i)
323              endif
324              xloadact(1,i)=xbounact(ilboun(id))
325              xloaddiff(1,i)=xloadact(1,i)-xloaddiff(1,i)
326            endif
327          endif
328        endif
329!
330        if(ipresboun.eq.0) then
331          if(nam.gt.0) then
332            iamloadi1=iamload(1,i)
333            iamloadi2=iamload(2,i)
334          else
335            iamloadi1=0
336            iamloadi2=0
337          endif
338!
339          if(iabsload.eq.0) then
340            xloaddiff(1,i)=xloadact(1,i)
341          else
342            xloaddiff(1,i)=xloadact(1,i)-xloaddiff(1,i)
343          endif
344          if(iamloadi1.gt.0) then
345            if(amname(iamloadi1).eq.'RAMP12357111317') then
346              xloadact(1,i)=xloadold(1,i)+
347     &             (xload(1,i)-xloadold(1,i))*reltime
348            else
349              xloadact(1,i)=xload(1,i)*ampli(iamloadi1)
350            endif
351          elseif(nmethod.eq.1) then
352            xloadact(1,i)=xloadold(1,i)+
353     &           (xload(1,i)-xloadold(1,i))*reltime
354          else
355            xloadact(1,i)=xload(1,i)
356          endif
357          xloaddiff(1,i)=xloadact(1,i)-xloaddiff(1,i)
358!
359          if(iabsload.eq.0) then
360            xloaddiff(2,i)=xloadact(1,i)
361          else
362            xloaddiff(2,i)=xloadact(2,i)-xloaddiff(2,i)
363          endif
364          if(iamloadi2.gt.0) then
365            if(amname(iamloadi2).eq.'RAMP12357111317') then
366              xloadact(2,i)=xloadold(2,i)+
367     &             (xload(2,i)-xloadold(2,i))*reltime
368            else
369              xloadact(2,i)=xload(2,i)*ampli(iamloadi2)
370            endif
371          elseif(nmethod.eq.1) then
372            xloadact(2,i)=xload(2,i)
373          else
374            xloadact(2,i)=xload(2,i)
375          endif
376          xloaddiff(2,i)=xloadact(2,i)-xloaddiff(2,i)
377        endif
378      enddo
379!
380      do i=1,nbody
381        if(nam.gt.0) then
382          iambodyi=ibody(2,i)
383        else
384          iambodyi=0
385        endif
386!
387        if(iabsload.eq.0) then
388          xbodydiff(1,i)=xbodyact(1,i)
389        else
390          xbodydiff(1,i)=xbodyact(1,i)-xbodydiff(1,i)
391        endif
392        if(iambodyi.gt.0) then
393          if(amname(iambodyi).eq.'RAMP12357111317') then
394            xbodyact(1,i)=xbodyold(1,i)+
395     &           (xbody(1,i)-xbodyold(1,i))*reltime
396          else
397            xbodyact(1,i)=xbody(1,i)*ampli(iambodyi)
398          endif
399        elseif(nmethod.eq.1) then
400          xbodyact(1,i)=xbodyold(1,i)+
401     &         (xbody(1,i)-xbodyold(1,i))*reltime
402        else
403          xbodyact(1,i)=xbody(1,i)
404        endif
405        xbodydiff(1,i)=xbodyact(1,i)-xbodydiff(1,i)
406      enddo
407!
408      return
409      end
410