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      real*8 function fumid(n,x,cotet,kontet,ipoeln,ieln,node,iedge,
20     &     ipoeled,ieled,iedgmid,iedtet)
21!
22!     determination of the quality of the shell of quadratic tetrahedron
23!     elements around a node; consists of a combination of:
24!     - the quality of the linear tetrahedrons in which it can be
25!       subdivided (P.L. George)
26!     - a penalty function for negative Jacobians at the integration
27!       points
28!
29      implicit none
30!
31      integer n,kontet(4,*),ipoeln(*),ieln(2,*),node,ine(2,6),indexe,
32     &     ielem,nodes(10),n1,n2,i,j,k,m,iedge,ipoeled(*),ieled(2,*),
33     &     iedgmid(*),iedtet(6,*),i6(4,8),nodesedge(4),iflag
34!
35      real*8 cotet(3,*),alpha,x(n),volume,surface(4),totsurface,
36     &     edgelength(6),radius,hmax,quality,xi,et,ze,shp(4,10),xsj,
37     &     xl(3,10)
38!
39!     nodes belonging to the linear sub-tetrahedra
40!
41      data i6 /8,9,10,4,1,5,7,8,7,6,3,10,9,8,10,7,
42     &     8,9,5,7,9,10,6,7,5,6,7,9,5,2,6,9/
43!
44!     edges of a tetrahedron
45!
46      data ine /1,2,2,3,1,3,1,4,2,4,3,4/
47!
48      data iflag /2/
49!
50      include "gauss.f"
51!
52      fumid=0.d0
53!
54!     alpha is the proporionality factor
55!
56      alpha=dsqrt(6.d0)/12.d0
57!
58      indexe=ipoeled(iedge)
59!
60      do
61        if(indexe.eq.0) exit
62        ielem=ieled(1,indexe)
63!
64!     determine the coordinates of the nodes belonging to the element
65!
66!     vertex nodes
67!
68        do i=1,4
69          nodes(i)=kontet(i,ielem)
70          do k=1,3
71            xl(k,i)=cotet(k,nodes(i))
72          enddo
73        enddo
74!
75!     middle nodes
76!
77        do i=1,6
78          nodes(i+4)=iedgmid(iedtet(i,ielem))
79          do k=1,3
80            xl(k,i+4)=cotet(k,nodes(i+4))
81          enddo
82!
83!     replace the coordinates by the optimization variables
84!
85          if(nodes(i+4).eq.node) then
86            do m=1,3
87              cotet(m,node)=x(m)
88              xl(m,i+4)=x(m)
89            enddo
90          endif
91        enddo
92!
93!     loop over the linear sub-tetrahedrons
94!
95        do i=1,8
96!
97          call calcvol(nodes(i6(1,i)),nodes(i6(2,i)),nodes(i6(3,i)),
98     &         nodes(i6(4,i)),cotet,volume)
99          if(volume.le.1.d-15) volume=1.d-15
100!
101!     calculating the area of each face of the element
102!
103          call calcsurf(nodes(i6(1,i)),nodes(i6(2,i)),nodes(i6(3,i)),
104     &         cotet,surface(1))
105          call calcsurf(nodes(i6(2,i)),nodes(i6(3,i)),nodes(i6(4,i)),
106     &         cotet,surface(2))
107          call calcsurf(nodes(i6(3,i)),nodes(i6(4,i)),nodes(i6(1,i)),
108     &         cotet,surface(3))
109          call calcsurf(nodes(i6(4,i)),nodes(i6(1,i)),nodes(i6(2,i)),
110     &         cotet,surface(4))
111!
112!     calculating the total surface
113!
114          totsurface=surface(1)+surface(2)+surface(3)+surface(4)
115!
116!     radius of the inscribed sphere
117!
118          radius=3.d0*volume/totsurface
119!
120!     length of each edge
121!
122          nodesedge(1)=nodes(i6(1,i))
123          nodesedge(2)=nodes(i6(2,i))
124          nodesedge(3)=nodes(i6(3,i))
125          nodesedge(4)=nodes(i6(4,i))
126!
127          do j=1,6
128            n1=nodesedge(ine(1,j))
129            n2=nodesedge(ine(2,j))
130            edgelength(j)=dsqrt((cotet(1,n1)-cotet(1,n2))**2+
131     &           (cotet(2,n1)-cotet(2,n2))**2+
132     &           (cotet(3,n1)-cotet(3,n2))**2)
133          enddo
134!
135!     maximum edge length
136!
137          hmax=maxval(edgelength)
138!
139!     quality of the element
140!
141          quality=alpha*hmax/radius
142!
143          fumid=max(fumid,quality)
144        enddo
145!
146!     determine the Jacobian in each integration point
147!     notice: the penalty due to a negative Jacobian must be
148!     bigger than the size of quality (volume >= 1.d-15)
149!
150        do j=1,4
151          xi=gauss3d5(1,j)
152          et=gauss3d5(2,j)
153          ze=gauss3d5(3,j)
154          call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
155          if(xsj.le.0.d0) then
156            fumid=fumid+1.d30*(1.d0-xsj)
157            exit
158          endif
159        enddo
160!
161        indexe=ieled(2,indexe)
162      enddo
163!
164      return
165      end
166