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 meshqualitycavity(no1,no2,no3,no4,cotet,quality,volume)
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!     version in which a) the element is described by nodes no1..no4
29!     b) the volume is calculated too
30!     used in cavity_refine
31!
32      implicit none
33!
34      integer j,nodes(4),n1,n2,ine(2,6),no1,no2,no3,no4
35!
36      real*8 cotet(3,*),quality,volume,surface(4),totsurface,alpha,
37     &     edgelength(6),hmax,radius
38!
39      data ine /1,2,2,3,1,3,1,4,2,4,3,4/
40!
41!
42!
43!     alpha is the proporionality factor
44!
45      alpha=dsqrt(6.d0)/12.d0
46!
47!     calculating the quality for just one element
48!
49      nodes(1)=no1
50      nodes(2)=no2
51      nodes(3)=no3
52      nodes(4)=no4
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=alpha*hmax/radius
96!
97!     are the next lines really needed?
98!
99c      if(quality.lt.1.d0) then
100c        quality=1.d0/quality
101c      endif
102!
103      return
104      end
105