45 void Foam::pairPatchAgglomeration::compactLevels(
const label nCreatedLevels)
53 bool Foam::pairPatchAgglomeration::continueAgglomerating
56 const label nLocalOld,
57 const label nMarkedEdges
66 label nGlobalOld =
returnReduce(nLocalOld, sumOp<label>());
67 label nGlobalMarked =
returnReduce(nMarkedEdges, sumOp<label>());
72 || nGlobal > nGlobalFacesInCoarsestLevel_
74 && (nGlobal != nGlobalOld || nGlobalMarked > 0);
78 void Foam::pairPatchAgglomeration::setLevel0EdgeWeights()
80 const bPatch& coarsePatch = patchLevels_[0];
81 const auto& coarseEdges = coarsePatch.edges();
87 for (label i = 0; i < coarsePatch.nInternalEdges(); i++)
89 scalar edgeLength = coarseEdges[i].mag(coarsePatch.localPoints());
91 const labelList& eFaces = coarsePatch.edgeFaces()[i];
93 if (eFaces.size() == 2)
96 coarsePatch.faceNormals()[eFaces[0]]
97 & coarsePatch.faceNormals()[eFaces[1]];
99 const edge edgeCommon = edge(eFaces[0], eFaces[1]);
101 if (facePairWeight_.found(edgeCommon))
103 facePairWeight_[edgeCommon] += edgeLength;
107 facePairWeight_.insert(edgeCommon, edgeLength);
112 facePairWeight_[edgeCommon] = -1.0;
120 for (label
k = j+1;
k<eFaces.size();
k++)
122 facePairWeight_.insert
124 edge(eFaces[j], eFaces[
k]),
136 <<
" nEdges:" << coarsePatch.nEdges() <<
" of which:" <<
nl 137 <<
" boundary:" << coarsePatch.nBoundaryEdges() <<
nl 138 <<
" non-manifold:" << nNonManif <<
nl 139 <<
" feature (angle < " << featureAngle_ <<
"):" << nFeat <<
nl 145 void Foam::pairPatchAgglomeration::setEdgeWeights
147 const label fineLevelIndex
150 const bPatch& coarsePatch = patchLevels_[fineLevelIndex];
151 const auto& coarseEdges = coarsePatch.edges();
152 const labelList& fineToCoarse = restrictAddressing_[fineLevelIndex];
154 edgeHashSet fineFeaturedFaces(coarsePatch.nEdges()/10);
161 const edge
e = iter.key();
162 const edge edgeFeatured
167 fineFeaturedFaces.insert(edgeFeatured);
177 facePairWeight_.clear();
178 facePairWeight_.reserve(coarsePatch.nEdges());
180 for (label i = 0; i < coarsePatch.nInternalEdges(); i++)
182 scalar edgeLength = coarseEdges[i].mag(coarsePatch.localPoints());
184 const labelList& eFaces = coarsePatch.edgeFaces()[i];
186 if (eFaces.size() == 2)
188 const edge edgeCommon(eFaces[0], eFaces[1]);
191 if (fineFeaturedFaces.found(edgeCommon))
193 auto w = facePairWeight_.find(edgeCommon);
196 facePairWeight_.insert(edgeCommon, -1.0);
199 else if (w() != -1.0)
208 auto w = facePairWeight_.find(edgeCommon);
218 facePairWeight_.insert(edgeCommon, edgeLength);
227 for (label
k = j+1;
k<eFaces.size();
k++)
229 facePairWeight_.insert
231 edge(eFaces[j], eFaces[
k]),
242 Pout<<
"Level:" << fineLevelIndex
243 <<
" nEdges:" << coarsePatch.nEdges() <<
" of which:" <<
nl 244 <<
" boundary:" << coarsePatch.nBoundaryEdges() <<
nl 245 <<
" non-manifold:" << nNonManif <<
nl 246 <<
" feature (angle < " << featureAngle_ <<
"):" << nFeat <<
nl 252 bool Foam::pairPatchAgglomeration::isSingleEdgeLoop
262 SubList<label>(allFaces, faceIDs.size()) = faceIDs;
263 allFaces.
last() = facei;
268 IndirectList<face>(
patch, allFaces),
272 return (upp.edgeLoops().size() == 1);
276 Foam::label Foam::pairPatchAgglomeration::maxValidNeighbour
278 const bool addToCluster,
288 const auto& fFaces =
patch.faceFaces()[facei];
290 label matchFaceNeibNo = -1;
291 scalar maxFaceWeight = -0.5;
297 for (
const label faceNeig : fFaces)
299 const label coarsei = fineToCoarse[faceNeig];
303 const edge edgeCommon = edge(facei, faceNeig);
304 const auto& weight = facePairWeight_[edgeCommon];
307 (weight > maxFaceWeight)
311 maxFaceWeight = weight;
312 matchFaceNeibNo = faceNeig;
320 for (
const label faceNeig : fFaces)
322 const label coarsei = fineToCoarse[faceNeig];
326 const edge edgeCommon = edge(facei, faceNeig);
327 const auto& weight = facePairWeight_[edgeCommon];
328 if (weight > maxFaceWeight)
330 maxFaceWeight = weight;
331 matchFaceNeibNo = faceNeig;
337 return matchFaceNeibNo;
343 Foam::pairPatchAgglomeration::pairPatchAgglomeration
355 nFacesInCoarsestLevel_
359 nGlobalFacesInCoarsestLevel_(
labelMax),
365 controlDict.getOrDefault<scalar>(
"featureAngle", 0)
368 restrictAddressing_(maxLevels_),
369 restrictTopBottomAddressing_(
identity(faces.size())),
370 patchLevels_(maxLevels_),
371 facePairWeight_(faces.size())
380 setLevel0EdgeWeights();
384 Foam::pairPatchAgglomeration::pairPatchAgglomeration
388 const label mergeLevels,
389 const label maxLevels,
390 const label nFacesInCoarsestLevel,
391 const label nGlobalFacesInCoarsestLevel,
392 const scalar featureAngle
395 mergeLevels_(mergeLevels),
396 maxLevels_(maxLevels),
397 nFacesInCoarsestLevel_(nFacesInCoarsestLevel),
398 nGlobalFacesInCoarsestLevel_(nGlobalFacesInCoarsestLevel),
399 featureAngle_(featureAngle),
401 restrictAddressing_(maxLevels_),
402 restrictTopBottomAddressing_(
identity(faces.size())),
403 patchLevels_(maxLevels_),
404 facePairWeight_(faces.size())
413 setLevel0EdgeWeights();
431 return patchLevels_[i];
435 void Foam::pairPatchAgglomeration::mapBaseToTopAgglom
437 const label fineLevelIndex
440 const labelList& fineToCoarse = restrictAddressing_[fineLevelIndex];
441 forAll(restrictTopBottomAddressing_, i)
443 restrictTopBottomAddressing_[i] =
444 fineToCoarse[restrictTopBottomAddressing_[i]];
449 bool Foam::pairPatchAgglomeration::agglomeratePatch
453 const label fineLevelIndex,
457 if (
min(fineToCoarse) == -1)
463 if (fineToCoarse.size() == 0)
468 if (fineToCoarse.size() !=
patch.size())
471 <<
"restrict map does not correspond to fine level. " <<
endl 472 <<
" Sizes: restrictMap: " << fineToCoarse.size()
473 <<
" nEqns: " <<
patch.size()
477 const label nCoarseI =
max(fineToCoarse) + 1;
478 List<face> patchFaces(nCoarseI);
487 for (label coarseI = 0; coarseI < nCoarseI; coarseI++)
489 const labelList& fineFaces = coarseToFine[coarseI];
494 IndirectList<face>(
patch, fineFaces),
498 if (upp.edgeLoops().size() != 1)
520 if (fineFaces.size() >= 2)
524 for (label
k = j+1;
k<fineFaces.size();
k++)
526 const edge
e(fineFaces[j], fineFaces[
k]);
528 auto w = facePairWeight_.find(
e);
531 facePairWeight_.insert(
e, -1.0);
534 else if (w() != -1.0)
548 patchFaces[coarseI] = face
563 SubList<face>(patchFaces, nCoarseI, 0),
574 label nPairLevels = 0;
575 label nCreatedLevels = 1;
577 label nCoarseFaces = 0;
578 label nCoarseFacesOld = 0;
579 label nMarkedEdges = 0;
581 while (nCreatedLevels < maxLevels_)
583 const bPatch&
patch = patchLevels_[nCreatedLevels - 1];
588 bool createdLevel =
false;
589 while (!createdLevel)
593 tfinalAgglom = agglomerateOneLevel(nCoarseFaces,
patch);
595 if (nCoarseFaces == 0)
604 createdLevel = agglomeratePatch
618 const auto& agglomPatch = patchLevels_[nCreatedLevels];
620 Pout<<
"Writing new patch at level:" << nCreatedLevels
622 os.
write(agglomPatch, agglomPatch.points(),
true);
625 restrictAddressing_.set(nCreatedLevels, tfinalAgglom);
627 mapBaseToTopAgglom(nCreatedLevels);
629 setEdgeWeights(nCreatedLevels);
631 if (nPairLevels % mergeLevels_)
633 combineLevels(nCreatedLevels);
642 nFaces_[nCreatedLevels] = nCoarseFaces;
647 if (!continueAgglomerating(nCoarseFaces, nCoarseFacesOld, nMarkedEdges))
652 nCoarseFacesOld = nCoarseFaces;
655 compactLevels(nCreatedLevels);
665 const label nFineFaces =
patch.size();
668 auto& coarseCellMap = tcoarseCellMap.ref();
677 if (coarseCellMap[facei] < 0)
679 const label matchFaceNeibNo = maxValidNeighbour
688 if (matchFaceNeibNo >= 0)
691 coarseCellMap[facei] = nCoarseFaces;
692 coarseCellMap[matchFaceNeibNo] = nCoarseFaces;
701 const label clusterMatchFaceNo = maxValidNeighbour
710 if (clusterMatchFaceNo >= 0)
713 const label coarsei = coarseCellMap[clusterMatchFaceNo];
714 coarseCellMap[facei] = coarsei;
720 coarseCellMap[facei] = nCoarseFaces;
731 for (label facei=0; facei<nFineFaces; facei++)
733 if (coarseCellMap[facei] < 0)
737 <<
" is not part of a cluster" 742 return tcoarseCellMap;
746 void Foam::pairPatchAgglomeration::combineLevels(
const label curLevel)
748 label prevLevel = curLevel - 1;
751 nFaces_[prevLevel] = nFaces_[curLevel];
756 const labelList& curResAddr = restrictAddressing_[curLevel];
757 labelList& prevResAddr = restrictAddressing_[prevLevel];
761 prevResAddr[i] = curResAddr[prevResAddr[i]];
765 restrictAddressing_.set(curLevel,
nullptr);
767 patchLevels_.set(prevLevel, patchLevels_.set(curLevel,
nullptr));
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
List of values generated by applying the access operation to each list item.
PrimitivePatch< List< face >, const pointField > bPatch
void size(const label n)
Older name for setAddressableSize.
PrimitivePatch<::Foam::List< face >, const pointField > bPatch
errorManipArg< error, int > exit(error &err, const int errNo=1)
virtual Ostream & write(const char c) override
Write character.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
virtual const fileName & name() const override
Read/write access to the name of the stream.
PtrList< labelField > restrictAddressing_
Cell restriction addressing array.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
IntListType renumber(const labelUList &oldToNew, const IntListType &input)
Renumber the values within a list.
Unit conversion functions.
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
HashSet< edge, Hash< edge > > edgeHashSet
A HashSet with edge for its key. Hashing (and ==) on an edge is symmetric.
label k
Boltzmann constant.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field...
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.
PtrList< bPatch > patchLevels_
Hierarchy of patch addressing.
#define forAll(list, i)
Loop across all elements in list.
virtual ~pairPatchAgglomeration()
A list of faces which address into the list of points.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
dimensionedScalar cos(const dimensionedScalar &ds)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
labelListList invertOneToMany(const label len, const labelUList &map)
Invert one-to-many map. Unmapped elements will be size 0.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
const bPatch & patchLevel(const label leveli) const
Return primitivePatch of given level.
void agglomerate()
Agglomerate patch.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
runTime controlDict().readEntry("adjustTimeStep"
The central control dictionary, the contents of which are either taken directly from the FOAM_CONTROL...
errorManip< error > abort(error &err)
An OFstream that keeps track of vertices and provides convenience output methods for OBJ files...
int debug
Static debugging option.
OBJstream os(runTime.globalPath()/outputName)
defineTypeNameAndDebug(combustionModel, 0)
T & last()
Access last element of the list, position [size()-1].
const std::string patch
OpenFOAM patch number as a std::string.
labelList nFaces_
The number of faces in each level.
List< label > labelList
A List of labels.
A class for managing temporary objects.
bool returnReduceOr(const bool value, const label comm=UPstream::worldComm)
Perform logical (or) MPI Allreduce on a copy. Uses UPstream::reduceOr.
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
forAllConstIters(mixture.phases(), phase)