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

Go to the source code of this file.

Functions/Subroutines

subroutine objective_mass_dx (co, kon, ipkon, lakon, nelcon, rhcon, ielmat, ielorien, norien, ntmat_, matname, mi, thicke, mortar, nea, neb, ielprop, prop, distmin, ndesi, nodedesi, nobject, g0, dgdx, iobject, xmass, istartdesi, ialdesi, xdesi, idesvar)
 

Function/Subroutine Documentation

◆ objective_mass_dx()

subroutine objective_mass_dx ( real*8, dimension(3,*)  co,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
character*8, dimension(*)  lakon,
integer, dimension(2,*)  nelcon,
real*8, dimension(0:1,ntmat_,*)  rhcon,
integer, dimension(mi(3),*)  ielmat,
integer, dimension(mi(3),*)  ielorien,
integer  norien,
integer  ntmat_,
character*80, dimension(*)  matname,
integer, dimension(*)  mi,
real*8, dimension(mi(3),*)  thicke,
integer  mortar,
integer  nea,
integer  neb,
integer, dimension(*)  ielprop,
real*8, dimension(*)  prop,
real*8  distmin,
integer  ndesi,
integer, dimension(*)  nodedesi,
integer  nobject,
real*8, dimension(*)  g0,
real*8, dimension(ndesi,nobject)  dgdx,
integer  iobject,
real*8, dimension(*)  xmass,
integer, dimension(*)  istartdesi,
integer, dimension(*)  ialdesi,
real*8, dimension(3,*)  xdesi,
integer  idesvar 
)
24 !
25 ! calculates the mass and its derivative w.r.t. the coordinates
26 ! of the mesh
27 !
28  implicit none
29 !
30  character*8 lakon(*),lakonl
31  character*80 amat,matname(*)
32 !
33  integer kon(*),konl(26),nea,neb,mi(*),mint2d,nopes,nelcon(2,*),
34  & ielmat(mi(3),*),ielorien(mi(3),*),ntmat_,ipkon(*),iflag,null,
35  & mt,i,ii,j,k,jj,kk,indexe,nope,norien,ihyper,iactive,
36  & imat,mint3d,iorien,ilayer,nlayer,ki,kl,istartdesi(*),
37  & ielprop(*),mortar,idesvar,node,ndesi,nodedesi(*),
38  & nobject,iobject,ialdesi(*),ij
39 !
40  real*8 co(3,*),shp(4,26),xl(3,26),vl(0:mi(2),26),
41  & prop(*),rhcon(0:1,ntmat_,*),xs2(3,7),thickness,rho,xl2(3,8),
42  & xsj2(3),shp2(7,8),xi,et,ze,
43  & xsj,weight,gs(8,4),a,tlayer(4),dlayer(4),xlayer(mi(3),4),
44  & thicke(mi(3),*),distmin,xmassel,g0(*),xmass(*),xdesi(3,*),
45  & dgdx(ndesi,nobject)
46 !
47  include "gauss.f"
48 !
49  iflag=3
50  null=0
51 !
52  mt=mi(2)+1
53 !
54 ! -------------------------------------------------------------
55 ! Initialisation of the loop for the variation of
56 ! the internal forces
57 ! -------------------------------------------------------------
58 !
59 !
60 ! Loop over all elements in thread
61  do ij=nea,neb
62  if(idesvar.gt.0) then
63  i=ialdesi(ij)
64  else
65  i=ij
66  endif
67 
68 ! initialisation of mass
69 !
70  xmassel=0.d0
71  lakonl=lakon(i)
72 !
73  if(ipkon(i).lt.0) cycle
74 !
75 ! no 3D-fluid elements
76 !
77  if(lakonl(1:1).eq.'F') cycle
78  if(lakonl(1:7).eq.'DCOUP3D') cycle
79 !
80  if(lakonl(7:8).ne.'LC') then
81 !
82  imat=ielmat(1,i)
83  amat=matname(imat)
84  if(norien.gt.0) then
85  iorien=ielorien(1,i)
86  else
87  iorien=0
88  endif
89 !
90  if(nelcon(1,imat).lt.0) then
91  ihyper=1
92  else
93  ihyper=0
94  endif
95  else
96 !
97 ! composite materials
98 !
99 ! determining the number of layers
100 !
101  nlayer=0
102  do k=1,mi(3)
103  if(ielmat(k,i).ne.0) then
104  nlayer=nlayer+1
105  endif
106  enddo
107 !
108 ! determining the layer thickness and global thickness
109 ! at the shell integration points
110 !
111  if(lakonl(4:5).eq.'20') then
112 !
113  mint2d=4
114  nopes=8
115 !
116  iflag=1
117  indexe=ipkon(i)
118  do kk=1,mint2d
119  xi=gauss3d2(1,kk)
120  et=gauss3d2(2,kk)
121  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
122  tlayer(kk)=0.d0
123  do k=1,nlayer
124  thickness=0.d0
125  do j=1,nopes
126  thickness=thickness+thicke(k,indexe+j)*
127  & shp2(4,j)
128  enddo
129  tlayer(kk)=tlayer(kk)+thickness
130  xlayer(k,kk)=thickness
131  enddo
132  enddo
133  iflag=3
134 !
135  ilayer=0
136  do k=1,4
137  dlayer(k)=0.d0
138  enddo
139  elseif(lakonl(4:5).eq.'15') then
140 !
141  mint2d=3
142  nopes=6
143 !
144  iflag=1
145  indexe=ipkon(i)
146  do kk=1,mint2d
147  xi=gauss3d10(1,kk)
148  et=gauss3d10(2,kk)
149  call shape6tri(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,3
165  dlayer(k)=0.d0
166  enddo
167  endif
168 !
169  endif
170 !
171  indexe=ipkon(i)
172 c Bernhardi start
173  if(lakonl(1:5).eq.'C3D8I') then
174  nope=11
175  elseif(lakonl(4:5).eq.'20') then
176 c Bernhardi end
177  nope=20
178  elseif(lakonl(4:4).eq.'2') then
179  nope=26
180  elseif(lakonl(4:4).eq.'8') then
181  nope=8
182  elseif(lakonl(4:5).eq.'10') then
183  nope=10
184  elseif(lakonl(4:5).eq.'14') then
185  nope=14
186  elseif(lakonl(4:4).eq.'4') then
187  nope=4
188  elseif(lakonl(4:5).eq.'15') then
189  nope=15
190  elseif(lakonl(4:4).eq.'6') then
191  nope=6
192  elseif((lakonl(1:1).eq.'E').and.(lakonl(7:7).ne.'F')) then
193 !
194 ! spring elements, contact spring elements and
195 ! dashpot elements
196 !
197  if(lakonl(7:7).eq.'C') then
198 !
199 ! contact spring elements
200 !
201  if(mortar.eq.1) then
202 !
203 ! face-to-face penalty
204 !
205  nope=kon(ipkon(i))
206  elseif(mortar.eq.0) then
207 !
208 ! node-to-face penalty
209 !
210  nope=ichar(lakonl(8:8))-47
211  konl(nope+1)=kon(indexe+nope+1)
212  endif
213  else
214 !
215 ! genuine spring elements and dashpot elements
216 !
217  nope=ichar(lakonl(8:8))-47
218  endif
219  else
220  cycle
221  endif
222 !
223 ! check whether a design variable belongs to the element
224 !
225  if(idesvar.gt.0) then
226  do ii=1,nope
227  node=kon(indexe+ii)
228  if(node.eq.nodedesi(idesvar)) then
229  iactive=ii
230  exit
231  endif
232  enddo
233  endif
234 !
235  if(lakonl(4:5).eq.'8R') then
236  mint3d=1
237  elseif(lakonl(4:7).eq.'20RB') then
238  if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
239  mint3d=50
240  else
241  call beamintscheme(lakonl,mint3d,ielprop(i),prop,
242  & null,xi,et,ze,weight)
243  endif
244  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'26R').or.
245  & (lakonl(4:6).eq.'20R')) then
246  if(lakonl(7:8).eq.'LC') then
247  mint3d=8*nlayer
248  else
249  mint3d=8
250  endif
251  elseif(lakonl(4:4).eq.'2') then
252  mint3d=27
253  elseif((lakonl(4:5).eq.'10').or.(lakonl(4:5).eq.'14')) then
254  mint3d=4
255  elseif(lakonl(4:4).eq.'4') then
256  mint3d=1
257  elseif(lakonl(4:5).eq.'15') then
258  if(lakonl(7:8).eq.'LC') then
259  mint3d=6*nlayer
260  else
261  mint3d=9
262  endif
263  elseif(lakonl(4:4).eq.'6') then
264  mint3d=2
265  elseif(lakonl(1:1).eq.'E') then
266  mint3d=0
267  endif
268 
269  do j=1,nope
270  konl(j)=kon(indexe+j)
271  do k=1,3
272  xl(k,j)=co(k,konl(j))
273  enddo
274  enddo
275 !
276 ! computation of the objective function and the derivate
277 ! if the designnode belongs to the considered element
278 ! if not, next element in loop will be taken
279 !
280  if(idesvar.gt.0) then
281  do j=1,3
282  xl(j,iactive)=xl(j,iactive)+xdesi(j,idesvar)
283  enddo
284  endif
285 !
286  do jj=1,mint3d
287  if(lakonl(4:5).eq.'8R') then
288  xi=gauss3d1(1,jj)
289  et=gauss3d1(2,jj)
290  ze=gauss3d1(3,jj)
291  weight=weight3d1(jj)
292  elseif(lakonl(4:7).eq.'20RB') then
293  if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
294  xi=gauss3d13(1,jj)
295  et=gauss3d13(2,jj)
296  ze=gauss3d13(3,jj)
297  weight=weight3d13(jj)
298  else
299  call beamintscheme(lakonl,mint3d,ielprop(i),prop,
300  & jj,xi,et,ze,weight)
301  endif
302  elseif((lakonl(4:4).eq.'8').or.
303  & (lakonl(4:6).eq.'20R').or.(lakonl(4:6).eq.'26R'))
304  & then
305  if(lakonl(7:8).ne.'LC') then
306  xi=gauss3d2(1,jj)
307  et=gauss3d2(2,jj)
308  ze=gauss3d2(3,jj)
309  weight=weight3d2(jj)
310  else
311  kl=mod(jj,8)
312  if(kl.eq.0) kl=8
313 !
314  xi=gauss3d2(1,kl)
315  et=gauss3d2(2,kl)
316  ze=gauss3d2(3,kl)
317  weight=weight3d2(kl)
318 !
319  ki=mod(jj,4)
320  if(ki.eq.0) ki=4
321 !
322  if(kl.eq.1) then
323  ilayer=ilayer+1
324  if(ilayer.gt.1) then
325  do k=1,4
326  dlayer(k)=dlayer(k)+xlayer(ilayer-1,k)
327  enddo
328  endif
329  endif
330  ze=2.d0*(dlayer(ki)+(ze+1.d0)/
331  & 2.d0*xlayer(ilayer,ki))/tlayer(ki)-1.d0
332  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
333 !
334 ! material and orientation
335 !
336  imat=ielmat(ilayer,i)
337  amat=matname(imat)
338  if(norien.gt.0) then
339  iorien=ielorien(ilayer,i)
340  else
341  iorien=0
342  endif
343 !
344  if(nelcon(1,imat).lt.0) then
345  ihyper=1
346  else
347  ihyper=0
348  endif
349  endif
350  elseif(lakonl(4:4).eq.'2') then
351  xi=gauss3d3(1,jj)
352  et=gauss3d3(2,jj)
353  ze=gauss3d3(3,jj)
354  weight=weight3d3(jj)
355  elseif((lakonl(4:5).eq.'10').or.(lakonl(4:5).eq.'14'))
356  & then
357  xi=gauss3d5(1,jj)
358  et=gauss3d5(2,jj)
359  ze=gauss3d5(3,jj)
360  weight=weight3d5(jj)
361  elseif(lakonl(4:4).eq.'4') then
362  xi=gauss3d4(1,jj)
363  et=gauss3d4(2,jj)
364  ze=gauss3d4(3,jj)
365  weight=weight3d4(jj)
366  elseif(lakonl(4:5).eq.'15') then
367  if(lakonl(7:8).ne.'LC') then
368  xi=gauss3d8(1,jj)
369  et=gauss3d8(2,jj)
370  ze=gauss3d8(3,jj)
371  weight=weight3d8(jj)
372  else
373  kl=mod(jj,6)
374  if(kl.eq.0) kl=6
375 !
376  xi=gauss3d10(1,kl)
377  et=gauss3d10(2,kl)
378  ze=gauss3d10(3,kl)
379  weight=weight3d10(kl)
380 !
381  ki=mod(jj,3)
382  if(ki.eq.0) ki=3
383 !
384  if(kl.eq.1) then
385  ilayer=ilayer+1
386  if(ilayer.gt.1) then
387  do k=1,3
388  dlayer(k)=dlayer(k)+xlayer(ilayer-1,k)
389  enddo
390  endif
391  endif
392  ze=2.d0*(dlayer(ki)+(ze+1.d0)/
393  & 2.d0*xlayer(ilayer,ki))/tlayer(ki)-1.d0
394  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
395 !
396 ! material and orientation
397 !
398  imat=ielmat(ilayer,i)
399  amat=matname(imat)
400  if(norien.gt.0) then
401  iorien=ielorien(ilayer,i)
402  else
403  iorien=0
404  endif
405 !
406  if(nelcon(1,imat).lt.0) then
407  ihyper=1
408  else
409  ihyper=0
410  endif
411  endif
412  elseif(lakonl(4:4).eq.'6') then
413  xi=gauss3d7(1,jj)
414  et=gauss3d7(2,jj)
415  ze=gauss3d7(3,jj)
416  weight=weight3d7(jj)
417  endif
418 !
419 c Bernhardi start
420  if(lakonl(1:5).eq.'C3D8R') then
421  call shape8hr(xl,xsj,shp,gs,a)
422  elseif(lakonl(1:5).eq.'C3D8I') then
423  call shape8hu(xi,et,ze,xl,xsj,shp,iflag)
424  elseif(nope.eq.20) then
425 c Bernhardi end
426  if(lakonl(7:7).eq.'A') then
427  call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag)
428  elseif((lakonl(7:7).eq.'E').or.
429  & (lakonl(7:7).eq.'S')) then
430  call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag)
431  else
432  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
433  endif
434  elseif(nope.eq.26) then
435  call shape26h(xi,et,ze,xl,xsj,shp,iflag,konl)
436  elseif(nope.eq.8) then
437  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
438  elseif(nope.eq.10) then
439  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
440  elseif(nope.eq.14) then
441  call shape14tet(xi,et,ze,xl,xsj,shp,iflag,konl)
442  elseif(nope.eq.4) then
443  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
444  elseif(nope.eq.15) then
445  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
446  else
447  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
448  endif
449 !
450 ! calculation of the objective function
451 !
452  rho=rhcon(1,1,imat)
453  xmassel=xmassel+weight*xsj*rho
454 !
455 ! end of loop over all integration points
456  enddo
457 !
458 ! summation of the objective function over all elements
459 !
460  if(idesvar.eq.0) then
461  xmass(i)=xmassel
462  g0(iobject)=g0(iobject)+xmassel
463  else
464  dgdx(idesvar,iobject)=dgdx(idesvar,iobject)+
465  & (xmassel-xmass(i))/distmin
466  endif
467 !
468 ! end of loop over all elements in thread
469 !
470  enddo
471 !
472 !
473  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
shape8q
subroutine shape8q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape8q.f:20
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
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
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)