50 const List<label>& toProc
61 label proci = toProc[i];
71 sendMap[proci].resize_nocopy(nSend[proci]);
78 label proci = toProc[i];
79 sendMap[proci][nSend[proci]++] = i;
88 void Foam::backgroundMeshDecomposition::initialRefinement()
95 mesh_.time().timeName(),
105 const conformationSurfaces& geometry = geometryToConformTo_;
107 decompositionMethod& decomposer =
114 List<volumeType> volumeStatus
125 forAll(volumeStatus, celli)
129 treeBoundBox cellBb(mesh_.cellBb(celli));
131 if (geometry.overlaps(cellBb))
135 else if (geometry.inside(cellBb.centre()))
147 labelList refCells = selectRefinementCells
156 meshCutter_.consistentRefinement
163 forAll(newCellsToRefine, nCTRI)
165 label celli = newCellsToRefine[nCTRI];
172 icellWeights[celli] =
max 175 icellWeights[celli]/8.0
185 polyTopoChange meshMod(mesh_);
188 meshCutter_.setRefinement(newCellsToRefine, meshMod);
191 autoPtr<mapPolyMesh> map = meshMod.changeMesh
201 mesh_.updateMesh(map());
204 meshCutter_.updateMesh(map());
209 const labelList& cellMap = map().cellMap();
211 List<volumeType> newVolumeStatus(cellMap.size());
215 label oldCelli = cellMap[newCelli];
223 newVolumeStatus[newCelli] = volumeStatus[oldCelli];
227 volumeStatus.transfer(newVolumeStatus);
230 Info<<
" Background mesh refined from " 232 <<
" to " << mesh_.globalData().nTotalCells()
233 <<
" cells." <<
endl;
237 forAll(volumeStatus, celli)
241 treeBoundBox cellBb(mesh_.cellBb(celli));
243 if (geometry.overlaps(cellBb))
247 else if (geometry.inside(cellBb.centre()))
259 bool removeOutsideCells =
false;
261 if (removeOutsideCells)
263 DynamicList<label> cellsToRemove;
265 forAll(volumeStatus, celli)
269 cellsToRemove.append(celli);
273 removeCells cellRemover(mesh_);
276 polyTopoChange meshMod(mesh_);
278 labelList exposedFaces = cellRemover.getExposedFaces
284 cellRemover.setRefinement
293 autoPtr<mapPolyMesh> map = meshMod.changeMesh
303 mesh_.updateMesh(map());
306 meshCutter_.updateMesh(map());
307 cellRemover.updateMesh(map());
312 const labelList& cellMap = map().cellMap();
314 List<volumeType> newVolumeStatus(cellMap.size());
318 label oldCelli = cellMap[newCelli];
326 newVolumeStatus[newCelli] = volumeStatus[oldCelli];
330 volumeStatus.transfer(newVolumeStatus);
335 - mesh_.globalData().nTotalCells()
336 <<
" cells." <<
endl;
348 labelList newDecomp = decomposer.decompose
355 fvMeshDistribute distributor(mesh_);
357 autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute
362 meshCutter_.distribute(mapDist());
364 mapDist().distributeCellData(volumeStatus);
368 printMeshData(mesh_);
391 void Foam::backgroundMeshDecomposition::printMeshData
423 Info<<
"Processor " << proci <<
" " 424 <<
"Number of cells = " << globalCells.localSize(proci)
448 bool Foam::backgroundMeshDecomposition::refineCell
452 scalar& weightEstimate
460 treeBoundBox cellBb(mesh_.cellBb(celli));
462 weightEstimate = 1.0;
565 Foam::labelList Foam::backgroundMeshDecomposition::selectRefinementCells
567 List<volumeType>& volumeStatus,
576 forAll(volumeStatus, celli)
580 if (meshCutter_.cellLevel()[celli] < minLevels_)
582 cellsToRefine.insert(celli);
598 cellsToRefine.insert(celli);
603 return cellsToRefine.toc();
607 void Foam::backgroundMeshDecomposition::buildPatchAndTree()
614 mesh_.nBoundaryFaces(),
615 mesh_.nInternalFaces()
620 boundaryFacesPtr_.reset
624 tmpBoundaryFaces.localFaces(),
625 tmpBoundaryFaces.localPoints()
630 treeBoundBox overallBb(boundaryFacesPtr_().localPoints());
632 Random& rnd = rndGen_;
636 new indexedOctree<treeDataBPatch>
644 overallBb.extend(rnd, 1
e-4),
657 globalBackgroundBounds_.reset();
658 forAll(allBackgroundMeshBounds_, proci)
660 globalBackgroundBounds_.add(allBackgroundMeshBounds_[proci]);
668 /
"backgroundMeshDecomposition_proc_" 670 +
"_boundaryFaces.obj" 673 const faceList& faces = boundaryFacesPtr_().localFaces();
674 const List<point>&
points = boundaryFacesPtr_().localPoints();
682 const face&
f = faces[i];
686 if (foamToObj.insert(
f[fPI], vertI))
697 fStr<<
' ' << foamToObj[
f[fPI]] + 1;
708 Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
712 const conformationSurfaces& geometryToConformTo,
713 const dictionary& coeffsDict,
714 const fileName& decompDictFile
718 geometryToConformTo_(geometryToConformTo),
724 "backgroundMeshDecomposition",
728 IOobject::AUTO_WRITE,
729 IOobject::NO_REGISTER
740 allBackgroundMeshBounds_(Pstream::nProcs()),
741 globalBackgroundBounds_(),
742 mergeDist_(1
e-6*mesh_.bounds().
mag()),
743 spanScale_(coeffsDict.
get<scalar>(
"spanScale")),
746 coeffsDict.getOrDefault<scalar>(
"minCellSizeLimit", 0)
748 minLevels_(coeffsDict.
get<label>(
"minLevels")),
749 volRes_(coeffsDict.
get<label>(
"sampleResolution")),
750 maxCellWeightCoeff_(coeffsDict.
get<scalar>(
"maxCellWeightCoeff"))
755 <<
"This cannot be used when not running in parallel." 759 const decompositionMethod& decomposer =
762 if (!decomposer.parallelAware())
765 <<
"You have selected decomposition method " 766 << decomposer.typeName
767 <<
" which is not parallel aware." <<
endl 771 Info<<
nl <<
"Building initial background mesh decomposition" <<
endl;
799 label nOccupiedCells = 0;
803 if (icellWeights[cI] > 1 - SMALL)
813 scalar cellWeightLimit =
max 816 *
sum(cellWeights).value()
823 Info<<
" cellWeightLimit " << cellWeightLimit <<
endl;
825 Pout<<
" sum(cellWeights) " <<
sum(cellWeights.primitiveField())
826 <<
" max(cellWeights) " <<
max(cellWeights.primitiveField())
834 if (icellWeights[cWI] > cellWeightLimit)
836 cellsToRefine.insert(cWI);
848 meshCutter_.consistentRefinement
855 if (
debug && !cellsToRefine.empty())
857 Pout<<
" cellWeights too large in " << cellsToRefine.size()
861 forAll(newCellsToRefine, nCTRI)
863 label celli = newCellsToRefine[nCTRI];
865 icellWeights[celli] /= 8.0;
869 polyTopoChange meshMod(mesh_);
872 meshCutter_.setRefinement(newCellsToRefine, meshMod);
875 autoPtr<mapPolyMesh> map = meshMod.changeMesh
885 mesh_.updateMesh(map());
888 meshCutter_.updateMesh(map());
890 Info<<
" Background mesh refined from " 892 <<
" to " << mesh_.globalData().nTotalCells()
893 <<
" cells." <<
endl;
906 printMeshData(mesh_);
908 Pout<<
" Pre distribute sum(cellWeights) " 910 <<
" max(cellWeights) " 915 decompositionMethod& decomposer =
918 labelList newDecomp = decomposer.decompose
925 Info<<
" Redistributing background mesh cells" <<
endl;
927 fvMeshDistribute distributor(mesh_);
929 autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute(newDecomp);
931 meshCutter_.distribute(mapDist());
935 printMeshData(mesh_);
937 Pout<<
" Post distribute sum(cellWeights) " 939 <<
" max(cellWeights) " 968 const List<point>&
pts 975 posProc[pI] = positionOnThisProcessor(
pts[pI]);
984 const treeBoundBox& box
988 return !bFTreePtr_().findBox(box).empty();
995 const scalar radiusSqr
1000 return bFTreePtr_().findNearest(centre, radiusSqr).hit();
1010 return bFTreePtr_().findLine(start,
end);
1020 return bFTreePtr_().findLineAny(start,
end);
1026 const List<point>&
pts 1029 DynamicList<label> toCandidateProc;
1030 DynamicList<point> testPoints;
1034 label nTotalCandidates = 0;
1040 label nCandidates = 0;
1042 forAll(allBackgroundMeshBounds_, proci)
1046 if (allBackgroundMeshBounds_[proci].overlaps(pt,
sqr(SMALL*100)))
1048 toCandidateProc.append(proci);
1049 testPoints.append(pt);
1055 ptBlockStart[pI] = nTotalCandidates;
1056 ptBlockSize[pI] = nCandidates;
1058 nTotalCandidates += nCandidates;
1062 label preDistributionToCandidateProcSize = toCandidateProc.size();
1064 autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1066 map().distribute(testPoints);
1068 List<scalar> distanceSqrToCandidate(testPoints.size(),
sqr(GREAT));
1081 distanceSqrToCandidate[tPI] = info.point().distSqr(testPoints[tPI]);
1085 map().reverseDistribute
1087 preDistributionToCandidateProcSize,
1088 distanceSqrToCandidate
1097 SubList<scalar> ptNearestProcResults
1099 distanceSqrToCandidate,
1104 scalar nearestProcDistSqr = GREAT;
1106 forAll(ptNearestProcResults, pPRI)
1108 if (ptNearestProcResults[pPRI] < nearestProcDistSqr)
1110 nearestProcDistSqr = ptNearestProcResults[pPRI];
1112 ptNearestProc[pI] = toCandidateProc[ptBlockStart[pI] + pPRI];
1118 Pout<<
pts[pI] <<
" nearestProcDistSqr " << nearestProcDistSqr
1119 <<
" ptNearestProc[pI] " << ptNearestProc[pI] <<
endl;
1122 if (ptNearestProc[pI] < 0)
1125 <<
"The position " <<
pts[pI]
1126 <<
" did not find a nearest point on the background mesh." 1131 return ptNearestProc;
1139 const List<point>& starts,
1140 const List<point>& ends,
1141 bool includeOwnProcessor
1144 DynamicList<label> toCandidateProc;
1145 DynamicList<point> testStarts;
1146 DynamicList<point> testEnds;
1147 labelList segmentBlockStart(starts.size(), -1);
1148 labelList segmentBlockSize(starts.size(), -1);
1150 label nTotalCandidates = 0;
1154 const point&
s = starts[sI];
1155 const point&
e = ends[sI];
1160 label nCandidates = 0;
1162 forAll(allBackgroundMeshBounds_, proci)
1169 && allBackgroundMeshBounds_[proci].intersects(
s,
e,
p)
1172 toCandidateProc.append(proci);
1173 testStarts.append(
s);
1180 segmentBlockStart[sI] = nTotalCandidates;
1181 segmentBlockSize[sI] = nCandidates;
1183 nTotalCandidates += nCandidates;
1187 label preDistributionToCandidateProcSize = toCandidateProc.size();
1189 autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1191 map().distribute(testStarts);
1192 map().distribute(testEnds);
1194 List<pointIndexHit> segmentIntersectsCandidate(testStarts.size());
1199 const point&
s = testStarts[sI];
1200 const point&
e = testEnds[sI];
1203 segmentIntersectsCandidate[sI] = bFTreePtr_().findLine(
s,
e);
1206 map().reverseDistribute
1208 preDistributionToCandidateProcSize,
1209 segmentIntersectsCandidate
1212 List<List<pointIndexHit>> segmentHitProcs(starts.size());
1215 DynamicList<pointIndexHit> tmpProcHits;
1219 tmpProcHits.clear();
1223 SubList<pointIndexHit> segmentProcResults
1225 segmentIntersectsCandidate,
1226 segmentBlockSize[sI],
1227 segmentBlockStart[sI]
1230 forAll(segmentProcResults, sPRI)
1232 if (segmentProcResults[sPRI].hit())
1234 tmpProcHits.append(segmentProcResults[sPRI]);
1236 tmpProcHits.last().setIndex
1238 toCandidateProc[segmentBlockStart[sI] + sPRI]
1243 segmentHitProcs[sI] = tmpProcHits;
1246 return segmentHitProcs;
1252 const point& centre,
1253 const scalar& radiusSqr
1256 forAll(allBackgroundMeshBounds_, proci)
1258 if (bFTreePtr_().findNearest(centre, radiusSqr).hit())
1270 const point& centre,
1271 const scalar radiusSqr
1276 forAll(allBackgroundMeshBounds_, proci)
1282 && allBackgroundMeshBounds_[proci].overlaps(centre, radiusSqr)
1288 toProc.append(proci);
static const word & zeroGradientType() noexcept
The type name for zeroGradient patch fields.
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
List of values generated by applying the access operation to each list item.
static scalar & perturbTol() noexcept
Get the perturbation tolerance.
void size(const label n)
Older name for setAddressableSize.
bool positionOnThisProcessor(const point &pt) const
Is the given position inside the domain of this decomposition.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
PrimitivePatch<::Foam::List< face >, const pointField > bPatch
errorManipArg< error, int > exit(error &err, const int errNo=1)
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.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static rangeType allProcs(const label communicator=worldComm)
Range of process indices for all processes.
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
PointIndexHit< point > pointIndexHit
A PointIndexHit with a 3D point.
static bool & parRun() noexcept
Test if this a parallel run.
This class describes the interaction of an object (often a face) and a point. It carries the info of ...
labelList overlapProcessors(const point ¢re, const scalar radiusSqr) const
Ignore writing from objectRegistry::writeObject()
const dimensionSet dimless
Dimensionless.
static int myProcNo(const label communicator=worldComm)
Rank of this process in the communicator (starting from masterNo()). Can be negative if the process i...
static const decompositionModel & New(const polyMesh &mesh, const fileName &decompDictFile="", const dictionary *fallback=nullptr)
Read and register on mesh, optionally with alternative decomposeParDict path/name or with fallback co...
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.
List< labelList > labelListList
List of labelList.
DimensionedField< scalar, volMesh > Internal
The internal field type from which this GeometricField is derived.
decompositionMethod & decomposer() const
Return demand-driven decomposition method.
static autoPtr< mapDistribute > buildMap(const List< label > &toProc)
Build a mapDistribute for the supplied destination processor data.
#define forAll(list, i)
Loop across all elements in list.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
static void allGatherList(List< T > &values, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Gather data, but keep individual values separate. Uses linear/tree communication. ...
List< face > faceList
List of faces.
treeDataPrimitivePatch< bPatch > treeDataBPatch
labelList processorNearestPosition(const List< point > &pts) const
What is the nearest processor to the given position?
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator). It is 1 for serial run. ...
bool returnReduceAnd(const bool value, const label comm=UPstream::worldComm)
Perform logical (and) MPI Allreduce on a copy. Uses UPstream::reduceAnd.
const dimensionedScalar e
Elementary charge.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
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.
static const word null
An empty word.
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.
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
A location inside the volume.
A location outside the volume.
A location that is partly inside and outside.
int debug
Static debugging option.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
defineTypeNameAndDebug(combustionModel, 0)
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end, (exposing.
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field...
vector point
Point is a vector.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label nCells() const noexcept
Number of mesh cells.
messageStream Info
Information stream (stdout output on master, null elsewhere)
autoPtr< mapDistributePolyMesh > distribute(volScalarField &cellWeights)
Redistribute the background mesh based on a supplied weight field,.
List< label > labelList
A List of labels.
static autoPtr< T > New(Args &&... args)
Construct autoPtr with forwarding arguments.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
List< bool > boolList
A List of bools.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
bool overlapsThisProcessor(const treeBoundBox &box) const
Does the given box overlap the faces of the boundary of this.
static constexpr const zero Zero
Global zero (0)