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)
129 void 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_);
446 void Foam::backgroundMeshDecomposition::printMeshData
478 Info<<
"Processor " << proci <<
" "
479 <<
"Number of cells = " << globalCells.localSize(proci)
503 bool Foam::backgroundMeshDecomposition::refineCell
507 scalar& weightEstimate
517 mesh_.cells()[celli].points
524 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_.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),
721 forAll(allBackgroundMeshBounds_, proci)
723 globalBackgroundBounds_.add(allBackgroundMeshBounds_[proci]);
731 /
"backgroundMeshDecomposition_proc_"
733 +
"_boundaryFaces.obj"
736 const faceList& faces = boundaryFacesPtr_().localFaces();
737 const List<point>&
points = boundaryFacesPtr_().localPoints();
739 Map<label> foamToObj(
points.size());
745 const face&
f = faces[i];
749 if (foamToObj.insert(
f[fPI], vertI))
760 fStr<<
' ' << foamToObj[
f[fPI]] + 1;
771 Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
775 const conformationSurfaces& geometryToConformTo,
776 const dictionary& coeffsDict,
777 const fileName& decompDictFile
781 geometryToConformTo_(geometryToConformTo),
787 "backgroundMeshDecomposition",
791 IOobject::AUTO_WRITE,
803 allBackgroundMeshBounds_(Pstream::nProcs()),
804 globalBackgroundBounds_(),
805 mergeDist_(1
e-6*mesh_.bounds().
mag()),
806 spanScale_(coeffsDict.
get<scalar>(
"spanScale")),
809 coeffsDict.getOrDefault<scalar>(
"minCellSizeLimit", 0)
811 minLevels_(coeffsDict.
get<label>(
"minLevels")),
812 volRes_(coeffsDict.
get<label>(
"sampleResolution")),
813 maxCellWeightCoeff_(coeffsDict.
get<scalar>(
"maxCellWeightCoeff"))
818 <<
"This cannot be used when not running in parallel."
822 const decompositionMethod& decomposer =
825 if (!decomposer.parallelAware())
828 <<
"You have selected decomposition method "
829 << decomposer.typeName
830 <<
" which is not parallel aware." <<
endl
834 Info<<
nl <<
"Building initial background mesh decomposition" <<
endl;
862 label nOccupiedCells = 0;
866 if (icellWeights[cI] > 1 - SMALL)
876 scalar cellWeightLimit =
max
879 *
sum(cellWeights).value()
886 Info<<
" cellWeightLimit " << cellWeightLimit <<
endl;
888 Pout<<
" sum(cellWeights) " <<
sum(cellWeights.primitiveField())
889 <<
" max(cellWeights) " <<
max(cellWeights.primitiveField())
897 if (icellWeights[cWI] > cellWeightLimit)
899 cellsToRefine.insert(cWI);
903 if (
returnReduce(cellsToRefine.size(), sumOp<label>()) == 0)
911 meshCutter_.consistentRefinement
918 if (
debug && !cellsToRefine.empty())
920 Pout<<
" cellWeights too large in " << cellsToRefine.size()
924 forAll(newCellsToRefine, nCTRI)
926 label celli = newCellsToRefine[nCTRI];
928 icellWeights[celli] /= 8.0;
932 polyTopoChange meshMod(mesh_);
935 meshCutter_.setRefinement(newCellsToRefine, meshMod);
938 autoPtr<mapPolyMesh> map = meshMod.changeMesh
948 mesh_.updateMesh(map());
951 meshCutter_.updateMesh(map());
953 Info<<
" Background mesh refined from "
955 <<
" to " << mesh_.globalData().nTotalCells()
956 <<
" cells." <<
endl;
969 printMeshData(mesh_);
971 Pout<<
" Pre distribute sum(cellWeights) "
973 <<
" max(cellWeights) "
978 decompositionMethod& decomposer =
981 labelList newDecomp = decomposer.decompose
988 Info<<
" Redistributing background mesh cells" <<
endl;
990 fvMeshDistribute distributor(mesh_);
992 autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute(newDecomp);
994 meshCutter_.distribute(mapDist());
998 printMeshData(mesh_);
1000 Pout<<
" Post distribute sum(cellWeights) "
1001 <<
sum(icellWeights)
1002 <<
" max(cellWeights) "
1003 <<
max(icellWeights)
1009 cellWeights.write();
1012 buildPatchAndTree();
1031 const List<point>& pts
1034 boolList posProc(pts.size(),
true);
1038 posProc[pI] = positionOnThisProcessor(pts[pI]);
1047 const treeBoundBox& box
1051 return !bFTreePtr_().findBox(box).empty();
1057 const point& centre,
1058 const scalar radiusSqr
1063 return bFTreePtr_().findNearest(centre, radiusSqr).hit();
1073 return bFTreePtr_().findLine(start,
end);
1083 return bFTreePtr_().findLineAny(start,
end);
1089 const List<point>& pts
1092 DynamicList<label> toCandidateProc;
1093 DynamicList<point> testPoints;
1097 label nTotalCandidates = 0;
1101 const point& pt = pts[pI];
1103 label nCandidates = 0;
1105 forAll(allBackgroundMeshBounds_, proci)
1109 if (allBackgroundMeshBounds_[proci].overlaps(pt,
sqr(SMALL*100)))
1111 toCandidateProc.append(proci);
1112 testPoints.append(pt);
1118 ptBlockStart[pI] = nTotalCandidates;
1119 ptBlockSize[pI] = nCandidates;
1121 nTotalCandidates += nCandidates;
1125 label preDistributionToCandidateProcSize = toCandidateProc.size();
1127 autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1129 map().distribute(testPoints);
1131 List<scalar> distanceSqrToCandidate(testPoints.size(),
sqr(GREAT));
1144 distanceSqrToCandidate[tPI] =
magSqr
1146 testPoints[tPI] - info.hitPoint()
1151 map().reverseDistribute
1153 preDistributionToCandidateProcSize,
1154 distanceSqrToCandidate
1157 labelList ptNearestProc(pts.size(), -1);
1163 SubList<scalar> ptNearestProcResults
1165 distanceSqrToCandidate,
1170 scalar nearestProcDistSqr = GREAT;
1172 forAll(ptNearestProcResults, pPRI)
1174 if (ptNearestProcResults[pPRI] < nearestProcDistSqr)
1176 nearestProcDistSqr = ptNearestProcResults[pPRI];
1178 ptNearestProc[pI] = toCandidateProc[ptBlockStart[pI] + pPRI];
1184 Pout<< pts[pI] <<
" nearestProcDistSqr " << nearestProcDistSqr
1185 <<
" ptNearestProc[pI] " << ptNearestProc[pI] <<
endl;
1188 if (ptNearestProc[pI] < 0)
1191 <<
"The position " << pts[pI]
1192 <<
" did not find a nearest point on the background mesh."
1197 return ptNearestProc;
1205 const List<point>& starts,
1206 const List<point>& ends,
1207 bool includeOwnProcessor
1210 DynamicList<label> toCandidateProc;
1211 DynamicList<point> testStarts;
1212 DynamicList<point> testEnds;
1213 labelList segmentBlockStart(starts.size(), -1);
1214 labelList segmentBlockSize(starts.size(), -1);
1216 label nTotalCandidates = 0;
1220 const point&
s = starts[sI];
1221 const point&
e = ends[sI];
1226 label nCandidates = 0;
1228 forAll(allBackgroundMeshBounds_, proci)
1235 && allBackgroundMeshBounds_[proci].intersects(
s,
e,
p)
1238 toCandidateProc.append(proci);
1239 testStarts.append(
s);
1246 segmentBlockStart[sI] = nTotalCandidates;
1247 segmentBlockSize[sI] = nCandidates;
1249 nTotalCandidates += nCandidates;
1253 label preDistributionToCandidateProcSize = toCandidateProc.size();
1255 autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1257 map().distribute(testStarts);
1258 map().distribute(testEnds);
1260 List<pointIndexHit> segmentIntersectsCandidate(testStarts.size());
1265 const point&
s = testStarts[sI];
1266 const point&
e = testEnds[sI];
1269 segmentIntersectsCandidate[sI] = bFTreePtr_().findLine(
s,
e);
1272 map().reverseDistribute
1274 preDistributionToCandidateProcSize,
1275 segmentIntersectsCandidate
1278 List<List<pointIndexHit>> segmentHitProcs(starts.size());
1281 DynamicList<pointIndexHit> tmpProcHits;
1285 tmpProcHits.clear();
1289 SubList<pointIndexHit> segmentProcResults
1291 segmentIntersectsCandidate,
1292 segmentBlockSize[sI],
1293 segmentBlockStart[sI]
1296 forAll(segmentProcResults, sPRI)
1298 if (segmentProcResults[sPRI].hit())
1300 tmpProcHits.append(segmentProcResults[sPRI]);
1302 tmpProcHits.last().setIndex
1304 toCandidateProc[segmentBlockStart[sI] + sPRI]
1309 segmentHitProcs[sI] = tmpProcHits;
1312 return segmentHitProcs;
1318 const point& centre,
1319 const scalar& radiusSqr
1322 forAll(allBackgroundMeshBounds_, proci)
1324 if (bFTreePtr_().findNearest(centre, radiusSqr).hit())
1336 const point& centre,
1337 const scalar radiusSqr
1342 forAll(allBackgroundMeshBounds_, proci)
1348 && allBackgroundMeshBounds_[proci].overlaps(centre, radiusSqr)
1354 toProc.append(proci);