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

Go to the source code of this file.

Functions/Subroutines

subroutine resultsem (co, kon, ipkon, lakon, v, elcon, nelcon, ielmat, ntmat_, vini, dtime, matname, mi, ncmat_, nea, neb, sti, alcon, nalcon, h0, istartset, iendset, ialset, iactive, fn)
 

Function/Subroutine Documentation

◆ resultsem()

subroutine resultsem ( real*8, dimension(3,*)  co,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
character*8, dimension(*)  lakon,
real*8, dimension(0:mi(2),*)  v,
real*8, dimension(0:ncmat_,ntmat_,*)  elcon,
integer, dimension(2,*)  nelcon,
integer, dimension(mi(3),*)  ielmat,
integer  ntmat_,
real*8, dimension(0:mi(2),*)  vini,
real*8  dtime,
character*80, dimension(*)  matname,
integer, dimension(*)  mi,
integer  ncmat_,
integer  nea,
integer  neb,
real*8, dimension(6,mi(1),*)  sti,
real*8, dimension(0:6,ntmat_,*)  alcon,
integer, dimension(2,*)  nalcon,
real*8, dimension(3,*)  h0,
integer, dimension(*)  istartset,
integer, dimension(*)  iendset,
integer, dimension(*)  ialset,
integer, dimension(3)  iactive,
real*8, dimension(0:mi(2),*)  fn 
)
22 !
23 ! calculates the heat flux and the material tangent at the integration
24 ! points and the internal concentrated flux at the nodes
25 !
26  implicit none
27 !
28  character*8 lakon(*),lakonl
29  character*80 amat,matname(*)
30 !
31  integer kon(*),konl(26),mi(*),nelcon(2,*),ielmat(mi(3),*),
32  & ntmat_,ipkon(*),null,three,iflag,mt,i,j,k,m1,kk,i1,m3,indexe,
33  & nope,imat,mint3d,ncmat_,nea,neb,nalcon(2,*),mm,l,istart,iset,
34  & isurf,ilength,istartset(*),iendset(*),ialset(*),nopes,m,
35  & m2,ig,id,iactive(3),konl2(9),ifaceq(8,6),ifacet(6,4),iel,
36  & ifacew(8,5),mint2d,nfaces,one
37 !
38  real*8 co(3,*),v(0:mi(2),*),shp(4,26),xl(3,26),vl(0:mi(2),26),
39  & elcon(0:ncmat_,ntmat_,*),vkl(0:mi(2),3),vini(0:mi(2),*),c1,
40  & elconloc(21),xi,et,ze,xsj,t1l,dtime,weight,xsj2(3),shp2(7,9),
41  & vinikl(0:mi(2),3),alpha(6),vl2(0:mi(2),9),xl2(1:3,9),
42  & h0l(3),al(3),ainil(3),sti(6,mi(1),*),xs2(3,7),phi,
43  & um,alcon(0:6,ntmat_,*),h0(3,*),vinil(0:mi(2),26),
44  & fn(0:mi(2),*)
45 !
46  include "gauss.f"
47 !
48  data ifaceq /4,3,2,1,11,10,9,12,
49  & 5,6,7,8,13,14,15,16,
50  & 1,2,6,5,9,18,13,17,
51  & 2,3,7,6,10,19,14,18,
52  & 3,4,8,7,11,20,15,19,
53  & 4,1,5,8,12,17,16,20/
54  data ifacet /1,3,2,7,6,5,
55  & 1,2,4,5,9,8,
56  & 2,3,4,6,10,9,
57  & 1,4,3,8,10,7/
58  data ifacew /1,3,2,9,8,7,0,0,
59  & 4,5,6,10,11,12,0,0,
60  & 1,2,5,4,7,14,10,13,
61  & 2,3,6,5,8,15,11,14,
62  & 4,6,3,1,12,15,9,13/
63 !
64  iflag=3
65  null=0
66  one=1
67  three=3
68 !
69  mt=mi(2)+1
70 !
71 ! calculation of temperatures and thermal flux
72 !
73  do i=nea,neb
74 !
75  lakonl(1:8)=lakon(i)(1:8)
76 !
77  if(ipkon(i).lt.0) cycle
78 !
79  imat=ielmat(1,i)
80  amat=matname(imat)
81 !
82  indexe=ipkon(i)
83  if(lakonl(4:5).eq.'20') then
84  nope=20
85  nopes=8
86  elseif(lakonl(4:4).eq.'8') then
87  nope=8
88  nopes=4
89  elseif(lakonl(4:5).eq.'10') then
90  nope=10
91  nopes=6
92  elseif(lakonl(4:4).eq.'4') then
93  nope=4
94  nopes=3
95  elseif(lakonl(4:5).eq.'15') then
96  nope=15
97  elseif(lakonl(4:4).eq.'6') then
98  nope=6
99  else
100  cycle
101  endif
102 !
103  if(lakonl(4:5).eq.'8R') then
104  mint3d=1
105  mint2d=1
106  elseif((lakonl(4:4).eq.'8').or.
107  & (lakonl(4:6).eq.'20R')) then
108  mint3d=8
109  mint2d=4
110  elseif(lakonl(4:4).eq.'2') then
111  mint3d=27
112  mint2d=9
113  elseif(lakonl(4:5).eq.'10') then
114  mint3d=4
115  mint2d=3
116  elseif(lakonl(4:4).eq.'4') then
117  mint3d=1
118  mint2d=1
119  elseif(lakonl(4:5).eq.'15') then
120  mint3d=9
121  elseif(lakonl(4:4).eq.'6') then
122  mint3d=2
123  endif
124 !
125  do j=1,nope
126  konl(j)=kon(indexe+j)
127  do k=1,3
128  xl(k,j)=co(k,konl(j))
129  enddo
130  do k=1,5
131  vl(k,j)=v(k,konl(j))
132  enddo
133  vinil(4,j)=vini(4,konl(j))
134  enddo
135 !
136  do kk=1,mint3d
137  if(lakonl(4:5).eq.'8R') then
138  xi=gauss3d1(1,kk)
139  et=gauss3d1(2,kk)
140  ze=gauss3d1(3,kk)
141  weight=weight3d1(kk)
142  elseif((lakonl(4:4).eq.'8').or.
143  & (lakonl(4:6).eq.'20R'))
144  & then
145  xi=gauss3d2(1,kk)
146  et=gauss3d2(2,kk)
147  ze=gauss3d2(3,kk)
148  weight=weight3d2(kk)
149  elseif(lakonl(4:4).eq.'2') then
150  xi=gauss3d3(1,kk)
151  et=gauss3d3(2,kk)
152  ze=gauss3d3(3,kk)
153  weight=weight3d3(kk)
154  elseif(lakonl(4:5).eq.'10') then
155  xi=gauss3d5(1,kk)
156  et=gauss3d5(2,kk)
157  ze=gauss3d5(3,kk)
158  weight=weight3d5(kk)
159  elseif(lakonl(4:4).eq.'4') then
160  xi=gauss3d4(1,kk)
161  et=gauss3d4(2,kk)
162  ze=gauss3d4(3,kk)
163  weight=weight3d4(kk)
164  elseif(lakonl(4:5).eq.'15') then
165  xi=gauss3d8(1,kk)
166  et=gauss3d8(2,kk)
167  ze=gauss3d8(3,kk)
168  weight=weight3d8(kk)
169  elseif(lakonl(4:4).eq.'6') then
170  xi=gauss3d7(1,kk)
171  et=gauss3d7(2,kk)
172  ze=gauss3d7(3,kk)
173  weight=weight3d7(kk)
174  endif
175 !
176  if(nope.eq.20) then
177  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
178  elseif(nope.eq.8) then
179  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
180  elseif(nope.eq.10) then
181  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
182  elseif(nope.eq.4) then
183  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
184  elseif(nope.eq.15) then
185  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
186  else
187  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
188  endif
189 !
190  c1=xsj*weight
191 !
192 ! vkl(m2,m3) contains the derivative of the m2-
193 ! component of the displacement with respect to
194 ! direction m3
195 !
196  do k=1,5
197  do m3=1,3
198  vkl(k,m3)=0.d0
199  enddo
200 !
201  do m1=1,nope
202  do m3=1,3
203  vkl(k,m3)=vkl(k,m3)+shp(m3,m1)*vl(k,m1)
204  enddo
205  enddo
206  enddo
207 !
208  do m3=1,3
209  vinikl(4,m3)=0.d0
210  enddo
211  do m1=1,nope
212  do m3=1,3
213  vinikl(4,m3)=vinikl(4,m3)+shp(m3,m1)*vinil(4,m1)
214  enddo
215  enddo
216 !
217 ! calculating the temperature difference in
218 ! the integration point
219 !
220  t1l=0.d0
221  do j=1,3
222  h0l(j)=0.d0
223  al(j)=0.d0
224  ainil(j)=0.d0
225  enddo
226  if(lakonl(4:5).eq.'8 ') then
227  do i1=1,8
228  do j=1,3
229  h0l(j)=h0l(j)+h0(j,konl(i1))/8.d0
230  al(j)=al(j)+v(j,konl(i1))/8.d0
231  ainil(j)=ainil(j)+vini(j,konl(i1))/8.d0
232  enddo
233  t1l=t1l+v(0,konl(i1))/8.d0
234  enddo
235  elseif(lakonl(4:6).eq.'20 ') then
236  call linscal(v,konl,nope,kk,t1l,mi(2))
237  call linvec(h0,konl,nope,kk,h0l,one,three)
238  call linvec(v,konl,nope,kk,al,null,mi(2))
239  call linvec(vini,konl,nope,kk,ainil,null,mi(2))
240  elseif(lakonl(4:6).eq.'10T') then
241  call linscal10(v,konl,t1l,mi(2),shp)
242  call linvec10(h0,konl,h0l,one,three,shp)
243  call linvec10(v,konl,al,null,mi(2),shp)
244  call linvec10(vini,konl,ainil,null,mi(2),shp)
245  else
246  do i1=1,nope
247  t1l=t1l+shp(4,i1)*v(0,konl(i1))
248  do j=1,3
249  h0l(j)=h0l(j)+shp(4,i1)*h0(j,konl(i1))
250  al(j)=al(j)+shp(4,i1)*v(j,konl(i1))
251  ainil(j)=ainil(j)+shp(4,i1)*vini(j,konl(i1))
252  enddo
253  enddo
254  endif
255 !
256 ! material data (permeability)
257 !
258  call materialdata_em(elcon,nelcon,alcon,nalcon,
259  & imat,ntmat_,t1l,elconloc,ncmat_,alpha)
260 !
261  um=elconloc(1)
262 !
263  if(int(elconloc(2)).eq.1) then
264 !
265 ! magnetic field B in phi-domain
266 !
267  do k=1,3
268  sti(k+3,kk,i)=um*(h0l(k)-vkl(5,k))
269  enddo
270 !
271 ! calculating the electromagnetic force K_phiphi
272 !
273  do m1=1,nope
274  do m3=1,3
275  fn(5,konl(m1))=fn(5,konl(m1))-c1*
276  & um*shp(m3,m1)*vkl(5,m3)
277  enddo
278  enddo
279  else
280 !
281 ! magnetic field B in A and A-V domain
282 !
283  sti(4,kk,i)=vkl(3,2)-vkl(2,3)
284  sti(5,kk,i)=vkl(1,3)-vkl(3,1)
285  sti(6,kk,i)=vkl(2,1)-vkl(1,2)
286 !
287 ! electric field E in A-V domain
288 !
289  if(int(elconloc(2)).eq.2) then
290  do k=1,3
291  sti(k,kk,i)=(ainil(k)-al(k)+
292  & vinikl(4,k)-vkl(4,k))/dtime
293  enddo
294  endif
295 !
296 ! calculating the electromagnetic force K_AA
297 !
298  do m1=1,nope
299  do m2=1,3
300  do m3=1,3
301  fn(m2,konl(m1))=fn(m2,konl(m1))+c1*
302  & (shp(m3,m1)*(vkl(m2,m3)-vkl(m3,m2))+
303  & shp(m2,m1)*vkl(m3,m3))/um
304  enddo
305  enddo
306  enddo
307 !
308  endif
309 !
310  enddo
311 !
312 ! surface integrals
313 !
314 ! determining the number of faces per element
315 !
316  if((lakonl(4:4).eq.'8').or.(lakonl(4:4).eq.'2')) then
317  nfaces=6
318  elseif((lakonl(4:4).eq.'6').or.(lakonl(4:5).eq.'15')) then
319  nfaces=5
320  elseif((lakonl(4:4).eq.'4').or.(lakonl(4:5).eq.'10')) then
321  nfaces=4
322  endif
323 !
324  m=int(elconloc(2))
325  if(iactive(m).eq.0) cycle
326  iset=iactive(m)
327  istart=istartset(iset)
328  ilength=iendset(iset)-istart+1
329 !
330  isurf=10*i+nfaces
331  call nident(ialset(istart),isurf,ilength,id)
332 !
333  do
334  if(id.eq.0) exit
335  isurf=ialset(istart+id-1)
336  iel=int(isurf/10.d0)
337  if(iel.ne.i) exit
338  ig=isurf-10*iel
339 !
340 ! treatment of wedge faces
341 !
342  if(lakonl(4:4).eq.'6') then
343  mint2d=1
344  if(ig.le.2) then
345  nopes=3
346  else
347  nopes=4
348  endif
349  endif
350  if(lakonl(4:5).eq.'15') then
351  if(ig.le.2) then
352  mint2d=3
353  nopes=6
354  else
355  mint2d=4
356  nopes=8
357  endif
358  endif
359 !
360 ! face connectivity is stored in konl2
361 !
362  if((nope.eq.20).or.(nope.eq.8)) then
363  do j=1,nopes
364  konl2(j)=konl(ifaceq(j,ig))
365  enddo
366  elseif((nope.eq.10).or.(nope.eq.4)) then
367  do j=1,nopes
368  konl2(j)=konl(ifacet(j,ig))
369  enddo
370  else
371  do j=1,nopes
372  konl2(j)=konl(ifacew(j,ig))
373  enddo
374  endif
375 !
376 ! face coordinates and solution fields
377 !
378  do j=1,nopes
379  do k=1,3
380  xl2(k,j)=co(k,konl2(j))
381  vl2(k,j)=v(k,konl2(j))
382  enddo
383  vl2(5,j)=v(5,konl2(j))
384  enddo
385 !
386  do kk=1,mint2d
387 !
388  if((lakonl(4:5).eq.'8R').or.
389  & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
390  xi=gauss2d1(1,kk)
391  et=gauss2d1(2,kk)
392  weight=weight2d1(kk)
393  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')
394  & .or.((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
395  xi=gauss2d2(1,kk)
396  et=gauss2d2(2,kk)
397  weight=weight2d2(kk)
398  elseif(lakonl(4:4).eq.'2') then
399  xi=gauss2d3(1,kk)
400  et=gauss2d3(2,kk)
401  weight=weight2d3(kk)
402  elseif((lakonl(4:5).eq.'10').or.
403  & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
404  xi=gauss2d5(1,kk)
405  et=gauss2d5(2,kk)
406  weight=weight2d5(kk)
407  elseif((lakonl(4:4).eq.'4').or.
408  & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
409  xi=gauss2d4(1,kk)
410  et=gauss2d4(2,kk)
411  weight=weight2d4(kk)
412  endif
413 !
414  if(nopes.eq.8) then
415  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
416  elseif(nopes.eq.4) then
417  call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
418  elseif(nopes.eq.6) then
419  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
420  else
421  call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
422  endif
423 !
424 ! determining the electromagnetic force
425 !
426  if(m.gt.1) then
427 !
428 ! determining the values of phi
429 !
430  phi=0.d0
431  do m1=1,nopes
432  phi=phi+shp2(4,m1)*vl2(5,m1)
433  enddo
434 !
435  do m1=1,nopes
436  do k=1,3
437  l=k+1
438  if(l.gt.3) l=1
439  mm=l+1
440  if(mm.gt.3) mm=1
441  fn(k,konl2(m1))=fn(k,konl2(m1))-
442  & (shp2(l,m1)*xsj2(mm)-shp2(mm,m1)*xsj2(l))
443  & *phi*weight
444  enddo
445  enddo
446  elseif(m.eq.1) then
447 !
448 ! determining the derivative of A w.r.t. x, y and z
449 !
450  do k=1,3
451  do m3=1,3
452  vkl(k,m3)=0.d0
453  enddo
454  do m1=1,nopes
455  do m3=1,3
456  vkl(k,m3)=vkl(k,m3)+shp2(m3,m1)*vl2(k,m1)
457  enddo
458  enddo
459  enddo
460 !
461  do m1=1,nopes
462  do k=1,3
463  l=k+1
464  if(l.gt.3) l=1
465  mm=l+1
466  if(mm.gt.3) mm=1
467  fn(5,konl2(m1))=fn(5,konl2(m1))+
468  & shp2(4,m1)*(vkl(mm,l)-vkl(l,mm))*xsj2(k)
469  & *weight
470  enddo
471  enddo
472  endif
473 !
474  enddo
475  id=id-1
476  enddo
477  enddo
478 !
479  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
c1
static double * c1
Definition: mafillvcompmain.c:30
linscal
subroutine linscal(scal, konl, nope, jj, scall, idim)
Definition: linscal.f:20
linvec10
subroutine linvec10(vec, konl, vecl, istart, iend, shp)
Definition: linvec10.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)