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

Go to the source code of this file.

Functions/Subroutines

subroutine printoutface (co, rhcon, nrhcon, ntmat_, vold, shcon, nshcon, cocon, ncocon, icompressible, istartset, iendset, ipkonf, lakonf, konf, ialset, prset, ttime, nset, set, nprint, prlab, ielmat, mi, ithermal, nactdoh, icfd, time, stn)
 

Function/Subroutine Documentation

◆ printoutface()

subroutine printoutface ( real*8, dimension(3,*)  co,
real*8, dimension(0:1,ntmat_,*)  rhcon,
integer, dimension(*)  nrhcon,
integer  ntmat_,
real*8, dimension(0:mi(2),*)  vold,
real*8, dimension(0:3,ntmat_,*)  shcon,
integer, dimension(*)  nshcon,
real*8, dimension(0:6,ntmat_,*)  cocon,
integer, dimension(2,*)  ncocon,
integer  icompressible,
integer, dimension(*)  istartset,
integer, dimension(*)  iendset,
integer, dimension(*)  ipkonf,
character*8, dimension(*)  lakonf,
integer, dimension(*)  konf,
integer, dimension(*)  ialset,
character*81, dimension(*)  prset,
real*8  ttime,
integer  nset,
character*81, dimension(*)  set,
integer  nprint,
character*6, dimension(*)  prlab,
integer, dimension(mi(3),*)  ielmat,
integer, dimension(*)  mi,
integer  ithermal,
integer, dimension(*)  nactdoh,
integer  icfd,
real*8  time,
real*8, dimension(6,*)  stn 
)
23 !
24 ! calculation and printout of the lift and drag forces
25 !
26  implicit none
27 !
28  integer icompressible
29 !
30  character*8 lakonl,lakonf(*)
31  character*6 prlab(*)
32  character*80 faset
33  character*81 set(*),prset(*)
34 !
35  integer konl(20),ifaceq(8,6),nelem,ii,nprint,i,j,i1,i2,j1,
36  & ncocon(2,*),k1,jj,ig,nrhcon(*),nshcon(*),ntmat_,nope,nopes,imat,
37  & mint2d,ifacet(6,4),ifacew(8,5),iflag,indexe,jface,istartset(*),
38  & iendset(*),ipkonf(*),konf(*),iset,ialset(*),nset,ipos,ithermal,
39  & mi(*),ielmat(mi(3),*),nactdoh(*),icfd,nelemcfd
40 !
41  real*8 co(3,*),xl(3,20),shp(4,20),xs2(3,7),dvi,f(0:3),time,
42  & vkl(0:3,3),rhcon(0:1,ntmat_,*),t(3,3),div,shcon(0:3,ntmat_,*),
43  & voldl(0:mi(2),20),cocon(0:6,ntmat_,*),xl2(3,8),xsj2(3),
44  & shp2(7,8),vold(0:mi(2),*),xi,et,xsj,temp,xi3d,et3d,ze3d,weight,
45  & xlocal20(3,9,6),xlocal4(3,1,4),xlocal10(3,3,4),xlocal6(3,1,5),
46  & xlocal15(3,4,5),xlocal8(3,4,6),xlocal8r(3,1,6),ttime,pres,
47  & tf(0:3),tn,tt,dd,coords(3),cond,stn(6,*),xm(3),df(3),cg(3),
48  & area,xn(3),xnormforc,shearforc
49 !
50  include "gauss.f"
51  include "xlocal.f"
52 !
53  data ifaceq /4,3,2,1,11,10,9,12,
54  & 5,6,7,8,13,14,15,16,
55  & 1,2,6,5,9,18,13,17,
56  & 2,3,7,6,10,19,14,18,
57  & 3,4,8,7,11,20,15,19,
58  & 4,1,5,8,12,17,16,20/
59  data ifacet /1,3,2,7,6,5,
60  & 1,2,4,5,9,8,
61  & 2,3,4,6,10,9,
62  & 1,4,3,8,10,7/
63  data ifacew /1,3,2,9,8,7,0,0,
64  & 4,5,6,10,11,12,0,0,
65  & 1,2,5,4,7,14,10,13,
66  & 2,3,6,5,8,15,11,14,
67  & 4,6,3,1,12,15,9,13/
68  data iflag /3/
69 !
70 ! flag to check whether sof and/or som output was already
71 ! printed
72 !
73  do ii=1,nprint
74 !
75 ! check whether there are facial print requests
76 !
77 ! DRAG: drag forces in cfd-calculations
78 ! FLUX: heat flux
79 ! SOF: forces and moments on a section
80 ! (= an internal or external surface)
81 !
82  if((prlab(ii)(1:4).eq.'DRAG').or.(prlab(ii)(1:4).eq.'FLUX').or.
83  & (prlab(ii)(1:3).eq.'SOF'))
84  & then
85 !
86  ipos=index(prset(ii),' ')
87  faset=' '
88  faset(1:ipos-1)=prset(ii)(1:ipos-1)
89 !
90 ! printing the header
91 !
92  write(5,*)
93  if(prlab(ii)(1:4).eq.'DRAG') then
94 !
95 ! initialisierung forces
96 !
97  do i=1,3
98  f(i)=0.d0
99  enddo
100 !
101  write(5,120) faset(1:ipos-2),ttime+time
102  120 format(
103  & ' surface stress at the integration points for set ',
104  & a,' and time ',e14.7)
105  write(5,*)
106  write(5,124)
107  124 format(' el fa int tx ty
108  & tz normal stress shear stress and coordinates
109  &')
110  elseif(prlab(ii)(1:4).eq.'FLUX') then
111 !
112 ! initialisierung of the flux
113 !
114  f(0)=0.d0
115  write(5,121) faset(1:ipos-2),ttime+time
116  121 format(' heat flux at the integration points for set ',
117  & a,' and time ',e14.7)
118  write(5,*)
119  write(5,125)
120  125 format(' el fa int heat flux q and c
121  &oordinates')
122  else
123 !
124 ! initialize total force and total moment
125 !
126  do i=1,3
127  f(i)=0.d0
128  xm(i)=0.d0
129  cg(i)=0.d0
130  xn(i)=0.d0
131  enddo
132  area=0.d0
133  endif
134  write(5,*)
135 !
136 ! printing the data
137 !
138  do iset=1,nset
139  if(set(iset).eq.prset(ii)) exit
140  enddo
141 !
142  do jj=istartset(iset),iendset(iset)
143 !
144  jface=ialset(jj)
145 !
146  nelem=int(jface/10.d0)
147  ig=jface-10*nelem
148 c write(*,*) 'printoutface elem ig ',nelem,ig
149 !
150 ! for CFD calculations the elements were renumbered
151 !
152  if(icfd.eq.1) then
153  nelemcfd=nactdoh(nelem)
154  indexe=ipkonf(nelemcfd)
155  lakonl=lakonf(nelemcfd)
156  imat=ielmat(1,nelemcfd)
157  else
158  indexe=ipkonf(nelem)
159  lakonl=lakonf(nelem)
160  imat=ielmat(1,nelem)
161  endif
162 !
163  if(lakonl(4:4).eq.'2') then
164  nope=20
165  nopes=8
166  elseif(lakonl(4:4).eq.'8') then
167  nope=8
168  nopes=4
169  elseif(lakonl(4:5).eq.'10') then
170  nope=10
171  nopes=6
172  elseif(lakonl(4:4).eq.'4') then
173  nope=4
174  nopes=3
175  elseif(lakonl(4:5).eq.'15') then
176  nope=15
177  elseif(lakonl(4:4).eq.'6') then
178  nope=6
179  endif
180 !
181  if(lakonl(4:5).eq.'8R') then
182  mint2d=1
183  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R'))
184  & then
185  if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'S').or.
186  & (lakonl(7:7).eq.'E')) then
187  mint2d=2
188  else
189  mint2d=4
190  endif
191  elseif(lakonl(4:4).eq.'2') then
192  mint2d=9
193  elseif(lakonl(4:5).eq.'10') then
194  mint2d=3
195  elseif(lakonl(4:4).eq.'4') then
196  mint2d=1
197  endif
198 !
199 ! local topology
200 !
201  do i=1,nope
202  konl(i)=konf(indexe+i)
203  enddo
204 !
205 ! computation of the coordinates of the local nodes
206 !
207  do i=1,nope
208  do j=1,3
209  xl(j,i)=co(j,konl(i))
210  enddo
211  enddo
212 !
213 ! temperature, displacement (structures) or
214 ! temperature, velocities and pressure (fluids)
215 !
216  do i1=1,nope
217  do i2=0,mi(2)
218  voldl(i2,i1)=vold(i2,konl(i1))
219  enddo
220  enddo
221 !
222 ! for structural calculations: adding the displacements to the
223 ! coordinates
224 !
225  if(icfd.eq.0) then
226  do i=1,nope
227  do j=1,3
228  xl(j,i)=xl(j,i)+voldl(j,i)
229  enddo
230  enddo
231  endif
232 !
233 ! treatment of wedge faces
234 !
235  if(lakonl(4:4).eq.'6') then
236  mint2d=1
237  if(ig.le.2) then
238  nopes=3
239  else
240  nopes=4
241  endif
242  endif
243  if(lakonl(4:5).eq.'15') then
244  if(ig.le.2) then
245  mint2d=3
246  nopes=6
247  else
248  mint2d=4
249  nopes=8
250  endif
251  endif
252 !
253  if(icfd.ne.0) then
254 !
255 ! CFD: undeformed structure
256 !
257  if((nope.eq.20).or.(nope.eq.8)) then
258  do i=1,nopes
259  do j=1,3
260  xl2(j,i)=co(j,konl(ifaceq(i,ig)))
261  enddo
262  enddo
263  elseif((nope.eq.10).or.(nope.eq.4)) then
264  do i=1,nopes
265  do j=1,3
266  xl2(j,i)=co(j,konl(ifacet(i,ig)))
267  enddo
268  enddo
269  else
270  do i=1,nopes
271  do j=1,3
272  xl2(j,i)=co(j,konl(ifacew(i,ig)))
273  enddo
274  enddo
275  endif
276  else
277 !
278 ! no CFD: deformed structure
279 !
280  if((nope.eq.20).or.(nope.eq.8)) then
281  do i=1,nopes
282  do j=1,3
283  xl2(j,i)=co(j,konl(ifaceq(i,ig)))
284  & +vold(j,konl(ifaceq(i,ig)))
285  enddo
286  enddo
287  elseif((nope.eq.10).or.(nope.eq.4)) then
288  do i=1,nopes
289  do j=1,3
290  xl2(j,i)=co(j,konl(ifacet(i,ig)))
291  & +vold(j,konl(ifacet(i,ig)))
292  enddo
293  enddo
294  else
295  do i=1,nopes
296  do j=1,3
297  xl2(j,i)=co(j,konl(ifacew(i,ig)))+
298  & vold(j,konl(ifacew(i,ig)))
299  enddo
300  enddo
301  endif
302  endif
303 !
304  do i=1,mint2d
305 !
306 ! local coordinates of the surface integration
307 ! point within the surface local coordinate system
308 !
309  if((lakonl(4:5).eq.'8R').or.
310  & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
311  xi=gauss2d1(1,i)
312  et=gauss2d1(2,i)
313  weight=weight2d1(i)
314  elseif((lakonl(4:4).eq.'8').or.
315  & (lakonl(4:6).eq.'20R').or.
316  & ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
317  xi=gauss2d2(1,i)
318  et=gauss2d2(2,i)
319  weight=weight2d2(i)
320  elseif(lakonl(4:4).eq.'2') then
321  xi=gauss2d3(1,i)
322  et=gauss2d3(2,i)
323  weight=weight2d3(i)
324  elseif((lakonl(4:5).eq.'10').or.
325  & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
326  xi=gauss2d5(1,i)
327  et=gauss2d5(2,i)
328  weight=weight2d5(i)
329  elseif((lakonl(4:4).eq.'4').or.
330  & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
331  xi=gauss2d4(1,i)
332  et=gauss2d4(2,i)
333  weight=weight2d4(i)
334  endif
335 !
336 ! local surface normal
337 !
338  if(nopes.eq.8) then
339  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
340  elseif(nopes.eq.4) then
341  call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
342  elseif(nopes.eq.6) then
343  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
344  else
345  call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
346  endif
347 !
348 ! global coordinates of the integration point
349 !
350  do j1=1,3
351  coords(j1)=0.d0
352  do i1=1,nopes
353  coords(j1)=coords(j1)+shp2(4,i1)*xl2(j1,i1)
354  enddo
355  enddo
356 !
357 ! local coordinates of the surface integration
358 ! point within the element local coordinate system
359 !
360  if((prlab(ii)(1:4).eq.'DRAG').or.
361  & (prlab(ii)(1:4).eq.'FLUX')) then
362 !
363 ! deformation gradient is only needed for
364 ! DRAG and FLUX applications
365 !
366  if(lakonl(4:5).eq.'8R') then
367  xi3d=xlocal8r(1,i,ig)
368  et3d=xlocal8r(2,i,ig)
369  ze3d=xlocal8r(3,i,ig)
370  call shape8h(xi3d,et3d,ze3d,xl,xsj,shp,iflag)
371  elseif(lakonl(4:4).eq.'8') then
372  xi3d=xlocal8(1,i,ig)
373  et3d=xlocal8(2,i,ig)
374  ze3d=xlocal8(3,i,ig)
375  call shape8h(xi3d,et3d,ze3d,xl,xsj,shp,iflag)
376  elseif(lakonl(4:6).eq.'20R') then
377  xi3d=xlocal8(1,i,ig)
378  et3d=xlocal8(2,i,ig)
379  ze3d=xlocal8(3,i,ig)
380  call shape20h(xi3d,et3d,ze3d,xl,xsj,shp,iflag)
381  elseif(lakonl(4:4).eq.'2') then
382  xi3d=xlocal20(1,i,ig)
383  et3d=xlocal20(2,i,ig)
384  ze3d=xlocal20(3,i,ig)
385  call shape20h(xi3d,et3d,ze3d,xl,xsj,shp,iflag)
386  elseif(lakonl(4:5).eq.'10') then
387  xi3d=xlocal10(1,i,ig)
388  et3d=xlocal10(2,i,ig)
389  ze3d=xlocal10(3,i,ig)
390  call shape10tet(xi3d,et3d,ze3d,xl,xsj,shp,iflag)
391  elseif(lakonl(4:4).eq.'4') then
392  xi3d=xlocal4(1,i,ig)
393  et3d=xlocal4(2,i,ig)
394  ze3d=xlocal4(3,i,ig)
395  call shape4tet(xi3d,et3d,ze3d,xl,xsj,shp,iflag)
396  elseif(lakonl(4:5).eq.'15') then
397  xi3d=xlocal15(1,i,ig)
398  et3d=xlocal15(2,i,ig)
399  ze3d=xlocal15(3,i,ig)
400  call shape15w(xi3d,et3d,ze3d,xl,xsj,shp,iflag)
401  elseif(lakonl(4:4).eq.'6') then
402  xi3d=xlocal6(1,i,ig)
403  et3d=xlocal6(2,i,ig)
404  ze3d=xlocal6(3,i,ig)
405  call shape6w(xi3d,et3d,ze3d,xl,xsj,shp,iflag)
406  endif
407  endif
408 !
409 ! calculating of
410 ! the temperature temp
411 ! the static pressure pres
412 ! the velocity gradient vkl
413 ! in the integration point
414 !
415  if(prlab(ii)(1:4).eq.'DRAG') then
416  temp=0.d0
417  pres=0.d0
418  do i1=1,3
419  do j1=1,3
420  vkl(i1,j1)=0.d0
421  enddo
422  enddo
423  do i1=1,nope
424  temp=temp+shp(4,i1)*voldl(0,i1)
425  pres=pres+shp(4,i1)*voldl(4,i1)
426  do j1=1,3
427  do k1=1,3
428  vkl(j1,k1)=vkl(j1,k1)+
429  & shp(k1,i1)*voldl(j1,i1)
430  enddo
431  enddo
432  enddo
433  if(icompressible.eq.1)
434  & div=vkl(1,1)+vkl(2,2)+vkl(3,3)
435 !
436 ! material data (density, dynamic viscosity, heat capacity and
437 ! conductivity)
438 !
439  call materialdata_dvi(shcon,nshcon,imat,dvi,temp,
440  & ntmat_,ithermal)
441 !
442 ! determining the stress
443 !
444  do i1=1,3
445  do j1=1,3
446  t(i1,j1)=vkl(i1,j1)+vkl(j1,i1)
447  enddo
448  if(icompressible.eq.1)
449  & t(i1,i1)=t(i1,i1)-2.d0*div/3.d0
450  enddo
451 !
452  dd=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
453  & xsj2(3)*xsj2(3))
454  do i1=1,3
455  tf(i1)=dvi*(t(i1,1)*xsj2(1)+t(i1,2)*xsj2(2)+
456  & t(i1,3)*xsj2(3))-pres*xsj2(i1)
457  f(i1)=f(i1)+tf(i1)*weight
458  tf(i1)=tf(i1)/dd
459  enddo
460  tn=(tf(1)*xsj2(1)+tf(2)*xsj2(2)+tf(3)*xsj2(3))/dd
461  tt=dsqrt((tf(1)-tn*xsj2(1)/dd)**2+
462  & (tf(2)-tn*xsj2(2)/dd)**2+
463  & (tf(3)-tn*xsj2(3)/dd)**2)
464  write(5,'(i10,1x,i3,1x,i3,1p,8(1x,e13.6))')nelem,
465  & ig,i,(tf(i1),i1=1,3),tn,tt,(coords(i1),i1=1,3)
466 !
467  elseif(prlab(ii)(1:4).eq.'FLUX') then
468  temp=0.d0
469  do j1=1,3
470  vkl(0,j1)=0.d0
471  enddo
472  do i1=1,nope
473  temp=temp+shp(4,i1)*voldl(0,i1)
474  do k1=1,3
475  vkl(0,k1)=vkl(0,k1)+shp(k1,i1)*voldl(0,i1)
476  enddo
477  enddo
478 !
479 ! material data (conductivity)
480 !
481  call materialdata_cond(imat,ntmat_,temp,cocon,
482  & ncocon,cond)
483 !
484 ! determining the stress
485 !
486  dd=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
487  & xsj2(3)*xsj2(3))
488 !
489  tf(0)=-cond*(vkl(0,1)*xsj2(1)+
490  & vkl(0,2)*xsj2(2)+
491  & vkl(0,3)*xsj2(3))
492  f(0)=f(0)+tf(0)*weight
493  tf(0)=tf(0)/dd
494  write(5,'(i10,1x,i3,1x,i3,1p,4(1x,e13.6))')nelem,
495  & ig,i,tf(0),(coords(i1),i1=1,3)
496 !
497  else
498 !
499 ! forces and moments in sections
500 ! These are obtained by integrating the stress tensor;
501 ! the stresses at the integration points are interpolated
502 ! from the values at the nodes of the face; these values
503 ! were extrapolated from the integration point values within
504 ! the element. This approach is different from the one
505 ! under "DRAG" because here the stress-strain relationship
506 ! can be highly nonlinear (plasticity..); this is not the
507 ! case in cfd.
508 !
509 ! interpolation of the stress tensor at the integration
510 ! point
511 !
512  do i1=1,3
513  do j1=i1,3
514  t(i1,j1)=0.d0
515  enddo
516  enddo
517 !
518  if((nope.eq.20).or.(nope.eq.8)) then
519  do i1=1,nopes
520  t(1,1)=t(1,1)+
521  & shp2(4,i1)*stn(1,konl(ifaceq(i1,ig)))
522  t(2,2)=t(2,2)+
523  & shp2(4,i1)*stn(2,konl(ifaceq(i1,ig)))
524  t(3,3)=t(3,3)+
525  & shp2(4,i1)*stn(3,konl(ifaceq(i1,ig)))
526  t(1,2)=t(1,2)+
527  & shp2(4,i1)*stn(4,konl(ifaceq(i1,ig)))
528  t(1,3)=t(1,3)+
529  & shp2(4,i1)*stn(5,konl(ifaceq(i1,ig)))
530  t(2,3)=t(2,3)+
531  & shp2(4,i1)*stn(6,konl(ifaceq(i1,ig)))
532  enddo
533  elseif((nope.eq.10).or.(nope.eq.4)) then
534  do i1=1,nopes
535  t(1,1)=t(1,1)+
536  & shp2(4,i1)*stn(1,konl(ifacet(i1,ig)))
537  t(2,2)=t(2,2)+
538  & shp2(4,i1)*stn(2,konl(ifacet(i1,ig)))
539  t(3,3)=t(3,3)+
540  & shp2(4,i1)*stn(3,konl(ifacet(i1,ig)))
541  t(1,2)=t(1,2)+
542  & shp2(4,i1)*stn(4,konl(ifacet(i1,ig)))
543  t(1,3)=t(1,3)+
544  & shp2(4,i1)*stn(5,konl(ifacet(i1,ig)))
545  t(2,3)=t(2,3)+
546  & shp2(4,i1)*stn(6,konl(ifacet(i1,ig)))
547  enddo
548  else
549  do i1=1,nopes
550  t(1,1)=t(1,1)+
551  & shp2(4,i1)*stn(1,konl(ifacew(i1,ig)))
552  t(2,2)=t(2,2)+
553  & shp2(4,i1)*stn(2,konl(ifacew(i1,ig)))
554  t(3,3)=t(3,3)+
555  & shp2(4,i1)*stn(3,konl(ifacew(i1,ig)))
556  t(1,2)=t(1,2)+
557  & shp2(4,i1)*stn(4,konl(ifacew(i1,ig)))
558  t(1,3)=t(1,3)+
559  & shp2(4,i1)*stn(5,konl(ifacew(i1,ig)))
560  t(2,3)=t(2,3)+
561  & shp2(4,i1)*stn(6,konl(ifacew(i1,ig)))
562  enddo
563  endif
564  t(2,1)=t(1,2)
565  t(3,1)=t(1,3)
566  t(3,2)=t(2,3)
567 !
568 ! calculating the force contribution
569 !
570  do i1=1,3
571  df(i1)=(t(i1,1)*xsj2(1)+
572  & t(i1,2)*xsj2(2)+
573  & t(i1,3)*xsj2(3))*weight
574  enddo
575 !
576 ! update total force and total moment
577 !
578  do i1=1,3
579  f(i1)=f(i1)+df(i1)
580  enddo
581  xm(1)=xm(1)+coords(2)*df(3)-coords(3)*df(2)
582  xm(2)=xm(2)+coords(3)*df(1)-coords(1)*df(3)
583  xm(3)=xm(3)+coords(1)*df(2)-coords(2)*df(1)
584 !
585 ! update area, center of gravity and mean normal
586 !
587  dd=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
588  & xsj2(3)*xsj2(3))
589  area=area+dd
590  do i1=1,3
591  cg(i1)=cg(i1)+coords(i1)*dd
592  xn(i1)=xn(i1)+xsj2(i1)
593  enddo
594  endif
595  enddo
596  enddo
597 !
598  if(prlab(ii)(1:4).eq.'DRAG') then
599  write(5,*)
600  write(5,122) faset(1:ipos-2),ttime+time
601  122 format(' total surface force (fx,fy,fz) for set ',a,
602  & ' and time ',e14.7)
603  write(5,*)
604  write(5,'(1p,3(1x,e13.6))') (f(j),j=1,3)
605  elseif(prlab(ii)(1:4).eq.'FLUX') then
606  write(5,*)
607  write(5,123) faset(1:ipos-2),ttime+time
608  123 format(' total surface flux (q) for set ',a,
609  & ' and time ',e14.7)
610  write(5,*)
611  write(5,'(1p,1x,e13.6)') f(0)
612  else
613 !
614 ! surface set statistics
615 !
616  write(5,*)
617  write(5,130) faset(1:ipos-2),ttime+time
618  130 format(' statistics for surface set ',a,
619  & ' and time ',e14.7)
620 !
621 ! total force and moment about the origin
622 !
623  write(5,*)
624  write(5,126)
625 !
626  126 format(' total surface force (fx,fy,fz) ',
627  & 'and moment about the origin(mx,my,mz)')
628  write(5,*)
629  write(5,'(2x,1p,6(1x,e13.6))') (f(j),j=1,3),(xm(j),j=1,3)
630 !
631 ! center of gravity and mean normal
632 !
633  do i1=1,3
634  cg(i1)=cg(i1)/area
635  xn(i1)=xn(i1)/area
636  enddo
637  write(5,*)
638  write(5,127)
639  127 format(' center of gravity and mean normal')
640  write(5,*)
641  write(5,'(2x,1p,6(1x,e13.6))')(cg(j),j=1,3),(xn(j),j=1,3)
642 !
643 ! moment about the center of gravity
644 !
645  write(5,*)
646  write(5,129)
647  129 format(
648  & ' moment about the center of gravity(mx,my,mz)')
649  write(5,*)
650  write(5,'(2x,1p,6(1x,e13.6))')
651  & xm(1)-cg(2)*f(3)+cg(3)*f(2),
652  & xm(2)-cg(3)*f(1)+cg(1)*f(3),
653  & xm(3)-cg(1)*f(2)+cg(2)*f(1)
654 !
655 ! area, normal and shear force
656 !
657  xnormforc=f(1)*xn(1)+f(2)*xn(2)+f(3)*xn(3)
658  shearforc=sqrt((f(1)-xnormforc*xn(1))**2+
659  & (f(2)-xnormforc*xn(2))**2+
660  & (f(3)-xnormforc*xn(3))**2)
661  write(5,*)
662  write(5,128)
663  128 format(' area, ',
664  & ' normal force (+ = tension) and shear force (size)')
665  write(5,*)
666  write(5,'(2x,1p,3(1x,e13.6))') area,xnormforc,shearforc
667  endif
668 !
669  endif
670  enddo
671 !
672  return
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
df
subroutine df(x, u, uprime, rpar, nev)
Definition: subspace.f:133
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
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
materialdata_cond
subroutine materialdata_cond(imat, ntmat_, t1l, cocon, ncocon, cond)
Definition: materialdata_cond.f:19
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_dvi
subroutine materialdata_dvi(shcon, nshcon, imat, dvi, t1l, ntmat_, ithermal)
Definition: materialdata_dvi.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)