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