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