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