1!
2!     CalculiX - A 3-dimensional finite element program
3!              Copyright (C) 1998-2021 Guido Dhondt
4!
5!     This program is free software; you can redistribute it and/or
6!     modify it under the terms of the GNU General Public License as
7!     published by the Free Software Foundation(version 2);
8!
9!
10!     This program is distributed in the hope that it will be useful,
11!     but WITHOUT ANY WARRANTY; without even the implied warranty of
12!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!     GNU General Public License for more details.
14!
15!     You should have received a copy of the GNU General Public License
16!     along with this program; if not, write to the Free Software
17!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18!
19      subroutine zienzhu(co,nk,kon,ipkon,lakon,ne,stn,
20     & ipneigh,neigh,sti,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
194c               write(*,*) 'Node',nodebase,' is not valid as',
195c     &              ' a patch base node. Seeking another',
196c     &              ' 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!
411c               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
624c                  write(*,*) 'conical patch found, node',nodebase,
625c     &                 ' 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
666c               write(*,*) 'node',nodebase,' has',nelemv,
667c     &              ' 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
785c               write(*,*) 'now node',nodebase,' has',nelemv,
786c     &              ' 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!
813c            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
878      end
879