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