50 const List<label>& toProc
61 label proci = toProc[i];
72 sendMap[proci].setSize(nSend[proci]);
80 label proci = toProc[i];
82 sendMap[proci][nSend[proci]++] = i;
103 forAll(constructMap, proci)
107 label nRecv = recvSizes[proci];
109 constructMap[proci].setSize(nRecv);
111 for (label i = 0; i < nRecv; i++)
113 constructMap[proci][i] = constructSize++;
122 std::move(constructMap)
129void Foam::backgroundMeshDecomposition::initialRefinement()
136 mesh_.time().timeName(),
143 zeroGradientFvPatchScalarField::typeName
146 const conformationSurfaces& geometry = geometryToConformTo_;
148 decompositionMethod& decomposer =
155 List<volumeType> volumeStatus
166 forAll(volumeStatus, celli)
172 mesh_.cells()[celli].points
179 if (geometry.overlaps(cellBb))
183 else if (geometry.inside(cellBb.centre()))
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.centre()))
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_);
412 autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute
417 meshCutter_.distribute(mapDist());
419 mapDist().distributeCellData(volumeStatus);
423 printMeshData(mesh_);
446void Foam::backgroundMeshDecomposition::printMeshData
478 Info<<
"Processor " << proci <<
" "
479 <<
"Number of cells = " << globalCells.localSize(proci)
507 scalar& weightEstimate
517 mesh_.cells()[celli].points
524 weightEstimate = 1.0;
627Foam::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();
669void Foam::backgroundMeshDecomposition::buildPatchAndTree()
676 mesh_.nBoundaryFaces(),
677 mesh_.nInternalFaces()
682 boundaryFacesPtr_.reset
686 tmpBoundaryFaces.localFaces(),
687 tmpBoundaryFaces.localPoints()
692 treeBoundBox overallBb(boundaryFacesPtr_().localPoints());
694 Random& rnd = rndGen_;
698 new indexedOctree<treeDataBPatch>
706 overallBb.extend(rnd, 1
e-4),
720 forAll(allBackgroundMeshBounds_, proci)
722 globalBackgroundBounds_.add(allBackgroundMeshBounds_[proci]);
730 /
"backgroundMeshDecomposition_proc_"
732 +
"_boundaryFaces.obj"
735 const faceList& faces = boundaryFacesPtr_().localFaces();
736 const List<point>&
points = boundaryFacesPtr_().localPoints();
744 const face&
f = faces[i];
748 if (foamToObj.insert(
f[fPI], vertI))
759 fStr<<
' ' << foamToObj[
f[fPI]] + 1;
774 const conformationSurfaces& geometryToConformTo,
775 const dictionary& coeffsDict,
776 const fileName& decompDictFile
780 geometryToConformTo_(geometryToConformTo),
786 "backgroundMeshDecomposition",
790 IOobject::AUTO_WRITE,
802 allBackgroundMeshBounds_(Pstream::nProcs()),
803 globalBackgroundBounds_(),
804 mergeDist_(1
e-6*mesh_.bounds().
mag()),
805 spanScale_(coeffsDict.
get<scalar>(
"spanScale")),
808 coeffsDict.getOrDefault<scalar>(
"minCellSizeLimit", 0)
810 minLevels_(coeffsDict.
get<label>(
"minLevels")),
811 volRes_(coeffsDict.
get<label>(
"sampleResolution")),
812 maxCellWeightCoeff_(coeffsDict.
get<scalar>(
"maxCellWeightCoeff"))
814 if (!Pstream::parRun())
817 <<
"This cannot be used when not running in parallel."
821 const decompositionMethod& decomposer =
822 decompositionModel::New(mesh_, decompDictFile).decomposer();
824 if (!decomposer.parallelAware())
827 <<
"You have selected decomposition method "
828 << decomposer.typeName
829 <<
" which is not parallel aware." <<
endl
833 Info<<
nl <<
"Building initial background mesh decomposition" <<
endl;
844 volScalarField& cellWeights
861 label nOccupiedCells = 0;
865 if (icellWeights[cI] > 1 - SMALL)
875 scalar cellWeightLimit =
max
878 *
sum(cellWeights).value()
885 Info<<
" cellWeightLimit " << cellWeightLimit <<
endl;
887 Pout<<
" sum(cellWeights) " <<
sum(cellWeights.primitiveField())
888 <<
" max(cellWeights) " <<
max(cellWeights.primitiveField())
896 if (icellWeights[cWI] > cellWeightLimit)
898 cellsToRefine.insert(cWI);
902 if (
returnReduce(cellsToRefine.size(), sumOp<label>()) == 0)
917 if (debug && !cellsToRefine.empty())
919 Pout<<
" cellWeights too large in " << cellsToRefine.size()
923 forAll(newCellsToRefine, nCTRI)
925 label celli = newCellsToRefine[nCTRI];
927 icellWeights[celli] /= 8.0;
931 polyTopoChange meshMod(mesh_);
937 autoPtr<mapPolyMesh> map = meshMod.changeMesh
952 Info<<
" Background mesh refined from "
955 <<
" cells." <<
endl;
968 printMeshData(mesh_);
970 Pout<<
" Pre distribute sum(cellWeights) "
972 <<
" max(cellWeights) "
977 decompositionMethod& decomposer =
980 labelList newDecomp = decomposer.decompose
987 Info<<
" Redistributing background mesh cells" <<
endl;
989 fvMeshDistribute distributor(mesh_);
991 autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute(newDecomp);
997 printMeshData(mesh_);
999 Pout<<
" Post distribute sum(cellWeights) "
1000 <<
sum(icellWeights)
1001 <<
" max(cellWeights) "
1002 <<
max(icellWeights)
1008 cellWeights.write();
1011 buildPatchAndTree();
1030 const List<point>& pts
1033 boolList posProc(pts.size(),
true);
1037 posProc[pI] = positionOnThisProcessor(pts[pI]);
1046 const treeBoundBox& box
1050 return !bFTreePtr_().findBox(box).empty();
1056 const point& centre,
1057 const scalar radiusSqr
1062 return bFTreePtr_().findNearest(centre, radiusSqr).hit();
1072 return bFTreePtr_().findLine(start, end);
1082 return bFTreePtr_().findLineAny(start, end);
1088 const List<point>& pts
1091 DynamicList<label> toCandidateProc;
1092 DynamicList<point> testPoints;
1096 label nTotalCandidates = 0;
1100 const point& pt = pts[pI];
1102 label nCandidates = 0;
1104 forAll(allBackgroundMeshBounds_, proci)
1108 if (allBackgroundMeshBounds_[proci].overlaps(pt,
sqr(SMALL*100)))
1110 toCandidateProc.append(proci);
1111 testPoints.append(pt);
1117 ptBlockStart[pI] = nTotalCandidates;
1118 ptBlockSize[pI] = nCandidates;
1120 nTotalCandidates += nCandidates;
1124 label preDistributionToCandidateProcSize = toCandidateProc.size();
1126 autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1128 map().distribute(testPoints);
1130 List<scalar> distanceSqrToCandidate(testPoints.size(),
sqr(GREAT));
1143 distanceSqrToCandidate[tPI] =
magSqr
1145 testPoints[tPI] - info.hitPoint()
1150 map().reverseDistribute
1152 preDistributionToCandidateProcSize,
1153 distanceSqrToCandidate
1156 labelList ptNearestProc(pts.size(), -1);
1162 SubList<scalar> ptNearestProcResults
1164 distanceSqrToCandidate,
1169 scalar nearestProcDistSqr = GREAT;
1171 forAll(ptNearestProcResults, pPRI)
1173 if (ptNearestProcResults[pPRI] < nearestProcDistSqr)
1175 nearestProcDistSqr = ptNearestProcResults[pPRI];
1177 ptNearestProc[pI] = toCandidateProc[ptBlockStart[pI] + pPRI];
1183 Pout<< pts[pI] <<
" nearestProcDistSqr " << nearestProcDistSqr
1184 <<
" ptNearestProc[pI] " << ptNearestProc[pI] <<
endl;
1187 if (ptNearestProc[pI] < 0)
1190 <<
"The position " << pts[pI]
1191 <<
" did not find a nearest point on the background mesh."
1196 return ptNearestProc;
1204 const List<point>& starts,
1205 const List<point>& ends,
1206 bool includeOwnProcessor
1209 DynamicList<label> toCandidateProc;
1210 DynamicList<point> testStarts;
1211 DynamicList<point> testEnds;
1212 labelList segmentBlockStart(starts.size(), -1);
1213 labelList segmentBlockSize(starts.size(), -1);
1215 label nTotalCandidates = 0;
1219 const point&
s = starts[sI];
1220 const point&
e = ends[sI];
1225 label nCandidates = 0;
1227 forAll(allBackgroundMeshBounds_, proci)
1234 && allBackgroundMeshBounds_[proci].intersects(
s,
e,
p)
1237 toCandidateProc.append(proci);
1238 testStarts.append(
s);
1245 segmentBlockStart[sI] = nTotalCandidates;
1246 segmentBlockSize[sI] = nCandidates;
1248 nTotalCandidates += nCandidates;
1252 label preDistributionToCandidateProcSize = toCandidateProc.size();
1254 autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1256 map().distribute(testStarts);
1257 map().distribute(testEnds);
1259 List<pointIndexHit> segmentIntersectsCandidate(testStarts.size());
1264 const point&
s = testStarts[sI];
1265 const point&
e = testEnds[sI];
1268 segmentIntersectsCandidate[sI] = bFTreePtr_().findLine(
s,
e);
1271 map().reverseDistribute
1273 preDistributionToCandidateProcSize,
1274 segmentIntersectsCandidate
1277 List<List<pointIndexHit>> segmentHitProcs(starts.size());
1280 DynamicList<pointIndexHit> tmpProcHits;
1284 tmpProcHits.clear();
1288 SubList<pointIndexHit> segmentProcResults
1290 segmentIntersectsCandidate,
1291 segmentBlockSize[sI],
1292 segmentBlockStart[sI]
1295 forAll(segmentProcResults, sPRI)
1297 if (segmentProcResults[sPRI].hit())
1299 tmpProcHits.append(segmentProcResults[sPRI]);
1301 tmpProcHits.last().setIndex
1303 toCandidateProc[segmentBlockStart[sI] + sPRI]
1308 segmentHitProcs[sI] = tmpProcHits;
1311 return segmentHitProcs;
1317 const point& centre,
1318 const scalar& radiusSqr
1321 forAll(allBackgroundMeshBounds_, proci)
1323 if (bFTreePtr_().findNearest(centre, radiusSqr).hit())
1335 const point& centre,
1336 const scalar radiusSqr
1341 forAll(allBackgroundMeshBounds_, proci)
1347 && allBackgroundMeshBounds_[proci].overlaps(centre, radiusSqr)
1353 toProc.append(proci);
DimensionedField< scalar, volMesh > Internal
The internal field type from which this GeometricField is derived.
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
UPstream::rangeType allProcs() const noexcept
Range of ranks indices associated with PstreamBuffers.
label nProcs() const noexcept
Number of ranks associated with PstreamBuffers.
static void allGatherList(const List< commsStruct > &comms, List< T > &values, const int tag, const label comm)
static void exchangeSizes(const labelUList &sendProcs, const labelUList &recvProcs, const Container &sendData, labelList &sizes, const label tag=UPstream::msgType(), const label comm=UPstream::worldComm)
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
void size(const label n)
Older name for setAddressableSize.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Store a background polyMesh to use for the decomposition of space and queries for parallel conformalV...
labelList processorNearestPosition(const List< point > &pts) const
What is the nearest processor to the given position?
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end, (exposing.
bool overlapsThisProcessor(const treeBoundBox &box) const
Does the given box overlap the faces of the boundary of this.
bool positionOnThisProcessor(const point &pt) const
Is the given position inside the domain of this decomposition.
autoPtr< mapDistributePolyMesh > distribute(volScalarField &cellWeights)
Redistribute the background mesh based on a supplied weight field,.
static autoPtr< mapDistribute > buildMap(const List< label > &toProc)
Build a mapDistribute for the supplied destination processor data.
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.
pointIndexHit findLine(const point &start, const point &end) const
Find nearest intersection of line between start and end, (exposing.
bool overlapsOtherProcessors(const point ¢re, const scalar &radiusSqr) const
labelList overlapProcessors(const point ¢re, const scalar radiusSqr) const
static const boundBox invertedBox
A large inverted boundBox: min/max == +/- ROOTVGREAT.
virtual bool write(const bool valid=true) const
Write mesh using IO settings from time.
virtual void updateMesh(const mapPolyMesh &mpm)
Update mesh corresponding to the given map.
label nTotalCells() const noexcept
Return total number of cells in decomposed mesh.
labelListList setRefinement(const labelList &cells, polyTopoChange &)
Insert refinement. All selected cells will be split into 8.
labelList consistentRefinement(const labelUList &cellLevel, const labelList &cellsToRefine, const bool maxSet) const
Given valid mesh and current cell level and proposed.
void distribute(const mapDistributePolyMesh &)
Update local numbering for mesh redistribution.
void updateMesh(const mapPolyMesh &)
Update local numbering for changed mesh.
static scalar & perturbTol()
Get the perturbation tolerance.
const globalMeshData & globalData() const
Return parallel info.
const vectorField & cellCentres() const
label nCells() const noexcept
Number of mesh cells.
int myProcNo() const noexcept
Return processor number.
Container with cells to refine. Refinement given as single direction.
@ OUTSIDE
A location outside the volume.
@ MIXED
A location that is partly inside and outside.
@ INSIDE
A location inside the volume.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
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< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
const dimensionSet dimless
Dimensionless.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
List< label > labelList
A List of labels.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
messageStream Info
Information stream (stdout output on master, null elsewhere)
vector point
Point is a vector.
treeDataPrimitivePatch< bPatch > treeDataBPatch
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
List< bool > boolList
A List of bools.
static constexpr const zero Zero
Global zero (0)
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
PrimitivePatch<::Foam::List< face >, const pointField > bPatch
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
List< face > faceList
A List of faces.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.