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),
786 "backgroundMeshDecomposition",
790 IOobject::AUTO_WRITE,
802 allBackgroundMeshBounds_(Pstream::nProcs()),
803 globalBackgroundBounds_(),
811 IOobject::MUST_READ_IF_MODIFIED,
815 decomposerPtr_(decompositionMethod::
New(decomposeDict_)),
816 mergeDist_(1e-6*mesh_.bounds().
mag()),
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_, mergeDist_);
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
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.
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,.
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.