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

Go to the source code of this file.

Functions/Subroutines

subroutine patch (iterms, node, sti, scpav, mi, kon, ipkon, ipoints, members, linpatch, co, lakon, iavflag)
 

Function/Subroutine Documentation

◆ patch()

subroutine patch ( integer  iterms,
integer  node,
real*8, dimension(6,mi(1),*)  sti,
real*8, dimension(6,*)  scpav,
integer, dimension(*)  mi,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
integer  ipoints,
integer, dimension(*)  members,
integer  linpatch,
real*8, dimension(3,*)  co,
character*8, dimension(*)  lakon,
integer  iavflag 
)
21 !
22 ! computes the smoothed nodal stresses for an element patch
23 !
24 ! author: Sascha Merz
25 !
26  implicit none
27 !
28  integer i,j,nope,mint3d,indexe,ielem,kon(*),ipkon(*),
29  & iterms,k,iflag,nrhs,info,node,ipnt,
30  & mi(*),irow,ipoints,members(*),linpatch,ielidx,iavflag
31 !
32  real*8 xl(3,20),co(3,*),shp(4,20),
33  & pgauss(3),pol(20),pntdist,w,scpav(6,*),xsj,
34  & sti(6,mi(1),*),xi,et,ze,rv1(ipoints),pp(ipoints,iterms),
35  & pdat(ipoints,6),pfit(6),pwrk(iterms),pre(ipoints,iterms),
36  & z(ipoints,ipoints)
37 !
38  real*8 tmpstr(6),gauss3d5e(3,4)
39 !
40  character*8 lakon(*)
41 !
42  logical matu,matv
43 !
44  include 'gauss.f'
45 !
46 ! initialize
47 !
48 ! iavflag: if 1, then build average of integration point values
49 !
50  if(iavflag.eq.0) then
51  do i=1,6
52  pfit(i)=0.d0
53  enddo
54  iflag=1
55 !
56 ! irow: row number of the rectangular matrix of the overdetermined
57 ! system of equations
58 !
59  irow=0
60 !
61 ! loop over patch elements
62 ! linpatch: number of elements in patch
63 !
64  do ielidx=1,linpatch
65 !
66  ielem=members(ielidx)
67  if(lakon(ielem)(1:5).eq.'C3D20') then
68 !
69 ! nope: nodes per element
70 !
71  nope=20
72  if(lakon(ielem)(6:6).eq.'R') then
73  mint3d=8
74  else
75  mint3d=27
76  endif
77  elseif(lakon(ielem)(1:4).eq.'C3D8') then
78  nope=8
79  if(lakon(ielem)(5:5).eq.'R') then
80  mint3d=1
81  else
82  mint3d=8
83  endif
84  elseif(lakon(ielem)(1:5).eq.'C3D10') then
85  nope=10
86  mint3d=4
87  endif
88 !
89  indexe=ipkon(ielem)
90 !
91 ! coordinates of the nodes belonging to the element
92 !
93  do j=1,nope
94  do k=1,3
95  xl(k,j)=co(k,kon(indexe+j))
96  enddo
97  enddo
98 !
99 ! loop over the integration points in one element
100 !
101  do ipnt=1,mint3d
102  irow=irow+1
103  if(nope.eq.10) then
104  xi=gauss3d5(1,ipnt)
105  et=gauss3d5(2,ipnt)
106  ze=gauss3d5(3,ipnt)
107  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
108  elseif(nope.eq.20) then
109  if(mint3d.eq.8) then
110  xi=gauss3d2(1,ipnt)
111  et=gauss3d2(2,ipnt)
112  ze=gauss3d2(3,ipnt)
113  elseif(mint3d.eq.27) then
114  xi=gauss3d3(1,ipnt)
115  et=gauss3d3(2,ipnt)
116  ze=gauss3d3(3,ipnt)
117  endif
118  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
119  elseif(nope.eq.8) then
120  if(mint3d.eq.1) then
121  xi=gauss3d1(1,ipnt)
122  et=gauss3d1(2,ipnt)
123  ze=gauss3d1(3,ipnt)
124  elseif(mint3d.eq.8) then
125  xi=gauss3d2(1,ipnt)
126  et=gauss3d2(2,ipnt)
127  ze=gauss3d2(3,ipnt)
128  endif
129  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
130  endif
131 !
132 ! in pgauss the global coordinates of the integration
133 ! point are saved.
134 !
135  do j=1,3
136  pgauss(j)=0.d0
137  do k=1,nope
138  pgauss(j)=pgauss(j)+shp(4,k)*xl(j,k)
139  enddo
140 !
141 ! the origin of the coordinate system is moved to the
142 ! evaluated node for higher numerical stability
143 !
144  pgauss(j)=pgauss(j)-co(j,node)
145  enddo
146 !
147 ! evaluate patch polynomial for the integration point
148 !
149  pol(1)=1.d0
150  pol(2)=pgauss(1)
151  pol(3)=pgauss(2)
152  pol(4)=pgauss(3)
153  pol(5)=pgauss(1)*pgauss(2)
154  pol(6)=pgauss(1)*pgauss(3)
155  pol(7)=pgauss(2)*pgauss(3)
156  pol(8)=pgauss(1)*pgauss(1)
157  pol(9)=pgauss(2)*pgauss(2)
158  pol(10)=pgauss(3)*pgauss(3)
159  pol(11)=pgauss(1)*pgauss(2)*pgauss(3)
160  if(iterms.gt.11) then
161  pol(12)=pgauss(1)*pgauss(1)*pgauss(2)
162  pol(13)=pgauss(1)*pgauss(2)*pgauss(2)
163  pol(14)=pgauss(1)*pgauss(1)*pgauss(3)
164  pol(15)=pgauss(1)*pgauss(3)*pgauss(3)
165  pol(16)=pgauss(2)*pgauss(2)*pgauss(3)
166  pol(17)=pgauss(2)*pgauss(3)*pgauss(3)
167  if(iterms.gt.17) then
168  pol(18)=pgauss(1)*pgauss(1)*pgauss(1)
169  pol(19)=pgauss(2)*pgauss(2)*pgauss(2)
170  pol(20)=pgauss(3)*pgauss(3)*pgauss(3)
171  endif
172  endif
173 !
174 ! weighting for integration point
175 !
176  pntdist=dsqrt(pgauss(1)*pgauss(1)
177  & +pgauss(2)*pgauss(2)
178  & +pgauss(3)*pgauss(3))
179  w=pntdist**-1.5d0
180 !
181  do j=1,6
182  pdat(irow,j)=sti(j,ipnt,ielem)*w
183  enddo
184 !
185  do j=1,iterms
186  pp(irow,j)=pol(j)*w
187  enddo
188  enddo
189  enddo
190 !
191 ! using singular value decomposition for the least squares fit
192 !
193  matu=.false.
194  matv=.true.
195  nrhs=6
196 !
197  call hybsvd(ipoints,ipoints,ipoints,ipoints,ipoints,
198  & ipoints,iterms,pp,pwrk,matu,pp,matv,
199  & pre,z,pdat,nrhs,info,rv1)
200  if(info.ne.0) then
201  write(*,*) '*ERROR in patch: Bad conditioned matrix,',
202  & ' using average of sampling point values.'
203  iavflag=1
204  endif
205  endif
206 !
207 ! matrix multiplication. only the first value of the
208 ! solution vector is needed. the singular values are manipulated
209 ! to increase the numerical stability
210 !
211  if(iavflag.eq.0) then
212  do j=1,iterms
213  if(pwrk(j).lt.1.e-22) then
214  pwrk(j)=0.d0
215  else
216  pwrk(j)=1.d0/pwrk(j)
217  endif
218  enddo
219  do j=1,iterms
220  pre(1,j)=pre(1,j)*pwrk(j)
221  enddo
222  do j=1,nrhs
223  do k=1,iterms
224  pfit(j)=pfit(j)+pre(1,k)*pdat(k,j)
225  enddo
226  enddo
227 !
228 ! pfit is an array containing the coefficients for the polynom
229 ! for the six stress components
230 !
231 ! solution in the node
232 !
233  do j=1,6
234  scpav(j,node)=scpav(j,node)+pfit(j)
235  enddo
236  endif
237 !
238 ! if there are not enough elements to fit a polynomial,
239 ! build average value of the sampling point values
240 !
241  if(iavflag.eq.1) then
242  do j=1,6
243  tmpstr(j)=0.d0
244  enddo
245  irow=0
246  do ielidx=1,linpatch
247  ielem=members(ielidx)
248  if(lakon(ielem)(1:5).eq.'C3D20') then
249  if(lakon(ielem)(6:6).eq.'R') then
250  mint3d=8
251  else
252  mint3d=27
253  endif
254  elseif(lakon(ielem)(1:5).eq.'C3D10') then
255  mint3d=4
256  elseif(lakon(ielem)(1:4).eq.'C3D8') then
257  if(lakon(ielem)(5:5).eq.'R') then
258  mint3d=1
259  else
260  mint3d=8
261  endif
262  endif
263  do ipnt=1,mint3d
264  irow=irow+1
265  do j=1,6
266  tmpstr(j)=tmpstr(j)+sti(j,ipnt,ielem)
267  enddo
268  enddo
269  enddo
270  do j=1,6
271  scpav(j,node)=scpav(j,node)+tmpstr(j)/irow
272  enddo
273  endif
274 !
275  return
hybsvd
subroutine hybsvd(NA, NU, NV, NZ, NB, M, N, A, W, MATU, U, MATV, V, Z, B, IRHS, IERR, RV1)
Definition: hybsvd.f:3
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
shape20h
subroutine shape20h(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape20h.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)