47 const List<label>& toProc
58 label proci = toProc[i];
69 sendMap[proci].setSize(nSend[proci]);
77 label proci = toProc[i];
79 sendMap[proci][nSend[proci]++] = i;
100 forAll(constructMap, proci)
104 label nRecv = recvSizes[proci];
106 constructMap[proci].setSize(nRecv);
108 for (
label i = 0; i < nRecv; i++)
110 constructMap[proci][i] = constructSize++;
115 return autoPtr<mapDistribute>
129 void Foam::backgroundMeshDecomposition::initialRefinement()
136 mesh_.time().timeName(),
143 zeroGradientFvPatchScalarField::typeName
146 const conformationSurfaces& geometry = geometryToConformTo_;
148 decompositionMethod& decomposer = decomposerPtr_();
154 List<volumeType> volumeStatus
165 forAll(volumeStatus, celli)
171 mesh_.cells()[celli].points
178 if (geometry.overlaps(cellBb))
182 else if (geometry.inside(cellBb.midpoint()))
194 labelList refCells = selectRefinementCells
203 meshCutter_.consistentRefinement
210 forAll(newCellsToRefine, nCTRI)
212 label celli = newCellsToRefine[nCTRI];
219 icellWeights[celli] =
max 222 icellWeights[celli]/8.0
226 if (
returnReduce(newCellsToRefine.size(), sumOp<label>()) == 0)
232 polyTopoChange meshMod(mesh_);
235 meshCutter_.setRefinement(newCellsToRefine, meshMod);
238 autoPtr<mapPolyMesh> map = meshMod.changeMesh
248 mesh_.updateMesh(map);
251 meshCutter_.updateMesh(map);
256 const labelList& cellMap = map().cellMap();
258 List<volumeType> newVolumeStatus(cellMap.size());
262 label oldCelli = cellMap[newCelli];
270 newVolumeStatus[newCelli] = volumeStatus[oldCelli];
274 volumeStatus.transfer(newVolumeStatus);
277 Info<<
" Background mesh refined from " 279 <<
" to " << mesh_.globalData().nTotalCells()
280 <<
" cells." <<
endl;
284 forAll(volumeStatus, celli)
290 mesh_.cells()[celli].points
297 if (geometry.overlaps(cellBb))
301 else if (geometry.inside(cellBb.midpoint()))
313 bool removeOutsideCells =
false;
315 if (removeOutsideCells)
317 DynamicList<label> cellsToRemove;
319 forAll(volumeStatus, celli)
323 cellsToRemove.append(celli);
327 removeCells cellRemover(mesh_);
330 polyTopoChange meshMod(mesh_);
332 labelList exposedFaces = cellRemover.getExposedFaces
338 cellRemover.setRefinement
347 autoPtr<mapPolyMesh> map = meshMod.changeMesh
357 mesh_.updateMesh(map);
360 meshCutter_.updateMesh(map);
361 cellRemover.updateMesh(map);
366 const labelList& cellMap = map().cellMap();
368 List<volumeType> newVolumeStatus(cellMap.size());
372 label oldCelli = cellMap[newCelli];
380 newVolumeStatus[newCelli] = volumeStatus[oldCelli];
384 volumeStatus.transfer(newVolumeStatus);
389 - mesh_.globalData().nTotalCells()
390 <<
" cells." <<
endl;
402 labelList newDecomp = decomposer.decompose
409 fvMeshDistribute distributor(mesh_);
411 autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute
416 meshCutter_.distribute(mapDist);
418 mapDist().distributeCellData(volumeStatus);
422 printMeshData(mesh_);
445 void Foam::backgroundMeshDecomposition::printMeshData
452 globalIndex globalCells(mesh.nCells());
477 Info<<
"Processor " << proci <<
" " 478 <<
"Number of cells = " << globalCells.localSize(proci)
502 bool Foam::backgroundMeshDecomposition::refineCell
506 scalar& weightEstimate
516 mesh_.cells()[celli].points
523 weightEstimate = 1.0;
627 Foam::labelList Foam::backgroundMeshDecomposition::selectRefinementCells
629 List<volumeType>& volumeStatus,
638 forAll(volumeStatus, celli)
642 if (meshCutter_.cellLevel()[celli] < minLevels_)
644 cellsToRefine.insert(celli);
660 cellsToRefine.insert(celli);
665 return cellsToRefine.toc();
669 void Foam::backgroundMeshDecomposition::buildPatchAndTree()
676 mesh_.nFaces() - mesh_.nInternalFaces(),
677 mesh_.nInternalFaces()
682 boundaryFacesPtr_.reset
686 tmpBoundaryFaces.localFaces(),
687 tmpBoundaryFaces.localPoints()
692 treeBoundBox overallBb(boundaryFacesPtr_().localPoints());
696 new indexedOctree<treeDataBPatch>
704 overallBb.extend(1e-4),
717 point bbMin(great, great, great);
718 point bbMax(-great, -great, -great);
720 forAll(allBackgroundMeshBounds_, proci)
722 bbMin =
min(bbMin, allBackgroundMeshBounds_[proci].
min());
723 bbMax =
max(bbMax, allBackgroundMeshBounds_[proci].
max());
726 globalBackgroundBounds_ = treeBoundBox(bbMin, bbMax);
733 /
"backgroundMeshDecomposition_proc_" 735 +
"_boundaryFaces.obj" 738 const faceList& faces = boundaryFacesPtr_().localFaces();
739 const List<point>& points = boundaryFacesPtr_().localPoints();
741 Map<label> foamToObj(points.size());
747 const face&
f = faces[i];
751 if (foamToObj.insert(f[fPI], vertI))
762 fStr<<
' ' << foamToObj[f[fPI]] + 1;
777 const conformationSurfaces& geometryToConformTo,
778 const dictionary& coeffsDict
782 geometryToConformTo_(geometryToConformTo),
787 "backgroundMeshDecomposition",
791 IOobject::AUTO_WRITE,
803 allBackgroundMeshBounds_(Pstream::nProcs()),
804 globalBackgroundBounds_(),
812 IOobject::MUST_READ_IF_MODIFIED,
816 decomposerPtr_(decompositionMethod::
New(decomposeDict_)),
817 spanScale_(coeffsDict.
lookup<scalar>(
"spanScale")),
820 coeffsDict.lookupOrDefault<scalar>(
"minCellSizeLimit", 0.0)
823 volRes_(coeffsDict.
lookup<
label>(
"sampleResolution")),
824 maxCellWeightCoeff_(coeffsDict.
lookup<scalar>(
"maxCellWeightCoeff"))
829 <<
"This cannot be used when not running in parallel." 833 if (!decomposerPtr_().parallelAware())
836 <<
"You have selected decomposition method " 837 << decomposerPtr_().typeName
838 <<
" which is not parallel aware." << endl
842 Info<<
nl <<
"Building initial background mesh decomposition" <<
endl;
876 label nOccupiedCells = 0;
880 if (icellWeights[cI] > 1 - small)
890 scalar cellWeightLimit =
max 893 *
sum(cellWeights).value()
900 Info<<
" cellWeightLimit " << cellWeightLimit <<
endl;
902 Pout<<
" sum(cellWeights) " <<
sum(cellWeights.primitiveField())
903 <<
" max(cellWeights) " <<
max(cellWeights.primitiveField())
911 if (icellWeights[cWI] > cellWeightLimit)
913 cellsToRefine.insert(cWI);
917 if (
returnReduce(cellsToRefine.size(), sumOp<label>()) == 0)
925 meshCutter_.consistentRefinement
932 if (debug && !cellsToRefine.empty())
934 Pout<<
" cellWeights too large in " << cellsToRefine.size()
938 forAll(newCellsToRefine, nCTRI)
940 label celli = newCellsToRefine[nCTRI];
942 icellWeights[celli] /= 8.0;
946 polyTopoChange meshMod(mesh_);
949 meshCutter_.setRefinement(newCellsToRefine, meshMod);
952 autoPtr<mapPolyMesh> map = meshMod.changeMesh
962 mesh_.updateMesh(map);
965 meshCutter_.updateMesh(map);
967 Info<<
" Background mesh refined from " 969 <<
" to " << mesh_.globalData().nTotalCells()
970 <<
" cells." <<
endl;
983 printMeshData(mesh_);
985 Pout<<
" Pre distribute sum(cellWeights) " 987 <<
" max(cellWeights) " 992 labelList newDecomp = decomposerPtr_().decompose
999 Info<<
" Redistributing background mesh cells" <<
endl;
1001 fvMeshDistribute distributor(mesh_);
1003 autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute(newDecomp);
1005 meshCutter_.distribute(mapDist);
1009 printMeshData(mesh_);
1011 Pout<<
" Post distribute sum(cellWeights) " 1012 <<
sum(icellWeights)
1013 <<
" max(cellWeights) " 1014 <<
max(icellWeights)
1020 cellWeights.
write();
1023 buildPatchAndTree();
1042 const List<point>& pts
1045 boolList posProc(pts.size(),
true);
1049 posProc[pI] = positionOnThisProcessor(pts[pI]);
1058 const treeBoundBox& box
1062 return !bFTreePtr_().findBox(box).empty();
1068 const point& centre,
1069 const scalar radiusSqr
1074 return bFTreePtr_().findNearest(centre, radiusSqr).hit();
1084 return bFTreePtr_().findLine(start, end);
1094 return bFTreePtr_().findLineAny(start, end);
1100 const List<point>& pts
1103 DynamicList<label> toCandidateProc;
1104 DynamicList<point> testPoints;
1108 label nTotalCandidates = 0;
1112 const point& pt = pts[pI];
1114 label nCandidates = 0;
1116 forAll(allBackgroundMeshBounds_, proci)
1120 if (allBackgroundMeshBounds_[proci].overlaps(pt,
sqr(small*100)))
1122 toCandidateProc.append(proci);
1123 testPoints.append(pt);
1129 ptBlockStart[pI] = nTotalCandidates;
1130 ptBlockSize[pI] = nCandidates;
1132 nTotalCandidates += nCandidates;
1136 label preDistributionToCandidateProcSize = toCandidateProc.size();
1138 autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1140 map().distribute(testPoints);
1142 List<scalar> distanceSqrToCandidate(testPoints.size(),
sqr(great));
1155 distanceSqrToCandidate[tPI] =
magSqr 1157 testPoints[tPI] - info.hitPoint()
1162 map().reverseDistribute
1164 preDistributionToCandidateProcSize,
1165 distanceSqrToCandidate
1168 labelList ptNearestProc(pts.size(), -1);
1174 SubList<scalar> ptNearestProcResults
1176 distanceSqrToCandidate,
1181 scalar nearestProcDistSqr = great;
1183 forAll(ptNearestProcResults, pPRI)
1185 if (ptNearestProcResults[pPRI] < nearestProcDistSqr)
1187 nearestProcDistSqr = ptNearestProcResults[pPRI];
1189 ptNearestProc[pI] = toCandidateProc[ptBlockStart[pI] + pPRI];
1195 Pout<< pts[pI] <<
" nearestProcDistSqr " << nearestProcDistSqr
1196 <<
" ptNearestProc[pI] " << ptNearestProc[pI] <<
endl;
1199 if (ptNearestProc[pI] < 0)
1202 <<
"The position " << pts[pI]
1203 <<
" did not find a nearest point on the background mesh." 1208 return ptNearestProc;
1216 const List<point>& starts,
1217 const List<point>& ends,
1218 bool includeOwnProcessor
1221 DynamicList<label> toCandidateProc;
1222 DynamicList<point> testStarts;
1223 DynamicList<point> testEnds;
1224 labelList segmentBlockStart(starts.size(), -1);
1225 labelList segmentBlockSize(starts.size(), -1);
1227 label nTotalCandidates = 0;
1231 const point& s = starts[sI];
1232 const point& e = ends[sI];
1237 label nCandidates = 0;
1239 forAll(allBackgroundMeshBounds_, proci)
1246 && allBackgroundMeshBounds_[proci].intersects(s, e,
p)
1249 toCandidateProc.append(proci);
1250 testStarts.append(s);
1257 segmentBlockStart[sI] = nTotalCandidates;
1258 segmentBlockSize[sI] = nCandidates;
1260 nTotalCandidates += nCandidates;
1264 label preDistributionToCandidateProcSize = toCandidateProc.size();
1266 autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1268 map().distribute(testStarts);
1269 map().distribute(testEnds);
1271 List<pointIndexHit> segmentIntersectsCandidate(testStarts.size());
1276 const point& s = testStarts[sI];
1277 const point& e = testEnds[sI];
1280 segmentIntersectsCandidate[sI] = bFTreePtr_().findLine(s, e);
1283 map().reverseDistribute
1285 preDistributionToCandidateProcSize,
1286 segmentIntersectsCandidate
1289 List<List<pointIndexHit>> segmentHitProcs(starts.size());
1292 DynamicList<pointIndexHit> tmpProcHits;
1296 tmpProcHits.clear();
1300 SubList<pointIndexHit> segmentProcResults
1302 segmentIntersectsCandidate,
1303 segmentBlockSize[sI],
1304 segmentBlockStart[sI]
1307 forAll(segmentProcResults, sPRI)
1309 if (segmentProcResults[sPRI].hit())
1311 tmpProcHits.append(segmentProcResults[sPRI]);
1313 tmpProcHits.last().setIndex
1315 toCandidateProc[segmentBlockStart[sI] + sPRI]
1320 segmentHitProcs[sI] = tmpProcHits;
1323 return segmentHitProcs;
1329 const point& centre,
1330 const scalar& radiusSqr
1333 forAll(allBackgroundMeshBounds_, proci)
1335 if (bFTreePtr_().findNearest(centre, radiusSqr).hit())
1347 const point& centre,
1348 const scalar radiusSqr
1353 forAll(allBackgroundMeshBounds_, proci)
1359 && allBackgroundMeshBounds_[proci].overlaps(centre, radiusSqr)
1365 toProc.append(proci);
1370 return Foam::move(toProc);
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.
virtual Ostream & write(const char)=0
Write character.
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
const dimensionSet dimless
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.
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< SubList< face >, 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.
vector point
Point is a vector.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
prefixOSstream Pout(cout, "Pout")
PrimitivePatch< faceList, const pointField > bPatch
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,.
backgroundMeshDecomposition(const Time &runTime, Random &rndGen, const conformationSurfaces &geometryToConformTo, const dictionary &coeffsDict)
Construct from components in foamyHexMesh operation.
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.