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_corio(co,nk,konl,lakonl,p1,p2,omx,bodyfx,nbody,s,sm,
20     &  ff,nelem,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
21     &  ielmat,ielorien,norien,orab,ntmat_,
22     &  t0,t1,ithermal,vold,iperturb,nelemload,
23     &  sideload,xload,nload,idist,sti,stx,iexpl,plicon,
24     &  nplicon,plkcon,nplkcon,xstiff,npmat_,dtime,
25     &  matname,mi,ncmat_,ttime,time,istep,iinc,nmethod,ielprop,prop)
26!
27!     computation of the element matrix and rhs for the element with
28!     the topology in konl
29!
30!     ff: rhs without temperature and eigenstress contribution
31!
32!     nmethod=0: check for positive Jacobian
33!     nmethod=1: stiffness matrix + right hand side
34!     nmethod=2: stiffness matrix + mass matrix
35!     nmethod=3: static stiffness + buckling stiffness
36!     nmethod=4: stiffness matrix + mass matrix
37!
38      implicit none
39!
40      logical mass,stiffness,buckling,rhsi
41!
42      character*8 lakonl
43      character*20 sideload(*)
44      character*80 matname(*),amat
45!
46      integer konl(20),ifaceq(8,6),nelemload(2,*),nk,nbody,nelem,
47     &  mi(*),mattyp,ithermal(*),iperturb(*),nload,idist,i,j,k,l,i1,i2,
48     &  j1,nmethod,k1,l1,ii,jj,ii1,jj1,id,ipointer,ig,m1,m2,m3,m4,kk,
49     &  nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),
50     &  ielorien(mi(3),*),null,ielprop(*),
51     &  ntmat_,nope,nopes,norien,ihyper,iexpl,kode,imat,mint2d,
52     &  mint3d,ifacet(7,4),nopev,iorien,istiff,ncmat_,
53     &  ifacew(8,5),intscheme,n,ipointeri,ipointerj,istep,iinc,
54     &  layer,kspt,jltyp,iflag,iperm(60),m,iscale,ne0,
55     &  nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_
56!
57      real*8 co(3,*),xl(3,20),shp(4,20),xs2(3,7),prop(*),
58     &  s(60,60),w(3,3),p1(3),p2(3),bodyf(3),bodyfx(3),ff(60),
59     &  bf(3),q(3),shpj(4,20),elcon(0:ncmat_,ntmat_,*),
60     &  rhcon(0:1,ntmat_,*),xkl(3,3),eknlsign,reltime,
61     &  alcon(0:6,ntmat_,*),alzero(*),orab(7,*),t0(*),t1(*),
62     &  anisox(3,3,3,3),voldl(0:mi(2),20),vo(3,3),
63     &  xl2(3,8),xsj2(3),shp2(7,8),vold(0:mi(2),*),xload(2,*),
64     &  v(3,3,3,3),
65     &  om,omx,e,un,al,um,xi,et,ze,tt,const,xsj,xsjj,sm(60,60),
66     &  sti(6,mi(1),*),stx(6,mi(1),*),s11,s22,s33,s12,s13,s23,s11b,
67     &  s22b,s33b,s12b,s13b,s23b,t0l,t1l,
68     &  senergy,senergyb,rho,elas(21),
69     &  sume,factorm,factore,alp,elconloc(21),eth(6),
70     &  weight,coords(3),dmass,xl1(3,8),term
71!
72      real*8 plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
73     &  xstiff(27,mi(1),*),plconloc(802),dtime,ttime,time,
74     &  sax(60,60),ffax(60),gs(8,4),a
75!
76      data iflag /3/
77      data iperm /13,14,-15,16,17,-18,19,20,-21,22,23,-24,
78     &            1,2,-3,4,5,-6,7,8,-9,10,11,-12,
79     &            37,38,-39,40,41,-42,43,44,-45,46,47,-48,
80     &            25,26,-27,28,29,-30,31,32,-33,34,35,-36,
81     &            49,50,-51,52,53,-54,55,56,-57,58,59,-60/
82!
83      include "gauss.f"
84!
85      null=0
86!
87      imat=ielmat(1,nelem)
88      amat=matname(imat)
89!
90c     Bernhardi start
91      if(lakonl(1:5).eq.'C3D8I') then
92         nope=11
93      elseif(lakonl(4:4).eq.'2') then
94c     Bernhardi end
95         nope=20
96      elseif(lakonl(4:4).eq.'8') then
97         nope=8
98      elseif(lakonl(4:5).eq.'10') then
99         nope=10
100      elseif(lakonl(4:4).eq.'4') then
101         nope=4
102      elseif(lakonl(4:5).eq.'15') then
103         nope=15
104      elseif(lakonl(4:4).eq.'6') then
105         nope=6
106      elseif(lakonl(1:2).eq.'ES') then
107         read(lakonl(8:8),'(i1)') nope
108         nope=nope+1
109      endif
110!
111      if(lakonl(4:5).eq.'8R') then
112         mint3d=1
113      elseif(lakonl(4:7).eq.'20RB') then
114         if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
115            mint3d=50
116         else
117            call beamintscheme(lakonl,mint3d,ielprop(nelem),prop,
118     &           null,xi,et,ze,weight)
119         endif
120      elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) then
121         if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'S').or.
122     &        (lakonl(7:7).eq.'E')) then
123            mint3d=4
124         else
125            mint3d=8
126         endif
127      elseif(lakonl(4:4).eq.'2') then
128         mint3d=27
129      elseif(lakonl(4:5).eq.'10') then
130         mint3d=4
131      elseif(lakonl(4:4).eq.'4') then
132         mint3d=1
133      elseif(lakonl(4:5).eq.'15') then
134         mint3d=9
135      elseif(lakonl(4:4).eq.'6') then
136         mint3d=2
137      else
138         mint3d=0
139      endif
140!
141!     computation of the coordinates of the local nodes
142!
143      do i=1,nope
144         do j=1,3
145            xl(j,i)=co(j,konl(i))
146         enddo
147      enddo
148!
149!     initialisation of s
150!
151      do i=1,3*nope
152        do j=1,3*nope
153          s(i,j)=0.d0
154        enddo
155      enddo
156!
157!     computation of the matrix: loop over the Gauss points
158!
159      do kk=1,mint3d
160         if(lakonl(4:5).eq.'8R') then
161            xi=gauss3d1(1,kk)
162            et=gauss3d1(2,kk)
163            ze=gauss3d1(3,kk)
164            weight=weight3d1(kk)
165         elseif(lakonl(4:7).eq.'20RB') then
166            if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
167               xi=gauss3d13(1,kk)
168               et=gauss3d13(2,kk)
169               ze=gauss3d13(3,kk)
170               weight=weight3d13(kk)
171            else
172               call beamintscheme(lakonl,mint3d,ielprop(nelem),prop,
173     &              kk,xi,et,ze,weight)
174            endif
175         elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R'))
176     &           then
177            xi=gauss3d2(1,kk)
178            et=gauss3d2(2,kk)
179            ze=gauss3d2(3,kk)
180            weight=weight3d2(kk)
181         elseif(lakonl(4:4).eq.'2') then
182            xi=gauss3d3(1,kk)
183            et=gauss3d3(2,kk)
184            ze=gauss3d3(3,kk)
185            weight=weight3d3(kk)
186         elseif(lakonl(4:5).eq.'10') then
187            xi=gauss3d5(1,kk)
188            et=gauss3d5(2,kk)
189            ze=gauss3d5(3,kk)
190            weight=weight3d5(kk)
191         elseif(lakonl(4:4).eq.'4') then
192            xi=gauss3d4(1,kk)
193            et=gauss3d4(2,kk)
194            ze=gauss3d4(3,kk)
195            weight=weight3d4(kk)
196         elseif(lakonl(4:5).eq.'15') then
197            xi=gauss3d8(1,kk)
198            et=gauss3d8(2,kk)
199            ze=gauss3d8(3,kk)
200            weight=weight3d8(kk)
201         elseif(lakonl(4:4).eq.'6') then
202            xi=gauss3d7(1,kk)
203            et=gauss3d7(2,kk)
204            ze=gauss3d7(3,kk)
205            weight=weight3d7(kk)
206         endif
207!
208!     calculation of the shape functions and their derivatives
209!     in the gauss point
210!
211c     Bernhardi start
212         if(lakonl(1:5).eq.'C3D8R') then
213            call shape8hr(xl,xsj,shp,gs,a)
214         elseif(lakonl(1:5).eq.'C3D8I') then
215            call shape8hu(xi,et,ze,xl,xsj,shp,iflag)
216         else if(nope.eq.20) then
217c     Bernhardi end
218            if(lakonl(7:7).eq.'A') then
219               call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag)
220            elseif((lakonl(7:7).eq.'E').or.(lakonl(7:7).eq.'S')) then
221               call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag)
222            else
223               call shape20h(xi,et,ze,xl,xsj,shp,iflag)
224            endif
225         elseif(nope.eq.8) then
226            call shape8h(xi,et,ze,xl,xsj,shp,iflag)
227         elseif(nope.eq.10) then
228            call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
229         elseif(nope.eq.4) then
230            call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
231         elseif(nope.eq.15) then
232            call shape15w(xi,et,ze,xl,xsj,shp,iflag)
233         else
234            call shape6w(xi,et,ze,xl,xsj,shp,iflag)
235         endif
236!
237!           check the jacobian determinant
238!
239         if(xsj.lt.1.d-20) then
240            write(*,*) '*ERROR in e_c3d: nonpositive jacobian'
241            write(*,*) '       determinant in element',nelem
242            write(*,*)
243            xsj=dabs(xsj)
244            nmethod=0
245         endif
246!
247!           material data and local stiffness matrix
248!
249         istiff=1
250         call materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
251     &        imat,amat,iorien,coords,orab,ntmat_,elas,rho,
252     &        nelem,ithermal,alzero,mattyp,t0l,t1l,
253     &        ihyper,istiff,elconloc,eth,kode,plicon,
254     &        nplicon,plkcon,nplkcon,npmat_,
255     &        plconloc,mi(1),dtime,kk,
256     &        xstiff,ncmat_)
257!
258!
259!           initialisation for the body forces
260!
261         om=2.d0*rho*dsqrt(omx)*weight
262!
263!           incorporating the jacobian determinant in the shape
264!           functions
265!
266         xsjj=dsqrt(xsj)
267         do i1=1,nope
268            shpj(1,i1)=shp(1,i1)*xsjj
269            shpj(2,i1)=shp(2,i1)*xsjj
270            shpj(3,i1)=shp(3,i1)*xsjj
271            shpj(4,i1)=shp(4,i1)*xsj
272         enddo
273!
274         jj1=1
275         do jj=1,nope
276!
277            ii1=1
278            do ii=1,jj
279!
280!     Coriolis matrix
281!
282               dmass=
283     &              om*shpj(4,ii)*shp(4,jj)
284               s(ii1,jj1+1)=s(ii1,jj1+1)-p2(3)*dmass
285               s(ii1,jj1+2)=s(ii1,jj1+2)+p2(2)*dmass
286               s(ii1+1,jj1)=s(ii1+1,jj1)+p2(3)*dmass
287               s(ii1+1,jj1+2)=s(ii1+1,jj1+2)-p2(1)*dmass
288               s(ii1+2,jj1)=s(ii1+2,jj1)-p2(2)*dmass
289               s(ii1+2,jj1+1)=s(ii1+2,jj1+1)+p2(1)*dmass
290!
291               ii1=ii1+3
292            enddo
293            jj1=jj1+3
294         enddo
295      enddo
296!
297!     for axially symmetric and plane stress/strain elements:
298!     complete s and sm
299!
300      if((lakonl(6:7).eq.'RA').or.(lakonl(6:7).eq.'RS').or.
301     &     (lakonl(6:7).eq.'RE')) then
302         do i=1,60
303            do j=i,60
304               k=abs(iperm(i))
305               l=abs(iperm(j))
306               if(k.gt.l) then
307                  m=k
308                  k=l
309                  l=m
310               endif
311               sax(i,j)=s(k,l)*iperm(i)*iperm(j)/(k*l)
312            enddo
313         enddo
314         do i=1,60
315            do j=i,60
316               s(i,j)=s(i,j)+sax(i,j)
317            enddo
318         enddo
319!
320      endif
321!
322      return
323      end
324
325