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 smoothingmidnodes(cotet,ipoed,kontet,iedtet,iedgmid, 20 & ipoeled,ieled,qualityjac,iponn,inn,h,iexternedg,netet_, 21 & nktet_) 22! 23! smoothing a quadratic tetrahedral mesh by moving the midnodes 24! 25 implicit none 26! 27 integer iedge,ipoed(*),iedtet(6,*),kontet(4,*),iedgmid(*), 28 & ipoeled(*),ieled(2,*),iter,iponn(*),inn(2,*),indexe,ielemold, 29 & indexeold,i,j,k,node,ielem,index,iexternedg(*),neigh,netet_, 30 & nktet_ 31! 32 real*8 cotet(3,*),qualityjac(*),relaxw,qmax,sumalpha,cpycotet(3), 33 & haverage,alphaj,h(*),sumalphap(3) 34! 35! relaxation factor 36! 37 relaxw=0.5d0 38! 39! number of smoothing iterations 40! 41 iter=10 42! 43 do k=1,iter 44! 45! subsurface nodes 46! 47 loop1: do i=1,nktet_ 48! 49 iedge=ipoed(i) 50! 51 if(iedge.eq.0) cycle 52 if(iexternedg(iedge).ne.0) cycle 53! 54 node=iedgmid(iedge) 55! 56 if(iponn(node).eq.0) cycle 57! 58! calculate the quality of the shell (all elements containing 59! the node), it is the worst quality of all elements belonging 60! to the shell (a shell is the set of all elements containing 61! a specific edge 62! 63 qmax=0.d0 64 indexe=ipoeled(iedge) 65! 66 do 67 if(indexe.eq.0) exit 68 ielem=ieled(1,indexe) 69! 70 if(qualityjac(ielem).gt.qmax) qmax=qualityjac(ielem) 71 indexe=ieled(2,indexe) 72 enddo 73! 74! initializing the numerator (sumalphap) and denominator 75! (sumalpha) of the new position 76! 77 sumalpha=0.d0 78 do j=1,3 79 sumalphap(j)=0.d0 80 enddo 81! 82! calculating the new position: loop over all neighboring nodes 83! 84 index=iponn(node) 85 do 86 if(index.eq.0) exit 87! 88! neighboring node 89! 90 neigh=inn(1,index) 91 haverage=(h(node)+h(neigh))/2.d0 92! 93! weight = inverse of the square of the desired edge length 94! 95 alphaj=1.d0/(haverage*haverage) 96 sumalpha=sumalpha+alphaj 97 do j=1,3 98 sumalphap(j)=sumalphap(j)+alphaj*cotet(j,neigh) 99 enddo 100 index=inn(2,index) 101 enddo 102! 103! storing the old position of the node; determining the new 104! position 105! 106 do j=1,3 107 cpycotet(j)=cotet(j,node) 108 cotet(j,node)=(1.d0-relaxw)*cotet(j,node)+ 109 & relaxw*sumalphap(j)/sumalpha 110 enddo 111! 112! check the quality of the tetrahedral elements in the shell for 113! the new position of node; as soon as a tetrahedron is detected 114! with a quality exceeding qmax the coordinates of node are reverted 115! to the ones before smoothing 116! 117 indexe=ipoeled(iedge) 118 do 119 if(indexe.eq.0) exit 120 ielem=ieled(1,indexe) 121! 122 call quadmeshquality(netet_,cotet,kontet,iedtet, 123 & iedgmid,qualityjac,ielem) 124! 125 if(qualityjac(ielem).gt.qmax) then 126! 127! 1. restore the original coordinates 128! 2. recalculate the quality of the tetrahedra in the shell 129! covered so far in the current loop 130! 131 do j=1,3 132 cotet(j,node)=cpycotet(j) 133 enddo 134 indexeold=ipoeled(iedge) 135 do 136 ielemold=ieled(1,indexeold) 137! 138 call quadmeshquality(netet_,cotet,kontet,iedtet, 139 & iedgmid,qualityjac,ielemold) 140 if(ielemold.eq.ielem) exit 141 indexeold=ieled(2,indexeold) 142c write(*,*) '*INFO in smoothingmidnodes: reverting to' 143c write(*,*) ' old coordinate for node ',node 144c write(*,*) ' and element ',ielem 145! 146 enddo 147 cycle loop1 148 endif 149 indexe=ieled(2,indexe) 150 enddo 151 enddo loop1 152! 153 enddo 154! 155 return 156 end 157