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

Go to the source code of this file.

Functions/Subroutines

subroutine e_c3d_cs_se (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, distmin, ndesi, nodedesi, dfl, icoordinate, dxstiff, ne, xdesi, istartelem, ialelem, v, sigma, nk)
 

Function/Subroutine Documentation

◆ e_c3d_cs_se()

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