26 const wedgePolyPatch& wpp
29 const polyBoundaryMesh&
patches =
mesh.boundaryMesh();
31 scalar wppCosAngle = wpp.cosAngle();
39 && isA<wedgePolyPatch>(
patches[patchi])
42 const wedgePolyPatch&
pp =
43 refCast<const wedgePolyPatch>(
patches[patchi]);
46 scalar ppCosAngle = wpp.centreNormal() &
pp.n();
50 pp.size() == wpp.size()
51 &&
mag(
pp.axis() & wpp.axis()) >= (1-1
e-3)
52 &&
mag(ppCosAngle - wppCosAngle) >= 1
e-3
67 const Vector<label>& directions,
72 EdgeMap<label> edgesInError;
78 const polyBoundaryMesh&
patches =
mesh.boundaryMesh();
81 if (
patches[patchi].size() && isA<wedgePolyPatch>(
patches[patchi]))
83 const wedgePolyPatch&
pp =
84 refCast<const wedgePolyPatch>(
patches[patchi]);
86 scalar wedgeAngle =
acos(
pp.cosAngle());
90 Info<<
" Wedge " <<
pp.name() <<
" with angle " 91 <<
radToDeg(wedgeAngle) <<
" degrees" 98 if (oppositePatchi == -1)
102 Info<<
" ***Cannot find opposite wedge for wedge " 108 const wedgePolyPatch& opp =
109 refCast<const wedgePolyPatch>(
patches[oppositePatchi]);
112 if (
mag(opp.axis() &
pp.axis()) < (1-1
e-3))
116 Info<<
" ***Wedges do not have the same axis." 117 <<
" Encountered " <<
pp.axis()
118 <<
" on patch " <<
pp.name()
119 <<
" which differs from " << opp.axis()
120 <<
" on opposite wedge patch" << opp.axis()
131 const face&
f =
pp[i];
135 label p1 =
f.nextLabel(fp);
136 edgesInError.insert(edge(
p0, p1), -1);
145 const point& pt =
p[
pp.meshPoints()[i]];
146 scalar d =
mag((pt -
p0) &
pp.n());
152 Info<<
" ***Wedge patch " <<
pp.name() <<
" not planar." 153 <<
" Point " << pt <<
" is not in patch plane by " 166 label nEdgesInError = 0;
170 const face&
f = fcs[facei];
175 label p1 =
f.nextLabel(fp);
179 scalar magD =
mag(d);
181 if (magD > ROOTVSMALL)
186 label nEmptyDirs = 0;
187 label nNonEmptyDirs = 0;
190 if (
mag(d[cmpt]) > 1
e-6)
192 if (directions[cmpt] == 0)
207 else if (nEmptyDirs == 1)
210 if (nNonEmptyDirs > 0)
212 if (edgesInError.insert(edge(
p0, p1), facei))
218 else if (nEmptyDirs > 1)
221 if (edgesInError.insert(edge(
p0, p1), facei))
231 label nErrorEdges =
returnReduce(nEdgesInError, sumOp<label>());
237 Info<<
" ***Number of edges not aligned with or perpendicular to " 238 <<
"non-empty directions: " << nErrorEdges <<
endl;
243 setPtr->reserve(2*nEdgesInError);
248 setPtr->insert(iter.key().first());
249 setPtr->insert(iter.key().second());
259 Info<<
" All edges aligned with or perpendicular to " 260 <<
"non-empty directions." <<
endl;
270 class transformPositionList
277 const coupledPolyPatch& cpp,
278 List<pointField>&
pts 284 List<pointField> newPts(
pts.size());
287 newPts[facei].setSize(
pts[facei].size());
300 if (facePts.size() > index)
302 ptsAtIndex[facei] = facePts[index];
314 cpp.transformPosition(ptsAtIndex);
320 if (facePts.size() > index)
322 facePts[index] = ptsAtIndex[facei];
329 pts.transfer(newPts);
337 const polyMesh&
mesh,
344 const polyBoundaryMesh&
patches =
mesh.boundaryMesh();
348 List<pointField> nbrPoints(fcs.size() -
mesh.nInternalFaces());
355 const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
362 label bFacei = cpp.start() + i -
mesh.nInternalFaces();
363 const face&
f = cpp[i];
364 nbrPoints[bFacei].setSize(
f.size());
368 nbrPoints[bFacei][fp] =
p0;
373 syncTools::syncBoundaryFaceList
378 transformPositionList()
382 label nErrorFaces = 0;
383 scalar avgMismatch = 0;
384 label nCoupledPoints = 0;
390 const coupledPolyPatch& cpp =
391 refCast<const coupledPolyPatch>(
patches[patchi]);
408 label bFacei = cpp.start() + i -
mesh.nInternalFaces();
409 const face&
f = cpp[i];
411 if (
f.size() != nbrPoints[bFacei].size())
414 <<
"Local face size : " <<
f.size()
415 <<
" does not equal neighbour face size : " 416 << nbrPoints[bFacei].size()
424 scalar d =
mag(
p0 - nbrPoints[bFacei][j]);
426 if (d > smallDist[i])
430 setPtr->insert(cpp.start()+i);
447 reduce(nErrorFaces, sumOp<label>());
448 reduce(avgMismatch, maxOp<scalar>());
449 reduce(nCoupledPoints, sumOp<label>());
451 if (nCoupledPoints > 0)
453 avgMismatch /= nCoupledPoints;
460 Info<<
" **Error in coupled point location: " 462 <<
" faces have their 0th or consecutive vertex not opposite" 463 <<
" their coupled equivalent. Average mismatch " 464 << avgMismatch <<
"." 473 Info<<
" Coupled point location match (average " 474 << avgMismatch <<
") OK." <<
endl;
483 const polyMesh&
mesh,
485 const fileName& fName,
489 const Map<label>& meshPointMap,
494 autoPtr<globalIndex>& globalFaces,
495 autoPtr<globalIndex>& globalPoints
508 uniqueMeshPointLabels,
517 globalFaces().gather(weights, mergedWeights);
519 if (Pstream::master())
529 wr.write(
"weightsSum", mergedWeights);
537 const polyMesh&
mesh,
538 const bool allGeometry,
539 autoPtr<surfaceWriter>& surfWriter,
540 autoPtr<coordSetWriter>& setWriter
543 label noFailedChecks = 0;
545 Info<<
"\nChecking geometry..." <<
endl;
548 const boundBox& globalBb =
mesh.bounds();
550 Info<<
" Overall domain bounding box " 551 << globalBb.min() <<
" " << globalBb.max() <<
endl;
555 scalar minDistSqr =
magSqr(1
e-6 * globalBb.span());
559 Info<<
" Mesh has " <<
mesh.nGeometricD()
560 <<
" geometric (non-empty/wedge) directions " << validDirs <<
endl;
564 Info<<
" Mesh has " <<
mesh.nSolutionD()
565 <<
" solution (non-empty) directions " << solDirs <<
endl;
567 if (
mesh.nGeometricD() < 3)
569 pointSet nonAlignedPoints(
mesh,
"nonAlignedEdges",
mesh.nPoints()/100);
579 &&
mesh.checkEdgeAlignment(
true, validDirs, &nonAlignedPoints)
586 nonAlignedPoints.size(),
592 Info<<
" <<Writing " << nNonAligned
593 <<
" points on non-aligned edges to set " 594 << nonAlignedPoints.name() <<
endl;
595 nonAlignedPoints.instance() =
mesh.pointsInstance();
596 nonAlignedPoints.write();
597 if (setWriter && setWriter->enabled())
605 if (
mesh.checkClosedBoundary(
true)) noFailedChecks++;
609 cellSet aspectCells(
mesh,
"highAspectRatioCells",
mesh.nCells()/100+1);
612 mesh.checkClosedCells
627 Info<<
" <<Writing " << nNonClosed
628 <<
" non closed cells to set " <<
cells.name() <<
endl;
629 cells.instance() =
mesh.pointsInstance();
631 if (surfWriter && surfWriter->enabled())
638 label nHighAspect =
returnReduce(aspectCells.size(), sumOp<label>());
642 Info<<
" <<Writing " << nHighAspect
643 <<
" cells with high aspect ratio to set " 644 << aspectCells.name() <<
endl;
645 aspectCells.instance() =
mesh.pointsInstance();
647 if (surfWriter && surfWriter->enabled())
655 faceSet faces(
mesh,
"zeroAreaFaces",
mesh.nFaces()/100+1);
656 if (
mesh.checkFaceAreas(
true, &faces))
660 label nFaces =
returnReduce(faces.size(), sumOp<label>());
664 Info<<
" <<Writing " << nFaces
665 <<
" zero area faces to set " << faces.name() <<
endl;
666 faces.instance() =
mesh.pointsInstance();
668 if (surfWriter && surfWriter->enabled())
686 Info<<
" <<Writing " << nCells
687 <<
" zero volume cells to set " <<
cells.name() <<
endl;
688 cells.instance() =
mesh.pointsInstance();
690 if (surfWriter && surfWriter->enabled())
699 faceSet faces(
mesh,
"nonOrthoFaces",
mesh.nFaces()/100+1);
700 if (
mesh.checkFaceOrthogonality(
true, &faces))
705 label nFaces =
returnReduce(faces.size(), sumOp<label>());
709 Info<<
" <<Writing " << nFaces
710 <<
" non-orthogonal faces to set " << faces.name() <<
endl;
711 faces.instance() =
mesh.pointsInstance();
713 if (surfWriter && surfWriter->enabled())
721 faceSet faces(
mesh,
"wrongOrientedFaces",
mesh.nFaces()/100 + 1);
722 if (
mesh.checkFacePyramids(
true, -SMALL, &faces))
726 label nFaces =
returnReduce(faces.size(), sumOp<label>());
730 Info<<
" <<Writing " << nFaces
731 <<
" faces with incorrect orientation to set " 732 << faces.name() <<
endl;
733 faces.instance() =
mesh.pointsInstance();
735 if (surfWriter && surfWriter->enabled())
744 faceSet faces(
mesh,
"skewFaces",
mesh.nFaces()/100+1);
745 if (
mesh.checkFaceSkewness(
true, &faces))
749 label nFaces =
returnReduce(faces.size(), sumOp<label>());
753 Info<<
" <<Writing " << nFaces
754 <<
" skew faces to set " << faces.name() <<
endl;
755 faces.instance() =
mesh.pointsInstance();
757 if (surfWriter && surfWriter->enabled())
766 faceSet faces(
mesh,
"coupledFaces",
mesh.nFaces()/100 + 1);
771 label nFaces =
returnReduce(faces.size(), sumOp<label>());
775 Info<<
" <<Writing " << nFaces
776 <<
" faces with incorrectly matched 0th (or consecutive)" 778 << faces.name() <<
endl;
779 faces.instance() =
mesh.pointsInstance();
781 if (surfWriter && surfWriter->enabled())
791 faceSet faces(
mesh,
"lowQualityTetFaces",
mesh.nFaces()/100+1);
794 polyMeshTetDecomposition::checkFaceTets
797 polyMeshTetDecomposition::minTetQuality,
805 label nFaces =
returnReduce(faces.size(), sumOp<label>());
809 Info<<
" <<Writing " << nFaces
810 <<
" faces with low quality or negative volume " 811 <<
"decomposition tets to set " << faces.name() <<
endl;
812 faces.instance() =
mesh.pointsInstance();
814 if (surfWriter && surfWriter->enabled())
826 if (
mesh.checkEdgeLength(
true, minDistSqr, &
points))
835 <<
" points on short edges to set " <<
points.name()
839 if (setWriter && setWriter->enabled())
848 if (
mesh.checkPointNearness(
false, minDistSqr, &
points))
856 pointSet nearPoints(
mesh,
"nearPoints",
points);
858 <<
" near (closer than " <<
Foam::sqrt(minDistSqr)
859 <<
" apart) points to set " << nearPoints.
name() <<
endl;
860 nearPoints.instance() =
mesh.pointsInstance();
862 if (setWriter && setWriter->enabled())
872 faceSet faces(
mesh,
"concaveFaces",
mesh.nFaces()/100 + 1);
873 if (
mesh.checkFaceAngles(
true, 10, &faces))
877 label nFaces =
returnReduce(faces.size(), sumOp<label>());
881 Info<<
" <<Writing " << nFaces
882 <<
" faces with concave angles to set " << faces.name()
884 faces.instance() =
mesh.pointsInstance();
886 if (surfWriter && surfWriter->enabled())
896 faceSet faces(
mesh,
"warpedFaces",
mesh.nFaces()/100 + 1);
897 if (
mesh.checkFaceFlatness(
true, 0.8, &faces))
901 label nFaces =
returnReduce(faces.size(), sumOp<label>());
905 Info<<
" <<Writing " << nFaces
906 <<
" warped faces to set " << faces.name() <<
endl;
907 faces.instance() =
mesh.pointsInstance();
909 if (surfWriter && surfWriter->enabled())
919 cellSet
cells(
mesh,
"underdeterminedCells",
mesh.nCells()/100);
920 if (
mesh.checkCellDeterminant(
true, &
cells))
926 Info<<
" <<Writing " << nCells
927 <<
" under-determined cells to set " <<
cells.name() <<
endl;
928 cells.instance() =
mesh.pointsInstance();
930 if (surfWriter && surfWriter->enabled())
940 if (
mesh.checkConcaveCells(
true, &
cells))
946 Info<<
" <<Writing " << nCells
947 <<
" concave cells to set " <<
cells.name() <<
endl;
948 cells.instance() =
mesh.pointsInstance();
950 if (surfWriter && surfWriter->enabled())
959 faceSet faces(
mesh,
"lowWeightFaces",
mesh.nFaces()/100);
960 if (
mesh.checkFaceWeight(
true, 0.05, &faces))
964 label nFaces =
returnReduce(faces.size(), sumOp<label>());
966 Info<<
" <<Writing " << nFaces
967 <<
" faces with low interpolation weights to set " 968 << faces.name() <<
endl;
969 faces.instance() =
mesh.pointsInstance();
971 if (surfWriter && surfWriter->enabled())
980 faceSet faces(
mesh,
"lowVolRatioFaces",
mesh.nFaces()/100);
981 if (
mesh.checkVolRatio(
true, 0.01, &faces))
985 label nFaces =
returnReduce(faces.size(), sumOp<label>());
987 Info<<
" <<Writing " << nFaces
988 <<
" faces with low volume ratio cells to set " 989 << faces.name() <<
endl;
990 faces.instance() =
mesh.pointsInstance();
992 if (surfWriter && surfWriter->enabled())
1001 const polyBoundaryMesh&
pbm =
mesh.boundaryMesh();
1003 const word tmName(
mesh.time().timeName());
1004 const word procAndTime(
Foam::name(Pstream::myProcNo()) +
"_" + tmName);
1006 autoPtr<surfaceWriter> patchWriter;
1009 patchWriter.reset(
new surfaceWriters::vtkWriter());
1012 surfaceWriter& wr = (surfWriter ? *surfWriter : *patchWriter);
1016 const fileName outputDir
1018 mesh.time().globalPath()/functionObject::outputPrefix
1019 /
mesh.regionName()/
"checkMesh" 1024 if (isA<cyclicAMIPolyPatch>(
pbm[patchi]))
1026 const cyclicAMIPolyPatch& cpp =
1027 refCast<const cyclicAMIPolyPatch>(
pbm[patchi]);
1031 Info<<
"Calculating AMI weights between owner patch: " 1032 << cpp.name() <<
" and neighbour patch: " 1033 << cpp.neighbPatch().name() <<
endl;
1041 autoPtr<globalIndex> globalFaces;
1042 autoPtr<globalIndex> globalPoints;
1047 outputDir / cpp.name() +
"-src_" + tmName,
1048 ami.srcWeightsSum(),
1059 if (isA<cyclicACMIPolyPatch>(
pbm[patchi]))
1061 const cyclicACMIPolyPatch&
pp =
1062 refCast<const cyclicACMIPolyPatch>(
pbm[patchi]);
1065 globalFaces().gather(
pp.mask(), mergedMask);
1067 if (Pstream::master())
1073 (outputDir / cpp.name() +
"-src_" + tmName),
1077 wr.write(
"mask", mergedMask);
1083 const auto& nbrPp = cpp.neighbPatch();
1088 autoPtr<globalIndex> globalFaces;
1089 autoPtr<globalIndex> globalPoints;
1096 / nbrPp.name() +
"-tgt_" + tmName
1098 ami.tgtWeightsSum(),
1101 nbrPp.meshPointMap(),
1109 if (isA<cyclicACMIPolyPatch>(
pbm[patchi]))
1111 const cyclicACMIPolyPatch&
pp =
1112 refCast<const cyclicACMIPolyPatch>(
pbm[patchi]);
1114 globalFaces().gather
1116 pp.neighbPatch().mask(),
1120 if (Pstream::master())
1128 / nbrPp.name() +
"-tgt_" + tmName
1133 wr.write(
"mask", mergedMask);
1140 else if (isA<mappedPatchBase>(
pbm[patchi]))
1142 const auto&
pp =
pbm[patchi];
1143 const auto& cpp = refCast<const mappedPatchBase>(
pp);
1149 autoPtr<globalIndex> globalFaces;
1150 autoPtr<globalIndex> globalPoints;
1155 outputDir /
pp.name() +
"-src_" + tmName,
1156 ami.srcWeightsSum(),
1167 if (cpp.sameWorld())
1170 const polyPatch& nbrPp = cpp.samplePolyPatch();
1175 autoPtr<globalIndex> globalFaces;
1176 autoPtr<globalIndex> globalPoints;
1182 outputDir / nbrPp.name() +
"-tgt_" + tmName,
1183 ami.tgtWeightsSum(),
1186 nbrPp.meshPointMap(),
1198 return noFailedChecks;
const polyBoundaryMesh & pbm
dimensionedScalar acos(const dimensionedScalar &ds)
void collectAndWriteAMIWeights(const polyMesh &mesh, surfaceWriter &wr, const fileName &fName, const scalarField &weights, const faceList &localFaces, const labelList &meshPoints, const Map< label > &meshPointMap, faceList &mergedFaces, pointField &mergedPoints, autoPtr< globalIndex > &globalFaces, autoPtr< globalIndex > &globalPoints)
Collect AMI weights to master and write.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
AMIInterpolation AMIPatchToPatchInterpolation
Patch-to-patch interpolation == Foam::AMIInterpolation.
Unit conversion functions.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
::Foam::direction nComponents(const expressions::valueTypeCode) noexcept
The number of components associated with given valueTypeCode.
void mergeAndWrite(const polyMesh &mesh, surfaceWriter &writer, const word &name, const indirectPrimitivePatch &setPatch, const fileName &outputDir)
Generate merged surface on master and write. Needs input patch.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
#define forAll(list, i)
Loop across all elements in list.
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
List< face > faceList
List of faces.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
void clear()
Clear the list, i.e. set size to zero.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
bool checkWedges(const polyMesh &, const bool report, const Vector< label > &, labelHashSet *)
Check wedge orientation.
errorManip< error > abort(error &err)
label findOppositeWedge(const polyMesh &, const wedgePolyPatch &)
bool checkCoupledPoints(const polyMesh &, const bool report, labelHashSet *)
Check 0th vertex on coupled faces.
const word & name() const noexcept
Return const reference to name.
vector point
Point is a vector.
constexpr scalar radToDeg(const scalar rad) noexcept
Conversion from radians to degrees.
const polyBoundaryMesh & patches
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
static const Vector< label > one
messageStream Info
Information stream (stdout output on master, null elsewhere)
List< label > labelList
A List of labels.
const volScalarField & p0
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
forAllConstIters(mixture.phases(), phase)
label checkGeometry(const polyMesh &mesh, const bool allGeometry, autoPtr< surfaceWriter > &surfWriter, autoPtr< coordSetWriter > &setWriter)
static constexpr const zero Zero
Global zero (0)