1#ifndef DUNE_ALU3DGRID_MESHQUALITY_HH
2#define DUNE_ALU3DGRID_MESHQUALITY_HH
11#include <dune/grid/common/rangegenerators.hh>
12#include <dune/geometry/referenceelements.hh>
16 template <
class Gr
idView>
17 static void meshQuality(
const GridView& gridView,
const double time, std::ostream& out )
19 static const int dim = GridView::dimension;
21 double minVolShortestEdgeRatio = 1e308;
22 double maxVolShortestEdgeRatio = 0;
23 double avgVolShortestEdgeRatio = 0;
25 double minVolLongestEdgeRatio = 1e308;
26 double maxVolLongestEdgeRatio = 0;
27 double avgVolLongestEdgeRatio = 0;
29 double minVolSmallestFaceRatio = 1e308;
30 double maxVolSmallestFaceRatio = 0;
31 double avgVolSmallestFaceRatio = 0;
33 double minVolBiggestFaceRatio = 1e308;
34 double maxVolBiggestFaceRatio = 0;
35 double avgVolBiggestFaceRatio = 0;
37 double minAreaShortestEdgeRatio = 1e308;
38 double maxAreaShortestEdgeRatio = 0;
39 double avgAreaShortestEdgeRatio = 0;
41 double minAreaLongestEdgeRatio = 1e308;
42 double maxAreaLongestEdgeRatio = 0;
43 double avgAreaLongestEdgeRatio = 0;
45 double minSine = 1e308;
56 const double factorEdgeTet = std::sqrt(2.) / double(12);
57 const double factorFaceEdgeTet = std::sqrt( 3. ) / 4.;
58 const double factorFaceTet = factorEdgeTet * std::pow( ( 1. / factorFaceEdgeTet ), 1.5 ) ;
60 std::vector< double > faceVols;
62 const auto& indexSet = gridView.indexSet();
64 std::vector< int > nElementsPerEdge( indexSet.size( dim-1 ), 0 );
65 std::vector< int > nElementsPerVertex( indexSet.size( dim ), 0 );
68 for(
const auto& element : elements( gridView, Dune::Partitions::interiorBorder ) )
70 const double vol = element.geometry().volume();
71 const Dune::GeometryType geomType = element.type();
73 if( ! geomType.isSimplex() )
75 if( gridView.comm().rank() == 0 )
77 std::cout <<
"MeshQuality check only works for simplex grids, skipping check!" << std::endl;
82 const double factorEdge = geomType.isCube() ? 1.0 : factorEdgeTet;
83 const double factorFaceEdge = geomType.isCube() ? 1.0 : factorFaceEdgeTet;
84 const double factorFace = geomType.isCube() ? 1.0 : factorFaceTet;
86 double shortestEdge = 1e308 ;
87 double longestEdge = 0;
89 const int edges = element.subEntities( dim-1 );
90 for(
int e = 0 ; e < edges ; ++ e )
92 const auto& edge = element.template subEntity<dim-1>( e );
93 const auto& geo = edge.geometry();
94 assert( geo.corners() == 2 );
96 double edgeLength = geo.volume();
97 shortestEdge = std::min( shortestEdge, edgeLength );
98 longestEdge = std::max( longestEdge, edgeLength );
100 const int idx = indexSet.subIndex( element, e, dim-1 );
101 ++(nElementsPerEdge[ idx ]);
104 double smallestFace = 1e308 ;
105 double biggestFace = 0;
107 const int faces = element.subEntities( dim-2 );
109 faceVols.reserve( faces );
110 for(
int f = 0 ; f < faces ; ++ f )
112 const auto& face = element.template subEntity<dim-2>( f );
113 const auto& geo = face.geometry();
116 double faceSize = geo.volume();
117 faceVols.push_back( faceSize );
119 smallestFace = std::min( smallestFace, faceSize );
120 biggestFace = std::max( biggestFace, faceSize );
122 double shortestFaceEdge = 1e308 ;
123 double longestFaceEdge = 0;
125 const int faceEdges = face.subEntities( dim-1 );
126 for(
int e = 0 ; e < faceEdges ; ++ e )
128 const int fe = Dune::ReferenceElements<double, dim>::simplex().subEntity(f,1,e,2);
129 const auto& edge = element.template subEntity<dim-1>( fe );
130 const auto& geo = edge.geometry();
131 assert( geo.corners() == 2 );
133 double edgeLength = geo.volume();
134 shortestFaceEdge = std::min( shortestFaceEdge, edgeLength );
135 longestFaceEdge = std::max( longestFaceEdge, edgeLength );
140 shortestFaceEdge = factorFaceEdge * std::pow( shortestFaceEdge, 2 );
141 longestFaceEdge = factorFaceEdge * std::pow( longestFaceEdge, 2 );
144 const double areaShortest = faceSize / shortestFaceEdge ;
145 minAreaShortestEdgeRatio = std::min( minAreaShortestEdgeRatio, areaShortest );
146 maxAreaShortestEdgeRatio = std::max( maxAreaShortestEdgeRatio, areaShortest );
147 avgAreaShortestEdgeRatio += areaShortest;
149 const double areaLongest = faceSize / longestFaceEdge ;
150 minAreaLongestEdgeRatio = std::min( minAreaLongestEdgeRatio, areaLongest );
151 maxAreaLongestEdgeRatio = std::max( maxAreaLongestEdgeRatio, areaLongest );
152 avgAreaLongestEdgeRatio += areaLongest;
155 const int vertices = element.subEntities( dim );
158 for (
int i=0; i<vertices; ++i )
160 const int idx = indexSet.subIndex( element, i, dim );
161 ++(nElementsPerVertex[ idx ]);
163 double sine = (vol * 6.0);
165 for(
int f=0; f<faces; ++f )
167 if( f == 3-i ) continue ;
168 sine /= (faceVols[ f ] * 2.0);
171 minSine = std::min( minSine, sine );
172 maxSine = std::max( maxSine, sine );
176 sumSine /= double(vertices);
181 shortestEdge = factorEdge * std::pow( shortestEdge, 3 );
182 longestEdge = factorEdge * std::pow( longestEdge, 3 );
184 const double volShortest = vol / shortestEdge;
185 minVolShortestEdgeRatio = std::min( minVolShortestEdgeRatio, volShortest );
186 maxVolShortestEdgeRatio = std::max( maxVolShortestEdgeRatio, volShortest );
187 avgVolShortestEdgeRatio += volShortest;
189 const double volLongest = vol / longestEdge ;
190 minVolLongestEdgeRatio = std::min( minVolLongestEdgeRatio, volLongest );
191 maxVolLongestEdgeRatio = std::max( maxVolLongestEdgeRatio, volLongest );
192 avgVolLongestEdgeRatio += volLongest;
197 smallestFace = factorFace * std::pow( smallestFace, 1.5 );
198 biggestFace = factorFace * std::pow( biggestFace, 1.5 );
200 const double volSmallest = vol / smallestFace;
201 minVolSmallestFaceRatio = std::min( minVolSmallestFaceRatio, volSmallest );
202 maxVolSmallestFaceRatio = std::max( maxVolSmallestFaceRatio, volSmallest );
203 avgVolSmallestFaceRatio += volSmallest;
205 const double volBiggest = vol / biggestFace ;
206 minVolBiggestFaceRatio = std::min( minVolBiggestFaceRatio, volBiggest );
207 maxVolBiggestFaceRatio = std::max( maxVolBiggestFaceRatio, volBiggest );
208 avgVolBiggestFaceRatio += volBiggest;
213 int maxElementsEdge = 0;
214 int minElementsEdge = 100000000;
215 double avgElementsEdge = 0;
216 double nEdges = nElementsPerEdge.size();
217 for(
const auto& nElem : nElementsPerEdge )
219 maxElementsEdge = std::max( maxElementsEdge, nElem );
220 minElementsEdge = std::min( minElementsEdge, nElem );
221 avgElementsEdge += nElem;
224 int maxElementsVertex = 0;
225 int minElementsVertex = 100000000;
226 double avgElementsVertex = 0;
227 double nVertices = nElementsPerVertex.size();
228 for(
const auto& nElem : nElementsPerVertex )
230 minElementsVertex = std::min( minElementsVertex, nElem );
231 maxElementsVertex = std::max( maxElementsVertex, nElem );
232 avgElementsVertex += nElem;
235 nElements = gridView.grid().comm().sum( nElements );
237 minVolShortestEdgeRatio = gridView.grid().comm().min( minVolShortestEdgeRatio );
238 maxVolShortestEdgeRatio = gridView.grid().comm().max( maxVolShortestEdgeRatio );
239 avgVolShortestEdgeRatio = gridView.grid().comm().sum( avgVolShortestEdgeRatio );
241 minVolLongestEdgeRatio = gridView.grid().comm().min( minVolLongestEdgeRatio );
242 maxVolLongestEdgeRatio = gridView.grid().comm().max( maxVolLongestEdgeRatio );
243 avgVolLongestEdgeRatio = gridView.grid().comm().sum( avgVolLongestEdgeRatio );
245 avgVolShortestEdgeRatio /= double(nElements);
246 avgVolLongestEdgeRatio /= double(nElements);
248 minVolSmallestFaceRatio = gridView.grid().comm().min( minVolSmallestFaceRatio );
249 maxVolSmallestFaceRatio = gridView.grid().comm().max( maxVolSmallestFaceRatio );
250 avgVolSmallestFaceRatio = gridView.grid().comm().sum( avgVolSmallestFaceRatio );
252 minVolBiggestFaceRatio = gridView.grid().comm().min( minVolBiggestFaceRatio );
253 maxVolBiggestFaceRatio = gridView.grid().comm().max( maxVolBiggestFaceRatio );
254 avgVolBiggestFaceRatio = gridView.grid().comm().sum( avgVolBiggestFaceRatio );
256 avgVolSmallestFaceRatio /= double(nElements);
257 avgVolBiggestFaceRatio /= double(nElements);
259 minAreaShortestEdgeRatio = gridView.grid().comm().min( minAreaShortestEdgeRatio );
260 maxAreaShortestEdgeRatio = gridView.grid().comm().max( maxAreaShortestEdgeRatio );
261 avgAreaShortestEdgeRatio = gridView.grid().comm().sum( avgAreaShortestEdgeRatio );
263 minAreaLongestEdgeRatio = gridView.grid().comm().min( minAreaLongestEdgeRatio );
264 maxAreaLongestEdgeRatio = gridView.grid().comm().max( maxAreaLongestEdgeRatio );
265 avgAreaLongestEdgeRatio = gridView.grid().comm().sum( avgAreaLongestEdgeRatio );
267 avgSine /= double(nElements);
269 minSine = gridView.grid().comm().min( minSine );
270 maxSine = gridView.grid().comm().max( maxSine );
271 avgSine = gridView.grid().comm().sum( avgSine );
273 minElementsEdge = gridView.grid().comm().min( minElementsEdge );
274 minElementsVertex = gridView.grid().comm().min( minElementsVertex );
275 maxElementsEdge = gridView.grid().comm().max( maxElementsEdge );
276 maxElementsVertex = gridView.grid().comm().max( maxElementsVertex );
278 nEdges = gridView.grid().comm().sum( nEdges );
279 nVertices = gridView.grid().comm().sum( nVertices );
281 avgElementsVertex = gridView.grid().comm().sum( avgElementsVertex );
282 avgElementsEdge = gridView.grid().comm().sum( avgElementsEdge );
284 avgElementsEdge /= nEdges;
285 avgElementsVertex /= nVertices;
287 avgAreaShortestEdgeRatio /= double(4 * nElements);
288 avgAreaLongestEdgeRatio /= double(4 * nElements);
290 if( gridView.grid().comm().rank() == 0 )
293 out <<
"# MeshQuality: time=1 V/LE(min=2,max=3,avg=4) V/SE(min=5,max=6,avg=7) V/LF(min=8,max=9,avg=10) V/SF(min=11,max=12,avg=13) A/LE(min=14,max=15,avg=16) A/SE(min=17,max=18,avg=19) " << std::endl;
295 out << std::setw(space) << minVolLongestEdgeRatio <<
" " << maxVolLongestEdgeRatio <<
" " << avgVolLongestEdgeRatio <<
" ";
296 out << std::setw(space) << minVolShortestEdgeRatio <<
" " << maxVolShortestEdgeRatio <<
" " << avgVolShortestEdgeRatio <<
" " ;
297 out << std::setw(space) << minVolBiggestFaceRatio <<
" " << maxVolBiggestFaceRatio <<
" " << avgVolBiggestFaceRatio <<
" ";
298 out << std::setw(space) << minVolSmallestFaceRatio <<
" " << maxVolSmallestFaceRatio <<
" " << avgVolSmallestFaceRatio <<
" ";
299 out << std::setw(space) << minAreaLongestEdgeRatio <<
" " << maxAreaLongestEdgeRatio <<
" " << avgAreaLongestEdgeRatio <<
" ";
300 out << std::setw(space) << minAreaShortestEdgeRatio <<
" " << maxAreaShortestEdgeRatio <<
" " << avgAreaShortestEdgeRatio <<
" ";
301 out << std::setw(space) << minSine <<
" " << maxSine <<
" " << avgSine <<
" ";
302 out << std::setw(space) << minElementsEdge <<
" " << maxElementsEdge <<
" " << avgElementsEdge <<
" ";
303 out << std::setw(space) << minElementsVertex <<
" " << maxElementsVertex <<
" " << avgElementsVertex <<
" ";
Definition: alu3dinclude.hh:63
static void meshQuality(const GridView &gridView, const double time, std::ostream &out)
Definition: meshquality.hh:17