46 const List<label>& toProc
57 label proci = toProc[i];
68 sendMap[proci].setSize(nSend[proci]);
76 label proci = toProc[i];
78 sendMap[proci][nSend[proci]++] = i;
99 forAll(constructMap, proci)
103 label nRecv = recvSizes[proci];
105 constructMap[proci].setSize(nRecv);
107 for (
label i = 0; i < nRecv; i++)
109 constructMap[proci][i] = constructSize++;
114 return autoPtr<mapDistribute>
128 void Foam::backgroundMeshDecomposition::initialRefinement()
135 mesh_.time().timeName(),
142 zeroGradientFvPatchScalarField::typeName
145 const conformationSurfaces& geometry = geometryToConformTo_;
147 decompositionMethod& decomposer = decomposerPtr_();
153 List<volumeType> volumeStatus
164 forAll(volumeStatus, celli)
170 mesh_.cells()[celli].points
177 if (geometry.overlaps(cellBb))
181 else if (geometry.inside(cellBb.midpoint()))
193 labelList refCells = selectRefinementCells
202 meshCutter_.consistentRefinement
209 forAll(newCellsToRefine, nCTRI)
211 label celli = newCellsToRefine[nCTRI];
218 icellWeights[celli] =
max 221 icellWeights[celli]/8.0
225 if (
returnReduce(newCellsToRefine.size(), sumOp<label>()) == 0)
231 polyTopoChange meshMod(mesh_);
234 meshCutter_.setRefinement(newCellsToRefine, meshMod);
237 autoPtr<mapPolyMesh> map = meshMod.changeMesh
247 mesh_.updateMesh(map);
250 meshCutter_.updateMesh(map);
255 const labelList& cellMap = map().cellMap();
257 List<volumeType> newVolumeStatus(cellMap.size());
261 label oldCelli = cellMap[newCelli];
269 newVolumeStatus[newCelli] = volumeStatus[oldCelli];
273 volumeStatus.transfer(newVolumeStatus);
276 Info<<
" Background mesh refined from " 278 <<
" to " << mesh_.globalData().nTotalCells()
279 <<
" cells." <<
endl;
283 forAll(volumeStatus, celli)
289 mesh_.cells()[celli].points
296 if (geometry.overlaps(cellBb))
300 else if (geometry.inside(cellBb.midpoint()))
312 bool removeOutsideCells =
false;
314 if (removeOutsideCells)
316 DynamicList<label> cellsToRemove;
318 forAll(volumeStatus, celli)
322 cellsToRemove.append(celli);
326 removeCells cellRemover(mesh_);
329 polyTopoChange meshMod(mesh_);
331 labelList exposedFaces = cellRemover.getExposedFaces
337 cellRemover.setRefinement
346 autoPtr<mapPolyMesh> map = meshMod.changeMesh
356 mesh_.updateMesh(map);
359 meshCutter_.updateMesh(map);
360 cellRemover.updateMesh(map);
365 const labelList& cellMap = map().cellMap();
367 List<volumeType> newVolumeStatus(cellMap.size());
371 label oldCelli = cellMap[newCelli];
379 newVolumeStatus[newCelli] = volumeStatus[oldCelli];
383 volumeStatus.transfer(newVolumeStatus);
388 - mesh_.globalData().nTotalCells()
389 <<
" cells." <<
endl;
401 labelList newDecomp = decomposer.decompose
408 fvMeshDistribute distributor(mesh_, mergeDist_);
410 autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute
415 meshCutter_.distribute(mapDist);
417 mapDist().distributeCellData(volumeStatus);
421 printMeshData(mesh_);
444 void Foam::backgroundMeshDecomposition::printMeshData
451 globalIndex globalCells(mesh.nCells());
476 Info<<
"Processor " << proci <<
" " 477 <<
"Number of cells = " << globalCells.localSize(proci)
501 bool Foam::backgroundMeshDecomposition::refineCell
505 scalar& weightEstimate
515 mesh_.cells()[celli].points
522 weightEstimate = 1.0;
626 Foam::labelList Foam::backgroundMeshDecomposition::selectRefinementCells
628 List<volumeType>& volumeStatus,
637 forAll(volumeStatus, celli)
641 if (meshCutter_.cellLevel()[celli] < minLevels_)
643 cellsToRefine.insert(celli);
659 cellsToRefine.insert(celli);
664 return cellsToRefine.toc();
668 void Foam::backgroundMeshDecomposition::buildPatchAndTree()
675 mesh_.nFaces() - mesh_.nInternalFaces(),
676 mesh_.nInternalFaces()
681 boundaryFacesPtr_.reset
685 tmpBoundaryFaces.localFaces(),
686 tmpBoundaryFaces.localPoints()
691 treeBoundBox overallBb(boundaryFacesPtr_().localPoints());
695 new indexedOctree<treeDataBPatch>
703 overallBb.extend(1e-4),
716 point bbMin(great, great, great);
717 point bbMax(-great, -great, -great);
719 forAll(allBackgroundMeshBounds_, proci)
721 bbMin =
min(bbMin, allBackgroundMeshBounds_[proci].
min());
722 bbMax =
max(bbMax, allBackgroundMeshBounds_[proci].
max());
725 globalBackgroundBounds_ = treeBoundBox(bbMin, bbMax);
732 /
"backgroundMeshDecomposition_proc_" 734 +
"_boundaryFaces.obj" 737 const faceList& faces = boundaryFacesPtr_().localFaces();
738 const List<point>& points = boundaryFacesPtr_().localPoints();
740 Map<label> foamToObj(points.size());
746 const face&
f = faces[i];
750 if (foamToObj.insert(f[fPI], vertI))
761 fStr<<
' ' << foamToObj[f[fPI]] + 1;
776 const conformationSurfaces& geometryToConformTo,
777 const dictionary& coeffsDict
781 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 mergeDist_(1e-6*mesh_.bounds().
mag()),
821 coeffsDict.lookupOrDefault<scalar>(
"minCellSizeLimit", 0.0)
830 <<
"This cannot be used when not running in parallel." 834 if (!decomposerPtr_().parallelAware())
837 <<
"You have selected decomposition method " 838 << decomposerPtr_().typeName
839 <<
" which is not parallel aware." << endl
843 Info<<
nl <<
"Building initial background mesh decomposition" <<
endl;
877 label nOccupiedCells = 0;
881 if (icellWeights[cI] > 1 - small)
891 scalar cellWeightLimit =
max 894 *
sum(cellWeights).value()
901 Info<<
" cellWeightLimit " << cellWeightLimit <<
endl;
903 Pout<<
" sum(cellWeights) " <<
sum(cellWeights.primitiveField())
904 <<
" max(cellWeights) " <<
max(cellWeights.primitiveField())
912 if (icellWeights[cWI] > cellWeightLimit)
914 cellsToRefine.insert(cWI);
918 if (
returnReduce(cellsToRefine.size(), sumOp<label>()) == 0)
926 meshCutter_.consistentRefinement
933 if (debug && !cellsToRefine.empty())
935 Pout<<
" cellWeights too large in " << cellsToRefine.size()
939 forAll(newCellsToRefine, nCTRI)
941 label celli = newCellsToRefine[nCTRI];
943 icellWeights[celli] /= 8.0;
947 polyTopoChange meshMod(mesh_);
950 meshCutter_.setRefinement(newCellsToRefine, meshMod);
953 autoPtr<mapPolyMesh> map = meshMod.changeMesh
963 mesh_.updateMesh(map);
966 meshCutter_.updateMesh(map);
968 Info<<
" Background mesh refined from " 970 <<
" to " << mesh_.globalData().nTotalCells()
971 <<
" cells." <<
endl;
984 printMeshData(mesh_);
986 Pout<<
" Pre distribute sum(cellWeights) " 988 <<
" max(cellWeights) " 993 labelList newDecomp = decomposerPtr_().decompose
1000 Info<<
" Redistributing background mesh cells" <<
endl;
1002 fvMeshDistribute distributor(mesh_, mergeDist_);
1004 autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute(newDecomp);
1006 meshCutter_.distribute(mapDist);
1010 printMeshData(mesh_);
1012 Pout<<
" Post distribute sum(cellWeights) " 1013 <<
sum(icellWeights)
1014 <<
" max(cellWeights) " 1015 <<
max(icellWeights)
1021 cellWeights.
write();
1024 buildPatchAndTree();
1043 const List<point>& pts
1046 boolList posProc(pts.size(),
true);
1050 posProc[pI] = positionOnThisProcessor(pts[pI]);
1059 const treeBoundBox& box
1063 return !bFTreePtr_().findBox(box).empty();
1069 const point& centre,
1070 const scalar radiusSqr
1075 return bFTreePtr_().findNearest(centre, radiusSqr).hit();
1085 return bFTreePtr_().findLine(start, end);
1095 return bFTreePtr_().findLineAny(start, end);
1101 const List<point>& pts
1104 DynamicList<label> toCandidateProc;
1105 DynamicList<point> testPoints;
1109 label nTotalCandidates = 0;
1113 const point& pt = pts[pI];
1115 label nCandidates = 0;
1117 forAll(allBackgroundMeshBounds_, proci)
1121 if (allBackgroundMeshBounds_[proci].overlaps(pt,
sqr(small*100)))
1123 toCandidateProc.append(proci);
1124 testPoints.append(pt);
1130 ptBlockStart[pI] = nTotalCandidates;
1131 ptBlockSize[pI] = nCandidates;
1133 nTotalCandidates += nCandidates;
1137 label preDistributionToCandidateProcSize = toCandidateProc.size();
1139 autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1141 map().distribute(testPoints);
1143 List<scalar> distanceSqrToCandidate(testPoints.size(),
sqr(great));
1156 distanceSqrToCandidate[tPI] =
magSqr 1158 testPoints[tPI] - info.hitPoint()
1163 map().reverseDistribute
1165 preDistributionToCandidateProcSize,
1166 distanceSqrToCandidate
1169 labelList ptNearestProc(pts.size(), -1);
1175 SubList<scalar> ptNearestProcResults
1177 distanceSqrToCandidate,
1182 scalar nearestProcDistSqr = great;
1184 forAll(ptNearestProcResults, pPRI)
1186 if (ptNearestProcResults[pPRI] < nearestProcDistSqr)
1188 nearestProcDistSqr = ptNearestProcResults[pPRI];
1190 ptNearestProc[pI] = toCandidateProc[ptBlockStart[pI] + pPRI];
1196 Pout<< pts[pI] <<
" nearestProcDistSqr " << nearestProcDistSqr
1197 <<
" ptNearestProc[pI] " << ptNearestProc[pI] <<
endl;
1200 if (ptNearestProc[pI] < 0)
1203 <<
"The position " << pts[pI]
1204 <<
" did not find a nearest point on the background mesh." 1209 return ptNearestProc;
1217 const List<point>& starts,
1218 const List<point>& ends,
1219 bool includeOwnProcessor
1222 DynamicList<label> toCandidateProc;
1223 DynamicList<point> testStarts;
1224 DynamicList<point> testEnds;
1225 labelList segmentBlockStart(starts.size(), -1);
1226 labelList segmentBlockSize(starts.size(), -1);
1228 label nTotalCandidates = 0;
1232 const point& s = starts[sI];
1233 const point& e = ends[sI];
1238 label nCandidates = 0;
1240 forAll(allBackgroundMeshBounds_, proci)
1247 && allBackgroundMeshBounds_[proci].intersects(s, e,
p)
1250 toCandidateProc.append(proci);
1251 testStarts.append(s);
1258 segmentBlockStart[sI] = nTotalCandidates;
1259 segmentBlockSize[sI] = nCandidates;
1261 nTotalCandidates += nCandidates;
1265 label preDistributionToCandidateProcSize = toCandidateProc.size();
1267 autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1269 map().distribute(testStarts);
1270 map().distribute(testEnds);
1272 List<pointIndexHit> segmentIntersectsCandidate(testStarts.size());
1277 const point& s = testStarts[sI];
1278 const point& e = testEnds[sI];
1281 segmentIntersectsCandidate[sI] = bFTreePtr_().findLine(s, e);
1284 map().reverseDistribute
1286 preDistributionToCandidateProcSize,
1287 segmentIntersectsCandidate
1290 List<List<pointIndexHit>> segmentHitProcs(starts.size());
1293 DynamicList<pointIndexHit> tmpProcHits;
1297 tmpProcHits.clear();
1301 SubList<pointIndexHit> segmentProcResults
1303 segmentIntersectsCandidate,
1304 segmentBlockSize[sI],
1305 segmentBlockStart[sI]
1308 forAll(segmentProcResults, sPRI)
1310 if (segmentProcResults[sPRI].hit())
1312 tmpProcHits.append(segmentProcResults[sPRI]);
1314 tmpProcHits.last().setIndex
1316 toCandidateProc[segmentBlockStart[sI] + sPRI]
1321 segmentHitProcs[sI] = tmpProcHits;
1324 return segmentHitProcs;
1330 const point& centre,
1331 const scalar& radiusSqr
1334 forAll(allBackgroundMeshBounds_, proci)
1336 if (bFTreePtr_().findNearest(centre, radiusSqr).hit())
1348 const point& centre,
1349 const scalar radiusSqr
1354 forAll(allBackgroundMeshBounds_, proci)
1360 && allBackgroundMeshBounds_[proci].overlaps(centre, radiusSqr)
1366 toProc.append(proci);
1371 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.
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.
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< 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.
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 > &)
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,.
virtual Ostream & write(const token &)=0
Write next token to stream.
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.