CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
e_c3d.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine e_c3d (co, kon, lakonl, p1, p2, omx, bodyfx, nbody, s, sm, ff, nelem, nmethod, elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, ielorien, norien, orab, ntmat_, t0, t1, ithermal, vold, iperturb, nelemload, sideload, xload, nload, idist, sti, stx, iexpl, plicon, nplicon, plkcon, nplkcon, xstiff, npmat_, dtime, matname, mi, ncmat_, mass, stiffness, buckling, rhsi, intscheme, ttime, time, istep, iinc, coriolis, xloadold, reltime, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, veold, springarea, nstate_, xstateini, xstate, ne0, ipkon, thicke, integerglob, doubleglob, tieset, istartset, iendset, ialset, ntie, nasym, pslavsurf, pmastsurf, mortar, clearini, ielprop, prop, kscale)
 

Function/Subroutine Documentation

◆ e_c3d()

subroutine e_c3d ( real*8, dimension(3,*), intent(in)  co,
integer, dimension(*), intent(in)  kon,
character*8, intent(in)  lakonl,
real*8, dimension(3), intent(in)  p1,
real*8, dimension(3), intent(in)  p2,
real*8, intent(in)  omx,
real*8, dimension(3), intent(in)  bodyfx,
integer, intent(in)  nbody,
real*8, dimension(60,60), intent(inout)  s,
real*8, dimension(60,60), intent(inout)  sm,
real*8, dimension(60), intent(inout)  ff,
integer, intent(in)  nelem,
integer, intent(inout)  nmethod,
real*8, dimension(0:ncmat_,ntmat_,*), intent(in)  elcon,
integer, dimension(2,*), intent(in)  nelcon,
real*8, dimension(0:1,ntmat_,*), intent(in)  rhcon,
integer, dimension(*), intent(in)  nrhcon,
real*8, dimension(0:6,ntmat_,*), intent(in)  alcon,
integer, dimension(2,*), intent(in)  nalcon,
real*8, dimension(*), intent(in)  alzero,
integer, dimension(mi(3),*), intent(in)  ielmat,
integer, dimension(mi(3),*), intent(in)  ielorien,
integer, intent(in)  norien,
real*8, dimension(7,*), intent(in)  orab,
integer, intent(in)  ntmat_,
real*8, dimension(*), intent(in)  t0,
real*8, dimension(*), intent(in)  t1,
integer, intent(in)  ithermal,
real*8, dimension(0:mi(2),*), intent(in)  vold,
integer, dimension(*), intent(in)  iperturb,
integer, dimension(2,*), intent(in)  nelemload,
character*20, dimension(*), intent(in)  sideload,
real*8, dimension(2,*), intent(inout)  xload,
integer, intent(in)  nload,
integer, intent(in)  idist,
real*8, dimension(6,mi(1),*), intent(in)  sti,
real*8, dimension(6,mi(1),*), intent(in)  stx,
integer, intent(in)  iexpl,
real*8, dimension(0:2*npmat_,ntmat_,*), intent(in)  plicon,
integer, dimension(0:ntmat_,*), intent(in)  nplicon,
real*8, dimension(0:2*npmat_,ntmat_,*), intent(in)  plkcon,
integer, dimension(0:ntmat_,*), intent(in)  nplkcon,
real*8, dimension(27,mi(1),*), intent(in)  xstiff,
integer, intent(in)  npmat_,
real*8, intent(in)  dtime,
character*80, dimension(*), intent(in)  matname,
integer, dimension(*), intent(in)  mi,
integer, intent(in)  ncmat_,
integer, intent(in)  mass,
integer, intent(in)  stiffness,
integer, intent(in)  buckling,
integer, intent(in)  rhsi,
integer, intent(in)  intscheme,
real*8, intent(in)  ttime,
real*8, intent(in)  time,
integer, intent(in)  istep,
integer, intent(in)  iinc,
integer, intent(in)  coriolis,
real*8, dimension(2,*), intent(in)  xloadold,
real*8, intent(in)  reltime,
integer, dimension(*), intent(in)  ipompc,
integer, dimension(3,*), intent(in)  nodempc,
real*8, dimension(*), intent(in)  coefmpc,
integer, intent(in)  nmpc,
integer, dimension(*), intent(in)  ikmpc,
integer, dimension(*), intent(in)  ilmpc,
real*8, dimension(0:mi(2),*), intent(in)  veold,
real*8, dimension(2,*), intent(inout)  springarea,
integer, intent(in)  nstate_,
real*8, dimension(nstate_,mi(1),*), intent(in)  xstateini,
real*8, dimension(nstate_,mi(1),*), intent(inout)  xstate,
integer, intent(in)  ne0,
integer, dimension(*), intent(in)  ipkon,
real*8, dimension(mi(3),*), intent(in)  thicke,
integer, dimension(*), intent(in)  integerglob,
real*8, dimension(*), intent(in)  doubleglob,
character*81, dimension(3,*), intent(in)  tieset,
integer, dimension(*), intent(in)  istartset,
integer, dimension(*), intent(in)  iendset,
integer, dimension(*), intent(in)  ialset,
integer, intent(in)  ntie,
integer, intent(in)  nasym,
real*8, dimension(3,*), intent(in)  pslavsurf,
real*8, dimension(6,*), intent(in)  pmastsurf,
integer, intent(in)  mortar,
real*8, dimension(3,9,*), intent(in)  clearini,
integer, dimension(*), intent(in)  ielprop,
real*8, dimension(*), intent(in)  prop,
integer  kscale 
)
31 !
32 ! computation of the element matrix and rhs for the element with
33 ! the topology in konl
34 !
35 ! ff: rhs without temperature and eigenstress contribution
36 !
37 ! nmethod=0: check for positive Jacobian
38 ! nmethod=1: stiffness matrix + right hand side
39 ! nmethod=2: stiffness matrix + mass matrix
40 ! nmethod=3: static stiffness + buckling stiffness
41 ! nmethod=4: stiffness matrix + mass matrix
42 !
43  implicit none
44 !
45  integer mass,stiffness,buckling,rhsi,coriolis
46 !
47  character*1 entity
48  character*8 lakonl
49  character*20 sideload(*)
50  character*80 matname(*),amat
51  character*81 tieset(3,*)
52 !
53  integer konl(26),ifaceq(8,6),nelemload(2,*),nbody,nelem,
54  & mi(*),iloc,jfaces,igauss,mortar,kon(*),ielprop(*),null,
55  & mattyp,ithermal,iperturb(*),nload,idist,i,j,k,l,i1,i2,j1,
56  & nmethod,k1,l1,ii,jj,ii1,jj1,id,ipointer,ig,m1,m2,m3,m4,kk,
57  & nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),six,
58  & ielorien(mi(3),*),ilayer,nlayer,ki,kl,ipkon(*),indexe,
59  & ntmat_,nope,nopes,norien,ihyper,iexpl,kode,imat,mint2d,
60  & mint3d,ifacet(6,4),nopev,iorien,istiff,ncmat_,iface,
61  & ifacew(8,5),intscheme,n,ipointeri,ipointerj,istep,iinc,
62  & layer,kspt,jltyp,iflag,iperm(60),m,ipompc(*),nodempc(3,*),
63  & nmpc,ikmpc(*),ilmpc(*),iscale,nstate_,ne0,iselect(6),kscale,
64  & istartset(*),iendset(*),ialset(*),ntie,integerglob(*),nasym,
65  & nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,nopered
66 !
67  real*8 co(3,*),xl(3,26),shp(4,26),xs2(3,7),veold(0:mi(2),*),
68  & s(60,60),w(3,3),p1(3),p2(3),bodyf(3),bodyfx(3),ff(60),
69  & bf(3),q(3),shpj(4,26),elcon(0:ncmat_,ntmat_,*),t(3),
70  & rhcon(0:1,ntmat_,*),xkl(3,3),eknlsign,reltime,prop(*),
71  & alcon(0:6,ntmat_,*),alzero(*),orab(7,*),t0(*),t1(*),
72  & anisox(3,3,3,3),voldl(0:mi(2),26),vo(3,3),xloadold(2,*),
73  & xl2(3,9),xsj2(3),shp2(7,9),vold(0:mi(2),*),xload(2,*),
74  & xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*),
75  & v(3,3,3,3),springarea(2,*),thickness,tlayer(4),dlayer(4),
76  & om,omx,e,un,al,um,xi,et,ze,tt,const,xsj,xsjj,sm(60,60),
77  & sti(6,mi(1),*),stx(6,mi(1),*),s11,s22,s33,s12,s13,s23,s11b,
78  & s22b,s33b,s12b,s13b,s23b,t0l,t1l,coefmpc(*),xlayer(mi(3),4),
79  & senergy,senergyb,rho,elas(21),summass,summ,thicke(mi(3),*),
80  & sume,factorm,factore,alp,elconloc(21),eth(6),doubleglob(*),
81  & weight,coords(3),dmass,xl1(3,9),term,clearini(3,9,*),
82  & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
83  & xstiff(27,mi(1),*),plconloc(802),dtime,ttime,time,tvar(2),
84  & sax(60,60),ffax(60),gs(8,4),a,stress(6),stre(3,3),
85  & pslavsurf(3,*),pmastsurf(6,*),xmass
86 !
87  intent(in) co,kon,lakonl,p1,p2,omx,bodyfx,nbody,
88  & nelem,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
89  & ielmat,ielorien,norien,orab,ntmat_,
90  & t0,t1,ithermal,vold,iperturb,nelemload,
91  & sideload,nload,idist,sti,stx,iexpl,plicon,
92  & nplicon,plkcon,nplkcon,xstiff,npmat_,dtime,
93  & matname,mi,ncmat_,mass,stiffness,buckling,rhsi,intscheme,
94  & ttime,time,istep,iinc,coriolis,xloadold,reltime,
95  & ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,veold,
96  & nstate_,xstateini,ne0,ipkon,thicke,
97  & integerglob,doubleglob,tieset,istartset,iendset,ialset,ntie,
98  & nasym,pslavsurf,pmastsurf,mortar,clearini,ielprop,prop
99 !
100  intent(inout) s,sm,ff,xload,nmethod,springarea,xstate
101 !
102  include "gauss.f"
103 !
104  ifaceq=reshape((/4,3,2,1,11,10,9,12,
105  & 5,6,7,8,13,14,15,16,
106  & 1,2,6,5,9,18,13,17,
107  & 2,3,7,6,10,19,14,18,
108  & 3,4,8,7,11,20,15,19,
109  & 4,1,5,8,12,17,16,20/),(/8,6/))
110  ifacet=reshape((/1,3,2,7,6,5,
111  & 1,2,4,5,9,8,
112  & 2,3,4,6,10,9,
113  & 1,4,3,8,10,7/),(/6,4/))
114  ifacew=reshape((/1,3,2,9,8,7,0,0,
115  & 4,5,6,10,11,12,0,0,
116  & 1,2,5,4,7,14,10,13,
117  & 2,3,6,5,8,15,11,14,
118  & 4,6,3,1,12,15,9,13/),(/8,5/))
119  iflag=3
120  null=0
121  iperm=(/13,14,-15,16,17,-18,19,20,-21,22,23,-24,
122  & 1,2,-3,4,5,-6,7,8,-9,10,11,-12,
123  & 37,38,-39,40,41,-42,43,44,-45,46,47,-48,
124  & 25,26,-27,28,29,-30,31,32,-33,34,35,-36,
125  & 49,50,-51,52,53,-54,55,56,-57,58,59,-60/)
126 !
127  tvar(1)=time
128  tvar(2)=ttime+time
129 !
130  summass=0.d0
131 !
132  indexe=ipkon(nelem)
133 c Bernhardi start
134  if(lakonl(1:5).eq.'C3D8I') then
135  nope=11
136  nopev=8
137  nopes=4
138  elseif(lakonl(4:5).eq.'20') then
139 c Bernhardi end
140  nope=20
141  nopev=8
142  nopes=8
143  elseif(lakonl(4:4).eq.'8') then
144  nope=8
145  nopev=8
146  nopes=4
147  elseif(lakonl(4:5).eq.'10') then
148  nope=10
149  nopev=4
150  nopes=6
151  elseif(lakonl(4:4).eq.'4') then
152  nope=4
153  nopev=4
154  nopes=3
155  elseif(lakonl(4:5).eq.'15') then
156  nope=15
157  nopev=6
158  elseif(lakonl(4:4).eq.'6') then
159  nope=6
160  nopev=6
161  elseif(lakonl(1:2).eq.'ES') then
162  if(lakonl(7:7).eq.'C') then
163  if(mortar.eq.0) then
164  nope=ichar(lakonl(8:8))-47
165  konl(nope+1)=kon(indexe+nope+1)
166  elseif(mortar.eq.1) then
167  nope=kon(indexe)
168  endif
169  else
170  nope=ichar(lakonl(8:8))-47
171  endif
172  elseif(lakonl(1:4).eq.'MASS') then
173  nope=1
174  endif
175 !
176 ! material and orientation
177 !
178  if(lakonl(7:8).ne.'LC') then
179 !
180  imat=ielmat(1,nelem)
181  amat=matname(imat)
182  if(norien.gt.0) then
183  iorien=ielorien(1,nelem)
184  else
185  iorien=0
186  endif
187 !
188  if(nelcon(1,imat).lt.0) then
189  ihyper=1
190  else
191  ihyper=0
192  endif
193  elseif(lakonl(4:5).eq.'20') then
194 !
195 ! composite materials
196 !
197 ! determining the number of layers
198 !
199  nlayer=0
200  do k=1,mi(3)
201  if(ielmat(k,nelem).ne.0) then
202  nlayer=nlayer+1
203  endif
204  enddo
205  mint2d=4
206 !
207 ! determining the layer thickness and global thickness
208 ! at the shell integration points
209 !
210  iflag=1
211  do kk=1,mint2d
212  xi=gauss3d2(1,kk)
213  et=gauss3d2(2,kk)
214  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
215  tlayer(kk)=0.d0
216  do i=1,nlayer
217  thickness=0.d0
218  do j=1,nopes
219  thickness=thickness+thicke(i,indexe+j)*shp2(4,j)
220  enddo
221  tlayer(kk)=tlayer(kk)+thickness
222  xlayer(i,kk)=thickness
223  enddo
224  enddo
225  iflag=3
226 !
227  ilayer=0
228  do i=1,4
229  dlayer(i)=0.d0
230  enddo
231 !
232  elseif(lakonl(4:5).eq.'15') then
233 !
234 ! composite materials
235 !
236 ! determining the number of layers
237 !
238  nlayer=0
239  do k=1,mi(3)
240  if(ielmat(k,nelem).ne.0) then
241  nlayer=nlayer+1
242  endif
243  enddo
244  mint2d=3
245 !
246 ! determining the layer thickness and global thickness
247 ! at the shell integration points
248 !
249  iflag=1
250  do kk=1,mint2d
251  xi=gauss3d10(1,kk)
252  et=gauss3d10(2,kk)
253  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
254  tlayer(kk)=0.d0
255  do i=1,nlayer
256  thickness=0.d0
257  do j=1,6
258  thickness=thickness+thicke(i,indexe+j)*shp2(4,j)
259  enddo
260  tlayer(kk)=tlayer(kk)+thickness
261  xlayer(i,kk)=thickness
262  enddo
263  enddo
264  iflag=3
265 !
266  ilayer=0
267  do i=1,3
268  dlayer(i)=0.d0
269  enddo
270 !
271  endif
272 !
273  if(intscheme.eq.0) then
274  if(lakonl(4:5).eq.'8R') then
275  mint2d=1
276  mint3d=1
277  elseif(lakonl(4:7).eq.'20RB') then
278  if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
279  mint2d=4
280  mint3d=50
281  else
282  mint2d=4
283  call beamintscheme(lakonl,mint3d,ielprop(nelem),prop,
284  & null,xi,et,ze,weight)
285  endif
286  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) then
287  if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'S').or.
288  & (lakonl(7:7).eq.'E')) then
289  mint2d=2
290  mint3d=4
291  else
292  mint2d=4
293  if(lakonl(7:8).eq.'LC') then
294  mint3d=8*nlayer
295  else
296  mint3d=8
297  endif
298  endif
299  elseif(lakonl(4:4).eq.'2') then
300  mint2d=9
301  mint3d=27
302  elseif(lakonl(4:5).eq.'10') then
303  mint2d=3
304  mint3d=4
305  elseif(lakonl(4:4).eq.'4') then
306  mint2d=1
307  mint3d=1
308  elseif(lakonl(4:5).eq.'15') then
309  if(lakonl(7:8).eq.'LC') then
310  mint3d=6*nlayer
311  else
312  mint3d=9
313  endif
314  elseif(lakonl(4:4).eq.'6') then
315  mint3d=2
316  else
317  mint3d=0
318  endif
319  else
320  if((lakonl(4:4).eq.'8').or.(lakonl(4:4).eq.'2')) then
321  mint3d=27
322  if(lakonl(4:4).eq.'8') then
323  if(lakonl(5:5).eq.'R') then
324  mint2d=1
325  else
326  mint2d=4
327  endif
328  else
329  if(lakonl(6:6).eq.'R') then
330  mint2d=4
331  else
332  mint2d=9
333  endif
334  endif
335  elseif((lakonl(4:5).eq.'10').or.(lakonl(4:4).eq.'4')) then
336  mint3d=15
337  if(lakonl(4:5).eq.'10') then
338  mint2d=3
339  else
340  mint2d=1
341  endif
342  elseif((lakonl(4:5).eq.'15').or.(lakonl(4:4).eq.'6')) then
343  mint3d=9
344  else
345  mint3d=0
346  endif
347  endif
348 !
349 ! computation of the coordinates of the local nodes
350 !
351  do i=1,nope
352  konl(i)=kon(indexe+i)
353  do j=1,3
354  xl(j,i)=co(j,konl(i))
355  enddo
356  enddo
357 !
358 ! initialisation for distributed forces
359 !
360  if(rhsi.eq.1) then
361  if(idist.ne.0) then
362  do i=1,3*nope
363  ff(i)=0.d0
364  enddo
365  endif
366  endif
367 !
368 ! displacements for 2nd order static and modal theory
369 !
370  if(((iperturb(1).eq.1).or.(iperturb(2).eq.1)).and.
371  & (stiffness.eq.1).and.(buckling.eq.0)) then
372  do i1=1,nope
373  do i2=1,3
374  voldl(i2,i1)=vold(i2,konl(i1))
375  enddo
376  enddo
377  endif
378 !
379 ! initialisation of sm
380 !
381  if((mass.eq.1).or.(buckling.eq.1).or.(coriolis.eq.1)) then
382  do i=1,3*nope
383  do j=1,3*nope
384  sm(i,j)=0.d0
385  enddo
386  enddo
387  endif
388 !
389 ! initialisation of s
390 !
391  do i=1,3*nope
392  do j=1,3*nope
393  s(i,j)=0.d0
394  enddo
395  enddo
396 !
397 ! calculating the stiffness matrix for the (contact) spring elements
398 ! and the mass matrix for the mass elements
399 !
400  if(mint3d.eq.0) then
401 !
402 ! for a
403 ! first step in a displacement-driven geometrically nonlinear
404 ! calculation or a nonlinear calculation with linear strains
405 ! the next block has to be evaluated
406 !
407  if(iperturb(2).eq.0) then
408  do i1=1,nope
409  do i2=1,3
410  voldl(i2,i1)=vold(i2,konl(i1))
411  enddo
412  enddo
413  endif
414 !
415  kode=nelcon(1,imat)
416  if(lakonl(1:2).eq.'ES') then
417  if(lakonl(7:7).ne.'C') then
418  t0l=0.d0
419  t1l=0.d0
420  if(ithermal.eq.1) then
421  t0l=(t0(konl(1))+t0(konl(2)))/2.d0
422  t1l=(t1(konl(1))+t1(konl(2)))/2.d0
423  elseif(ithermal.ge.2) then
424  t0l=(t0(konl(1))+t0(konl(2)))/2.d0
425  t1l=(vold(0,konl(1))+vold(0,konl(2)))/2.d0
426  endif
427  endif
428 !
429  if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'1').or.
430  & (lakonl(7:7).eq.'2').or.(mortar.eq.0)) then
431  call springstiff_n2f(xl,elas,konl,voldl,s,imat,elcon,
432  & nelcon,ncmat_,ntmat_,nope,lakonl,t1l,kode,elconloc,
433  & plicon,nplicon,npmat_,iperturb,
434  & springarea(1,konl(nope+1)),nmethod,mi,ne0,nstate_,
435  & xstateini,xstate,reltime,nasym,ielorien,orab,norien,
436  & nelem)
437  elseif(mortar.eq.1) then
438  iloc=kon(indexe+nope+1)
439  jfaces=kon(indexe+nope+2)
440  igauss=kon(indexe+nope+1)
441  call springstiff_f2f(xl,elas,voldl,s,imat,elcon,nelcon,
442  & ncmat_,ntmat_,nope,lakonl,t1l,kode,elconloc,plicon,
443  & nplicon,npmat_,iperturb,springarea(1,iloc),nmethod,
444  & mi,ne0,nstate_,xstateini,xstate,reltime,
445  & nasym,iloc,jfaces,igauss,pslavsurf,
446  & pmastsurf,clearini,kscale)
447  endif
448  elseif(lakonl(1:4).eq.'MASS') then
449 !
450 ! mass matrix contribution
451 !
452  if(mass.eq.1) then
453  do i1=1,3
454  sm(i1,i1)=rhcon(1,1,imat)
455  enddo
456  endif
457 !
458 ! contribution to the rhs
459 !
460  if(rhsi.eq.1) then
461  if(nbody.ne.0) then
462 !
463  xmass=rhcon(1,1,imat)
464 !
465  do i1=1,3
466  ff(i1)=bodyfx(i1)*xmass
467  enddo
468 !
469  if(omx.gt.0.d0) then
470 !
471 ! omega**2 * mass
472 !
473  om=omx*xmass
474  do i1=1,3
475 !
476 ! computation of the global coordinates of the
477 ! mass node
478 !
479  if((iperturb(1).ne.1).and.(iperturb(2).ne.1))
480  & then
481  q(i1)=xl(i1,1)
482  else
483  q(i1)=xl(i1,1)+voldl(i1,1)
484  endif
485 !
486  q(i1)=q(i1)-p1(i1)
487  enddo
488  const=q(1)*p2(1)+q(2)*p2(2)+q(3)*p2(3)
489 !
490 ! inclusion of the centrifugal force into the body force
491 !
492  do i1=1,3
493  ff(i1)=ff(i1)+(q(i1)-const*p2(i1))*om
494  enddo
495  endif
496  endif
497  endif
498 !
499  endif
500  return
501  endif
502 !
503 ! computation of the matrix: loop over the Gauss points
504 !
505  do kk=1,mint3d
506  if(intscheme.eq.0) then
507  if(lakonl(4:5).eq.'8R') then
508  xi=gauss3d1(1,kk)
509  et=gauss3d1(2,kk)
510  ze=gauss3d1(3,kk)
511  weight=weight3d1(kk)
512  elseif(lakonl(4:7).eq.'20RB') then
513  if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
514  xi=gauss3d13(1,kk)
515  et=gauss3d13(2,kk)
516  ze=gauss3d13(3,kk)
517  weight=weight3d13(kk)
518  else
519  call beamintscheme(lakonl,mint3d,ielprop(nelem),prop,
520  & kk,xi,et,ze,weight)
521  endif
522  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R'))
523  & then
524  if(lakonl(7:8).ne.'LC') then
525  xi=gauss3d2(1,kk)
526  et=gauss3d2(2,kk)
527  ze=gauss3d2(3,kk)
528  weight=weight3d2(kk)
529  else
530  kl=mod(kk,8)
531  if(kl.eq.0) kl=8
532 !
533  xi=gauss3d2(1,kl)
534  et=gauss3d2(2,kl)
535  ze=gauss3d2(3,kl)
536  weight=weight3d2(kl)
537 !
538  ki=mod(kk,4)
539  if(ki.eq.0) ki=4
540 !
541  if(kl.eq.1) then
542  ilayer=ilayer+1
543  if(ilayer.gt.1) then
544  do i=1,4
545  dlayer(i)=dlayer(i)+xlayer(ilayer-1,i)
546  enddo
547  endif
548  endif
549  ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/
550  & tlayer(ki)-1.d0
551  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
552 !
553 ! material and orientation
554 !
555  imat=ielmat(ilayer,nelem)
556  amat=matname(imat)
557  if(norien.gt.0) then
558  iorien=ielorien(ilayer,nelem)
559  else
560  iorien=0
561  endif
562 !
563  if(nelcon(1,imat).lt.0) then
564  ihyper=1
565  else
566  ihyper=0
567  endif
568  endif
569  elseif(lakonl(4:4).eq.'2') then
570  xi=gauss3d3(1,kk)
571  et=gauss3d3(2,kk)
572  ze=gauss3d3(3,kk)
573  weight=weight3d3(kk)
574  elseif(lakonl(4:5).eq.'10') then
575  xi=gauss3d5(1,kk)
576  et=gauss3d5(2,kk)
577  ze=gauss3d5(3,kk)
578  weight=weight3d5(kk)
579  elseif(lakonl(4:4).eq.'4') then
580  xi=gauss3d4(1,kk)
581  et=gauss3d4(2,kk)
582  ze=gauss3d4(3,kk)
583  weight=weight3d4(kk)
584  elseif(lakonl(4:5).eq.'15') then
585  if(lakonl(7:8).ne.'LC') then
586  xi=gauss3d8(1,kk)
587  et=gauss3d8(2,kk)
588  ze=gauss3d8(3,kk)
589  weight=weight3d8(kk)
590  else
591  kl=mod(kk,6)
592  if(kl.eq.0) kl=6
593 !
594  xi=gauss3d10(1,kl)
595  et=gauss3d10(2,kl)
596  ze=gauss3d10(3,kl)
597  weight=weight3d10(kl)
598 !
599  ki=mod(kk,3)
600  if(ki.eq.0) ki=3
601 !
602  if(kl.eq.1) then
603  ilayer=ilayer+1
604  if(ilayer.gt.1) then
605  do i=1,3
606  dlayer(i)=dlayer(i)+xlayer(ilayer-1,i)
607  enddo
608  endif
609  endif
610  ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/
611  & tlayer(ki)-1.d0
612  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
613 !
614 ! material and orientation
615 !
616  imat=ielmat(ilayer,nelem)
617  amat=matname(imat)
618  if(norien.gt.0) then
619  iorien=ielorien(ilayer,nelem)
620  else
621  iorien=0
622  endif
623 !
624  if(nelcon(1,imat).lt.0) then
625  ihyper=1
626  else
627  ihyper=0
628  endif
629  endif
630  elseif(lakonl(4:4).eq.'6') then
631  xi=gauss3d7(1,kk)
632  et=gauss3d7(2,kk)
633  ze=gauss3d7(3,kk)
634  weight=weight3d7(kk)
635  endif
636  else
637  if((lakonl(4:4).eq.'8').or.(lakonl(4:4).eq.'2')) then
638  xi=gauss3d3(1,kk)
639  et=gauss3d3(2,kk)
640  ze=gauss3d3(3,kk)
641  weight=weight3d3(kk)
642  elseif((lakonl(4:5).eq.'10').or.(lakonl(4:4).eq.'4')) then
643  xi=gauss3d6(1,kk)
644  et=gauss3d6(2,kk)
645  ze=gauss3d6(3,kk)
646  weight=weight3d6(kk)
647  else
648  xi=gauss3d8(1,kk)
649  et=gauss3d8(2,kk)
650  ze=gauss3d8(3,kk)
651  weight=weight3d8(kk)
652  endif
653  endif
654 c if(nelem.eq.1) then
655 c write(*,*) 'kk', kk
656 c write(*,*) 'coords',xi,et,ze
657 c write(*,*) 'weight',weight
658 c write(*,*) 'dlayer',dlayer(ki)
659 c endif
660 !
661 ! calculation of the shape functions and their derivatives
662 ! in the gauss point
663 !
664 c Bernhardi start
665  if(lakonl(1:5).eq.'C3D8R') then
666  call shape8hr(xl,xsj,shp,gs,a)
667  elseif(lakonl(1:5).eq.'C3D8I') then
668  call shape8hu(xi,et,ze,xl,xsj,shp,iflag)
669  elseif(nope.eq.20) then
670 c Bernhardi end
671  if(lakonl(7:7).eq.'A') then
672  call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag)
673  elseif((lakonl(7:7).eq.'E').or.(lakonl(7:7).eq.'S')) then
674  call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag)
675  else
676  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
677  endif
678  elseif(nope.eq.8) then
679  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
680  elseif(nope.eq.10) then
681  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
682  elseif(nope.eq.4) then
683  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
684  elseif(nope.eq.15) then
685  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
686  else
687  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
688  endif
689 !
690 ! check the jacobian determinant
691 !
692  if(xsj.lt.1.d-20) then
693  write(*,*) '*ERROR in e_c3d: nonpositive jacobian'
694  write(*,*) ' determinant in element',nelem
695  write(*,*)
696  xsj=dabs(xsj)
697  nmethod=0
698  endif
699 !
700 c if((iperturb(1).ne.0).and.stiffness.and.(.not.buckling))
701  if(((iperturb(1).eq.1).or.(iperturb(2).eq.1)).and.
702  & (stiffness.eq.1).and.(buckling.eq.0))then
703 !
704 ! stresses for 2nd order static and modal theory
705 !
706  s11=sti(1,kk,nelem)
707  s22=sti(2,kk,nelem)
708  s33=sti(3,kk,nelem)
709  s12=sti(4,kk,nelem)
710  s13=sti(5,kk,nelem)
711  s23=sti(6,kk,nelem)
712  endif
713 !
714 ! calculating the temperature in the integration
715 ! point
716 !
717  t0l=0.d0
718  t1l=0.d0
719  if(ithermal.eq.1) then
720  if(lakonl(4:5).eq.'8 ') then
721  do i1=1,nope
722  t0l=t0l+t0(konl(i1))/8.d0
723  t1l=t1l+t1(konl(i1))/8.d0
724  enddo
725  elseif(lakonl(4:6).eq.'20 ')then
726  nopered=20
727  call lintemp(t0,t1,konl,nopered,kk,t0l,t1l)
728  elseif(lakonl(4:6).eq.'10T') then
729  call linscal10(t0,konl,t0l,null,shp)
730  call linscal10(t1,konl,t1l,null,shp)
731  else
732  do i1=1,nope
733  t0l=t0l+shp(4,i1)*t0(konl(i1))
734  t1l=t1l+shp(4,i1)*t1(konl(i1))
735  enddo
736  endif
737  elseif(ithermal.ge.2) then
738  if(lakonl(4:5).eq.'8 ') then
739  do i1=1,nope
740  t0l=t0l+t0(konl(i1))/8.d0
741  t1l=t1l+vold(0,konl(i1))/8.d0
742  enddo
743  elseif(lakonl(4:6).eq.'20 ')then
744  nopered=20
745  call lintemp_th(t0,vold,konl,nopered,kk,t0l,t1l,mi)
746  elseif(lakonl(4:6).eq.'10T') then
747  call linscal10(t0,konl,t0l,null,shp)
748  call linscal10(vold,konl,t1l,mi(2),shp)
749  else
750  do i1=1,nope
751  t0l=t0l+shp(4,i1)*t0(konl(i1))
752  t1l=t1l+shp(4,i1)*vold(0,konl(i1))
753  enddo
754  endif
755  endif
756  tt=t1l-t0l
757 !
758 ! calculating the coordinates of the integration point
759 ! for material orientation purposes (for cylindrical
760 ! coordinate systems)
761 !
762  if(iorien.gt.0) then
763  do j=1,3
764  coords(j)=0.d0
765  do i1=1,nope
766  coords(j)=coords(j)+shp(4,i1)*co(j,konl(i1))
767  enddo
768  enddo
769  endif
770 !
771 ! for deformation plasticity: calculating the Jacobian
772 ! and the inverse of the deformation gradient
773 ! needed to convert the stiffness matrix in the spatial
774 ! frame of reference to the material frame
775 !
776  kode=nelcon(1,imat)
777 !
778 ! material data and local stiffness matrix
779 !
780  istiff=1
781  call materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
782  & imat,amat,iorien,coords,orab,ntmat_,elas,rho,
783  & nelem,ithermal,alzero,mattyp,t0l,t1l,
784  & ihyper,istiff,elconloc,eth,kode,plicon,
785  & nplicon,plkcon,nplkcon,npmat_,
786  & plconloc,mi(1),dtime,kk,
787  & xstiff,ncmat_)
788 !
789  if(mattyp.eq.1) then
790 c write(*,*) 'elastic co', elas(1),elas(2)
791  e=elas(1)
792  un=elas(2)
793  um=e/(1.d0+un)
794  al=un*um/(1.d0-2.d0*un)
795  um=um/2.d0
796  elseif(mattyp.eq.2) then
797 c call orthotropic(elas,anisox)
798  else
799  call anisotropic(elas,anisox)
800  endif
801 !
802 ! initialisation for the body forces
803 !
804  om=omx*rho
805  if(rhsi.eq.1) then
806  if(nbody.ne.0) then
807  do ii=1,3
808  bodyf(ii)=bodyfx(ii)*rho
809  enddo
810  endif
811  endif
812 !
813  if(buckling.eq.1) then
814 !
815 ! buckling stresses
816 !
817  s11b=stx(1,kk,nelem)
818  s22b=stx(2,kk,nelem)
819  s33b=stx(3,kk,nelem)
820  s12b=stx(4,kk,nelem)
821  s13b=stx(5,kk,nelem)
822  s23b=stx(6,kk,nelem)
823 !
824  endif
825 !
826 ! incorporating the jacobian determinant in the shape
827 ! functions
828 !
829  xsjj=dsqrt(xsj)
830  do i1=1,nope
831  shpj(1,i1)=shp(1,i1)*xsjj
832  shpj(2,i1)=shp(2,i1)*xsjj
833  shpj(3,i1)=shp(3,i1)*xsjj
834  shpj(4,i1)=shp(4,i1)*xsj
835  enddo
836 !
837 ! determination of the stiffness, and/or mass and/or
838 ! buckling matrix
839 !
840  if((stiffness.eq.1).or.(mass.eq.1).or.(buckling.eq.1).or.
841  & (coriolis.eq.1)) then
842 !
843  if(((iperturb(1).ne.1).and.(iperturb(2).ne.1)).or.
844  & (buckling.eq.1)) then
845 !
846 ! calculating the total mass of the element for
847 ! lumping purposes: only for explicit
848 ! dynamic calculations
849 !
850  if((mass.eq.1).and.(iexpl.gt.1)) then
851  summass=summass+rho*xsj
852  endif
853 !
854  jj1=1
855  do jj=1,nope
856 !
857  ii1=1
858  do ii=1,jj
859 !
860 ! all products of the shape functions for a given ii
861 ! and jj
862 !
863  do i1=1,3
864  do j1=1,3
865  w(i1,j1)=shpj(i1,ii)*shpj(j1,jj)
866  enddo
867  enddo
868 !
869 ! the following section calculates the static
870 ! part of the stiffness matrix which, for buckling
871 ! calculations, is done in a preliminary static
872 ! call
873 !
874  if(buckling.eq.0) then
875 !
876  if(mattyp.eq.1) then
877 !
878  s(ii1,jj1)=s(ii1,jj1)+(al*w(1,1)+
879  & um*(2.d0*w(1,1)+w(2,2)+w(3,3)))*weight
880 !
881  s(ii1,jj1+1)=s(ii1,jj1+1)+(al*w(1,2)+
882  & um*w(2,1))*weight
883  s(ii1,jj1+2)=s(ii1,jj1+2)+(al*w(1,3)+
884  & um*w(3,1))*weight
885  s(ii1+1,jj1)=s(ii1+1,jj1)+(al*w(2,1)+
886  & um*w(1,2))*weight
887  s(ii1+1,jj1+1)=s(ii1+1,jj1+1)+(al*w(2,2)+
888  & um*(2.d0*w(2,2)+w(1,1)+w(3,3)))*weight
889  s(ii1+1,jj1+2)=s(ii1+1,jj1+2)+(al*w(2,3)+
890  & um*w(3,2))*weight
891  s(ii1+2,jj1)=s(ii1+2,jj1)+(al*w(3,1)+
892  & um*w(1,3))*weight
893  s(ii1+2,jj1+1)=s(ii1+2,jj1+1)+(al*w(3,2)+
894  & um*w(2,3))*weight
895  s(ii1+2,jj1+2)=s(ii1+2,jj1+2)+(al*w(3,3)+
896  & um*(2.d0*w(3,3)+w(2,2)+w(1,1)))*weight
897 !
898  elseif(mattyp.eq.2) then
899 !
900  s(ii1,jj1)=s(ii1,jj1)+(elas(1)*w(1,1)+
901  & elas(7)*w(2,2)+elas(8)*w(3,3))*weight
902  s(ii1,jj1+1)=s(ii1,jj1+1)+(elas(2)*w(1,2)+
903  & elas(7)*w(2,1))*weight
904  s(ii1,jj1+2)=s(ii1,jj1+2)+(elas(4)*w(1,3)+
905  & elas(8)*w(3,1))*weight
906  s(ii1+1,jj1)=s(ii1+1,jj1)+(elas(7)*w(1,2)+
907  & elas(2)*w(2,1))*weight
908  s(ii1+1,jj1+1)=s(ii1+1,jj1+1)+
909  & (elas(7)*w(1,1)+
910  & elas(3)*w(2,2)+elas(9)*w(3,3))*weight
911  s(ii1+1,jj1+2)=s(ii1+1,jj1+2)+
912  & (elas(5)*w(2,3)+
913  & elas(9)*w(3,2))*weight
914  s(ii1+2,jj1)=s(ii1+2,jj1)+
915  & (elas(8)*w(1,3)+
916  & elas(4)*w(3,1))*weight
917  s(ii1+2,jj1+1)=s(ii1+2,jj1+1)+
918  & (elas(9)*w(2,3)+
919  & elas(5)*w(3,2))*weight
920  s(ii1+2,jj1+2)=s(ii1+2,jj1+2)+
921  & (elas(8)*w(1,1)+
922  & elas(9)*w(2,2)+elas(6)*w(3,3))*weight
923 !
924  else
925 !
926  do i1=1,3
927  do j1=1,3
928  do k1=1,3
929  do l1=1,3
930  s(ii1+i1-1,jj1+j1-1)=
931  & s(ii1+i1-1,jj1+j1-1)
932  & +anisox(i1,k1,j1,l1)
933  & *w(k1,l1)*weight
934  enddo
935  enddo
936  enddo
937  enddo
938 !
939  endif
940 !
941 ! mass matrix
942 !
943  if(mass.eq.1) then
944  sm(ii1,jj1)=sm(ii1,jj1)
945  & +rho*shpj(4,ii)*shp(4,jj)*weight
946  sm(ii1+1,jj1+1)=sm(ii1,jj1)
947  sm(ii1+2,jj1+2)=sm(ii1,jj1)
948  endif
949 !
950 ! Coriolis matrix
951 !
952  if(coriolis.eq.1) then
953  dmass=2.d0*
954  & rho*shpj(4,ii)*shp(4,jj)*weight*dsqrt(omx)
955  sm(ii1,jj1+1)=sm(ii1,jj1+1)-p2(3)*dmass
956  sm(ii1,jj1+2)=sm(ii1,jj1+2)+p2(2)*dmass
957  sm(ii1+1,jj1)=sm(ii1+1,jj1)+p2(3)*dmass
958  sm(ii1+1,jj1+2)=sm(ii1+1,jj1+2)-p2(1)*dmass
959  sm(ii1+2,jj1)=sm(ii1+2,jj1)-p2(2)*dmass
960  sm(ii1+2,jj1+1)=sm(ii1+2,jj1+1)+p2(1)*dmass
961  endif
962 !
963  else
964 !
965 ! buckling matrix
966 !
967  senergyb=
968  & (s11b*w(1,1)+s12b*(w(1,2)+w(2,1))
969  & +s13b*(w(1,3)+w(3,1))+s22b*w(2,2)
970  & +s23b*(w(2,3)+w(3,2))+s33b*w(3,3))*weight
971  sm(ii1,jj1)=sm(ii1,jj1)-senergyb
972  sm(ii1+1,jj1+1)=sm(ii1+1,jj1+1)-senergyb
973  sm(ii1+2,jj1+2)=sm(ii1+2,jj1+2)-senergyb
974 !
975  endif
976 !
977  ii1=ii1+3
978  enddo
979  jj1=jj1+3
980  enddo
981  else
982 !
983 ! stiffness matrix for static and modal
984 ! 2nd order calculations
985 !
986 ! large displacement stiffness
987 !
988  do i1=1,3
989  do j1=1,3
990  vo(i1,j1)=0.d0
991  do k1=1,nope
992  vo(i1,j1)=vo(i1,j1)+shp(j1,k1)*voldl(i1,k1)
993  enddo
994  enddo
995  enddo
996 !
997  if(mattyp.eq.1) then
998  call wcoef(v,vo,al,um)
999  endif
1000 !
1001 ! calculating the total mass of the element for
1002 ! lumping purposes: only for explicit
1003 ! dynamic calculations
1004 !
1005  if((mass.eq.1).and.(iexpl.gt.1)) then
1006  summass=summass+rho*xsj
1007  endif
1008 !
1009  jj1=1
1010  do jj=1,nope
1011 !
1012  ii1=1
1013  do ii=1,jj
1014 !
1015 ! all products of the shape functions for a given ii
1016 ! and jj
1017 !
1018  do i1=1,3
1019  do j1=1,3
1020  w(i1,j1)=shpj(i1,ii)*shpj(j1,jj)
1021  enddo
1022  enddo
1023 !
1024  if(mattyp.eq.1) then
1025 !
1026  do m1=1,3
1027  do m2=1,3
1028  do m3=1,3
1029  do m4=1,3
1030  s(ii1+m2-1,jj1+m1-1)=
1031  & s(ii1+m2-1,jj1+m1-1)
1032  & +v(m4,m3,m2,m1)*w(m4,m3)*weight
1033  enddo
1034  enddo
1035  enddo
1036  enddo
1037 !
1038  elseif(mattyp.eq.2) then
1039 !
1040  call orthonl(w,vo,elas,s,ii1,jj1,weight)
1041 !
1042  else
1043 !
1044  call anisonl(w,vo,elas,s,ii1,jj1,weight)
1045 !
1046  endif
1047 !
1048 ! stress stiffness
1049 !
1050  senergy=
1051  & (s11*w(1,1)+s12*(w(1,2)+w(2,1))
1052  & +s13*(w(1,3)+w(3,1))+s22*w(2,2)
1053  & +s23*(w(2,3)+w(3,2))+s33*w(3,3))*weight
1054  s(ii1,jj1)=s(ii1,jj1)+senergy
1055  s(ii1+1,jj1+1)=s(ii1+1,jj1+1)+senergy
1056  s(ii1+2,jj1+2)=s(ii1+2,jj1+2)+senergy
1057 !
1058 ! stiffness contribution of centrifugal forces
1059 !
1060  if((mass.eq.1).and.(om.gt.0.d0)) then
1061  dmass=shpj(4,ii)*shp(4,jj)*weight*om
1062  do m1=1,3
1063  s(ii1+m1-1,jj1+m1-1)=s(ii1+m1-1,jj1+m1-1)-
1064  & dmass
1065  do m2=1,3
1066  s(ii1+m1-1,jj1+m2-1)=s(ii1+m1-1,jj1+m2-1)+
1067  & dmass*p2(m1)*p2(m2)
1068  enddo
1069  enddo
1070  endif
1071 !
1072 ! mass matrix
1073 !
1074  if(mass.eq.1) then
1075  sm(ii1,jj1)=sm(ii1,jj1)
1076  & +rho*shpj(4,ii)*shp(4,jj)*weight
1077  sm(ii1+1,jj1+1)=sm(ii1,jj1)
1078  sm(ii1+2,jj1+2)=sm(ii1,jj1)
1079  endif
1080 !
1081 ! Coriolis matrix
1082 !
1083  if(coriolis.eq.1) then
1084  dmass=2.d0*
1085  & rho*shpj(4,ii)*shp(4,jj)*weight*dsqrt(omx)
1086  sm(ii1,jj1+1)=sm(ii1,jj1+1)-p2(3)*dmass
1087  sm(ii1,jj1+2)=sm(ii1,jj1+2)+p2(2)*dmass
1088  sm(ii1+1,jj1)=sm(ii1+1,jj1)+p2(3)*dmass
1089  sm(ii1+1,jj1+2)=sm(ii1+1,jj1+2)-p2(1)*dmass
1090  sm(ii1+2,jj1)=sm(ii1+2,jj1)-p2(2)*dmass
1091  sm(ii1+2,jj1+1)=sm(ii1+2,jj1+1)+p2(1)*dmass
1092  endif
1093 !
1094  ii1=ii1+3
1095  enddo
1096  jj1=jj1+3
1097  enddo
1098  endif
1099 !
1100  endif
1101 !
1102 ! add hourglass control stiffnesses: C3D8R only.
1103  if(lakonl(1:5).eq.'C3D8R') then
1104  call hgstiffness(s,elas,a,gs)
1105  endif
1106 !
1107 ! computation of the right hand side
1108 !
1109  if(rhsi.eq.1) then
1110 !
1111 ! distributed body flux
1112 !
1113  if(nload.gt.0) then
1114  call nident2(nelemload,nelem,nload,id)
1115  do
1116  if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
1117  if((sideload(id)(1:2).ne.'BX').and.
1118  & (sideload(id)(1:2).ne.'BY').and.
1119  & (sideload(id)(1:2).ne.'BZ')) then
1120  id=id-1
1121  cycle
1122  endif
1123  if(sideload(id)(3:4).eq.'NU') then
1124  do j=1,3
1125  coords(j)=0.d0
1126  do i1=1,nope
1127  coords(j)=coords(j)+
1128  & shp(4,i1)*xl(j,i1)
1129  enddo
1130  enddo
1131  if(sideload(id)(1:2).eq.'BX') then
1132  jltyp=1
1133  elseif(sideload(id)(1:2).eq.'BY') then
1134  jltyp=2
1135  elseif(sideload(id)(1:2).eq.'BZ') then
1136  jltyp=3
1137  endif
1138  iscale=1
1139  call dload(xload(1,id),istep,iinc,tvar,nelem,i,
1140  & layer,kspt,coords,jltyp,sideload(id),vold,co,
1141  & lakonl,konl,ipompc,nodempc,coefmpc,nmpc,ikmpc,
1142  & ilmpc,iscale,veold,rho,amat,mi)
1143  if((nmethod.eq.1).and.(iscale.ne.0))
1144  & xload(1,id)=xloadold(1,id)+
1145  & (xload(1,id)-xloadold(1,id))*reltime
1146  endif
1147  jj1=1
1148  do jj=1,nope
1149  if(sideload(id)(1:2).eq.'BX')
1150  & ff(jj1)=ff(jj1)+xload(1,id)*shpj(4,jj)*weight
1151  if(sideload(id)(1:2).eq.'BY')
1152  & ff(jj1+1)=ff(jj1+1)+xload(1,id)*shpj(4,jj)*weight
1153  if(sideload(id)(1:2).eq.'BZ')
1154  & ff(jj1+2)=ff(jj1+2)+xload(1,id)*shpj(4,jj)*weight
1155  jj1=jj1+3
1156  enddo
1157  id=id-1
1158  enddo
1159  endif
1160 !
1161 ! body forces
1162 !
1163  if(nbody.ne.0) then
1164  if(om.gt.0.d0) then
1165  do i1=1,3
1166 !
1167 ! computation of the global coordinates of the gauss
1168 ! point
1169 !
1170  q(i1)=0.d0
1171 c if(iperturb(1).eq.0) then
1172  if((iperturb(1).ne.1).and.(iperturb(2).ne.1)) then
1173  do j1=1,nope
1174  q(i1)=q(i1)+shp(4,j1)*xl(i1,j1)
1175  enddo
1176  else
1177  do j1=1,nope
1178  q(i1)=q(i1)+shp(4,j1)*
1179  & (xl(i1,j1)+voldl(i1,j1))
1180  enddo
1181  endif
1182 !
1183  q(i1)=q(i1)-p1(i1)
1184  enddo
1185  const=q(1)*p2(1)+q(2)*p2(2)+q(3)*p2(3)
1186 !
1187 ! inclusion of the centrifugal force into the body force
1188 !
1189  do i1=1,3
1190  bf(i1)=bodyf(i1)+(q(i1)-const*p2(i1))*om
1191  enddo
1192  else
1193  do i1=1,3
1194  bf(i1)=bodyf(i1)
1195  enddo
1196  endif
1197  jj1=1
1198  do jj=1,nope
1199  ff(jj1)=ff(jj1)+bf(1)*shpj(4,jj)*weight
1200  ff(jj1+1)=ff(jj1+1)+bf(2)*shpj(4,jj)*weight
1201  ff(jj1+2)=ff(jj1+2)+bf(3)*shpj(4,jj)*weight
1202  jj1=jj1+3
1203  enddo
1204  endif
1205 !
1206  endif
1207 !
1208  enddo
1209 !
1210 c write(*,*) nelem
1211 c write(*,'(6(1x,e11.4))') ((s(i1,j1),i1=1,j1),j1=1,60)
1212 c write(*,*)
1213 c
1214  if((buckling.eq.0).and.(nload.ne.0)) then
1215 !
1216 ! distributed loads
1217 !
1218  call nident2(nelemload,nelem,nload,id)
1219  do
1220  if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
1221  if(sideload(id)(1:1).ne.'P') then
1222  id=id-1
1223  cycle
1224  endif
1225 c read(sideload(id)(2:2),'(i1)') ig
1226  ig=ichar(sideload(id)(2:2))-48
1227 !
1228 !
1229 ! treatment of wedge faces
1230 !
1231  if(lakonl(4:4).eq.'6') then
1232  mint2d=1
1233  if(ig.le.2) then
1234  nopes=3
1235  else
1236  nopes=4
1237  endif
1238  endif
1239  if(lakonl(4:5).eq.'15') then
1240  if(ig.le.2) then
1241  mint2d=3
1242  nopes=6
1243  else
1244  mint2d=4
1245  nopes=8
1246  endif
1247  endif
1248 !
1249 c Bernhardi start
1250  if((nope.eq.20).or.(nope.eq.8).or.
1251  & (nope.eq.11)) then
1252 c Bernhardi end
1253 c if(iperturb(1).eq.0) then
1254  if((iperturb(1).ne.1).and.(iperturb(2).ne.1)) then
1255  do i=1,nopes
1256  do j=1,3
1257  xl2(j,i)=co(j,konl(ifaceq(i,ig)))
1258  enddo
1259  enddo
1260  else
1261  if(mass.eq.1) then
1262  do i=1,nopes
1263  do j=1,3
1264  xl1(j,i)=co(j,konl(ifaceq(i,ig)))
1265  enddo
1266  enddo
1267  endif
1268  do i=1,nopes
1269  do j=1,3
1270  xl2(j,i)=co(j,konl(ifaceq(i,ig)))+
1271  & vold(j,konl(ifaceq(i,ig)))
1272  enddo
1273  enddo
1274  endif
1275  elseif((nope.eq.10).or.(nope.eq.4)) then
1276 c if(iperturb(1).eq.0) then
1277  if((iperturb(1).ne.1).and.(iperturb(2).ne.1)) then
1278  do i=1,nopes
1279  do j=1,3
1280  xl2(j,i)=co(j,konl(ifacet(i,ig)))
1281  enddo
1282  enddo
1283  else
1284  if(mass.eq.1) then
1285  do i=1,nopes
1286  do j=1,3
1287  xl1(j,i)=co(j,konl(ifacet(i,ig)))
1288  enddo
1289  enddo
1290  endif
1291  do i=1,nopes
1292  do j=1,3
1293  xl2(j,i)=co(j,konl(ifacet(i,ig)))+
1294  & vold(j,konl(ifacet(i,ig)))
1295  enddo
1296  enddo
1297  endif
1298  else
1299 c if(iperturb(1).eq.0) then
1300  if((iperturb(1).ne.1).and.(iperturb(2).ne.1)) then
1301  do i=1,nopes
1302  do j=1,3
1303  xl2(j,i)=co(j,konl(ifacew(i,ig)))
1304  enddo
1305  enddo
1306  else
1307  if(mass.eq.1) then
1308  do i=1,nopes
1309  do j=1,3
1310  xl1(j,i)=co(j,konl(ifacew(i,ig)))
1311  enddo
1312  enddo
1313  endif
1314  do i=1,nopes
1315  do j=1,3
1316  xl2(j,i)=co(j,konl(ifacew(i,ig)))+
1317  & vold(j,konl(ifacew(i,ig)))
1318  enddo
1319  enddo
1320  endif
1321  endif
1322 !
1323  do i=1,mint2d
1324  if((lakonl(4:5).eq.'8R').or.
1325  & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
1326  xi=gauss2d1(1,i)
1327  et=gauss2d1(2,i)
1328  weight=weight2d1(i)
1329  elseif((lakonl(4:4).eq.'8').or.
1330  & (lakonl(4:6).eq.'20R').or.
1331  & ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
1332  xi=gauss2d2(1,i)
1333  et=gauss2d2(2,i)
1334  weight=weight2d2(i)
1335  elseif(lakonl(4:4).eq.'2') then
1336  xi=gauss2d3(1,i)
1337  et=gauss2d3(2,i)
1338  weight=weight2d3(i)
1339  elseif((lakonl(4:5).eq.'10').or.
1340  & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
1341  xi=gauss2d5(1,i)
1342  et=gauss2d5(2,i)
1343  weight=weight2d5(i)
1344  elseif((lakonl(4:4).eq.'4').or.
1345  & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
1346  xi=gauss2d4(1,i)
1347  et=gauss2d4(2,i)
1348  weight=weight2d4(i)
1349  endif
1350 !
1351  if(rhsi.eq.1) then
1352  if(nopes.eq.9) then
1353  call shape9q(xi,et,xl2,xsj2,xs2,shp2,iflag)
1354  elseif(nopes.eq.8) then
1355  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
1356  elseif(nopes.eq.4) then
1357  call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
1358  elseif(nopes.eq.6) then
1359  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
1360  elseif(nopes.eq.7) then
1361  call shape7tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
1362  else
1363  call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
1364  endif
1365 !
1366 ! for nonuniform load: determine the coordinates of the
1367 ! point (transferred into the user subroutine)
1368 !
1369  if(sideload(id)(3:4).eq.'NU') then
1370  do k=1,3
1371  coords(k)=0.d0
1372  do j=1,nopes
1373  coords(k)=coords(k)+xl2(k,j)*shp2(4,j)
1374  enddo
1375  enddo
1376 c read(sideload(id)(2:2),'(i1)') jltyp
1377  jltyp=ichar(sideload(id)(2:2))-48
1378  jltyp=jltyp+20
1379  iscale=1
1380  call dload(xload(1,id),istep,iinc,tvar,nelem,i,layer,
1381  & kspt,coords,jltyp,sideload(id),vold,co,lakonl,
1382  & konl,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,
1383  & iscale,veold,rho,amat,mi)
1384  if((nmethod.eq.1).and.(iscale.ne.0))
1385  & xload(1,id)=xloadold(1,id)+
1386  & (xload(1,id)-xloadold(1,id))*reltime
1387  elseif(sideload(id)(3:4).eq.'SM') then
1388 !
1389 ! submodel boundary: interpolation from the
1390 ! global model
1391 !
1392  do k=1,3
1393  coords(k)=0.d0
1394  do j=1,nopes
1395  coords(k)=coords(k)+xl2(k,j)*shp2(4,j)
1396  enddo
1397  enddo
1398 c read(sideload(id)(2:2),'(i1)') jltyp
1399  jltyp=ichar(sideload(id)(2:2))-48
1400 !
1401  entity='T'
1402  six=6
1403  do k=1,6
1404  iselect(k)=k+4
1405  enddo
1406  iface=10*nelem+jltyp
1407  call interpolsubmodel(integerglob,doubleglob,stress,
1408  & coords,iselect,six,iface,tieset,istartset,
1409  & iendset,ialset,ntie,entity)
1410 c write(*,*) 'e_c3d ',(stress(k),k=1,6)
1411 !
1412 ! cave: stress order of cgx: xx,yy,zz,xy,yz,xz
1413 !
1414  t(1)=stress(1)*xsj2(1)+stress(4)*xsj2(2)+
1415  & stress(6)*xsj2(3)
1416  t(2)=stress(4)*xsj2(1)+stress(2)*xsj2(2)+
1417  & stress(5)*xsj2(3)
1418  t(3)=stress(6)*xsj2(1)+stress(5)*xsj2(2)+
1419  & stress(3)*xsj2(3)
1420 !
1421  xload(1,id)=-1.d0
1422  do k=1,3
1423  xsj2(k)=t(k)
1424  enddo
1425  endif
1426 !
1427  do k=1,nopes
1428 c Bernhardi start
1429  if((nope.eq.20).or.(nope.eq.8).or.
1430  & (nope.eq.11)) then
1431 c Bernhardi end
1432  ipointer=(ifaceq(k,ig)-1)*3
1433  elseif((nope.eq.10).or.(nope.eq.4)) then
1434  ipointer=(ifacet(k,ig)-1)*3
1435  else
1436  ipointer=(ifacew(k,ig)-1)*3
1437  endif
1438  ff(ipointer+1)=ff(ipointer+1)-shp2(4,k)*xload(1,id)
1439  & *xsj2(1)*weight
1440  ff(ipointer+2)=ff(ipointer+2)-shp2(4,k)*xload(1,id)
1441  & *xsj2(2)*weight
1442  ff(ipointer+3)=ff(ipointer+3)-shp2(4,k)*xload(1,id)
1443  & *xsj2(3)*weight
1444  enddo
1445 !
1446 ! stiffness< contribution of the distributed load
1447 ! reference: Dhondt G., The Finite Element Method for
1448 ! three-dimensional thermomechanical Applications,
1449 ! Wiley, 2004, p 153, eqn. (3.54).
1450 !
1451 c elseif((mass.eq.1).and.(iperturb(1).ne.0)) then
1452  elseif((mass.eq.1).and.
1453  & ((iperturb(1).eq.1).or.(iperturb(2).eq.1))) then
1454  if(nopes.eq.9) then
1455  call shape9q(xi,et,xl1,xsj2,xs2,shp2,iflag)
1456  elseif(nopes.eq.8) then
1457  call shape8q(xi,et,xl1,xsj2,xs2,shp2,iflag)
1458  elseif(nopes.eq.4) then
1459  call shape4q(xi,et,xl1,xsj2,xs2,shp2,iflag)
1460  elseif(nopes.eq.6) then
1461  call shape6tri(xi,et,xl1,xsj2,xs2,shp2,iflag)
1462  elseif(nopes.eq.7) then
1463  call shape7tri(xi,et,xl1,xsj2,xs2,shp2,iflag)
1464  else
1465  call shape3tri(xi,et,xl1,xsj2,xs2,shp2,iflag)
1466  endif
1467 !
1468 ! for nonuniform load: determine the coordinates of the
1469 ! point (transferred into the user subroutine)
1470 !
1471  if(sideload(id)(3:4).eq.'NU') then
1472  do k=1,3
1473  coords(k)=0.d0
1474  do j=1,nopes
1475  coords(k)=coords(k)+xl1(k,j)*shp2(4,j)
1476  enddo
1477  enddo
1478 c read(sideload(id)(2:2),'(i1)') jltyp
1479  jltyp=ichar(sideload(id)(2:2))-48
1480  jltyp=jltyp+20
1481  iscale=1
1482  call dload(xload(1,id),istep,iinc,tvar,nelem,i,layer,
1483  & kspt,coords,jltyp,sideload(id),vold,co,lakonl,
1484  & konl,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,
1485  & iscale,veold,rho,amat,mi)
1486  if((nmethod.eq.1).and.(iscale.ne.0))
1487  & xload(1,id)=xloadold(1,id)+
1488  & (xload(1,id)-xloadold(1,id))*reltime
1489  elseif(sideload(id)(3:4).eq.'SM') then
1490 !
1491 ! submodel boundary: interpolation from the
1492 ! global model
1493 !
1494  do k=1,3
1495  coords(k)=0.d0
1496  do j=1,nopes
1497  coords(k)=coords(k)+xl2(k,j)*shp2(4,j)
1498  enddo
1499  enddo
1500 c read(sideload(id)(2:2),'(i1)') jltyp
1501  jltyp=ichar(sideload(id)(2:2))-48
1502 !
1503  entity='T'
1504  six=6
1505  do k=1,6
1506  iselect(k)=k+4
1507  enddo
1508  iface=10*nelem+jltyp
1509  call interpolsubmodel(integerglob,doubleglob,stress,
1510  & coords,iselect,six,iface,tieset,istartset,
1511  & iendset,ialset,ntie,entity)
1512 c write(*,*) 'e_c3d ',(stress(k),k=1,6)
1513 !
1514  stre(1,1)=stress(1)
1515  stre(1,2)=stress(4)
1516  stre(1,3)=stress(6)
1517  stre(2,1)=stress(4)
1518  stre(2,2)=stress(2)
1519  stre(2,3)=stress(5)
1520  stre(3,1)=stress(6)
1521  stre(3,2)=stress(5)
1522  stre(3,3)=stress(3)
1523  endif
1524 !
1525 ! calculation of the deformation gradient
1526 !
1527  do k=1,3
1528  do l=1,3
1529  xkl(k,l)=0.d0
1530  do ii=1,nopes
1531  xkl(k,l)=xkl(k,l)+shp2(l,ii)*xl2(k,ii)
1532  enddo
1533  enddo
1534  enddo
1535 !
1536  do ii=1,nopes
1537 c Bernhardi start
1538  if((nope.eq.20).or.(nope.eq.8).or.
1539  & (nope.eq.11)) then
1540 c Bernhardi end
1541  ipointeri=(ifaceq(ii,ig)-1)*3
1542  elseif((nope.eq.10).or.(nope.eq.4))then
1543  ipointeri=(ifacet(ii,ig)-1)*3
1544  else
1545  ipointeri=(ifacew(ii,ig)-1)*3
1546  endif
1547  do jj=1,nopes
1548 c Bernhardi start
1549  if((nope.eq.20).or.(nope.eq.8)
1550  & .or.(nope.eq.11)) then
1551 c Bernhardi end
1552  ipointerj=(ifaceq(jj,ig)-1)*3
1553  elseif((nope.eq.10).or.(nope.eq.4)) then
1554  ipointerj=(ifacet(jj,ig)-1)*3
1555  else
1556  ipointerj=(ifacew(jj,ig)-1)*3
1557  endif
1558 !
1559 ! if no submodel: only pressure
1560 ! else: complete stress vector
1561 !
1562  if(sideload(id)(3:4).ne.'SM') then
1563  do k=1,3
1564  do l=1,3
1565  if(k.eq.l) cycle
1566  eknlsign=1.d0
1567  if(k*l.eq.2) then
1568  n=3
1569  if(k.lt.l) eknlsign=-1.d0
1570  elseif(k*l.eq.3) then
1571  n=2
1572  if(k.gt.l) eknlsign=-1.d0
1573  else
1574  n=1
1575  if(k.lt.l) eknlsign=-1.d0
1576  endif
1577  term=weight*xload(1,id)*shp2(4,ii)*
1578  & eknlsign*(xsj2(1)*
1579  & (xkl(n,2)*shp2(3,jj)-xkl(n,3)*
1580  & shp2(2,jj))+xsj2(2)*
1581  & (xkl(n,3)*shp2(1,jj)-xkl(n,1)*
1582  & shp2(3,jj))+xsj2(3)*
1583  & (xkl(n,1)*shp2(2,jj)-xkl(n,2)*
1584  & shp2(1,jj)))
1585  s(ipointeri+k,ipointerj+l)=
1586  & s(ipointeri+k,ipointerj+l)+term/2.d0
1587  s(ipointerj+l,ipointeri+k)=
1588  & s(ipointerj+l,ipointeri+k)+term/2.d0
1589  enddo
1590  enddo
1591  else
1592  do kk=1,3
1593  do k=1,3
1594  do l=1,3
1595  if(k.eq.l) cycle
1596  eknlsign=1.d0
1597  if(k*l.eq.2) then
1598  n=3
1599  if(k.lt.l) eknlsign=-1.d0
1600  elseif(k*l.eq.3) then
1601  n=2
1602  if(k.gt.l) eknlsign=-1.d0
1603  else
1604  n=1
1605  if(k.lt.l) eknlsign=-1.d0
1606  endif
1607  term=-weight*stre(kk,k)*shp2(4,ii)*
1608  & eknlsign*(xsj2(1)*
1609  & (xkl(n,2)*shp2(3,jj)-xkl(n,3)*
1610  & shp2(2,jj))+xsj2(2)*
1611  & (xkl(n,3)*shp2(1,jj)-xkl(n,1)*
1612  & shp2(3,jj))+xsj2(3)*
1613  & (xkl(n,1)*shp2(2,jj)-xkl(n,2)*
1614  & shp2(1,jj)))
1615  s(ipointeri+kk,ipointerj+l)=
1616  & s(ipointeri+kk,ipointerj+l)+term/2.d0
1617  s(ipointerj+l,ipointeri+kk)=
1618  & s(ipointerj+l,ipointeri+kk)+term/2.d0
1619  enddo
1620  enddo
1621  enddo
1622  endif
1623  enddo
1624  enddo
1625 !
1626  endif
1627  enddo
1628 !
1629  id=id-1
1630  enddo
1631  endif
1632 !
1633 ! for axially symmetric and plane stress/strain elements:
1634 ! complete s and sm
1635 !
1636  if(((lakonl(4:5).eq.'8 ').or.
1637  & ((lakonl(4:6).eq.'20R').and.(lakonl(7:8).ne.'BR'))).and.
1638  & ((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'S').or.
1639  & (lakonl(7:7).eq.'E'))) then
1640  do i=1,60
1641  do j=i,60
1642  k=abs(iperm(i))
1643  l=abs(iperm(j))
1644  if(k.gt.l) then
1645  m=k
1646  k=l
1647  l=m
1648  endif
1649  sax(i,j)=s(k,l)*iperm(i)*iperm(j)/(k*l)
1650  enddo
1651  enddo
1652  do i=1,60
1653  do j=i,60
1654  s(i,j)=s(i,j)+sax(i,j)
1655  enddo
1656  enddo
1657 !
1658  if((nload.ne.0).or.(nbody.ne.0)) then
1659  do i=1,60
1660  k=abs(iperm(i))
1661  ffax(i)=ff(k)*iperm(i)/k
1662  enddo
1663  do i=1,60
1664  ff(i)=ff(i)+ffax(i)
1665  enddo
1666  endif
1667 !
1668  if(mass.eq.1) then
1669  summass=2.d0*summass
1670  do i=1,60
1671  do j=i,60
1672  k=abs(iperm(i))
1673  l=abs(iperm(j))
1674  if(k.gt.l) then
1675  m=k
1676  k=l
1677  l=m
1678  endif
1679  sax(i,j)=sm(k,l)*iperm(i)*iperm(j)/(k*l)
1680  enddo
1681  enddo
1682  do i=1,60
1683  do j=i,60
1684  sm(i,j)=sm(i,j)+sax(i,j)
1685  enddo
1686  enddo
1687  endif
1688  endif
1689 !
1690  if((mass.eq.1).and.(iexpl.gt.1)) then
1691 !
1692 ! scaling the diagonal terms of the mass matrix such that the total mass
1693 ! is right (LUMPING; for explicit dynamic calculations)
1694 !
1695  sume=0.d0
1696  summ=0.d0
1697  do i=1,3*nopev,3
1698  sume=sume+sm(i,i)
1699  enddo
1700  do i=3*nopev+1,3*nope,3
1701  summ=summ+sm(i,i)
1702  enddo
1703 !
1704  if(nope.eq.20) then
1705 c alp=.2215d0
1706  alp=.2917d0
1707 ! maybe alp=.2917d0 is better??
1708  elseif(nope.eq.10) then
1709  alp=0.1203d0
1710  elseif(nope.eq.15) then
1711  alp=0.2141d0
1712  endif
1713 !
1714  if((nope.eq.20).or.(nope.eq.10).or.
1715  & (nope.eq.15)) then
1716  factore=summass*alp/(1.d0+alp)/sume
1717  factorm=summass/(1.d0+alp)/summ
1718  else
1719  factore=summass/sume
1720  endif
1721 !
1722  do i=1,3*nopev,3
1723  sm(i,i)=sm(i,i)*factore
1724  sm(i+1,i+1)=sm(i,i)
1725  sm(i+2,i+2)=sm(i,i)
1726  enddo
1727  do i=3*nopev+1,3*nope,3
1728  sm(i,i)=sm(i,i)*factorm
1729  sm(i+1,i+1)=sm(i,i)
1730  sm(i+2,i+2)=sm(i,i)
1731  enddo
1732 !
1733  endif
1734 !
1735  return
shape20h_pl
subroutine shape20h_pl(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape20h_pl.f:20
shape6w
subroutine shape6w(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape6w.f:20
shape9q
subroutine shape9q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape9q.f:20
linscal10
subroutine linscal10(scal, konl, scall, idim, shp)
Definition: linscal10.f:20
shape8q
subroutine shape8q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape8q.f:20
shape3tri
subroutine shape3tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape3tri.f:20
shape7tri
subroutine shape7tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape7tri.f:20
shape10tet
subroutine shape10tet(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape10tet.f:20
shape8h
subroutine shape8h(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape8h.f:20
nident2
subroutine nident2(x, px, n, id)
Definition: nident2.f:27
anisonl
subroutine anisonl(w, vo, elas, s, ii1, jj1, weight)
Definition: anisonl.f:20
shape15w
subroutine shape15w(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape15w.f:20
shape20h_ax
subroutine shape20h_ax(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape20h_ax.f:20
hgstiffness
subroutine hgstiffness(s, elas, a, gs)
Definition: hgstiffness.f:20
idist
static ITG * idist
Definition: radflowload.c:39
orthonl
subroutine orthonl(w, vo, elas, s, ii1, jj1, weight)
Definition: orthonl.f:20
lintemp_th
subroutine lintemp_th(t0, vold, konl, nope, jj, t0l, t1l, mi)
Definition: lintemp_th.f:20
shape4q
subroutine shape4q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape4q.f:20
materialdata_me
subroutine materialdata_me(elcon, nelcon, rhcon, nrhcon, alcon, nalcon, imat, amat, iorien, pgauss, orab, ntmat_, elas, rho, iel, ithermal, alzero, mattyp, t0l, t1l, ihyper, istiff, elconloc, eth, kode, plicon, nplicon, plkcon, nplkcon, npmat_, plconloc, mi, dtime, iint, xstiff, ncmat_)
Definition: materialdata_me.f:23
lintemp
subroutine lintemp(t0, t1, konl, nope, jj, t0l, t1l)
Definition: lintemp.f:20
shape20h
subroutine shape20h(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape20h.f:20
shape8hr
subroutine shape8hr(xl, xsj, shp, gs, a)
Definition: shape8hr.f:20
beamintscheme
subroutine beamintscheme(lakonl, mint3d, npropstart, prop, kk, xi, et, ze, weight)
Definition: beamintscheme.f:21
shape4tet
subroutine shape4tet(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape4tet.f:20
shape8hu
subroutine shape8hu(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape8hu.f:20
shape6tri
subroutine shape6tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape6tri.f:20
anisotropic
subroutine anisotropic(anisol, anisox)
Definition: anisotropic.f:20
springstiff_n2f
subroutine springstiff_n2f(xl, elas, konl, voldl, s, imat, elcon, nelcon, ncmat_, ntmat_, nope, lakonl, t1l, kode, elconloc, plicon, nplicon, npmat_, iperturb, springarea, nmethod, mi, ne0, nstate_, xstateini, xstate, reltime, nasym, ielorien, orab, norien, nelem)
Definition: springstiff_n2f.f:24
springstiff_f2f
subroutine springstiff_f2f(xl, elas, voldl, s, imat, elcon, nelcon, ncmat_, ntmat_, nope, lakonl, t1l, kode, elconloc, plicon, nplicon, npmat_, iperturb, springarea, nmethod, mi, ne0, nstate_, xstateini, xstate, reltime, nasym, iloc, jfaces, igauss, pslavsurf, pmastsurf, clearini, kscale)
Definition: springstiff_f2f.f:24
thickness
subroutine thickness(dgdx, nobject, nodedesiboun, ndesiboun, objectset, xo, yo, zo, x, y, z, nx, ny, nz, co, ifree, ndesia, ndesib, iobject, ndesi, dgdxglob, nk)
Definition: thickness.f:22
wcoef
subroutine wcoef(v, vo, al, um)
Definition: wcoef.f:20
interpolsubmodel
subroutine interpolsubmodel(integerglob, doubleglob, value, coords, iselect, nselect, nodeface, tieset, istartset, iendset, ialset, ntie, entity)
Definition: interpolsubmodel.f:22
dload
subroutine dload(f, kstep, kinc, time, noel, npt, layer, kspt, coords, jltyp, loadtype, vold, co, lakonl, konl, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, iscale, veold, rho, amat, mi)
Definition: dload.f:23
Hosted by OpenAircraft.com, (Michigan UAV, LLC)