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

Go to the source code of this file.

Functions/Subroutines

subroutine calcmass (ipkon, lakon, kon, co, mi, nelem, ne, thicke, ielmat, nope, t0, t1, rhcon, nrhcon, ntmat_, ithermal, csmasstot, ielprop, prop)
 

Function/Subroutine Documentation

◆ calcmass()

subroutine calcmass ( integer, dimension(*)  ipkon,
character*8, dimension(*)  lakon,
integer, dimension(*)  kon,
real*8, dimension(3,*)  co,
integer, dimension(*)  mi,
integer  nelem,
integer  ne,
real*8, dimension(mi(3),*)  thicke,
integer, dimension(mi(3),*)  ielmat,
integer  nope,
real*8, dimension(*)  t0,
real*8, dimension(*)  t1,
real*8, dimension(0:1,ntmat_,*)  rhcon,
integer, dimension(*)  nrhcon,
integer  ntmat_,
integer  ithermal,
real*8  csmasstot,
integer, dimension(*)  ielprop,
real*8, dimension(*)  prop 
)
22 !
23 ! calculates the mass of element "nelem"
24 !
25  implicit none
26 !
27  character*8 lakon(*)
28 !
29  integer ipkon(*),nelem,kon(*),mi(*),nope,indexe,i,j,k,konl(20),
30  & mint3d,jj,iflag,ne,ki,kl,ilayer,nlayer,kk,ntmat_,i1,imat,
31  & nopes,ielmat(mi(3),*),mint2d,nrhcon(*),ithermal,null,ielprop(*)
32 !
33  real*8 csmasstot,csmass,co(3,*),prop(*),
34  & xl(3,20),xi,et,ze,xsj,shp(4,20),weight,
35  & a,gs(8,4),dlayer(4),tlayer(4),thickness,t0l,t1l,rho,
36  & thicke(mi(3),*),xlayer(mi(3),4),shp2(7,8),xs2(3,7),xsj2(3),
37  & xl2(3,8),t0(*),t1(*),rhcon(0:1,ntmat_,*)
38 !
39  include "gauss.f"
40 !
41  data iflag /2/
42  null=0
43 !
44  if(ipkon(nelem).lt.0) return
45 !
46  indexe=ipkon(nelem)
47 !
48  if(lakon(nelem)(1:5).eq.'C3D8I') then
49  nope=11
50  elseif(lakon(nelem)(4:4).eq.'2') then
51  nope=20
52  elseif(lakon(nelem)(4:4).eq.'8') then
53  nope=8
54  elseif(lakon(nelem)(4:5).eq.'10') then
55  nope=10
56  elseif(lakon(nelem)(4:4).eq.'4') then
57  nope=4
58  elseif(lakon(nelem)(4:5).eq.'15') then
59  nope=15
60  elseif(lakon(nelem)(4:5).eq.'6') then
61  nope=6
62  else
63  nope=0
64  endif
65 !
66 ! composite materials
67 !
68  if(lakon(nelem)(7:8).ne.'LC') then
69  imat=ielmat(1,nelem)
70  else
71 !
72 ! determining the number of layers
73 !
74  nlayer=0
75  do k=1,mi(3)
76  if(ielmat(k,nelem).ne.0) then
77  nlayer=nlayer+1
78  endif
79  enddo
80 !
81  if(lakon(nelem)(4:4).eq.'2') then
82  mint2d=4
83  nopes=8
84 !
85 ! determining the layer thickness and global thickness
86 ! at the shell integration points
87 !
88  iflag=1
89  indexe=ipkon(nelem)
90  do kk=1,mint2d
91  xi=gauss3d2(1,kk)
92  et=gauss3d2(2,kk)
93  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
94  tlayer(kk)=0.d0
95  do i=1,nlayer
96  thickness=0.d0
97  do j=1,nopes
98  thickness=thickness+thicke(i,indexe+j)*shp2(4,j)
99  enddo
100  tlayer(kk)=tlayer(kk)+thickness
101  xlayer(i,kk)=thickness
102  enddo
103  enddo
104  iflag=2
105 !
106  ilayer=0
107  do i=1,4
108  dlayer(i)=0.d0
109  enddo
110  elseif(lakon(nelem)(4:5).eq.'15') then
111  mint2d=3
112  nopes=6
113 !
114 ! determining the layer thickness and global thickness
115 ! at the shell integration points
116 !
117  iflag=1
118  indexe=ipkon(nelem)
119  do kk=1,mint2d
120  xi=gauss3d10(1,kk)
121  et=gauss3d10(2,kk)
122  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
123  tlayer(kk)=0.d0
124  do i=1,nlayer
125  thickness=0.d0
126  do j=1,nopes
127  thickness=thickness+thicke(i,indexe+j)*shp2(4,j)
128  enddo
129  tlayer(kk)=tlayer(kk)+thickness
130  xlayer(i,kk)=thickness
131  enddo
132  enddo
133  iflag=2
134 !
135  ilayer=0
136  do i=1,3
137  dlayer(i)=0.d0
138  enddo
139  endif
140 !
141  endif
142 !
143  do j=1,nope
144  konl(j)=kon(indexe+j)
145  do k=1,3
146  xl(k,j)=co(k,konl(j))
147  enddo
148  enddo
149 !
150  csmass=0.d0
151 !
152  if(lakon(nelem)(4:5).eq.'8R') then
153  mint3d=1
154  elseif(lakon(nelem)(4:7).eq.'20RB') then
155  if((lakon(nelem)(8:8).eq.'R').or.
156  & (lakon(nelem)(8:8).eq.'C')) then
157  mint3d=50
158  else
159  call beamintscheme(lakon(nelem),mint3d,ielprop(nelem),prop,
160  & null,xi,et,ze,weight)
161  endif
162  elseif((lakon(nelem)(4:4).eq.'8').or.
163  & (lakon(nelem)(4:6).eq.'20R')) then
164  if(lakon(nelem)(7:8).eq.'LC') then
165  mint3d=8*nlayer
166  else
167  mint3d=8
168  endif
169  elseif(lakon(nelem)(4:4).eq.'2') then
170  mint3d=27
171  elseif(lakon(nelem)(4:5).eq.'10') then
172  mint3d=4
173  elseif(lakon(nelem)(4:4).eq.'4') then
174  mint3d=1
175  elseif(lakon(nelem)(4:5).eq.'15') then
176  if(lakon(nelem)(7:8).eq.'LC') then
177  mint3d=6*nlayer
178  else
179  mint3d=9
180  endif
181  elseif(lakon(nelem)(4:5).eq.'6') then
182  mint3d=2
183  else
184  mint3d=0
185  endif
186 !
187  do jj=1,mint3d
188  if(lakon(nelem)(4:5).eq.'8R') then
189  xi=gauss3d1(1,jj)
190  et=gauss3d1(2,jj)
191  ze=gauss3d1(3,jj)
192  weight=weight3d1(jj)
193  elseif(lakon(nelem)(4:7).eq.'20RB') then
194  if((lakon(nelem)(8:8).eq.'R').or.
195  & (lakon(nelem)(8:8).eq.'C')) then
196  xi=gauss3d13(1,kk)
197  et=gauss3d13(2,kk)
198  ze=gauss3d13(3,kk)
199  weight=weight3d13(kk)
200  else
201  call beamintscheme(lakon(nelem),mint3d,ielprop(nelem),
202  & prop,kk,xi,et,ze,weight)
203  endif
204  elseif((lakon(nelem)(4:4).eq.'8').or.
205  & (lakon(nelem)(4:6).eq.'20R'))
206  & then
207  if(lakon(nelem)(7:8).ne.'LC') then
208  xi=gauss3d2(1,jj)
209  et=gauss3d2(2,jj)
210  ze=gauss3d2(3,jj)
211  weight=weight3d2(jj)
212  else
213  kl=mod(jj,8)
214  if(kl.eq.0) kl=8
215 !
216  xi=gauss3d2(1,kl)
217  et=gauss3d2(2,kl)
218  ze=gauss3d2(3,kl)
219  weight=weight3d2(kl)
220 !
221  ki=mod(jj,4)
222  if(ki.eq.0) ki=4
223 !
224  if(kl.eq.1) then
225  ilayer=ilayer+1
226  if(ilayer.gt.1) then
227  do i=1,4
228  dlayer(i)=dlayer(i)+xlayer(ilayer-1,i)
229  enddo
230  endif
231  endif
232  ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/
233  & tlayer(ki)-1.d0
234  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
235  imat=ielmat(ilayer,nelem)
236  endif
237  elseif(lakon(nelem)(4:4).eq.'2') then
238  xi=gauss3d3(1,jj)
239  et=gauss3d3(2,jj)
240  ze=gauss3d3(3,jj)
241  weight=weight3d3(jj)
242  elseif(lakon(nelem)(4:5).eq.'10') then
243  xi=gauss3d5(1,jj)
244  et=gauss3d5(2,jj)
245  ze=gauss3d5(3,jj)
246  weight=weight3d5(jj)
247  elseif(lakon(nelem)(4:4).eq.'4') then
248  xi=gauss3d4(1,jj)
249  et=gauss3d4(2,jj)
250  ze=gauss3d4(3,jj)
251  weight=weight3d4(jj)
252  elseif(lakon(nelem)(4:5).eq.'15') then
253  if(lakon(nelem)(7:8).ne.'LC') then
254  xi=gauss3d8(1,jj)
255  et=gauss3d8(2,jj)
256  ze=gauss3d8(3,jj)
257  weight=weight3d8(jj)
258  else
259  kl=mod(jj,6)
260  if(kl.eq.0) kl=6
261 !
262  xi=gauss3d10(1,kl)
263  et=gauss3d10(2,kl)
264  ze=gauss3d10(3,kl)
265  weight=weight3d10(kl)
266 !
267  ki=mod(jj,3)
268  if(ki.eq.0) ki=3
269 !
270  if(kl.eq.1) then
271  ilayer=ilayer+1
272  if(ilayer.gt.1) then
273  do i=1,3
274  dlayer(i)=dlayer(i)+xlayer(ilayer-1,i)
275  enddo
276  endif
277  endif
278  ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/
279  & tlayer(ki)-1.d0
280  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
281  imat=ielmat(ilayer,nelem)
282  endif
283  else
284  xi=gauss3d7(1,jj)
285  et=gauss3d7(2,jj)
286  ze=gauss3d7(3,jj)
287  weight=weight3d7(jj)
288  endif
289 !
290  if(lakon(nelem)(1:5).eq.'C3D8R') then
291  call shape8hr(xl,xsj,shp,gs,a)
292  elseif(lakon(nelem)(1:5).eq.'C3D8I') then
293  call shape8hu(xi,et,ze,xl,xsj,shp,iflag)
294  elseif(nope.eq.20) then
295  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
296  elseif(nope.eq.8) then
297  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
298  elseif(nope.eq.10) then
299  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
300  elseif(nope.eq.4) then
301  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
302  elseif(nope.eq.15) then
303  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
304  else
305  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
306  endif
307 !
308 ! calculating the initial temperature in the integration
309 ! point
310 !
311  t0l=0.d0
312  t1l=0.d0
313  if(ithermal.eq.1) then
314  if(lakon(nelem)(4:5).eq.'8 ') then
315  do i1=1,nope
316  t0l=t0l+t0(konl(i1))/8.d0
317  enddo
318  elseif(lakon(nelem)(4:6).eq.'20 ') then
319  call lintemp(t0,t1,konl,nope,jj,t0l,t1l)
320  elseif(lakon(nelem)(4:6).eq.'10T') then
321  call linscal10(t0,konl,t0l,null,shp)
322  else
323  do i1=1,nope
324  t0l=t0l+shp(4,i1)*t0(konl(i1))
325  enddo
326  endif
327  endif
328 !
329 ! calculating the density
330 !
331  call materialdata_rho(rhcon,nrhcon,imat,rho,
332  & t0l,ntmat_,ithermal)
333 !
334  csmass=csmass+weight*xsj*rho
335  enddo
336 !
337  csmasstot=csmasstot+csmass
338 c write(*,*) 'calcmass ',nelem,csmass,csmasstot
339 !
340  return
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
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
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
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)