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 fuvertex(n,x,cotet,kontet,ipoeln,ieln,node,iedge,
20     &     ipoeled,ieled,iedgmid,iedtet)
21!
22!     determination of the quality of the ball of linear tetrahedron
23!     elements around a node (P.L. George)
24!
25      implicit none
26!
27      integer n,kontet(4,*),ipoeln(*),ieln(2,*),node,ine(2,6),indexe,
28     &     ielem,nodes(4),n1,n2,i,j,iedge,ipoeled(*),ieled(2,*),
29     &     iedgmid(*),iedtet(6,*)
30!
31      real*8 cotet(3,*),alpha,x(n),volume,surface(4),totsurface,
32     &     edgelength(6),radius,hmax,quality
33!
34      data ine /1,2,2,3,1,3,1,4,2,4,3,4/
35!
36      fuvertex=0.d0
37!
38!     alpha is the proporionality factor
39!
40      alpha=dsqrt(6.d0)/12.d0
41!
42      indexe=ipoeln(node)
43!
44      do
45        if(indexe.eq.0) exit
46!
47!       an element belonging to the ball of node
48!
49        ielem=ieln(1,indexe)
50!
51        do i=1,4
52          nodes(i)=kontet(i,ielem)
53          if(nodes(i).eq.node) then
54            do j=1,3
55              cotet(j,node)=x(j)
56            enddo
57          endif
58        enddo
59!
60!     calculating the volume of the element
61!
62        call calcvol(nodes(1),nodes(2),nodes(3),nodes(4),cotet,
63     &       volume)
64        if(volume.le.0.d0) volume=1.d-30
65!
66!     calculating area of each face in the element
67!
68        call calcsurf(nodes(1),nodes(2),nodes(3),cotet,
69     &       surface(1))
70        call calcsurf(nodes(2),nodes(3),nodes(4),cotet,
71     &       surface(2))
72        call calcsurf(nodes(3),nodes(4),nodes(1),cotet,
73     &       surface(3))
74        call calcsurf(nodes(4),nodes(1),nodes(2),cotet,
75     &       surface(4))
76!
77!     calculating the total surface
78!
79        totsurface=surface(1)+surface(2)+surface(3)+surface(4)
80!
81!     radius of the inscribed sphere
82!
83        radius=3.d0*volume/totsurface
84!
85!     length of each edge
86!
87        do j=1,6
88          n1=nodes(ine(1,j))
89          n2=nodes(ine(2,j))
90          edgelength(j)=dsqrt((cotet(1,n1)-cotet(1,n2))**2+
91     &         (cotet(2,n1)-cotet(2,n2))**2+
92     &         (cotet(3,n1)-cotet(3,n2))**2)
93        enddo
94!
95!     maximum edge length
96!
97        hmax=maxval(edgelength)
98!
99!     quality of the element
100!
101        quality=alpha*hmax/radius
102!
103!     worst quality of elements treated so far
104!
105        fuvertex=max(fuvertex,quality)
106!
107        indexe=ieln(2,indexe)
108      enddo
109      write(*,100) x(1),x(2),x(3),fuvertex
110 100  format('fuvertex ',4(1x,f15.8))
111!
112      return
113      end
114