FreeFOAM The Cross-Platform CFD Toolkit
checkGeometry.C
Go to the documentation of this file.
1 #include "checkGeometry.H"
2 #include <OpenFOAM/polyMesh.H>
3 #include <meshTools/cellSet.H>
4 #include <meshTools/faceSet.H>
5 #include <meshTools/pointSet.H>
6 #include <OpenFOAM/EdgeMap.H>
9 
10 
11 // Find wedge with opposite orientation. Note: does not actually check that
12 // it is opposite, only that it has opposite normal and same axis
13 Foam::label Foam::findOppositeWedge
14 (
15  const polyMesh& mesh,
16  const wedgePolyPatch& wpp
17 )
18 {
19  const polyBoundaryMesh& patches = mesh.boundaryMesh();
20 
21  scalar wppCosAngle = wpp.centreNormal()&wpp.patchNormal();
22 
23  forAll(patches, patchI)
24  {
25  if
26  (
27  patchI != wpp.index()
28  && patches[patchI].size()
29  && isA<wedgePolyPatch>(patches[patchI])
30  )
31  {
32  const wedgePolyPatch& pp = refCast<const wedgePolyPatch>
33  (
34  patches[patchI]
35  );
36 
37  // Calculate (cos of) angle to wpp (not pp!) centre normal
38  scalar ppCosAngle = wpp.centreNormal()&pp.patchNormal();
39 
40  if
41  (
42  pp.size() == wpp.size()
43  && mag(pp.axis() & wpp.axis()) >= (1-1E-3)
44  && mag(ppCosAngle - wppCosAngle) >= 1E-3
45  )
46  {
47  return patchI;
48  }
49  }
50  }
51  return -1;
52 }
53 
54 
56 (
57  const polyMesh& mesh,
58  const bool report,
59  const Vector<label>& directions,
60  labelHashSet* setPtr
61 )
62 {
63  // To mark edges without calculating edge addressing
64  EdgeMap<label> edgesInError;
65 
66  const pointField& p = mesh.points();
67  const faceList& fcs = mesh.faces();
68 
69 
70  const polyBoundaryMesh& patches = mesh.boundaryMesh();
71  forAll(patches, patchI)
72  {
73  if (patches[patchI].size() && isA<wedgePolyPatch>(patches[patchI]))
74  {
75  const wedgePolyPatch& pp = refCast<const wedgePolyPatch>
76  (
77  patches[patchI]
78  );
79 
80  scalar wedgeAngle = acos(pp.centreNormal()&pp.patchNormal());
81 
82  if (report)
83  {
84  Info<< " Wedge " << pp.name() << " with angle "
85  << 180/mathematicalConstant::pi*wedgeAngle << " degrees"
86  << endl;
87  }
88 
89  // Find opposite
90  label oppositePatchI = findOppositeWedge(mesh, pp);
91 
92  if (oppositePatchI == -1)
93  {
94  if (report)
95  {
96  Info<< " ***Cannot find opposite wedge for wedge "
97  << pp.name() << endl;
98  }
99  return true;
100  }
101 
102  const wedgePolyPatch& opp = refCast<const wedgePolyPatch>
103  (
104  patches[oppositePatchI]
105  );
106 
107 
108  if (mag(opp.axis() & pp.axis()) < (1-1E-3))
109  {
110  if (report)
111  {
112  Info<< " ***Wedges do not have the same axis."
113  << " Encountered " << pp.axis()
114  << " on patch " << pp.name()
115  << " which differs from " << opp.axis()
116  << " on opposite wedge patch" << opp.axis()
117  << endl;
118  }
119  return true;
120  }
121 
122 
123 
124  // Mark edges on wedgePatches
125  forAll(pp, i)
126  {
127  const face& f = pp[i];
128  forAll(f, fp)
129  {
130  label p0 = f[fp];
131  label p1 = f.nextLabel(fp);
132  edgesInError.insert(edge(p0, p1), -1); // non-error value
133  }
134  }
135 
136 
137  // Check that wedge patch is flat
138  const point& p0 = p[pp.meshPoints()[0]];
139  forAll(pp.meshPoints(), i)
140  {
141  const point& pt = p[pp.meshPoints()[i]];
142  scalar d = mag((pt-p0) & pp.patchNormal());
143 
144  if (d > sqrt(SMALL))
145  {
146  if (report)
147  {
148  Info<< " ***Wedge patch " << pp.name() << " not planar."
149  << " Point " << pt << " is not in patch plane by "
150  << d << " meter."
151  << endl;
152  }
153  return true;
154  }
155  }
156  }
157  }
158 
159 
160 
161  // Check all non-wedge faces
162  label nEdgesInError = 0;
163 
164  forAll(fcs, faceI)
165  {
166  const face& f = fcs[faceI];
167 
168  forAll(f, fp)
169  {
170  label p0 = f[fp];
171  label p1 = f.nextLabel(fp);
172  if (p0 < p1)
173  {
174  vector d(p[p1]-p[p0]);
175  scalar magD = mag(d);
176 
177  if (magD > ROOTVSMALL)
178  {
179  d /= magD;
180 
181  // Check how many empty directions are used by the edge.
182  label nEmptyDirs = 0;
183  label nNonEmptyDirs = 0;
184  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
185  {
186  if (mag(d[cmpt]) > 1e-6)
187  {
188  if (directions[cmpt] == 0)
189  {
190  nEmptyDirs++;
191  }
192  else
193  {
194  nNonEmptyDirs++;
195  }
196  }
197  }
198 
199  if (nEmptyDirs == 0)
200  {
201  // Purely in ok directions.
202  }
203  else if (nEmptyDirs == 1)
204  {
205  // Ok if purely in empty directions.
206  if (nNonEmptyDirs > 0)
207  {
208  if (edgesInError.insert(edge(p0, p1), faceI))
209  {
210  nEdgesInError++;
211  }
212  }
213  }
214  else if (nEmptyDirs > 1)
215  {
216  // Always an error
217  if (edgesInError.insert(edge(p0, p1), faceI))
218  {
219  nEdgesInError++;
220  }
221  }
222  }
223  }
224  }
225  }
226 
227  label nErrorEdges = returnReduce(nEdgesInError, sumOp<label>());
228 
229  if (nErrorEdges > 0)
230  {
231  if (report)
232  {
233  Info<< " ***Number of edges not aligned with or perpendicular to "
234  << "non-empty directions: " << nErrorEdges << endl;
235  }
236 
237  if (setPtr)
238  {
239  setPtr->resize(2*nEdgesInError);
240  forAllConstIter(EdgeMap<label>, edgesInError, iter)
241  {
242  if (iter() >= 0)
243  {
244  setPtr->insert(iter.key()[0]);
245  setPtr->insert(iter.key()[1]);
246  }
247  }
248  }
249 
250  return true;
251  }
252  else
253  {
254  if (report)
255  {
256  Info<< " All edges aligned with or perpendicular to "
257  << "non-empty directions." << endl;
258  }
259  return false;
260  }
261 }
262 
263 
264 Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
265 {
266  label noFailedChecks = 0;
267 
268  Info<< "\nChecking geometry..." << endl;
269 
270  // Get a small relative length from the bounding box
271  const boundBox& globalBb = mesh.bounds();
272 
273  Info<< " Overall domain bounding box "
274  << globalBb.min() << " " << globalBb.max() << endl;
275 
276 
277  // Min length
278  scalar minDistSqr = magSqr(1e-6 * globalBb.span());
279 
280  // Non-empty directions
281  const Vector<label> validDirs = (mesh.geometricD() + Vector<label>::one)/2;
282  Info<< " Mesh (non-empty, non-wedge) directions " << validDirs << endl;
283 
284  const Vector<label> solDirs = (mesh.solutionD() + Vector<label>::one)/2;
285  Info<< " Mesh (non-empty) directions " << solDirs << endl;
286 
287  if (mesh.nGeometricD() < 3)
288  {
289  pointSet nonAlignedPoints(mesh, "nonAlignedEdges", mesh.nPoints()/100);
290 
291  if
292  (
293  (
294  validDirs != solDirs
295  && checkWedges(mesh, true, validDirs, &nonAlignedPoints)
296  )
297  || (
298  validDirs == solDirs
299  && mesh.checkEdgeAlignment(true, validDirs, &nonAlignedPoints)
300  )
301  )
302  {
303  noFailedChecks++;
304  label nNonAligned = returnReduce
305  (
306  nonAlignedPoints.size(),
307  sumOp<label>()
308  );
309 
310  if (nNonAligned > 0)
311  {
312  Info<< " <<Writing " << nNonAligned
313  << " points on non-aligned edges to set "
314  << nonAlignedPoints.name() << endl;
315  nonAlignedPoints.write();
316  }
317  }
318  }
319 
320  if (mesh.checkClosedBoundary(true)) noFailedChecks++;
321 
322  {
323  cellSet cells(mesh, "nonClosedCells", mesh.nCells()/100+1);
324  cellSet aspectCells(mesh, "highAspectRatioCells", mesh.nCells()/100+1);
325  if
326  (
327  mesh.checkClosedCells
328  (
329  true,
330  &cells,
331  &aspectCells,
332  mesh.geometricD()
333  )
334  )
335  {
336  noFailedChecks++;
337 
338  label nNonClosed = returnReduce(cells.size(), sumOp<label>());
339 
340  if (nNonClosed > 0)
341  {
342  Info<< " <<Writing " << nNonClosed
343  << " non closed cells to set " << cells.name() << endl;
344  cells.write();
345  }
346  }
347 
348  label nHighAspect = returnReduce(aspectCells.size(), sumOp<label>());
349 
350  if (nHighAspect > 0)
351  {
352  Info<< " <<Writing " << nHighAspect
353  << " cells with high aspect ratio to set "
354  << aspectCells.name() << endl;
355  aspectCells.write();
356  }
357  }
358 
359  {
360  faceSet faces(mesh, "zeroAreaFaces", mesh.nFaces()/100 + 1);
361  if (mesh.checkFaceAreas(true, &faces))
362  {
363  noFailedChecks++;
364 
365  label nFaces = returnReduce(faces.size(), sumOp<label>());
366 
367  if (nFaces > 0)
368  {
369  Info<< " <<Writing " << nFaces
370  << " zero area faces to set " << faces.name() << endl;
371  faces.write();
372  }
373  }
374  }
375 
376  {
377  cellSet cells(mesh, "zeroVolumeCells", mesh.nCells()/100 + 1);
378  if (mesh.checkCellVolumes(true, &cells))
379  {
380  noFailedChecks++;
381 
382  label nCells = returnReduce(cells.size(), sumOp<label>());
383 
384  if (nCells > 0)
385  {
386  Info<< " <<Writing " << nCells
387  << " zero volume cells to set " << cells.name() << endl;
388  cells.write();
389  }
390  }
391  }
392 
393  {
394  faceSet faces(mesh, "nonOrthoFaces", mesh.nFaces()/100 + 1);
395  if (mesh.checkFaceOrthogonality(true, &faces))
396  {
397  noFailedChecks++;
398  }
399 
400  label nFaces = returnReduce(faces.size(), sumOp<label>());
401 
402  if (nFaces > 0)
403  {
404  Info<< " <<Writing " << nFaces
405  << " non-orthogonal faces to set " << faces.name() << endl;
406  faces.write();
407  }
408  }
409 
410 
411  {
412  faceSet faces(mesh, "wrongOrientedFaces", mesh.nFaces()/100 + 1);
413  if (mesh.checkFacePyramids(true, -SMALL, &faces))
414  {
415  noFailedChecks++;
416 
417  label nFaces = returnReduce(faces.size(), sumOp<label>());
418 
419  if (nFaces > 0)
420  {
421  Info<< " <<Writing " << nFaces
422  << " faces with incorrect orientation to set "
423  << faces.name() << endl;
424  faces.write();
425  }
426  }
427  }
428 
429  {
430  faceSet faces(mesh, "skewFaces", mesh.nFaces()/100 + 1);
431  if (mesh.checkFaceSkewness(true, &faces))
432  {
433  noFailedChecks++;
434 
435  label nFaces = returnReduce(faces.size(), sumOp<label>());
436 
437  if (nFaces > 0)
438  {
439  Info<< " <<Writing " << nFaces
440  << " skew faces to set " << faces.name() << endl;
441  faces.write();
442  }
443  }
444  }
445 
446  if (allGeometry)
447  {
448  // Note use of nPoints since don't want edge construction.
449  pointSet points(mesh, "shortEdges", mesh.nPoints()/1000 + 1);
450  if (mesh.checkEdgeLength(true, minDistSqr, &points))
451  {
452  //noFailedChecks++;
453 
454  label nPoints = returnReduce(points.size(), sumOp<label>());
455 
456  if (nPoints > 0)
457  {
458  Info<< " <<Writing " << nPoints
459  << " points on short edges to set " << points.name()
460  << endl;
461  points.write();
462  }
463  }
464 
465  label nEdgeClose = returnReduce(points.size(), sumOp<label>());
466 
467  if (mesh.checkPointNearness(false, minDistSqr, &points))
468  {
469  //noFailedChecks++;
470 
471  label nPoints = returnReduce(points.size(), sumOp<label>());
472 
473  if (nPoints > nEdgeClose)
474  {
475  pointSet nearPoints(mesh, "nearPoints", points);
476  Info<< " <<Writing " << nPoints
477  << " near (closer than " << Foam::sqrt(minDistSqr)
478  << " apart) points to set " << nearPoints.name() << endl;
479  nearPoints.write();
480  }
481  }
482  }
483 
484  if (allGeometry)
485  {
486  faceSet faces(mesh, "concaveFaces", mesh.nFaces()/100 + 1);
487  if (mesh.checkFaceAngles(true, 10, &faces))
488  {
489  //noFailedChecks++;
490 
491  label nFaces = returnReduce(faces.size(), sumOp<label>());
492 
493  if (nFaces > 0)
494  {
495  Info<< " <<Writing " << nFaces
496  << " faces with concave angles to set " << faces.name()
497  << endl;
498  faces.write();
499  }
500  }
501  }
502 
503  if (allGeometry)
504  {
505  faceSet faces(mesh, "warpedFaces", mesh.nFaces()/100 + 1);
506  if (mesh.checkFaceFlatness(true, 0.8, &faces))
507  {
508  //noFailedChecks++;
509 
510  label nFaces = returnReduce(faces.size(), sumOp<label>());
511 
512  if (nFaces > 0)
513  {
514  Info<< " <<Writing " << nFaces
515  << " warped faces to set " << faces.name() << endl;
516  faces.write();
517  }
518  }
519  }
520 
521  if (allGeometry)
522  {
523  cellSet cells(mesh, "underdeterminedCells", mesh.nCells()/100);
524  if (mesh.checkCellDeterminant(true, &cells, mesh.geometricD()))
525  {
526  noFailedChecks++;
527 
528  label nCells = returnReduce(cells.size(), sumOp<label>());
529 
530  Info<< " <<Writing " << nCells
531  << " under-determined cells to set " << cells.name() << endl;
532  cells.write();
533  }
534  }
535 
536 
537  return noFailedChecks;
538 }
539 
540 // ************************ vim: set sw=4 sts=4 et: ************************ //