24 const wedgePolyPatch& wpp
27 const polyBoundaryMesh& patches = mesh.boundaryMesh();
29 scalar wppCosAngle = wpp.cosAngle();
37 && isA<wedgePolyPatch>(patches[
patchi])
40 const wedgePolyPatch& pp =
41 refCast<const wedgePolyPatch>(patches[
patchi]);
44 scalar ppCosAngle = wpp.centreNormal() & pp.n();
48 pp.size() == wpp.size()
49 &&
mag(pp.axis() & wpp.axis()) >= (1-1
e-3)
50 &&
mag(ppCosAngle - wppCosAngle) >= 1
e-3
65 const Vector<label>& directions,
70 EdgeMap<label> edgesInError;
76 const polyBoundaryMesh& patches = mesh.boundaryMesh();
79 if (patches[patchi].size() && isA<wedgePolyPatch>(patches[patchi]))
81 const wedgePolyPatch& pp =
82 refCast<const wedgePolyPatch>(patches[
patchi]);
84 scalar wedgeAngle =
acos(pp.cosAngle());
88 Info<<
" Wedge " << pp.name() <<
" with angle " 89 <<
radToDeg(wedgeAngle) <<
" degrees" 96 if (oppositePatchi == -1)
100 Info<<
" ***Cannot find opposite wedge for wedge " 101 << pp.name() <<
endl;
106 const wedgePolyPatch& opp =
107 refCast<const wedgePolyPatch>(patches[oppositePatchi]);
110 if (
mag(opp.axis() & pp.axis()) < (1-1
e-3))
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()
129 const face& f = pp[i];
133 label p1 = f.nextLabel(fp);
134 edgesInError.insert(edge(p0, p1), -1);
140 const point& p0 = p[pp.meshPoints()[0]];
141 forAll(pp.meshPoints(), i)
143 const point& pt = p[pp.meshPoints()[i]];
144 scalar d =
mag((pt - p0) & pp.n());
150 Info<<
" ***Wedge patch " << pp.name() <<
" not planar." 151 <<
" Point " << pt <<
" is not in patch plane by " 164 label nEdgesInError = 0;
168 const face& f = fcs[facei];
173 label p1 = f.nextLabel(fp);
177 scalar magD =
mag(d);
179 if (magD > ROOTVSMALL)
184 label nEmptyDirs = 0;
185 label nNonEmptyDirs = 0;
186 for (
direction cmpt=0; cmpt<vector::nComponents; cmpt++)
188 if (
mag(d[cmpt]) > 1
e-6)
190 if (directions[cmpt] == 0)
205 else if (nEmptyDirs == 1)
208 if (nNonEmptyDirs > 0)
210 if (edgesInError.insert(edge(p0, p1), facei))
216 else if (nEmptyDirs > 1)
219 if (edgesInError.insert(edge(p0, p1), facei))
235 Info<<
" ***Number of edges not aligned with or perpendicular to " 236 <<
"non-empty directions: " << nErrorEdges <<
endl;
241 setPtr->resize(2*nEdgesInError);
246 setPtr->insert(iter.key()[0]);
247 setPtr->insert(iter.key()[1]);
258 Info<<
" All edges aligned with or perpendicular to " 259 <<
"non-empty directions." <<
endl;
269 class transformPositionList
276 const coupledPolyPatch& cpp,
277 List<pointField>& pts
283 List<pointField> newPts(pts.size());
286 newPts[facei].setSize(pts[facei].size());
299 if (facePts.size() > index)
301 ptsAtIndex[facei] = facePts[index];
313 cpp.transformPosition(ptsAtIndex);
319 if (facePts.size() > index)
321 facePts[index] = ptsAtIndex[facei];
336 const polyMesh& mesh,
343 const polyBoundaryMesh& patches = mesh.boundaryMesh();
347 List<pointField> nbrPoints(fcs.size() - mesh.nInternalFaces());
352 if (patches[patchi].coupled())
354 const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
361 label bFacei = cpp.start() + i - mesh.nInternalFaces();
362 const face& f = cpp[i];
363 nbrPoints[bFacei].setSize(f.size());
366 const point& p0 = p[f[fp]];
367 nbrPoints[bFacei][fp] = p0;
377 transformPositionList()
381 label nErrorFaces = 0;
382 scalar avgMismatch = 0;
383 label nCoupledPoints = 0;
387 if (patches[patchi].coupled())
389 const coupledPolyPatch& cpp =
390 refCast<const coupledPolyPatch>(patches[
patchi]);
407 label bFacei = cpp.start() + i - mesh.nInternalFaces();
408 const face& f = cpp[i];
410 if (f.size() != nbrPoints[bFacei].size())
413 <<
"Local face size : " << f.size()
414 <<
" does not equal neighbour face size : " 415 << nbrPoints[bFacei].size()
422 const point& p0 = p[f[fp]];
423 scalar d =
mag(p0 - nbrPoints[bFacei][j]);
425 if (d > smallDist[i])
429 setPtr->insert(cpp.start()+i);
446 reduce(nErrorFaces, sumOp<label>());
447 reduce(avgMismatch, maxOp<scalar>());
448 reduce(nCoupledPoints, sumOp<label>());
450 if (nCoupledPoints > 0)
452 avgMismatch /= nCoupledPoints;
459 Info<<
" **Error in coupled point location: " 461 <<
" faces have their 0th or consecutive vertex not opposite" 462 <<
" their coupled equivalent. Average mismatch " 463 << avgMismatch <<
"." 473 Info<<
" Coupled point location match (average " 474 << avgMismatch <<
") OK." <<
endl;
484 const polyMesh& mesh,
485 const bool allGeometry,
486 const autoPtr<surfaceWriter>& surfWriter,
487 const autoPtr<writer<scalar>>& setWriter
490 label noFailedChecks = 0;
492 Info<<
"\nChecking geometry..." <<
endl;
495 const boundBox& globalBb = mesh.bounds();
497 Info<<
" Overall domain bounding box " 498 << globalBb.min() <<
" " << globalBb.max() <<
endl;
502 scalar minDistSqr =
magSqr(1
e-6 * globalBb.span());
506 Info<<
" Mesh has " << mesh.nGeometricD()
507 <<
" geometric (non-empty/wedge) directions " << validDirs <<
endl;
511 Info<<
" Mesh has " << mesh.nSolutionD()
512 <<
" solution (non-empty) directions " << solDirs <<
endl;
514 if (mesh.nGeometricD() < 3)
516 pointSet nonAlignedPoints(mesh,
"nonAlignedEdges", mesh.nPoints()/100);
522 &&
checkWedges(mesh,
true, validDirs, &nonAlignedPoints)
526 && mesh.checkEdgeAlignment(
true, validDirs, &nonAlignedPoints)
533 nonAlignedPoints.size(),
539 Info<<
" <<Writing " << nNonAligned
540 <<
" points on non-aligned edges to set " 541 << nonAlignedPoints.name() <<
endl;
542 nonAlignedPoints.instance() = mesh.pointsInstance();
543 nonAlignedPoints.
write();
544 if (setWriter.valid())
552 if (mesh.checkClosedBoundary(
true)) noFailedChecks++;
555 cellSet
cells(mesh,
"nonClosedCells", mesh.nCells()/100+1);
556 cellSet aspectCells(mesh,
"highAspectRatioCells", mesh.nCells()/100+1);
559 mesh.checkClosedCells
574 Info<<
" <<Writing " << nNonClosed
575 <<
" non closed cells to set " <<
cells.name() <<
endl;
576 cells.instance() = mesh.pointsInstance();
578 if (surfWriter.valid())
589 Info<<
" <<Writing " << nHighAspect
590 <<
" cells with high aspect ratio to set " 591 << aspectCells.name() <<
endl;
592 aspectCells.instance() = mesh.pointsInstance();
594 if (surfWriter.valid())
602 faceSet faces(mesh,
"zeroAreaFaces", mesh.nFaces()/100+1);
603 if (mesh.checkFaceAreas(
true, &faces))
611 Info<<
" <<Writing " << nFaces
612 <<
" zero area faces to set " << faces.name() <<
endl;
613 faces.instance() = mesh.pointsInstance();
615 if (surfWriter.valid())
624 cellSet
cells(mesh,
"zeroVolumeCells", mesh.nCells()/100+1);
625 if (mesh.checkCellVolumes(
true, &
cells))
633 Info<<
" <<Writing " << nCells
634 <<
" zero volume cells to set " <<
cells.name() <<
endl;
635 cells.instance() = mesh.pointsInstance();
637 if (surfWriter.valid())
646 faceSet faces(mesh,
"nonOrthoFaces", mesh.nFaces()/100+1);
647 if (mesh.checkFaceOrthogonality(
true, &faces))
656 Info<<
" <<Writing " << nFaces
657 <<
" non-orthogonal faces to set " << faces.name() <<
endl;
658 faces.instance() = mesh.pointsInstance();
660 if (surfWriter.valid())
668 faceSet faces(mesh,
"wrongOrientedFaces", mesh.nFaces()/100 + 1);
669 if (mesh.checkFacePyramids(
true, -SMALL, &faces))
677 Info<<
" <<Writing " << nFaces
678 <<
" faces with incorrect orientation to set " 679 << faces.name() <<
endl;
680 faces.instance() = mesh.pointsInstance();
682 if (surfWriter.valid())
691 faceSet faces(mesh,
"skewFaces", mesh.nFaces()/100+1);
692 if (mesh.checkFaceSkewness(
true, &faces))
700 Info<<
" <<Writing " << nFaces
701 <<
" skew faces to set " << faces.name() <<
endl;
702 faces.instance() = mesh.pointsInstance();
704 if (surfWriter.valid())
713 faceSet faces(mesh,
"coupledFaces", mesh.nFaces()/100 + 1);
722 Info<<
" <<Writing " << nFaces
723 <<
" faces with incorrectly matched 0th (or consecutive)" 725 << faces.name() <<
endl;
726 faces.instance() = mesh.pointsInstance();
728 if (surfWriter.valid())
738 faceSet faces(mesh,
"lowQualityTetFaces", mesh.nFaces()/100+1);
756 Info<<
" <<Writing " << nFaces
757 <<
" faces with low quality or negative volume " 758 <<
"decomposition tets to set " << faces.name() <<
endl;
759 faces.instance() = mesh.pointsInstance();
761 if (surfWriter.valid())
772 pointSet
points(mesh,
"shortEdges", mesh.nPoints()/1000 + 1);
773 if (mesh.checkEdgeLength(
true, minDistSqr, &
points))
781 Info<<
" <<Writing " << nPoints
782 <<
" points on short edges to set " <<
points.name()
784 points.instance() = mesh.pointsInstance();
786 if (setWriter.valid())
795 if (mesh.checkPointNearness(
false, minDistSqr, &
points))
801 if (nPoints > nEdgeClose)
803 pointSet nearPoints(mesh,
"nearPoints",
points);
804 Info<<
" <<Writing " << nPoints
805 <<
" near (closer than " <<
Foam::sqrt(minDistSqr)
806 <<
" apart) points to set " << nearPoints.
name() <<
endl;
807 nearPoints.instance() = mesh.pointsInstance();
809 if (setWriter.valid())
819 faceSet faces(mesh,
"concaveFaces", mesh.nFaces()/100 + 1);
820 if (mesh.checkFaceAngles(
true, 10, &faces))
828 Info<<
" <<Writing " << nFaces
829 <<
" faces with concave angles to set " << faces.name()
831 faces.instance() = mesh.pointsInstance();
833 if (surfWriter.valid())
843 faceSet faces(mesh,
"warpedFaces", mesh.nFaces()/100 + 1);
844 if (mesh.checkFaceFlatness(
true, 0.8, &faces))
852 Info<<
" <<Writing " << nFaces
853 <<
" warped faces to set " << faces.name() <<
endl;
854 faces.instance() = mesh.pointsInstance();
856 if (surfWriter.valid())
866 cellSet
cells(mesh,
"underdeterminedCells", mesh.nCells()/100);
867 if (mesh.checkCellDeterminant(
true, &
cells))
873 Info<<
" <<Writing " << nCells
874 <<
" under-determined cells to set " <<
cells.name() <<
endl;
875 cells.instance() = mesh.pointsInstance();
877 if (surfWriter.valid())
886 cellSet
cells(mesh,
"concaveCells", mesh.nCells()/100);
887 if (mesh.checkConcaveCells(
true, &
cells))
893 Info<<
" <<Writing " << nCells
894 <<
" concave cells to set " <<
cells.name() <<
endl;
895 cells.instance() = mesh.pointsInstance();
897 if (surfWriter.valid())
906 faceSet faces(mesh,
"lowWeightFaces", mesh.nFaces()/100);
907 if (mesh.checkFaceWeight(
true, 0.05, &faces))
913 Info<<
" <<Writing " << nFaces
914 <<
" faces with low interpolation weights to set " 915 << faces.name() <<
endl;
916 faces.instance() = mesh.pointsInstance();
918 if (surfWriter.valid())
927 faceSet faces(mesh,
"lowVolRatioFaces", mesh.nFaces()/100);
928 if (mesh.checkVolRatio(
true, 0.01, &faces))
934 Info<<
" <<Writing " << nFaces
935 <<
" faces with low volume ratio cells to set " 936 << faces.name() <<
endl;
937 faces.instance() = mesh.pointsInstance();
939 if (surfWriter.valid())
948 const polyBoundaryMesh& pbm = mesh.boundaryMesh();
950 const word tmName(mesh.time().timeName());
953 autoPtr<surfaceWriter> patchWriter;
954 if (!surfWriter.valid())
956 patchWriter.reset(
new vtkSurfaceWriter());
958 const surfaceWriter& wr =
967 if (isA<cyclicAMIPolyPatch>(pbm[patchi]))
969 const cyclicAMIPolyPatch& cpp =
970 refCast<const cyclicAMIPolyPatch>(pbm[
patchi]);
974 Info<<
"Calculating AMI weights between owner patch: " 975 << cpp.name() <<
" and neighbour patch: " 976 << cpp.neighbPatch().name() <<
endl;
985 autoPtr<globalIndex> globalPoints;
986 autoPtr<globalIndex> globalFaces;
997 uniqueMeshPointLabels,
1006 globalFaces().gather
1010 ami.srcWeightsSum(),
1032 autoPtr<globalIndex> globalPoints;
1033 autoPtr<globalIndex> globalFaces;
1039 cpp.neighbPatch().localFaces(),
1040 cpp.neighbPatch().meshPoints(),
1041 cpp.neighbPatch().meshPointMap(),
1044 uniqueMeshPointLabels,
1053 globalFaces().gather
1057 ami.tgtWeightsSum(),
1080 return noFailedChecks;
dimensionedScalar acos(const dimensionedScalar &ds)
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
const double e
Elementary charge.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Unit conversion functions.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
static bool master(const label communicator=0)
Am I the master process.
Vector< scalar > vector
A scalar version of the templated Vector.
static label worldComm
Default communicator (all processors)
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
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.
const dimensionedScalar e
Elementary charge.
label checkGeometry(const polyMesh &mesh, const bool allGeometry, const autoPtr< surfaceWriter > &, const autoPtr< writer< scalar >> &)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
bool checkWedges(const polyMesh &, const bool report, const Vector< label > &, labelHashSet *)
Check wedge orientation.
List< label > labelList
A List of labels.
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)
errorManip< error > abort(error &err)
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.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
vector point
Point is a vector.
static const Vector< label > one
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.
AMIInterpolation< PrimitivePatch< face, SubList, const pointField & >, PrimitivePatch< face, SubList, const pointField & > > AMIPatchToPatchInterpolation
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
static const scalar minTetQuality
Minimum tetrahedron quality.
static List< int > & procID(label communicator)
Process ID of given process index.