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

Go to the source code of this file.

Functions/Subroutines

subroutine resultsmech (co, kon, ipkon, lakon, ne, v, stx, elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, ielorien, norien, orab, ntmat_, t0, t1, ithermal, prestr, iprestr, eme, iperturb, fn, iout, qa, vold, nmethod, veold, dtime, time, ttime, plicon, nplicon, plkcon, nplkcon, xstateini, xstiff, xstate, npmat_, matname, mi, ielas, icmd, ncmat_, nstate_, stiini, vini, ener, eei, enerini, istep, iinc, springarea, reltime, calcul_fn, calcul_qa, calcul_cauchy, nener, ikin, nal, ne0, thicke, emeini, pslavsurf, pmastsurf, mortar, clearini, nea, neb, ielprop, prop, kscale)
 

Function/Subroutine Documentation

◆ resultsmech()

subroutine resultsmech ( real*8, dimension(3,*), intent(in)  co,
integer, dimension(*), intent(in)  kon,
integer, dimension(*), intent(in)  ipkon,
character*8, dimension(*), intent(in)  lakon,
integer, intent(in)  ne,
real*8, dimension(0:mi(2),*), intent(in)  v,
real*8, dimension(6,mi(1),*), intent(inout)  stx,
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, dimension(2), intent(in)  ithermal,
real*8, dimension(6,mi(1),*), intent(in)  prestr,
integer, intent(in)  iprestr,
real*8, dimension(6,mi(1),*), intent(inout)  eme,
integer, dimension(*), intent(in)  iperturb,
real*8, dimension(0:mi(2),*), intent(inout)  fn,
integer, intent(in)  iout,
real*8, dimension(*), intent(inout)  qa,
real*8, dimension(0:mi(2),*), intent(in)  vold,
integer, intent(in)  nmethod,
real*8, dimension(0:mi(2),*), intent(in)  veold,
real*8, intent(in)  dtime,
real*8, intent(in)  time,
real*8, intent(in)  ttime,
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(nstate_,mi(1),*), intent(in)  xstateini,
real*8, dimension(27,mi(1),*), intent(inout)  xstiff,
real*8, dimension(nstate_,mi(1),*), intent(in)  xstate,
integer, intent(in)  npmat_,
character*80, dimension(*), intent(in)  matname,
integer, dimension(*), intent(in)  mi,
integer, intent(in)  ielas,
integer, intent(in)  icmd,
integer, intent(in)  ncmat_,
integer, intent(in)  nstate_,
real*8, dimension(6,mi(1),*), intent(in)  stiini,
real*8, dimension(0:mi(2),*), intent(in)  vini,
real*8, dimension(mi(1),*), intent(inout)  ener,
real*8, dimension(6,mi(1),*), intent(inout)  eei,
real*8, dimension(mi(1),*), intent(in)  enerini,
integer, intent(in)  istep,
integer, intent(in)  iinc,
real*8, dimension(2,*), intent(in)  springarea,
real*8, intent(in)  reltime,
integer, intent(in)  calcul_fn,
integer, intent(in)  calcul_qa,
integer, intent(in)  calcul_cauchy,
integer, intent(in)  nener,
integer, intent(in)  ikin,
integer, intent(inout)  nal,
integer, intent(in)  ne0,
real*8, dimension(mi(3),*), intent(in)  thicke,
real*8, dimension(6,mi(1),*), intent(in)  emeini,
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, intent(in)  nea,
integer, intent(in)  neb,
integer, dimension(*), intent(in)  ielprop,
real*8, dimension(*), intent(in)  prop,
integer, intent(in)  kscale 
)
29 !
30 ! calculates stresses and the material tangent at the integration
31 ! points and the internal forces at the nodes
32 !
33  implicit none
34 !
35  integer cauchy
36 !
37  character*8 lakon(*),lakonl
38  character*80 amat,matname(*)
39 !
40  integer kon(*),konl(26),nea,neb,mi(*),mint2d,nopes,
41  & nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),
42  & ielorien(mi(3),*),ntmat_,ipkon(*),ne0,iflag,null,kscale,
43  & istep,iinc,mt,ne,mattyp,ithermal(2),iprestr,i,j,k,m1,m2,jj,
44  & i1,m3,m4,kk,nener,indexe,nope,norien,iperturb(*),iout,
45  & nal,icmd,ihyper,nmethod,kode,imat,mint3d,iorien,ielas,
46  & istiff,ncmat_,nstate_,ikin,ilayer,nlayer,ki,kl,ielprop(*),
47  & nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,calcul_fn,
48  & calcul_cauchy,calcul_qa,nopered,mortar,jfaces,iloc,igauss
49 !
50  real*8 co(3,*),v(0:mi(2),*),shp(4,26),stiini(6,mi(1),*),
51  & stx(6,mi(1),*),xl(3,26),vl(0:mi(2),26),stre(6),prop(*),
52  & elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),xs2(3,7),
53  & alcon(0:6,ntmat_,*),vini(0:mi(2),*),thickness,
54  & alzero(*),orab(7,*),elas(21),rho,fn(0:mi(2),*),
55  & fnl(3,10),skl(3,3),beta(6),q(0:mi(2),26),xl2(3,8),
56  & vkl(0:3,3),t0(*),t1(*),prestr(6,mi(1),*),eme(6,mi(1),*),
57  & ckl(3,3),vold(0:mi(2),*),eloc(9),veold(0:mi(2),*),
58  & springarea(2,*),elconloc(21),eth(6),xkl(3,3),voldl(0:mi(2),26),
59  & xikl(3,3),ener(mi(1),*),emec(6),eei(6,mi(1),*),enerini(mi(1),*),
60  & emec0(6),vel(1:3,26),veoldl(0:mi(2),26),xsj2(3),shp2(7,8),
61  & e,un,al,um,am1,xi,et,ze,tt,exx,eyy,ezz,exy,exz,eyz,
62  & xsj,qa(*),vj,t0l,t1l,dtime,weight,pgauss(3),vij,time,ttime,
63  & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
64  & xstiff(27,mi(1),*),xstate(nstate_,mi(1),*),plconloc(802),
65  & vokl(3,3),xstateini(nstate_,mi(1),*),vikl(3,3),
66  & gs(8,4),a,reltime,tlayer(4),dlayer(4),xlayer(mi(3),4),
67  & thicke(mi(3),*),emeini(6,mi(1),*),clearini(3,9,*),
68  & pslavsurf(3,*),pmastsurf(6,*)
69 !
70  intent(in) co,kon,ipkon,lakon,ne,v,
71  & elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
72  & ielmat,ielorien,norien,orab,ntmat_,t0,t1,ithermal,prestr,
73  & iprestr,iperturb,iout,vold,nmethod,
74  & veold,dtime,time,ttime,plicon,nplicon,plkcon,nplkcon,
75  & xstateini,xstate,npmat_,matname,mi,ielas,icmd,
76  & ncmat_,nstate_,stiini,vini,enerini,istep,iinc,
77  & springarea,reltime,calcul_fn,calcul_qa,calcul_cauchy,nener,
78  & ikin,ne0,thicke,emeini,pslavsurf,
79  & pmastsurf,mortar,clearini,nea,neb,ielprop,prop,kscale
80 !
81  intent(inout) nal,qa,fn,xstiff,ener,eme,eei,stx
82 !
83  include "gauss.f"
84 !
85  iflag=3
86  null=0
87 !
88  mt=mi(2)+1
89  nal=0
90  qa(3)=-1.
91  qa(4)=0.
92 !
93  do i=nea,neb
94 !
95  lakonl=lakon(i)
96 !
97  if(ipkon(i).lt.0) cycle
98 !
99 ! no 3D-fluid elements
100 !
101  if(lakonl(1:1).eq.'F') cycle
102  if(lakonl(1:7).eq.'DCOUP3D') cycle
103 !
104 ! user elements
105 !
106  if(lakonl(1:1).eq.'U') then
107  call resultsmech_u(co,kon,ipkon,lakon,ne,v,
108  & stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
109  & ielmat,ielorien,norien,orab,ntmat_,t0,t1,ithermal,prestr,
110  & iprestr,eme,iperturb,fn,iout,qa,vold,nmethod,
111  & veold,dtime,time,ttime,plicon,nplicon,plkcon,nplkcon,
112  & xstateini,xstiff,xstate,npmat_,matname,mi,ielas,icmd,
113  & ncmat_,nstate_,stiini,vini,ener,eei,enerini,istep,iinc,
114  & reltime,calcul_fn,calcul_qa,calcul_cauchy,nener,
115  & ikin,nal,ne0,thicke,emeini,i,ielprop,prop)
116  cycle
117  endif
118 !
119  if(lakonl(7:8).ne.'LC') then
120 !
121  imat=ielmat(1,i)
122  amat=matname(imat)
123  if(norien.gt.0) then
124  iorien=ielorien(1,i)
125  else
126  iorien=0
127  endif
128 !
129  if(nelcon(1,imat).lt.0) then
130  ihyper=1
131  else
132  ihyper=0
133  endif
134  elseif(lakonl(4:5).eq.'20') then
135 !
136 ! composite materials
137 ! S8R - composite element
138 !
139  mint2d=4
140  nopes=8
141 !
142 ! determining the number of layers
143 !
144  nlayer=0
145  do k=1,mi(3)
146  if(ielmat(k,i).ne.0) then
147  nlayer=nlayer+1
148  endif
149  enddo
150 !
151 ! determining the layer thickness and global thickness
152 ! at the shell integration points
153 !
154  iflag=1
155  indexe=ipkon(i)
156  do kk=1,mint2d
157  xi=gauss3d2(1,kk)
158  et=gauss3d2(2,kk)
159  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
160  tlayer(kk)=0.d0
161  do k=1,nlayer
162  thickness=0.d0
163  do j=1,nopes
164  thickness=thickness+thicke(k,indexe+j)*shp2(4,j)
165  enddo
166  tlayer(kk)=tlayer(kk)+thickness
167  xlayer(k,kk)=thickness
168  enddo
169  enddo
170  iflag=3
171 !
172  ilayer=0
173  do k=1,4
174  dlayer(k)=0.d0
175  enddo
176 !
177 !
178 ! S6 - composite element
179 !
180  elseif(lakonl(4:5).eq.'15') then
181  mint2d=3
182  nopes=6
183 ! determining the number of layers
184 !
185  nlayer=0
186  do k=1,mi(3)
187  if(ielmat(k,i).ne.0) then
188  nlayer=nlayer+1
189  endif
190  enddo
191 !
192 ! determining the layer thickness and global thickness
193 ! at the shell integration points
194 !
195  iflag=1
196  indexe=ipkon(i)
197  do kk=1,mint2d
198  xi=gauss3d10(1,kk)
199  et=gauss3d10(2,kk)
200  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
201  tlayer(kk)=0.d0
202  do k=1,nlayer
203  thickness=0.d0
204  do j=1,nopes
205  thickness=thickness+thicke(k,indexe+j)*shp2(4,j)
206  enddo
207  tlayer(kk)=tlayer(kk)+thickness
208  xlayer(k,kk)=thickness
209  enddo
210  enddo
211  iflag=3
212 !
213  ilayer=0
214  do k=1,3
215  dlayer(k)=0.d0
216  enddo
217 !
218  endif
219 !
220  indexe=ipkon(i)
221 c Bernhardi start
222  if(lakonl(1:5).eq.'C3D8I') then
223  nope=11
224  elseif(lakonl(4:5).eq.'20') then
225 c Bernhardi end
226  nope=20
227  elseif(lakonl(4:4).eq.'8') then
228  nope=8
229  elseif(lakonl(4:5).eq.'10') then
230  nope=10
231  elseif(lakonl(4:4).eq.'4') then
232  nope=4
233  elseif(lakonl(4:5).eq.'15') then
234  nope=15
235  elseif(lakonl(4:4).eq.'6') then
236  nope=6
237  elseif((lakonl(1:1).eq.'E').and.(lakonl(7:7).ne.'F')) then
238 !
239 ! spring elements, contact spring elements and
240 ! dashpot elements
241 !
242  if(lakonl(7:7).eq.'C') then
243 !
244 ! contact spring elements
245 !
246  if(mortar.eq.1) then
247 !
248 ! face-to-face penalty
249 !
250  nope=kon(ipkon(i))
251  elseif(mortar.eq.0) then
252 !
253 ! node-to-face penalty
254 !
255  nope=ichar(lakonl(8:8))-47
256  konl(nope+1)=kon(indexe+nope+1)
257  endif
258  else
259 !
260 ! genuine spring elements and dashpot elements
261 !
262  nope=ichar(lakonl(8:8))-47
263  endif
264  else
265  cycle
266  endif
267 !
268  if(lakonl(4:5).eq.'8R') then
269  mint3d=1
270  elseif(lakonl(4:7).eq.'20RB') then
271  if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
272  mint3d=50
273  else
274  call beamintscheme(lakonl,mint3d,ielprop(i),prop,
275  & null,xi,et,ze,weight)
276  endif
277  elseif((lakonl(4:4).eq.'8').or.
278  & (lakonl(4:6).eq.'20R')) then
279  if(lakonl(7:8).eq.'LC') then
280  mint3d=8*nlayer
281  else
282  mint3d=8
283  endif
284  elseif(lakonl(4:4).eq.'2') then
285  mint3d=27
286  elseif(lakonl(4:5).eq.'10') then
287  mint3d=4
288  elseif(lakonl(4:4).eq.'4') then
289  mint3d=1
290  elseif(lakonl(4:5).eq.'15') then
291  if(lakonl(7:8).eq.'LC') then
292  mint3d=6*nlayer
293  else
294  mint3d=9
295  endif
296  elseif(lakonl(4:4).eq.'6') then
297  mint3d=2
298  elseif(lakonl(1:1).eq.'E') then
299  mint3d=0
300  endif
301 !
302  do j=1,nope
303  konl(j)=kon(indexe+j)
304  do k=1,3
305  xl(k,j)=co(k,konl(j))
306  vl(k,j)=v(k,konl(j))
307  voldl(k,j)=vold(k,konl(j))
308  enddo
309  enddo
310 !
311 ! q contains the nodal forces per element; initialisation of q
312 !
313  if((iperturb(1).ge.2).or.((iperturb(1).le.0).and.(iout.lt.1)))
314  & then
315  do m1=1,nope
316  do m2=0,mi(2)
317  q(m2,m1)=fn(m2,konl(m1))
318  enddo
319  enddo
320  endif
321 !
322 ! calculating the forces for the contact elements
323 !
324 c write(*,*) 'resultsmech ',i,lakonl,mint3d
325  if(mint3d.eq.0) then
326 !
327 ! "normal" spring and dashpot elements
328 !
329  kode=nelcon(1,imat)
330  if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'1').or.
331  & (lakonl(7:7).eq.'2')) then
332  t0l=0.d0
333  t1l=0.d0
334  if(ithermal(1).eq.1) then
335  t0l=(t0(konl(1))+t0(konl(2)))/2.d0
336  t1l=(t1(konl(1))+t1(konl(2)))/2.d0
337  elseif(ithermal(1).ge.2) then
338  t0l=(t0(konl(1))+t0(konl(2)))/2.d0
339  t1l=(vold(0,konl(1))+vold(0,konl(2)))/2.d0
340  endif
341  endif
342 !
343 ! spring elements (including contact springs)
344 !
345  if(lakonl(2:2).eq.'S') then
346  if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'1').or.
347  & (lakonl(7:7).eq.'2').or.((mortar.eq.0).and.
348  & ((nmethod.ne.1).or.(iperturb(1).ge.2).or.(iout.ne.-1))))
349  & then
350  call springforc_n2f(xl,konl,vl,imat,elcon,nelcon,elas,
351  & fnl,ncmat_,ntmat_,nope,lakonl,t1l,kode,elconloc,
352  & plicon,nplicon,npmat_,ener(1,i),nener,
353  & stx(1,1,i),mi,springarea(1,konl(nope+1)),nmethod,
354  & ne0,nstate_,xstateini,xstate,reltime,
355  & ielas,ener(1,i+ne),ielorien,orab,norien,i)
356  elseif((mortar.eq.1).and.
357  & ((nmethod.ne.1).or.(iperturb(1).ge.2).or.(iout.ne.-1)))
358  & then
359  iloc=kon(indexe+nope+1)
360  jfaces=kon(indexe+nope+2)
361  igauss=kon(indexe+nope+1)
362  call springforc_f2f(xl,vl,imat,elcon,nelcon,elas,
363  & fnl,ncmat_,ntmat_,nope,lakonl,t1l,kode,elconloc,
364  & plicon,nplicon,npmat_,ener(1,i),nener,
365  & stx(1,1,i),mi,springarea(1,iloc),nmethod,
366  & ne0,nstate_,xstateini,xstate,reltime,
367  & ielas,iloc,jfaces,igauss,pslavsurf,pmastsurf,
368  & clearini,ener(1,i+ne),kscale,konl,iout,i)
369  endif
370 !
371 ! next lines are not executed in linstatic.c before the
372 ! setup of the stiffness matrix (i.e. nmethod=1 and
373 ! iperturb(1)<1 and iout=-1).
374 !
375  if((lakonl(7:7).eq.'A').or.
376  & ((nmethod.ne.1).or.(iperturb(1).ge.2).or.(iout.ne.-1)))
377  & then
378  do j=1,nope
379  do k=1,3
380  fn(k,konl(j))=fn(k,konl(j))+fnl(k,j)
381  enddo
382  enddo
383  endif
384  endif
385  elseif(ikin.eq.1) then
386  do j=1,nope
387  do k=1,3
388  veoldl(k,j)=veold(k,konl(j))
389  enddo
390  enddo
391  endif
392 !
393  do jj=1,mint3d
394  if(lakonl(4:5).eq.'8R') then
395  xi=gauss3d1(1,jj)
396  et=gauss3d1(2,jj)
397  ze=gauss3d1(3,jj)
398  weight=weight3d1(jj)
399  elseif(lakonl(4:7).eq.'20RB') then
400  if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
401  xi=gauss3d13(1,jj)
402  et=gauss3d13(2,jj)
403  ze=gauss3d13(3,jj)
404  weight=weight3d13(jj)
405  else
406  call beamintscheme(lakonl,mint3d,ielprop(i),prop,
407  & jj,xi,et,ze,weight)
408  endif
409  elseif((lakonl(4:4).eq.'8').or.
410  & (lakonl(4:6).eq.'20R'))
411  & then
412  if(lakonl(7:8).ne.'LC') then
413  xi=gauss3d2(1,jj)
414  et=gauss3d2(2,jj)
415  ze=gauss3d2(3,jj)
416  weight=weight3d2(jj)
417  else
418  kl=mod(jj,8)
419  if(kl.eq.0) kl=8
420 !
421  xi=gauss3d2(1,kl)
422  et=gauss3d2(2,kl)
423  ze=gauss3d2(3,kl)
424  weight=weight3d2(kl)
425 !
426  ki=mod(jj,4)
427  if(ki.eq.0) ki=4
428 !
429  if(kl.eq.1) then
430  ilayer=ilayer+1
431  if(ilayer.gt.1) then
432  do k=1,4
433  dlayer(k)=dlayer(k)+xlayer(ilayer-1,k)
434  enddo
435  endif
436  endif
437  ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/
438  & tlayer(ki)-1.d0
439  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
440 !
441 ! material and orientation
442 !
443  imat=ielmat(ilayer,i)
444  amat=matname(imat)
445  if(norien.gt.0) then
446  iorien=ielorien(ilayer,i)
447  else
448  iorien=0
449  endif
450 !
451  if(nelcon(1,imat).lt.0) then
452  ihyper=1
453  else
454  ihyper=0
455  endif
456  endif
457  elseif(lakonl(4:4).eq.'2') then
458  xi=gauss3d3(1,jj)
459  et=gauss3d3(2,jj)
460  ze=gauss3d3(3,jj)
461  weight=weight3d3(jj)
462  elseif(lakonl(4:5).eq.'10') then
463  xi=gauss3d5(1,jj)
464  et=gauss3d5(2,jj)
465  ze=gauss3d5(3,jj)
466  weight=weight3d5(jj)
467  elseif(lakonl(4:4).eq.'4') then
468  xi=gauss3d4(1,jj)
469  et=gauss3d4(2,jj)
470  ze=gauss3d4(3,jj)
471  weight=weight3d4(jj)
472  elseif(lakonl(4:5).eq.'15') then
473  if(lakonl(7:8).ne.'LC') then
474  xi=gauss3d8(1,jj)
475  et=gauss3d8(2,jj)
476  ze=gauss3d8(3,jj)
477  weight=weight3d8(jj)
478  else
479  kl=mod(jj,6)
480  if(kl.eq.0) kl=6
481 !
482  xi=gauss3d10(1,kl)
483  et=gauss3d10(2,kl)
484  ze=gauss3d10(3,kl)
485  weight=weight3d10(kl)
486 !
487  ki=mod(jj,3)
488  if(ki.eq.0) ki=3
489 !
490  if(kl.eq.1) then
491  ilayer=ilayer+1
492  if(ilayer.gt.1) then
493  do k=1,3
494  dlayer(k)=dlayer(k)+xlayer(ilayer-1,k)
495  enddo
496  endif
497  endif
498  ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/
499  & tlayer(ki)-1.d0
500  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
501 !
502 ! material and orientation
503 !
504  imat=ielmat(ilayer,i)
505  amat=matname(imat)
506  if(norien.gt.0) then
507  iorien=ielorien(ilayer,i)
508  else
509  iorien=0
510  endif
511 !
512  if(nelcon(1,imat).lt.0) then
513  ihyper=1
514  else
515  ihyper=0
516  endif
517  endif
518  elseif(lakonl(4:4).eq.'6') then
519  xi=gauss3d7(1,jj)
520  et=gauss3d7(2,jj)
521  ze=gauss3d7(3,jj)
522  weight=weight3d7(jj)
523  endif
524 !
525 c Bernhardi start
526  if(lakonl(1:5).eq.'C3D8R') then
527  call shape8hr(xl,xsj,shp,gs,a)
528  elseif(lakonl(1:5).eq.'C3D8I') then
529  call shape8hu(xi,et,ze,xl,xsj,shp,iflag)
530  elseif(nope.eq.20) then
531 c Bernhardi end
532  if(lakonl(7:7).eq.'A') then
533  call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag)
534  elseif((lakonl(7:7).eq.'E').or.
535  & (lakonl(7:7).eq.'S')) then
536  call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag)
537  else
538  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
539  endif
540  elseif(nope.eq.8) then
541  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
542  elseif(nope.eq.10) then
543  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
544  elseif(nope.eq.4) then
545  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
546  elseif(nope.eq.15) then
547  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
548  else
549  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
550  endif
551 !
552 ! vkl(m2,m3) contains the derivative of the m2-
553 ! component of the displacement with respect to
554 ! direction m3
555 !
556  do m2=1,3
557  do m3=1,3
558  vkl(m2,m3)=0.d0
559  enddo
560  enddo
561 !
562  do m1=1,nope
563  do m2=1,3
564  do m3=1,3
565  vkl(m2,m3)=vkl(m2,m3)+shp(m3,m1)*vl(m2,m1)
566  enddo
567 c write(*,*) 'vnoeie',i,konl(m1),(vkl(m2,k),k=1,3)
568  enddo
569  enddo
570 !
571 ! for frequency analysis or buckling with preload the
572 ! strains are calculated with respect to the deformed
573 ! configuration
574 ! for a linear iteration within a nonlinear increment:
575 ! the tangent matrix is calculated at strain at the end
576 ! of the previous increment
577 !
578  if((iperturb(1).eq.1).or.(iperturb(1).eq.-1))then
579  do m2=1,3
580  do m3=1,3
581  vokl(m2,m3)=0.d0
582  enddo
583  enddo
584 !
585  do m1=1,nope
586  do m2=1,3
587  do m3=1,3
588  vokl(m2,m3)=vokl(m2,m3)+
589  & shp(m3,m1)*voldl(m2,m1)
590  enddo
591  enddo
592  enddo
593  endif
594 !
595  kode=nelcon(1,imat)
596 !
597 ! calculating the strain
598 !
599 ! attention! exy,exz and eyz are engineering strains!
600 !
601  exx=vkl(1,1)
602  eyy=vkl(2,2)
603  ezz=vkl(3,3)
604  exy=vkl(1,2)+vkl(2,1)
605  exz=vkl(1,3)+vkl(3,1)
606  eyz=vkl(2,3)+vkl(3,2)
607 !
608  if(iperturb(2).eq.1) then
609 !
610 ! Lagrangian strain
611 !
612  exx=exx+(vkl(1,1)**2+vkl(2,1)**2+vkl(3,1)**2)/2.d0
613  eyy=eyy+(vkl(1,2)**2+vkl(2,2)**2+vkl(3,2)**2)/2.d0
614  ezz=ezz+(vkl(1,3)**2+vkl(2,3)**2+vkl(3,3)**2)/2.d0
615  exy=exy+vkl(1,1)*vkl(1,2)+vkl(2,1)*vkl(2,2)+
616  & vkl(3,1)*vkl(3,2)
617  exz=exz+vkl(1,1)*vkl(1,3)+vkl(2,1)*vkl(2,3)+
618  & vkl(3,1)*vkl(3,3)
619  eyz=eyz+vkl(1,2)*vkl(1,3)+vkl(2,2)*vkl(2,3)+
620  & vkl(3,2)*vkl(3,3)
621 !
622 ! for frequency analysis or buckling with preload the
623 ! strains are calculated with respect to the deformed
624 ! configuration
625 !
626  elseif(iperturb(1).eq.1) then
627  exx=exx+vokl(1,1)*vkl(1,1)+vokl(2,1)*vkl(2,1)+
628  & vokl(3,1)*vkl(3,1)
629  eyy=eyy+vokl(1,2)*vkl(1,2)+vokl(2,2)*vkl(2,2)+
630  & vokl(3,2)*vkl(3,2)
631  ezz=ezz+vokl(1,3)*vkl(1,3)+vokl(2,3)*vkl(2,3)+
632  & vokl(3,3)*vkl(3,3)
633  exy=exy+vokl(1,1)*vkl(1,2)+vokl(1,2)*vkl(1,1)+
634  & vokl(2,1)*vkl(2,2)+vokl(2,2)*vkl(2,1)+
635  & vokl(3,1)*vkl(3,2)+vokl(3,2)*vkl(3,1)
636  exz=exz+vokl(1,1)*vkl(1,3)+vokl(1,3)*vkl(1,1)+
637  & vokl(2,1)*vkl(2,3)+vokl(2,3)*vkl(2,1)+
638  & vokl(3,1)*vkl(3,3)+vokl(3,3)*vkl(3,1)
639  eyz=eyz+vokl(1,2)*vkl(1,3)+vokl(1,3)*vkl(1,2)+
640  & vokl(2,2)*vkl(2,3)+vokl(2,3)*vkl(2,2)+
641  & vokl(3,2)*vkl(3,3)+vokl(3,3)*vkl(3,2)
642  endif
643 !
644 ! storing the local strains
645 !
646  if(iperturb(1).ne.-1) then
647  eloc(1)=exx
648  eloc(2)=eyy
649  eloc(3)=ezz
650  eloc(4)=exy/2.d0
651  eloc(5)=exz/2.d0
652  eloc(6)=eyz/2.d0
653  else
654 !
655 ! linear iteration within a nonlinear increment:
656 !
657  eloc(1)=vokl(1,1)+
658  & (vokl(1,1)**2+vokl(2,1)**2+vokl(3,1)**2)/2.d0
659  eloc(2)=vokl(2,2)+
660  & (vokl(1,2)**2+vokl(2,2)**2+vokl(3,2)**2)/2.d0
661  eloc(3)=vokl(3,3)+
662  & (vokl(1,3)**2+vokl(2,3)**2+vokl(3,3)**2)/2.d0
663  eloc(4)=(vokl(1,2)+vokl(2,1)+vokl(1,1)*vokl(1,2)+
664  & vokl(2,1)*vokl(2,2)+vokl(3,1)*vokl(3,2))/2.d0
665  eloc(5)=(vokl(1,3)+vokl(3,1)+vokl(1,1)*vokl(1,3)+
666  & vokl(2,1)*vokl(2,3)+vokl(3,1)*vokl(3,3))/2.d0
667  eloc(6)=(vokl(2,3)+vokl(3,2)+vokl(1,2)*vokl(1,3)+
668  & vokl(2,2)*vokl(2,3)+vokl(3,2)*vokl(3,3))/2.d0
669  endif
670 !
671 ! calculating the deformation gradient (needed to
672 ! convert the element stiffness matrix from spatial
673 ! coordinates to material coordinates
674 ! deformation plasticity)
675 !
676  if((kode.eq.-50).or.(kode.le.-100)) then
677 !
678 ! calculating the deformation gradient
679 !
680 c Bernhardi start
681  xkl(1,1)=vkl(1,1)+1.0d0
682  xkl(2,2)=vkl(2,2)+1.0d0
683  xkl(3,3)=vkl(3,3)+1.0d0
684 c Bernhardi end
685  xkl(1,2)=vkl(1,2)
686  xkl(1,3)=vkl(1,3)
687  xkl(2,3)=vkl(2,3)
688  xkl(2,1)=vkl(2,1)
689  xkl(3,1)=vkl(3,1)
690  xkl(3,2)=vkl(3,2)
691 !
692 ! calculating the Jacobian
693 !
694  vj=xkl(1,1)*(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2))
695  & -xkl(1,2)*(xkl(2,1)*xkl(3,3)-xkl(2,3)*xkl(3,1))
696  & +xkl(1,3)*(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1))
697 !
698 ! inversion of the deformation gradient (only for
699 ! deformation plasticity)
700 !
701  if(kode.eq.-50) then
702 !
703  ckl(1,1)=(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2))/vj
704  ckl(2,2)=(xkl(1,1)*xkl(3,3)-xkl(1,3)*xkl(3,1))/vj
705  ckl(3,3)=(xkl(1,1)*xkl(2,2)-xkl(1,2)*xkl(2,1))/vj
706  ckl(1,2)=(xkl(1,3)*xkl(3,2)-xkl(1,2)*xkl(3,3))/vj
707  ckl(1,3)=(xkl(1,2)*xkl(2,3)-xkl(2,2)*xkl(1,3))/vj
708  ckl(2,3)=(xkl(2,1)*xkl(1,3)-xkl(1,1)*xkl(2,3))/vj
709  ckl(2,1)=(xkl(3,1)*xkl(2,3)-xkl(2,1)*xkl(3,3))/vj
710  ckl(3,1)=(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1))/vj
711  ckl(3,2)=(xkl(3,1)*xkl(1,2)-xkl(1,1)*xkl(3,2))/vj
712 !
713 ! converting the Lagrangian strain into Eulerian
714 ! strain (only for deformation plasticity)
715 !
716  cauchy=0
717  call str2mat(eloc,ckl,vj,cauchy)
718  endif
719 !
720  endif
721 !
722 ! calculating fields for incremental plasticity
723 !
724  if(kode.le.-100) then
725 !
726 ! calculating the deformation gradient at the
727 ! start of the increment
728 !
729 ! calculating the displacement gradient at the
730 ! start of the increment
731 !
732  do m2=1,3
733  do m3=1,3
734  vikl(m2,m3)=0.d0
735  enddo
736  enddo
737 !
738  do m1=1,nope
739  do m2=1,3
740  do m3=1,3
741  vikl(m2,m3)=vikl(m2,m3)
742  & +shp(m3,m1)*vini(m2,konl(m1))
743  enddo
744  enddo
745  enddo
746 !
747 ! calculating the deformation gradient of the old
748 ! fields
749 !
750  xikl(1,1)=vikl(1,1)+1
751  xikl(2,2)=vikl(2,2)+1.
752  xikl(3,3)=vikl(3,3)+1.
753  xikl(1,2)=vikl(1,2)
754  xikl(1,3)=vikl(1,3)
755  xikl(2,3)=vikl(2,3)
756  xikl(2,1)=vikl(2,1)
757  xikl(3,1)=vikl(3,1)
758  xikl(3,2)=vikl(3,2)
759 !
760 ! calculating the Jacobian
761 !
762  vij=xikl(1,1)*(xikl(2,2)*xikl(3,3)
763  & -xikl(2,3)*xikl(3,2))
764  & -xikl(1,2)*(xikl(2,1)*xikl(3,3)
765  & -xikl(2,3)*xikl(3,1))
766  & +xikl(1,3)*(xikl(2,1)*xikl(3,2)
767  & -xikl(2,2)*xikl(3,1))
768 !
769 ! stresses at the start of the increment
770 !
771  do m1=1,6
772  stre(m1)=stiini(m1,jj,i)
773  enddo
774 !
775  endif
776 !
777 ! prestress values
778 !
779  if(iprestr.eq.1) then
780  do kk=1,6
781  beta(kk)=-prestr(kk,jj,i)
782  enddo
783  else
784  do kk=1,6
785  beta(kk)=0.d0
786  enddo
787  endif
788 !
789  if(ithermal(1).ge.1) then
790 !
791 ! calculating the temperature difference in
792 ! the integration point
793 !
794  t0l=0.d0
795  t1l=0.d0
796  if(ithermal(1).eq.1) then
797  if((lakonl(4:5).eq.'8 ').or.
798  & (lakonl(4:5).eq.'8I')) then
799  do i1=1,8
800  t0l=t0l+t0(konl(i1))/8.d0
801  t1l=t1l+t1(konl(i1))/8.d0
802  enddo
803  elseif(lakonl(4:6).eq.'20 ') then
804  nopered=20
805  call lintemp(t0,t1,konl,nopered,jj,t0l,t1l)
806  elseif(lakonl(4:6).eq.'10T') then
807  call linscal10(t0,konl,t0l,null,shp)
808  call linscal10(t1,konl,t1l,null,shp)
809  else
810  do i1=1,nope
811  t0l=t0l+shp(4,i1)*t0(konl(i1))
812  t1l=t1l+shp(4,i1)*t1(konl(i1))
813  enddo
814  endif
815  elseif(ithermal(1).ge.2) then
816  if((lakonl(4:5).eq.'8 ').or.
817  & (lakonl(4:5).eq.'8I')) then
818  do i1=1,8
819  t0l=t0l+t0(konl(i1))/8.d0
820  t1l=t1l+vold(0,konl(i1))/8.d0
821  enddo
822  elseif(lakonl(4:6).eq.'20 ') then
823  nopered=20
824  call lintemp_th(t0,vold,konl,nopered,jj,t0l,t1l,mi)
825  elseif(lakonl(4:6).eq.'10T') then
826  call linscal10(t0,konl,t0l,null,shp)
827  call linscal10(vold,konl,t1l,mi(2),shp)
828  else
829  do i1=1,nope
830  t0l=t0l+shp(4,i1)*t0(konl(i1))
831  t1l=t1l+shp(4,i1)*vold(0,konl(i1))
832  enddo
833  endif
834  endif
835  tt=t1l-t0l
836  endif
837 !
838 ! calculating the coordinates of the integration point
839 ! for material orientation purposes (for cylindrical
840 ! coordinate systems)
841 !
842  if((iorien.gt.0).or.(kode.le.-100)) then
843  do j=1,3
844  pgauss(j)=0.d0
845  do i1=1,nope
846  pgauss(j)=pgauss(j)+shp(4,i1)*co(j,konl(i1))
847  enddo
848  enddo
849  endif
850 !
851 ! material data; for linear elastic materials
852 ! this includes the calculation of the stiffness
853 ! matrix
854 !
855  istiff=0
856 !
857  call materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon,
858  & nalcon,imat,amat,iorien,pgauss,orab,ntmat_,
859  & elas,rho,i,ithermal,alzero,mattyp,t0l,t1l,ihyper,
860  & istiff,elconloc,eth,kode,plicon,nplicon,
861  & plkcon,nplkcon,npmat_,plconloc,mi(1),dtime,jj,
862  & xstiff,ncmat_)
863 !
864 ! determining the mechanical strain
865 !
866  if(ithermal(1).ne.0) then
867  do m1=1,6
868  emec(m1)=eloc(m1)-eth(m1)
869 c emec0(m1)=emeini(m1,jj,i)
870  enddo
871  else
872  do m1=1,6
873  emec(m1)=eloc(m1)
874 c emec0(m1)=emeini(m1,jj,i)
875  enddo
876  endif
877  if(kode.le.-100) then
878  do m1=1,6
879  emec0(m1)=emeini(m1,jj,i)
880  enddo
881  endif
882 !
883 ! subtracting the plastic initial strains
884 !
885  if(iprestr.eq.2) then
886  do m1=1,6
887  emec(m1)=emec(m1)-prestr(m1,jj,i)
888  enddo
889  endif
890 !
891 ! calculating the local stiffness and stress
892 !
893  call mechmodel(elconloc,elas,emec,kode,emec0,ithermal,
894  & icmd,beta,stre,xkl,ckl,vj,xikl,vij,
895  & plconloc,xstate,xstateini,ielas,
896  & amat,t1l,dtime,time,ttime,i,jj,nstate_,mi(1),
897  & iorien,pgauss,orab,eloc,mattyp,qa(3),istep,iinc,
898  & ipkon,nmethod,iperturb,qa(4))
899 !
900  if(((nmethod.ne.4).or.(iperturb(1).ne.0)).and.
901  & (nmethod.ne.5).and.(icmd.ne.3)) then
902  do m1=1,21
903  xstiff(m1,jj,i)=elas(m1)
904  enddo
905  endif
906 !
907  if(iperturb(1).eq.-1) then
908 !
909 ! if the forced displacements were changed at
910 ! the start of a nonlinear step, the nodal
911 ! forces due do this displacements are
912 ! calculated in a purely linear way, and
913 ! the first iteration is purely linear in order
914 ! to allow the displacements to redistribute
915 ! in a quasi-static way (only applies to
916 ! quasi-static analyses (*STATIC))
917 !
918  eloc(1)=exx-vokl(1,1)
919  eloc(2)=eyy-vokl(2,2)
920  eloc(3)=ezz-vokl(3,3)
921  eloc(4)=exy-(vokl(1,2)+vokl(2,1))
922  eloc(5)=exz-(vokl(1,3)+vokl(3,1))
923  eloc(6)=eyz-(vokl(2,3)+vokl(3,2))
924 !
925  if(mattyp.eq.1) then
926  e=elas(1)
927  un=elas(2)
928  um=e/(1.d0+un)
929  al=un*um/(1.d0-2.d0*un)
930  um=um/2.d0
931  am1=al*(eloc(1)+eloc(2)+eloc(3))
932  stre(1)=am1+2.d0*um*eloc(1)
933  stre(2)=am1+2.d0*um*eloc(2)
934  stre(3)=am1+2.d0*um*eloc(3)
935  stre(4)=um*eloc(4)
936  stre(5)=um*eloc(5)
937  stre(6)=um*eloc(6)
938  elseif(mattyp.eq.2) then
939  stre(1)=eloc(1)*elas(1)+eloc(2)*elas(2)
940  & +eloc(3)*elas(4)
941  stre(2)=eloc(1)*elas(2)+eloc(2)*elas(3)
942  & +eloc(3)*elas(5)
943  stre(3)=eloc(1)*elas(4)+eloc(2)*elas(5)
944  & +eloc(3)*elas(6)
945  stre(4)=eloc(4)*elas(7)
946  stre(5)=eloc(5)*elas(8)
947  stre(6)=eloc(6)*elas(9)
948  elseif(mattyp.eq.3) then
949  stre(1)=eloc(1)*elas(1)+eloc(2)*elas(2)+
950  & eloc(3)*elas(4)+eloc(4)*elas(7)+
951  & eloc(5)*elas(11)+eloc(6)*elas(16)
952  stre(2)=eloc(1)*elas(2)+eloc(2)*elas(3)+
953  & eloc(3)*elas(5)+eloc(4)*elas(8)+
954  & eloc(5)*elas(12)+eloc(6)*elas(17)
955  stre(3)=eloc(1)*elas(4)+eloc(2)*elas(5)+
956  & eloc(3)*elas(6)+eloc(4)*elas(9)+
957  & eloc(5)*elas(13)+eloc(6)*elas(18)
958  stre(4)=eloc(1)*elas(7)+eloc(2)*elas(8)+
959  & eloc(3)*elas(9)+eloc(4)*elas(10)+
960  & eloc(5)*elas(14)+eloc(6)*elas(19)
961  stre(5)=eloc(1)*elas(11)+eloc(2)*elas(12)+
962  & eloc(3)*elas(13)+eloc(4)*elas(14)+
963  & eloc(5)*elas(15)+eloc(6)*elas(20)
964  stre(6)=eloc(1)*elas(16)+eloc(2)*elas(17)+
965  & eloc(3)*elas(18)+eloc(4)*elas(19)+
966  & eloc(5)*elas(20)+eloc(6)*elas(21)
967  endif
968  endif
969 !
970 ! updating the internal energy and mechanical strain
971 ! for user materials (kode<=-100) the mechanical strain has to
972 ! be updated at the end of each increment (also if no output
973 ! is requested), since it is input to the umat routine
974 !
975  if((iout.gt.0).or.(iout.eq.-2).or.(kode.le.-100).or.
976  & ((nmethod.eq.4).and.(iperturb(1).gt.1).and.
977  & (ithermal(1).le.1))) then
978  if(ithermal(1).eq.0) then
979  do m1=1,6
980  eth(m1)=0.d0
981  enddo
982  endif
983  if(nener.eq.1) then
984  ener(jj,i)=enerini(jj,i)+
985  & ((eloc(1)-eth(1)-emeini(1,jj,i))*
986  & (stre(1)+stiini(1,jj,i))+
987  & (eloc(2)-eth(2)-emeini(2,jj,i))*
988  & (stre(2)+stiini(2,jj,i))+
989  & (eloc(3)-eth(3)-emeini(3,jj,i))*
990  & (stre(3)+stiini(3,jj,i)))/2.d0+
991  & (eloc(4)-eth(4)-emeini(4,jj,i))*(stre(4)+stiini(4,jj,i))+
992  & (eloc(5)-eth(5)-emeini(5,jj,i))*(stre(5)+stiini(5,jj,i))+
993  & (eloc(6)-eth(6)-emeini(6,jj,i))*(stre(6)+stiini(6,jj,i))
994  endif
995  eme(1,jj,i)=eloc(1)-eth(1)
996  eme(2,jj,i)=eloc(2)-eth(2)
997  eme(3,jj,i)=eloc(3)-eth(3)
998  eme(4,jj,i)=eloc(4)-eth(4)
999  eme(5,jj,i)=eloc(5)-eth(5)
1000  eme(6,jj,i)=eloc(6)-eth(6)
1001  endif
1002 !
1003  if((iout.gt.0).or.(iout.eq.-2).or.(kode.le.-100)) then
1004 !
1005  eei(1,jj,i)=eloc(1)
1006  eei(2,jj,i)=eloc(2)
1007  eei(3,jj,i)=eloc(3)
1008  eei(4,jj,i)=eloc(4)
1009  eei(5,jj,i)=eloc(5)
1010  eei(6,jj,i)=eloc(6)
1011  endif
1012 !
1013 ! updating the kinetic energy
1014 !
1015  if(ikin.eq.1) then
1016 !
1017  call materialdata_rho(rhcon,nrhcon,imat,rho,t1l,
1018  & ntmat_,ithermal)
1019  do m1=1,3
1020  vel(m1,1)=0.d0
1021  do i1= 1,nope
1022  vel(m1,1)=vel(m1,1)+shp(4,i1)*veoldl(m1,i1)
1023  enddo
1024  enddo
1025  ener(jj,i+ne)=rho*(vel(1,1)*vel(1,1)+
1026  & vel(2,1)*vel(2,1)+ vel(3,1)*vel(3,1))/2.d0
1027  endif
1028 !
1029  skl(1,1)=stre(1)
1030  skl(2,2)=stre(2)
1031  skl(3,3)=stre(3)
1032  skl(2,1)=stre(4)
1033  skl(3,1)=stre(5)
1034  skl(3,2)=stre(6)
1035 !
1036  stx(1,jj,i)=skl(1,1)
1037  stx(2,jj,i)=skl(2,2)
1038  stx(3,jj,i)=skl(3,3)
1039  stx(4,jj,i)=skl(2,1)
1040  stx(5,jj,i)=skl(3,1)
1041  stx(6,jj,i)=skl(3,2)
1042 !
1043  skl(1,2)=skl(2,1)
1044  skl(1,3)=skl(3,1)
1045  skl(2,3)=skl(3,2)
1046 !
1047 ! calculation of the nodal forces
1048 !
1049  if(calcul_fn.eq.1)then
1050 !
1051 ! calculating fn using skl
1052 !
1053  do m1=1,nope
1054  do m2=1,3
1055 !
1056 ! linear elastic part
1057 !
1058  do m3=1,3
1059  fn(m2,konl(m1))=fn(m2,konl(m1))+
1060  & xsj*skl(m2,m3)*shp(m3,m1)*weight
1061  enddo
1062 !
1063 ! nonlinear geometric part
1064 !
1065  if(iperturb(2).eq.1) then
1066  do m3=1,3
1067  do m4=1,3
1068  fn(m2,konl(m1))=fn(m2,konl(m1))+
1069  & xsj*skl(m4,m3)*weight*
1070  & (vkl(m2,m4)*shp(m3,m1)+
1071  & vkl(m2,m3)*shp(m4,m1))/2.d0
1072  enddo
1073  enddo
1074  endif
1075 !
1076  enddo
1077  enddo
1078 c Bernhardi start
1079  if(lakonl(1:5).eq.'C3D8R') then
1080  call hgforce (fn,elas,a,gs,vl,mi,konl)
1081  endif
1082 c Bernhardi end
1083  endif
1084 !
1085 ! calculation of the Cauchy stresses
1086 !
1087  if(calcul_cauchy.eq.1) then
1088 !
1089 ! changing the displacement gradients into
1090 ! deformation gradients
1091 !
1092 c if(kode.ne.-50) then
1093  if((kode.ne.-50).and.(kode.gt.-100)) then
1094 c Bernhardi start
1095  xkl(1,1)=vkl(1,1)+1.0d0
1096  xkl(2,2)=vkl(2,2)+1.0d0
1097  xkl(3,3)=vkl(3,3)+1.0d0
1098 c Bernhardi end
1099  xkl(1,2)=vkl(1,2)
1100  xkl(1,3)=vkl(1,3)
1101  xkl(2,3)=vkl(2,3)
1102  xkl(2,1)=vkl(2,1)
1103  xkl(3,1)=vkl(3,1)
1104  xkl(3,2)=vkl(3,2)
1105 !
1106  vj=xkl(1,1)*(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2))
1107  & -xkl(1,2)*(xkl(2,1)*xkl(3,3)-xkl(2,3)*xkl(3,1))
1108  & +xkl(1,3)*(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1))
1109  endif
1110 !
1111  do m1=1,3
1112  do m2=1,m1
1113  ckl(m1,m2)=0.d0
1114  do m3=1,3
1115  do m4=1,3
1116  ckl(m1,m2)=ckl(m1,m2)+
1117  & skl(m3,m4)*xkl(m1,m3)*xkl(m2,m4)
1118  enddo
1119  enddo
1120  ckl(m1,m2)=ckl(m1,m2)/vj
1121  enddo
1122  enddo
1123 !
1124  stx(1,jj,i)=ckl(1,1)
1125  stx(2,jj,i)=ckl(2,2)
1126  stx(3,jj,i)=ckl(3,3)
1127  stx(4,jj,i)=ckl(2,1)
1128  stx(5,jj,i)=ckl(3,1)
1129  stx(6,jj,i)=ckl(3,2)
1130  endif
1131 !
1132  enddo
1133 !
1134 ! q contains the contributions to the nodal force in the nodes
1135 ! belonging to the element at stake from other elements (elements
1136 ! already treated). These contributions have to be
1137 ! subtracted to get the contributions attributable to the element
1138 ! at stake only
1139 !
1140  if(calcul_qa.eq.1) then
1141  do m1=1,nope
1142  do m2=1,3
1143  qa(1)=qa(1)+dabs(fn(m2,konl(m1))-q(m2,m1))
1144  enddo
1145  enddo
1146  nal=nal+3*nope
1147  endif
1148  enddo
1149 !
1150 c write(*,*) 'resultsmech ',qa(4)
1151 !
1152  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
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
springforc_n2f
subroutine springforc_n2f(xl, konl, vl, imat, elcon, nelcon, elas, fnl, ncmat_, ntmat_, nope, lakonl, t1l, kode, elconloc, plicon, nplicon, npmat_, senergy, nener, cstr, mi, springarea, nmethod, ne0, nstate_, xstateini, xstate, reltime, ielas, venergy, ielorien, orab, norien, nelem)
Definition: springforc_n2f.f:24
mechmodel
subroutine mechmodel(elconloc, elas, emec, kode, emec0, ithermal, icmd, beta, stre, xkl, ckl, vj, xikl, vij, plconloc, xstate, xstateini, ielas, amat, t1l, dtime, time, ttime, iel, iint, nstate_, mi, iorien, pgauss, orab, eloc, mattyp, pnewdt, istep, iinc, ipkon, nmethod, iperturb, depvisc)
Definition: mechmodel.f:24
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
shape15w
subroutine shape15w(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape15w.f:20
hgforce
subroutine hgforce(fn, elas, a, gs, vl, mi, konl)
Definition: hgforce.f:20
shape20h_ax
subroutine shape20h_ax(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape20h_ax.f:20
lintemp_th
subroutine lintemp_th(t0, vold, konl, nope, jj, t0l, t1l, mi)
Definition: lintemp_th.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
str2mat
subroutine str2mat(str, ckl, vj, cauchy)
Definition: str2mat.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
cauchy
subroutine cauchy(n, x, l, u, nbd, g, iorder, iwhere, t, d, xcp, m, wy, ws, sy, wt, theta, col, head, p, c, wbp, v, nseg, iprint, sbgnrm, info, epsmch)
Definition: lbfgsb.f:1225
shape8hu
subroutine shape8hu(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape8hu.f:20
nal
static ITG * nal
Definition: results.c:31
shape6tri
subroutine shape6tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape6tri.f:20
resultsmech_u
subroutine resultsmech_u(co, kon, ipkon, lakon, ne, v, stx, elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, ielorien, norien, orab, ntmat_, t0, t1, ithermal, prestr, iprestr, eme, iperturb, fn, iout, qa, vold, nmethod, veold, dtime, time, ttime, plicon, nplicon, plkcon, nplkcon, xstateini, xstiff, xstate, npmat_, matname, mi, ielas, icmd, ncmat_, nstate_, stiini, vini, ener, eei, enerini, istep, iinc, reltime, calcul_fn, calcul_qa, calcul_cauchy, nener, ikin, nal, ne0, thicke, emeini, nelem, ielprop, prop)
Definition: resultsmech_u.f:28
springforc_f2f
subroutine springforc_f2f(xl, vl, imat, elcon, nelcon, elas, fnl, ncmat_, ntmat_, nope, lakonl, t1l, kode, elconloc, plicon, nplicon, npmat_, senergy, nener, cstr, mi, springarea, nmethod, ne0, nstate_, xstateini, xstate, reltime, ielas, iloc, jfaces, igauss, pslavsurf, pmastsurf, clearini, venergy, kscale, konl, iout, nelem)
Definition: springforc_f2f.f:26
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
materialdata_rho
subroutine materialdata_rho(rhcon, nrhcon, imat, rho, t1l, ntmat_, ithermal)
Definition: materialdata_rho.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)