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,
493 autoPtr<globalIndex> globalPoints;
494 autoPtr<globalIndex> globalFaces;
502 patch.meshPointMap(),
505 uniqueMeshPointLabels,
526 vtkSurfaceWriter().write
542 const polyBoundaryMesh& pbm = mesh.boundaryMesh();
544 const word tmName(mesh.time().timeName());
548 if (isA<cyclicAMIPolyPatch>(pbm[patchi]))
550 const cyclicAMIPolyPatch& cpp =
551 refCast<const cyclicAMIPolyPatch>(pbm[
patchi]);
555 Info<<
"Calculating AMI weights between owner patch: " 556 << cpp.name() <<
" and neighbour patch: " 557 << cpp.neighbPatch().name() <<
endl;
564 fileName(
"postProcessing") /
"src_" + tmName
570 cpp.neighbWeightsSum(),
571 fileName(
"postProcessing") /
"tgt_" + tmName
581 const polyMesh& mesh,
582 const bool allGeometry,
583 const autoPtr<surfaceWriter>& surfWriter,
584 const autoPtr<writer<scalar>>& setWriter
587 label noFailedChecks = 0;
589 Info<<
"\nChecking geometry..." <<
endl;
592 const boundBox& globalBb = mesh.bounds();
594 Info<<
" Overall domain bounding box " 595 << globalBb.min() <<
" " << globalBb.max() <<
endl;
599 scalar minDistSqr =
magSqr(1
e-6 * globalBb.span());
603 Info<<
" Mesh has " << mesh.nGeometricD()
604 <<
" geometric (non-empty/wedge) directions " << validDirs <<
endl;
608 Info<<
" Mesh has " << mesh.nSolutionD()
609 <<
" solution (non-empty) directions " << solDirs <<
endl;
611 if (mesh.nGeometricD() < 3)
613 pointSet nonAlignedPoints(mesh,
"nonAlignedEdges", mesh.nPoints()/100);
619 &&
checkWedges(mesh,
true, validDirs, &nonAlignedPoints)
623 && mesh.checkEdgeAlignment(
true, validDirs, &nonAlignedPoints)
630 nonAlignedPoints.size(),
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())
649 if (mesh.checkClosedBoundary(
true)) noFailedChecks++;
652 cellSet
cells(mesh,
"nonClosedCells", mesh.nCells()/100+1);
653 cellSet aspectCells(mesh,
"highAspectRatioCells", mesh.nCells()/100+1);
656 mesh.checkClosedCells
671 Info<<
" <<Writing " << nNonClosed
672 <<
" non closed cells to set " <<
cells.name() <<
endl;
673 cells.instance() = mesh.pointsInstance();
675 if (surfWriter.valid())
686 Info<<
" <<Writing " << nHighAspect
687 <<
" cells with high aspect ratio to set " 688 << aspectCells.name() <<
endl;
689 aspectCells.instance() = mesh.pointsInstance();
691 if (surfWriter.valid())
699 faceSet faces(mesh,
"zeroAreaFaces", mesh.nFaces()/100+1);
700 if (mesh.checkFaceAreas(
true, &faces))
708 Info<<
" <<Writing " << nFaces
709 <<
" zero area faces to set " << faces.name() <<
endl;
710 faces.instance() = mesh.pointsInstance();
712 if (surfWriter.valid())
721 cellSet
cells(mesh,
"zeroVolumeCells", mesh.nCells()/100+1);
722 if (mesh.checkCellVolumes(
true, &
cells))
730 Info<<
" <<Writing " << nCells
731 <<
" zero volume cells to set " <<
cells.name() <<
endl;
732 cells.instance() = mesh.pointsInstance();
734 if (surfWriter.valid())
743 faceSet faces(mesh,
"nonOrthoFaces", mesh.nFaces()/100+1);
744 if (mesh.checkFaceOrthogonality(
true, &faces))
753 Info<<
" <<Writing " << nFaces
754 <<
" non-orthogonal faces to set " << faces.name() <<
endl;
755 faces.instance() = mesh.pointsInstance();
757 if (surfWriter.valid())
765 faceSet faces(mesh,
"wrongOrientedFaces", mesh.nFaces()/100 + 1);
766 if (mesh.checkFacePyramids(
true, -small, &faces))
774 Info<<
" <<Writing " << nFaces
775 <<
" faces with incorrect orientation to set " 776 << faces.name() <<
endl;
777 faces.instance() = mesh.pointsInstance();
779 if (surfWriter.valid())
788 faceSet faces(mesh,
"skewFaces", mesh.nFaces()/100+1);
789 if (mesh.checkFaceSkewness(
true, &faces))
797 Info<<
" <<Writing " << nFaces
798 <<
" skew faces to set " << faces.name() <<
endl;
799 faces.instance() = mesh.pointsInstance();
801 if (surfWriter.valid())
810 faceSet faces(mesh,
"coupledFaces", mesh.nFaces()/100 + 1);
819 Info<<
" <<Writing " << nFaces
820 <<
" faces with incorrectly matched 0th (or consecutive)" 822 << faces.name() <<
endl;
823 faces.instance() = mesh.pointsInstance();
825 if (surfWriter.valid())
835 faceSet faces(mesh,
"lowQualityTetFaces", mesh.nFaces()/100+1);
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();
858 if (surfWriter.valid())
869 pointSet
points(mesh,
"shortEdges", mesh.nPoints()/1000 + 1);
870 if (mesh.checkEdgeLength(
true, minDistSqr, &
points))
878 Info<<
" <<Writing " << nPoints
879 <<
" points on short edges to set " <<
points.name()
881 points.instance() = mesh.pointsInstance();
883 if (setWriter.valid())
892 if (mesh.checkPointNearness(
false, minDistSqr, &
points))
898 if (nPoints > nEdgeClose)
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();
906 if (setWriter.valid())
916 faceSet faces(mesh,
"concaveFaces", mesh.nFaces()/100 + 1);
917 if (mesh.checkFaceAngles(
true, 10, &faces))
925 Info<<
" <<Writing " << nFaces
926 <<
" faces with concave angles to set " << faces.name()
928 faces.instance() = mesh.pointsInstance();
930 if (surfWriter.valid())
940 faceSet faces(mesh,
"warpedFaces", mesh.nFaces()/100 + 1);
941 if (mesh.checkFaceFlatness(
true, 0.8, &faces))
949 Info<<
" <<Writing " << nFaces
950 <<
" warped faces to set " << faces.name() <<
endl;
951 faces.instance() = mesh.pointsInstance();
953 if (surfWriter.valid())
963 cellSet
cells(mesh,
"underdeterminedCells", mesh.nCells()/100);
964 if (mesh.checkCellDeterminant(
true, &
cells))
970 Info<<
" <<Writing " << nCells
971 <<
" under-determined cells to set " <<
cells.name() <<
endl;
972 cells.instance() = mesh.pointsInstance();
974 if (surfWriter.valid())
983 cellSet
cells(mesh,
"concaveCells", mesh.nCells()/100);
984 if (mesh.checkConcaveCells(
true, &
cells))
990 Info<<
" <<Writing " << nCells
991 <<
" concave cells to set " <<
cells.name() <<
endl;
992 cells.instance() = mesh.pointsInstance();
994 if (surfWriter.valid())
1003 faceSet faces(mesh,
"lowWeightFaces", mesh.nFaces()/100);
1004 if (mesh.checkFaceWeight(
true, 0.05, &faces))
1010 Info<<
" <<Writing " << nFaces
1011 <<
" faces with low interpolation weights to set " 1012 << faces.name() <<
endl;
1013 faces.instance() = mesh.pointsInstance();
1015 if (surfWriter.valid())
1024 faceSet faces(mesh,
"lowVolRatioFaces", mesh.nFaces()/100);
1025 if (mesh.checkVolRatio(
true, 0.01, &faces))
1031 Info<<
" <<Writing " << nFaces
1032 <<
" faces with low volume ratio cells to set " 1033 << faces.name() <<
endl;
1034 faces.instance() = mesh.pointsInstance();
1036 if (surfWriter.valid())
1048 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.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Unit conversion functions.
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.
label checkGeometry(const polyMesh &mesh, const bool allGeometry, const autoPtr< surfaceWriter > &, const autoPtr< writer< scalar >> &)
volScalarField & e
Elementary charge.
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.
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)
vector point
Point is a vector.
void writeAMIWeightsSums(const polyMesh &)
Write out the weights-sums on all the AMI patches.
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.
const doubleScalar e
Elementary charge.
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.
static const scalar minTetQuality
Minimum tetrahedron quality.
static List< int > & procID(label communicator)
Process ID of given process index.