checkGeometry.C
Go to the documentation of this file.
1 #include "PatchTools.H"
2 #include "checkGeometry.H"
3 #include "polyMesh.H"
4 #include "cellSet.H"
5 #include "faceSet.H"
6 #include "pointSet.H"
7 #include "EdgeMap.H"
8 #include "wedgePolyPatch.H"
9 #include "unitConversion.H"
11 
12 #include "vtkSurfaceWriter.H"
13 #include "writer.H"
14 
15 #include "checkTools.H"
16 #include "cyclicAMIPolyPatch.H"
17 #include "Time.H"
18 
19 // Find wedge with opposite orientation. Note: does not actually check that
20 // it is opposite, only that it has opposite normal and same axis
22 (
23  const polyMesh& mesh,
24  const wedgePolyPatch& wpp
25 )
26 {
27  const polyBoundaryMesh& patches = mesh.boundaryMesh();
28 
29  scalar wppCosAngle = wpp.cosAngle();
30 
31  forAll(patches, patchi)
32  {
33  if
34  (
35  patchi != wpp.index()
36  && patches[patchi].size()
37  && isA<wedgePolyPatch>(patches[patchi])
38  )
39  {
40  const wedgePolyPatch& pp =
41  refCast<const wedgePolyPatch>(patches[patchi]);
42 
43  // Calculate (cos of) angle to wpp (not pp!) centre normal
44  scalar ppCosAngle = wpp.centreNormal() & pp.n();
45 
46  if
47  (
48  pp.size() == wpp.size()
49  && mag(pp.axis() & wpp.axis()) >= (1-1e-3)
50  && mag(ppCosAngle - wppCosAngle) >= 1e-3
51  )
52  {
53  return patchi;
54  }
55  }
56  }
57  return -1;
58 }
59 
60 
62 (
63  const polyMesh& mesh,
64  const bool report,
65  const Vector<label>& directions,
66  labelHashSet* setPtr
67 )
68 {
69  // To mark edges without calculating edge addressing
70  EdgeMap<label> edgesInError;
71 
72  const pointField& p = mesh.points();
73  const faceList& fcs = mesh.faces();
74 
75 
76  const polyBoundaryMesh& patches = mesh.boundaryMesh();
77  forAll(patches, patchi)
78  {
79  if (patches[patchi].size() && isA<wedgePolyPatch>(patches[patchi]))
80  {
81  const wedgePolyPatch& pp =
82  refCast<const wedgePolyPatch>(patches[patchi]);
83 
84  scalar wedgeAngle = acos(pp.cosAngle());
85 
86  if (report)
87  {
88  Info<< " Wedge " << pp.name() << " with angle "
89  << radToDeg(wedgeAngle) << " degrees"
90  << endl;
91  }
92 
93  // Find opposite
94  label oppositePatchi = findOppositeWedge(mesh, pp);
95 
96  if (oppositePatchi == -1)
97  {
98  if (report)
99  {
100  Info<< " ***Cannot find opposite wedge for wedge "
101  << pp.name() << endl;
102  }
103  return true;
104  }
105 
106  const wedgePolyPatch& opp =
107  refCast<const wedgePolyPatch>(patches[oppositePatchi]);
108 
109 
110  if (mag(opp.axis() & pp.axis()) < (1-1e-3))
111  {
112  if (report)
113  {
114  Info<< " ***Wedges do not have the same axis."
115  << " Encountered " << pp.axis()
116  << " on patch " << pp.name()
117  << " which differs from " << opp.axis()
118  << " on opposite wedge patch" << opp.axis()
119  << endl;
120  }
121  return true;
122  }
123 
124 
125 
126  // Mark edges on wedgePatches
127  forAll(pp, i)
128  {
129  const face& f = pp[i];
130  forAll(f, fp)
131  {
132  label p0 = f[fp];
133  label p1 = f.nextLabel(fp);
134  edgesInError.insert(edge(p0, p1), -1); // non-error value
135  }
136  }
137 
138 
139  // Check that wedge patch is flat
140  const point& p0 = p[pp.meshPoints()[0]];
141  forAll(pp.meshPoints(), i)
142  {
143  const point& pt = p[pp.meshPoints()[i]];
144  scalar d = mag((pt - p0) & pp.n());
145 
146  if (d > rootSmall)
147  {
148  if (report)
149  {
150  Info<< " ***Wedge patch " << pp.name() << " not planar."
151  << " Point " << pt << " is not in patch plane by "
152  << d << " metre."
153  << endl;
154  }
155  return true;
156  }
157  }
158  }
159  }
160 
161 
162 
163  // Check all non-wedge faces
164  label nEdgesInError = 0;
165 
166  forAll(fcs, facei)
167  {
168  const face& f = fcs[facei];
169 
170  forAll(f, fp)
171  {
172  label p0 = f[fp];
173  label p1 = f.nextLabel(fp);
174  if (p0 < p1)
175  {
176  vector d(p[p1]-p[p0]);
177  scalar magD = mag(d);
178 
179  if (magD > rootVSmall)
180  {
181  d /= magD;
182 
183  // Check how many empty directions are used by the edge.
184  label nEmptyDirs = 0;
185  label nNonEmptyDirs = 0;
186  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
187  {
188  if (mag(d[cmpt]) > 1e-6)
189  {
190  if (directions[cmpt] == 0)
191  {
192  nEmptyDirs++;
193  }
194  else
195  {
196  nNonEmptyDirs++;
197  }
198  }
199  }
200 
201  if (nEmptyDirs == 0)
202  {
203  // Purely in ok directions.
204  }
205  else if (nEmptyDirs == 1)
206  {
207  // Ok if purely in empty directions.
208  if (nNonEmptyDirs > 0)
209  {
210  if (edgesInError.insert(edge(p0, p1), facei))
211  {
212  nEdgesInError++;
213  }
214  }
215  }
216  else if (nEmptyDirs > 1)
217  {
218  // Always an error
219  if (edgesInError.insert(edge(p0, p1), facei))
220  {
221  nEdgesInError++;
222  }
223  }
224  }
225  }
226  }
227  }
228 
229  label nErrorEdges = returnReduce(nEdgesInError, sumOp<label>());
230 
231  if (nErrorEdges > 0)
232  {
233  if (report)
234  {
235  Info<< " ***Number of edges not aligned with or perpendicular to "
236  << "non-empty directions: " << nErrorEdges << endl;
237  }
238 
239  if (setPtr)
240  {
241  setPtr->resize(2*nEdgesInError);
242  forAllConstIter(EdgeMap<label>, edgesInError, iter)
243  {
244  if (iter() >= 0)
245  {
246  setPtr->insert(iter.key()[0]);
247  setPtr->insert(iter.key()[1]);
248  }
249  }
250  }
251 
252  return true;
253  }
254  else
255  {
256  if (report)
257  {
258  Info<< " All edges aligned with or perpendicular to "
259  << "non-empty directions." << endl;
260  }
261  return false;
262  }
263 }
264 
265 
266 namespace Foam
267 {
268  //- Default transformation behaviour for position
269  class transformPositionList
270  {
271  public:
272 
273  //- Transform patch-based field
274  void operator()
275  (
276  const coupledPolyPatch& cpp,
277  List<pointField>& pts
278  ) const
279  {
280  // Each element of pts is all the points in the face. Convert into
281  // lists of size cpp to transform.
282 
283  List<pointField> newPts(pts.size());
284  forAll(pts, facei)
285  {
286  newPts[facei].setSize(pts[facei].size());
287  }
288 
289  label index = 0;
290  while (true)
291  {
292  label n = 0;
293 
294  // Extract for every face the i'th position
295  pointField ptsAtIndex(pts.size(), Zero);
296  forAll(cpp, facei)
297  {
298  const pointField& facePts = pts[facei];
299  if (facePts.size() > index)
300  {
301  ptsAtIndex[facei] = facePts[index];
302  n++;
303  }
304  }
305 
306  if (n == 0)
307  {
308  break;
309  }
310 
311  // Now ptsAtIndex will have for every face either zero or
312  // the position of the i'th vertex. Transform.
313  cpp.transformPosition(ptsAtIndex);
314 
315  // Extract back from ptsAtIndex into newPts
316  forAll(cpp, facei)
317  {
318  pointField& facePts = newPts[facei];
319  if (facePts.size() > index)
320  {
321  facePts[index] = ptsAtIndex[facei];
322  }
323  }
324 
325  index++;
326  }
327 
328  pts.transfer(newPts);
329  }
330  };
331 }
332 
333 
335 (
336  const polyMesh& mesh,
337  const bool report,
338  labelHashSet* setPtr
339 )
340 {
341  const pointField& p = mesh.points();
342  const faceList& fcs = mesh.faces();
343  const polyBoundaryMesh& patches = mesh.boundaryMesh();
344 
345  // Zero'th point on coupled faces
346  // pointField nbrZeroPoint(fcs.size()-mesh.nInternalFaces(), vector::max);
347  List<pointField> nbrPoints(fcs.size() - mesh.nInternalFaces());
348 
349  // Exchange zero point
350  forAll(patches, patchi)
351  {
352  if (patches[patchi].coupled())
353  {
354  const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
355  (
356  patches[patchi]
357  );
358 
359  forAll(cpp, i)
360  {
361  label bFacei = cpp.start() + i - mesh.nInternalFaces();
362  const face& f = cpp[i];
363  nbrPoints[bFacei].setSize(f.size());
364  forAll(f, fp)
365  {
366  const point& p0 = p[f[fp]];
367  nbrPoints[bFacei][fp] = p0;
368  }
369  }
370  }
371  }
373  (
374  mesh,
375  nbrPoints,
376  eqOp<pointField>(),
377  transformPositionList()
378  );
379 
380  // Compare to local ones. Use same tolerance as for matching
381  label nErrorFaces = 0;
382  scalar avgMismatch = 0;
383  label nCoupledPoints = 0;
384 
385  forAll(patches, patchi)
386  {
387  if (patches[patchi].coupled())
388  {
389  const coupledPolyPatch& cpp =
390  refCast<const coupledPolyPatch>(patches[patchi]);
391 
392  if (cpp.owner())
393  {
394  scalarField smallDist
395  (
396  cpp.calcFaceTol
397  (
398  // cpp.matchTolerance(),
399  cpp,
400  cpp.points(),
401  cpp.faceCentres()
402  )
403  );
404 
405  forAll(cpp, i)
406  {
407  label bFacei = cpp.start() + i - mesh.nInternalFaces();
408  const face& f = cpp[i];
409 
410  if (f.size() != nbrPoints[bFacei].size())
411  {
413  << "Local face size : " << f.size()
414  << " does not equal neighbour face size : "
415  << nbrPoints[bFacei].size()
416  << abort(FatalError);
417  }
418 
419  label fp = 0;
420  forAll(f, j)
421  {
422  const point& p0 = p[f[fp]];
423  scalar d = mag(p0 - nbrPoints[bFacei][j]);
424 
425  if (d > smallDist[i])
426  {
427  if (setPtr)
428  {
429  setPtr->insert(cpp.start()+i);
430  }
431  nErrorFaces++;
432 
433  break;
434  }
435 
436  avgMismatch += d;
437  nCoupledPoints++;
438 
439  fp = f.rcIndex(fp);
440  }
441  }
442  }
443  }
444  }
445 
446  reduce(nErrorFaces, sumOp<label>());
447  reduce(avgMismatch, maxOp<scalar>());
448  reduce(nCoupledPoints, sumOp<label>());
449 
450  if (nCoupledPoints > 0)
451  {
452  avgMismatch /= nCoupledPoints;
453  }
454 
455  if (nErrorFaces > 0)
456  {
457  if (report)
458  {
459  Info<< " **Error in coupled point location: "
460  << nErrorFaces
461  << " faces have their 0th or consecutive vertex not opposite"
462  << " their coupled equivalent. Average mismatch "
463  << avgMismatch << "."
464  << endl;
465  }
466 
467  return true;
468  }
469  else
470  {
471  if (report)
472  {
473  Info<< " Coupled point location match (average "
474  << avgMismatch << ") OK." << endl;
475  }
476 
477  return false;
478  }
479 }
480 
481 
483 (
484  const polyMesh& mesh,
485  const primitivePatch& patch,
486  const scalarField& wghtSum,
487  const fileName& file
488 )
489 {
490  // Collect geometry
491  labelList pointToGlobal;
492  labelList uniqueMeshPointLabels;
493  autoPtr<globalIndex> globalPoints;
494  autoPtr<globalIndex> globalFaces;
495  faceList mergedFaces;
496  pointField mergedPoints;
498  (
499  mesh,
500  patch.localFaces(),
501  patch.meshPoints(),
502  patch.meshPointMap(),
503 
504  pointToGlobal,
505  uniqueMeshPointLabels,
506  globalPoints,
507  globalFaces,
508 
509  mergedFaces,
510  mergedPoints
511  );
512 
513  // Collect field
514  scalarField mergedWeights;
515  globalFaces().gather
516  (
519  wghtSum,
520  mergedWeights
521  );
522 
523  // Write the surface
524  if (Pstream::master())
525  {
526  vtkSurfaceWriter().write
527  (
528  file.path(),
529  file.name(),
530  mergedPoints,
531  mergedFaces,
532  "weightsSum",
533  mergedWeights,
534  false
535  );
536  }
537 }
538 
539 
540 void Foam::writeAMIWeightsSums(const polyMesh& mesh)
541 {
542  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
543 
544  const word tmName(mesh.time().timeName());
545 
546  forAll(pbm, patchi)
547  {
548  if (isA<cyclicAMIPolyPatch>(pbm[patchi]))
549  {
550  const cyclicAMIPolyPatch& cpp =
551  refCast<const cyclicAMIPolyPatch>(pbm[patchi]);
552 
553  if (cpp.owner())
554  {
555  Info<< "Calculating AMI weights between owner patch: "
556  << cpp.name() << " and neighbour patch: "
557  << cpp.neighbPatch().name() << endl;
558 
560  (
561  mesh,
562  cpp,
563  cpp.weightsSum(),
564  fileName("postProcessing") / "src_" + tmName
565  );
567  (
568  mesh,
569  cpp.neighbPatch(),
570  cpp.neighbWeightsSum(),
571  fileName("postProcessing") / "tgt_" + tmName
572  );
573  }
574  }
575  }
576 }
577 
578 
580 (
581  const polyMesh& mesh,
582  const bool allGeometry,
583  const autoPtr<surfaceWriter>& surfWriter,
584  const autoPtr<writer<scalar>>& setWriter
585 )
586 {
587  label noFailedChecks = 0;
588 
589  Info<< "\nChecking geometry..." << endl;
590 
591  // Get a small relative length from the bounding box
592  const boundBox& globalBb = mesh.bounds();
593 
594  Info<< " Overall domain bounding box "
595  << globalBb.min() << " " << globalBb.max() << endl;
596 
597 
598  // Min length
599  scalar minDistSqr = magSqr(1e-6 * globalBb.span());
600 
601  // Geometric directions
602  const Vector<label> validDirs = (mesh.geometricD() + Vector<label>::one)/2;
603  Info<< " Mesh has " << mesh.nGeometricD()
604  << " geometric (non-empty/wedge) directions " << validDirs << endl;
605 
606  // Solution directions
607  const Vector<label> solDirs = (mesh.solutionD() + Vector<label>::one)/2;
608  Info<< " Mesh has " << mesh.nSolutionD()
609  << " solution (non-empty) directions " << solDirs << endl;
610 
611  if (mesh.nGeometricD() < 3)
612  {
613  pointSet nonAlignedPoints(mesh, "nonAlignedEdges", mesh.nPoints()/100);
614 
615  if
616  (
617  (
618  validDirs != solDirs
619  && checkWedges(mesh, true, validDirs, &nonAlignedPoints)
620  )
621  || (
622  validDirs == solDirs
623  && mesh.checkEdgeAlignment(true, validDirs, &nonAlignedPoints)
624  )
625  )
626  {
627  noFailedChecks++;
628  label nNonAligned = returnReduce
629  (
630  nonAlignedPoints.size(),
631  sumOp<label>()
632  );
633 
634  if (nNonAligned > 0)
635  {
636  Info<< " <<Writing " << nNonAligned
637  << " points on non-aligned edges to set "
638  << nonAlignedPoints.name() << endl;
639  nonAlignedPoints.instance() = mesh.pointsInstance();
640  nonAlignedPoints.write();
641  if (setWriter.valid())
642  {
643  mergeAndWrite(setWriter, nonAlignedPoints);
644  }
645  }
646  }
647  }
648 
649  if (mesh.checkClosedBoundary(true)) noFailedChecks++;
650 
651  {
652  cellSet cells(mesh, "nonClosedCells", mesh.nCells()/100+1);
653  cellSet aspectCells(mesh, "highAspectRatioCells", mesh.nCells()/100+1);
654  if
655  (
656  mesh.checkClosedCells
657  (
658  true,
659  &cells,
660  &aspectCells,
661  mesh.geometricD()
662  )
663  )
664  {
665  noFailedChecks++;
666 
667  label nNonClosed = returnReduce(cells.size(), sumOp<label>());
668 
669  if (nNonClosed > 0)
670  {
671  Info<< " <<Writing " << nNonClosed
672  << " non closed cells to set " << cells.name() << endl;
673  cells.instance() = mesh.pointsInstance();
674  cells.write();
675  if (surfWriter.valid())
676  {
677  mergeAndWrite(surfWriter(), cells);
678  }
679  }
680  }
681 
682  label nHighAspect = returnReduce(aspectCells.size(), sumOp<label>());
683 
684  if (nHighAspect > 0)
685  {
686  Info<< " <<Writing " << nHighAspect
687  << " cells with high aspect ratio to set "
688  << aspectCells.name() << endl;
689  aspectCells.instance() = mesh.pointsInstance();
690  aspectCells.write();
691  if (surfWriter.valid())
692  {
693  mergeAndWrite(surfWriter(), aspectCells);
694  }
695  }
696  }
697 
698  {
699  faceSet faces(mesh, "zeroAreaFaces", mesh.nFaces()/100+1);
700  if (mesh.checkFaceAreas(true, &faces))
701  {
702  noFailedChecks++;
703 
704  label nFaces = returnReduce(faces.size(), sumOp<label>());
705 
706  if (nFaces > 0)
707  {
708  Info<< " <<Writing " << nFaces
709  << " zero area faces to set " << faces.name() << endl;
710  faces.instance() = mesh.pointsInstance();
711  faces.write();
712  if (surfWriter.valid())
713  {
714  mergeAndWrite(surfWriter(), faces);
715  }
716  }
717  }
718  }
719 
720  {
721  cellSet cells(mesh, "zeroVolumeCells", mesh.nCells()/100+1);
722  if (mesh.checkCellVolumes(true, &cells))
723  {
724  noFailedChecks++;
725 
726  label nCells = returnReduce(cells.size(), sumOp<label>());
727 
728  if (nCells > 0)
729  {
730  Info<< " <<Writing " << nCells
731  << " zero volume cells to set " << cells.name() << endl;
732  cells.instance() = mesh.pointsInstance();
733  cells.write();
734  if (surfWriter.valid())
735  {
736  mergeAndWrite(surfWriter(), cells);
737  }
738  }
739  }
740  }
741 
742  {
743  faceSet faces(mesh, "nonOrthoFaces", mesh.nFaces()/100+1);
744  if (mesh.checkFaceOrthogonality(true, &faces))
745  {
746  noFailedChecks++;
747  }
748 
749  label nFaces = returnReduce(faces.size(), sumOp<label>());
750 
751  if (nFaces > 0)
752  {
753  Info<< " <<Writing " << nFaces
754  << " non-orthogonal faces to set " << faces.name() << endl;
755  faces.instance() = mesh.pointsInstance();
756  faces.write();
757  if (surfWriter.valid())
758  {
759  mergeAndWrite(surfWriter(), faces);
760  }
761  }
762  }
763 
764  {
765  faceSet faces(mesh, "wrongOrientedFaces", mesh.nFaces()/100 + 1);
766  if (mesh.checkFacePyramids(true, -small, &faces))
767  {
768  noFailedChecks++;
769 
770  label nFaces = returnReduce(faces.size(), sumOp<label>());
771 
772  if (nFaces > 0)
773  {
774  Info<< " <<Writing " << nFaces
775  << " faces with incorrect orientation to set "
776  << faces.name() << endl;
777  faces.instance() = mesh.pointsInstance();
778  faces.write();
779  if (surfWriter.valid())
780  {
781  mergeAndWrite(surfWriter(), faces);
782  }
783  }
784  }
785  }
786 
787  {
788  faceSet faces(mesh, "skewFaces", mesh.nFaces()/100+1);
789  if (mesh.checkFaceSkewness(true, &faces))
790  {
791  noFailedChecks++;
792 
793  label nFaces = returnReduce(faces.size(), sumOp<label>());
794 
795  if (nFaces > 0)
796  {
797  Info<< " <<Writing " << nFaces
798  << " skew faces to set " << faces.name() << endl;
799  faces.instance() = mesh.pointsInstance();
800  faces.write();
801  if (surfWriter.valid())
802  {
803  mergeAndWrite(surfWriter(), faces);
804  }
805  }
806  }
807  }
808 
809  {
810  faceSet faces(mesh, "coupledFaces", mesh.nFaces()/100 + 1);
811  if (checkCoupledPoints(mesh, true, &faces))
812  {
813  noFailedChecks++;
814 
815  label nFaces = returnReduce(faces.size(), sumOp<label>());
816 
817  if (nFaces > 0)
818  {
819  Info<< " <<Writing " << nFaces
820  << " faces with incorrectly matched 0th (or consecutive)"
821  << " vertex to set "
822  << faces.name() << endl;
823  faces.instance() = mesh.pointsInstance();
824  faces.write();
825  if (surfWriter.valid())
826  {
827  mergeAndWrite(surfWriter(), faces);
828  }
829  }
830  }
831  }
832 
833  if (allGeometry)
834  {
835  faceSet faces(mesh, "lowQualityTetFaces", mesh.nFaces()/100+1);
836  if
837  (
839  (
840  mesh,
842  true,
843  &faces
844  )
845  )
846  {
847  noFailedChecks++;
848 
849  label nFaces = returnReduce(faces.size(), sumOp<label>());
850 
851  if (nFaces > 0)
852  {
853  Info<< " <<Writing " << nFaces
854  << " faces with low quality or negative volume "
855  << "decomposition tets to set " << faces.name() << endl;
856  faces.instance() = mesh.pointsInstance();
857  faces.write();
858  if (surfWriter.valid())
859  {
860  mergeAndWrite(surfWriter(), faces);
861  }
862  }
863  }
864  }
865 
866  if (allGeometry)
867  {
868  // Note use of nPoints since don't want edge construction.
869  pointSet points(mesh, "shortEdges", mesh.nPoints()/1000 + 1);
870  if (mesh.checkEdgeLength(true, minDistSqr, &points))
871  {
872  // noFailedChecks++;
873 
874  label nPoints = returnReduce(points.size(), sumOp<label>());
875 
876  if (nPoints > 0)
877  {
878  Info<< " <<Writing " << nPoints
879  << " points on short edges to set " << points.name()
880  << endl;
881  points.instance() = mesh.pointsInstance();
882  points.write();
883  if (setWriter.valid())
884  {
885  mergeAndWrite(setWriter, points);
886  }
887  }
888  }
889 
890  label nEdgeClose = returnReduce(points.size(), sumOp<label>());
891 
892  if (mesh.checkPointNearness(false, minDistSqr, &points))
893  {
894  // noFailedChecks++;
895 
896  label nPoints = returnReduce(points.size(), sumOp<label>());
897 
898  if (nPoints > nEdgeClose)
899  {
900  pointSet nearPoints(mesh, "nearPoints", points);
901  Info<< " <<Writing " << nPoints
902  << " near (closer than " << Foam::sqrt(minDistSqr)
903  << " apart) points to set " << nearPoints.name() << endl;
904  nearPoints.instance() = mesh.pointsInstance();
905  nearPoints.write();
906  if (setWriter.valid())
907  {
908  mergeAndWrite(setWriter, nearPoints);
909  }
910  }
911  }
912  }
913 
914  if (allGeometry)
915  {
916  faceSet faces(mesh, "concaveFaces", mesh.nFaces()/100 + 1);
917  if (mesh.checkFaceAngles(true, 10, &faces))
918  {
919  // noFailedChecks++;
920 
921  label nFaces = returnReduce(faces.size(), sumOp<label>());
922 
923  if (nFaces > 0)
924  {
925  Info<< " <<Writing " << nFaces
926  << " faces with concave angles to set " << faces.name()
927  << endl;
928  faces.instance() = mesh.pointsInstance();
929  faces.write();
930  if (surfWriter.valid())
931  {
932  mergeAndWrite(surfWriter(), faces);
933  }
934  }
935  }
936  }
937 
938  if (allGeometry)
939  {
940  faceSet faces(mesh, "warpedFaces", mesh.nFaces()/100 + 1);
941  if (mesh.checkFaceFlatness(true, 0.8, &faces))
942  {
943  // noFailedChecks++;
944 
945  label nFaces = returnReduce(faces.size(), sumOp<label>());
946 
947  if (nFaces > 0)
948  {
949  Info<< " <<Writing " << nFaces
950  << " warped faces to set " << faces.name() << endl;
951  faces.instance() = mesh.pointsInstance();
952  faces.write();
953  if (surfWriter.valid())
954  {
955  mergeAndWrite(surfWriter(), faces);
956  }
957  }
958  }
959  }
960 
961  if (allGeometry)
962  {
963  cellSet cells(mesh, "underdeterminedCells", mesh.nCells()/100);
964  if (mesh.checkCellDeterminant(true, &cells))
965  {
966  noFailedChecks++;
967 
968  label nCells = returnReduce(cells.size(), sumOp<label>());
969 
970  Info<< " <<Writing " << nCells
971  << " under-determined cells to set " << cells.name() << endl;
972  cells.instance() = mesh.pointsInstance();
973  cells.write();
974  if (surfWriter.valid())
975  {
976  mergeAndWrite(surfWriter(), cells);
977  }
978  }
979  }
980 
981  if (allGeometry)
982  {
983  cellSet cells(mesh, "concaveCells", mesh.nCells()/100);
984  if (mesh.checkConcaveCells(true, &cells))
985  {
986  noFailedChecks++;
987 
988  label nCells = returnReduce(cells.size(), sumOp<label>());
989 
990  Info<< " <<Writing " << nCells
991  << " concave cells to set " << cells.name() << endl;
992  cells.instance() = mesh.pointsInstance();
993  cells.write();
994  if (surfWriter.valid())
995  {
996  mergeAndWrite(surfWriter(), cells);
997  }
998  }
999  }
1000 
1001  if (allGeometry)
1002  {
1003  faceSet faces(mesh, "lowWeightFaces", mesh.nFaces()/100);
1004  if (mesh.checkFaceWeight(true, 0.05, &faces))
1005  {
1006  noFailedChecks++;
1007 
1008  label nFaces = returnReduce(faces.size(), sumOp<label>());
1009 
1010  Info<< " <<Writing " << nFaces
1011  << " faces with low interpolation weights to set "
1012  << faces.name() << endl;
1013  faces.instance() = mesh.pointsInstance();
1014  faces.write();
1015  if (surfWriter.valid())
1016  {
1017  mergeAndWrite(surfWriter(), faces);
1018  }
1019  }
1020  }
1021 
1022  if (allGeometry)
1023  {
1024  faceSet faces(mesh, "lowVolRatioFaces", mesh.nFaces()/100);
1025  if (mesh.checkVolRatio(true, 0.01, &faces))
1026  {
1027  noFailedChecks++;
1028 
1029  label nFaces = returnReduce(faces.size(), sumOp<label>());
1030 
1031  Info<< " <<Writing " << nFaces
1032  << " faces with low volume ratio cells to set "
1033  << faces.name() << endl;
1034  faces.instance() = mesh.pointsInstance();
1035  faces.write();
1036  if (surfWriter.valid())
1037  {
1038  mergeAndWrite(surfWriter(), faces);
1039  }
1040  }
1041  }
1042 
1043  if (allGeometry)
1044  {
1045  writeAMIWeightsSums(mesh);
1046  }
1047 
1048  return noFailedChecks;
1049 }
static void syncBoundaryFaceList(const polyMesh &, UList< T > &, const CombineOp &cop, const TransformOp &top, const bool parRun=Pstream::parRun())
Synchronize values on boundary faces only.
dimensionedScalar acos(const dimensionedScalar &ds)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Unit conversion functions.
uint8_t direction
Definition: direction.H:45
List< face > faceList
Definition: faceListFwd.H:43
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:423
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
static label worldComm
Default communicator (all processors)
Definition: UPstream.H:278
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:208
static bool checkFaceTets(const polyMesh &mesh, scalar tol=minTetQuality, const bool report=false, labelHashSet *setPtr=nullptr)
Check face-decomposition tet volume.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
label checkGeometry(const polyMesh &mesh, const bool allGeometry, const autoPtr< surfaceWriter > &, const autoPtr< writer< scalar >> &)
const cellShapeList & cells
const pointField & points
volScalarField & e
Elementary charge.
Definition: createFields.H:11
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static void gatherAndMerge(const scalar mergeDist, const PrimitivePatch< FaceList, PointField > &p, Field< typename PrimitivePatch< FaceList, PointField >::PointType > &mergedPoints, List< typename PrimitivePatch< FaceList, PointField >::FaceType > &mergedFaces, labelList &pointMergeMap)
Gather points and faces onto master and merge into single patch.
bool checkWedges(const polyMesh &, const bool report, const Vector< label > &, labelHashSet *)
Check wedge orientation.
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
void mergeAndWrite(const polyMesh &mesh, const surfaceWriter &writer, const word &name, const indirectPrimitivePatch setPatch, const fileName &outputDir)
Generate merged surface on master and write. Needs input patch.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
label findOppositeWedge(const polyMesh &, const wedgePolyPatch &)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
bool checkCoupledPoints(const polyMesh &, const bool report, labelHashSet *)
Check 0th vertex on coupled faces.
const word & name() const
Return const reference to name.
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
Addressing for a faceList slice.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
label patchi
vector point
Point is a vector.
Definition: point.H:41
void writeAMIWeightsSums(const polyMesh &)
Write out the weights-sums on all the AMI patches.
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
virtual Ostream & write(const token &)=0
Write next token to stream.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:98
void writeAMIWeightsSum(const polyMesh &, const primitivePatch &, const scalarField &, const fileName &)
Write out the weights-sum on the given AMI patch.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
static const scalar minTetQuality
Minimum tetrahedron quality.
static List< int > & procID(label communicator)
Process ID of given process index.
Definition: UPstream.H:440
Namespace for OpenFOAM.