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

Go to the source code of this file.

Functions/Subroutines

subroutine e_c3d_rhs_th (co, nk, konl, lakonl, ff, nelem, nmethod, t0, t1, vold, nelemload, sideload, xload, nload, idist, dtime, ttime, time, istep, iinc, xloadold, reltime, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, mi, ielprop, prop, sti, xstateini, xstate, nstate_)
 

Function/Subroutine Documentation

◆ e_c3d_rhs_th()

subroutine e_c3d_rhs_th ( real*8, dimension(3,*)  co,
integer  nk,
integer, dimension(20)  konl,
character*8  lakonl,
real*8, dimension(60)  ff,
integer  nelem,
integer  nmethod,
real*8, dimension(*)  t0,
real*8, dimension(*)  t1,
real*8, dimension(0:mi(2),*)  vold,
integer, dimension(2,*)  nelemload,
character*20, dimension(*)  sideload,
real*8, dimension(2,*)  xload,
integer  nload,
integer  idist,
real*8  dtime,
real*8  ttime,
real*8  time,
integer  istep,
integer  iinc,
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,
integer, dimension(*)  mi,
integer, dimension(*)  ielprop,
real*8, dimension(*)  prop,
real*8, dimension(6,mi(1),*)  sti,
real*8, dimension(nstate_,mi(1),*)  xstateini,
real*8, dimension(nstate_,mi(1),*)  xstate,
integer  nstate_ 
)
25 !
26 ! computation of the rhs for the element with
27 ! the topology in konl
28 !
29 ! ff: rhs without temperature and eigenstress contribution
30 !
31  implicit none
32 !
33  logical ivolumeforce
34 !
35  character*8 lakonl
36  character*20 sideload(*)
37 !
38  integer konl(20),ifaceq(8,6),nelemload(2,*),nk,nelem,nmethod,
39  & nload,idist,i,j,k,i1,iflag,ipompc(*),nodempc(3,*),nmpc,
40  & jj,id,ipointer,ig,kk,nope,nopes,mint2d,ikmpc(*),ilmpc(*),
41  & mint3d,ifacet(6,4),nopev,ifacew(8,5),iinc,istep,jltyp,
42  & iscale,mi(*),ielprop(*),null,nstate_
43 !
44  real*8 co(3,*),xl(3,20),shp(4,20),xs2(3,7),xloadold(2,*),
45  & ff(60),shpj(4,20),dxsj2,temp,press,t0(*),t1(*),coords(3),
46  & xl2(3,8),xsj2(3),shp2(7,8),vold(0:mi(2),*),xload(2,*),
47  & xi,et,ze,xsj,xsjj,t1l,ttime,time,weight,pgauss(3),tvar(2),
48  & reltime,areaj,coefmpc(*),tl2(8),prop(*),sti(6,mi(1),*),
49  & xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*)
50 !
51  real*8 dtime
52 !
53  include "gauss.f"
54 !
55  data ifaceq /4,3,2,1,11,10,9,12,
56  & 5,6,7,8,13,14,15,16,
57  & 1,2,6,5,9,18,13,17,
58  & 2,3,7,6,10,19,14,18,
59  & 3,4,8,7,11,20,15,19,
60  & 4,1,5,8,12,17,16,20/
61  data ifacet /1,3,2,7,6,5,
62  & 1,2,4,5,9,8,
63  & 2,3,4,6,10,9,
64  & 1,4,3,8,10,7/
65  data ifacew /1,3,2,9,8,7,0,0,
66  & 4,5,6,10,11,12,0,0,
67  & 1,2,5,4,7,14,10,13,
68  & 2,3,6,5,8,15,11,14,
69  & 4,6,3,1,12,15,9,13/
70  data iflag /3/
71  data null /0/
72 !
73  tvar(1)=time
74  tvar(2)=ttime+time
75 !
76  if(lakonl(4:4).eq.'2') then
77  nope=20
78  nopev=8
79  nopes=8
80  elseif(lakonl(4:4).eq.'8') then
81  nope=8
82  nopev=8
83  nopes=4
84  elseif(lakonl(4:5).eq.'10') then
85  nope=10
86  nopev=4
87  nopes=6
88  elseif(lakonl(4:4).eq.'4') then
89  nope=4
90  nopev=4
91  nopes=3
92  elseif(lakonl(4:5).eq.'15') then
93  nope=15
94  nopev=6
95  else
96  nope=6
97  nopev=6
98  endif
99 !
100  if(lakonl(4:5).eq.'8R') then
101  mint2d=1
102  mint3d=1
103  elseif(lakonl(4:7).eq.'20RB') then
104  if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
105  mint2d=4
106  mint3d=50
107  else
108  mint2d=4
109  call beamintscheme(lakonl,mint3d,ielprop(nelem),prop,
110  & null,xi,et,ze,weight)
111  endif
112  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) then
113  mint2d=4
114  mint3d=8
115  elseif(lakonl(4:4).eq.'2') then
116  mint2d=9
117  mint3d=27
118  elseif(lakonl(4:5).eq.'10') then
119  mint2d=3
120  mint3d=4
121  elseif(lakonl(4:4).eq.'4') then
122  mint2d=1
123  mint3d=1
124  elseif(lakonl(4:5).eq.'15') then
125  mint3d=9
126  else
127  mint3d=2
128  endif
129 !
130 ! computation of the coordinates of the local nodes
131 !
132  do i=1,nope
133  do j=1,3
134  xl(j,i)=co(j,konl(i))
135  enddo
136  enddo
137 !
138 ! initialisation for distributed forces
139 !
140 c if(idist.ne.0) then
141  do i=1,nope
142  ff(i)=0.d0
143  enddo
144 c endif
145 !
146 ! computation of the body forces
147 !
148  ivolumeforce=.false.
149 c if(nload.gt.0) then
150  call nident2(nelemload,nelem,nload,id)
151  do
152  if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
153  if(sideload(id)(1:2).ne.'BF') then
154  id=id-1
155  cycle
156  else
157  ivolumeforce=.true.
158  exit
159  endif
160  enddo
161 c endif
162 !
163 ! computation of the matrix: loop over the Gauss points
164 !
165  if(ivolumeforce) then
166  do kk=1,mint3d
167  if(lakonl(4:5).eq.'8R') then
168  xi=gauss3d1(1,kk)
169  et=gauss3d1(2,kk)
170  ze=gauss3d1(3,kk)
171  weight=weight3d1(kk)
172  elseif(lakonl(4:7).eq.'20RB') then
173  if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
174  xi=gauss3d13(1,kk)
175  et=gauss3d13(2,kk)
176  ze=gauss3d13(3,kk)
177  weight=weight3d13(kk)
178  else
179  call beamintscheme(lakonl,mint3d,ielprop(nelem),prop,
180  & kk,xi,et,ze,weight)
181  endif
182  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R'))
183  & then
184  xi=gauss3d2(1,kk)
185  et=gauss3d2(2,kk)
186  ze=gauss3d2(3,kk)
187  weight=weight3d2(kk)
188  elseif(lakonl(4:4).eq.'2') then
189  xi=gauss3d3(1,kk)
190  et=gauss3d3(2,kk)
191  ze=gauss3d3(3,kk)
192  weight=weight3d3(kk)
193  elseif(lakonl(4:5).eq.'10') then
194  xi=gauss3d5(1,kk)
195  et=gauss3d5(2,kk)
196  ze=gauss3d5(3,kk)
197  weight=weight3d5(kk)
198  elseif(lakonl(4:4).eq.'4') then
199  xi=gauss3d4(1,kk)
200  et=gauss3d4(2,kk)
201  ze=gauss3d4(3,kk)
202  weight=weight3d4(kk)
203  elseif(lakonl(4:5).eq.'15') then
204  xi=gauss3d8(1,kk)
205  et=gauss3d8(2,kk)
206  ze=gauss3d8(3,kk)
207  weight=weight3d8(kk)
208  else
209  xi=gauss3d7(1,kk)
210  et=gauss3d7(2,kk)
211  ze=gauss3d7(3,kk)
212  weight=weight3d7(kk)
213  endif
214 !
215 ! calculation of the shape functions and their derivatives
216 ! in the gauss point
217 !
218  if(nope.eq.20) then
219  if(lakonl(7:7).eq.'A') then
220  call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag)
221  elseif((lakonl(7:7).eq.'E').or.(lakonl(7:7).eq.'S')) then
222  call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag)
223  else
224  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
225  endif
226  elseif(nope.eq.8) then
227  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
228  elseif(nope.eq.10) then
229  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
230  elseif(nope.eq.4) then
231  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
232  elseif(nope.eq.15) then
233  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
234  else
235  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
236  endif
237 !
238 ! check the jacobian determinant
239 !
240  if(xsj.lt.1.d-20) then
241  write(*,*) '*ERROR in e_c3d_rhs_th: nonpositive jacobian'
242  write(*,*) ' determinant in element',nelem
243  write(*,*)
244  xsj=dabs(xsj)
245  nmethod=0
246  endif
247 !
248 ! calculating the temperature in the integration
249 ! point
250 !
251  t1l=0.d0
252 !
253  do i1=1,nope
254  t1l=t1l+shp(4,i1)*vold(0,konl(i1))
255  enddo
256 !
257 ! incorporating the jacobian determinant in the shape
258 ! functions
259 !
260  xsjj=dsqrt(xsj)
261  do i1=1,nope
262  shpj(1,i1)=shp(1,i1)*xsjj
263  shpj(2,i1)=shp(2,i1)*xsjj
264  shpj(3,i1)=shp(3,i1)*xsjj
265  shpj(4,i1)=shp(4,i1)*xsj
266  enddo
267 !
268 ! computation of the right hand side
269 !
270 ! distributed heat flux
271 !
272 c if(nload.gt.0) then
273  call nident2(nelemload,nelem,nload,id)
274  areaj=xsj*weight
275  do
276  if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
277  if(sideload(id)(1:2).ne.'BF') then
278  id=id-1
279  cycle
280  endif
281  if(sideload(id)(3:4).eq.'NU') then
282  do j=1,3
283  pgauss(j)=0.d0
284  do i1=1,nope
285  pgauss(j)=pgauss(j)+
286  & shp(4,i1)*co(j,konl(i1))
287  enddo
288  enddo
289  jltyp=1
290  iscale=1
291  call dflux(xload(1,id),t1l,istep,iinc,tvar,
292  & nelem,kk,pgauss,jltyp,temp,press,sideload(id),
293  & areaj,vold,co,lakonl,konl,ipompc,nodempc,coefmpc,
294  & nmpc,ikmpc,ilmpc,iscale,mi)
295  if((nmethod.eq.1).and.(iscale.ne.0))
296  & xload(1,id)=xloadold(1,id)+
297  & (xload(1,id)-xloadold(1,id))*reltime
298  endif
299  do jj=1,nope
300  ff(jj)=ff(jj)+xload(1,id)*shpj(4,jj)*weight
301  enddo
302  exit
303  enddo
304 c endif
305 !
306  enddo
307  endif
308 !
309 ! distributed loads
310 !
311 c if(nload.eq.0) then
312 c return
313 c endif
314  call nident2(nelemload,nelem,nload,id)
315  do
316  if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
317  if((sideload(id)(1:1).ne.'F').and.
318  & (sideload(id)(1:1).ne.'R').and.
319  & (sideload(id)(1:1).ne.'S')) then
320  id=id-1
321  cycle
322  endif
323  read(sideload(id)(2:2),'(i1)') ig
324 !
325 ! treatment of wedge faces
326 !
327  if(lakonl(4:4).eq.'6') then
328  mint2d=1
329  if(ig.le.2) then
330  nopes=3
331  else
332  nopes=4
333  endif
334  endif
335  if(lakonl(4:5).eq.'15') then
336  if(ig.le.2) then
337  mint2d=3
338  nopes=6
339  else
340  mint2d=4
341  nopes=8
342  endif
343  endif
344 !
345  if((nope.eq.20).or.(nope.eq.8)) then
346  do i=1,nopes
347  tl2(i)=vold(0,konl(ifaceq(i,ig)))
348  do j=1,3
349  xl2(j,i)=co(j,konl(ifaceq(i,ig)))+
350  & vold(j,konl(ifaceq(i,ig)))
351  enddo
352  enddo
353  elseif((nope.eq.10).or.(nope.eq.4)) then
354  do i=1,nopes
355  tl2(i)=vold(0,konl(ifacet(i,ig)))
356  do j=1,3
357  xl2(j,i)=co(j,konl(ifacet(i,ig)))+
358  & vold(j,konl(ifacet(i,ig)))
359  enddo
360  enddo
361  else
362  do i=1,nopes
363  tl2(i)=vold(0,konl(ifacew(i,ig)))
364  do j=1,3
365  xl2(j,i)=co(j,konl(ifacew(i,ig)))+
366  & vold(j,konl(ifacew(i,ig)))
367  enddo
368  enddo
369  endif
370 !
371  do i=1,mint2d
372  if((lakonl(4:5).eq.'8R').or.
373  & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
374  xi=gauss2d1(1,i)
375  et=gauss2d1(2,i)
376  weight=weight2d1(i)
377  elseif((lakonl(4:4).eq.'8').or.
378  & (lakonl(4:6).eq.'20R').or.
379  & ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
380  xi=gauss2d2(1,i)
381  et=gauss2d2(2,i)
382  weight=weight2d2(i)
383  elseif(lakonl(4:4).eq.'2') then
384  xi=gauss2d3(1,i)
385  et=gauss2d3(2,i)
386  weight=weight2d3(i)
387  elseif((lakonl(4:5).eq.'10').or.
388  & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
389  xi=gauss2d5(1,i)
390  et=gauss2d5(2,i)
391  weight=weight2d5(i)
392  elseif((lakonl(4:4).eq.'4').or.
393  & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
394  xi=gauss2d4(1,i)
395  et=gauss2d4(2,i)
396  weight=weight2d4(i)
397  endif
398 !
399  if(nopes.eq.8) then
400  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
401  elseif(nopes.eq.4) then
402  call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
403  elseif(nopes.eq.6) then
404  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
405  else
406  call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
407  endif
408 !
409  dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
410  & xsj2(3)*xsj2(3))
411  areaj=dxsj2*weight
412 !
413  temp=0.d0
414  do j=1,nopes
415  temp=temp+tl2(j)*shp2(4,j)
416  enddo
417 !
418 ! for nonuniform load: determine the coordinates of the
419 ! point (transferred into the user subroutine)
420 !
421  if(sideload(id)(3:4).eq.'NU') then
422  do k=1,3
423  coords(k)=0.d0
424  do j=1,nopes
425  coords(k)=coords(k)+xl2(k,j)*shp2(4,j)
426  enddo
427  enddo
428  read(sideload(id)(2:2),'(i1)') jltyp
429  jltyp=jltyp+10
430  if(sideload(id)(1:1).eq.'S') then
431  iscale=1
432  call dflux(xload(1,id),temp,istep,iinc,tvar,
433  & nelem,i,coords,jltyp,temp,press,sideload(id),
434  & areaj,vold,co,lakonl,konl,ipompc,nodempc,
435  & coefmpc,nmpc,ikmpc,ilmpc,iscale,mi)
436  if((nmethod.eq.1).and.(iscale.ne.0))
437  & xload(1,id)=xloadold(1,id)+
438  & (xload(1,id)-xloadold(1,id))*reltime
439  endif
440  endif
441 !
442  do k=1,nopes
443  if((nope.eq.20).or.(nope.eq.8)) then
444  ipointer=ifaceq(k,ig)
445  elseif((nope.eq.10).or.(nope.eq.4)) then
446  ipointer=ifacet(k,ig)
447  else
448  ipointer=ifacew(k,ig)
449  endif
450  if(sideload(id)(1:1).eq.'S') then
451 !
452 ! flux INTO the face is positive (input deck convention)
453 ! this is different from the convention in the theory
454 !
455  ff(ipointer)=ff(ipointer)+shp2(4,k)*xload(1,id)
456  & *dxsj2*weight
457  elseif(sideload(id)(1:1).eq.'F') then
458  write(*,*) '*ERROR in e_c3d_rhs_th.f: no'
459  write(*,*) ' film conditions allowed'
460  write(*,*) ' in an modal dynamic calculation'
461  call exit(201)
462  elseif(sideload(id)(1:1).eq.'R') then
463  write(*,*) '*ERROR in e_c3d_rhs_th.f: no'
464  write(*,*) ' radiation conditions allowed'
465  write(*,*) ' in an modal dynamic calculation'
466  call exit(201)
467  endif
468  enddo
469 !
470  enddo
471 !
472  id=id-1
473  enddo
474 !
475  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
dflux
subroutine dflux(flux, sol, kstep, kinc, time, noel, npt, coords, jltyp, temp, press, loadtype, area, vold, co, lakonl, konl, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, iscale, mi, sti, xstateini, xstate, nstate_, dtime)
Definition: dflux.f:23
shape6tri
subroutine shape6tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape6tri.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)