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