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 meshquality(netet_,kontet,cotet,quality,ielem)
20!
21!     calculate the element quality. The measure used is proportional
22!     to the ratio of the longest edge divided by the radius of the
23!     inscribed sphere. The proportionality constant is such that
24!     the quality is 1 for an equilateral tetrahedron. For all other
25!     elements it exceeds 1. The bigger this number, the worse the
26!     quality
27!
28!     if ielem>0 the quality for this element is calculated, if
29!     ielem=0 the quality for all elements in the mesh is determined
30!
31      implicit none
32!
33      integer netet_,kontet(4,*),ielem,i,j,nodes(4),n1,n2,ine(2,6)
34!
35      real*8 cotet(3,*),quality(*),volume,surface(4),totsurface,alpha,
36     &     edgelength(6),hmax,radius
37!
38      data ine /1,2,2,3,1,3,1,4,2,4,3,4/
39!
40!     alpha is the proporionality factor
41!
42      alpha=dsqrt(6.d0)/12.d0
43!
44      if(ielem.eq.0) then
45!
46!     calculating the quality for all elements
47!
48        do i=1,netet_
49          if(kontet(1,i).ne.0) then
50            do j=1,4
51              nodes(j)=kontet(j,i)
52            enddo
53!
54!     calculating the volume of the element
55!
56            call calcvol(nodes(1),nodes(2),nodes(3),nodes(4),cotet,
57     &           volume)
58            if(volume.le.0.d0) volume=1.d-30
59!
60!     calculating area of each face in the element
61!
62            call calcsurf(nodes(1),nodes(2),nodes(3),cotet,
63     &           surface(1))
64            call calcsurf(nodes(2),nodes(3),nodes(4),cotet,
65     &           surface(2))
66            call calcsurf(nodes(3),nodes(4),nodes(1),cotet,
67     &           surface(3))
68            call calcsurf(nodes(4),nodes(1),nodes(2),cotet,
69     &           surface(4))
70!
71!     calculating the total surface
72!
73            totsurface=surface(1)+surface(2)+surface(3)+surface(4)
74!
75!     radius of the inscribed sphere
76!
77            radius=3.d0*volume/totsurface
78!
79!     length of each edge
80!
81            do j=1,6
82              n1=nodes(ine(1,j))
83              n2=nodes(ine(2,j))
84              edgelength(j)=dsqrt((cotet(1,n1)-cotet(1,n2))**2+
85     &             (cotet(2,n1)-cotet(2,n2))**2+
86     &             (cotet(3,n1)-cotet(3,n2))**2)
87            enddo
88!
89!     maximum edge length
90!
91            hmax=maxval(edgelength)
92!
93!     quality
94!
95            quality(i)=alpha*hmax/radius
96!
97!     are the next lines really needed?
98!
99c            if(quality(i).lt.1.d0) then
100c              quality(i)=1.d0/quality(i)
101c           endif
102          endif
103        enddo
104      else
105!
106!     calculating the quality for just one element
107!
108        i=ielem
109!
110        if(kontet(1,i).ne.0) then
111          do j=1,4
112            nodes(j)=kontet(j,i)
113          enddo
114!
115!     calculating the volume of the element
116!
117          call calcvol(nodes(1),nodes(2),nodes(3),nodes(4),cotet,
118     &         volume)
119          if(volume.le.0.d0) volume=1.d-30
120!
121!     calculating area of each face in the element
122!
123          call calcsurf(nodes(1),nodes(2),nodes(3),cotet,
124     &         surface(1))
125          call calcsurf(nodes(2),nodes(3),nodes(4),cotet,
126     &         surface(2))
127          call calcsurf(nodes(3),nodes(4),nodes(1),cotet,
128     &         surface(3))
129          call calcsurf(nodes(4),nodes(1),nodes(2),cotet,
130     &         surface(4))
131!
132!     calculating the total surface
133!
134          totsurface=surface(1)+surface(2)+surface(3)+surface(4)
135!
136!     radius of the inscribed sphere
137!
138          radius=3.d0*volume/totsurface
139!
140!     length of each edge
141!
142          do j=1,6
143            n1=nodes(ine(1,j))
144            n2=nodes(ine(2,j))
145            edgelength(j)=dsqrt((cotet(1,n1)-cotet(1,n2))**2+
146     &           (cotet(2,n1)-cotet(2,n2))**2+
147     &           (cotet(3,n1)-cotet(3,n2))**2)
148          enddo
149!
150!     maximum edge length
151!
152          hmax=maxval(edgelength)
153!
154!     quality
155!
156          quality(i)=alpha*hmax/radius
157!
158!     are the next lines really needed?
159!
160c          if(quality(i).lt.1.d0) then
161c            quality(i)=1.d0/quality(i)
162c          endif
163        endif
164      endif
165!
166      return
167      end
168