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

Go to the source code of this file.

Functions/Subroutines

subroutine zienzhu (co, nk, kon, ipkon, lakon, ne, stn, ipneigh, neigh, sti, mi)
 

Function/Subroutine Documentation

◆ zienzhu()

subroutine zienzhu ( real*8, dimension(3,*)  co,
integer  nk,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
character*8, dimension(*)  lakon,
integer  ne,
real*8, dimension(6,*)  stn,
integer, dimension(*)  ipneigh,
integer, dimension(2,*)  neigh,
real*8, dimension(6,mi(1),*)  sti,
integer, dimension(*)  mi 
)
21 !
22 ! modified zienkiewicz zhu pointwise error estimator
23 !
24 ! author: Sascha Merz
25 !
26  implicit none
27 !
28  character*8 lakon(*)
29 !
30  integer maxmid
31 !
32 ! maximal number of midnodes surrounding a patch base node
33 !
34  parameter(maxmid=400)
35 !
36  integer kon(*),nk,ne,i,j,ipkon(*),indexe,nodebase,
37  & k,ipneigh(*),neigh(2,*),ifree,node,ielem,ielem1,index,index1,m,
38  & jj,mi(*),ii,ncount,nope,itypflag,inum(nk),nenei20(3,8),maxcommon,
39  & icommon,idxnode,iecount,index2,lneigh8a(7,8),ipoints,icont,
40  & nmids(maxmid),nelem(3),nelemv,iavflag,members(ne),ielem2,
41  & linpatch,iterms,inodes(4),iaddelem,ielidx,nenei10(3,4),iscount,
42  & idxnode1,maxcommon1
43 !
44  real*8 co(3,*),stn(6,*),sti(6,mi(1),*),angle,scpav(6,nk),
45  & angmax
46 !
47  data lneigh8a /2,3,4,5,6,7,8,1,3,4,5,6,7,8,1,2,4,5,6,7,8,
48  & 1,2,3,5,6,7,8,1,2,3,4,6,7,8,1,2,3,4,5,7,8,
49  & 1,2,3,4,5,6,8,1,2,3,4,5,6,7/
50 !
51  data nenei10 /5,7,8,5,6,9,6,7,10,8,9,10/
52 !
53  data nenei20 /9,12,17,9,10,18,10,11,19,11,12,20,
54  & 13,16,17,13,14,18,14,15,19,15,16,20/
55 !
56  write(*,*) 'Estimating the stress errors'
57  write(*,*)
58 !
59 ! initialization
60 !
61  ifree=0
62  do i=1,nk
63  ipneigh(i)=0
64  inum(i)=0
65  do j=1,6
66  scpav(j,i)=0.d0
67  enddo
68  enddo
69 !
70 ! build neighbour list
71 !
72  do i=1,ne
73  indexe=ipkon(i)
74  if(lakon(i)(1:4).eq.'C3D4') then
75  nope=4
76  elseif(lakon(i)(1:4).eq.'C3D6') then
77  nope=6
78  elseif(lakon(i)(1:4).eq.'C3D8') then
79  nope=8
80  elseif(lakon(i)(1:5).eq.'C3D10') then
81  nope=4
82  elseif(lakon(i)(1:5).eq.'C3D15') then
83  nope=6
84  elseif(lakon(i)(1:5).eq.'C3D20') then
85  nope=8
86  else
87  cycle
88  endif
89  do j=1,nope
90  node=kon(indexe+j)
91  ifree=ifree+1
92  neigh(2,ifree)=ipneigh(node)
93  ipneigh(node)=ifree
94  neigh(1,ifree)=i
95  enddo
96  enddo
97 !
98  patches:do nodebase=1,nk
99  if(ipneigh(nodebase).eq.0) cycle patches
100  index=ipneigh(nodebase)
101 !
102 ! initialization
103 !
104  do j=1,maxmid
105  nmids(j)=0
106  enddo
107  do j=1,3
108  nelem(j)=0
109  enddo
110  idxnode=nodebase
111  linpatch=0;ipoints=0;iavflag=0;itypflag=0
112 !
113 ! analyze neighbour structure
114 !
115  do
116  ielem=neigh(1,index)
117 !
118  if(lakon(ielem)(1:5).eq.'C3D20') then
119  nelem(1)=nelem(1)+1
120  elseif(lakon(ielem)(1:5).eq.'C3D10') then
121  nelem(2)=nelem(2)+1
122  elseif(lakon(ielem)(1:4).eq.'C3D8') then
123  nelem(3)=nelem(3)+1
124  endif
125 !
126  if(neigh(2,index).eq.0)exit
127  index=neigh(2,index)
128  enddo
129 !
130 ! itypflag 1: using hex20 estimator
131 ! itypflag 2: using tet10 estimator
132 ! itypflag 3: using hex8 estimator
133 !
134  if(nelem(1).gt.0) then
135  itypflag=1
136  elseif(nelem(1).eq.0) then
137  if(nelem(2).gt.0) then
138  itypflag=2
139  elseif(nelem(2).eq.0) then
140  if(nelem(3).gt.0) then
141  itypflag=3
142  else
143  itypflag=0
144  endif
145  endif
146  endif
147 !
148 ! case 1: element type not supported
149 !
150  if(itypflag.eq.0) then
151  write(*,*) '*WARINING in estimator: Elements of node',
152  & nodebase,' cannot be used for error estimation.'
153  do j=1,6
154  scpav(j,nodebase)=stn(j,nodebase)
155  enddo
156 !
157 !................ case 2: using hex estimator...........................
158 ! !
159 ! !
160 !
161  elseif(itypflag.eq.1.or.itypflag.eq.3) then
162 !
163 ! get spaceangle
164 !
165  call angsum(lakon,kon,ipkon,neigh,ipneigh,
166  & co,nodebase,itypflag,angle)
167 !
168 ! determine surface structure, if surface node
169 !
170  if(angle.lt.12.56535d0) then
171  call chksurf(lakon,kon,ipkon,neigh,ipneigh,co,
172  & itypflag,nodebase,icont,iscount,angmax)
173  endif
174 !
175  if(itypflag.eq.1) then
176  nelemv=nelem(1)
177  elseif(itypflag.eq.3) then
178  nelemv=nelem(3)
179  endif
180 !
181 ! if the elements neighbouring the node are not appropriate
182 ! to form a valid patch, then look for an alternative patch
183 ! base node.
184 !
185  if(
186  & angle.lt.8.64d0
187  & .and.
188  & iscount.eq.nelemv
189  & .and.
190  & angmax.lt.3.d1
191  & .or.
192  & nelemv.lt.5
193  & ) then
194 c write(*,*) 'Node',nodebase,' is not valid as',
195 c & ' a patch base node. Seeking another',
196 c & ' patch base node.'
197 !
198 ! node is on surface. seeking a node within volume having
199 ! the most elements in common
200 !
201  index=ipneigh(nodebase)
202 !
203 ! index points to the location in neigh, where the element
204 ! number is given
205 !
206  maxcommon=0
207  elements: do
208  if(index.eq.0) exit elements
209  ielem=neigh(1,index)
210 !
211  if(.not.((lakon(ielem)(1:5).eq.'C3D20'
212  & .and.itypflag.eq.1)
213  & .or.(lakon(ielem)(1:4).eq.'C3D8'
214  & .and.itypflag.eq.3))) then
215  index=neigh(2,index)
216  cycle elements
217  endif
218 !
219 ! ielem: element number
220 !
221  indexe=ipkon(ielem)
222 !
223 ! find element # corresp. 2 global node #
224 !
225  do m=1,8
226  if(kon(indexe+m).eq.nodebase) exit
227  enddo
228 !
229 ! for every corner node of a neighbouring
230 ! element count the common elements
231 !
232  nodes: do j=1,7
233  k=lneigh8a(j,m)
234 !
235 ! get global node # for neighbouring node
236 ! and count the elements in common
237 !
238  node=kon(ipkon(ielem)+k)
239 !
240 ! skip, if this node is on surface as well
241 !
242  call angsum(lakon,kon,ipkon,neigh,ipneigh,
243  & co,node,itypflag,angle)
244  if(angle.lt.12.56535d0) cycle nodes
245 !
246 ! go through neighbouring elements to find common
247 ! elements
248 !
249  icommon=0
250  index1=ipneigh(nodebase)
251  outer: do
252  if(index1.eq.0) exit outer
253  ielem1=neigh(1,index1)
254  if(.not.((lakon(ielem1)(1:5).eq.'C3D20'
255  & .and.itypflag.eq.1)
256  & .or.(lakon(ielem1)(1:4).eq.'C3D8'
257  & .and.itypflag.eq.3))) then
258  index1=neigh(2,index1)
259  cycle outer
260  endif
261  index2=ipneigh(node)
262  inner: do
263  if(index2.eq.0) exit inner
264  ielem2=neigh(1,index2)
265  if(.not.((lakon(ielem2)(1:5).eq.'C3D20'
266  & .and.itypflag.eq.1)
267  & .or.(lakon(ielem2)(1:4).eq.'C3D8'
268  & .and.itypflag.eq.3))) then
269  index2=neigh(2,index2)
270  cycle inner
271  endif
272 !
273 ! check, if element is common
274 !
275  if(ielem1.eq.ielem2) then
276  icommon=icommon+1
277  endif
278  index2=neigh(2,index2)
279  enddo inner
280  index1=neigh(2,index1)
281  enddo outer
282 !
283 ! check, if the last two nodes have the most
284 ! common elements
285 !
286  if(icommon.gt.maxcommon) then
287  maxcommon=icommon
288  idxnode=node
289  endif
290  enddo nodes
291 !
292  index=neigh(2,index)
293  enddo elements
294 !
295 ! however, if there is a node on the surface, having more
296 ! elements in common and the original base node is not on
297 ! a free surface, then use this node as basenode.
298 !
299  if(icont.eq.0) then
300  index=ipneigh(nodebase)
301 !
302  maxcommon1=0
303  elements1: do
304  if(index.eq.0) exit elements1
305  ielem=neigh(1,index)
306 !
307  if(.not.((lakon(ielem)(1:5).eq.'C3D20'
308  & .and.itypflag.eq.1)
309  & .or.(lakon(ielem)(1:4).eq.'C3D8'
310  & .and.itypflag.eq.3))) then
311  index=neigh(2,index)
312  cycle elements1
313  endif
314 !
315  indexe=ipkon(ielem)
316 !
317  do m=1,8
318  if(kon(indexe+m).eq.nodebase) exit
319  enddo
320 !
321  nodes1: do j=1,7
322  k=lneigh8a(j,m)
323 !
324  node=kon(ipkon(ielem)+k)
325 !
326  icommon=0
327  index1=ipneigh(nodebase)
328  outer1: do
329  if(index1.eq.0) exit outer1
330  ielem1=neigh(1,index1)
331  if(.not.((lakon(ielem1)(1:5).eq.'C3D20'
332  & .and.itypflag.eq.1)
333  & .or.(lakon(ielem1)(1:4).eq.'C3D8'
334  & .and.itypflag.eq.3))) then
335  index1=neigh(2,index1)
336  cycle outer1
337  endif
338  index2=ipneigh(node)
339  inner1: do
340  if(index2.eq.0) exit inner1
341  ielem2=neigh(1,index2)
342  if(.not.((lakon(ielem2)(1:5).eq.'C3D20'
343  & .and.itypflag.eq.1)
344  & .or.(lakon(ielem2)(1:4).eq.'C3D8'
345  & .and.itypflag.eq.3))) then
346  index2=neigh(2,index2)
347  cycle inner1
348  endif
349 !
350  if(ielem1.eq.ielem2) then
351  icommon=icommon+1
352  endif
353  index2=neigh(2,index2)
354  enddo inner1
355  index1=neigh(2,index1)
356  enddo outer1
357 !
358  if(icommon.gt.maxcommon1) then
359  maxcommon1=icommon
360  idxnode1=node
361  endif
362  enddo nodes1
363 !
364  index=neigh(2,index)
365  enddo elements1
366 !
367  if(maxcommon1.gt.maxcommon) then
368  idxnode=idxnode1
369  endif
370  endif
371  index=ipneigh(idxnode)
372  nelemv=0
373  do
374  if(index.eq.0) exit
375  ielem=neigh(1,index)
376  if(.not.((lakon(ielem)(1:5).eq.'C3D20'
377  & .and.itypflag.eq.1)
378  & .or.(lakon(ielem)(1:4).eq.'C3D8'
379  & .and.itypflag.eq.3))) then
380  index=neigh(2,index)
381  cycle
382  endif
383  nelemv=nelemv+1
384  index=neigh(2,index)
385  enddo
386 !
387 ! verify
388 !
389  call angsum(lakon,kon,ipkon,neigh,ipneigh,
390  & co,idxnode,itypflag,angle)
391  if(angle.lt.12.56535d0) then
392  call chksurf(lakon,kon,ipkon,neigh,ipneigh,co,
393  & itypflag,idxnode,icont,iscount,angmax)
394  endif
395  if(
396  & angle.lt.8.64d0
397  & .and.
398  & iscount.eq.nelemv
399  & .and.
400  & angmax.lt.3.d1
401  & .or.
402  & nelemv.lt.5
403  & ) then
404  write(*,*) '*WARNING in estimator: Patch not',
405  & ' appropriate for patch recovery,'
406  write(*,*) ' using',
407  & ' average of sampling point values', nodebase
408  iavflag=1
409  endif
410 !
411 c write(*,*) 'Alternative node is',idxnode
412  endif
413 !
414  index=ipneigh(idxnode)
415 !
416 ! determine patch elements (members), the number of
417 ! patch elements (linpatch) and the number of
418 ! sampling points (ipoints) in the patch
419 !
420  do
421  if(index.eq.0) exit
422  ielem=neigh(1,index)
423  if(.not.((lakon(ielem)(1:5).eq.'C3D20'
424  & .and.itypflag.eq.1)
425  & .or.(lakon(ielem)(1:4).eq.'C3D8'
426  & .and.itypflag.eq.3))) then
427  index=neigh(2,index)
428  cycle
429  endif
430  if(lakon(ielem)(1:4).eq.'C3D8') then
431  if(lakon(ielem)(5:5).eq.'R') then
432  ipoints=ipoints+1
433  else
434  ipoints=ipoints+8
435  endif
436  elseif(lakon(ielem)(1:5).eq.'C3D20') then
437  if(lakon(ielem)(6:6).eq.'R') then
438  ipoints=ipoints+8
439  else
440  ipoints=ipoints+27
441  endif
442  endif
443  linpatch=linpatch+1
444  members(linpatch)=ielem
445  index=neigh(2,index)
446  enddo
447 !
448  if(itypflag.eq.1) iterms=20
449  if(itypflag.eq.3) iterms=4
450 !
451 ! evaluate patch for patch base node (nodebase)
452 !
453  call patch(iterms,nodebase,sti,scpav,mi(1),kon,ipkon,
454  & ipoints,members,linpatch,co,lakon,iavflag)
455  inum(nodebase)=1
456 !
457 ! midnodes
458 !
459  if(itypflag.eq.1) then
460  k=1
461  hexelements: do ielidx=1,linpatch
462  ielem=members(ielidx)
463 !
464 ! find out index of base node in kon
465 !
466  do m=1,8
467  if(nodebase.eq.kon(ipkon(ielem)+m)) exit
468  if(m.eq.8) then
469  cycle hexelements
470  endif
471  enddo
472 !
473  hexnodes: do j=1,3
474 !
475 ! if node already in patch, skip
476 !
477  if(k.gt.1) then
478  do ii=1,k-1
479  if(nmids(ii)
480  & .eq.kon(ipkon(ielem)+nenei20(j,m))) then
481  cycle hexnodes
482  endif
483  enddo
484  endif
485  nmids(k)=kon(ipkon(ielem)+nenei20(j,m))
486  inum(nmids(k))=inum(nmids(k))+1
487  call patch(iterms,nmids(k),sti,scpav,mi(1),kon,
488  & ipkon,ipoints,members,linpatch,co,lakon,
489  & iavflag)
490  k=k+1
491  if(k.gt.maxmid) then
492  write(*,*) '*ERROR in estimator: array size',
493  & ' for midnodes exceeded'
494  exit hexelements
495  endif
496  enddo hexnodes
497  enddo hexelements
498  endif
499 !
500 !................ case 3: using tet estimator..........................
501 ! !
502 ! !
503  elseif(itypflag.eq.2) then
504 !
505 !
506 ! for the tetrahedral elements it could happen, that additional
507 ! elements have to be added to the patch.
508 ! the information already obtained is the number of elements
509 ! neighbouring the node.
510 ! now a check should be undertaken, if the shape of the initial
511 ! patch is useful.
512 !
513 !
514 ! case 1: the patch has a conical shape. this occurs, if
515 ! if all neighbouring elements of the
516 ! patch base node have another vertice node in common
517 ! (integration points are then not reasonably distributed).
518 ! the patch looks then like a cone or pyramid:
519 !
520 ! x
521 ! /|\
522 ! / | \
523 ! / | \
524 ! / | \
525 ! / | \
526 ! .........x_____o_____x.................
527 !
528 ! ...... surface
529 !
530 ! x patch member node
531 ! o patch base node
532 !
533 ! a alternative patch base node has to be found, if the following
534 ! conditions occur:
535 !
536 ! 1. the node is not a volume node and there is a conical
537 ! patch with a even number of patch members
538 !
539 ! 2. the node has only one neighbouring element. then it is
540 ! most likely a corner node and a good patch will be created
541 ! if another edgenode of this element is used as patch base node
542 !
543 ! 3. the node is located at a free surface and the number of
544 ! neighbouring elements is odd and the elements have conical
545 ! shape
546 !
547  call angsum(lakon,kon,ipkon,neigh,ipneigh,co,nodebase,
548  & itypflag,angle)
549 !
550 ! if the node is on the surface, then check, if the surfaces
551 ! adjacent to the node have normal vectors with an angle of
552 ! maximal 10 degree (free surface is assumed, otherwise edge
553 ! or corner)
554 !
555  if(angle.lt.12.56535d0) then
556  call chksurf(lakon,kon,ipkon,neigh,ipneigh,co,
557  & itypflag,nodebase,icont,iscount,angmax)
558  endif
559 !
560  nelemv=nelem(2)
561 !
562  icommon=0
563 !
564  if(
565  & angle.lt.12.56535d0
566  & .and.(
567  & mod(nelemv,2).eq.0
568  & .or.
569  & nelemv.eq.1
570  & )
571  & .or.(
572  & icont.eq.1
573  & .and.
574  & mod(nelemv,2).ne.0
575  & )
576  & ) then
577 !
578 ! searching for another common vertex node
579 !
580  index=ipneigh(nodebase)
581  idxnode=0
582  node=0
583  conical: do
584  if(index.eq.0) exit
585  ielem=neigh(1,index)
586  if(.not.(lakon(ielem)(1:5).eq.'C3D10'
587  & .and.itypflag.eq.2)) then
588  index=neigh(2,index)
589  cycle
590  endif
591  do j=1,4
592  node=kon(ipkon(ielem)+j)
593  if(node.eq.nodebase) cycle
594 !
595 ! is this node a member of all patch elements?
596 !
597  index1=ipneigh(nodebase)
598  ncount=0
599  do
600  if(index1.eq.0) exit
601  ielem1=neigh(1,index1)
602  if(.not.(lakon(ielem1)(1:5).eq.'C3D10'
603  & .and.itypflag.eq.2)) then
604  index1=neigh(2,index1)
605  cycle
606  endif
607  do jj=1,4
608  if(idxnode.eq.kon(ipkon(ielem1)+jj)
609  & .and.idxnode.ne.0) cycle
610  if(node.eq.kon(ipkon(ielem1)+jj))
611  & ncount=ncount+1
612  enddo
613  index1=neigh(2,index1)
614  enddo
615  if(ncount.eq.nelemv) then
616  icommon=icommon+1
617  idxnode=node
618  endif
619  enddo
620  index=neigh(2,index)
621  enddo conical
622 !
623  if(icommon.gt.0) then
624 c write(*,*) 'conical patch found, node',nodebase,
625 c & ' using node',idxnode,' as base node instead'
626 !
627 ! count elements
628 !
629  index=ipneigh(idxnode)
630  nelemv=0
631  do
632  if(index.eq.0) exit
633  ielem=neigh(1,index)
634  if(.not.(lakon(ielem)(1:5).eq.'C3D10'
635  & .and.itypflag.eq.2)) then
636  index=neigh(2,index)
637  cycle
638  endif
639  nelemv=nelemv+1
640  index=neigh(2,index)
641  enddo
642 !
643  else
644  idxnode=nodebase
645  endif
646  endif
647 !
648 ! add neighbouring elements to patch firstly
649 !
650  linpatch=0
651  index=ipneigh(idxnode)
652  do
653  if(index.eq.0) exit
654  ielem=neigh(1,index)
655  if(.not.(lakon(ielem)(1:5).eq.'C3D10'
656  & .and.itypflag.eq.2)) then
657  index=neigh(2,index)
658  cycle
659  endif
660  linpatch=linpatch+1
661  members(linpatch)=ielem
662  index=neigh(2,index)
663  enddo
664 !
665  if(nelemv.gt.0.and.nelemv.lt.4) then
666 c write(*,*) 'node',nodebase,' has',nelemv,
667 c & ' neighbouring elements'
668 !
669 ! patch has to be extended
670 !
671 ! add more elements to patch, where the elements
672 ! have to share a face with the initial patch elements
673 !
674  index=ipneigh(idxnode)
675  iecount=0
676  do
677  if(index.eq.0) exit
678  ielem=neigh(1,index)
679  if(.not.(lakon(ielem)(1:5).eq.'C3D10'
680  & .and.itypflag.eq.2)) then
681  index=neigh(2,index)
682  cycle
683  endif
684 !
685 ! check every element in the neighbour list
686 ! if there is any element in the element list
687 ! having three nodes in common (is a common
688 ! surface).
689 !
690 ! save element nodes
691 !
692  do j=1,4
693  inodes(j)=kon(ipkon(ielem)+j)
694  enddo
695 !
696 ! loop over each element in the model and
697 ! check against patch element
698 !
699  tetloop: do k=1,ne
700  if(.not.(lakon(ielem)(1:5).eq.'C3D10'
701  & .and.itypflag.eq.2)) then
702  cycle
703  endif
704 !
705 ! skip counting, if element is already in patch
706 !
707  index1=ipneigh(idxnode)
708  do
709  if(index1.eq.0) exit
710  ielem1=neigh(1,index1)
711  if(.not.(lakon(ielem1)(1:5).eq.'C3D10'
712  & .and.itypflag.eq.2)) then
713  index1=neigh(2,index1)
714  cycle
715  endif
716  if(ielem1.eq.k) cycle tetloop
717  index1=neigh(2,index1)
718  enddo
719 !
720  ncount=0
721  do j=1,4
722  do jj=1,4
723  if(inodes(jj).eq.kon(ipkon(k)+j)) then
724  ncount=ncount+1
725  endif
726  if(ncount.gt.2) then
727  linpatch=linpatch+1
728  iecount=iecount+1
729  members(linpatch)=k
730 !
731 ! keep number of 2nd found element in mind
732 !
733  if(nelemv+iecount.eq.2)
734  & iaddelem=k
735  cycle tetloop
736  endif
737  enddo
738  enddo
739  enddo tetloop
740  index=neigh(2,index)
741  enddo
742 !
743 ! if there are still just two elements, treat
744 ! those elements as initial patch and extend
745 !
746  if(nelemv+iecount.eq.2) then
747 !
748 ! search again all elements
749 !
750  do j=1,4
751  inodes(j)=kon(ipkon(iaddelem)+j)
752  enddo
753  loop3: do k=1,ne
754  if(.not.(lakon(k)(1:5).eq.'C3D10'
755  & .and.itypflag.eq.2)) then
756  cycle loop3
757  endif
758  index=ipneigh(idxnode)
759  do
760  if(index.eq.0) exit
761  ielem=neigh(1,index)
762  if(.not.(lakon(ielem)(1:5).eq.'C3D10'
763  & .and.itypflag.eq.2)) then
764  exit
765  endif
766  index=neigh(2,index)
767  enddo
768  if(ielem.eq.k.or.iaddelem.eq.k) cycle loop3
769  ncount=0
770  do j=1,4
771  do jj=1,4
772  if(inodes(jj).eq.kon(ipkon(k)+j))
773  & ncount=ncount+1
774  if(ncount.gt.2) then
775  linpatch=linpatch+1
776  iecount=iecount+1
777  members(linpatch)=k
778  cycle loop3
779  endif
780  enddo
781  enddo
782  enddo loop3
783  endif
784  nelemv=nelemv+iecount
785 c write(*,*) 'now node',nodebase,' has',nelemv,
786 c & ' elements'
787 !
788 ! verify
789 !
790  if(nelemv.lt.5) then
791  write(*,*) '*WARNING in estimator: Patch not',
792  & ' appropriate for patch recovery,'
793  write(*,*) ' using',
794  & ' average of sampling point values:', nodebase
795  iavflag=1
796  endif
797  endif
798 !
799 ! nodal stresses for patch
800 !
801 ! determine degree of patch polynomial
802 !
803  ipoints=linpatch*4
804  iterms=10
805  if(ipoints.ge.17.and.ipoints.lt.21) then
806  iterms=10
807  elseif(ipoints.ge.21.and.ipoints.lt.63) then
808  iterms=11
809  elseif(ipoints.ge.63) then
810  iterms=17
811  endif
812 !
813 c write(*,*) 'using',iterms,' terms for patch',nodebase
814 !
815 ! patch for patch base node
816 !
817  call patch(iterms,nodebase,sti,scpav,mi(1),kon,ipkon,
818  & ipoints,members,linpatch,co,lakon,iavflag)
819  inum(nodebase)=1
820 !
821 ! midnodes
822 !
823  k=1
824  tetelements: do ielidx=1,linpatch
825  ielem=members(ielidx)
826 !
827 ! find out index of base node in kon
828 !
829  do m=1,4
830  if(nodebase.eq.kon(ipkon(ielem)+m)) exit
831  if(m.eq.4) then
832  cycle tetelements
833  endif
834  enddo
835 !
836  tetnodes: do j=1,3
837 !
838 ! if node already in patch, skip
839 !
840  if(k.gt.1) then
841  do ii=1,k-1
842  if(nmids(ii)
843  & .eq.kon(ipkon(ielem)+nenei10(j,m))) then
844  cycle tetnodes
845  endif
846  enddo
847  endif
848  nmids(k)=kon(ipkon(ielem)+nenei10(j,m))
849  inum(nmids(k))=inum(nmids(k))+1
850  call patch(iterms,nmids(k),sti,scpav,mi(1),kon,
851  & ipkon,ipoints,members,linpatch,co,lakon,
852  & iavflag)
853  k=k+1
854  if(k.gt.maxmid) then
855  write(*,*) '*ERROR in estimator: array size',
856  & ' for midnodes exceeded'
857  exit tetelements
858  endif
859  enddo tetnodes
860  enddo tetelements
861  endif
862  enddo patches
863 !.......................................................................
864 !
865  do i=1,nk
866  if(inum(i).gt.0) then
867  do j=1,6
868  stn(j,i)=scpav(j,i)/inum(i)
869  enddo
870  else
871  do j=1,6
872  stn(j,i)=0.d0
873  enddo
874  endif
875  enddo
876 !
877  return
chksurf
subroutine chksurf(lakon, kon, ipkon, neigh, ipneigh, co, itypflag, node, icont, iscount, angmax)
Definition: chksurf.f:21
angsum
subroutine angsum(lakon, kon, ipkon, neigh, ipneigh, co, node, itypflag, angle)
Definition: angsum.f:21
nodes
subroutine nodes(inpc, textpart, co, nk, nk_, set, istartset, iendset, ialset, nset, nset_, nalset, nalset_, istep, istat, n, iline, ipol, inl, ipoinp, inp, ipoinpc)
Definition: nodes.f:22
patch
subroutine patch(iterms, node, sti, scpav, mi, kon, ipkon, ipoints, members, linpatch, co, lakon, iavflag)
Definition: patch.f:21
elements
subroutine elements(inpc, textpart, kon, ipkon, lakon, nkon, ne, ne_, set, istartset, iendset, ialset, nset, nset_, nalset, nalset_, mi, ixfree, iponor, xnor, istep, istat, n, iline, ipol, inl, ipoinp, inp, iaxial, ipoinpc, solid, cfd, network, filab, nlabel, out3d, iuel, nuel_)
Definition: elements.f:24
Hosted by OpenAircraft.com, (Michigan UAV, LLC)