dune-alugrid  2.8-git
meshquality.hh
Go to the documentation of this file.
1 #ifndef DUNE_ALU3DGRID_MESHQUALITY_HH
2 #define DUNE_ALU3DGRID_MESHQUALITY_HH
3 
4 #include <iomanip>
5 #include <string>
6 #include <sstream>
7 #include <fstream>
8 #include <vector>
9 #include <cmath>
10 
11 #include <dune/grid/common/rangegenerators.hh>
12 #include <dune/geometry/referenceelements.hh>
13 
14 namespace Dune {
15 
16  template <class GridView>
17  static void meshQuality( const GridView& gridView, const double time, std::ostream& out )
18  {
19  static const int dim = GridView::dimension;
20 
21  double minVolShortestEdgeRatio = 1e308;
22  double maxVolShortestEdgeRatio = 0;
23  double avgVolShortestEdgeRatio = 0;
24 
25  double minVolLongestEdgeRatio = 1e308;
26  double maxVolLongestEdgeRatio = 0;
27  double avgVolLongestEdgeRatio = 0;
28 
29  double minVolSmallestFaceRatio = 1e308;
30  double maxVolSmallestFaceRatio = 0;
31  double avgVolSmallestFaceRatio = 0;
32 
33  double minVolBiggestFaceRatio = 1e308;
34  double maxVolBiggestFaceRatio = 0;
35  double avgVolBiggestFaceRatio = 0;
36 
37  double minAreaShortestEdgeRatio = 1e308;
38  double maxAreaShortestEdgeRatio = 0;
39  double avgAreaShortestEdgeRatio = 0;
40 
41  double minAreaLongestEdgeRatio = 1e308;
42  double maxAreaLongestEdgeRatio = 0;
43  double avgAreaLongestEdgeRatio = 0;
44 
45  double minSine = 1e308;
46  double maxSine = 0;
47  double avgSine = 0;
48 
49  // TODO: number of elements connected to one edge
50  // TODO: number of elements connected to one vertex
51 
52  //in a regular tetrahedron, we have
53  // volume = (sqrt(2)/12)*edge^3
54  // faceArea = (sqrt(3)/4) *edge^2
55 
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 ) ;
59 
60  std::vector< double > faceVols;
61 
62  const auto& indexSet = gridView.indexSet();
63 
64  std::vector< int > nElementsPerEdge( indexSet.size( dim-1 ), 0 );
65  std::vector< int > nElementsPerVertex( indexSet.size( dim ), 0 );
66 
67  size_t nElements = 0;
68  for( const auto& element : elements( gridView, Dune::Partitions::interiorBorder ) )
69  {
70  const double vol = element.geometry().volume();
71  const Dune::GeometryType geomType = element.type();
72 
73  if( ! geomType.isSimplex() )
74  {
75  if( gridView.comm().rank() == 0 )
76  {
77  std::cout << "MeshQuality check only works for simplex grids, skipping check!" << std::endl;
78  }
79  return ;
80  }
81 
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;
85 
86  double shortestEdge = 1e308 ;
87  double longestEdge = 0;
88 
89  const int edges = element.subEntities( dim-1 );
90  for( int e = 0 ; e < edges ; ++ e )
91  {
92  const auto& edge = element.template subEntity<dim-1>( e );
93  const auto& geo = edge.geometry();
94  assert( geo.corners() == 2 );
95  //( geo.corner( 1 ) - geo.corner( 0 ) ).two_norm();
96  double edgeLength = geo.volume();
97  shortestEdge = std::min( shortestEdge, edgeLength );
98  longestEdge = std::max( longestEdge, edgeLength );
99 
100  const int idx = indexSet.subIndex( element, e, dim-1 );
101  ++(nElementsPerEdge[ idx ]);
102  }
103 
104  double smallestFace = 1e308 ;
105  double biggestFace = 0;
106 
107  const int faces = element.subEntities( dim-2 );
108  faceVols.clear();
109  faceVols.reserve( faces );
110  for( int f = 0 ; f < faces ; ++ f )
111  {
112  const auto& face = element.template subEntity<dim-2>( f );
113  const auto& geo = face.geometry();
114  //assert( geo.corners() == 3 );
115  //( geo.corner( 1 ) - geo.corner( 0 ) ).two_norm();
116  double faceSize = geo.volume();
117  faceVols.push_back( faceSize );
118 
119  smallestFace = std::min( smallestFace, faceSize );
120  biggestFace = std::max( biggestFace, faceSize );
121 
122  double shortestFaceEdge = 1e308 ;
123  double longestFaceEdge = 0;
124 
125  const int faceEdges = face.subEntities( dim-1 );
126  for( int e = 0 ; e < faceEdges ; ++ e )
127  {
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 );
132  //( geo.corner( 1 ) - geo.corner( 0 ) ).two_norm();
133  double edgeLength = geo.volume();
134  shortestFaceEdge = std::min( shortestFaceEdge, edgeLength );
135  longestFaceEdge = std::max( longestFaceEdge, edgeLength );
136  }
137 
138  //in a regular tetrahedron, we have
139  // volume = (sqrt(2)/12)*edge^3
140  shortestFaceEdge = factorFaceEdge * std::pow( shortestFaceEdge, 2 );
141  longestFaceEdge = factorFaceEdge * std::pow( longestFaceEdge, 2 );
142 
143 
144  const double areaShortest = faceSize / shortestFaceEdge ;
145  minAreaShortestEdgeRatio = std::min( minAreaShortestEdgeRatio, areaShortest );
146  maxAreaShortestEdgeRatio = std::max( maxAreaShortestEdgeRatio, areaShortest );
147  avgAreaShortestEdgeRatio += areaShortest;
148 
149  const double areaLongest = faceSize / longestFaceEdge ;
150  minAreaLongestEdgeRatio = std::min( minAreaLongestEdgeRatio, areaLongest );
151  maxAreaLongestEdgeRatio = std::max( maxAreaLongestEdgeRatio, areaLongest );
152  avgAreaLongestEdgeRatio += areaLongest;
153  }
154 
155  const int vertices = element.subEntities( dim );
156 
157  double sumSine = 0 ;
158  for ( int i=0; i<vertices; ++i )
159  {
160  const int idx = indexSet.subIndex( element, i, dim );
161  ++(nElementsPerVertex[ idx ]);
162 
163  double sine = (vol * 6.0);
164  sine *= sine;
165  for( int f=0; f<faces; ++f )
166  {
167  if( f == 3-i ) continue ;
168  sine /= (faceVols[ f ] * 2.0);
169  }
170 
171  minSine = std::min( minSine, sine );
172  maxSine = std::max( maxSine, sine );
173  sumSine += sine ;
174  }
175 
176  sumSine /= double(vertices);
177  avgSine += sumSine ;
178 
179  //in a regular tetrahedron, we have
180  // volume = (sqrt(2)/12)*edge^3
181  shortestEdge = factorEdge * std::pow( shortestEdge, 3 );
182  longestEdge = factorEdge * std::pow( longestEdge, 3 );
183 
184  const double volShortest = vol / shortestEdge;
185  minVolShortestEdgeRatio = std::min( minVolShortestEdgeRatio, volShortest );
186  maxVolShortestEdgeRatio = std::max( maxVolShortestEdgeRatio, volShortest );
187  avgVolShortestEdgeRatio += volShortest;
188 
189  const double volLongest = vol / longestEdge ;
190  minVolLongestEdgeRatio = std::min( minVolLongestEdgeRatio, volLongest );
191  maxVolLongestEdgeRatio = std::max( maxVolLongestEdgeRatio, volLongest );
192  avgVolLongestEdgeRatio += volLongest;
193 
194  //in a regular tetrahedron, we have
195  // volume = (sqrt(2)/12)*edge^3
196  // faceArea = (sqrt(3)/4) *edge^2
197  smallestFace = factorFace * std::pow( smallestFace, 1.5 );
198  biggestFace = factorFace * std::pow( biggestFace, 1.5 );
199 
200  const double volSmallest = vol / smallestFace;
201  minVolSmallestFaceRatio = std::min( minVolSmallestFaceRatio, volSmallest );
202  maxVolSmallestFaceRatio = std::max( maxVolSmallestFaceRatio, volSmallest );
203  avgVolSmallestFaceRatio += volSmallest;
204 
205  const double volBiggest = vol / biggestFace ;
206  minVolBiggestFaceRatio = std::min( minVolBiggestFaceRatio, volBiggest );
207  maxVolBiggestFaceRatio = std::max( maxVolBiggestFaceRatio, volBiggest );
208  avgVolBiggestFaceRatio += volBiggest;
209 
210  ++ nElements;
211  }
212 
213  int maxElementsEdge = 0;
214  int minElementsEdge = 100000000;
215  double avgElementsEdge = 0;
216  double nEdges = nElementsPerEdge.size();
217  for( const auto& nElem : nElementsPerEdge )
218  {
219  maxElementsEdge = std::max( maxElementsEdge, nElem );
220  minElementsEdge = std::min( minElementsEdge, nElem );
221  avgElementsEdge += nElem;
222  }
223 
224  int maxElementsVertex = 0;
225  int minElementsVertex = 100000000;
226  double avgElementsVertex = 0;
227  double nVertices = nElementsPerVertex.size();
228  for( const auto& nElem : nElementsPerVertex )
229  {
230  minElementsVertex = std::min( minElementsVertex, nElem );
231  maxElementsVertex = std::max( maxElementsVertex, nElem );
232  avgElementsVertex += nElem;
233  }
234 
235  nElements = gridView.grid().comm().sum( nElements );
236 
237  minVolShortestEdgeRatio = gridView.grid().comm().min( minVolShortestEdgeRatio );
238  maxVolShortestEdgeRatio = gridView.grid().comm().max( maxVolShortestEdgeRatio );
239  avgVolShortestEdgeRatio = gridView.grid().comm().sum( avgVolShortestEdgeRatio );
240 
241  minVolLongestEdgeRatio = gridView.grid().comm().min( minVolLongestEdgeRatio );
242  maxVolLongestEdgeRatio = gridView.grid().comm().max( maxVolLongestEdgeRatio );
243  avgVolLongestEdgeRatio = gridView.grid().comm().sum( avgVolLongestEdgeRatio );
244 
245  avgVolShortestEdgeRatio /= double(nElements);
246  avgVolLongestEdgeRatio /= double(nElements);
247 
248  minVolSmallestFaceRatio = gridView.grid().comm().min( minVolSmallestFaceRatio );
249  maxVolSmallestFaceRatio = gridView.grid().comm().max( maxVolSmallestFaceRatio );
250  avgVolSmallestFaceRatio = gridView.grid().comm().sum( avgVolSmallestFaceRatio );
251 
252  minVolBiggestFaceRatio = gridView.grid().comm().min( minVolBiggestFaceRatio );
253  maxVolBiggestFaceRatio = gridView.grid().comm().max( maxVolBiggestFaceRatio );
254  avgVolBiggestFaceRatio = gridView.grid().comm().sum( avgVolBiggestFaceRatio );
255 
256  avgVolSmallestFaceRatio /= double(nElements);
257  avgVolBiggestFaceRatio /= double(nElements);
258 
259  minAreaShortestEdgeRatio = gridView.grid().comm().min( minAreaShortestEdgeRatio );
260  maxAreaShortestEdgeRatio = gridView.grid().comm().max( maxAreaShortestEdgeRatio );
261  avgAreaShortestEdgeRatio = gridView.grid().comm().sum( avgAreaShortestEdgeRatio );
262 
263  minAreaLongestEdgeRatio = gridView.grid().comm().min( minAreaLongestEdgeRatio );
264  maxAreaLongestEdgeRatio = gridView.grid().comm().max( maxAreaLongestEdgeRatio );
265  avgAreaLongestEdgeRatio = gridView.grid().comm().sum( avgAreaLongestEdgeRatio );
266 
267  avgSine /= double(nElements);
268 
269  minSine = gridView.grid().comm().min( minSine );
270  maxSine = gridView.grid().comm().max( maxSine );
271  avgSine = gridView.grid().comm().sum( avgSine );
272 
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 );
277 
278  nEdges = gridView.grid().comm().sum( nEdges );
279  nVertices = gridView.grid().comm().sum( nVertices );
280 
281  avgElementsVertex = gridView.grid().comm().sum( avgElementsVertex );
282  avgElementsEdge = gridView.grid().comm().sum( avgElementsEdge );
283 
284  avgElementsEdge /= nEdges;
285  avgElementsVertex /= nVertices;
286 
287  avgAreaShortestEdgeRatio /= double(4 * nElements);
288  avgAreaLongestEdgeRatio /= double(4 * nElements);
289 
290  if( gridView.grid().comm().rank() == 0 )
291  {
292  int space = 18;
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;
294  out << time << " ";
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 << " ";
304  out << std::endl;
305  }
306  }
307 
308 } // end namespace Dune
309 
310 #endif // #ifndef DUNE_ALU3DGRID_MESHQUALITY_HH
Definition: alu3dinclude.hh:63
static void meshQuality(const GridView &gridView, const double time, std::ostream &out)
Definition: meshquality.hh:17
@ simplex
use only simplex elements (i.e., triangles or tetrahedra)
Definition: declaration.hh:18