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! ------------------------------------------------------------- 20 subroutine e_c3d_us45(co,kon,lakonl,p1,p2,omx,bodyfx,nbody,s,sm, 21 & ff,nelem,nmethod,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero, 22 & ielmat,ielorien,norien,orab,ntmat_, 23 & t0,t1,ithermal,vold,iperturb,nelemload, 24 & sideload,xload,nload,idist,sti,stx,iexpl,plicon, 25 & nplicon,plkcon,nplkcon,xstiff,npmat_,dtime, 26 & matname,mi,ncmat_,mass,stiffness,buckling,rhsi,intscheme, 27 & ttime,time,istep,iinc,coriolis,xloadold,reltime, 28 & ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,veold, 29 & ne0,ipkon,thicke, 30 & integerglob,doubleglob,tieset,istartset,iendset,ialset,ntie, 31 & nasym,ielprop,prop) 32! 33! 34! 35! INPUT: 36! 37! co(1..3,i) coordinates of node i 38! kon(*) contains the topology of all elements. The 39! topology of element i starts at kon(ipkon(i)+1) 40! and continues until all nodes are covered. The 41! number of nodes depends on the element label 42! lakonl label of current element nelem (character*8) 43! p1(1..3) coordinates of a point on the rotation axis 44! (if applicable) 45! p2(1..3) unit vector on rotation axis (if applicable) 46! omx rotational speed square 47! bodyfx(1..3) acceleration vector 48! nbody number of body loads 49! nelem element number 50! nmethod procedure: 51! 1: static analysis 52! 2: frequency analysis 53! 3: buckling analysis 54! 4: (linear or nonlinear) dynamic analysis 55! 5: steady state dynamics analysis 56! 6: Coriolis frequency calculation 57! 7: flutter frequency calculation 58! 8: magnetostatics 59! 9: magnetodynamics 60! 10: electromagnetic eigenvalue problems 61! 11: superelement creation or Green function 62! calculation 63! 12: sensitivity analysis 64! elcon elastic constants (cf. List of variables and 65! their meaning in the User's Manual) 66! nelcon integers describing the elastic constant fields 67! (cf. User's Manual) 68! rhcon density constants (cf. User's Manual) 69! nrhcon integers describing the density constant fields 70! (cf. User's Manual) 71! alcon thermal expansion constants (cf. User's Manual) 72! nalcon integers describing the thermal expansion 73! constants (cf. User's Manual) 74! alzero thermal expansion reference values (cf. User's Manual) 75! ielmat(i) material for element i 76! ielorien(i) orientation for element i 77! norien number of orientations 78! orab(7,*) description of all local coordinate systems. 79! (cf. List of variables and their meaning in the 80! User's manual) 81! ntmat_ maximum number of material temperature points 82! t0(i) temperature in node i at start of calculation 83! t1(i) temperature in node i at the end of the current 84! increment 85! ithermal cf. ithermal(1) in the List of variables and 86! their meaning in the User's Manual 87! vold(0..mi(2),i) value of the variables in node i at the end 88! of the previous iteration 89! iperturb(*) describes the kind of nonlinearity of the 90! calculation, cf. User's Manual 91! nelemload field describing facial distributed load (cf. 92! User's Manual) 93! sideload field describing facial distributed load (cf. 94! User's Manual) 95! xload facial distributed load values (cf. 96! User's Manual) 97! nload number of facial distributed loads 98! idist 0: no distributed forces in the model 99! 1: distributed forces are present in the model 100! (body forces or thermal loads or residual 101! stresses or distributed facial loads) 102! sti(i,j,k) stress component i in integration point j 103! of element k at the end of the previous step 104! stx(i,j,k) stress component i in integration point j 105! of element k at the end of a static step 106! describing a reference buckling load 107! iexpl 0: structure: implicit dynamics, 108! fluid: incompressible 109! iexpl 1: structure: implicit dynamics, 110! fluid: compressible 111! iexpl 2: structure: explicit dynamics, 112! fluid: incompressible 113! iexpl 3: structure: explicit dynamics, 114! fluid: compressible 115! plicon,nplicon fields describing isotropic hardening of 116! a plastic material or spring constants of 117! a nonlinear spring (cf. User's Manual) 118! plkcon,nplkcon fields describing kinematic hardening of 119! a plastic material or gap conductance 120! constants (cf. User's Manual) 121! xstiff(i,j,k) stiffness (i=1...21) and conductivity 122! (i=22..27) constants in integration point j 123! of element k 124! npmat_ maximum number of plastic constants 125! dtime length of present time increment 126! matname(i) name of material i 127! mi(1) max # of integration points per element (max 128! over all elements) 129! mi(2) max degree of freedom per node (max over all 130! nodes) in fields like v(0:mi(2))... 131! mi(3) max number of layers in the structure 132! ncmat_ max number of elastic constants 133! mass(1) if 1: mass matrix needed, else not 134! mass(2) if 1: heat capacity matrix needed, else not 135! stiffness if 1: stiffness matrix needed, else not 136! buckling if 1: linear buckling calculation, else not 137! rhsi if 1: right hand side needed, else not 138! intscheme 1: integration point scheme for the calculation 139! of the initial acceleration in a dynamic 140! step 141! 0: any other case 142! ttime total time at the start of the present increment 143! time step time at the end of the present increment 144! istep current step number 145! iinc current increment number within the actual step 146! coriolis if 1: Coriolis forces are taken into account, 147! else not 148! xloadold load at the end of the previous increment 149! (cf. User's Manual) 150! reltime relative step time (between 0 and 1) 151! ipompc(i) pointer into field nodempc and coefmpc for 152! MPC i 153! nodempc field containing the nodes and dof's of all 154! MPC's (cf. User's Manual) 155! coefmpc field containing the coefficients of all 156! MPC's (cf. User's Manual) 157! nmpc number of MPC's in the model 158! ikmpc(i) field containing an integer code number 159! describing the dependent dof of some MPC: 160! ikmpc(i)=8*(node-1)+dir, where the dependent 161! dof is applied in direction dir of node "node" 162! ikmpc is sorted in ascending order 163! ilmpc(i) MPC number corresponding to ikmpc(i) (cf. 164! User's manual) 165! veold(j,i) time rate of variable j in node i at the end 166! of the previous iteration 167! ne0 element number before the first contact element; 168! contact elements are appended after ne0 169! ipkon(i) points to the location in field kon preceding 170! the topology of element i 171! thicke(j,i) thickness of layer j in node i 172! integerglob integer field needed for determining the 173! interface loading of a submodel 174! doubleglob real field needed for determining the 175! interface loading of a submodel 176! tieset(1..3,*) tie character information for tie i (cf. 177! User's Manual) 178! istartset integer describing set information (cf. 179! User's Manual) 180! iendset integer describing set information (cf. 181! User's Manual) 182! ialset integer describing set information (cf. 183! User's Manual) 184! ntie number of ties in the model 185! nasym 1: asymmetric contributions, else purely 186! symmetric 187! ielprop(i) points to the location in field prop preceding 188! the properties of element i 189! prop(*) contains the properties and some beam 190! elements (cf. User's Manual) 191! 192! 193! OUTPUT: 194! 195! s element stiffness matrix 196! sm element mass matrix 197! ff element load vector 198! xload facial distributed load values (cf. 199! User's Manual) 200! nmethod may be changed by the user, e.g. a 201! change to 0 for negative Jacobians leads 202! to a program exit in the calling program 203! for the other value: cf. above 204! 205 206! 207 integer mass,stiffness,buckling,rhsi,coriolis 208! 209 character*8 lakonl 210 character*20 sideload(*) 211 character*80 matname(*),amat 212 character*81 tieset(3,*) 213! 214 integer konl(26),nelemload(2,*),nbody,nelem,mi(*),kon(*), 215 & ielprop(*),index,mattyp,ithermal,iperturb(*),nload,idist, 216 & i,j,k,i1,nmethod,kk,nelcon(2,*),nrhcon(*),nalcon(2,*), 217 & ielmat(mi(3),*),ielorien(mi(3),*),ipkon(*),indexe, 218 & ntmat_,nope,norien,ihyper,iexpl,kode,imat,iorien,istiff, 219 & ncmat_,intscheme,istep,iinc,ipompc(*),nodempc(3,*), 220 & nmpc,ikmpc(*),ilmpc(*),ne0,ndof,istartset(*),iendset(*), 221 & ialset(*),ntie,integerglob(*),nasym,nplicon(0:ntmat_,*), 222 & nplkcon(0:ntmat_,*),npmat_,jjj 223! 224 real*8 co(3,*),xl(3,26),veold(0:mi(2),*),rho,s(60,60),bodyfx(3), 225 & ff(60),elconloc(21),coords(3),p1(3),elcon(0:ncmat_,ntmat_,*), 226 & p2(3),eth(6),rhcon(0:1,ntmat_,*),reltime,prop(*),tm(3,3), 227 & alcon(0:6,ntmat_,*),alzero(*),orab(7,*),t0(*),t1(*), 228 & xloadold(2,*),vold(0:mi(2),*),xload(2,*),omx,e,un,um,tt, 229 & sm(60,60),sti(6,mi(1),*),stx(6,mi(1),*),t0l,t1l,coefmpc(*), 230 & elas(21),thicke(mi(3),*),doubleglob(*),dl,e2(3),e3(3), 231 & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*), 232 & xstiff(27,mi(1),*),plconloc(802),dtime,ttime,time,tmg(24,24), 233 & a,xi11,xi12,xi22,xk,e1(3),offset1,offset2,y1,y2,y3,z1,z2,z3, 234 & sg(24,24), h, Dm(3,3),q1,Nrs(4), dNr(4), dNs(4),x(4,3), xc(3), 235 & Jm(2,2), xg(4,3), invJm(2,2), detinvJm, detJm, dNx(4), dNy(4), 236 & bm(3,24), g_p(4,3), Kmem(24,24),ri,si,Kb(24,24), bb(3,24), 237 & bs(2,24), Ks(24,24), Ds(2,2), Db(3,3), Gt,Kshell(24,24),kdmax, 238 & m_3t(6,6), N_u(6,24),Mshell(24,24),gpthick(3,2),dett,dettt, 239 & Qin(3,3),Qs(2,2) 240! 241 intent(in) co,kon,lakonl,p1,p2,omx,bodyfx,nbody, 242 & nelem,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero, 243 & ielmat,ielorien,norien,orab,ntmat_, 244 & t0,t1,ithermal,vold,iperturb,nelemload, 245 & sideload,nload,idist,sti,stx,iexpl,plicon, 246 & nplicon,plkcon,nplkcon,xstiff,npmat_,dtime, 247 & matname,mi,ncmat_,mass,stiffness,buckling,rhsi,intscheme, 248 & ttime,time,istep,iinc,coriolis,xloadold,reltime, 249 & ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,veold, 250 & ne0,ipkon,thicke, 251 & integerglob,doubleglob,tieset,istartset,iendset,ialset,ntie, 252 & nasym,ielprop,prop 253 254 intent(inout) s,sm,ff,xload,nmethod 255 indexe=ipkon(nelem) 256 ! 257 ! Gauß points(4): 258 g_p(1,1) = -0.577350269189626 259 g_p(2,1) = +0.577350269189626 260 g_p(3,1) = +0.577350269189626 261 g_p(4,1) = -0.577350269189626 262 ! 263 g_p(1,2) = -0.577350269189626 264 g_p(2,2) = -0.577350269189626 265 g_p(3,2) = +0.577350269189626 266 g_p(4,2) = +0.577350269189626 267 ! 268 g_p(1,3) = +1.000000000000000 269 g_p(2,3) = +1.000000000000000 270 g_p(3,3) = +1.000000000000000 271 g_p(4,3) = +1.000000000000000 272 ! 273 gpthick(1,1) = +1.d0 274 gpthick(2,1) = 0.d0 275 gpthick(3,1) = -1.d0 276 gpthick(1,2) = 2.d0/6.d0 277 gpthick(2,2) = 8.d0/6.d0 278 gpthick(3,2) = 2.d0/6.d0 279 280! 281! material and orientation 282! 283 imat=ielmat(1,nelem) 284 amat=matname(imat) 285 if(norien.gt.0) then 286 iorien=ielorien(1,nelem) 287 else 288 iorien=0 289 endif 290! 291 if(nelcon(1,imat).lt.0) then 292 ihyper=1 293 else 294 ihyper=0 295 endif 296! 297! properties of the cross section 298! 299 index=ielprop(nelem) 300 h = prop(index+1) ! general beam section shell section?! 301 dett = h/2.d0 302 dettt = h**3/8.d0 303 ! ********************************************************** 304 nope = 4 ! number of element nodes 305 ndof = 6 ! u,v,w,rx,ry,rz 306 ! ********************************************************** 307 ! 308 ! initialisation for distributed forces 309 ! 310 if(rhsi.eq.1) then 311 if(idist.ne.0) then 312 ff(:) = 0.d0 313 endif 314 endif 315 ! 316 ! displacements for 2nd order static and modal theory 317 ! 318 if(((iperturb(1).eq.1).or.(iperturb(2).eq.1)).and. 319 & (stiffness.eq.1).and.(buckling.eq.0)) then 320 write(*,*) '*ERROR in e_c3d_US3: no second order' 321 write(*,*) ' calculation for this type of element' 322 call exit(201) 323 endif 324 ! 325 if((buckling.eq.1).or.(coriolis.eq.1)) then 326 write(*,*) '*ERROR in e_c3d_US3: no buckling' 327 write(*,*) ' calculation for this type of element' 328 call exit(201) 329 endif 330 ! 331 ! initialization for the body forces 332 ! 333 if(rhsi.eq.1) then 334 if(nbody.ne.0) then 335 write(*,*) '*ERROR in e_c3d_u45: no body forces' 336 write(*,*) ' for this type of element' 337 call exit(201) 338 endif 339 endif 340 ! 341 if(buckling.eq.1) then 342 write(*,*) '*ERROR in e_c3d_u45: no buckling ' 343 write(*,*) ' calculation for this type of element' 344 call exit(201) 345 endif 346 ! 347 ! computation of the coordinates of the local nodes 348 ! 349 do i=1,nope 350 konl(i)=kon(indexe+i) 351 do j=1,3 352 xl(i,j) = co(j,konl(i)) 353 xg(i,j) = co(j,konl(i)) 354 enddo 355 enddo 356 ! 357 ! element frame & tranformation matrix tm (e1,e2,e3)T 358 call us4_csys(xg,tm,tmg) 359 ! nodal coordinates in element frame: 360 x(1,:) = matmul(tm,xg(1,:)) 361 x(2,:) = matmul(tm,xg(2,:)) 362 x(3,:) = matmul(tm,xg(3,:)) 363 x(4,:) = matmul(tm,xg(4,:)) 364 ! 365 Dm(:,:) = 0.d0 366 Db(:,:) = 0.d0 367 Ds(:,:) = 0.d0 368 ! 369 kode=nelcon(1,imat) 370 ! 371 if(kode.eq.2) then 372 mattyp=1 373 else 374 mattyp=0 375 endif 376 ! 377 ! material data and local stiffness matrix 378 ! 379 istiff=0 380 do jjj = 1,3 ! workaround for thermals 381 call materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon,nalcon, 382 & imat,amat,iorien,coords,orab,ntmat_,elas,rho, 383 & nelem,ithermal,alzero,mattyp,t0l,t1l, 384 & ihyper,istiff,elconloc,eth,kode,plicon, 385 & nplicon,plkcon,nplkcon,npmat_, 386 & plconloc,mi(1),dtime,kk, 387 & xstiff,ncmat_) 388 ! 389 e = xstiff(1,jjj,1) 390 un = xstiff(2,jjj,1) 391 rho = rhcon(1,1,imat) 392 ! 393 call us3_linel_Qi(e,un,Qin,Qs) 394 ! 395 Dm = Dm+Qin*gpthick(jjj,2)*dett 396 Db = Db+Qin*((gpthick(jjj,1))**2)*gpthick(jjj,2)*dettt 397 Ds = Ds+Qs*gpthick(jjj,2)*dett 398 enddo 399 ! 400 ! stiffness and mass matrix 401 ! 402 call us4_Kb(x,Db,Kb) 403 call us4_Ks(x,Ds,Ks) 404 call us4_Km(x,Dm,Kmem) 405 call us4_M(x,h,rho,Mshell) 406 ! 407 Kshell = Kmem + Kb + Ks 408 ! artifical drilling stiffness (Krotz) in orede to avoid singularities 409 kdmax = 0.d0 410 do k = 1,24 411 if(kdmax.LT.abs(Kshell(k,k))) then 412 kdmax = abs(Kshell(k,k)) 413 endif 414 enddo 415 Kshell(6,6) = kdmax/10000.d0 416 Kshell(12,12) = kdmax/10000.d0 417 Kshell(18,18) = kdmax/10000.d0 418 Kshell(24,24) = kdmax/10000.d0 419 ! 420 Kshell = matmul(matmul(transpose(tmg),Kshell),(tmg)) 421 Mshell = matmul(matmul(transpose(tmg),Mshell),(tmg)) 422 ! 423 do i=1,ndof*nope 424 do j=1,ndof*nope 425 s(i,j) = Kshell(i,j) 426 sm(i,j) = Mshell(i,j) 427 enddo 428 enddo 429 ! 430 return 431 end 432 433