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

Go to the source code of this file.

Functions/Subroutines

subroutine calcenergy (ipkon, lakon, kon, co, ener, mi, ne, thicke, ielmat, energyini, energy, ielprop, prop)
 

Function/Subroutine Documentation

◆ calcenergy()

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