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

Go to the source code of this file.

Functions/Subroutines

subroutine e_c3d_rhs (co, nk, konl, lakonl, p1, p2, omx, bodyfx, nbody, ff, nelem, nmethod, rhcon, ielmat, ntmat_, vold, iperturb, nelemload, sideload, xload, nload, idist, ttime, time, istep, iinc, dtime, xloadold, reltime, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, veold, matname, mi, ielprop, prop)
 

Function/Subroutine Documentation

◆ e_c3d_rhs()

subroutine e_c3d_rhs ( real*8, dimension(3,*)  co,
integer  nk,
integer, dimension(20)  konl,
character*8  lakonl,
real*8, dimension(3,2)  p1,
real*8, dimension(3,2)  p2,
real*8, dimension(2)  omx,
real*8, dimension(3)  bodyfx,
integer  nbody,
real*8, dimension(60)  ff,
integer  nelem,
integer  nmethod,
real*8, dimension(0:1,ntmat_,*)  rhcon,
integer, dimension(mi(3),*)  ielmat,
integer  ntmat_,
real*8, dimension(0:mi(2),*)  vold,
integer  iperturb,
integer, dimension(2,*)  nelemload,
character*20, dimension(*)  sideload,
real*8, dimension(2,*)  xload,
integer  nload,
integer  idist,
real*8  ttime,
real*8  time,
integer  istep,
integer  iinc,
real*8  dtime,
real*8, dimension(2,*)  xloadold,
real*8  reltime,
integer, dimension(*)  ipompc,
integer, dimension(3,*)  nodempc,
real*8, dimension(*)  coefmpc,
integer  nmpc,
integer, dimension(*)  ikmpc,
integer, dimension(*)  ilmpc,
real*8, dimension(0:mi(2),*)  veold,
character*80, dimension(*)  matname,
integer, dimension(*)  mi,
integer, dimension(*)  ielprop,
real*8, dimension(*)  prop 
)
24 !
25 ! computation of the rhs for the element with
26 ! the topology in konl: only for nonlinear calculations (i.e.
27 ! the field ff contains no temperature and eigenstress
28 ! contributions)
29 !
30 ! RESTRICTION: the only material parameter occurring in the
31 ! loading is the density. In the present subroutine the density
32 ! is assumed to be temperature INDEPENDENT.
33 !
34  implicit none
35 !
36  logical ivolumeforce
37 !
38  character*8 lakonl
39  character*20 sideload(*)
40  character*80 matname(*),amat
41 !
42  integer mi(*),nk,konl(20),ielmat(mi(3),*),nbody,ifaceq(8,6),
43  & nelemload(2,*),ielprop(*),null,
44  & nelem,nmethod,iperturb,nload,idist,i,i1,j1,jj,jj1,id,kk,
45  & ipointer,nope,nopes,j,k,ntmat_,i2,imat,ii,ig,mint2d,mint3d,
46  & ifacet(6,4),ifacew(8,5),istep,iinc,layer,kspt,jltyp,iflag,
47  & ipompc(*),nodempc(3,*),nmpc,ikmpc(*),ilmpc(*),iscale
48 !
49  real*8 co(3,*),p1(3,2),p2(3,2),omx(2),bodyfx(3),veold(0:mi(2),*),
50  & rhcon(0:1,ntmat_,*),xs2(3,7),xloadold(2,*),reltime,
51  & bodyf(3),om(2),rho,bf(3),q(3),shpj(4,20),xl(3,20),
52  & shp(4,20),voldl(3,20),xl2(3,8),xsj2(3),shp2(7,8),
53  & vold(0:mi(2),*),prop(*),
54  & xload(2,*),xi,et,ze,const,xsj,ff(60),weight,ttime,time,tvar(2),
55  & coords(3),dtime,coefmpc(*)
56 !
57  include "gauss.f"
58 !
59  data ifaceq /4,3,2,1,11,10,9,12,
60  & 5,6,7,8,13,14,15,16,
61  & 1,2,6,5,9,18,13,17,
62  & 2,3,7,6,10,19,14,18,
63  & 3,4,8,7,11,20,15,19,
64  & 4,1,5,8,12,17,16,20/
65  data ifacet /1,3,2,7,6,5,
66  & 1,2,4,5,9,8,
67  & 2,3,4,6,10,9,
68  & 1,4,3,8,10,7/
69  data ifacew /1,3,2,9,8,7,0,0,
70  & 4,5,6,10,11,12,0,0,
71  & 1,2,5,4,7,14,10,13,
72  & 2,3,6,5,8,15,11,14,
73  & 4,6,3,1,12,15,9,13/
74  data iflag /3/
75  data null /0/
76 !
77  tvar(1)=time
78  tvar(2)=ttime+time
79 !
80  if(lakonl(4:4).eq.'2') then
81  nope=20
82  nopes=8
83  elseif(lakonl(4:4).eq.'8') then
84  nope=8
85  nopes=4
86  elseif(lakonl(4:5).eq.'10') then
87  nope=10
88  nopes=6
89  elseif(lakonl(4:4).eq.'4') then
90  nope=4
91  nopes=3
92  elseif(lakonl(4:5).eq.'15') then
93  nope=15
94  else
95  nope=6
96  endif
97 !
98  if(lakonl(4:5).eq.'8R') then
99  mint2d=1
100  mint3d=1
101  elseif(lakonl(4:7).eq.'20RB') then
102  if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
103  mint2d=4
104  mint3d=50
105  else
106  mint2d=4
107  call beamintscheme(lakonl,mint3d,ielprop(nelem),prop,
108  & null,xi,et,ze,weight)
109  endif
110  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) then
111  mint2d=4
112  mint3d=8
113  elseif(lakonl(4:4).eq.'2') then
114  mint2d=9
115  mint3d=27
116  elseif(lakonl(4:5).eq.'10') then
117  mint2d=3
118  mint3d=4
119  elseif(lakonl(4:4).eq.'4') then
120  mint2d=1
121  mint3d=1
122  elseif(lakonl(4:5).eq.'15') then
123  mint3d=9
124  else
125  mint3d=2
126  endif
127 !
128 ! computation of the coordinates of the local nodes
129 !
130  do i=1,nope
131  do j=1,3
132  xl(j,i)=co(j,konl(i))
133  enddo
134  enddo
135 !
136 ! initialisation for distributed forces
137 !
138 c if(idist.ne.0) then
139  do i=1,3*nope
140  ff(i)=0.d0
141  enddo
142 c endif
143 !
144 ! displacements for 2nd order static and modal theory
145 !
146  if(iperturb.ne.0) then
147  do i1=1,nope
148  do i2=1,3
149  voldl(i2,i1)=vold(i2,konl(i1))
150  enddo
151  enddo
152  endif
153 !
154 ! calculating the density: no temperature dependence assumed!
155 ! otherwise cfr. e_c3d.f
156 !
157  imat=ielmat(1,nelem)
158  amat=matname(imat)
159  rho=rhcon(1,1,imat)
160 !
161 ! computation of the body forces
162 !
163  ivolumeforce=.false.
164  if(nbody.ne.0) then
165  ivolumeforce=.true.
166  do ii=1,2
167  om(ii)=omx(ii)*rho
168  enddo
169  do ii=1,3
170  bodyf(ii)=bodyfx(ii)*rho
171  enddo
172  elseif(nload.gt.0) then
173  call nident2(nelemload,nelem,nload,id)
174  do
175  if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
176  if((sideload(id)(1:2).ne.'BX').and.
177  & (sideload(id)(1:2).ne.'BY').and.
178  & (sideload(id)(1:2).ne.'BZ')) then
179  id=id-1
180  cycle
181  else
182  ivolumeforce=.true.
183  exit
184  endif
185  enddo
186  endif
187 !
188 ! computation of the rhs: loop over the Gauss points
189 !
190  if(ivolumeforce) then
191  do kk=1,mint3d
192  if(lakonl(4:5).eq.'8R') then
193  xi=gauss3d1(1,kk)
194  et=gauss3d1(2,kk)
195  ze=gauss3d1(3,kk)
196  weight=weight3d1(kk)
197  elseif(lakonl(4:7).eq.'20RB') then
198  if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
199  xi=gauss3d13(1,kk)
200  et=gauss3d13(2,kk)
201  ze=gauss3d13(3,kk)
202  weight=weight3d13(kk)
203  else
204  call beamintscheme(lakonl,mint3d,ielprop(nelem),prop,
205  & kk,xi,et,ze,weight)
206  endif
207  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R'))
208  & then
209  xi=gauss3d2(1,kk)
210  et=gauss3d2(2,kk)
211  ze=gauss3d2(3,kk)
212  weight=weight3d2(kk)
213  elseif(lakonl(4:4).eq.'2') then
214  xi=gauss3d3(1,kk)
215  et=gauss3d3(2,kk)
216  ze=gauss3d3(3,kk)
217  weight=weight3d3(kk)
218  elseif(lakonl(4:5).eq.'10') then
219  xi=gauss3d5(1,kk)
220  et=gauss3d5(2,kk)
221  ze=gauss3d5(3,kk)
222  weight=weight3d5(kk)
223  elseif(lakonl(4:4).eq.'4') then
224  xi=gauss3d4(1,kk)
225  et=gauss3d4(2,kk)
226  ze=gauss3d4(3,kk)
227  weight=weight3d4(kk)
228  elseif(lakonl(4:5).eq.'15') then
229  xi=gauss3d8(1,kk)
230  et=gauss3d8(2,kk)
231  ze=gauss3d8(3,kk)
232  weight=weight3d8(kk)
233  else
234  xi=gauss3d7(1,kk)
235  et=gauss3d7(2,kk)
236  ze=gauss3d7(3,kk)
237  weight=weight3d7(kk)
238  endif
239 !
240 ! calculation of the shape functions and their derivatives
241 ! in the gauss point
242 !
243  if(nope.eq.20) then
244  if(lakonl(7:7).eq.'A') then
245  call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag)
246  elseif((lakonl(7:7).eq.'E').or.(lakonl(7:7).eq.'S')) then
247  call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag)
248  else
249  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
250  endif
251  elseif(nope.eq.8) then
252  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
253  elseif(nope.eq.10) then
254  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
255  elseif(nope.eq.4) then
256  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
257  elseif(nope.eq.15) then
258  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
259  else
260  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
261  endif
262 !
263 ! check the jacobian determinant
264 !
265  if(xsj.lt.1.d-20) then
266  write(*,*) '*ERROR in e_c3d_rhs: nonpositive jacobian'
267  write(*,*) ' determinant in element',nelem
268  write(*,*)
269  xsj=dabs(xsj)
270  nmethod=0
271  endif
272 !
273 ! computation of the right hand side
274 !
275 ! incorporating the jacobian determinant in the shape
276 ! functions
277 !
278  if((nload.gt.0).or.(nbody.ne.0)) then
279  do i1=1,nope
280  shpj(4,i1)=shp(4,i1)*xsj
281  enddo
282  endif
283 !
284 ! distributed body flux
285 !
286  if(nload.gt.0) then
287  call nident2(nelemload,nelem,nload,id)
288  do
289  if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
290  if((sideload(id)(1:2).ne.'BX').and.
291  & (sideload(id)(1:2).ne.'BY').and.
292  & (sideload(id)(1:2).ne.'BZ')) then
293  id=id-1
294  cycle
295  endif
296  if(sideload(id)(3:4).eq.'NU') then
297  do j=1,3
298  coords(j)=0.d0
299  do i1=1,nope
300  coords(j)=coords(j)+
301  & shp(4,i1)*xl(j,i1)
302  enddo
303  enddo
304  if(sideload(id)(1:2).eq.'BX') then
305  jltyp=1
306  elseif(sideload(id)(1:2).eq.'BY') then
307  jltyp=2
308  elseif(sideload(id)(1:2).eq.'BZ') then
309  jltyp=3
310  endif
311  iscale=1
312  call dload(xload(1,id),istep,iinc,tvar,nelem,i,
313  & layer,kspt,coords,jltyp,sideload(id),vold,co,
314  & lakonl,konl,ipompc,nodempc,coefmpc,nmpc,ikmpc,
315  & ilmpc,iscale,veold,rho,amat,mi)
316  if((nmethod.eq.1).and.(iscale.ne.0))
317  & xload(1,id)=xloadold(1,id)+
318  & (xload(1,id)-xloadold(1,id))*reltime
319  endif
320  jj1=1
321  do jj=1,nope
322  if(sideload(id)(1:2).eq.'BX')
323  & ff(jj1)=ff(jj1)+xload(1,id)*shpj(4,jj)*weight
324  if(sideload(id)(1:2).eq.'BY')
325  & ff(jj1+1)=ff(jj1+1)+xload(1,id)*shpj(4,jj)
326  & *weight
327  if(sideload(id)(1:2).eq.'BZ')
328  & ff(jj1+2)=ff(jj1+2)+xload(1,id)*shpj(4,jj)
329  & *weight
330  jj1=jj1+3
331  enddo
332  id=id-1
333  enddo
334  endif
335 !
336 ! body forces
337 !
338  if(nbody.ne.0) then
339 !
340  do i1=1,3
341  bf(i1)=0.d0
342  enddo
343  do jj1=1,2
344  if(dabs(om(jj1)).lt.1.d-20) cycle
345  do i1=1,3
346 !
347 ! computation of the global coordinates of the gauss
348 ! point
349 !
350  q(i1)=0.d0
351  if(iperturb.eq.0) then
352  do j1=1,nope
353  q(i1)=q(i1)+shp(4,j1)*xl(i1,j1)
354  enddo
355  else
356  do j1=1,nope
357  q(i1)=q(i1)+shp(4,j1)*
358  & (xl(i1,j1)+voldl(i1,j1))
359  enddo
360  endif
361 !
362  q(i1)=q(i1)-p1(i1,jj1)
363  enddo
364  const=q(1)*p2(1,jj1)+q(2)*p2(2,jj1)+q(3)*p2(3,jj1)
365 !
366 ! inclusion of the centrifugal force into the body force
367 !
368  do i1=1,3
369  bf(i1)=bf(i1)+(q(i1)-const*p2(i1,jj1))*om(jj1)
370  enddo
371  enddo
372 !
373  do i1=1,3
374  bf(i1)=bf(i1)+bodyf(i1)
375  enddo
376 !
377  jj1=1
378  do jj=1,nope
379  ff(jj1)=ff(jj1)+bf(1)*shpj(4,jj)*weight
380  ff(jj1+1)=ff(jj1+1)+bf(2)*shpj(4,jj)*weight
381  ff(jj1+2)=ff(jj1+2)+bf(3)*shpj(4,jj)*weight
382  jj1=jj1+3
383  enddo
384  endif
385 !
386  enddo
387  endif
388 !
389 ! distributed loads
390 !
391  if(nload.eq.0) then
392  return
393  endif
394  call nident2(nelemload,nelem,nload,id)
395  do
396  if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
397  if(sideload(id)(1:1).ne.'P') then
398  id=id-1
399  cycle
400  endif
401  read(sideload(id)(2:2),'(i1)') ig
402 !
403 ! treatment of wedge faces
404 !
405  if(lakonl(4:4).eq.'6') then
406  mint2d=1
407  if(ig.le.2) then
408  nopes=3
409  else
410  nopes=4
411  endif
412  endif
413  if(lakonl(4:5).eq.'15') then
414  if(ig.le.2) then
415  mint2d=3
416  nopes=6
417  else
418  mint2d=4
419  nopes=8
420  endif
421  endif
422 !
423  if((nope.eq.20).or.(nope.eq.8)) then
424  if(iperturb.eq.0) then
425  do i=1,nopes
426  do j=1,3
427  xl2(j,i)=co(j,konl(ifaceq(i,ig)))
428  enddo
429  enddo
430  else
431  do i=1,nopes
432  do j=1,3
433  xl2(j,i)=co(j,konl(ifaceq(i,ig)))+
434  & vold(j,konl(ifaceq(i,ig)))
435  enddo
436  enddo
437  endif
438  elseif((nope.eq.10).or.(nope.eq.4)) then
439  if(iperturb.eq.0) then
440  do i=1,nopes
441  do j=1,3
442  xl2(j,i)=co(j,konl(ifacet(i,ig)))
443  enddo
444  enddo
445  else
446  do i=1,nopes
447  do j=1,3
448  xl2(j,i)=co(j,konl(ifacet(i,ig)))+
449  & vold(j,konl(ifacet(i,ig)))
450  enddo
451  enddo
452  endif
453  else
454  if(iperturb.eq.0) then
455  do i=1,nopes
456  do j=1,3
457  xl2(j,i)=co(j,konl(ifacew(i,ig)))
458  enddo
459  enddo
460  else
461  do i=1,nopes
462  do j=1,3
463  xl2(j,i)=co(j,konl(ifacew(i,ig)))+
464  & vold(j,konl(ifacew(i,ig)))
465  enddo
466  enddo
467  endif
468  endif
469 !
470  do i=1,mint2d
471  if((lakonl(4:5).eq.'8R').or.
472  & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
473  xi=gauss2d1(1,i)
474  et=gauss2d1(2,i)
475  weight=weight2d1(i)
476  elseif((lakonl(4:4).eq.'8').or.
477  & (lakonl(4:6).eq.'20R').or.
478  & ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
479  xi=gauss2d2(1,i)
480  et=gauss2d2(2,i)
481  weight=weight2d2(i)
482  elseif(lakonl(4:4).eq.'2') then
483  xi=gauss2d3(1,i)
484  et=gauss2d3(2,i)
485  weight=weight2d3(i)
486  elseif((lakonl(4:5).eq.'10').or.
487  & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
488  xi=gauss2d5(1,i)
489  et=gauss2d5(2,i)
490  weight=weight2d5(i)
491  elseif((lakonl(4:4).eq.'4').or.
492  & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
493  xi=gauss2d4(1,i)
494  et=gauss2d4(2,i)
495  weight=weight2d4(i)
496  endif
497 !
498  if(nopes.eq.8) then
499  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
500  elseif(nopes.eq.4) then
501  call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
502  elseif(nopes.eq.6) then
503  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
504  else
505  call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
506  endif
507 !
508 ! for nonuniform load: determine the coordinates of the
509 ! point (transferred into the user subroutine)
510 !
511  if(sideload(id)(3:4).eq.'NU') then
512  do k=1,3
513  coords(k)=0.d0
514  do j=1,nopes
515  coords(k)=coords(k)+xl2(k,j)*shp2(4,j)
516  enddo
517  enddo
518  read(sideload(id)(2:2),'(i1)') jltyp
519  jltyp=jltyp+20
520  iscale=1
521  call dload(xload(1,id),istep,iinc,tvar,nelem,i,layer,
522  & kspt,coords,jltyp,sideload(id),vold,co,lakonl,
523  & konl,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,
524  & iscale,veold,rho,amat,mi)
525  if((nmethod.eq.1).and.(iscale.ne.0))
526  & xload(1,id)=xloadold(1,id)+
527  & (xload(1,id)-xloadold(1,id))*reltime
528  endif
529 !
530  do k=1,nopes
531  if((nope.eq.20).or.(nope.eq.8)) then
532  ipointer=(ifaceq(k,ig)-1)*3
533  elseif((nope.eq.10).or.(nope.eq.4)) then
534  ipointer=(ifacet(k,ig)-1)*3
535  else
536  ipointer=(ifacew(k,ig)-1)*3
537  endif
538  ff(ipointer+1)=ff(ipointer+1)-shp2(4,k)*xload(1,id)
539  & *xsj2(1)*weight
540  ff(ipointer+2)=ff(ipointer+2)-shp2(4,k)*xload(1,id)
541  & *xsj2(2)*weight
542  ff(ipointer+3)=ff(ipointer+3)-shp2(4,k)*xload(1,id)
543  & *xsj2(3)*weight
544  enddo
545  enddo
546 !
547  id=id-1
548  enddo
549 !
550  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
shape3tri
subroutine shape3tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape3tri.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
nident2
subroutine nident2(x, px, n, id)
Definition: nident2.f:27
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
idist
static ITG * idist
Definition: radflowload.c:39
shape4q
subroutine shape4q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape4q.f:20
shape20h
subroutine shape20h(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape20h.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
shape6tri
subroutine shape6tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape6tri.f:20
dload
subroutine dload(f, kstep, kinc, time, noel, npt, layer, kspt, coords, jltyp, loadtype, vold, co, lakonl, konl, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, iscale, veold, rho, amat, mi)
Definition: dload.f:23
Hosted by OpenAircraft.com, (Michigan UAV, LLC)