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(co,nk,konl,lakonl,p1,p2,omx,bodyfx,nbody,
20     &  ff,nelem,nmethod,rhcon,ielmat,ntmat_,vold,iperturb,nelemload,
21     &  sideload,xload,nload,idist,ttime,time,istep,iinc,dtime,
22     &  xloadold,reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,
23     &  veold,matname,mi,ielprop,prop)
24!
25!     computation of the rhs for the element with
26!     the topology in konl: only for nonlinear calculations (i.e.
27!     the field ff contains no temperature and eigenstress
28!     contributions)
29!
30!     RESTRICTION: the only material parameter occurring in the
31!     loading is the density. In the present subroutine the density
32!     is assumed to be temperature INDEPENDENT.
33!
34      implicit none
35!
36      logical ivolumeforce
37!
38      character*8 lakonl
39      character*20 sideload(*)
40      character*80 matname(*),amat
41!
42      integer mi(*),nk,konl(20),ielmat(mi(3),*),nbody,ifaceq(8,6),
43     &  nelemload(2,*),ielprop(*),null,
44     &  nelem,nmethod,iperturb(*),nload,idist,i,i1,j1,jj,jj1,id,kk,
45     &  ipointer,nope,nopes,j,k,ntmat_,i2,imat,ii,ig,mint2d,mint3d,
46     &  ifacet(6,4),ifacew(8,5),istep,iinc,layer,kspt,jltyp,iflag,
47     &  ipompc(*),nodempc(3,*),nmpc,ikmpc(*),ilmpc(*),iscale
48!
49      real*8 co(3,*),p1(3,2),p2(3,2),omx(2),bodyfx(3),veold(0:mi(2),*),
50     &  rhcon(0:1,ntmat_,*),xs2(3,7),xloadold(2,*),reltime,
51     &  bodyf(3),om(2),rho,bf(3),q(3),shpj(4,20),xl(3,20),
52     &  shp(4,20),voldl(3,20),xl2(3,8),xsj2(3),shp2(7,8),
53     &  vold(0:mi(2),*),prop(*),
54     &  xload(2,*),xi,et,ze,const,xsj,ff(60),weight,ttime,time,tvar(2),
55     &  coords(3),dtime,coefmpc(*)
56!
57      include "gauss.f"
58!
59      data ifaceq /4,3,2,1,11,10,9,12,
60     &            5,6,7,8,13,14,15,16,
61     &            1,2,6,5,9,18,13,17,
62     &            2,3,7,6,10,19,14,18,
63     &            3,4,8,7,11,20,15,19,
64     &            4,1,5,8,12,17,16,20/
65      data ifacet /1,3,2,7,6,5,
66     &             1,2,4,5,9,8,
67     &             2,3,4,6,10,9,
68     &             1,4,3,8,10,7/
69      data ifacew /1,3,2,9,8,7,0,0,
70     &             4,5,6,10,11,12,0,0,
71     &             1,2,5,4,7,14,10,13,
72     &             2,3,6,5,8,15,11,14,
73     &             4,6,3,1,12,15,9,13/
74      data iflag /3/
75      data null /0/
76!
77      tvar(1)=time
78      tvar(2)=ttime+time
79!
80      if(lakonl(4:4).eq.'2') then
81         nope=20
82         nopes=8
83      elseif(lakonl(4:4).eq.'8') then
84         nope=8
85         nopes=4
86      elseif(lakonl(4:5).eq.'10') then
87         nope=10
88         nopes=6
89      elseif(lakonl(4:4).eq.'4') then
90         nope=4
91         nopes=3
92      elseif(lakonl(4:5).eq.'15') then
93         nope=15
94      else
95         nope=6
96      endif
97!
98      if(lakonl(4:5).eq.'8R') then
99         mint2d=1
100         mint3d=1
101      elseif(lakonl(4:7).eq.'20RB') then
102         if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
103            mint2d=4
104            mint3d=50
105         else
106            mint2d=4
107            call beamintscheme(lakonl,mint3d,ielprop(nelem),prop,
108     &           null,xi,et,ze,weight)
109         endif
110      elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) then
111         mint2d=4
112         mint3d=8
113      elseif(lakonl(4:4).eq.'2') then
114         mint2d=9
115         mint3d=27
116      elseif(lakonl(4:5).eq.'10') then
117         mint2d=3
118         mint3d=4
119      elseif(lakonl(4:4).eq.'4') then
120         mint2d=1
121         mint3d=1
122      elseif(lakonl(4:5).eq.'15') then
123         mint3d=9
124      else
125         mint3d=2
126      endif
127!
128!     computation of the coordinates of the local nodes
129!
130      do i=1,nope
131        do j=1,3
132          xl(j,i)=co(j,konl(i))
133        enddo
134      enddo
135!
136!       initialisation for distributed forces
137!
138c      if(idist.ne.0) then
139         do i=1,3*nope
140            ff(i)=0.d0
141         enddo
142c      endif
143!
144!     displacements for 2nd order static and modal theory
145!
146c      if(iperturb.ne.0) then
147      if((iperturb(1).eq.1).or.(iperturb(2).eq.1)) then
148         do i1=1,nope
149            do i2=1,3
150               voldl(i2,i1)=vold(i2,konl(i1))
151            enddo
152         enddo
153      endif
154!
155!     calculating the density: no temperature dependence assumed!
156!     otherwise cfr. e_c3d.f
157!
158      imat=ielmat(1,nelem)
159      amat=matname(imat)
160      rho=rhcon(1,1,imat)
161!
162!     computation of the body forces
163!
164      ivolumeforce=.false.
165      if(nbody.ne.0) then
166         ivolumeforce=.true.
167         do ii=1,2
168            om(ii)=omx(ii)*rho
169         enddo
170         do ii=1,3
171            bodyf(ii)=bodyfx(ii)*rho
172         enddo
173      elseif(nload.gt.0) then
174         call nident2(nelemload,nelem,nload,id)
175         do
176            if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
177            if((sideload(id)(1:2).ne.'BX').and.
178     &           (sideload(id)(1:2).ne.'BY').and.
179     &           (sideload(id)(1:2).ne.'BZ')) then
180               id=id-1
181               cycle
182            else
183               ivolumeforce=.true.
184               exit
185            endif
186         enddo
187      endif
188!
189!     computation of the rhs: loop over the Gauss points
190!
191      if(ivolumeforce) then
192         do kk=1,mint3d
193            if(lakonl(4:5).eq.'8R') then
194               xi=gauss3d1(1,kk)
195               et=gauss3d1(2,kk)
196               ze=gauss3d1(3,kk)
197               weight=weight3d1(kk)
198            elseif(lakonl(4:7).eq.'20RB') then
199               if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
200                  xi=gauss3d13(1,kk)
201                  et=gauss3d13(2,kk)
202                  ze=gauss3d13(3,kk)
203                  weight=weight3d13(kk)
204               else
205                  call beamintscheme(lakonl,mint3d,ielprop(nelem),prop,
206     &                 kk,xi,et,ze,weight)
207               endif
208            elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R'))
209     &              then
210               xi=gauss3d2(1,kk)
211               et=gauss3d2(2,kk)
212               ze=gauss3d2(3,kk)
213               weight=weight3d2(kk)
214            elseif(lakonl(4:4).eq.'2') then
215               xi=gauss3d3(1,kk)
216               et=gauss3d3(2,kk)
217               ze=gauss3d3(3,kk)
218               weight=weight3d3(kk)
219            elseif(lakonl(4:5).eq.'10') then
220               xi=gauss3d5(1,kk)
221               et=gauss3d5(2,kk)
222               ze=gauss3d5(3,kk)
223               weight=weight3d5(kk)
224            elseif(lakonl(4:4).eq.'4') then
225               xi=gauss3d4(1,kk)
226               et=gauss3d4(2,kk)
227               ze=gauss3d4(3,kk)
228               weight=weight3d4(kk)
229            elseif(lakonl(4:5).eq.'15') then
230               xi=gauss3d8(1,kk)
231               et=gauss3d8(2,kk)
232               ze=gauss3d8(3,kk)
233               weight=weight3d8(kk)
234            else
235               xi=gauss3d7(1,kk)
236               et=gauss3d7(2,kk)
237               ze=gauss3d7(3,kk)
238               weight=weight3d7(kk)
239            endif
240!
241!     calculation of the shape functions and their derivatives
242!     in the gauss point
243!
244            if(nope.eq.20) then
245               if(lakonl(7:7).eq.'A') then
246                  call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag)
247               elseif((lakonl(7:7).eq.'E').or.(lakonl(7:7).eq.'S')) then
248                  call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag)
249               else
250                  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
251               endif
252            elseif(nope.eq.8) then
253               call shape8h(xi,et,ze,xl,xsj,shp,iflag)
254            elseif(nope.eq.10) then
255               call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
256            elseif(nope.eq.4) then
257               call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
258            elseif(nope.eq.15) then
259               call shape15w(xi,et,ze,xl,xsj,shp,iflag)
260            else
261               call shape6w(xi,et,ze,xl,xsj,shp,iflag)
262            endif
263!
264!     check the jacobian determinant
265!
266            if(xsj.lt.1.d-20) then
267               write(*,*) '*ERROR in e_c3d_rhs: nonpositive jacobian'
268               write(*,*) '         determinant in element',nelem
269               write(*,*)
270               xsj=dabs(xsj)
271               nmethod=0
272            endif
273!
274!     computation of the right hand side
275!
276!     incorporating the jacobian determinant in the shape
277!     functions
278!
279            if((nload.gt.0).or.(nbody.ne.0)) then
280               do i1=1,nope
281                  shpj(4,i1)=shp(4,i1)*xsj
282               enddo
283            endif
284!
285!     distributed body flux
286!
287            if(nload.gt.0) then
288               call nident2(nelemload,nelem,nload,id)
289               do
290                  if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
291                  if((sideload(id)(1:2).ne.'BX').and.
292     &                 (sideload(id)(1:2).ne.'BY').and.
293     &                 (sideload(id)(1:2).ne.'BZ')) then
294                     id=id-1
295                     cycle
296                  endif
297                  if(sideload(id)(3:4).eq.'NU') then
298                     do j=1,3
299                        coords(j)=0.d0
300                        do i1=1,nope
301                           coords(j)=coords(j)+
302     &                          shp(4,i1)*xl(j,i1)
303                        enddo
304                     enddo
305                     if(sideload(id)(1:2).eq.'BX') then
306                        jltyp=1
307                     elseif(sideload(id)(1:2).eq.'BY') then
308                        jltyp=2
309                     elseif(sideload(id)(1:2).eq.'BZ') then
310                        jltyp=3
311                     endif
312                     iscale=1
313                     call dload(xload(1,id),istep,iinc,tvar,nelem,i,
314     &                    layer,kspt,coords,jltyp,sideload(id),vold,co,
315     &                    lakonl,konl,ipompc,nodempc,coefmpc,nmpc,ikmpc,
316     &                    ilmpc,iscale,veold,rho,amat,mi)
317                     if((nmethod.eq.1).and.(iscale.ne.0))
318     &                    xload(1,id)=xloadold(1,id)+
319     &                    (xload(1,id)-xloadold(1,id))*reltime
320                  endif
321                  jj1=1
322                  do jj=1,nope
323                     if(sideload(id)(1:2).eq.'BX')
324     &                    ff(jj1)=ff(jj1)+xload(1,id)*shpj(4,jj)*weight
325                     if(sideload(id)(1:2).eq.'BY')
326     &                    ff(jj1+1)=ff(jj1+1)+xload(1,id)*shpj(4,jj)
327     &                              *weight
328                     if(sideload(id)(1:2).eq.'BZ')
329     &                    ff(jj1+2)=ff(jj1+2)+xload(1,id)*shpj(4,jj)
330     &                              *weight
331                     jj1=jj1+3
332                  enddo
333                  id=id-1
334               enddo
335            endif
336!
337!     body forces
338!
339            if(nbody.ne.0) then
340!
341               do i1=1,3
342                  bf(i1)=0.d0
343               enddo
344               do jj1=1,2
345                  if(dabs(om(jj1)).lt.1.d-20) cycle
346                  do i1=1,3
347!
348!     computation of the global coordinates of the gauss
349!     point
350!
351                     q(i1)=0.d0
352c                     if(iperturb.eq.0) then
353                     if((iperturb(1).ne.1).and.(iperturb(2).ne.1))
354     &                    then
355                        do j1=1,nope
356                           q(i1)=q(i1)+shp(4,j1)*xl(i1,j1)
357                        enddo
358                     else
359                        do j1=1,nope
360                           q(i1)=q(i1)+shp(4,j1)*
361     &                          (xl(i1,j1)+voldl(i1,j1))
362                        enddo
363                     endif
364!
365                     q(i1)=q(i1)-p1(i1,jj1)
366                  enddo
367                  const=q(1)*p2(1,jj1)+q(2)*p2(2,jj1)+q(3)*p2(3,jj1)
368!
369!     inclusion of the centrifugal force into the body force
370!
371                  do i1=1,3
372                     bf(i1)=bf(i1)+(q(i1)-const*p2(i1,jj1))*om(jj1)
373                  enddo
374               enddo
375!
376               do i1=1,3
377                  bf(i1)=bf(i1)+bodyf(i1)
378               enddo
379!
380               jj1=1
381               do jj=1,nope
382                  ff(jj1)=ff(jj1)+bf(1)*shpj(4,jj)*weight
383                  ff(jj1+1)=ff(jj1+1)+bf(2)*shpj(4,jj)*weight
384                  ff(jj1+2)=ff(jj1+2)+bf(3)*shpj(4,jj)*weight
385                  jj1=jj1+3
386               enddo
387            endif
388!
389         enddo
390      endif
391!
392!     distributed loads
393!
394      if(nload.eq.0) then
395         return
396      endif
397      call nident2(nelemload,nelem,nload,id)
398      do
399         if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
400         if(sideload(id)(1:1).ne.'P') then
401            id=id-1
402            cycle
403         endif
404         read(sideload(id)(2:2),'(i1)') ig
405!
406!     treatment of wedge faces
407!
408         if(lakonl(4:4).eq.'6') then
409            mint2d=1
410            if(ig.le.2) then
411               nopes=3
412            else
413               nopes=4
414            endif
415         endif
416         if(lakonl(4:5).eq.'15') then
417            if(ig.le.2) then
418               mint2d=3
419               nopes=6
420            else
421               mint2d=4
422               nopes=8
423            endif
424         endif
425!
426         if((nope.eq.20).or.(nope.eq.8)) then
427c            if(iperturb.eq.0) then
428            if((iperturb(1).ne.1).and.(iperturb(2).ne.1)) then
429               do i=1,nopes
430                  do j=1,3
431                     xl2(j,i)=co(j,konl(ifaceq(i,ig)))
432                  enddo
433               enddo
434            else
435               do i=1,nopes
436                  do j=1,3
437                     xl2(j,i)=co(j,konl(ifaceq(i,ig)))+
438     &                    vold(j,konl(ifaceq(i,ig)))
439                  enddo
440               enddo
441            endif
442         elseif((nope.eq.10).or.(nope.eq.4)) then
443c            if(iperturb.eq.0) then
444            if((iperturb(1).ne.1).and.(iperturb(2).ne.1)) then
445               do i=1,nopes
446                  do j=1,3
447                     xl2(j,i)=co(j,konl(ifacet(i,ig)))
448                  enddo
449               enddo
450            else
451               do i=1,nopes
452                  do j=1,3
453                     xl2(j,i)=co(j,konl(ifacet(i,ig)))+
454     &                    vold(j,konl(ifacet(i,ig)))
455                  enddo
456               enddo
457            endif
458         else
459c            if(iperturb.eq.0) then
460            if((iperturb(1).ne.1).and.(iperturb(2).ne.1)) then
461               do i=1,nopes
462                  do j=1,3
463                     xl2(j,i)=co(j,konl(ifacew(i,ig)))
464                  enddo
465               enddo
466            else
467               do i=1,nopes
468                  do j=1,3
469                     xl2(j,i)=co(j,konl(ifacew(i,ig)))+
470     &                    vold(j,konl(ifacew(i,ig)))
471                  enddo
472               enddo
473            endif
474         endif
475!
476         do i=1,mint2d
477            if((lakonl(4:5).eq.'8R').or.
478     &           ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
479               xi=gauss2d1(1,i)
480               et=gauss2d1(2,i)
481               weight=weight2d1(i)
482            elseif((lakonl(4:4).eq.'8').or.
483     &              (lakonl(4:6).eq.'20R').or.
484     &              ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
485               xi=gauss2d2(1,i)
486               et=gauss2d2(2,i)
487               weight=weight2d2(i)
488            elseif(lakonl(4:4).eq.'2') then
489               xi=gauss2d3(1,i)
490               et=gauss2d3(2,i)
491               weight=weight2d3(i)
492            elseif((lakonl(4:5).eq.'10').or.
493     &              ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
494               xi=gauss2d5(1,i)
495               et=gauss2d5(2,i)
496               weight=weight2d5(i)
497            elseif((lakonl(4:4).eq.'4').or.
498     &              ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
499               xi=gauss2d4(1,i)
500               et=gauss2d4(2,i)
501               weight=weight2d4(i)
502            endif
503!
504            if(nopes.eq.8) then
505               call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
506            elseif(nopes.eq.4) then
507               call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
508            elseif(nopes.eq.6) then
509               call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
510            else
511               call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
512            endif
513!
514!     for nonuniform load: determine the coordinates of the
515!     point (transferred into the user subroutine)
516!
517            if(sideload(id)(3:4).eq.'NU') then
518               do k=1,3
519                  coords(k)=0.d0
520                  do j=1,nopes
521                     coords(k)=coords(k)+xl2(k,j)*shp2(4,j)
522                  enddo
523               enddo
524               read(sideload(id)(2:2),'(i1)') jltyp
525               jltyp=jltyp+20
526               iscale=1
527               call dload(xload(1,id),istep,iinc,tvar,nelem,i,layer,
528     &              kspt,coords,jltyp,sideload(id),vold,co,lakonl,
529     &              konl,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,
530     &              iscale,veold,rho,amat,mi)
531               if((nmethod.eq.1).and.(iscale.ne.0))
532     &              xload(1,id)=xloadold(1,id)+
533     &              (xload(1,id)-xloadold(1,id))*reltime
534            endif
535!
536            do k=1,nopes
537               if((nope.eq.20).or.(nope.eq.8)) then
538                  ipointer=(ifaceq(k,ig)-1)*3
539               elseif((nope.eq.10).or.(nope.eq.4)) then
540                  ipointer=(ifacet(k,ig)-1)*3
541               else
542                  ipointer=(ifacew(k,ig)-1)*3
543               endif
544               ff(ipointer+1)=ff(ipointer+1)-shp2(4,k)*xload(1,id)
545     &              *xsj2(1)*weight
546               ff(ipointer+2)=ff(ipointer+2)-shp2(4,k)*xload(1,id)
547     &              *xsj2(2)*weight
548               ff(ipointer+3)=ff(ipointer+3)-shp2(4,k)*xload(1,id)
549     &              *xsj2(3)*weight
550            enddo
551         enddo
552!
553         id=id-1
554      enddo
555!
556      return
557      end
558
559