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

Go to the source code of this file.

Functions/Subroutines

subroutine e_c3d_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, ieigenfrequency)
 

Function/Subroutine Documentation

◆ e_c3d_se()

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