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

Go to the source code of this file.

Functions/Subroutines

subroutine e_c3d_em (co, konl, lakonl, s, sm, ff, nelem, nmethod, ielmat, ntmat_, t1, ithermal, vold, idist, matname, mi, mass, rhsi, ncmat_, elcon, nelcon, h0, iactive, alcon, nalcon, istartset, iendset, ialset)
 

Function/Subroutine Documentation

◆ e_c3d_em()

subroutine e_c3d_em ( real*8, dimension(3,*)  co,
integer, dimension(26)  konl,
character*8  lakonl,
real*8, dimension(100,100)  s,
real*8, dimension(100,100)  sm,
real*8, dimension(100)  ff,
integer  nelem,
integer  nmethod,
integer, dimension(mi(3),*)  ielmat,
integer  ntmat_,
real*8, dimension(*)  t1,
integer  ithermal,
real*8, dimension(0:mi(2),*)  vold,
integer  idist,
character*80, dimension(*)  matname,
integer, dimension(*)  mi,
logical  mass,
logical  rhsi,
integer  ncmat_,
real*8, dimension(0:ncmat_,ntmat_,*)  elcon,
integer, dimension(2,*)  nelcon,
real*8, dimension(3,*)  h0,
integer, dimension(3)  iactive,
real*8, dimension(0:6,ntmat_,*)  alcon,
integer, dimension(2,*)  nalcon,
integer, dimension(*)  istartset,
integer, dimension(*)  iendset,
integer, dimension(*)  ialset 
)
25 !
26 ! computation of the element matrix and rhs for the element with
27 ! the topology in konl
28 !
29 ! ff: rhs without temperature and eigenstress contribution
30 !
31 ! nmethod=0: check for positive Jacobian
32 ! nmethod=1: stiffness matrix + right hand side
33 ! nmethod=2: stiffness matrix + mass matrix
34 ! nmethod=3: static stiffness + buckling stiffness
35 ! nmethod=4: stiffness matrix + mass matrix
36 !
37  implicit none
38 !
39  logical mass,rhsi
40 !
41  character*8 lakonl
42  character*80 matname(*),amat
43 !
44  integer konl(26),ifaceq(8,6),nelem,nmethod,iactive(3),
45  & ithermal,idist,i,j,k,i1,m,one,ii,jj,id,ipointer,ig,kk,mi(*),
46  & ielmat(mi(3),*),ntmat_,nope,nopes,imat,mint2d,mint3d,
47  & ifacet(6,4),nopev,ifacew(8,5),ipointeri,ipointerj,iflag,
48  & nelcon(2,*),ncmat_,nalcon(2,*),iel,ii1,ilength,istart,iset,
49  & isurf,jj1,istartset(*),iendset(*),ialset(*),three,nfaces,
50  & null
51 !
52  real*8 co(3,*),xl(3,26),shp(4,26),s(100,100),ff(100),xs2(3,7),
53  & t1(*),h0(3,*),xl2(3,9),xsj2(3),shp2(7,9),vold(0:mi(2),*),
54  & xi,et,ze,xsj,sm(100,100),t1l,weight,
55  & elcon(0:ncmat_,ntmat_,*),elconloc(21),sigma,
56  & h0l(3),h0l2(3,8),alcon(0:6,ntmat_,*),diag,um,uminv,alpha(6),
57  & d(3,3)
58 !
59  include "gauss.f"
60 !
61  data ifaceq /4,3,2,1,11,10,9,12,
62  & 5,6,7,8,13,14,15,16,
63  & 1,2,6,5,9,18,13,17,
64  & 2,3,7,6,10,19,14,18,
65  & 3,4,8,7,11,20,15,19,
66  & 4,1,5,8,12,17,16,20/
67  data ifacet /1,3,2,7,6,5,
68  & 1,2,4,5,9,8,
69  & 2,3,4,6,10,9,
70  & 1,4,3,8,10,7/
71  data ifacew /1,3,2,9,8,7,0,0,
72  & 4,5,6,10,11,12,0,0,
73  & 1,2,5,4,7,14,10,13,
74  & 2,3,6,5,8,15,11,14,
75  & 4,6,3,1,12,15,9,13/
76  data d /1.d0,0.d0,0.d0,0.d0,1.d0,0.d0,0.d0,0.d0,1.d0/
77 !
78  null=0
79  one=1
80  three=3
81  iflag=3
82 !
83  imat=ielmat(1,nelem)
84  amat=matname(imat)
85 !
86  if(lakonl(4:5).eq.'20') then
87  nope=20
88  nopev=8
89  nopes=8
90  elseif(lakonl(4:4).eq.'8') then
91  nope=8
92  nopev=8
93  nopes=4
94  elseif(lakonl(4:5).eq.'10') then
95  nope=10
96  nopev=4
97  nopes=6
98  elseif(lakonl(4:4).eq.'4') then
99  nope=4
100  nopev=4
101  nopes=3
102  elseif(lakonl(4:5).eq.'15') then
103  nope=15
104  nopev=6
105  elseif(lakonl(4:4).eq.'6') then
106  nope=6
107  nopev=6
108  endif
109 !
110  if(lakonl(4:5).eq.'8R') then
111  mint2d=1
112  mint3d=1
113  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) then
114  mint2d=4
115  mint3d=8
116  elseif(lakonl(4:4).eq.'2') then
117  mint2d=9
118  mint3d=27
119  elseif(lakonl(4:5).eq.'10') then
120  mint2d=3
121  mint3d=4
122  elseif(lakonl(4:4).eq.'4') then
123  mint2d=1
124  mint3d=1
125  elseif(lakonl(4:5).eq.'15') then
126  mint3d=9
127  elseif(lakonl(4:4).eq.'6') then
128  mint3d=2
129  else
130  mint3d=0
131  endif
132 !
133 ! computation of the coordinates of the local nodes
134 !
135  do i=1,nope
136  do j=1,3
137  xl(j,i)=co(j,konl(i))
138  enddo
139  enddo
140 !
141 ! initialisation for distributed forces
142 !
143  if(rhsi) then
144  if(idist.ne.0) then
145  do i=1,5*nope
146  ff(i)=0.d0
147  enddo
148  endif
149  endif
150 !
151 ! initialisation of sm
152 !
153  if(mass) then
154  do i=1,5*nope
155  do j=i,5*nope
156  sm(i,j)=0.d0
157  enddo
158  enddo
159  endif
160 !
161 ! initialisation of s
162 !
163  do i=1,5*nope
164  do j=i,5*nope
165  s(i,j)=0.d0
166  enddo
167  enddo
168 !
169 ! computation of the matrix: loop over the Gauss points
170 !
171  do kk=1,mint3d
172  if(lakonl(4:5).eq.'8R') then
173  xi=gauss3d1(1,kk)
174  et=gauss3d1(2,kk)
175  ze=gauss3d1(3,kk)
176  weight=weight3d1(kk)
177  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R'))
178  & then
179  xi=gauss3d2(1,kk)
180  et=gauss3d2(2,kk)
181  ze=gauss3d2(3,kk)
182  weight=weight3d2(kk)
183  elseif(lakonl(4:4).eq.'2') then
184  xi=gauss3d3(1,kk)
185  et=gauss3d3(2,kk)
186  ze=gauss3d3(3,kk)
187  weight=weight3d3(kk)
188  elseif(lakonl(4:5).eq.'10') then
189  xi=gauss3d5(1,kk)
190  et=gauss3d5(2,kk)
191  ze=gauss3d5(3,kk)
192  weight=weight3d5(kk)
193  elseif(lakonl(4:4).eq.'4') then
194  xi=gauss3d4(1,kk)
195  et=gauss3d4(2,kk)
196  ze=gauss3d4(3,kk)
197  weight=weight3d4(kk)
198  elseif(lakonl(4:5).eq.'15') then
199  xi=gauss3d8(1,kk)
200  et=gauss3d8(2,kk)
201  ze=gauss3d8(3,kk)
202  weight=weight3d8(kk)
203  else
204  xi=gauss3d7(1,kk)
205  et=gauss3d7(2,kk)
206  ze=gauss3d7(3,kk)
207  weight=weight3d7(kk)
208  endif
209 !
210 ! calculation of the shape functions and their derivatives
211 ! in the gauss point
212 !
213  if(nope.eq.20) then
214  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
215  elseif(nope.eq.8) then
216  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
217  elseif(nope.eq.10) then
218  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
219  elseif(nope.eq.4) then
220  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
221  elseif(nope.eq.15) then
222  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
223  else
224  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
225  endif
226 !
227 ! check the jacobian determinant
228 !
229  if(xsj.lt.1.d-20) then
230  write(*,*) '*ERROR in e_c3d_th: nonpositive jacobian'
231  write(*,*) ' determinant in element',nelem
232  write(*,*)
233  xsj=dabs(xsj)
234  nmethod=0
235  endif
236 !
237 ! calculating the temperature and the magnetic intensity
238 ! in the integration point (one order lower than the
239 ! interpolation of the primary unknowns)
240 !
241  do j=1,3
242  h0l(j)=0.d0
243  enddo
244  t1l=0.d0
245 !
246 ! calculating the magnetic intensity
247 !
248  if(lakonl(4:5).eq.'8 ') then
249  do i1=1,nope
250  do j=1,3
251  h0l(j)=h0l(j)+h0(j,konl(i1))/8.d0
252  enddo
253  enddo
254  elseif(lakonl(4:6).eq.'20 ') then
255  call linvec(h0,konl,nope,kk,h0l,one,three)
256  elseif(lakonl(4:6).eq.'10T') then
257  call linvec(h0,konl,h0l,one,three,shp)
258  else
259  do i1=1,nope
260  do j=1,3
261  h0l(j)=h0l(j)+shp(4,i1)*h0(j,konl(i1))
262  enddo
263  enddo
264  endif
265 !
266 ! calculating the temperature
267 !
268  if(ithermal.eq.1) then
269  if(lakonl(4:5).eq.'8 ') then
270  do i1=1,nope
271  t1l=t1l+t1(konl(i1))/8.d0
272  enddo
273  elseif(lakonl(4:6).eq.'20 ')then
274  call linscal(t1,konl,nope,kk,t1l,one)
275  elseif(lakonl(4:6).eq.'10T') then
276  call linscal10(t1,konl,t1l,null,shp)
277  else
278  do i1=1,nope
279  t1l=t1l+shp(4,i1)*t1(konl(i1))
280  enddo
281  endif
282  elseif(ithermal.ge.2) then
283  if(lakonl(4:5).eq.'8 ') then
284  do i1=1,nope
285  t1l=t1l+vold(0,konl(i1))/8.d0
286  enddo
287  elseif(lakonl(4:6).eq.'20 ')then
288  call linscal(vold,konl,nope,kk,t1l,mi(2))
289  elseif(lakonl(4:6).eq.'10T') then
290  call linscal10(vold,konl,t1l,mi(2),shp)
291  else
292  do i1=1,nope
293  t1l=t1l+shp(4,i1)*vold(0,konl(i1))
294  enddo
295  endif
296  endif
297 !
298 ! material data (electric conductivity and
299 ! magnetic permeability)
300 !
301  call materialdata_em(elcon,nelcon,alcon,nalcon,
302  & imat,ntmat_,t1l,elconloc,ncmat_,alpha)
303 !
304  weight=weight*xsj
305 !
306  uminv=weight/elconloc(1)
307  um=weight*elconloc(1)
308  sigma=weight*alpha(1)
309 !
310  jj1=0
311  do jj=1,nope
312  ii1=0
313  do ii=1,jj
314 !
315  if((int(elconloc(2)).eq.2).or.(int(elconloc(2)).eq.3))
316  & then
317 !
318 ! K_AA matrix
319 !
320  diag=shp(1,ii)*shp(1,jj)+shp(2,ii)*shp(2,jj)+
321  & shp(3,ii)*shp(3,jj)
322  do k=1,3
323  do m=1,3
324  s(ii1+k,jj1+m)=s(ii1+k,jj1+m)+
325  & (diag*d(k,m)-shp(m,ii)*shp(k,jj)
326  & +shp(k,ii)*shp(m,jj))*uminv
327  enddo
328  enddo
329 !
330 ! M_AA matrix
331 !
332  if(mass) then
333  do k=1,3
334  sm(ii1+k,jj1+k)=sm(ii1+k,jj1+k)+
335  & sigma*shp(4,ii)*shp(4,jj)
336  enddo
337  endif
338  if((int(elconloc(2)).eq.2).and.mass) then
339 !
340 ! M_AV and M_VA matrix
341 !
342  do k=1,3
343  sm(ii1+k,jj1+4)=sm(ii1+k,jj1+4)+
344  & sigma*shp(4,ii)*shp(k,jj)
345  sm(ii1+4,jj1+k)=sm(ii1+4,jj1+k)+
346  & sigma*shp(k,ii)*shp(4,jj)
347  enddo
348 !
349 ! M_VV matrix
350 !
351  sm(ii1+4,jj1+4)=sm(ii1+4,jj1+4)+
352  & sigma*(shp(1,ii)*shp(1,jj)+
353  & shp(2,ii)*shp(2,jj)+
354  & shp(3,ii)*shp(3,jj))
355  endif
356  else if(int(elconloc(2)).eq.1) then
357 !
358 ! K_phiphi matrix
359 !
360  s(ii1+5,jj1+5)=s(ii1+5,jj1+5)-
361  & um*(shp(1,ii)*shp(1,jj)+
362  & shp(2,ii)*shp(2,jj)+
363  & shp(3,ii)*shp(3,jj))
364  endif
365 !
366  ii1=ii1+5
367  enddo
368 !
369 ! F_phi matrix
370 !
371  if((rhsi).and.(int(elconloc(2)).eq.1)) then
372  ff(jj1+5)=ff(jj1+5)-um*(shp(1,jj)*h0l(1)+
373  & shp(2,jj)*h0l(2)+
374  & shp(3,jj)*h0l(3))
375  endif
376 !
377  jj1=jj1+5
378  enddo
379  enddo
380 !
381 ! surface integrals
382 !
383 ! determining the number of faces per element
384 !
385  if((lakonl(4:4).eq.'8').or.(lakonl(4:4).eq.'2')) then
386  nfaces=6
387  elseif((lakonl(4:4).eq.'6').or.(lakonl(4:5).eq.'15')) then
388  nfaces=5
389  elseif((lakonl(4:4).eq.'4').or.(lakonl(4:5).eq.'10')) then
390  nfaces=4
391  endif
392 !
393  m=int(elconloc(2))
394  if(iactive(m).eq.0) return
395  iset=iactive(m)
396  istart=istartset(iset)
397  ilength=iendset(iset)-istart+1
398 !
399  isurf=10*nelem+nfaces
400  call nident(ialset(istart),isurf,ilength,id)
401 !
402  do
403  if(id.eq.0) exit
404  isurf=ialset(istart+id-1)
405  iel=int(isurf/10.d0)
406  if(iel.ne.nelem) exit
407  ig=isurf-10*iel
408 !
409 ! treatment of wedge faces
410 !
411  if(lakonl(4:4).eq.'6') then
412  mint2d=1
413  if(ig.le.2) then
414  nopes=3
415  else
416  nopes=4
417  endif
418  endif
419  if(lakonl(4:5).eq.'15') then
420  if(ig.le.2) then
421  mint2d=3
422  nopes=6
423  else
424  mint2d=4
425  nopes=8
426  endif
427  endif
428 !
429  if((nope.eq.20).or.(nope.eq.8)) then
430  do i=1,nopes
431  do j=1,3
432  h0l2(j,i)=h0(j,konl(ifaceq(i,ig)))
433  xl2(j,i)=co(j,konl(ifaceq(i,ig)))
434  enddo
435  enddo
436  elseif((nope.eq.10).or.(nope.eq.4)) then
437  do i=1,nopes
438  do j=1,3
439  h0l2(j,i)=h0(j,konl(ifacet(i,ig)))
440  xl2(j,i)=co(j,konl(ifacet(i,ig)))
441  enddo
442  enddo
443  else
444  do i=1,nopes
445  do j=1,3
446  h0l2(j,i)=h0(j,konl(ifacew(i,ig)))
447  xl2(j,i)=co(j,konl(ifacew(i,ig)))
448  enddo
449  enddo
450  endif
451 !
452  do i=1,mint2d
453 !
454  if((lakonl(4:5).eq.'8R').or.
455  & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
456  xi=gauss2d1(1,i)
457  et=gauss2d1(2,i)
458  weight=weight2d1(i)
459  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R').or.
460  & ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
461  xi=gauss2d2(1,i)
462  et=gauss2d2(2,i)
463  weight=weight2d2(i)
464  elseif(lakonl(4:4).eq.'2') then
465  xi=gauss2d3(1,i)
466  et=gauss2d3(2,i)
467  weight=weight2d3(i)
468  elseif((lakonl(4:5).eq.'10').or.
469  & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
470  xi=gauss2d5(1,i)
471  et=gauss2d5(2,i)
472  weight=weight2d5(i)
473  elseif((lakonl(4:4).eq.'4').or.
474  & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
475  xi=gauss2d4(1,i)
476  et=gauss2d4(2,i)
477  weight=weight2d4(i)
478  endif
479 !
480  if(nopes.eq.8) then
481  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
482  elseif(nopes.eq.4) then
483  call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
484  elseif(nopes.eq.6) then
485  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
486  else
487  call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
488  endif
489 !
490  do k=1,3
491  h0l(k)=0.d0
492  do j=1,nopes
493  h0l(k)=h0l(k)+h0l2(k,j)*shp2(4,j)
494  enddo
495  enddo
496 !
497  if((rhsi).and.(m.gt.1)) then
498  do k=1,nopes
499  if((nope.eq.20).or.(nope.eq.8)) then
500  ipointer=5*(ifaceq(k,ig)-1)
501  elseif((nope.eq.10).or.(nope.eq.4)) then
502  ipointer=5*(ifacet(k,ig)-1)
503  else
504  ipointer=5*(ifacew(k,ig)-1)
505  endif
506 !
507 ! F_A vector
508 !
509  ff(ipointer+1)=ff(ipointer+1)+
510  & shp2(4,k)*(h0l(2)*xsj2(3)-h0l(3)*xsj2(2))*weight
511  ff(ipointer+2)=ff(ipointer+2)+
512  & shp2(4,k)*(h0l(3)*xsj2(1)-h0l(1)*xsj2(3))*weight
513  ff(ipointer+3)=ff(ipointer+3)+
514  & shp2(4,k)*(h0l(1)*xsj2(2)-h0l(2)*xsj2(1))*weight
515 !
516  enddo
517  endif
518 !
519  do ii=1,nopes
520  if((nope.eq.20).or.(nope.eq.8)) then
521  ipointeri=5*(ifaceq(ii,ig)-1)
522  elseif((nope.eq.10).or.(nope.eq.4)) then
523  ipointeri=5*(ifacet(ii,ig)-1)
524  else
525  ipointeri=5*(ifacew(ii,ig)-1)
526  endif
527  do jj=1,nopes
528  if((nope.eq.20).or.(nope.eq.8)) then
529  ipointerj=5*(ifaceq(jj,ig)-1)
530  elseif((nope.eq.10).or.(nope.eq.4)) then
531  ipointerj=5*(ifacet(jj,ig)-1)
532  else
533  ipointerj=5*(ifacew(jj,ig)-1)
534  endif
535 !
536  if(m.gt.1) then
537 !
538 ! K_Aphi matrix
539 !
540  if((ipointeri+1).gt.ipointerj+5) cycle
541  s(ipointeri+1,ipointerj+5)=
542  & s(ipointeri+1,ipointerj+5)
543  & +shp2(4,jj)*(shp2(3,ii)*xsj2(2)
544  & -shp2(2,ii)*xsj2(3))
545  & *weight
546 !
547  if((ipointeri+2).gt.ipointerj+5) cycle
548  s(ipointeri+2,ipointerj+5)=
549  & s(ipointeri+2,ipointerj+5)
550  & +shp2(4,jj)*(shp2(1,ii)*xsj2(3)
551  & -shp2(3,ii)*xsj2(1))
552  & *weight
553 !
554  if((ipointeri+3).gt.ipointerj+5) cycle
555  s(ipointeri+3,ipointerj+5)=
556  & s(ipointeri+3,ipointerj+5)
557  & +shp2(4,jj)*(shp2(2,ii)*xsj2(1)
558  & -shp2(1,ii)*xsj2(2))
559  & *weight
560 !
561 ! K_phiA matrix
562 !
563  if(ipointeri+5.gt.(ipointerj+1)) cycle
564  s(ipointeri+5,ipointerj+1)=
565  & s(ipointeri+5,ipointerj+1)
566  & +shp2(4,ii)*(shp2(3,jj)*xsj2(2)
567  & -shp2(2,jj)*xsj2(3))
568  & *weight
569 !
570  if(ipointeri+5.gt.(ipointerj+2)) cycle
571  s(ipointeri+5,ipointerj+2)=
572  & s(ipointeri+5,ipointerj+2)
573  & +shp2(4,ii)*(shp2(1,jj)*xsj2(3)
574  & -shp2(3,jj)*xsj2(1))
575  & *weight
576 !
577  if(ipointeri+5.gt.(ipointerj+3)) cycle
578  s(ipointeri+5,ipointerj+3)=
579  & s(ipointeri+5,ipointerj+3)
580  & +shp2(4,ii)*(shp2(2,jj)*xsj2(1)
581  & -shp2(1,jj)*xsj2(2))
582  & *weight
583  endif
584  enddo
585  enddo
586 !
587  enddo
588 !
589  id=id-1
590  enddo
591 !
592  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
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
shape15w
subroutine shape15w(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape15w.f:20
idist
static ITG * idist
Definition: radflowload.c:39
linscal
subroutine linscal(scal, konl, nope, jj, scall, idim)
Definition: linscal.f:20
linvec
subroutine linvec(vec, konl, nope, jj, vecl, istart, iend)
Definition: linvec.f:20
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
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
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
materialdata_em
subroutine materialdata_em(elcon, nelcon, alcon, nalcon, imat, ntmat_, t1l, elconloc, ncmat_, alpha)
Definition: materialdata_em.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)