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