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

Go to the source code of this file.

Functions/Subroutines

subroutine errorestimator (yi, yn, ipkon, kon, lakon, nk, ne, mi, ielmat, nterms, inum, co, vold, cflag)
 

Function/Subroutine Documentation

◆ errorestimator()

subroutine errorestimator ( real*8, dimension(nterms,mi(1),*)  yi,
real*8, dimension(nterms,*)  yn,
integer, dimension(*)  ipkon,
integer, dimension(*)  kon,
character*8, dimension(*)  lakon,
integer  nk,
integer  ne,
integer, dimension(*)  mi,
integer, dimension(mi(3),*)  ielmat,
integer  nterms,
integer, dimension(*)  inum,
real*8, dimension(3,*)  co,
real*8, dimension(0:mi(2),*)  vold,
character*1  cflag 
)
21 !
22 ! the error in the node is calculated based on the maximum difference
23 ! between the max principal stress (mechanical calculations)
24 ! or heat flux size (thermal calculations) at the integration
25 ! points in the elements belonging to the node
26 !
27 ! for mechanical applications this error is converted into a
28 ! true stress error with heuristic equations
29 !
30  implicit none
31 !
32  logical force
33 !
34  character*1 cflag
35  character*8 lakon(*),lakonl
36 !
37  integer ipkon(*),kon(*),mi(*),ne,indexe,null,nonei20(3,12),
38  & nonei10(3,6),nk,i,j,k,node,nonei15(3,9),nopev,nterms,
39  & mint3d,ielmat(mi(3),*),inum(*)
40 !
41  real*8 yi(nterms,mi(1),*),yn(nterms,*),size,wpsmin,wpsmax,
42  & absdiff,reldiff,sizemax,al(3),sizemin,c(3,3),
43  & wpsmin1,wpsmax1,wpsmin3,wpsmax3,co(3,*),vold(0:mi(2),*)
44 !
45  data nonei10 /5,1,2,6,2,3,7,3,1,8,1,4,9,2,4,10,3,4/
46 !
47  data nonei15 /7,1,2,8,2,3,9,3,1,10,4,5,11,5,6,12,6,4,
48  & 13,1,4,14,2,5,15,3,6/
49 !
50  data nonei20 /9,1,2,10,2,3,11,3,4,12,4,1,
51  & 13,5,6,14,6,7,15,7,8,16,8,5,
52  & 17,1,5,18,2,6,19,3,7,20,4,8/
53 !
54  null=0
55 !
56  do i=1,nk
57  do j=1,nterms
58  yn(j,i)=0.d0
59  enddo
60  enddo
61 !
62  do i=1,ne
63 !
64  if(ipkon(i).lt.0) cycle
65  indexe=ipkon(i)
66  lakonl=lakon(i)
67 !
68  if(lakonl(7:8).eq.'LC') cycle
69 !
70  if(lakonl(1:1).eq.'F') then
71  cycle
72  elseif(lakonl(4:4).eq.'2') then
73  nopev=8
74  elseif(lakonl(4:4).eq.'8') then
75  nopev=8
76  elseif(lakonl(4:5).eq.'10') then
77  nopev=4
78  elseif(lakonl(4:4).eq.'4') then
79  nopev=4
80  elseif(lakonl(4:5).eq.'15') then
81  nopev=6
82  elseif(lakonl(4:4).eq.'6') then
83  nopev=6
84  else
85  cycle
86  endif
87 !
88  if(lakonl(4:5).eq.'8R') then
89  mint3d=1
90  elseif((lakonl(4:4).eq.'8').or.
91  & (lakonl(4:6).eq.'20R')) then
92  mint3d=8
93  elseif(lakonl(4:4).eq.'2') then
94  mint3d=27
95  elseif(lakonl(4:5).eq.'10') then
96  mint3d=4
97  elseif(lakonl(4:4).eq.'4') then
98  mint3d=1
99  elseif(lakonl(4:5).eq.'15') then
100  mint3d=9
101  elseif(lakonl(4:5).eq.'6') then
102  mint3d=2
103  elseif(lakonl(1:2).eq.'ES') then
104  cycle
105  endif
106 !
107 ! calculating the maximal differences of first principal
108 ! stress or heat flux size across the integration points
109 !
110  absdiff=0.d0
111  reldiff=0.d0
112 !
113  if(nterms.eq.6) then
114 !
115 ! mechanical calculation: max principal stress
116 !
117  wpsmin1=1.d30
118  wpsmax1=-1.e30
119 !
120  wpsmin3=1.d30
121  wpsmax3=-1.e30
122 !
123  do j=1,mint3d
124  c(1,1)=yi(1,j,i)
125  c(2,2)=yi(2,j,i)
126  c(3,3)=yi(3,j,i)
127  c(1,2)=yi(4,j,i)
128  c(1,3)=yi(5,j,i)
129  c(2,3)=yi(6,j,i)
130 !
131 ! calculate the eigenvalues
132 !
133 ! al(1): smallest eigenvalue
134 ! al(3): largest eigenvalue
135 !
136  call calceigenvalues(c,al)
137 !
138  wpsmin1=min(wpsmin1,al(1))
139  wpsmax1=max(wpsmax1,al(1))
140  wpsmin3=min(wpsmin3,al(3))
141  wpsmax3=max(wpsmax3,al(3))
142 !
143  enddo
144 !
145 ! check which eigenvalue is the largest in
146 ! absolute value
147 !
148  if(wpsmax3.ge.-wpsmin1) then
149  wpsmin=wpsmin3
150  wpsmax=wpsmax3
151  else
152  wpsmin=wpsmin1
153  wpsmax=wpsmax1
154  endif
155 !
156  absdiff=wpsmax-wpsmin
157  if(max(dabs(wpsmax),dabs(wpsmin)).lt.1.d-30) then
158  reldiff=0.d0
159  else
160  reldiff=absdiff/(max(dabs(wpsmax),dabs(wpsmin)))
161  endif
162  else
163 !
164 ! thermal calculation: size of heat flux
165 !
166  sizemin=1.e30
167  sizemax=0.d0
168 !
169  do j=1,mint3d
170  c(1,1)=yi(1,j,i)
171  c(2,2)=yi(2,j,i)
172  c(3,3)=yi(3,j,i)
173 !
174  size=dsqrt(c(1,1)**2+c(2,2)**2+c(3,3)**2)
175  sizemin=min(sizemin,size)
176  sizemax=max(sizemax,size)
177 !
178  enddo
179  absdiff=sizemax-sizemin
180  if(max(sizemax,sizemin).lt.1.d-30) then
181  reldiff=0.d0
182  else
183  reldiff=absdiff/(max(dabs(sizemax),dabs(sizemin)))
184  endif
185  endif
186 !
187 ! transferring the maximum to the nodes belonging to the
188 ! element
189 !
190  do j=1,nopev
191  yn(2,kon(indexe+j))=max(yn(2,kon(indexe+j)),reldiff)
192  enddo
193 !
194  enddo
195 !
196 ! converting the error estimator into a stress error (%)
197 ! through heuristic relationships
198 !
199 ! not covered: C3D8R* and C3D6* elements
200 !
201  if(nterms.eq.6) then
202  do i=1,ne
203 !
204  if(ipkon(i).lt.0) cycle
205  indexe=ipkon(i)
206  lakonl=lakon(i)
207 !
208  if(lakonl(7:8).eq.'LC') cycle
209 !
210  if(lakonl(1:1).eq.'F') then
211  cycle
212  elseif(lakonl(4:4).eq.'2') then
213  nopev=8
214  elseif(lakonl(4:4).eq.'8') then
215  nopev=8
216  elseif(lakonl(4:5).eq.'10') then
217  nopev=4
218  elseif(lakonl(4:4).eq.'4') then
219  nopev=4
220  elseif(lakonl(4:5).eq.'15') then
221  nopev=6
222  elseif(lakonl(4:4).eq.'6') then
223  nopev=6
224  else
225  cycle
226  endif
227 !
228  if(lakonl(4:5).eq.'10') then
229 !
230 ! 10-node tetrahedral element
231 !
232  do j=1,nopev
233  node=kon(indexe+j)
234  yn(1,node)=max(yn(1,node),30.000d0*yn(2,node))
235  enddo
236  elseif((lakonl(4:7).eq.'20 ').or.
237  & (lakonl(4:7).eq.'20 L').or.
238  & (lakonl(4:7).eq.'20 B')) then
239 !
240 ! true 20-node brick element or S8 or B32 (shell/beam)
241 !
242  do j=1,nopev
243  node=kon(indexe+j)
244  if(yn(2,node).le.0.26d0) then
245  yn(1,node)=max(yn(1,node),5.115d0*yn(2,node))
246  else
247  yn(1,node)=max(yn(1,node),
248  & 27.792d0*yn(2,node)-5.895d0)
249  endif
250  enddo
251  elseif(lakonl(4:6).eq.'20 ') then
252 !
253 ! expanded 20-node brick element (plane stress,
254 ! plane strain or axisymmetric)
255 !
256  do j=1,nopev
257  node=kon(indexe+j)
258  if(yn(2,node).le.0.325d0) then
259  yn(1,node)=max(yn(1,node),9.538d0*yn(2,node))
260  else
261  yn(1,node)=
262  & max(yn(1,node),53.695d0*yn(2,node)-14.351d0)
263  endif
264  enddo
265  elseif((lakonl(4:7).eq.'20R ').or.
266  & (lakonl(4:7).eq.'20RL').or.
267  & (lakonl(4:7).eq.'20RB')) then
268 !
269 ! true 20-node brick element with reduced integration
270 ! or S8R or B32R (shell/beam)
271 !
272  do j=1,nopev
273  node=kon(indexe+j)
274  if(yn(2,node).le.0.18d0) then
275  yn(1,node)=max(yn(1,node),20.278d0*yn(2,node))
276  else
277  yn(1,node)=max(yn(1,node),
278  & 74.318d0*yn(2,node)-9.727d0)
279  endif
280  enddo
281  elseif(lakonl(4:6).eq.'20R') then
282 !
283 ! expanded 20-node brick element with reduced integration
284 ! (plane stress, plane strain or axisymmetric)
285 !
286  do j=1,nopev
287  node=kon(indexe+j)
288  yn(1,node)=max(yn(1,node),54.054d0*yn(2,node))
289  enddo
290  elseif(lakonl(4:5).eq.'8I') then
291 !
292 ! true C3D8I-element or S4 or BE31 (shell/beam)
293 !
294  do j=1,nopev
295  node=kon(indexe+j)
296  if(yn(2,node).le.0.165d0) then
297  yn(1,node)=max(yn(1,node),30.303d0*yn(2,node))
298  else
299  yn(1,node)=
300  & max(yn(1,node),139.535d0*yn(2,node)-18.023d0)
301  endif
302  enddo
303  elseif(lakonl(4:7).eq.'8 ') then
304 !
305 ! true 8-node brick element
306 !
307  do j=1,nopev
308  node=kon(indexe+j)
309  if(yn(2,node).le.0.157d0) then
310  yn(1,node)=max(yn(1,node),31.847d0*yn(2,node))
311  else
312  yn(1,node)=max(yn(1,node),
313  & 85.324d0*yn(2,node)-8.396d0)
314  endif
315  enddo
316  elseif(lakonl(4:5).eq.'8 ') then
317 !
318 ! expanded 8-node brick element (plane stress, plane strain
319 ! or axisymmetric)
320 !
321  do j=1,nopev
322  node=kon(indexe+j)
323  yn(1,node)=max(yn(1,node),74.074d0*yn(2,node))
324  enddo
325  elseif(lakonl(4:5).eq.'15') then
326 !
327 ! true or expanded 15-node wedge element
328 !
329  do j=1,nopev
330  node=kon(indexe+j)
331  yn(1,node)=max(yn(1,node),46.189d0*yn(2,node))
332  enddo
333  endif
334  enddo
335  endif
336 !
337 ! determining the field values in the midside nodes
338 !
339  do i=1,ne
340 !
341  if(ipkon(i).lt.0) cycle
342  indexe=ipkon(i)
343  lakonl=lakon(i)
344 !
345  if(lakonl(7:8).eq.'LC') cycle
346 !
347  if(lakonl(4:5).eq.'20') then
348  do j=9,20
349  do k=1,2
350  yn(k,kon(indexe+j))=(
351  & yn(k,kon(indexe+nonei20(2,j-8)))+
352  & yn(k,kon(indexe+nonei20(3,j-8))))/2.d0
353  enddo
354  enddo
355  elseif(lakonl(4:5).eq.'10') then
356  do j=5,10
357  do k=1,2
358  yn(k,kon(indexe+j))=(yn(k,kon(indexe+nonei10(2,j-4)))+
359  & yn(k,kon(indexe+nonei10(3,j-4))))/2.d0
360  enddo
361  enddo
362  elseif(lakonl(4:5).eq.'15') then
363  do j=7,15
364  do k=1,2
365  yn(k,kon(indexe+j))=(
366  & yn(k,kon(indexe+nonei15(2,j-6)))+
367  & yn(k,kon(indexe+nonei15(3,j-6))))/2.d0
368  enddo
369  enddo
370  endif
371  enddo
372 !
373 ! mapping 3D on 1D/2D
374 !
375  if(cflag.eq.'I') then
376  force=.false.
377  call map3dto1d2d(yn,ipkon,inum,kon,lakon,nterms,nk,
378  & ne,cflag,co,vold,force,mi)
379  endif
380 !
381  return
max
#define max(a, b)
Definition: cascade.c:32
min
#define min(a, b)
Definition: cascade.c:31
calceigenvalues
subroutine calceigenvalues(c, al)
Definition: calceigenvalues.f:20
map3dto1d2d
subroutine map3dto1d2d(yn, ipkon, inum, kon, lakon, nfield, nk, ne, cflag, co, vold, force, mi)
Definition: map3dto1d2d.f:21
Hosted by OpenAircraft.com, (Michigan UAV, LLC)