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

Go to the source code of this file.

Functions/Subroutines

subroutine objective_shapeener_dx (co, kon, ipkon, lakon, ne, stx, elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, ielorien, norien, orab, ntmat_, t0, t1, ithermal, prestr, iprestr, iperturb, iout, vold, nmethod, veold, dtime, time, ttime, plicon, nplicon, plkcon, nplkcon, xstateini, xstiff, xstate, npmat_, matname, mi, ielas, icmd, ncmat_, nstate_, stiini, vini, ener, enerini, istep, iinc, springarea, reltime, calcul_qa, nener, ikin, ne0, thicke, emeini, pslavsurf, pmastsurf, mortar, clearini, nea, neb, ielprop, prop, distmin, ndesi, nodedesi, nobject, g0, dgdx, iobject, sti, xener, istartdesi, ialdesi, xdesi, idesvar)
 

Function/Subroutine Documentation

◆ objective_shapeener_dx()

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