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

Go to the source code of this file.

Functions/Subroutines

subroutine e_c3d_th (co, nk, kon, lakonl, s, sm, ff, nelem, nmethod, rhcon, nrhcon, ielmat, ielorien, norien, orab, ntmat_, t0, t1, ithermal, vold, iperturb, nelemload, sideload, xload, nload, idist, iexpl, dtime, matname, mi, mass, stiffness, buckling, rhsi, intscheme, physcon, shcon, nshcon, cocon, ncocon, ttime, time, istep, iinc, xstiff, xloadold, reltime, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, springarea, plkcon, nplkcon, npmat_, ncmat_, elcon, nelcon, lakon, pslavsurf, pmastsurf, mortar, clearini, plicon, nplicon, ipkon, ielprop, prop, iponoel, inoel, sti, xstateini, xstate, nstate_, network, ipobody, xbody, ibody)
 

Function/Subroutine Documentation

◆ e_c3d_th()

subroutine e_c3d_th ( real*8, dimension(3,*), intent(in)  co,
integer, intent(in)  nk,
integer, dimension(*), intent(in)  kon,
character*8, intent(in)  lakonl,
real*8, dimension(60,60), intent(inout)  s,
real*8, dimension(60,60), intent(inout)  sm,
real*8, dimension(60), intent(inout)  ff,
integer, intent(in)  nelem,
integer, intent(inout)  nmethod,
real*8, dimension(0:1,ntmat_,*), intent(in)  rhcon,
integer, dimension(*), intent(in)  nrhcon,
integer, dimension(mi(3),*), intent(in)  ielmat,
integer, dimension(mi(3),*), intent(in)  ielorien,
integer, intent(in)  norien,
real*8, dimension(7,*), intent(in)  orab,
integer, intent(in)  ntmat_,
real*8, dimension(*), intent(in)  t0,
real*8, dimension(*), intent(in)  t1,
integer, dimension(2), intent(in)  ithermal,
real*8, dimension(0:mi(2),*), intent(in)  vold,
integer, intent(in)  iperturb,
integer, dimension(2,*), intent(in)  nelemload,
character*20, dimension(*), intent(in)  sideload,
real*8, dimension(2,*), intent(inout)  xload,
integer, intent(in)  nload,
integer, intent(in)  idist,
integer, intent(in)  iexpl,
real*8, intent(in)  dtime,
character*80, dimension(*), intent(in)  matname,
integer, dimension(*), intent(in)  mi,
integer, intent(in)  mass,
integer, intent(in)  stiffness,
integer, intent(in)  buckling,
integer, intent(in)  rhsi,
integer, intent(in)  intscheme,
real*8, dimension(*), intent(in)  physcon,
real*8, dimension(0:3,ntmat_,*), intent(in)  shcon,
integer, dimension(*), intent(in)  nshcon,
real*8, dimension(0:6,ntmat_,*), intent(in)  cocon,
integer, dimension(2,*), intent(in)  ncocon,
real*8, intent(in)  ttime,
real*8, intent(in)  time,
integer, intent(in)  istep,
integer, intent(in)  iinc,
real*8, dimension(27,mi(1),*), intent(in)  xstiff,
real*8, dimension(2,*), intent(in)  xloadold,
real*8, intent(in)  reltime,
integer, dimension(*), intent(in)  ipompc,
integer, dimension(3,*), intent(in)  nodempc,
real*8, dimension(*), intent(in)  coefmpc,
integer, intent(in)  nmpc,
integer, dimension(*), intent(in)  ikmpc,
integer, dimension(*), intent(in)  ilmpc,
real*8, dimension(2,*), intent(inout)  springarea,
real*8, dimension(0:2*npmat_,ntmat_,*), intent(in)  plkcon,
integer, dimension(0:ntmat_,*), intent(in)  nplkcon,
integer, intent(in)  npmat_,
integer, intent(in)  ncmat_,
real*8, dimension(0:ncmat_,ntmat_,*), intent(in)  elcon,
integer, dimension(2,*), intent(in)  nelcon,
character*8, dimension(*), intent(in)  lakon,
real*8, dimension(3,*), intent(in)  pslavsurf,
real*8, dimension(2,*), intent(in)  pmastsurf,
integer, intent(in)  mortar,
real*8, dimension(3,9,*), intent(in)  clearini,
real*8, dimension(0:2*npmat_,ntmat_,*), intent(in)  plicon,
integer, dimension(0:ntmat_,*), intent(in)  nplicon,
integer, dimension(*), intent(in)  ipkon,
integer, dimension(*), intent(in)  ielprop,
real*8, dimension(*), intent(in)  prop,
integer, dimension(*)  iponoel,
integer, dimension(2,*)  inoel,
real*8, dimension(6,mi(1),*)  sti,
real*8, dimension(nstate_,mi(1),*)  xstateini,
real*8, dimension(nstate_,mi(1),*)  xstate,
integer  nstate_,
integer  network,
integer, dimension(2,*)  ipobody,
real*8, dimension(7,*)  xbody,
integer, dimension(3,*)  ibody 
)
31 !
32 ! computation of the element matrix and rhs for the element with
33 ! the topology in konl
34 !
35 ! ff: rhs without temperature and eigenstress contribution
36 !
37 ! nmethod=0: check for positive Jacobian
38 ! nmethod=1: stiffness matrix + right hand side
39 ! nmethod=2: stiffness matrix + mass matrix
40 ! nmethod=3: static stiffness + buckling stiffness
41 ! nmethod=4: stiffness matrix + mass matrix
42 !
43  implicit none
44 !
45  integer mass,stiffness,buckling,rhsi
46 !
47  character*8 lakonl,lakon(*)
48  character*20 sideload(*)
49  character*80 matname(*),amat
50 !
51  integer konl(26),ifaceq(8,6),nelemload(2,*),nk,nelem,nmethod,
52  & mattyp,ithermal(2),iperturb,nload,idist,i,j,k,i1,j1,l,m,
53  & ii,jj,id,ipointer,ig,kk,istiff,iperm(20),ipompc(*),mi(*),
54  & nrhcon(*),ielmat(mi(3),*),ielorien(mi(3),*),nodempc(3,*),nmpc,
55  & ntmat_,nope,nopes,norien,iexpl,imat,mint2d,ikmpc(*),iloc,
56  & mint3d,ifacet(6,4),nopev,iorien,ilmpc(*),kode,jfaces,null,
57  & ifacew(8,5),intscheme,ipointeri,ipointerj,ncocon(2,*),
58  & nshcon(*),iinc,istep,jltyp,nfield,node,iflag,iscale,ielprop(*),
59  & nplkcon(0:ntmat_,*),nelcon(2,*),npmat_,ncmat_,i2,ipkon(*),
60  & iemchange,kon(*),mortar,nplicon(0:ntmat_,*),indexe,igauss,
61  & iponoel(*),inoel(2,*),nstate_,network,ipobody(2,*),ibody(3,*)
62 !
63  real*8 co(3,*),xl(3,26),shp(4,26),xstiff(27,mi(1),*),
64  & s(60,60),w(3,3),ff(60),shpj(4,26),sinktemp,xs2(3,7),
65  & rhcon(0:1,ntmat_,*),dxsj2,temp,press,xloadold(2,*),
66  & orab(7,*),t0(*),t1(*),coords(3),c1,c2,reltime,prop(*),
67  & xl2(3,9),xsj2(3),shp2(7,9),vold(0:mi(2),*),xload(2,*),
68  & xi,et,ze,xsj,xsjj,sm(60,60),t1l,rho,summass,summ,ttime,time,
69  & sume,factorm,factore,alp,weight,pgauss(3),timeend(2),
70  & cocon(0:6,ntmat_,*),shcon(0:3,ntmat_,*),sph,coconloc(6),
71  & field,areaj,sax(60,60),ffax(60),coefmpc(*),tl2(8),
72  & voldl(0:mi(2),26),springarea(2,*),plkcon(0:2*npmat_,ntmat_,*),
73  & elcon(0:ncmat_,ntmat_,*),elconloc(21),pslavsurf(3,*),
74  & pmastsurf(2,*),clearini(3,9,*),plicon(0:2*npmat_,ntmat_,*),
75  & sti(6,mi(1),*),xstate(nstate_,mi(1),*),xbody(7,*),
76  & xstateini(nstate_,mi(1),*),heatnod,heatfac
77 !
78  real*8 dtime,physcon(*)
79 !
80  intent(in) co,nk,kon,lakonl,
81  & nelem,rhcon,nrhcon,ielmat,ielorien,norien,orab,
82  & ntmat_,t0,t1,ithermal,vold,iperturb,nelemload,
83  & sideload,nload,idist,iexpl,dtime,
84  & matname,mi,mass,stiffness,buckling,rhsi,intscheme,
85  & physcon,shcon,nshcon,cocon,ncocon,ttime,time,istep,iinc,
86  & xstiff,xloadold,reltime,
87  & ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,
88  & plkcon,nplkcon,npmat_,ncmat_,elcon,nelcon,lakon,
89  & pslavsurf,pmastsurf,mortar,clearini,plicon,nplicon,ipkon,
90  & ielprop,prop
91 !
92  intent(inout) s,sm,ff,xload,nmethod,springarea
93 !
94  include "gauss.f"
95 !
96  ifaceq=reshape((/4,3,2,1,11,10,9,12,
97  & 5,6,7,8,13,14,15,16,
98  & 1,2,6,5,9,18,13,17,
99  & 2,3,7,6,10,19,14,18,
100  & 3,4,8,7,11,20,15,19,
101  & 4,1,5,8,12,17,16,20/),(/8,6/))
102  ifacet=reshape((/1,3,2,7,6,5,
103  & 1,2,4,5,9,8,
104  & 2,3,4,6,10,9,
105  & 1,4,3,8,10,7/),(/6,4/))
106  ifacew=reshape((/1,3,2,9,8,7,0,0,
107  & 4,5,6,10,11,12,0,0,
108  & 1,2,5,4,7,14,10,13,
109  & 2,3,6,5,8,15,11,14,
110  & 4,6,3,1,12,15,9,13/),(/8,5/))
111  null=0
112  iperm=(/5,6,7,8,1,2,3,4,13,14,15,16,9,10,11,12,17,18,19,20/)
113 !
114  iflag=3
115 !
116  timeend(1)=time
117  timeend(2)=ttime+time
118 !
119  summass=0.d0
120  heatfac=0.d0
121 !
122  indexe=ipkon(nelem)
123  imat=ielmat(1,nelem)
124  amat=matname(imat)
125  if(norien.gt.0) then
126  iorien=ielorien(1,nelem)
127  else
128  iorien=0
129  endif
130 !
131  if(lakonl(4:5).eq.'20') then
132  nope=20
133  nopev=8
134  nopes=8
135  elseif(lakonl(4:4).eq.'8') then
136  nope=8
137  nopev=8
138  nopes=4
139  elseif(lakonl(4:5).eq.'10') then
140  nope=10
141  nopev=4
142  nopes=6
143  elseif(lakonl(4:4).eq.'4') then
144  nope=4
145  nopev=4
146  nopes=3
147  elseif(lakonl(4:5).eq.'15') then
148  nope=15
149  nopev=6
150  elseif(lakonl(4:4).eq.'6') then
151  nope=6
152  nopev=6
153  elseif(lakonl(1:2).eq.'ES') then
154  if(lakonl(7:7).eq.'C') then
155  if(mortar.eq.0) then
156 c read(lakonl(8:8),'(i1)') nope
157  nope=ichar(lakonl(8:8))-47
158 c nope=nope+1
159  konl(nope+1)=kon(indexe+nope+1)
160  elseif(mortar.eq.1) then
161  nope=kon(indexe)
162  endif
163  else
164 c read(lakonl(8:8),'(i1)') nope
165  nope=ichar(lakonl(8:8))-47
166 c nope=nope+1
167  endif
168  elseif((lakonl(1:2).eq.'D ').or.
169  & ((lakonl(1:1).eq.'D').and.(network.eq.1))) then
170  nope=3
171  endif
172 !
173  if(intscheme.eq.0) then
174 !
175 ! # of 2D and 3D integration points
176 !
177  if(lakonl(4:5).eq.'8R') then
178  mint2d=1
179  mint3d=1
180  elseif(lakonl(4:7).eq.'20RB') then
181  if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
182  mint2d=4
183  mint3d=50
184  else
185  mint2d=4
186  call beamintscheme(lakonl,mint3d,ielprop(nelem),prop,
187  & null,xi,et,ze,weight)
188  endif
189  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) then
190  if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'E')) then
191  mint2d=2
192  mint3d=4
193  else
194  mint2d=4
195  mint3d=8
196  endif
197  elseif(lakonl(4:4).eq.'2') then
198  mint2d=9
199  mint3d=27
200  elseif(lakonl(4:5).eq.'10') then
201  mint2d=3
202  mint3d=4
203  elseif(lakonl(4:4).eq.'4') then
204  mint2d=1
205  mint3d=1
206  elseif(lakonl(4:5).eq.'15') then
207  mint3d=9
208  elseif(lakonl(4:4).eq.'6') then
209  mint3d=2
210  else
211  mint3d=0
212  endif
213  else
214 !
215 ! # of 3D integration points
216 !
217  if((lakonl(4:4).eq.'8').or.(lakonl(4:4).eq.'2')) then
218  mint3d=27
219  elseif((lakonl(4:5).eq.'10').or.(lakonl(4:4).eq.'4')) then
220  mint3d=15
221  elseif((lakonl(4:5).eq.'15').or.(lakonl(4:4).eq.'6')) then
222  mint3d=9
223  else
224  mint3d=0
225  endif
226 !
227 ! # of 2D integration points
228 !
229  if(lakonl(4:5).eq.'8R') then
230  mint2d=1
231  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) then
232  if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'E')) then
233  mint2d=2
234  else
235  mint2d=4
236  endif
237  elseif(lakonl(4:4).eq.'2') then
238  mint2d=9
239  elseif(lakonl(4:5).eq.'10') then
240  mint2d=3
241  elseif(lakonl(4:4).eq.'4') then
242  mint2d=1
243  endif
244  endif
245 !
246 ! computation of the coordinates of the local nodes
247 !
248  do i=1,nope
249  konl(i)=kon(indexe+i)
250  do j=1,3
251  xl(j,i)=co(j,konl(i))
252  enddo
253  enddo
254 !
255 ! initialisation for distributed forces
256 !
257  if(rhsi.eq.1) then
258  if(idist.ne.0) then
259  do i=1,nope
260  ff(i)=0.d0
261  enddo
262  endif
263  endif
264 !
265 ! initialisation of sm
266 !
267  if((mass.eq.1).or.(buckling.eq.1)) then
268  do i=1,nope
269  do j=i,nope
270  sm(i,j)=0.d0
271  enddo
272  enddo
273  endif
274 !
275 ! initialisation of s
276 !
277  do i=1,nope
278  do j=i,nope
279  s(i,j)=0.d0
280  enddo
281  enddo
282 !
283 ! calculating the stiffness matrix for the contact spring elements
284 !
285  if(mint3d.eq.0) then
286  do i1=1,nope
287  do i2=0,3
288  voldl(i2,i1)=vold(i2,konl(i1))
289  enddo
290  enddo
291 !
292  if(lakonl(1:2).eq.'ES') then
293  if(lakonl(7:7).eq.'C') then
294 !
295 ! contact element
296 !
297  kode=nelcon(1,imat)
298  if(kode.eq.-51) then
299  if(mortar.eq.0) then
300  call springstiff_n2f_th(xl,voldl,s,imat,elcon,
301  & nelcon,ncmat_,ntmat_,nope,kode,plkcon,nplkcon,
302  & npmat_,iperturb,springarea(1,konl(nope+1)),mi,
303  & timeend,matname,konl(nope),nelem,istep,iinc)
304  elseif(mortar.eq.1) then
305  iloc=kon(indexe+nope+1)
306  jfaces=kon(indexe+nope+2)
307  igauss=kon(indexe+nope+1)
308  node=0
309  call springstiff_f2f_th(xl,voldl,s,imat,elcon,
310  & nelcon,ncmat_,ntmat_,nope,lakonl,kode,
311  & elconloc,plicon,nplicon,npmat_,
312  & springarea(1,iloc),
313  & nmethod,mi,reltime,jfaces,igauss,pslavsurf,
314  & pmastsurf,clearini,matname,plkcon,nplkcon,
315  & node,nelem,istep,iinc,timeend)
316  endif
317  endif
318  elseif(lakonl(7:7).eq.'F') then
319 !
320 ! advective element
321 !
322  call advecstiff(nope,voldl,ithermal,xl,nelemload,
323  & nelem,nload,lakon,xload,istep,time,ttime,dtime,
324  & sideload,vold,mi,xloadold,reltime,nmethod,s,
325  & iinc,iponoel,inoel,ielprop,prop,ielmat,shcon,
326  & nshcon,rhcon,nrhcon,ntmat_,ipkon,kon,cocon,ncocon,
327  & ipobody,xbody,ibody)
328  endif
329  elseif((lakonl(1:2).eq.'D ').or.
330  & ((lakonl(1:1).eq.'D').and.(network.eq.1))) then
331  do i=1,3
332  do j=1,i
333  s(i,j)=0.d0
334  enddo
335  enddo
336  call networkstiff(voldl,s,imat,konl,mi,ntmat_,shcon,
337  & nshcon,rhcon,nrhcon)
338  endif
339  return
340  endif
341 !
342 ! computation of the matrix: loop over the Gauss points
343 !
344  do kk=1,mint3d
345  if(intscheme.eq.0) then
346  if(lakonl(4:5).eq.'8R') then
347  xi=gauss3d1(1,kk)
348  et=gauss3d1(2,kk)
349  ze=gauss3d1(3,kk)
350  weight=weight3d1(kk)
351  elseif(lakonl(4:7).eq.'20RB') then
352  if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
353  xi=gauss3d13(1,kk)
354  et=gauss3d13(2,kk)
355  ze=gauss3d13(3,kk)
356  weight=weight3d13(kk)
357  else
358  call beamintscheme(lakonl,mint3d,ielprop(nelem),prop,
359  & kk,xi,et,ze,weight)
360  endif
361  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R'))
362  & then
363  xi=gauss3d2(1,kk)
364  et=gauss3d2(2,kk)
365  ze=gauss3d2(3,kk)
366  weight=weight3d2(kk)
367  elseif(lakonl(4:4).eq.'2') then
368  xi=gauss3d3(1,kk)
369  et=gauss3d3(2,kk)
370  ze=gauss3d3(3,kk)
371  weight=weight3d3(kk)
372  elseif(lakonl(4:5).eq.'10') then
373  xi=gauss3d5(1,kk)
374  et=gauss3d5(2,kk)
375  ze=gauss3d5(3,kk)
376  weight=weight3d5(kk)
377  elseif(lakonl(4:4).eq.'4') then
378  xi=gauss3d4(1,kk)
379  et=gauss3d4(2,kk)
380  ze=gauss3d4(3,kk)
381  weight=weight3d4(kk)
382  elseif(lakonl(4:5).eq.'15') then
383  xi=gauss3d8(1,kk)
384  et=gauss3d8(2,kk)
385  ze=gauss3d8(3,kk)
386  weight=weight3d8(kk)
387  else
388  xi=gauss3d7(1,kk)
389  et=gauss3d7(2,kk)
390  ze=gauss3d7(3,kk)
391  weight=weight3d7(kk)
392  endif
393  else
394  if((lakonl(4:4).eq.'8').or.(lakonl(4:4).eq.'2')) then
395  xi=gauss3d3(1,kk)
396  et=gauss3d3(2,kk)
397  ze=gauss3d3(3,kk)
398  weight=weight3d3(kk)
399  elseif((lakonl(4:5).eq.'10').or.(lakonl(4:4).eq.'4')) then
400  xi=gauss3d6(1,kk)
401  et=gauss3d6(2,kk)
402  ze=gauss3d6(3,kk)
403  weight=weight3d6(kk)
404  else
405  xi=gauss3d8(1,kk)
406  et=gauss3d8(2,kk)
407  ze=gauss3d8(3,kk)
408  weight=weight3d8(kk)
409  endif
410  endif
411 !
412 ! calculation of the shape functions and their derivatives
413 ! in the gauss point
414 !
415  if(nope.eq.20) then
416  if(lakonl(7:7).eq.'A') then
417  call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag)
418  elseif((lakonl(7:7).eq.'E').or.(lakonl(7:7).eq.'S')) then
419  call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag)
420  else
421  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
422  endif
423  elseif(nope.eq.8) then
424  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
425  elseif(nope.eq.10) then
426  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
427  elseif(nope.eq.4) then
428  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
429  elseif(nope.eq.15) then
430  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
431  else
432  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
433  endif
434 !
435 ! check the jacobian determinant
436 !
437  if(xsj.lt.1.d-20) then
438  write(*,*) '*ERROR in e_c3d_th: nonpositive jacobian'
439  write(*,*) ' determinant in element',nelem
440  write(*,*)
441  xsj=dabs(xsj)
442  nmethod=0
443  endif
444 !
445 ! calculating the temperature in the integration
446 ! point
447 !
448  t1l=0.d0
449 !
450  if((lakonl(4:7).eq.'20RA').or.(lakonl(4:7).eq.'20RS').or.
451  & (lakonl(4:7).eq.'20RE')) then
452  t1l=vold(0,konl(1))*(shp(4,1)+shp(4,5)+shp(4,17))
453  & +vold(0,konl(2))*(shp(4,2)+shp(4,6)+shp(4,18))
454  & +vold(0,konl(3))*(shp(4,3)+shp(4,7)+shp(4,19))
455  & +vold(0,konl(4))*(shp(4,4)+shp(4,8)+shp(4,20))
456  & +vold(0,konl(9))*(shp(4,9)+shp(4,13))
457  & +vold(0,konl(10))*(shp(4,10)+shp(4,14))
458  & +vold(0,konl(11))*(shp(4,11)+shp(4,15))
459  & +vold(0,konl(12))*(shp(4,12)+shp(4,16))
460  else
461  do i1=1,nope
462  t1l=t1l+shp(4,i1)*vold(0,konl(i1))
463  enddo
464  endif
465 !
466 ! calculating the coordinates of the integration point
467 ! for material orientation purposes (for cylindrical
468 ! coordinate systems)
469 !
470  if(iorien.gt.0) then
471  do j=1,3
472  pgauss(j)=0.d0
473  do i1=1,nope
474  pgauss(j)=pgauss(j)+shp(4,i1)*xl(j,i1)
475  enddo
476  enddo
477  endif
478 !
479 ! material data
480 !
481  istiff=1
482  call materialdata_th(cocon,ncocon,imat,iorien,pgauss,orab,
483  & ntmat_,coconloc,mattyp,t1l,rhcon,nrhcon,rho,shcon,
484  & nshcon,sph,xstiff,kk,nelem,istiff,mi(1))
485 !
486 ! incorporating the jacobian determinant in the shape
487 ! functions
488 !
489  xsjj=dsqrt(xsj)
490  do i1=1,nope
491  shpj(1,i1)=shp(1,i1)*xsjj
492  shpj(2,i1)=shp(2,i1)*xsjj
493  shpj(3,i1)=shp(3,i1)*xsjj
494  shpj(4,i1)=shp(4,i1)*xsj
495  enddo
496 !
497  c1=coconloc(1)*weight
498  c2=rho*sph*weight
499 !
500 ! determination of the stiffness, and/or mass and/or
501 ! buckling matrix
502 !
503  do jj=1,nope
504 !
505  do ii=1,jj
506 !
507 ! the following section calculates the static
508 ! part of the stiffness matrix which, for buckling
509 ! calculations, is done in a preliminary static
510 ! call
511 !
512  if(mattyp.eq.1) then
513 !
514  s(ii,jj)=s(ii,jj)+c1*
515  & (shpj(1,ii)*shpj(1,jj)+shpj(2,ii)*shpj(2,jj)
516  & +shpj(3,ii)*shpj(3,jj))
517 !
518  elseif(mattyp.eq.2) then
519 !
520  s(ii,jj)=s(ii,jj)+(coconloc(1)*shpj(1,ii)*shpj(1,jj)
521  & +coconloc(2)*shpj(2,ii)*shpj(2,jj)
522  & +coconloc(3)*shpj(3,ii)*shpj(3,jj))*weight
523 !
524  else
525 !
526  do i1=1,3
527  do j1=1,3
528  w(i1,j1)=shpj(i1,ii)*shpj(j1,jj)
529  enddo
530  enddo
531 !
532  s(ii,jj)=s(ii,jj)+
533  & (coconloc(1)*w(1,1)+
534  & coconloc(4)*(w(1,2)+w(2,1))+
535  & coconloc(2)*w(2,2)+
536  & coconloc(5)*(w(1,3)+w(3,1))+
537  & coconloc(6)*(w(2,3)+w(3,2))+
538  & coconloc(3)*w(3,3))*weight
539 !
540  endif
541 !
542 ! mass matrix
543 !
544  if(mass.eq.1) then
545  sm(ii,jj)=sm(ii,jj)
546  & +c2*shpj(4,ii)*shp(4,jj)
547  endif
548 !
549  enddo
550  enddo
551 !
552 ! computation of the right hand side
553 !
554  if(rhsi.eq.1) then
555 !
556 ! distributed heat flux
557 !
558  if(nload.gt.0) then
559  call nident2(nelemload,nelem,nload,id)
560  areaj=xsj*weight
561  do
562  if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
563  if(sideload(id)(1:2).ne.'BF') then
564  id=id-1
565  cycle
566  endif
567  if(sideload(id)(3:4).eq.'NU') then
568  if(ithermal(2).eq.2) then
569  do j=1,3
570  pgauss(j)=0.d0
571  do i1=1,nope
572  pgauss(j)=pgauss(j)+
573  & shp(4,i1)*xl(j,i1)
574  enddo
575  enddo
576  else
577  do j=1,3
578  pgauss(j)=0.d0
579  do i1=1,nope
580  pgauss(j)=pgauss(j)+
581  & shp(4,i1)*(xl(j,i1)+vold(j,konl(i1)))
582  enddo
583  enddo
584  endif
585  jltyp=1
586  iscale=1
587  call dflux(xload(1,id),t1l,istep,iinc,timeend,
588  & nelem,kk,pgauss,jltyp,temp,press,sideload(id),
589  & areaj,vold,co,lakonl,konl,ipompc,nodempc,coefmpc,
590  & nmpc,ikmpc,ilmpc,iscale,mi,sti,xstateini,xstate,
591  & nstate_,dtime)
592  if((nmethod.eq.1).and.(iscale.ne.0))
593  & xload(1,id)=xloadold(1,id)+
594  & (xload(1,id)-xloadold(1,id))*reltime
595  endif
596  do jj=1,nope
597  ff(jj)=ff(jj)+xload(1,id)*shpj(4,jj)*weight
598  enddo
599  exit
600  enddo
601  endif
602  endif
603 !
604  enddo
605 !
606  if((buckling.eq.0).and.(nload.ne.0)) then
607  iflag=2
608 !
609 ! distributed loads
610 !
611  call nident2(nelemload,nelem,nload,id)
612  do
613  if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
614  if((sideload(id)(1:1).ne.'F').and.
615  & (sideload(id)(1:1).ne.'R').and.
616  & (sideload(id)(1:1).ne.'S')) then
617  id=id-1
618  cycle
619  endif
620 c read(sideload(id)(2:2),'(i1)') ig
621  ig=ichar(sideload(id)(2:2))-48
622 !
623 !
624 ! treatment of wedge faces
625 !
626  if(lakonl(4:4).eq.'6') then
627  mint2d=1
628  if(ig.le.2) then
629  nopes=3
630  else
631  nopes=4
632  endif
633  endif
634  if(lakonl(4:5).eq.'15') then
635  if(ig.le.2) then
636  mint2d=3
637  nopes=6
638  else
639  mint2d=4
640  nopes=8
641  endif
642  endif
643 !
644  if((nope.eq.20).or.(nope.eq.8)) then
645  do i=1,nopes
646  tl2(i)=vold(0,konl(ifaceq(i,ig)))
647  if(ithermal(2).eq.2) then
648  do j=1,3
649  xl2(j,i)=co(j,konl(ifaceq(i,ig)))
650  enddo
651  else
652  do j=1,3
653  xl2(j,i)=co(j,konl(ifaceq(i,ig)))+
654  & vold(j,konl(ifaceq(i,ig)))
655  enddo
656  endif
657  enddo
658  elseif((nope.eq.10).or.(nope.eq.4)) then
659  do i=1,nopes
660  tl2(i)=vold(0,konl(ifacet(i,ig)))
661  if(ithermal(2).eq.2) then
662  do j=1,3
663  xl2(j,i)=co(j,konl(ifacet(i,ig)))
664  enddo
665  else
666  do j=1,3
667  xl2(j,i)=co(j,konl(ifacet(i,ig)))+
668  & vold(j,konl(ifacet(i,ig)))
669  enddo
670  endif
671  enddo
672  else
673  do i=1,nopes
674  tl2(i)=vold(0,konl(ifacew(i,ig)))
675  if(ithermal(2).eq.2) then
676  do j=1,3
677  xl2(j,i)=co(j,konl(ifacew(i,ig)))
678  enddo
679  else
680  do j=1,3
681  xl2(j,i)=co(j,konl(ifacew(i,ig)))+
682  & vold(j,konl(ifacew(i,ig)))
683  enddo
684  endif
685  enddo
686  endif
687 !
688 ! storing envnode values into xload for non-cavity
689 ! radiation
690 !
691  if((sideload(id)(1:1).eq.'R').and.
692  & (sideload(id)(3:4).ne.'CR').and.
693  & (nelemload(2,id).gt.0)) then
694  xload(2,id)=vold(0,nelemload(2,id))-physcon(1)
695  endif
696 !
697  do i=1,mint2d
698 !
699 ! copying the sink temperature to ensure the same
700 ! value in each integration point (sinktemp can be
701 ! changed in subroutine film: requirement from the
702 ! thermal people)
703 !
704  sinktemp=xload(2,id)
705 !
706  if((lakonl(4:5).eq.'8R').or.
707  & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
708  xi=gauss2d1(1,i)
709  et=gauss2d1(2,i)
710  weight=weight2d1(i)
711  elseif((lakonl(4:4).eq.'8').or.
712  & (lakonl(4:6).eq.'20R').or.
713  & ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
714  xi=gauss2d2(1,i)
715  et=gauss2d2(2,i)
716  weight=weight2d2(i)
717  elseif(lakonl(4:4).eq.'2') then
718  xi=gauss2d3(1,i)
719  et=gauss2d3(2,i)
720  weight=weight2d3(i)
721  elseif((lakonl(4:5).eq.'10').or.
722  & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
723  xi=gauss2d5(1,i)
724  et=gauss2d5(2,i)
725  weight=weight2d5(i)
726  elseif((lakonl(4:4).eq.'4').or.
727  & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
728  xi=gauss2d4(1,i)
729  et=gauss2d4(2,i)
730  weight=weight2d4(i)
731  endif
732 !
733  if(nopes.eq.9) then
734  call shape9q(xi,et,xl2,xsj2,xs2,shp2,iflag)
735  elseif(nopes.eq.8) then
736  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
737  elseif(nopes.eq.4) then
738  call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
739  elseif(nopes.eq.6) then
740  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
741  elseif(nopes.eq.7) then
742  call shape7tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
743  else
744  call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
745  endif
746 !
747  dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
748  & xsj2(3)*xsj2(3))
749  areaj=dxsj2*weight
750 !
751  temp=0.d0
752  do j=1,nopes
753  temp=temp+tl2(j)*shp2(4,j)
754  enddo
755 !
756 ! for nonuniform load: determine the coordinates of the
757 ! point (transferred into the user subroutine)
758 !
759  if((sideload(id)(3:4).eq.'NU').or.
760  & (sideload(id)(5:6).eq.'NU')) then
761  do k=1,3
762  coords(k)=0.d0
763  do j=1,nopes
764  coords(k)=coords(k)+xl2(k,j)*shp2(4,j)
765  enddo
766  enddo
767 c read(sideload(id)(2:2),'(i1)') jltyp
768  jltyp=ichar(sideload(id)(2:2))-48
769  jltyp=jltyp+10
770  if(sideload(id)(1:1).eq.'S') then
771  iscale=1
772  call dflux(xload(1,id),temp,istep,iinc,timeend,
773  & nelem,i,coords,jltyp,temp,press,sideload(id),
774  & areaj,vold,co,lakonl,konl,
775  & ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,iscale,mi,
776  & sti,xstateini,xstate,nstate_,dtime)
777 c if((nmethod.eq.1).and.(iscale.ne.0))
778 c & xload(1,id)=xloadold(1,id)+
779 c & (xload(1,id)-xloadold(1,id))*reltime
780  elseif(sideload(id)(1:1).eq.'F') then
781  node=nelemload(2,id)
782  call film(xload(1,id),sinktemp,temp,istep,
783  & iinc,timeend,nelem,i,coords,jltyp,field,nfield,
784  & sideload(id),node,areaj,vold,mi,
785  & ipkon,kon,lakon,iponoel,inoel,ielprop,prop,ielmat,
786  & shcon,nshcon,rhcon,nrhcon,ntmat_,cocon,ncocon,
787  & ipobody,xbody,ibody,heatnod,heatfac)
788 c if(nmethod.eq.1) xload(1,id)=xloadold(1,id)+
789 c & (xload(1,id)-xloadold(1,id))*reltime
790  elseif(sideload(id)(1:1).eq.'R') then
791  call radiate(xload(1,id),xload(2,id),temp,istep,
792  & iinc,timeend,nelem,i,coords,jltyp,field,nfield,
793  & sideload(id),node,areaj,vold,mi,iemchange)
794 c if(nmethod.eq.1) xload(1,id)=xloadold(1,id)+
795 c & (xload(1,id)-xloadold(1,id))*reltime
796  endif
797  endif
798 
799 !
800  do k=1,nopes
801  if((nope.eq.20).or.(nope.eq.8)) then
802  ipointer=ifaceq(k,ig)
803  elseif((nope.eq.10).or.(nope.eq.4)) then
804  ipointer=ifacet(k,ig)
805  else
806  ipointer=ifacew(k,ig)
807  endif
808  if(sideload(id)(1:1).eq.'S') then
809 !
810 ! flux INTO the face is positive (input deck convention)
811 ! this is different from the convention in the theory
812 !
813  ff(ipointer)=ff(ipointer)+shp2(4,k)*xload(1,id)
814  & *areaj
815  elseif(sideload(id)(1:1).eq.'F') then
816  ff(ipointer)=ff(ipointer)+shp2(4,k)*areaj*
817  & (heatfac-xload(1,id)*(temp-sinktemp))
818  elseif(sideload(id)(1:1).eq.'R') then
819  ff(ipointer)=ff(ipointer)-shp2(4,k)*physcon(2)*
820  & xload(1,id)*((temp-physcon(1))**4-
821  & (xload(2,id)-physcon(1))**4)*
822  & areaj
823  endif
824  enddo
825 !
826  do ii=1,nopes
827  if((nope.eq.20).or.(nope.eq.8)) then
828  ipointeri=ifaceq(ii,ig)
829  elseif((nope.eq.10).or.(nope.eq.4)) then
830  ipointeri=ifacet(ii,ig)
831  else
832  ipointeri=ifacew(ii,ig)
833  endif
834  do jj=1,nopes
835  if((nope.eq.20).or.(nope.eq.8)) then
836  ipointerj=ifaceq(jj,ig)
837  elseif((nope.eq.10).or.(nope.eq.4)) then
838  ipointerj=ifacet(jj,ig)
839  else
840  ipointerj=ifacew(jj,ig)
841  endif
842  if(ipointeri.gt.ipointerj) cycle
843  if(sideload(id)(1:1).eq.'F') then
844  s(ipointeri,ipointerj)=s(ipointeri,ipointerj)+
845  & shp2(4,ii)*shp2(4,jj)*xload(1,id)*areaj
846  elseif(sideload(id)(1:1).eq.'R') then
847  s(ipointeri,ipointerj)=s(ipointeri,ipointerj)+
848  & shp2(4,ii)*shp2(4,jj)*physcon(2)*xload(1,id)*
849  & areaj*4.d0*(temp-physcon(1))**3
850  endif
851  enddo
852  enddo
853 !
854  enddo
855 !
856  id=id-1
857  enddo
858  endif
859 !
860 ! for axially symmetric and plane stress/strain elements:
861 ! complete s and sm
862 !
863  if(((lakonl(4:5).eq.'8 ').or.
864  & ((lakonl(4:6).eq.'20R').and.(lakonl(7:8).ne.'BR'))).and.
865  & ((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'E'))) then
866  do i=1,20
867  do j=i,20
868  k=iperm(i)
869  l=iperm(j)
870  if(k.gt.l) then
871  m=k
872  k=l
873  l=m
874  endif
875  sax(i,j)=s(k,l)
876  enddo
877  enddo
878  do i=1,20
879  do j=i,20
880  s(i,j)=s(i,j)+sax(i,j)
881  enddo
882  enddo
883 !
884 ! special treatment of plane stress elements since lateral
885 ! heating is allowed (orthogonal to the plane)
886 !
887  if(nload.ne.0) then
888  do i=1,20
889  k=iperm(i)
890  ffax(i)=ff(k)
891  enddo
892  do i=1,20
893  ff(i)=ff(i)+ffax(i)
894  enddo
895  endif
896 !
897  if(mass.eq.1) then
898  do i=1,20
899  do j=i,20
900  k=iperm(i)
901  l=iperm(j)
902  if(k.gt.l) then
903  m=k
904  k=l
905  l=m
906  endif
907  sax(i,j)=sm(k,l)
908  enddo
909  enddo
910  do i=1,20
911  do j=i,20
912  sm(i,j)=sm(i,j)+sax(i,j)
913  enddo
914  enddo
915  endif
916  endif
917 !
918  if((mass.eq.1).and.(iexpl.gt.1)) then
919 !
920 ! scaling the diagonal terms of the mass matrix such that the total
921 ! mass is right (LUMPING; for explicit dynamic calculations)
922 !
923  sume=0.d0
924  summ=0.d0
925  do i=1,nopev
926  sume=sume+sm(i,i)
927  enddo
928  do i=nopev+1,nope
929  summ=summ+sm(i,i)
930  enddo
931 !
932  if(nope.eq.20) then
933  alp=.2917d0
934  elseif(nope.eq.10) then
935  alp=0.1203d0
936  elseif(nope.eq.15) then
937  alp=0.2141d0
938  endif
939 !
940  if((nope.eq.20).or.(nope.eq.10).or.
941  & (nope.eq.15)) then
942  factore=summass*alp/(1.d0+alp)/sume
943  factorm=summass/(1.d0+alp)/summ
944  else
945  factore=summass/sume
946  endif
947 !
948  do i=1,nopev
949  sm(i,i)=sm(i,i)*factore
950  enddo
951  do i=nopev+1,nope
952  sm(i,i)=sm(i,i)*factorm
953  enddo
954 !
955  endif
956 !
957  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
shape9q
subroutine shape9q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape9q.f:20
shape8q
subroutine shape8q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape8q.f:20
springstiff_n2f_th
subroutine springstiff_n2f_th(xl, voldl, s, imat, elcon, nelcon, ncmat_, ntmat_, nope, kode, plkcon, nplkcon, npmat_, iperturb, springarea, mi, timeend, matname, node, noel, istep, iinc)
Definition: springstiff_n2f_th.f:23
shape3tri
subroutine shape3tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape3tri.f:20
shape7tri
subroutine shape7tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape7tri.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
c1
static double * c1
Definition: mafillvcompmain.c:30
shape20h_ax
subroutine shape20h_ax(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape20h_ax.f:20
springstiff_f2f_th
subroutine springstiff_f2f_th(xl, voldl, s, imat, elcon, nelcon, ncmat_, ntmat_, nope, lakonl, kode, elconloc, plicon, nplicon, npmat_, springarea, nmethod, mi, reltime, jfaces, igauss, pslavsurf, pmastsurf, clearini, matname, plkcon, nplkcon, node, noel, istep, iinc, timeend)
Definition: springstiff_f2f_th.f:24
idist
static ITG * idist
Definition: radflowload.c:39
materialdata_th
subroutine materialdata_th(cocon, ncocon, imat, iorien, pgauss, orab, ntmat_, coconloc, mattyp, t1l, rhcon, nrhcon, rho, shcon, nshcon, sph, xstiff, iint, iel, istiff, mi)
Definition: materialdata_th.f:21
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
film
subroutine film(h, sink, temp, kstep, kinc, time, noel, npt, coords, jltyp, field, nfield, loadtype, node, area, vold, mi, ipkon, kon, lakon, iponoel, inoel, ielprop, prop, ielmat, shcon, nshcon, rhcon, nrhcon, ntmat_, cocon, ncocon, ipobody, xbody, ibody, heatnod, heatfac)
Definition: film.f:24
advecstiff
subroutine advecstiff(nope, voldl, ithermal, xl, nelemload, nelemadvec, nload, lakon, xload, istep, time, ttime, dtime, sideload, vold, mi, xloadold, reltime, nmethod, s, iinc, iponoel, inoel, ielprop, prop, ielmat, shcon, nshcon, rhcon, nrhcon, ntmat_, ipkon, kon, cocon, ncocon, ipobody, xbody, ibody)
Definition: advecstiff.f:24
radiate
subroutine radiate(e, sink, temp, kstep, kinc, time, noel, npt, coords, jltyp, field, nfield, loadtype, node, area, vold, mi, iemchange)
Definition: radiate.f:22
networkstiff
subroutine networkstiff(voldl, s, imat, konl, mi, ntmat_, shcon, nshcon, rhcon, nrhcon)
Definition: networkstiff.f:21
Hosted by OpenAircraft.com, (Michigan UAV, LLC)