48 const List<label>& toProc
59 label proci = toProc[i];
70 sendMap[proci].setSize(nSend[proci]);
78 label proci = toProc[i];
80 sendMap[proci][nSend[proci]++] = i;
101 forAll(constructMap, proci)
105 label nRecv = recvSizes[proci];
107 constructMap[proci].setSize(nRecv);
109 for (
label i = 0; i < nRecv; i++)
111 constructMap[proci][i] = constructSize++;
116 return autoPtr<mapDistribute>
130 void Foam::backgroundMeshDecomposition::initialRefinement()
137 mesh_.time().timeName(),
144 zeroGradientFvPatchScalarField::typeName
147 const conformationSurfaces& geometry = geometryToConformTo_;
149 decompositionMethod& decomposer = decomposerPtr_();
155 List<volumeType> volumeStatus
166 forAll(volumeStatus, celli)
172 mesh_.cells()[celli].points
179 if (geometry.overlaps(cellBb))
183 else if (geometry.inside(cellBb.midpoint()))
195 labelList refCells = selectRefinementCells
204 meshCutter_.consistentRefinement
211 forAll(newCellsToRefine, nCTRI)
213 label celli = newCellsToRefine[nCTRI];
220 icellWeights[celli] =
max 223 icellWeights[celli]/8.0
227 if (
returnReduce(newCellsToRefine.size(), sumOp<label>()) == 0)
233 polyTopoChange meshMod(mesh_);
236 meshCutter_.setRefinement(newCellsToRefine, meshMod);
239 autoPtr<mapPolyMesh> map = meshMod.changeMesh
249 mesh_.updateMesh(map);
252 meshCutter_.updateMesh(map);
257 const labelList& cellMap = map().cellMap();
259 List<volumeType> newVolumeStatus(cellMap.size());
263 label oldCelli = cellMap[newCelli];
271 newVolumeStatus[newCelli] = volumeStatus[oldCelli];
275 volumeStatus.transfer(newVolumeStatus);
278 Info<<
" Background mesh refined from " 280 <<
" to " << mesh_.globalData().nTotalCells()
281 <<
" cells." <<
endl;
285 forAll(volumeStatus, celli)
291 mesh_.cells()[celli].points
298 if (geometry.overlaps(cellBb))
302 else if (geometry.inside(cellBb.midpoint()))
314 bool removeOutsideCells =
false;
316 if (removeOutsideCells)
318 DynamicList<label> cellsToRemove;
320 forAll(volumeStatus, celli)
324 cellsToRemove.append(celli);
328 removeCells cellRemover(mesh_);
331 polyTopoChange meshMod(mesh_);
333 labelList exposedFaces = cellRemover.getExposedFaces
339 cellRemover.setRefinement
348 autoPtr<mapPolyMesh> map = meshMod.changeMesh
358 mesh_.updateMesh(map);
361 meshCutter_.updateMesh(map);
362 cellRemover.updateMesh(map);
367 const labelList& cellMap = map().cellMap();
369 List<volumeType> newVolumeStatus(cellMap.size());
373 label oldCelli = cellMap[newCelli];
381 newVolumeStatus[newCelli] = volumeStatus[oldCelli];
385 volumeStatus.transfer(newVolumeStatus);
390 - mesh_.globalData().nTotalCells()
391 <<
" cells." <<
endl;
403 labelList newDecomp = decomposer.decompose
410 fvMeshDistribute distributor(mesh_, mergeDist_);
412 autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute
417 meshCutter_.distribute(mapDist);
419 mapDist().distributeCellData(volumeStatus);
423 printMeshData(mesh_);
446 void Foam::backgroundMeshDecomposition::printMeshData
453 globalIndex globalCells(mesh.nCells());
478 Info<<
"Processor " << proci <<
" " 479 <<
"Number of cells = " << globalCells.localSize(proci)
503 bool Foam::backgroundMeshDecomposition::refineCell
507 scalar& weightEstimate
517 mesh_.cells()[celli].points
524 weightEstimate = 1.0;
628 Foam::labelList Foam::backgroundMeshDecomposition::selectRefinementCells
630 List<volumeType>& volumeStatus,
639 forAll(volumeStatus, celli)
643 if (meshCutter_.cellLevel()[celli] < minLevels_)
645 cellsToRefine.insert(celli);
661 cellsToRefine.insert(celli);
666 return cellsToRefine.toc();
670 void Foam::backgroundMeshDecomposition::buildPatchAndTree()
677 mesh_.nFaces() - mesh_.nInternalFaces(),
678 mesh_.nInternalFaces()
683 boundaryFacesPtr_.reset
687 tmpBoundaryFaces.localFaces(),
688 tmpBoundaryFaces.localPoints()
693 treeBoundBox overallBb(boundaryFacesPtr_().localPoints());
697 new indexedOctree<treeDataBPatch>
705 overallBb.extend(1e-4),
718 point bbMin(great, great, great);
719 point bbMax(-great, -great, -great);
721 forAll(allBackgroundMeshBounds_, proci)
723 bbMin =
min(bbMin, allBackgroundMeshBounds_[proci].
min());
724 bbMax =
max(bbMax, allBackgroundMeshBounds_[proci].
max());
727 globalBackgroundBounds_ = treeBoundBox(bbMin, bbMax);
734 /
"backgroundMeshDecomposition_proc_" 736 +
"_boundaryFaces.obj" 739 const faceList& faces = boundaryFacesPtr_().localFaces();
740 const List<point>& points = boundaryFacesPtr_().localPoints();
742 Map<label> foamToObj(points.size());
748 const face&
f = faces[i];
752 if (foamToObj.insert(f[fPI], vertI))
763 fStr<<
' ' << foamToObj[f[fPI]] + 1;
774 Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
778 const conformationSurfaces& geometryToConformTo,
779 const dictionary& coeffsDict
783 geometryToConformTo_(geometryToConformTo),
789 "backgroundMeshDecomposition",
793 IOobject::AUTO_WRITE,
805 allBackgroundMeshBounds_(Pstream::nProcs()),
806 globalBackgroundBounds_(),
814 IOobject::MUST_READ_IF_MODIFIED,
818 decomposerPtr_(decompositionMethod::
New(decomposeDict_)),
819 mergeDist_(1e-6*mesh_.bounds().
mag()),
823 coeffsDict.lookupOrDefault<scalar>(
"minCellSizeLimit", 0.0)
832 <<
"This cannot be used when not running in parallel." 836 if (!decomposerPtr_().parallelAware())
839 <<
"You have selected decomposition method " 840 << decomposerPtr_().typeName
841 <<
" which is not parallel aware." << endl
845 Info<<
nl <<
"Building initial background mesh decomposition" <<
endl;
879 label nOccupiedCells = 0;
883 if (icellWeights[cI] > 1 - small)
893 scalar cellWeightLimit =
max 896 *
sum(cellWeights).value()
903 Info<<
" cellWeightLimit " << cellWeightLimit <<
endl;
905 Pout<<
" sum(cellWeights) " <<
sum(cellWeights.primitiveField())
906 <<
" max(cellWeights) " <<
max(cellWeights.primitiveField())
914 if (icellWeights[cWI] > cellWeightLimit)
916 cellsToRefine.insert(cWI);
920 if (
returnReduce(cellsToRefine.size(), sumOp<label>()) == 0)
928 meshCutter_.consistentRefinement
935 if (debug && !cellsToRefine.empty())
937 Pout<<
" cellWeights too large in " << cellsToRefine.size()
941 forAll(newCellsToRefine, nCTRI)
943 label celli = newCellsToRefine[nCTRI];
945 icellWeights[celli] /= 8.0;
949 polyTopoChange meshMod(mesh_);
952 meshCutter_.setRefinement(newCellsToRefine, meshMod);
955 autoPtr<mapPolyMesh> map = meshMod.changeMesh
965 mesh_.updateMesh(map);
968 meshCutter_.updateMesh(map);
970 Info<<
" Background mesh refined from " 972 <<
" to " << mesh_.globalData().nTotalCells()
973 <<
" cells." <<
endl;
986 printMeshData(mesh_);
988 Pout<<
" Pre distribute sum(cellWeights) " 990 <<
" max(cellWeights) " 995 labelList newDecomp = decomposerPtr_().decompose
1002 Info<<
" Redistributing background mesh cells" <<
endl;
1004 fvMeshDistribute distributor(mesh_, mergeDist_);
1006 autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute(newDecomp);
1008 meshCutter_.distribute(mapDist);
1012 printMeshData(mesh_);
1014 Pout<<
" Post distribute sum(cellWeights) " 1015 <<
sum(icellWeights)
1016 <<
" max(cellWeights) " 1017 <<
max(icellWeights)
1023 cellWeights.
write();
1026 buildPatchAndTree();
1045 const List<point>& pts
1048 boolList posProc(pts.size(),
true);
1052 posProc[pI] = positionOnThisProcessor(pts[pI]);
1061 const treeBoundBox& box
1065 return !bFTreePtr_().findBox(box).empty();
1071 const point& centre,
1072 const scalar radiusSqr
1077 return bFTreePtr_().findNearest(centre, radiusSqr).hit();
1087 return bFTreePtr_().findLine(start, end);
1097 return bFTreePtr_().findLineAny(start, end);
1103 const List<point>& pts
1106 DynamicList<label> toCandidateProc;
1107 DynamicList<point> testPoints;
1111 label nTotalCandidates = 0;
1115 const point& pt = pts[pI];
1117 label nCandidates = 0;
1119 forAll(allBackgroundMeshBounds_, proci)
1123 if (allBackgroundMeshBounds_[proci].overlaps(pt,
sqr(small*100)))
1125 toCandidateProc.append(proci);
1126 testPoints.append(pt);
1132 ptBlockStart[pI] = nTotalCandidates;
1133 ptBlockSize[pI] = nCandidates;
1135 nTotalCandidates += nCandidates;
1139 label preDistributionToCandidateProcSize = toCandidateProc.size();
1141 autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1143 map().distribute(testPoints);
1145 List<scalar> distanceSqrToCandidate(testPoints.size(),
sqr(great));
1158 distanceSqrToCandidate[tPI] =
magSqr 1160 testPoints[tPI] - info.hitPoint()
1165 map().reverseDistribute
1167 preDistributionToCandidateProcSize,
1168 distanceSqrToCandidate
1171 labelList ptNearestProc(pts.size(), -1);
1177 SubList<scalar> ptNearestProcResults
1179 distanceSqrToCandidate,
1184 scalar nearestProcDistSqr = great;
1186 forAll(ptNearestProcResults, pPRI)
1188 if (ptNearestProcResults[pPRI] < nearestProcDistSqr)
1190 nearestProcDistSqr = ptNearestProcResults[pPRI];
1192 ptNearestProc[pI] = toCandidateProc[ptBlockStart[pI] + pPRI];
1198 Pout<< pts[pI] <<
" nearestProcDistSqr " << nearestProcDistSqr
1199 <<
" ptNearestProc[pI] " << ptNearestProc[pI] <<
endl;
1202 if (ptNearestProc[pI] < 0)
1205 <<
"The position " << pts[pI]
1206 <<
" did not find a nearest point on the background mesh." 1211 return ptNearestProc;
1219 const List<point>& starts,
1220 const List<point>& ends,
1221 bool includeOwnProcessor
1224 DynamicList<label> toCandidateProc;
1225 DynamicList<point> testStarts;
1226 DynamicList<point> testEnds;
1227 labelList segmentBlockStart(starts.size(), -1);
1228 labelList segmentBlockSize(starts.size(), -1);
1230 label nTotalCandidates = 0;
1234 const point& s = starts[sI];
1235 const point& e = ends[sI];
1240 label nCandidates = 0;
1242 forAll(allBackgroundMeshBounds_, proci)
1249 && allBackgroundMeshBounds_[proci].intersects(s, e,
p)
1252 toCandidateProc.append(proci);
1253 testStarts.append(s);
1260 segmentBlockStart[sI] = nTotalCandidates;
1261 segmentBlockSize[sI] = nCandidates;
1263 nTotalCandidates += nCandidates;
1267 label preDistributionToCandidateProcSize = toCandidateProc.size();
1269 autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1271 map().distribute(testStarts);
1272 map().distribute(testEnds);
1274 List<pointIndexHit> segmentIntersectsCandidate(testStarts.size());
1279 const point& s = testStarts[sI];
1280 const point& e = testEnds[sI];
1283 segmentIntersectsCandidate[sI] = bFTreePtr_().findLine(s, e);
1286 map().reverseDistribute
1288 preDistributionToCandidateProcSize,
1289 segmentIntersectsCandidate
1292 List<List<pointIndexHit>> segmentHitProcs(starts.size());
1295 DynamicList<pointIndexHit> tmpProcHits;
1299 tmpProcHits.clear();
1303 SubList<pointIndexHit> segmentProcResults
1305 segmentIntersectsCandidate,
1306 segmentBlockSize[sI],
1307 segmentBlockStart[sI]
1310 forAll(segmentProcResults, sPRI)
1312 if (segmentProcResults[sPRI].hit())
1314 tmpProcHits.append(segmentProcResults[sPRI]);
1316 tmpProcHits.last().setIndex
1318 toCandidateProc[segmentBlockStart[sI] + sPRI]
1323 segmentHitProcs[sI] = tmpProcHits;
1326 return segmentHitProcs;
1332 const point& centre,
1333 const scalar& radiusSqr
1336 forAll(allBackgroundMeshBounds_, proci)
1338 if (bFTreePtr_().findNearest(centre, radiusSqr).hit())
1350 const point& centre,
1351 const scalar radiusSqr
1356 forAll(allBackgroundMeshBounds_, proci)
1362 && allBackgroundMeshBounds_[proci].overlaps(centre, radiusSqr)
1368 toProc.append(proci);
List< labelList > labelListList
A List of labelList.
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
bool positionOnThisProcessor(const point &pt) const
Is the given position inside the domain of this decomposition.
#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.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Ostream & endl(Ostream &os)
Add newline and flush stream.
static scalar & perturbTol()
Get the perturbation tolerance.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
PointIndexHit< point > pointIndexHit
labelList overlapProcessors(const point ¢re, const scalar radiusSqr) const
DimensionedField< scalar, volMesh > Internal
Type of the internal field from which this GeometricField is derived.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
static autoPtr< mapDistribute > buildMap(const List< label > &toProc)
Build a mapDistribute for the supplied destination processor data.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
List< bool > boolList
Bool container classes.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
treeDataPrimitivePatch< bPatch > treeDataBPatch
static void exchangeSizes(const Container &sendData, labelList &sizes, const label comm=UPstream::worldComm)
Helper: exchange sizes of sendData. sendData is the data per.
labelList processorNearestPosition(const List< point > &pts) const
What is the nearest processor to the given position?
stressControl lookup("compactNormalStress") >> compactNormalStress
bool overlapsOtherProcessors(const point ¢re, const scalar &radiusSqr) const
pointIndexHit findLine(const point &start, const point &end) const
Find nearest intersection of line between start and end, (exposing.
List< List< pointIndexHit > > intersectsProcessors(const List< point > &starts, const List< point > &ends, bool includeOwnProcessor=false) const
Which processors are intersected by the line segment, returns all.
List< label > labelList
A List of labels.
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if successful.
PrimitivePatch< face, List, const pointField > bPatch
label readLabel(Istream &is)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end, (exposing.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
PrimitivePatch< face, SubList, const pointField & > primitivePatch
Addressing for a faceList slice.
word name(const complex &)
Return a string representation of a complex.
static bool & parRun()
Is this a parallel run?
static label nProcs(const label communicator=0)
Number of processes in parallel run.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
vector point
Point is a vector.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
prefixOSstream Pout(cout, "Pout")
dimensioned< scalar > mag(const dimensioned< Type > &)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
autoPtr< mapDistributePolyMesh > distribute(volScalarField &cellWeights)
Redistribute the background mesh based on a supplied weight field,.
virtual Ostream & write(const token &)=0
Write next token to stream.
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
~backgroundMeshDecomposition()
Destructor.
int system(const std::string &command)
Execute the specified command.
bool overlapsThisProcessor(const treeBoundBox &box) const
Does the given box overlap the faces of the boundary of this.