50 const wedgePolyPatch& wpp
53 const polyBoundaryMesh&
patches = mesh.boundaryMesh();
55 scalar wppCosAngle = wpp.cosAngle();
66 const wedgePolyPatch& pp =
70 scalar ppCosAngle = wpp.centreNormal() & pp.n();
74 pp.size() == wpp.size()
75 &&
mag(pp.axis() & wpp.axis()) >= (1-1
e-3)
76 &&
mag(ppCosAngle - wppCosAngle) >= 1
e-3
91 const Vector<label>& directions,
96 EdgeMap<label> edgesInError;
102 const polyBoundaryMesh&
patches = mesh.boundaryMesh();
107 const wedgePolyPatch& pp =
110 scalar wedgeAngle =
acos(pp.cosAngle());
114 Info<<
" Wedge " << pp.name() <<
" with angle "
115 <<
radToDeg(wedgeAngle) <<
" degrees"
122 if (oppositePatchi == -1)
126 Info<<
" ***Cannot find opposite wedge for wedge "
127 << pp.name() <<
endl;
132 const wedgePolyPatch& opp =
133 refCast<const wedgePolyPatch>(
patches[oppositePatchi]);
136 if (
mag(opp.axis() & pp.axis()) < (1-1
e-3))
140 Info<<
" ***Wedges do not have the same axis."
141 <<
" Encountered " << pp.axis()
142 <<
" on patch " << pp.name()
143 <<
" which differs from " << opp.axis()
144 <<
" on opposite wedge patch" << opp.axis()
155 const face&
f = pp[i];
159 label p1 =
f.nextLabel(fp);
160 edgesInError.insert(edge(p0, p1), -1);
166 const point& p0 =
p[pp.meshPoints()[0]];
167 forAll(pp.meshPoints(), i)
169 const point& pt =
p[pp.meshPoints()[i]];
170 scalar d =
mag((pt - p0) & pp.n());
176 Info<<
" ***Wedge patch " << pp.name() <<
" not planar."
177 <<
" Point " << pt <<
" is not in patch plane by "
190 label nEdgesInError = 0;
194 const face&
f = fcs[facei];
199 label p1 =
f.nextLabel(fp);
203 scalar magD =
mag(d);
205 if (magD > rootVSmall)
210 label nEmptyDirs = 0;
211 label nNonEmptyDirs = 0;
212 for (
direction cmpt=0; cmpt<vector::nComponents; cmpt++)
214 if (
mag(d[cmpt]) > 1
e-6)
216 if (directions[cmpt] == 0)
231 else if (nEmptyDirs == 1)
234 if (nNonEmptyDirs > 0)
236 if (edgesInError.insert(edge(p0, p1), facei))
242 else if (nEmptyDirs > 1)
245 if (edgesInError.insert(edge(p0, p1), facei))
261 Info<<
" ***Number of edges not aligned with or perpendicular to "
262 <<
"non-empty directions: " << nErrorEdges <<
endl;
267 setPtr->resize(2*nEdgesInError);
272 setPtr->insert(iter.key()[0]);
273 setPtr->insert(iter.key()[1]);
284 Info<<
" All edges aligned with or perpendicular to "
285 <<
"non-empty directions." <<
endl;
295 class transformPositionList
302 const coupledPolyPatch& cpp,
303 List<pointField>& pts
309 List<pointField> newPts(pts.size());
312 newPts[facei].setSize(pts[facei].size());
325 if (facePts.size() > index)
327 ptsAtIndex[facei] = facePts[index];
339 cpp.transform().transformPosition(ptsAtIndex, ptsAtIndex);
345 if (facePts.size() > index)
347 facePts[index] = ptsAtIndex[facei];
362 const polyMesh& mesh,
369 const polyBoundaryMesh&
patches = mesh.boundaryMesh();
373 List<pointField> nbrPoints(fcs.size() - mesh.nInternalFaces());
380 const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
387 label bFacei = cpp.start() + i - mesh.nInternalFaces();
388 const face&
f = cpp[i];
389 nbrPoints[bFacei].setSize(
f.size());
393 nbrPoints[bFacei][fp] = p0;
398 syncTools::syncBoundaryFaceList
403 transformPositionList()
407 label nErrorFaces = 0;
408 scalar avgMismatch = 0;
409 label nCoupledPoints = 0;
415 const coupledPolyPatch& cpp =
433 label bFacei = cpp.start() + i - mesh.nInternalFaces();
434 const face&
f = cpp[i];
436 if (
f.size() != nbrPoints[bFacei].size())
439 <<
"Local face size : " <<
f.size()
440 <<
" does not equal neighbour face size : "
441 << nbrPoints[bFacei].size()
449 scalar d =
mag(p0 - nbrPoints[bFacei][j]);
451 if (d > smallDist[i])
455 setPtr->insert(cpp.start()+i);
472 reduce(nErrorFaces, sumOp<label>());
473 reduce(avgMismatch, maxOp<scalar>());
474 reduce(nCoupledPoints, sumOp<label>());
476 if (nCoupledPoints > 0)
478 avgMismatch /= nCoupledPoints;
485 Info<<
" **Error in coupled point location: "
487 <<
" faces have their 0th or consecutive vertex not opposite"
488 <<
" their coupled equivalent. Average mismatch "
489 << avgMismatch <<
"."
499 Info<<
" Coupled point location match (average "
500 << avgMismatch <<
") OK." <<
endl;
510 const polyMesh& mesh,
511 const bool allGeometry,
512 const autoPtr<surfaceWriter>& surfWriter,
513 const autoPtr<Foam::setWriter>& setWriter
516 label noFailedChecks = 0;
518 Info<<
"\nChecking geometry..." <<
endl;
521 const boundBox& globalBb = mesh.bounds();
523 Info<<
" Overall domain bounding box "
524 << globalBb.min() <<
" " << globalBb.max() <<
endl;
528 scalar minDistSqr =
magSqr(1
e-6 * globalBb.span());
532 Info<<
" Mesh has " << mesh.nGeometricD()
533 <<
" geometric (non-empty/wedge) directions " << validDirs <<
endl;
537 Info<<
" Mesh has " << mesh.nSolutionD()
538 <<
" solution (non-empty) directions " << solDirs <<
endl;
540 if (mesh.nGeometricD() < 3)
542 pointSet nonAlignedPoints(mesh,
"nonAlignedEdges", mesh.nPoints()/100);
548 &&
checkWedges(mesh,
true, validDirs, &nonAlignedPoints)
552 && mesh.checkEdgeAlignment(
true, validDirs, &nonAlignedPoints)
559 nonAlignedPoints.size(),
565 Info<<
" <<Writing " << nNonAligned
566 <<
" points on non-aligned edges to set "
567 << nonAlignedPoints.name() <<
endl;
568 nonAlignedPoints.instance() = mesh.pointsInstance();
569 nonAlignedPoints.
write();
570 if (setWriter.valid())
578 if (mesh.checkClosedBoundary(
true)) noFailedChecks++;
581 cellSet
cells(mesh,
"nonClosedCells", mesh.nCells()/100+1);
582 cellSet aspectCells(mesh,
"highAspectRatioCells", mesh.nCells()/100+1);
585 mesh.checkClosedCells
600 Info<<
" <<Writing " << nNonClosed
601 <<
" non closed cells to set " <<
cells.name() <<
endl;
602 cells.instance() = mesh.pointsInstance();
604 if (surfWriter.valid())
615 Info<<
" <<Writing " << nHighAspect
616 <<
" cells with high aspect ratio to set "
617 << aspectCells.name() <<
endl;
618 aspectCells.instance() = mesh.pointsInstance();
620 if (surfWriter.valid())
628 faceSet faces(mesh,
"zeroAreaFaces", mesh.nFaces()/100+1);
629 if (mesh.checkFaceAreas(
true, &faces))
637 Info<<
" <<Writing " << nFaces
638 <<
" zero area faces to set " << faces.name() <<
endl;
639 faces.instance() = mesh.pointsInstance();
641 if (surfWriter.valid())
650 cellSet
cells(mesh,
"zeroVolumeCells", mesh.nCells()/100+1);
651 if (mesh.checkCellVolumes(
true, &
cells))
659 Info<<
" <<Writing " << nCells
660 <<
" zero volume cells to set " <<
cells.name() <<
endl;
661 cells.instance() = mesh.pointsInstance();
663 if (surfWriter.valid())
672 faceSet faces(mesh,
"nonOrthoFaces", mesh.nFaces()/100+1);
673 if (mesh.checkFaceOrthogonality(
true, &faces))
682 Info<<
" <<Writing " << nFaces
683 <<
" non-orthogonal faces to set " << faces.name() <<
endl;
684 faces.instance() = mesh.pointsInstance();
686 if (surfWriter.valid())
694 faceSet faces(mesh,
"wrongOrientedFaces", mesh.nFaces()/100 + 1);
695 if (mesh.checkFacePyramids(
true, -small, &faces))
703 Info<<
" <<Writing " << nFaces
704 <<
" faces with incorrect orientation to set "
705 << faces.name() <<
endl;
706 faces.instance() = mesh.pointsInstance();
708 if (surfWriter.valid())
717 faceSet faces(mesh,
"skewFaces", mesh.nFaces()/100+1);
718 if (mesh.checkFaceSkewness(
true, &faces))
726 Info<<
" <<Writing " << nFaces
727 <<
" skew faces to set " << faces.name() <<
endl;
728 faces.instance() = mesh.pointsInstance();
730 if (surfWriter.valid())
739 faceSet faces(mesh,
"coupledFaces", mesh.nFaces()/100 + 1);
748 Info<<
" <<Writing " << nFaces
749 <<
" faces with incorrectly matched 0th (or consecutive)"
751 << faces.name() <<
endl;
752 faces.instance() = mesh.pointsInstance();
754 if (surfWriter.valid())
764 faceSet faces(mesh,
"lowQualityTetFaces", mesh.nFaces()/100+1);
770 polyMeshTetDecomposition::minTetQuality,
782 Info<<
" <<Writing " << nFaces
783 <<
" faces with low quality or negative volume "
784 <<
"decomposition tets to set " << faces.name() <<
endl;
785 faces.instance() = mesh.pointsInstance();
787 if (surfWriter.valid())
798 pointSet
points(mesh,
"shortEdges", mesh.nPoints()/1000 + 1);
799 if (mesh.checkEdgeLength(
true, minDistSqr, &
points))
808 <<
" points on short edges to set " <<
points.name()
810 points.instance() = mesh.pointsInstance();
812 if (setWriter.valid())
821 if (mesh.checkPointNearness(
false, minDistSqr, &
points))
829 pointSet nearPoints(mesh,
"nearPoints",
points);
831 <<
" near (closer than " <<
Foam::sqrt(minDistSqr)
832 <<
" apart) points to set " << nearPoints.
name() <<
endl;
833 nearPoints.instance() = mesh.pointsInstance();
835 if (setWriter.valid())
845 faceSet faces(mesh,
"concaveFaces", mesh.nFaces()/100 + 1);
846 if (mesh.checkFaceAngles(
true, 10, &faces))
854 Info<<
" <<Writing " << nFaces
855 <<
" faces with concave angles to set " << faces.name()
857 faces.instance() = mesh.pointsInstance();
859 if (surfWriter.valid())
869 faceSet faces(mesh,
"warpedFaces", mesh.nFaces()/100 + 1);
870 if (mesh.checkFaceFlatness(
true, 0.8, &faces))
878 Info<<
" <<Writing " << nFaces
879 <<
" warped faces to set " << faces.name() <<
endl;
880 faces.instance() = mesh.pointsInstance();
882 if (surfWriter.valid())
892 cellSet
cells(mesh,
"underdeterminedCells", mesh.nCells()/100);
893 if (mesh.checkCellDeterminant(
true, &
cells))
899 Info<<
" <<Writing " << nCells
900 <<
" under-determined cells to set " <<
cells.name() <<
endl;
901 cells.instance() = mesh.pointsInstance();
903 if (surfWriter.valid())
912 cellSet
cells(mesh,
"concaveCells", mesh.nCells()/100);
913 if (mesh.checkConcaveCells(
true, &
cells))
919 Info<<
" <<Writing " << nCells
920 <<
" concave cells to set " <<
cells.name() <<
endl;
921 cells.instance() = mesh.pointsInstance();
923 if (surfWriter.valid())
932 faceSet faces(mesh,
"lowWeightFaces", mesh.nFaces()/100);
933 if (mesh.checkFaceWeight(
true, 0.05, &faces))
939 Info<<
" <<Writing " << nFaces
940 <<
" faces with low interpolation weights to set "
941 << faces.name() <<
endl;
942 faces.instance() = mesh.pointsInstance();
944 if (surfWriter.valid())
953 faceSet faces(mesh,
"lowVolRatioFaces", mesh.nFaces()/100);
954 if (mesh.checkVolRatio(
true, 0.01, &faces))
960 Info<<
" <<Writing " << nFaces
961 <<
" faces with low volume ratio cells to set "
962 << faces.name() <<
endl;
963 faces.instance() = mesh.pointsInstance();
965 if (surfWriter.valid())
974 const fileName outputPath =
975 mesh.time().globalPath()
976 /functionObjects::writeFile::outputPrefix
977 /(mesh.name() != polyMesh::defaultRegion ? mesh.name() : word())
981 const polyBoundaryMesh&
patches = mesh.boundaryMesh();
984 PtrList<scalarField> patchCoverage(
patches.size());
987 if (isA<nonConformalCyclicPolyPatch>(
patches[nccPatchi]))
989 const nonConformalCyclicPolyPatch& nccPp =
990 refCast<const nonConformalCyclicPolyPatch>
995 const polyPatch& origPp = nccPp.origPatch();
996 const polyPatch& nbrOrigPp = nccPp.nbrPatch().origPatch();
998 const patchToPatches::intersection& intersection =
999 nccPp.intersection();
1001 if (!patchCoverage.set(origPp.index()))
1010 patchCoverage[origPp.index()] +=
1011 intersection.srcCoverage();
1013 if (!patchCoverage.set(nbrOrigPp.index()))
1022 patchCoverage[nbrOrigPp.index()] +=
1023 intersection.tgtCoverage();
1031 if (patchCoverage.set(
patchi))
1038 autoPtr<globalIndex> globalPoints;
1039 autoPtr<globalIndex> globalFaces;
1047 patch.meshPointMap(),
1049 uniqueMeshPointLabels,
1058 globalFaces().gather
1060 UPstream::worldComm,
1061 labelList(UPstream::procID(UPstream::worldComm)),
1067 if (Pstream::master())
1071 mesh.time().writeFormat(),
1072 mesh.time().writeCompression()
1076 patch.name() +
"_coverage",
1088 return noFailedChecks;
#define forAll(list, i)
Loop across all elements in list.
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Routines for checking mesh geometry.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
virtual Ostream & write(const char)=0
Write character.
const word & name() const
Return const reference to name.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
volScalarField scalarField(fieldObject, mesh)
const fvPatchList & patches
const dimensionedScalar e
Elementary charge.
bool checkFaceTets(const bool report, const scalar minPyrVol, const polyMesh &, const vectorField &cellCentres, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *)
Check face tetrahedra volumes.
List< label > labelList
A List of labels.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
label findOppositeWedge(const polyMesh &, const wedgePolyPatch &)
Find wedge with opposite orientation. Note: does not actually check.
Ostream & endl(Ostream &os)
Add newline and flush stream.
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
errorManip< error > abort(error &err)
void mergeAndWrite(const polyMesh &mesh, const surfaceWriter &setWriter, const word &name, const indirectPrimitivePatch setPatch, const fileName &outputDir)
Generate merged surface on master and write. Needs input patch.
vectorField pointField
pointField is a vectorField.
vector point
Point is a vector.
label checkGeometry(const polyMesh &mesh, const bool allGeometry, const autoPtr< surfaceWriter > &, const autoPtr< setWriter > &)
Check the geometry.
Vector< scalar > vector
A scalar version of the templated Vector.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
dimensionedScalar sqrt(const dimensionedScalar &ds)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
dimensioned< scalar > mag(const dimensioned< Type > &)
bool checkWedges(const polyMesh &, const bool report, const Vector< label > &, labelHashSet *)
Check wedge orientation.
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
bool checkCoupledPoints(const polyMesh &, const bool report, labelHashSet *)
Check 0th vertex on coupled faces.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
dimensionedScalar acos(const dimensionedScalar &ds)
Unit conversion functions.