Compare commits
4 Commits
multiNodeD
...
hpathRenum
| Author | SHA1 | Date | |
|---|---|---|---|
| 44c76c5806 | |||
| aaf0245795 | |||
| 69df0ad468 | |||
| fdfd2d15a9 |
@ -9,7 +9,7 @@ EXE_INC = \
|
||||
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
|
||||
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
|
||||
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
|
||||
-I$(LIB_SRC)/regionFaModels/lnInclude
|
||||
-I$(LIB_SRC)/regionFaModels\lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-lfiniteVolume \
|
||||
|
||||
@ -7,7 +7,7 @@ EXE_INC = \
|
||||
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
|
||||
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
|
||||
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
|
||||
-I$(LIB_SRC)/regionFaModels/lnInclude
|
||||
-I$(LIB_SRC)/regionFaModels\lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-lfiniteVolume \
|
||||
|
||||
@ -10,7 +10,7 @@ EXE_INC = \
|
||||
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
|
||||
-I$(LIB_SRC)/dynamicMesh/lnInclude \
|
||||
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
|
||||
-I$(LIB_SRC)/regionFaModels/lnInclude
|
||||
-I$(LIB_SRC)/regionFaModels\lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-lfiniteVolume \
|
||||
|
||||
@ -7,7 +7,7 @@ EXE_INC = \
|
||||
-I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
|
||||
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
|
||||
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
|
||||
-I$(LIB_SRC)/regionFaModels/lnInclude
|
||||
-I$(LIB_SRC)/regionFaModels\lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-lfiniteVolume \
|
||||
|
||||
@ -19,7 +19,7 @@ EXE_INC = \
|
||||
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
|
||||
-I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
|
||||
-I$(LIB_SRC)/regionModels/regionModel/lnInclude \
|
||||
-I$(LIB_SRC)/regionFaModels/lnInclude
|
||||
-I$(LIB_SRC)/regionFaModels\lnInclude
|
||||
|
||||
|
||||
EXE_LIBS = \
|
||||
|
||||
@ -8,7 +8,7 @@ EXE_INC = \
|
||||
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
|
||||
-I$(LIB_SRC)/dynamicMesh/lnInclude \
|
||||
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
|
||||
-I$(LIB_SRC)/regionFaModels/lnInclude
|
||||
-I$(LIB_SRC)/regionFaModels\lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-lfiniteVolume \
|
||||
|
||||
@ -30,8 +30,8 @@ Description
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef Foam_CGALTriangulation3DKernel_H
|
||||
#define Foam_CGALTriangulation3DKernel_H
|
||||
#ifndef CGALTriangulation3DKernel_H
|
||||
#define CGALTriangulation3DKernel_H
|
||||
|
||||
// Silence boost bind deprecation warnings (before CGAL-5.2.1)
|
||||
#include "CGAL/version.h"
|
||||
@ -54,19 +54,9 @@ Description
|
||||
// #include "CGAL/Robust_circumcenter_traits_3.h"
|
||||
// typedef CGAL::Robust_circumcenter_traits_3<baseK> K;
|
||||
|
||||
#if defined(CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1050500000)
|
||||
// Prior to CGAL-5.5
|
||||
|
||||
#include "CGAL/Robust_circumcenter_filtered_traits_3.h"
|
||||
typedef CGAL::Robust_circumcenter_filtered_traits_3<baseK> K;
|
||||
|
||||
#else
|
||||
|
||||
#include "CGAL/Robust_weighted_circumcenter_filtered_traits_3.h"
|
||||
typedef CGAL::Robust_weighted_circumcenter_filtered_traits_3<baseK> K;
|
||||
|
||||
#endif
|
||||
|
||||
#else
|
||||
|
||||
// Very robust but expensive kernel
|
||||
|
||||
@ -4,6 +4,7 @@ EXE_INC = \
|
||||
-Wno-old-style-cast \
|
||||
$(COMP_FLAGS) \
|
||||
${CGAL_INC} \
|
||||
-DCGAL_HEADER_ONLY \
|
||||
-I$(LIB_SRC)/finiteVolume/lnInclude \
|
||||
-I$(LIB_SRC)/surfMesh/lnInclude \
|
||||
-I$(LIB_SRC)/meshTools/lnInclude \
|
||||
@ -11,6 +12,7 @@ EXE_INC = \
|
||||
|
||||
|
||||
EXE_LIBS = \
|
||||
/* ${CGAL_LIBS} */ \
|
||||
-lfiniteVolume \
|
||||
-lsurfMesh \
|
||||
-lmeshTools \
|
||||
|
||||
@ -79,8 +79,8 @@ bool Foam::laminarModels::generalizedNewtonianViscosityModels::Casson::read
|
||||
|
||||
coeffs.readEntry("m", m_);
|
||||
coeffs.readEntry("tau0", tau0_);
|
||||
coeffs.readEntry("nuMin", nuMin_);
|
||||
coeffs.readEntry("nuMax", nuMax_);
|
||||
coeffs.readEntry("nuMin_", nuMin_);
|
||||
coeffs.readEntry("nuMax_", nuMax_);
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
@ -389,7 +389,7 @@ bool Foam::functionObjects::fvExpressionField::read(const dictionary& dict)
|
||||
}
|
||||
|
||||
autowrite_ = dict.getOrDefault("autowrite", false);
|
||||
store_ = dict.getOrDefault("store", true);
|
||||
store_ = dict.getOrDefault("autowrite", true);
|
||||
|
||||
// "dimensions" is optional
|
||||
dimensions_.clear();
|
||||
|
||||
@ -103,7 +103,7 @@ bool Foam::functionObjects::ensightWrite::readSelection(const dictionary& dict)
|
||||
selectFields_.uniq();
|
||||
|
||||
blockFields_.clear();
|
||||
dict.readIfPresent("excludeFields", blockFields_);
|
||||
dict.readIfPresent("blockFields", blockFields_);
|
||||
blockFields_.uniq();
|
||||
|
||||
// Actions to define selection
|
||||
|
||||
@ -5,8 +5,8 @@
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2007-2023 PCOpt/NTUA
|
||||
Copyright (C) 2013-2023 FOSS GP
|
||||
Copyright (C) 2007-2021 PCOpt/NTUA
|
||||
Copyright (C) 2013-2021 FOSS GP
|
||||
Copyright (C) 2019 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -86,9 +86,10 @@ void ATCUaGradU::addATC(fvVectorMatrix& UaEqn)
|
||||
if (extraConvection_ > 0)
|
||||
{
|
||||
// Implicit part added to increase diagonal dominance
|
||||
UaEqn += ATClimiter_*extraConvection_*fvm::div(-phi, Ua);
|
||||
// Note: Maybe this needs to be multiplied with the ATClimiter ??
|
||||
UaEqn += extraConvection_*fvm::div(-phi, Ua);
|
||||
|
||||
// Correct rhs due to implicitly augmenting the adjoint convection
|
||||
// correct rhs due to implicitly augmenting the adjoint convection
|
||||
ATC_ += extraConvection_*(fvc::grad(UaForATC(), "gradUaATC")().T() & U);
|
||||
}
|
||||
|
||||
|
||||
@ -5,8 +5,8 @@
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2007-2023 PCOpt/NTUA
|
||||
Copyright (C) 2013-2023 FOSS GP
|
||||
Copyright (C) 2007-2021 PCOpt/NTUA
|
||||
Copyright (C) 2013-2021 FOSS GP
|
||||
Copyright (C) 2019 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -90,9 +90,10 @@ void ATCstandard::addATC(fvVectorMatrix& UaEqn)
|
||||
if (extraConvection_ > 0)
|
||||
{
|
||||
// Implicit part added to increase diagonal dominance
|
||||
UaEqn += ATClimiter_*extraConvection_*fvm::div(-phi, Ua);
|
||||
// Note: Maybe this needs to be multiplied with the ATClimiter ??
|
||||
UaEqn += extraConvection_*fvm::div(-phi, Ua);
|
||||
|
||||
// Correct rhs due to implicitly augmenting the adjoint convection
|
||||
// correct rhs due to implicitly augmenting the adjoint convection
|
||||
ATC_ += extraConvection_*(fvc::grad(Ua, "gradUaATC")().T() & U);
|
||||
}
|
||||
|
||||
|
||||
@ -4,7 +4,6 @@ simpleGeomDecomp/simpleGeomDecomp.C
|
||||
hierarchGeomDecomp/hierarchGeomDecomp.C
|
||||
manualDecomp/manualDecomp.C
|
||||
multiLevelDecomp/multiLevelDecomp.C
|
||||
multiNodeDecomp/multiNodeDecomp.C
|
||||
metisLikeDecomp/metisLikeDecomp.C
|
||||
structuredDecomp/structuredDecomp.C
|
||||
randomDecomp/randomDecomp.C
|
||||
|
||||
@ -1,58 +0,0 @@
|
||||
# New Multi-Level Decomposition
|
||||
The multi-node decomposition is an extension of the existing multi-level decomposition. It supports the syntax of the current multi-level decomposition, but allows to change the decomposition tree as you wish. For example, you may split into unbalanced nodes, set the weights of some nodes to be bigger than others, or perhaps use a different decomposition method for some nodes.
|
||||
You may set up the decomposition in two ways:
|
||||
1. Using a domains list and a default method:
|
||||
```
|
||||
numberOfSubdomains 8;
|
||||
multiNodeCoeffs {
|
||||
domains (2 4);
|
||||
method metis;
|
||||
}
|
||||
```
|
||||
2. Using a dictionary for each level:
|
||||
```
|
||||
numberOfSubdomains 8;
|
||||
multiLevelCoeffs {
|
||||
nodes {
|
||||
numberOfSubdomains 2;
|
||||
method metis;
|
||||
}
|
||||
cores {
|
||||
numberOfSubdomains 4;
|
||||
method scotch;
|
||||
}
|
||||
}
|
||||
```
|
||||
Note that if the total number of subdomains does not match the product of the number of subdomains at each level, but a default method is provided, a new level will be inferred in order to match the total number of subdomains.
|
||||
|
||||
This creates a "decomposition tree" - for example, the dictionaries above create a tree, where the root has two children, and each child has four children (who are the leaves of the tree). Every leaf in the tree is a subdomain in the final decomposition.
|
||||
After setting up the decomposition, we may edit specific nodes or ranges of nodes. For example, suppose we want to split into two nodes, the first one having four subdomains and the second having eight subdomains. We can use the above dictionaries, and then use:
|
||||
```
|
||||
domains[1] (8);
|
||||
```
|
||||
The squared brackets indicate which nodes in the tree should we edit - We want the second child of the root (the indexing starts from zero). If we wanted to change the first two children of the third child of the root, we would write:
|
||||
```
|
||||
domains[2][0-1] (8);
|
||||
```
|
||||
|
||||
Note that the total number of subdomains must match the number of subdomains declared after all modifications. In addition, note that the decomposition into two nodes will be done as if they were of the same size, hence the first four subdomains will be bigger than the other eight. In order to fix this, we may:
|
||||
1. Change the weight of the second node into twice the weight:
|
||||
```
|
||||
weight[1] 2;
|
||||
```
|
||||
2. Set the weights initialization into relative - this will cause the weights of the children to first be computed by the amount of leaves in their subtree. Note that this updates the whole subtree initialization, but using the `weight` parameter, we can override this initialization.
|
||||
```
|
||||
weightsInitialization[1] relative;
|
||||
```
|
||||
|
||||
|
||||
We may also set a special method dictionary that decomposes differently for some nodes:
|
||||
```
|
||||
method[2-4] {
|
||||
numberOfSubdomains 4;
|
||||
method metis;
|
||||
coeffs {
|
||||
...
|
||||
}
|
||||
}
|
||||
```
|
||||
@ -1,788 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2017 OpenFOAM Foundation
|
||||
Copyright (C) 2017-2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
||||
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
||||
for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "multiNodeDecomp.H"
|
||||
#include "addToRunTimeSelectionTable.H"
|
||||
#include "IFstream.H"
|
||||
#include "globalIndex.H"
|
||||
#include "mapDistribute.H"
|
||||
#include "DynamicList.H"
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
defineTypeNameAndDebug(multiNodeDecomp, 0);
|
||||
addToRunTimeSelectionTable
|
||||
(
|
||||
decompositionMethod,
|
||||
multiNodeDecomp,
|
||||
dictionary
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
namespace Foam {
|
||||
void multiNodeDecomp::initializeMetadata(const dictionary& coeffsDict) {
|
||||
word defaultMethod;
|
||||
dictionary defaultMethodDict;
|
||||
if(coeffsDict.readIfPresent("method", defaultMethod, keyType::LITERAL)) {
|
||||
defaultMethodDict.add("method",defaultMethod);
|
||||
const dictionary& subMethodCoeffsDict
|
||||
(
|
||||
findCoeffsDict
|
||||
(
|
||||
coeffsDict,
|
||||
defaultMethod + "Coeffs",
|
||||
selectionType::NULL_DICT
|
||||
)
|
||||
);
|
||||
if(subMethodCoeffsDict.size())
|
||||
defaultMethodDict.add(subMethodCoeffsDict.dictName(), subMethodCoeffsDict);
|
||||
}
|
||||
labelList domains;
|
||||
|
||||
label nTotal = 0;
|
||||
label nLevels = 0;
|
||||
|
||||
//Check if any meta argument is changed using the new syntax.
|
||||
//If they are, we cannot infer an additional level of decomposition,
|
||||
//as it may interfere with the indices.
|
||||
List<string> domainChanges = metaParser::getEntries(coeffsDict, "domains");
|
||||
List<string> methodChanges = metaParser::getEntries(coeffsDict, "method");
|
||||
List<string> weightChanges = metaParser::getEntries(coeffsDict, "weight");
|
||||
//We can parse weightMode without brackets too
|
||||
List<string> weightModeChanges = metaParser::getEntries(coeffsDict, "weightsInitialization", true);
|
||||
|
||||
bool bChangesDomains = !domainChanges.empty();
|
||||
bool bChangesArguments = bChangesDomains
|
||||
|| (!methodChanges.empty())
|
||||
|| (!weightChanges.empty())
|
||||
|| (!weightModeChanges.empty());
|
||||
|
||||
bool bMetadataInitialized = false;
|
||||
|
||||
// Found (non-recursive, no patterns) "method" and "domains" ?
|
||||
// Allow as quick short-cut entry
|
||||
if
|
||||
(
|
||||
// non-recursive, no patterns
|
||||
coeffsDict.readIfPresent("method", defaultMethod, keyType::LITERAL)
|
||||
// non-recursive, no patterns
|
||||
&& coeffsDict.readIfPresent("domains", domains, keyType::LITERAL)
|
||||
)
|
||||
{
|
||||
// Short-cut version specified by method, domains only
|
||||
|
||||
nTotal = (domains.empty() ? 0 : 1);
|
||||
for (const label n : domains)
|
||||
{
|
||||
nTotal *= n;
|
||||
++nLevels;
|
||||
}
|
||||
|
||||
//Update domains here
|
||||
if(nTotal != 0 && bChangesDomains) {
|
||||
rootMetadata_.initialize(
|
||||
domains,
|
||||
&defaultMethodDict
|
||||
);
|
||||
bMetadataInitialized = true;
|
||||
for(string key : domainChanges)
|
||||
rootMetadata_.updateDomains( key,
|
||||
coeffsDict.get<labelList>(key, keyType::LITERAL));
|
||||
|
||||
nTotal = rootMetadata_.getSize();
|
||||
}
|
||||
|
||||
if (nTotal == 1)
|
||||
{
|
||||
// Emit Warning
|
||||
nTotal = nDomains();
|
||||
nLevels = 1;
|
||||
domains.setSize(1);
|
||||
domains[0] = nTotal;
|
||||
}
|
||||
//If bChangesDomains is true, we do not want to add another dimension as this
|
||||
//may affect the user's assignments of domains/weights/methods later on.
|
||||
else if (nTotal > 0 && nTotal < nDomains() && !(nDomains() % nTotal) && !bChangesArguments)
|
||||
{
|
||||
// nTotal < nDomains, but with an integral factor,
|
||||
// which we insert as level 0
|
||||
++nLevels;
|
||||
|
||||
labelList old(std::move(domains));
|
||||
|
||||
domains.setSize(old.size()+1);
|
||||
|
||||
domains[0] = nDomains() / nTotal;
|
||||
forAll(old, i)
|
||||
{
|
||||
domains[i+1] = old[i];
|
||||
}
|
||||
nTotal *= domains[0];
|
||||
|
||||
|
||||
Info<<" inferred level 0 with " << domains[0]
|
||||
<< " domains" << nl << nl;
|
||||
}
|
||||
|
||||
if (!nLevels || nTotal != nDomains())
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Top level decomposition specifies " << nDomains()
|
||||
<< " domains which is not equal to the product of"
|
||||
<< " all sub domains " << nTotal
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
if(!bMetadataInitialized) {
|
||||
bMetadataInitialized = true;
|
||||
rootMetadata_.initialize(
|
||||
domains,
|
||||
&defaultMethodDict
|
||||
);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Specified by full dictionaries
|
||||
|
||||
// Create editable methods dictionaries
|
||||
// - Only consider sub-dictionaries with a "numberOfSubdomains" entry
|
||||
// This automatically filters out any coeffs dictionaries
|
||||
|
||||
label nTotal = 1;
|
||||
List<const dictionary*> methods;
|
||||
for (const entry& dEntry : coeffsDict)
|
||||
{
|
||||
word methodName;
|
||||
|
||||
if
|
||||
(
|
||||
dEntry.isDict()
|
||||
// non-recursive, no patterns
|
||||
&& dEntry.dict().found("numberOfSubdomains", keyType::LITERAL)
|
||||
)
|
||||
{
|
||||
domains.append(dEntry.dict().get<label>("numberOfSubdomains"));
|
||||
nTotal *= domains.last();
|
||||
// No method specified? can use a default method?
|
||||
|
||||
const bool addDefaultMethod
|
||||
(
|
||||
!(dEntry.dict().found("method", keyType::LITERAL))
|
||||
&& !defaultMethod.empty()
|
||||
);
|
||||
if(!(dEntry.dict().found("method",keyType::LITERAL)) && defaultMethod.empty()) {
|
||||
FatalErrorInFunction <<
|
||||
dEntry.keyword() <<
|
||||
" dictionary does not contain method, and no default method is specified."
|
||||
<< nl << exit(FatalError);
|
||||
}
|
||||
dictionary* levelDict = new dictionary(dEntry.dict());
|
||||
levelDict->remove("numberOfSubdomains");
|
||||
if(addDefaultMethod) levelDict->add("method", defaultMethod);
|
||||
methods.append(levelDict);
|
||||
}
|
||||
}
|
||||
if(domains.empty())
|
||||
nTotal = 0;
|
||||
|
||||
|
||||
rootMetadata_.initialize(domains, methods[0]);
|
||||
bMetadataInitialized = true;
|
||||
for(string key : domainChanges)
|
||||
rootMetadata_.updateDomains( key,
|
||||
coeffsDict.get<labelList>(key, keyType::LITERAL));
|
||||
|
||||
if(nTotal != nDomains()) {
|
||||
FatalErrorInFunction
|
||||
<< "Top level decomposition specifies " << nDomains()
|
||||
<< " domains which is not equal to the product of"
|
||||
<< " all sub domains " << nTotal << " manually defined by dictionaries. "
|
||||
<< exit(FatalError);
|
||||
}
|
||||
rootMetadata_.setLeveledDictionaries(methods);
|
||||
for(const dictionary* method : methods)
|
||||
delete method;
|
||||
}
|
||||
|
||||
|
||||
for(string key : methodChanges)
|
||||
rootMetadata_.updateMethod(key, coeffsDict.subDict(key, keyType::LITERAL));
|
||||
|
||||
for(string key : weightChanges)
|
||||
rootMetadata_.updateWeight(key, coeffsDict.get<label>(key, keyType::LITERAL));
|
||||
|
||||
for(string key : weightModeChanges) {
|
||||
word value = coeffsDict.get<word>(key, keyType::LITERAL);
|
||||
WeightsInitialization newValue = UNKNOWN;
|
||||
|
||||
if(value=="uniform")
|
||||
newValue = UNIFORM;
|
||||
else if(value == "relative")
|
||||
newValue = RELATIVE;
|
||||
else
|
||||
FatalErrorInFunction <<
|
||||
"unknown weights initialization (" << value << "). Must be one of: relative, uniform."
|
||||
<< nl << exit(FatalError);
|
||||
|
||||
rootMetadata_.updateWeightsInitialization(key, newValue);
|
||||
}
|
||||
|
||||
if(!rootMetadata_.isLeaf())
|
||||
rootMetadata_.constructMethods();
|
||||
}
|
||||
|
||||
|
||||
// Given a subset of cells determine the new global indices. The problem
|
||||
// is in the cells from neighbouring processors which need to be renumbered.
|
||||
void multiNodeDecomp::subsetGlobalCellCells
|
||||
(
|
||||
const label nDomains,
|
||||
const label domainI,
|
||||
const labelList& dist,
|
||||
|
||||
const labelListList& cellCells,
|
||||
const labelList& set,
|
||||
labelListList& subCellCells,
|
||||
labelList& cutConnections
|
||||
) const
|
||||
{
|
||||
// Determine new index for cells by inverting subset
|
||||
labelList oldToNew(invert(cellCells.size(), set));
|
||||
|
||||
globalIndex globalCells(cellCells.size());
|
||||
|
||||
// Subset locally the elements for which I have data
|
||||
subCellCells = UIndirectList<labelList>(cellCells, set);
|
||||
|
||||
// Get new indices for neighbouring processors
|
||||
List<Map<label>> compactMap;
|
||||
mapDistribute map(globalCells, subCellCells, compactMap);
|
||||
map.distribute(oldToNew);
|
||||
labelList allDist(dist);
|
||||
map.distribute(allDist);
|
||||
|
||||
// Now we have:
|
||||
// oldToNew : the locally-compact numbering of all our cellCells. -1 if
|
||||
// cellCell is not in set.
|
||||
// allDist : destination domain for all our cellCells
|
||||
// subCellCells : indexes into oldToNew and allDist
|
||||
|
||||
// Globally compact numbering for cells in set.
|
||||
globalIndex globalSubCells(set.size());
|
||||
|
||||
// Now subCellCells contains indices into oldToNew which are the
|
||||
// new locations of the neighbouring cells.
|
||||
|
||||
cutConnections.setSize(nDomains);
|
||||
cutConnections = 0;
|
||||
|
||||
forAll(subCellCells, subCelli)
|
||||
{
|
||||
labelList& cCells = subCellCells[subCelli];
|
||||
|
||||
// Keep the connections to valid mapped cells
|
||||
label newI = 0;
|
||||
forAll(cCells, i)
|
||||
{
|
||||
// Get locally-compact cell index of neighbouring cell
|
||||
const label nbrCelli = oldToNew[cCells[i]];
|
||||
if (nbrCelli == -1)
|
||||
{
|
||||
cutConnections[allDist[cCells[i]]]++;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Reconvert local cell index into global one
|
||||
|
||||
// Get original neighbour
|
||||
const label celli = set[subCelli];
|
||||
const label oldNbrCelli = cellCells[celli][i];
|
||||
// Get processor from original neighbour
|
||||
const label proci = globalCells.whichProcID(oldNbrCelli);
|
||||
// Convert into global compact numbering
|
||||
cCells[newI++] = globalSubCells.toGlobal(proci, nbrCelli);
|
||||
}
|
||||
}
|
||||
cCells.setSize(newI);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void multiNodeDecomp::decompose
|
||||
(
|
||||
const labelListList& pointPoints,
|
||||
const pointField& points,
|
||||
const scalarField& pointWeights,
|
||||
const labelUList& pointMap, // map back to original points
|
||||
const nodeMetadata& decomposeData,
|
||||
const label leafOffset,
|
||||
|
||||
labelList& finalDecomp
|
||||
) const
|
||||
{
|
||||
labelList dist
|
||||
(
|
||||
decomposeData.getMethod()->decompose
|
||||
(
|
||||
pointPoints,
|
||||
points,
|
||||
pointWeights
|
||||
)
|
||||
);
|
||||
|
||||
// Number of domains at the current level
|
||||
const label nCurrDomains = decomposeData.nDomains();
|
||||
|
||||
// Calculate the domain remapping.
|
||||
// The decompose() method delivers a distribution of [0..nDomains-1]
|
||||
// which we map to the final location according to the decomposition
|
||||
// leaf we are on.
|
||||
|
||||
labelList domainOffsets(nCurrDomains);
|
||||
domainOffsets[0] = leafOffset;
|
||||
for(label nDomain = 1; nDomain < nCurrDomains; ++nDomain) {
|
||||
domainOffsets[nDomain] = domainOffsets[nDomain-1] + decomposeData.getChild(nDomain-1)->getSize();
|
||||
}
|
||||
|
||||
// Extract processor+local index from point-point addressing
|
||||
forAll(pointMap, i)
|
||||
{
|
||||
finalDecomp[pointMap[i]] = domainOffsets[dist[i]];
|
||||
}
|
||||
|
||||
if (nCurrDomains > 0)
|
||||
{
|
||||
// Recurse
|
||||
|
||||
// Determine points per domain
|
||||
labelListList domainToPoints(invertOneToMany(nCurrDomains, dist));
|
||||
|
||||
for (label domainI = 0; domainI < nCurrDomains; ++domainI)
|
||||
{
|
||||
if(decomposeData.getChild(domainI)->isLeaf()) continue;
|
||||
// Extract elements for current domain
|
||||
const labelList domainPoints(findIndices(dist, domainI));
|
||||
|
||||
// Subset point-wise data.
|
||||
pointField subPoints(points, domainPoints);
|
||||
scalarField subWeights(pointWeights, domainPoints);
|
||||
labelList subPointMap(labelUIndList(pointMap, domainPoints));
|
||||
// Subset point-point addressing (adapt global numbering)
|
||||
labelListList subPointPoints;
|
||||
labelList nOutsideConnections;
|
||||
subsetGlobalCellCells
|
||||
(
|
||||
nCurrDomains,
|
||||
domainI,
|
||||
dist,
|
||||
|
||||
pointPoints,
|
||||
domainPoints,
|
||||
|
||||
subPointPoints,
|
||||
nOutsideConnections
|
||||
);
|
||||
|
||||
decompose
|
||||
(
|
||||
subPointPoints,
|
||||
subPoints,
|
||||
subWeights,
|
||||
subPointMap,
|
||||
*decomposeData.getChild(domainI),
|
||||
domainOffsets[domainI], // The offset for this level and leaf
|
||||
|
||||
finalDecomp
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
multiNodeDecomp::multiNodeDecomp
|
||||
(
|
||||
const dictionary& decompDict,
|
||||
const word& regionName
|
||||
)
|
||||
:
|
||||
decompositionMethod(decompDict, regionName),
|
||||
rootMetadata_()
|
||||
{
|
||||
const dictionary& coeffsDict(
|
||||
findCoeffsDict(
|
||||
typeName + "Coeffs",
|
||||
(selectionType::EXACT | selectionType::MANDATORY)
|
||||
)
|
||||
);
|
||||
initializeMetadata(coeffsDict);
|
||||
}
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
bool multiNodeDecomp::parallelAware() const
|
||||
{
|
||||
return rootMetadata_.parallelAware();
|
||||
}
|
||||
|
||||
|
||||
labelList multiNodeDecomp::decompose
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const pointField& cc,
|
||||
const scalarField& cWeights
|
||||
) const
|
||||
{
|
||||
CompactListList<label> cellCells;
|
||||
calcCellCells(mesh, identity(cc.size()), cc.size(), true, cellCells);
|
||||
|
||||
labelList finalDecomp(cc.size(), Zero);
|
||||
labelList cellMap(identity(cc.size()));
|
||||
|
||||
decompose
|
||||
(
|
||||
cellCells.unpack(),
|
||||
cc,
|
||||
cWeights,
|
||||
cellMap, // map back to original cells
|
||||
rootMetadata_,
|
||||
0,
|
||||
|
||||
finalDecomp
|
||||
);
|
||||
|
||||
return finalDecomp;
|
||||
}
|
||||
|
||||
|
||||
labelList multiNodeDecomp::decompose
|
||||
(
|
||||
const labelListList& globalPointPoints,
|
||||
const pointField& points,
|
||||
const scalarField& pointWeights
|
||||
) const
|
||||
{
|
||||
labelList finalDecomp(points.size(), Zero);
|
||||
labelList pointMap(identity(points.size()));
|
||||
|
||||
decompose
|
||||
(
|
||||
globalPointPoints,
|
||||
points,
|
||||
pointWeights,
|
||||
pointMap, // map back to original points
|
||||
rootMetadata_,
|
||||
0,
|
||||
|
||||
finalDecomp
|
||||
);
|
||||
|
||||
return finalDecomp;
|
||||
}
|
||||
|
||||
// * * * * * * * * * * * * Meta Parser Class * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
|
||||
List<string> multiNodeDecomp::metaParser::getEntries(const dictionary& dict, const string& argument, bool allowWithoutBrackets) {
|
||||
string argumentBracket = argument + "[";
|
||||
DynamicList<string, 4> Result;
|
||||
for(auto& dEntry : dict) {
|
||||
if(dEntry.keyword().starts_with(argumentBracket) || (allowWithoutBrackets && dEntry.keyword() == argument))
|
||||
Result.push_back(dEntry.keyword());
|
||||
}
|
||||
return Result;
|
||||
}
|
||||
|
||||
List<Pair<label>> multiNodeDecomp::metaParser::parseRanges(const string& key) {
|
||||
|
||||
// First, discard the argument and process the indices only.
|
||||
// The current syntax is argument[...]
|
||||
// Assuming that this key was returned in getEntries,
|
||||
// if there is no '[', it is OK and we use the
|
||||
// empty string (update the root).
|
||||
string indices = "";
|
||||
if(key.find_first_of('[') != key.npos) {
|
||||
// There is a '[' in the string.
|
||||
// We can substr from that location.
|
||||
label nFirstBracket = key.find('[');
|
||||
indices = key.substr(nFirstBracket);
|
||||
}
|
||||
|
||||
// All checks print an error message if failed, explaining why.
|
||||
|
||||
DynamicList<Pair<label>, 4> Result;
|
||||
label nCurPtr = 0, nIndicesLength = indices.size();
|
||||
// As long as there are more ranges to parse.
|
||||
while(nCurPtr != nIndicesLength) {
|
||||
// First, check if there is an opening bracket.
|
||||
if(indices[nCurPtr]!='[')
|
||||
FatalError
|
||||
<< "Error when parsing indices "
|
||||
<< indices << ": Expected '[', found "
|
||||
<< indices[nCurPtr] << ". Aborting\n"
|
||||
<< exit(FatalError);
|
||||
|
||||
// Then, find the matching close bracket.
|
||||
label nEndIndex = indices.find(']', nCurPtr);
|
||||
if(nEndIndex == nIndicesLength) {
|
||||
FatalError
|
||||
<< "Error when parsing indices "
|
||||
<< indices << ": Expected ']' after '['. Aborting\n"
|
||||
<< exit(FatalError);
|
||||
}
|
||||
// Read inside the brackets, mark the hyphen if it exists, and make sure
|
||||
// every character is either a digit or a hyphen.
|
||||
// Note that only one hyphen may exist.
|
||||
label nHyphenIdx=-1;
|
||||
for(label nCurIndex = nCurPtr+1; nCurIndex < nEndIndex; ++nCurIndex) {
|
||||
if(!isdigit(indices[nCurIndex])&&indices[nCurIndex]!='-') {
|
||||
FatalError
|
||||
<< "Error when parsing indices "
|
||||
<< indices << ": Expected digit/'-'/']', found "
|
||||
<< indices[nCurIndex] << ". Aborting\n"
|
||||
<< exit(FatalError);
|
||||
}
|
||||
if(indices[nCurIndex]=='-') {
|
||||
if(nHyphenIdx!=-1)
|
||||
FatalError
|
||||
<< "Error when parsing indices "
|
||||
<< indices << ": Found two hyphens(-) inside an index. Aborting\n"
|
||||
<< exit(FatalError);
|
||||
|
||||
nHyphenIdx = nCurIndex;
|
||||
}
|
||||
}
|
||||
label nLeft,nRight;
|
||||
if(nHyphenIdx == -1) {
|
||||
// Not a range - just a single index, or empty brackets (indicating to change the whole range).
|
||||
if(nCurPtr+1==nEndIndex) nLeft = 0, nRight = -1;
|
||||
else {
|
||||
string sNum = indices.substr(nCurPtr+1,nEndIndex-nCurPtr-1);
|
||||
nLeft = nRight = atoi(sNum.c_str());
|
||||
}
|
||||
} else {
|
||||
// A range of indices.
|
||||
// Assert that the hyphen is not right next to the brackets.
|
||||
if(nHyphenIdx+1==nEndIndex||nCurPtr+1==nHyphenIdx)
|
||||
FatalError
|
||||
<< "Error when parsing indices "
|
||||
<< indices << ": Expected number, found "
|
||||
<< (nCurPtr+1==nHyphenIdx?'-':']')
|
||||
<< ". Aborting\n"
|
||||
<< exit(FatalError);
|
||||
|
||||
// Parse the numbers
|
||||
string sLeftNum = indices.substr(nCurPtr+1,nHyphenIdx-nCurPtr-1);
|
||||
string sRightNum = indices.substr(nHyphenIdx+1,nEndIndex-nHyphenIdx-1);
|
||||
nLeft = atoi(sLeftNum.c_str());
|
||||
nRight = atoi(sRightNum.c_str());
|
||||
// Make sure left endpoint is at most the right endpoint
|
||||
if(nLeft>nRight)
|
||||
FatalError
|
||||
<< "Error when parsing indices "
|
||||
<< indices << ": right endpoint("<< nRight
|
||||
<< ") cannot be smaller than left endpoint("
|
||||
<< nLeft << "). Aborting\n"
|
||||
<< exit(FatalError);
|
||||
}
|
||||
// Move the pointer after the closing bracket and append to the result list.
|
||||
nCurPtr = nEndIndex + 1;
|
||||
Result.push_back({nLeft,nRight});
|
||||
}
|
||||
return Result;
|
||||
}
|
||||
|
||||
// * * * * * * * * * * * * Node Metadata Class * * * * * * * * * * * * * //
|
||||
|
||||
void multiNodeDecomp::nodeMetadata::setLeveledDictionaries(const List<const dictionary*>& dictionaries) {
|
||||
setLeveledDictionaries(dictionaries, 0);
|
||||
}
|
||||
|
||||
bool multiNodeDecomp::nodeMetadata::parallelAware() const {
|
||||
// The decomposition tree is parallel aware if and only if all methods used are parallel aware.
|
||||
// If this is a leaf, we are OK.
|
||||
if(children.empty())
|
||||
return true;
|
||||
|
||||
// Otherwise, check if the method used in this node is parallel aware.
|
||||
if(!method->parallelAware())
|
||||
return false;
|
||||
|
||||
// Check recursively, and if any child is not parallel aware - return false.
|
||||
for(auto& child : children)
|
||||
if(!child->parallelAware())
|
||||
return false;
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
void multiNodeDecomp::nodeMetadata::updateProcessorWeights() {
|
||||
label nDom = nDomains();
|
||||
word methodCoeffsName = coeffsDict->get<word>("method") + "Coeffs";
|
||||
// If processorWeights were set by the user, we do not modify them.
|
||||
if(
|
||||
// Check if the user did not specify processorWeights under the coeffs dictionary or the methodCoeffs dictionary
|
||||
!(coeffsDict->subDictOrAdd(methodCoeffsName).found("processorWeights", keyType::LITERAL)
|
||||
|| coeffsDict->subDictOrAdd("coeffs").found("processorWeights", keyType::LITERAL))) {
|
||||
// Then we should compute weights on our own
|
||||
Field<float> processorWeights(nDom);
|
||||
forAll(children, i) {
|
||||
if(children[i]->weight != 1)
|
||||
processorWeights[i] = children[i]->weight;
|
||||
else switch(weightsInitialization) {
|
||||
case RELATIVE:
|
||||
processorWeights[i] = children[i]->size;
|
||||
break;
|
||||
case UNIFORM:
|
||||
processorWeights[i] = 1;
|
||||
break;
|
||||
default:
|
||||
FatalError
|
||||
<< "Weights initialization is not handled in updateProcessorWeights. Aborting\n"
|
||||
<< exit(FatalError);
|
||||
}
|
||||
}
|
||||
|
||||
coeffsDict->subDictOrAdd(methodCoeffsName).add("processorWeights", processorWeights);
|
||||
}
|
||||
}
|
||||
void multiNodeDecomp::nodeMetadata::constructMethods() {
|
||||
// Special handling of nDomains = 1, because some decomposition methods crash when decomposing to one domain.
|
||||
label nDom = nDomains();
|
||||
if(nDom==1) {
|
||||
coeffsDict->clear();
|
||||
coeffsDict->add("method","none");
|
||||
} else
|
||||
updateProcessorWeights();
|
||||
coeffsDict->add("numberOfSubdomains",nDom);
|
||||
|
||||
// Non-verbose construction of decomposition methods would be nice
|
||||
method = decompositionMethod::New(*coeffsDict).release();
|
||||
// Cannot release coeffsDict from memory because method uses a reference that must stay alive
|
||||
|
||||
forAll(children, i) {
|
||||
if(!children[i]->isLeaf())
|
||||
children[i]->constructMethods();
|
||||
}
|
||||
}
|
||||
|
||||
// Recursively construct the decomposition tree, given the list of dimensions and a default method.
|
||||
void multiNodeDecomp::nodeMetadata::constructRecursive(const labelList& dims, const dictionary* defaultMethod) {
|
||||
if(!dims.empty()) {
|
||||
// The list of dimensions of the children is the current list without the first element.
|
||||
labelList newDims(dims.size() - 1);
|
||||
forAll(newDims, i)
|
||||
newDims[i] = dims[i+1];
|
||||
|
||||
// Construct children recursively
|
||||
// First, resize existing children
|
||||
// And delete the excess
|
||||
forAll(children, i) {
|
||||
if(i < dims[0])
|
||||
children[i]->constructRecursive(newDims, defaultMethod);
|
||||
else
|
||||
delete children[i];
|
||||
}
|
||||
label nOldSize = children.size();
|
||||
children.resize(dims[0]);
|
||||
// If the new array is bigger we will need to allocate new children.
|
||||
for(label i = nOldSize; i < dims[0]; ++i)
|
||||
children[i] = new nodeMetadata(newDims, defaultMethod);
|
||||
|
||||
// Compute size (number of leaves in subtree)
|
||||
size = dims[0];
|
||||
if(!children.empty())
|
||||
size *= children[0]->size;
|
||||
}
|
||||
}
|
||||
void multiNodeDecomp::nodeMetadata::updateNodes(const string& key, const std::function<void(nodeMetadata*)>& update) {
|
||||
List<Pair<label>> indicesList = metaParser::parseRanges(key);
|
||||
updateNodes(indicesList, update);
|
||||
}
|
||||
|
||||
// Parse the indices, and apply the update function to all matching nodes.
|
||||
// nCurPtr is used to indicate the index we are now parsing (instead of sending substrings of indices)
|
||||
void multiNodeDecomp::nodeMetadata::updateNodes(const List<Pair<label>>& indices, const std::function<void(nodeMetadata*)>& update, label nCurIdx) {
|
||||
if(nCurIdx == label(indices.size())) update(this);
|
||||
else {
|
||||
// Otherwise, call recursively.
|
||||
label nLeft, nRight, nChildren = children.size();
|
||||
nLeft = indices[nCurIdx].first();
|
||||
nRight = indices[nCurIdx].second();
|
||||
|
||||
// [0,-1] means the entire range.
|
||||
|
||||
if(nLeft==0 && nRight == -1)
|
||||
nRight = nChildren - 1;
|
||||
// Make sure that the indices do not exceed the number of children.
|
||||
if(nRight >= nChildren)
|
||||
FatalError
|
||||
<< "Error when parsing indices: The #" << (nCurIdx+1)
|
||||
<< " range ["<< nLeft <<"," << nRight<<"]:\n"
|
||||
<< " Cannot update indices bigger than number of children("
|
||||
<< nChildren << "). Aborting\n"
|
||||
<< exit(FatalError);
|
||||
|
||||
for(label nChildIdx = nLeft; nChildIdx <= nRight; ++nChildIdx)
|
||||
children[nChildIdx]->updateNodes(indices,update, nCurIdx+1);
|
||||
}
|
||||
// Recompute size assuming children are updated.
|
||||
if(!children.empty()) {
|
||||
size = 0;
|
||||
forAll(children, i)
|
||||
size += children[i]->size;
|
||||
}
|
||||
}
|
||||
|
||||
void multiNodeDecomp::nodeMetadata::setLeveledDictionaries(const List<const dictionary*>& dictionaries, label nLevel) {
|
||||
// Set the dictionary to this level, and to non-leaf children.
|
||||
setDict(*dictionaries[nLevel]);
|
||||
forAll(children, i) {
|
||||
if(children[i]->nDomains() > 0)
|
||||
children[i]->setLeveledDictionaries(dictionaries,nLevel+1);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,435 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2015 OpenFOAM Foundation
|
||||
Copyright (C) 2017-2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
||||
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
||||
for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Class
|
||||
Foam::multiNodeDecomp
|
||||
|
||||
Description
|
||||
Decompose given using consecutive application of decomposers,
|
||||
to perhaps uneven pieces.
|
||||
Note: If uneven pieces are required, the decomposition method
|
||||
used must support the processorWeights argument.
|
||||
|
||||
SourceFiles
|
||||
multiNodeDecomp.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef multiNodeDecomp_H
|
||||
#define multiNodeDecomp_H
|
||||
|
||||
#include "decompositionMethod.H"
|
||||
// #include "List.H"
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class multiNodeDecomp Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class multiNodeDecomp
|
||||
:
|
||||
public decompositionMethod
|
||||
{
|
||||
// Nested classes declarations
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class metaParser Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
// A class responsible for detecting and parsing metadata-related arguments.
|
||||
class metaParser {
|
||||
public:
|
||||
// Detect and return entries related to the given argument.
|
||||
// Input:
|
||||
// dict - the coeffs dictionary we are looking inside.
|
||||
// argument - the argument we're looking for.
|
||||
// allowWithoutBrackets - set to true if the argument can be detected without brackets
|
||||
// For example, domains should not be allowed without brackets, but weightsInitialization can.
|
||||
static List<string> getEntries(const dictionary& dict, const string& argument, bool allowWithoutBrackets = false);
|
||||
|
||||
// Given a key string of an entry returned from getEntries,
|
||||
// Parse and return all ranges described in the key.
|
||||
// Note that it is the user's responsibility to make sure that the right endpoint
|
||||
// does not exceed the number of children in each node.
|
||||
// The user may also write "[]" to specify all children of the node.
|
||||
// In this case, the range returned is [0,-1]. Otherwise, endpoints are always non-negative.
|
||||
// Input:
|
||||
// key - a key of an entry we are parsing, perhaps from an entry returned in getEntries.
|
||||
// Output:
|
||||
// A list of ranges, where the i-th range corresponds to the i-th level of the decomposition tree.
|
||||
// i.e. if two ranges are returned, we will traverse each child of the root in the first range, and recursively
|
||||
// each child in the second range, thus updating nodes at the third level.
|
||||
|
||||
static List<Pair<label>> parseRanges(const string& key);
|
||||
|
||||
};
|
||||
enum WeightsInitialization {
|
||||
RELATIVE,
|
||||
UNIFORM,
|
||||
UNKNOWN
|
||||
};
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class nodeMetadata Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
// A class holding all the information necessary for a multi-node decomposition, without building
|
||||
// the decompositionMethod objects.
|
||||
// The size indicates the number of processors in this subtree. It will be used
|
||||
// When computing the offset of the decomposition, and when using the relative weights initialization.
|
||||
// The weight is a multiplicative factor applied when decomposing. It is 1 by default and can be set by the user.
|
||||
// If the uniform weights initialization is used, all nodes will have the same weight. If the relative weights
|
||||
// initialization is used, each node's weight is set relatively to their size.
|
||||
// Then, the weight field can be used to change the weight of a specific node.
|
||||
// Note that if the coeffs dictionary contains a processorWeights field, it will not be overwritten.
|
||||
// We will then construct a new dictionary with the required numberOfSubdomains and processorWeights.
|
||||
|
||||
class nodeMetadata {
|
||||
public:
|
||||
|
||||
nodeMetadata() : weight(1), size(1), weightsInitialization(UNIFORM), children(0), coeffsDict(nullptr), method(nullptr) {}
|
||||
|
||||
// Constructs a decomposition data tree with dimensions dims and a default method.
|
||||
// Input: A list of domains for each level, and a default dictionary method.
|
||||
nodeMetadata(const labelList& dims, const dictionary* defaultMethod) : nodeMetadata() {
|
||||
initialize(dims, defaultMethod);
|
||||
}
|
||||
|
||||
// Initializes an existing nodeMetadata object.
|
||||
void initialize(const labelList& dims, const dictionary* defaultMethod) {
|
||||
setDict(*defaultMethod);
|
||||
constructRecursive(dims, defaultMethod);
|
||||
}
|
||||
|
||||
~nodeMetadata() {
|
||||
// Since this class represents a tree, we will need to destruct recursively.
|
||||
for(nodeMetadata* child : children)
|
||||
delete child;
|
||||
|
||||
// Only delete method and dict if they were assigned.
|
||||
if(method != nullptr)
|
||||
delete method;
|
||||
if(coeffsDict != nullptr)
|
||||
delete coeffsDict;
|
||||
}
|
||||
|
||||
|
||||
// Getters
|
||||
|
||||
// Get the weight of this node, with respect to the decomposition done in this node's parent
|
||||
label getWeight() const {
|
||||
return weight;
|
||||
}
|
||||
|
||||
// Get the coeffs dictionary for the decomposition of this node.
|
||||
const dictionary* getDict() const {
|
||||
return coeffsDict;
|
||||
}
|
||||
|
||||
// Get the modifiable coeffs dictionary for the decomposition of this node.
|
||||
dictionary* getMutableDict() {
|
||||
return coeffsDict;
|
||||
}
|
||||
|
||||
// Get the number of leaves in this subtree, i.e., number of processors
|
||||
// created under this node.
|
||||
label getSize() const {
|
||||
return size;
|
||||
}
|
||||
|
||||
// Get the decomposition method object of this node.
|
||||
// Note that construct methods must be called first, otherwise
|
||||
// a null pointer will be returned.
|
||||
const Foam::decompositionMethod* getMethod() const {
|
||||
return method;
|
||||
}
|
||||
|
||||
// Get the current weights initialization mode.
|
||||
bool getWeightsInitialization() const {
|
||||
return weightsInitialization;
|
||||
}
|
||||
|
||||
// Get a const pointer to a child of this node.
|
||||
const nodeMetadata* getChild(label index) const {
|
||||
return children[index];
|
||||
}
|
||||
|
||||
// Get a non-const pointer to a child of this node.
|
||||
nodeMetadata* getMutableChild(label index) {
|
||||
return children[index];
|
||||
}
|
||||
|
||||
// Returns the number of direct subdomains this node has.
|
||||
label nDomains() const {
|
||||
return children.size();
|
||||
}
|
||||
|
||||
// Returns whether this node represents a leaf (i.e., has no children)
|
||||
bool isLeaf() const {
|
||||
return children.empty();
|
||||
}
|
||||
// Setters
|
||||
|
||||
// Set the weight of this node, with respect to the decomposition done in this node's parent
|
||||
void setWeight(label weight) {
|
||||
this->weight = weight;
|
||||
}
|
||||
|
||||
// Set the coeffs dictionary for the decomposition of this node.
|
||||
// This creates a copy of the dictionary (and deletes the previous copy created)
|
||||
void setDict(const dictionary& dict) {
|
||||
if(coeffsDict != nullptr) {
|
||||
delete coeffsDict;
|
||||
}
|
||||
coeffsDict = new dictionary(dict);
|
||||
}
|
||||
|
||||
// Set the decomposition method object of this node.
|
||||
void setMethod(Foam::decompositionMethod* method) {
|
||||
this->method = method;
|
||||
}
|
||||
|
||||
// Sets the weights initialization mode. If setRecursive is true, propagate to the entire subtree (i.e., the root and all of the descendents)
|
||||
void setWeightsInitialization(WeightsInitialization newMode, bool setRecursive = true) {
|
||||
weightsInitialization = newMode;
|
||||
if(setRecursive) {
|
||||
for(nodeMetadata* child : children)
|
||||
child->setWeightsInitialization(newMode, true);
|
||||
}
|
||||
}
|
||||
|
||||
// Updates
|
||||
|
||||
// Update the weights of the nodes at the given indices to the given weight.
|
||||
// Input: A string indicating the indices of the nodes to be updated, and
|
||||
// the new weight of the nodes.
|
||||
void updateWeight(const string& indices, label newWeight) {
|
||||
updateNodes(indices, [newWeight](nodeMetadata* node) {
|
||||
node->setWeight(newWeight);
|
||||
});
|
||||
}
|
||||
|
||||
// Update the dimensions array of the nodes at the given indices to the given dimensions array.
|
||||
// Input: A string indicating the indices of the nodes to be updated, and
|
||||
// the new list of dimensions.
|
||||
void updateDomains(const string& indices,const labelList& dims) {
|
||||
updateNodes(indices, [dims](nodeMetadata* node) {
|
||||
// Reconstruct using this node's dict.
|
||||
// Note that first all domain changes are done,
|
||||
// And only then dictionaries are set.
|
||||
// So the descendents' dictionaries are not overwritten.
|
||||
node->constructRecursive(dims,node->getDict());
|
||||
});
|
||||
}
|
||||
|
||||
|
||||
// Update the method of the nodes at the given indices to the given method dictionary.
|
||||
// Input: A string indicating the indices of the nodes to be updated, and
|
||||
// the new method dictionary.
|
||||
void updateMethod(const string& indices, const dictionary& dict) {
|
||||
updateNodes(indices, [dict](nodeMetadata* node) {
|
||||
node->setDict(dict);
|
||||
});
|
||||
}
|
||||
|
||||
// Update the weight initialization mode of nodes at the given indices and their descendents to the new mode.
|
||||
// Input: A string indicating the indices of the nodes to be updated, and
|
||||
// the new weight mode.
|
||||
void updateWeightsInitialization(const string& indices, WeightsInitialization newMode) {
|
||||
updateNodes(indices, [newMode](nodeMetadata* node) {
|
||||
node->setWeightsInitialization(newMode);
|
||||
});
|
||||
}
|
||||
|
||||
// Given a list of dictionaries for each level, set the dictionaries accordingly.
|
||||
// Input: A list of dictionaries for each level.
|
||||
void setLeveledDictionaries(const List<const dictionary*>& dictionaries);
|
||||
|
||||
// To be used within the decompositionMethod's parallelAware function.
|
||||
// Returns whether all decompositions in this subtree are parallel aware
|
||||
// (i.e., synchronize domains across proc boundaries)
|
||||
bool parallelAware() const;
|
||||
|
||||
// Calculate (and add to the dictionary) the new processor weights if reqired,
|
||||
// Using the children's weights and the weight initialization mode.
|
||||
void updateProcessorWeights();
|
||||
|
||||
// Construct the decompositionMethod object for this node and all its descendents.
|
||||
void constructMethods();
|
||||
|
||||
private:
|
||||
|
||||
// The weight of this node in the parent's decomposition, relative to the other nodes.
|
||||
// Overrides weights set by the weights initialization.
|
||||
label weight;
|
||||
|
||||
// The size of the node indicates the total number of subdomains this node has.
|
||||
label size;
|
||||
|
||||
// An enum describing the weights initialization.
|
||||
WeightsInitialization weightsInitialization;
|
||||
|
||||
// The direct descendents.
|
||||
List<nodeMetadata*> children;
|
||||
|
||||
// The dictionary used to construct the decomposition method.
|
||||
dictionary* coeffsDict;
|
||||
|
||||
// The decomposition method of this node.
|
||||
const Foam::decompositionMethod* method;
|
||||
|
||||
// Recursively constructs the subtree rooted at this node
|
||||
// Input: A list of dimensions and the dictionary of the default method.
|
||||
void constructRecursive(const labelList& dims, const dictionary* defaultMethod);
|
||||
|
||||
// Update all nodes matching the given indices with the given updating function.
|
||||
// Input: A list of ranges for each level, and a function that receives a pointer to nodeMetadata,
|
||||
// updates it accordingly and returns nothing.
|
||||
void updateNodes(const string& key, const std::function<void(nodeMetadata*)>& update);
|
||||
|
||||
// Internal implementation of updateNodes.
|
||||
// The list of ranges are constructed by passing the key argument to the meta parser.
|
||||
// nCurIdx is an internal variable that indicates our location inside the indices array.
|
||||
void updateNodes(const List<Pair<label>>& indices, const std::function<void(nodeMetadata*)>& update, label nCurIdx = 0);
|
||||
|
||||
// This function is used inside the public setLeveledDictionaries function.
|
||||
void setLeveledDictionaries(const List<const dictionary*>& dictionaries, label nLevel);
|
||||
|
||||
// Parse the range of indices starting at the (string) index nStartIndex.
|
||||
// Input: The indices string, and the starting position of the range
|
||||
// (i.e, the position of the opening bracket)
|
||||
// Returns a pair representing the range if succeeded,
|
||||
// or crashes the program with an appropriate error message if failed to parse.
|
||||
Pair<label> parseRange(const string& indices, label nStartIndex) const;
|
||||
|
||||
};
|
||||
|
||||
|
||||
// Private Data
|
||||
|
||||
//- The decomposition metadata.
|
||||
nodeMetadata rootMetadata_;
|
||||
|
||||
|
||||
// Private Member Functions
|
||||
|
||||
//- Read coeffsDict and construct the decomposition metadata.
|
||||
void initializeMetadata(const dictionary& coeffsDict);
|
||||
|
||||
|
||||
//- Given connectivity across processors work out connectivity
|
||||
// for a (consistent) subset
|
||||
void subsetGlobalCellCells
|
||||
(
|
||||
const label nDomains,
|
||||
const label domainI,
|
||||
const labelList& dist,
|
||||
|
||||
const labelListList& cellCells,
|
||||
const labelList& set,
|
||||
labelListList& subCellCells,
|
||||
labelList& cutConnections
|
||||
) const;
|
||||
|
||||
//- Decompose at 'currLevel' without addressing
|
||||
void decompose
|
||||
(
|
||||
const labelListList& pointPoints,
|
||||
const pointField& points,
|
||||
const scalarField& pointWeights,
|
||||
const labelUList& pointMap, // map back to original points
|
||||
const nodeMetadata& decomposeData,
|
||||
const label leafOffset,
|
||||
|
||||
labelList& finalDecomp
|
||||
) const;
|
||||
|
||||
|
||||
//- No copy construct
|
||||
multiNodeDecomp(const multiNodeDecomp&) = delete;
|
||||
|
||||
//- No copy assignment
|
||||
void operator=(const multiNodeDecomp&) = delete;
|
||||
|
||||
public:
|
||||
|
||||
//- Runtime type information
|
||||
TypeName("multiNode");
|
||||
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct given decomposition dictionary and optional region name
|
||||
explicit multiNodeDecomp
|
||||
(
|
||||
const dictionary& decompDict,
|
||||
const word& regionName = ""
|
||||
);
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~multiNodeDecomp() = default;
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Is method parallel aware?
|
||||
// i.e. does it synchronize domains across proc boundaries
|
||||
virtual bool parallelAware() const;
|
||||
|
||||
//- Inherit decompose from decompositionMethod
|
||||
using decompositionMethod::decompose;
|
||||
|
||||
//- Return for every coordinate the wanted processor number.
|
||||
// Use the mesh connectivity (if needed)
|
||||
virtual labelList decompose
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const pointField& points,
|
||||
const scalarField& pointWeights
|
||||
) const;
|
||||
|
||||
//- Return for every coordinate the wanted processor number.
|
||||
// Explicitly provided connectivity - does not use mesh_.
|
||||
virtual labelList decompose
|
||||
(
|
||||
const labelListList& globalCellCells,
|
||||
const pointField& cc,
|
||||
const scalarField& cWeights
|
||||
) const;
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -5,5 +5,6 @@ randomRenumber/randomRenumber.C
|
||||
springRenumber/springRenumber.C
|
||||
structuredRenumber/structuredRenumber.C
|
||||
structuredRenumber/OppositeFaceCellWaveBase.C
|
||||
hpathRenumber/hpathRenumber.C
|
||||
|
||||
LIB = $(FOAM_LIBBIN)/librenumberMethods
|
||||
|
||||
976
src/renumber/renumberMethods/hpathRenumber/hpathRenumber.C
Normal file
976
src/renumber/renumberMethods/hpathRenumber/hpathRenumber.C
Normal file
@ -0,0 +1,976 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2023 Alon Zameret, Noam Manaker Morag
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
||||
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
||||
for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "hpathRenumber.H"
|
||||
#include "addToRunTimeSelectionTable.H"
|
||||
|
||||
#include "CircularBuffer.H"
|
||||
#include <iomanip>
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
defineTypeNameAndDebug(hpathRenumber, 0);
|
||||
|
||||
addToRunTimeSelectionTable
|
||||
(
|
||||
renumberMethod,
|
||||
hpathRenumber,
|
||||
dictionary
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Local Details * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class hpathRenumber::hpathFinder Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
// Class hpathFinder Declaration
|
||||
class hpathRenumber::hpathFinder
|
||||
{
|
||||
// Private Data
|
||||
|
||||
// The input mesh
|
||||
const Foam::polyMesh& mesh;
|
||||
|
||||
// Counter for the number of renumbered cells
|
||||
label nFoundCellCount;
|
||||
|
||||
// Marks cells that have been added to the renumbering as 'true'
|
||||
Foam::bitSet bIsRenumbered;
|
||||
|
||||
// For every data structure, I explain what it is used for and which
|
||||
// method is used to compute it
|
||||
|
||||
// For every cell, a list of all its point-neighbouring cells - getMeshGraph()
|
||||
Foam::labelListList nMeshGraph;
|
||||
|
||||
// For each cell its 'layer index': - getLayers()
|
||||
// - cells in the same layer will have the same layer index
|
||||
Foam::labelList nCellLayerIndex;
|
||||
|
||||
// For each cell its 'connected component index': - getConnectedComponents()
|
||||
// - cells in the same connected component will have the same connected component index
|
||||
Foam::labelList nCellConnnectedComponentIndex;
|
||||
|
||||
// For each cell its face-neighbours in the connected component - getConnnectedComponentGraph()
|
||||
Foam::labelListList nConnnectedComponentGraph;
|
||||
|
||||
// Marks cells that have already been by BFS - getStartingCellInConnnectedComponent()
|
||||
Foam::bitSet bBFSFoundCell;
|
||||
|
||||
// For each cell its point distance from the connected components start cell - reorderDistFromStart()
|
||||
Foam::labelList nCellPointDistFromStart;
|
||||
Foam::labelList nCellFaceDistFromStart;
|
||||
|
||||
// For each cell its DFS depth within the connected component - findPath()
|
||||
Foam::labelList nDFSCellDepth;
|
||||
|
||||
// For each cell its DFS parent within the connected component - findPath()
|
||||
Foam::labelList nDFSParentCell;
|
||||
|
||||
|
||||
public: // Public methods
|
||||
|
||||
//- Constructor
|
||||
explicit hpathFinder(const Foam::polyMesh& mesh);
|
||||
|
||||
// Member Functions
|
||||
|
||||
// Get Renumbering for the mesh
|
||||
void getRenumbering
|
||||
(
|
||||
Foam::labelList& cellOrder,
|
||||
bool bApplyLayerSeparation
|
||||
);
|
||||
|
||||
// Returns the accuracy of the renumbering:
|
||||
// Accuracy is defined as the percentage of consecutive cells that are also face-neighbours in the mesh
|
||||
// - Cells are face-neighbours if they have a common face
|
||||
float getAccuracy(const Foam::labelList& cellOrder) const;
|
||||
|
||||
|
||||
private: // Private methods
|
||||
|
||||
// Initialize data structures for later use
|
||||
void initialize();
|
||||
|
||||
// Given a cell and a face index, find its neighbour through the face
|
||||
// - If facing the boundary, returns -1
|
||||
// - Otherwise, returns FaceOwner/FaceNeighbour[nFaceIdx], the one that's different from nCellIdx
|
||||
label getNei(label nCellIdx, label nFaceIdx) const;
|
||||
|
||||
// Creates the 'Mesh-Graph': for every cell, a list of cells that are point-neighbours with it in the mesh
|
||||
// - Cells are point-neighbours if they have a common point
|
||||
void getMeshGraph();
|
||||
|
||||
// Separates the mesh into layers: each cell has its layer saved in nCellLayerIndex
|
||||
// Also returns for every layer a list of all cells in it
|
||||
void getLayerSeparation(Foam::DynamicList<Foam::DynamicList<label>>& nCellsByLayer);
|
||||
|
||||
// For every connected component, find the deepest cell and choose it as a starting cell
|
||||
// Returns a list with one starting cell per connected component
|
||||
void getStartingCells(Foam::DynamicList<label>& nStartingCells) const;
|
||||
|
||||
// Once the starting cells have been found, this method does the actual layer separation
|
||||
void getLayers(const Foam::labelList& nStartingCellList, Foam::DynamicList<Foam::DynamicList<label>>& nCellsByLayer);
|
||||
|
||||
// Renumber all cells in a layer
|
||||
void solveLayer(Foam::labelList nCellsInLayer, Foam::labelList& cellOrder);
|
||||
|
||||
// Seperates cells into face-connected components
|
||||
// This is done using a general DFS algorithm
|
||||
void getConnectedComponents(const Foam::labelList& nCellList, Foam::DynamicList<Foam::DynamicList<label>>& nCellsByConnnectedComponent);
|
||||
|
||||
// Find an approximate H-path through a connected component
|
||||
void solveConnectedComponent(const Foam::labelList& nCellsInConnectedComponent, Foam::labelList& cellOrder);
|
||||
|
||||
// Finds for each cell in the connected component a list of its face-neighbours within the component
|
||||
void getConnnectedComponentGraph(const Foam::labelList& nCellsInConnectedComponent);
|
||||
|
||||
// Finds a starting cell within the connected component
|
||||
label getStartingCellInConnnectedComponent(const Foam::labelList& nCellsInConnectedComponent);
|
||||
|
||||
// This method reorders the cells in the given connected component based on their distance from nStartCell
|
||||
void reorderDistFromStart(label nStartCell, const Foam::labelList& nCellsInConnectedComponent);
|
||||
|
||||
// Finds an H-path within the connected component and returns it in nResultHpath
|
||||
// H-path is guaranteed to start at nStartCell
|
||||
void findPath(label nStartCell, Foam::labelList& cellOrder);
|
||||
|
||||
// Resets data structures for the cells that weren't found
|
||||
void resetCells(Foam::labelList& nCellList);
|
||||
};
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::hpathRenumber::hpathRenumber(const dictionary& dict)
|
||||
:
|
||||
renumberMethod(dict),
|
||||
applyLayerSeparation_
|
||||
(
|
||||
dict.optionalSubDict(typeName + "Coeffs")
|
||||
.getOrDefault("layered", true)
|
||||
)
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::labelList Foam::hpathRenumber::renumber
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const pointField& points
|
||||
) const
|
||||
{
|
||||
Info<< nl
|
||||
<< "Starting Cell Renumbering: "
|
||||
<< mesh.nCells() << " Cells, "
|
||||
<< mesh.nFaces() << " Faces, "
|
||||
<< mesh.nPoints() << " Points" << nl << nl;
|
||||
|
||||
hpathRenumber::hpathFinder s(mesh);
|
||||
Foam::labelList cellOrder;
|
||||
s.getRenumbering(cellOrder, applyLayerSeparation_);
|
||||
|
||||
std::cout << std::endl << "Found path with accuracy of "
|
||||
<< std::fixed << std::setprecision(3)
|
||||
<< s.getAccuracy(cellOrder) << "%" << std::endl;
|
||||
|
||||
fprintf(stdout,"\n******************************************************\n\n\n");
|
||||
return cellOrder;
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * hpathRenumber::hpathFinder * * * * * * * * * * * //
|
||||
|
||||
// Public methods
|
||||
Foam::hpathRenumber::hpathFinder::hpathFinder(const polyMesh& mesh)
|
||||
:
|
||||
mesh(mesh)
|
||||
{}
|
||||
|
||||
|
||||
void Foam::hpathRenumber::hpathFinder::getRenumbering
|
||||
(
|
||||
Foam::labelList& cellOrder,
|
||||
bool bApplyLayerSeparation
|
||||
)
|
||||
{
|
||||
// Find a renumbering for the entire mesh
|
||||
|
||||
cellOrder.resize(mesh.nCells());
|
||||
// Counter for how many cells we have added to the renumbering so far
|
||||
nFoundCellCount = 0;
|
||||
|
||||
// Initialize the data structures:
|
||||
Info<< "Initializing Data Structures" << endl;
|
||||
initialize();
|
||||
|
||||
// Compute a graph to represent the mesh:
|
||||
// - Cells in the mesh will be connected in the graph if they have a *common point*
|
||||
getMeshGraph();
|
||||
|
||||
Foam::DynamicList<Foam::DynamicList<label>> nCellsByLayer;
|
||||
if (bApplyLayerSeparation)
|
||||
{
|
||||
getLayerSeparation(nCellsByLayer);
|
||||
}
|
||||
else
|
||||
{
|
||||
// If there is no layer separation, set the entire mesh as one 'layer'
|
||||
Foam::DynamicList<label> nAllCells(Foam::identity(mesh.nCells()));
|
||||
nCellsByLayer.append(nAllCells);
|
||||
}
|
||||
|
||||
std::cout << "Beginning Hpath Computation" << std::endl;
|
||||
// Find H-path for each layer separately
|
||||
for (const Foam::labelList& nCellsInLayer : nCellsByLayer)
|
||||
{
|
||||
// solveLayer() will find a renumbering for all the cells in the layer
|
||||
// Path will be appended into cellOrder
|
||||
solveLayer(nCellsInLayer, cellOrder);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
float Foam::hpathRenumber::hpathFinder::getAccuracy
|
||||
(
|
||||
const Foam::labelList& cellOrder
|
||||
) const
|
||||
{
|
||||
// Finding the number of "hits"
|
||||
// - A hit is a pair of consecutive cells in the renumbering that are also face-neighbours in the mesh
|
||||
// - Cells are face-neighbours if they have a common face
|
||||
label nHitCnt = 0;
|
||||
for (label nPathIdx = 0; nPathIdx < label(cellOrder.size()) - 1; nPathIdx++)
|
||||
{
|
||||
label nCurrCellIdx = cellOrder[nPathIdx];
|
||||
label nNextCellIdx = cellOrder[nPathIdx+1];
|
||||
|
||||
// For every pair of consecutive cells, we search for a common face between them
|
||||
for (label nFaceIdx : mesh.cells()[nCurrCellIdx])
|
||||
{
|
||||
label nNeiIdx = getNei(nCurrCellIdx, nFaceIdx);
|
||||
// If there is a common face between them, we add a hit!
|
||||
if (nNeiIdx == nNextCellIdx)
|
||||
{
|
||||
nHitCnt++;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
// The accuracy is the percentage of consecutive cells that were hits
|
||||
return 100.0f * float(nHitCnt) / float(label(cellOrder.size()) - 1);
|
||||
}
|
||||
|
||||
|
||||
// Private methods
|
||||
void Foam::hpathRenumber::hpathFinder::initialize()
|
||||
{
|
||||
// Data structure to keep track of cells already in the renumbering
|
||||
bIsRenumbered = Foam::bitSet(mesh.nCells(), false);
|
||||
|
||||
// List used to seperate cells into layers
|
||||
// Saves for every cell its layer index
|
||||
nCellLayerIndex = Foam::labelList(mesh.nCells(), -1);
|
||||
|
||||
// List used to seperate cells within the same layer into seperate connected components
|
||||
// Saves for every cell its connected component index
|
||||
nCellConnnectedComponentIndex = Foam::labelList(mesh.nCells(), -1);
|
||||
|
||||
// Marks cells that have already been by BFS
|
||||
bBFSFoundCell = Foam::bitSet(mesh.nCells(), false);
|
||||
// Saves for every cell its face-neighbours within the same connected component
|
||||
nConnnectedComponentGraph = Foam::labelListList(mesh.nCells());
|
||||
|
||||
// Saves for every cell within a layer its point-distance from the starting cell of that layer
|
||||
nCellPointDistFromStart = Foam::labelList(mesh.nCells(), -1);
|
||||
// Saves for every cell within a layer its face-distance from the starting cell of that layer
|
||||
nCellFaceDistFromStart = Foam::labelList(mesh.nCells(), -1);
|
||||
|
||||
// For each cell its DFS depth within the connected component
|
||||
nDFSCellDepth = Foam::labelList(mesh.nCells(), -1);
|
||||
// For each cell its DFS parent within the connected component
|
||||
nDFSParentCell = Foam::labelList(mesh.nCells(), -1);
|
||||
|
||||
// Now we can find the Mesh-Graph
|
||||
nMeshGraph = Foam::labelListList(mesh.nCells());
|
||||
}
|
||||
|
||||
|
||||
Foam::label Foam::hpathRenumber::hpathFinder::getNei(label nCell, label nFaceIdx) const
|
||||
{
|
||||
// If the face has no neighbor, return -1
|
||||
if (nFaceIdx >= mesh.faceNeighbour().size())
|
||||
{
|
||||
return -1;
|
||||
}
|
||||
|
||||
// Otherwise, the face connects 2 cells: the owner and the neighbor.
|
||||
// One of these should be nCell
|
||||
label nOwner = mesh.faceOwner()[nFaceIdx];
|
||||
label nNei = mesh.faceNeighbour()[nFaceIdx];
|
||||
|
||||
// Find which of these two cells is the input cell, return the other one
|
||||
return (nOwner == nCell) ? nNei : nOwner;
|
||||
}
|
||||
|
||||
|
||||
void Foam::hpathRenumber::hpathFinder::getMeshGraph()
|
||||
{
|
||||
// The purpose of the bitSet is to avoid adding the same cell as a neighbour more than once
|
||||
Foam::bitSet foundCells(mesh.nCells(), false);
|
||||
|
||||
for (label nCellIdx = 0; nCellIdx < mesh.nCells(); nCellIdx++)
|
||||
{
|
||||
// Dynamic list of point-neighbours of the current cell
|
||||
DynamicList<label> nNeiList;
|
||||
|
||||
for (label nPntIdx : mesh.cellPoints()[nCellIdx])
|
||||
{
|
||||
for (label nNeiCell : mesh.pointCells()[nPntIdx])
|
||||
{
|
||||
if (nNeiCell == nCellIdx) continue;
|
||||
if (foundCells[nNeiCell]) continue;
|
||||
|
||||
// We have found a cell with a common point to the current cell
|
||||
// Therefore, we can now add it as a neighbour in the Mesh-Graph
|
||||
nNeiList.append(nNeiCell);
|
||||
|
||||
foundCells[nNeiCell] = true;
|
||||
}
|
||||
}
|
||||
|
||||
// Convert from DynamicList to labelList
|
||||
nMeshGraph[nCellIdx] = nNeiList;
|
||||
|
||||
// Clear the bitSet before next iteration
|
||||
for (label neiCell : nMeshGraph[nCellIdx])
|
||||
{
|
||||
foundCells[neiCell] = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::hpathRenumber::hpathFinder::getLayerSeparation
|
||||
(
|
||||
Foam::DynamicList<Foam::DynamicList<label>>& nCellsByLayer
|
||||
)
|
||||
{
|
||||
// We want to separate the mesh into layers:
|
||||
// 1) The mesh is separated into connected components - getConnectedComponents()
|
||||
// 2) For each connected component we find the deepest cell and set it to be a starting cell - getStartingCells()
|
||||
// 3) Each connected component is separated into layers using a BFS from the starting - getLayers()
|
||||
|
||||
std::cout << "Finding Layer Separation" << std::endl;
|
||||
|
||||
// Step 1: Finding connected components
|
||||
Foam::labelList nAllCells(Foam::identity(mesh.nCells()));
|
||||
Foam::DynamicList<Foam::DynamicList<label>> nCellsByConnnectedComponent;
|
||||
getConnectedComponents(nAllCells, nCellsByConnnectedComponent);
|
||||
|
||||
// Step 2: Finding starting cells
|
||||
// Finds for each connected component its deepest cell and returns them
|
||||
// The starting cells are returned through nStartingCellList
|
||||
// There will be exactly one starting cell per connected component
|
||||
Foam::DynamicList<label> nStartingCellList;
|
||||
getStartingCells(nStartingCellList);
|
||||
|
||||
// Step 3: Layer separation
|
||||
// Layer separation is done in each connected component from the starting cell outwards
|
||||
getLayers(nStartingCellList, nCellsByLayer);
|
||||
|
||||
label nLayerCnt = nCellsByLayer.size();
|
||||
std::cout << "Mesh Separated into " << nLayerCnt << " layers" << std::endl;
|
||||
|
||||
// Reset connected component index list for solveLayer() to use
|
||||
nCellConnnectedComponentIndex = Foam::labelList(mesh.nCells(), -1);
|
||||
}
|
||||
|
||||
|
||||
void Foam::hpathRenumber::hpathFinder::getStartingCells
|
||||
(
|
||||
Foam::DynamicList<label>& nStartingCells)
|
||||
const
|
||||
{
|
||||
// The starting cell in each connected component should be the 'deepest' cell in that connected component
|
||||
// A cells 'depth' is its point-distance from any boundary cell
|
||||
|
||||
// The Algorithm:
|
||||
// 1) Find all boundary points
|
||||
// 2) Find all boundary cells
|
||||
// 3) Separate the boundary cells by connected component
|
||||
// 4) Find for each cell in the mesh its point-distance from the boundary
|
||||
// 5) Return through nStartingCells a list containing the deepest cell from each connected component
|
||||
|
||||
// nCellDepthList will contain for each cell its depth within its ConnectedComponent
|
||||
Foam::labelList nCellDepthList(mesh.nCells(), -1);
|
||||
|
||||
// Step 1: Finding all boundary points
|
||||
// - A boundary point is a point on a boundary face
|
||||
Foam::bitSet bIsBoundaryPts(mesh.nPoints(), false);
|
||||
|
||||
// Iterate over boundary faces and mark their points as boundary
|
||||
// - Ignore boundary faces of type "empty"
|
||||
const polyBoundaryMesh& bndMesh = mesh.boundaryMesh();
|
||||
for (label nBndType = 0; nBndType < bndMesh.size(); ++nBndType)
|
||||
{
|
||||
// If boundary is of type "empty" - skip it
|
||||
if (bndMesh[nBndType].type().compare("empty") == 0) continue;
|
||||
|
||||
labelRange range = bndMesh.patchRanges()[nBndType];
|
||||
// For every face in range set all of its points as boundary
|
||||
for (label nFaceIdx = range.min(); nFaceIdx <= range.max(); ++nFaceIdx)
|
||||
{
|
||||
for (label nPointIdx : mesh.faces()[nFaceIdx])
|
||||
{
|
||||
bIsBoundaryPts[nPointIdx] = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
// Step 2: Finding all boundary cells
|
||||
// - A boundary cell is a cell containing at least one boundary point
|
||||
Foam::bitSet bIsBoundaryCells(mesh.nCells(), false);
|
||||
|
||||
// For every boundary point: mark all cells that have it as boundary cells
|
||||
for (label nPntIdx = 0; nPntIdx < mesh.nPoints(); nPntIdx++)
|
||||
{
|
||||
if (!bIsBoundaryPts[nPntIdx]) continue;
|
||||
// Here we make use of the nPntCellList we found when computing the Mesh-Graph
|
||||
for (label nCellIdx : mesh.pointCells()[nPntIdx])
|
||||
{
|
||||
bIsBoundaryCells[nCellIdx] = true;
|
||||
}
|
||||
}
|
||||
|
||||
// Step 3: Separate the boundary cells based on which connected component they are a part of
|
||||
// - Find for each connected component a list of all its boundary cells
|
||||
label nConnnectedComponentCnt = *std::max_element(nCellConnnectedComponentIndex.begin(), nCellConnnectedComponentIndex.end()) + 1;
|
||||
Foam::List<Foam::DynamicList<label>> nBoundaryCellsByConnnectedComponent(nConnnectedComponentCnt);
|
||||
for (label nCellIdx = 0; nCellIdx < mesh.nCells(); nCellIdx++)
|
||||
{
|
||||
if (bIsBoundaryCells[nCellIdx])
|
||||
{
|
||||
nBoundaryCellsByConnnectedComponent[nCellConnnectedComponentIndex[nCellIdx]].append(nCellIdx);
|
||||
}
|
||||
}
|
||||
|
||||
// Steps 4+5:
|
||||
// - Find each cells distance from the boundary by running a BFS algorithm from the boundary in each connected component
|
||||
// - Along the way, find the deepest cell in each connected component and push them to the list
|
||||
|
||||
for (Foam::labelList& nBndCells : nBoundaryCellsByConnnectedComponent)
|
||||
{
|
||||
// We want to find the maximum depth cell within the connected component
|
||||
label nMaxDepthCell = nBndCells[0];
|
||||
|
||||
// Run a BFS algorithm from the boundary:
|
||||
// - All boundary cells are pushed to the queue with depth 0
|
||||
Foam::CircularBuffer<label> nBfsQueue;
|
||||
for (label nCellIdx : nBndCells)
|
||||
{
|
||||
nCellDepthList[nCellIdx] = 0;
|
||||
nBfsQueue.push_back(nCellIdx);
|
||||
}
|
||||
|
||||
while (!nBfsQueue.empty())
|
||||
{
|
||||
label nCurrCell = nBfsQueue.first();
|
||||
nBfsQueue.pop_front();
|
||||
|
||||
// Check if current cell is the deeper than the maximum depth cell found so far. If it is, we update nMaxDepthCell
|
||||
if (nCellDepthList[nCurrCell] > nCellDepthList[nMaxDepthCell])
|
||||
{
|
||||
nMaxDepthCell = nCurrCell;
|
||||
}
|
||||
|
||||
// The BFS algorithm needs to be based on point-neghbours, meaning by using the Mesh-Graph
|
||||
// Note that while cells in different connected components are never face-neighbours, but they may be point-neighbours (neighbours in the Mesh-Graph)
|
||||
// Therefore, when pushing all point-neighbouring cells we need to check that they are in the same connected component
|
||||
for (label nNeiCell : nMeshGraph[nCurrCell])
|
||||
{
|
||||
if (nCellDepthList[nNeiCell] != -1) continue;
|
||||
if (nCellConnnectedComponentIndex[nNeiCell] != nCellConnnectedComponentIndex[nCurrCell]) continue;
|
||||
|
||||
nCellDepthList[nNeiCell] = nCellDepthList[nCurrCell] + 1;
|
||||
nBfsQueue.push_back(nNeiCell);
|
||||
}
|
||||
}
|
||||
|
||||
// Finally, we can push the maximum depth cell we found
|
||||
nStartingCells.append(nMaxDepthCell);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::hpathRenumber::hpathFinder::getLayers
|
||||
(
|
||||
const Foam::labelList& nStartingCells,
|
||||
Foam::DynamicList<Foam::DynamicList<label>>& nCellsByLayer
|
||||
)
|
||||
{
|
||||
// Separates all cells in the mesh to layers
|
||||
// In each connected component:
|
||||
// 1) The starting cell is set as 'Layer 0'
|
||||
// 2) All cells that are point-neighbours of the starting cell are set as layer 1
|
||||
// 3) All cells that are point-neighbours of a cell in layer 1 are set as layer 2
|
||||
// 4) Repeat step (3) with increasing layer indices until all cells have been found
|
||||
// In this way, cells in each component are grouped in layers by their MINIMUM point-distance to their starting cell
|
||||
|
||||
// Note: In reality, the same BFS is run for all components at once. Because components are connected, this is equivalent
|
||||
|
||||
// Now we run BFS from starting cells
|
||||
Foam::CircularBuffer<label> nBfsQueue;
|
||||
|
||||
// Push all starting cells to queue
|
||||
for (label nCellIdx : nStartingCells)
|
||||
{
|
||||
nCellLayerIndex[nCellIdx] = 0;
|
||||
nBfsQueue.push_back(nCellIdx);
|
||||
}
|
||||
|
||||
// BFS algorithm will find for each cell its minimum distance in the Mesh-Graph from the corresponding starting cell
|
||||
while (!nBfsQueue.empty())
|
||||
{
|
||||
label nCurrCell = nBfsQueue.first();
|
||||
nBfsQueue.pop_front();
|
||||
|
||||
label nLayer = nCellLayerIndex[nCurrCell];
|
||||
if (label(nCellsByLayer.size()) <= nCellLayerIndex[nCurrCell])
|
||||
{
|
||||
nCellsByLayer.append(Foam::labelList());
|
||||
}
|
||||
nCellsByLayer[nLayer].append(nCurrCell);
|
||||
|
||||
for (label nNeiCell : nMeshGraph[nCurrCell])
|
||||
{
|
||||
// If neighbour is in a different connected component, ignore it
|
||||
if (nCellConnnectedComponentIndex[nNeiCell] != nCellConnnectedComponentIndex[nCurrCell]) continue;
|
||||
// If neighbour hasn't been visited yet, set it's layer and push it:
|
||||
if (nCellLayerIndex[nNeiCell] != -1) continue;
|
||||
nCellLayerIndex[nNeiCell] = nCellLayerIndex[nCurrCell] + 1;
|
||||
nBfsQueue.push_back(nNeiCell);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::hpathRenumber::hpathFinder::solveLayer
|
||||
(
|
||||
Foam::labelList nCellsInLayer,
|
||||
Foam::labelList& cellOrder
|
||||
)
|
||||
{
|
||||
// Note: this method assumes all cells in the nCellsInLayer have the same 'layer index' in nCellLayerIndex
|
||||
|
||||
// General rundown of the algorithm:
|
||||
// 1) Cells are separated into connected components
|
||||
// 2) Each connected component is solved separately and its path is added to the renumbering
|
||||
// 3) If there are still cells that have not been found, return to step (1)
|
||||
|
||||
while (!nCellsInLayer.empty())
|
||||
{
|
||||
// Step 1: separate the cells into connected components
|
||||
// - Connected components need to be connected by FACES (not points)
|
||||
Foam::DynamicList<Foam::DynamicList<label>> nCellsByConnnectedComponent;
|
||||
getConnectedComponents(nCellsInLayer, nCellsByConnnectedComponent);
|
||||
|
||||
label nOrigFoundCellCount = nFoundCellCount;
|
||||
|
||||
// Step 2: For each connected component we call getHpathinConnnectedComponent()
|
||||
for (const Foam::labelList& nCellsInConnectedComponent : nCellsByConnnectedComponent)
|
||||
{
|
||||
solveConnectedComponent(nCellsInConnectedComponent, cellOrder);
|
||||
}
|
||||
// Find how many cells were still not found
|
||||
label nRemainingCellCount = nCellsInLayer.size() - (nFoundCellCount - nOrigFoundCellCount);
|
||||
if (nRemainingCellCount == 0) break;
|
||||
|
||||
// Step 3: Find all the cells in the layer that were not found yet
|
||||
Foam::labelList nRemainingCells(nRemainingCellCount);
|
||||
label nIndex = 0;
|
||||
for (label nCellIdx : nCellsInLayer)
|
||||
{
|
||||
if (!bIsRenumbered[nCellIdx])
|
||||
{
|
||||
nRemainingCells[nIndex++] = nCellIdx;
|
||||
}
|
||||
}
|
||||
|
||||
resetCells(nRemainingCells);
|
||||
nCellsInLayer = nRemainingCells;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::hpathRenumber::hpathFinder::getConnectedComponents
|
||||
(
|
||||
const Foam::labelList& nCellList,
|
||||
Foam::DynamicList<Foam::DynamicList<label>>& nCellsByConnnectedComponent
|
||||
)
|
||||
{
|
||||
// This function separates cells in a certain layer into connected components
|
||||
// Each connected component will get a unique index
|
||||
// - nCellsByConnnectedComponent[i] will contain a list of all cells with connected component index 'i'
|
||||
|
||||
// We accomplish this using a generalized DFS algorithm:
|
||||
// While there are still cells that haven't been found:
|
||||
// 1) Choose some cell in the layer that hasn't been found yet
|
||||
// 2) Find all cells connected to that cell (using a dfs stack)
|
||||
// 3) Set all of these cells as a new connected component
|
||||
|
||||
// Notice: This method does not use the Mesh-Graph!
|
||||
// - This is because the Mesh-Graph is POINT-connected, but we are searching for FACE-connected components
|
||||
|
||||
// The index of the current connected component
|
||||
label nConnnectedComponentIndex = -1;
|
||||
|
||||
for (label nCellIdx : nCellList)
|
||||
{
|
||||
if (nCellConnnectedComponentIndex[nCellIdx] != -1) continue;
|
||||
|
||||
// This cell hasn't been found yet
|
||||
// So we set it to be the beginning of a new connected component
|
||||
// We increment nConnnectedComponentIndex, it is now the index of this new connected component
|
||||
nConnnectedComponentIndex++;
|
||||
nCellConnnectedComponentIndex[nCellIdx] = nConnnectedComponentIndex;
|
||||
nCellsByConnnectedComponent.append(Foam::labelList());
|
||||
|
||||
Foam::CircularBuffer<label> nDfsStack;
|
||||
nDfsStack.push_back(nCellIdx);
|
||||
|
||||
while(!nDfsStack.empty())
|
||||
{
|
||||
label nCurrCell = nDfsStack.last();
|
||||
nDfsStack.pop_back();
|
||||
|
||||
nCellsByConnnectedComponent[nConnnectedComponentIndex].append(nCurrCell);
|
||||
|
||||
// We push all face-neighbouring cells that:
|
||||
// 1) Are in the same layer
|
||||
// 2) Haven't been found yet
|
||||
for (label nFaceIdx : mesh.cells()[nCurrCell])
|
||||
{
|
||||
label nNeiCell = getNei(nCurrCell, nFaceIdx);
|
||||
if (nNeiCell < 0) continue;
|
||||
if (nCellLayerIndex[nNeiCell] != nCellLayerIndex[nCurrCell]) continue;
|
||||
if (nCellConnnectedComponentIndex[nNeiCell] != -1) continue;
|
||||
|
||||
nCellConnnectedComponentIndex[nNeiCell] = nConnnectedComponentIndex;
|
||||
nDfsStack.push_back(nNeiCell);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::hpathRenumber::hpathFinder::solveConnectedComponent
|
||||
(
|
||||
const Foam::labelList& nCellsInConnectedComponent,
|
||||
Foam::labelList& cellOrder
|
||||
)
|
||||
{
|
||||
// For small cases, find path manually (1-2 cells)
|
||||
if (nCellsInConnectedComponent.size() <= 2) {
|
||||
for (label nCellIdx : nCellsInConnectedComponent) {
|
||||
cellOrder[nFoundCellCount++] = nCellIdx;
|
||||
bIsRenumbered[nCellIdx] = true;
|
||||
}
|
||||
return;
|
||||
}
|
||||
// Get a graph representing the connected component
|
||||
getConnnectedComponentGraph(nCellsInConnectedComponent);
|
||||
|
||||
// Get the starting cell
|
||||
label nStartCell = getStartingCellInConnnectedComponent(nCellsInConnectedComponent);
|
||||
|
||||
// Reorder each cells neighbours in the ConnnectedComponent-Graph in descending order by distance from start cell
|
||||
reorderDistFromStart(nStartCell, nCellsInConnectedComponent);
|
||||
|
||||
// Get a path through the connected component, using the ConnnectedComponent-Graph
|
||||
findPath(nStartCell, cellOrder);
|
||||
}
|
||||
|
||||
|
||||
void Foam::hpathRenumber::hpathFinder::getConnnectedComponentGraph
|
||||
(
|
||||
const Foam::labelList& nCellsInConnectedComponent
|
||||
)
|
||||
{
|
||||
// This method computes the ConnnectedComponent-Graph
|
||||
// Cells are connected in the ConnnectedComponent-Graph if:
|
||||
// a) They are in the same layer
|
||||
// b) They are face-connected in the mesh (different from Mesh-Graph - there it was point-connected)
|
||||
// - note that this implies that they are also in the same connected component (hence the name)
|
||||
|
||||
for (label nCellIdx : nCellsInConnectedComponent)
|
||||
{
|
||||
Foam::DynamicList<label> nNeiList;
|
||||
|
||||
// For every face on the cell, find its neighbour through that face (if there is one)
|
||||
// If that neighbour:
|
||||
// 1) Hasn't been added to the reordering yet
|
||||
// 2) Is in the same layer
|
||||
// 3) Is in the same connected component
|
||||
// Then we add it as a neighbour in the ConnectedComponent-Graph
|
||||
for (label nFaceIdx : mesh.cells()[nCellIdx])
|
||||
{
|
||||
label nNeiCell = getNei(nCellIdx, nFaceIdx);
|
||||
if (nNeiCell < 0) continue;
|
||||
if (bIsRenumbered[nNeiCell]) continue;
|
||||
if (nCellLayerIndex[nNeiCell] != nCellLayerIndex[nCellIdx]) continue;
|
||||
|
||||
nNeiList.append(nNeiCell);
|
||||
}
|
||||
|
||||
nConnnectedComponentGraph[nCellIdx] = nNeiList;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Foam::label Foam::hpathRenumber::hpathFinder::getStartingCellInConnnectedComponent
|
||||
(
|
||||
const Foam::labelList& nCellsInConnectedComponent
|
||||
)
|
||||
{
|
||||
|
||||
// Strategy for choosing a starting cell in the connected component:
|
||||
// 1) Start from an arbitrary cell within the connected component
|
||||
// 2) Find the furthest cell from it: choose it as the starting cell
|
||||
|
||||
// The reasoning behind this strategy is as follows:
|
||||
// - In many cases the cells given may be of the general shape of a long straight path
|
||||
// - In these cases, by starting from the farthest cell we guarantee it will be at one of the edges of the path
|
||||
|
||||
// It does not matter which cell we start from, so arbitrarly start from first cell in the connected component
|
||||
label nCellIdx = nCellsInConnectedComponent[0];
|
||||
|
||||
// Find farthest cell from nCellIdx
|
||||
// This is done using a standard BFS algorithm
|
||||
Foam::CircularBuffer<label> nBfsQueue;
|
||||
|
||||
bBFSFoundCell[nCellIdx] = true;
|
||||
nBfsQueue.push_back(nCellIdx);
|
||||
|
||||
label nCurrCell = -1;
|
||||
|
||||
while(!nBfsQueue.empty())
|
||||
{
|
||||
nCurrCell = nBfsQueue.first();
|
||||
nBfsQueue.pop_front();
|
||||
|
||||
// Push all cells that have not yet been found by the BFS
|
||||
for (label nNeiCell : nConnnectedComponentGraph[nCurrCell])
|
||||
{
|
||||
if (bBFSFoundCell[nNeiCell]) continue;
|
||||
bBFSFoundCell[nNeiCell] = true;
|
||||
nBfsQueue.push_back(nNeiCell);
|
||||
}
|
||||
}
|
||||
|
||||
// nCurrCell should now be the farthest cell from the starting cell in the connected component
|
||||
return nCurrCell;
|
||||
}
|
||||
|
||||
|
||||
void Foam::hpathRenumber::hpathFinder::reorderDistFromStart
|
||||
(
|
||||
label nStartCell,
|
||||
const Foam::labelList& nCellsInConnectedComponent
|
||||
)
|
||||
{
|
||||
// findPath() tends to find much better results when each cell's neighbours are ordered in a specific way based on their face/point-distance to the starting cell
|
||||
|
||||
// First, this method finds for each cell in the connected component its face-distance from the start cell
|
||||
// Secondly, this method finds for each cell in the connected component its point-distance from the start cell
|
||||
// Finally, this method reorders each cells neighbours in the ConnnectedComponent-Graph accordingly
|
||||
|
||||
// - Both of these are done using a BFS algorithm: the first using face-neigbours and the second using point-neighbours
|
||||
|
||||
Foam::CircularBuffer<label> nBfsQueue;
|
||||
|
||||
// Find the face-distance of all cells in the connected component from the starting cell
|
||||
nCellFaceDistFromStart[nStartCell] = 0;
|
||||
// Push starting cell to queue
|
||||
nBfsQueue.push_back(nStartCell);
|
||||
|
||||
while (!nBfsQueue.empty())
|
||||
{
|
||||
label nCurrCell = nBfsQueue.first();
|
||||
nBfsQueue.pop_front();
|
||||
|
||||
// Push all face-neighbours that haven't already had their face distance found
|
||||
// Face-neighbours are saved in the ConnnectedComponent-Graph
|
||||
for (label nNeiCell : nConnnectedComponentGraph[nCurrCell])
|
||||
{
|
||||
if (nCellFaceDistFromStart[nNeiCell] != -1) continue;
|
||||
|
||||
nCellFaceDistFromStart[nNeiCell] = nCellFaceDistFromStart[nCurrCell] + 1;
|
||||
nBfsQueue.push_back(nNeiCell);
|
||||
}
|
||||
}
|
||||
// Find the *point*-distance of all cells in the connected component from the starting cell
|
||||
nCellPointDistFromStart[nStartCell] = 0;
|
||||
// Push starting cell to queue
|
||||
nBfsQueue.push_back(nStartCell);
|
||||
|
||||
while (!nBfsQueue.empty())
|
||||
{
|
||||
label nCurrCell = nBfsQueue.first();
|
||||
nBfsQueue.pop_front();
|
||||
|
||||
// Push all the point-neighbour that:
|
||||
// 1) Haven't already been pushed by the BFS previously
|
||||
// 2) Are in the same layer as the Starting Cell
|
||||
// 3) Are in the same connected component as the Starting Cell
|
||||
// Point neighbours are saved in the Mesh-Graph
|
||||
for (label nNeiCell : nMeshGraph[nCurrCell])
|
||||
{
|
||||
if (bIsRenumbered[nNeiCell]) continue;
|
||||
if (nCellPointDistFromStart[nNeiCell] != -1) continue;
|
||||
if (nCellLayerIndex[nNeiCell] != nCellLayerIndex[nStartCell]) continue;
|
||||
if (nCellConnnectedComponentIndex[nNeiCell] != nCellConnnectedComponentIndex[nStartCell]) continue;
|
||||
|
||||
nCellPointDistFromStart[nNeiCell] = nCellPointDistFromStart[nCurrCell] + 1;
|
||||
nBfsQueue.push_back(nNeiCell);
|
||||
}
|
||||
}
|
||||
|
||||
for (label nCellIdx : nCellsInConnectedComponent)
|
||||
{
|
||||
// Sorting of neighbours is done in two levels:
|
||||
// 1) Neigbours are sorted *descending*-order based on their *point*-distance from the starting cell
|
||||
// 2) Neigbours with the same *point*-distance are sorted in *ascending*-order based on their *face*-distance from the starting cell
|
||||
std::sort
|
||||
(
|
||||
nConnnectedComponentGraph[nCellIdx].begin(),
|
||||
nConnnectedComponentGraph[nCellIdx].end(),
|
||||
[this](label i, label j) {
|
||||
if (nCellPointDistFromStart[i] != nCellPointDistFromStart[j])
|
||||
return nCellPointDistFromStart[i] > nCellPointDistFromStart[j];
|
||||
else
|
||||
return nCellFaceDistFromStart[i] < nCellFaceDistFromStart[j];
|
||||
}
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::hpathRenumber::hpathFinder::findPath
|
||||
(
|
||||
label nStartCell, Foam::labelList& cellOrder
|
||||
)
|
||||
{
|
||||
// Tries to find the furthest cell from the starting cell within the ConnectedComponent-Graph
|
||||
// This is done using a DFS algorithm:
|
||||
// - Run DFS from starting cell
|
||||
// - Return path to deepest cell found by DFS
|
||||
|
||||
// Run a standard DFS algorithm from starting cell
|
||||
Foam::CircularBuffer<label> nDfsStack;
|
||||
nDfsStack.push_back(nStartCell);
|
||||
nDFSParentCell[nStartCell] = nStartCell;
|
||||
|
||||
// Used for finding and returning the best path
|
||||
label nBestCell = nStartCell;
|
||||
|
||||
while (!nDfsStack.empty())
|
||||
{
|
||||
label nCurrCell = nDfsStack.last();
|
||||
nDfsStack.pop_back();
|
||||
|
||||
// If the cell has already been popped previously, we can skip it
|
||||
if (nDFSCellDepth[nCurrCell] != -1) continue;
|
||||
|
||||
// Otherwise, we update its depth value. This means the cell will never be PUSHED again
|
||||
nDFSCellDepth[nCurrCell] = nDFSCellDepth[nDFSParentCell[nCurrCell]] + 1;
|
||||
|
||||
// If this is the new deepest cell, update best
|
||||
if (nDFSCellDepth[nCurrCell] > nDFSCellDepth[nBestCell]) {
|
||||
nBestCell = nCurrCell;
|
||||
}
|
||||
|
||||
// Cells will be popped from the stack in reverse order from how we pushed them
|
||||
// By iterating over the neighbors in reverse order, cells will be popped in the correct order
|
||||
for (label nNeiIdx = nConnnectedComponentGraph[nCurrCell].size() - 1; nNeiIdx >= 0; nNeiIdx--)
|
||||
{
|
||||
label nNeiCell = nConnnectedComponentGraph[nCurrCell][nNeiIdx];
|
||||
if (nDFSCellDepth[nNeiCell] == -1)
|
||||
{
|
||||
// Notice we only update the nCellDepth value of a cell when we POP it from the stack
|
||||
// This means some cells may be pushed many times, but they will only be popped once
|
||||
nDFSParentCell[nNeiCell] = nCurrCell;
|
||||
nDfsStack.push_back(nNeiCell);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
label nCellIdx = nBestCell;
|
||||
label nPrevCell = -1;
|
||||
while (nCellIdx != nPrevCell)
|
||||
{
|
||||
cellOrder[nFoundCellCount++] = nCellIdx;
|
||||
bIsRenumbered[nCellIdx] = true;
|
||||
|
||||
// The first cell in the path is always its own parent
|
||||
nPrevCell = nCellIdx;
|
||||
nCellIdx = nDFSParentCell[nCellIdx];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::hpathRenumber::hpathFinder::resetCells
|
||||
(
|
||||
Foam::labelList& nCellList
|
||||
)
|
||||
{
|
||||
// Reset the data structures for the cells that were not found in the renumbering
|
||||
for (label nCellIdx : nCellList)
|
||||
{
|
||||
nCellConnnectedComponentIndex[nCellIdx] = -1;
|
||||
nCellPointDistFromStart[nCellIdx] = -1;
|
||||
nCellFaceDistFromStart[nCellIdx] = -1;
|
||||
bBFSFoundCell[nCellIdx] = false;
|
||||
nDFSCellDepth[nCellIdx] = -1;
|
||||
nDFSParentCell[nCellIdx] = -1;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
142
src/renumber/renumberMethods/hpathRenumber/hpathRenumber.H
Normal file
142
src/renumber/renumberMethods/hpathRenumber/hpathRenumber.H
Normal file
@ -0,0 +1,142 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2023 Alon Zameret, Noam Manaker Morag
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
||||
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
||||
for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Class
|
||||
Foam::hpathRenumber
|
||||
|
||||
Description
|
||||
Renumbering based on Hamiltonian path.
|
||||
|
||||
SourceFiles
|
||||
hpathRenumber.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef Foam_hpathRenumber_H
|
||||
#define Foam_hpathRenumber_H
|
||||
|
||||
#include "renumberMethod.H"
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class hpathRenumber Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class hpathRenumber
|
||||
:
|
||||
public renumberMethod
|
||||
{
|
||||
// Private Data
|
||||
|
||||
//- Apply layer separation? (default: true)
|
||||
bool applyLayerSeparation_;
|
||||
|
||||
|
||||
// Private Classes
|
||||
|
||||
// Helper class
|
||||
class hpathFinder;
|
||||
|
||||
public:
|
||||
|
||||
// Generated Methods
|
||||
|
||||
//- No copy construct
|
||||
hpathRenumber(const hpathRenumber&) = delete;
|
||||
|
||||
//- No copy assignment
|
||||
void operator=(const hpathRenumber&) = delete;
|
||||
|
||||
|
||||
//- Runtime type information
|
||||
TypeName("hpath");
|
||||
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct given the renumber dictionary
|
||||
explicit hpathRenumber(const dictionary& dict);
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~hpathRenumber() = default;
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Return the order in which cells need to be visited
|
||||
//- (ie. from ordered back to original cell label).
|
||||
virtual labelList renumber(const pointField&) const
|
||||
{
|
||||
NotImplemented;
|
||||
return labelList();
|
||||
}
|
||||
|
||||
//- Return the order in which cells need to be visited
|
||||
//- (ie. from ordered back to original cell label).
|
||||
// Uses mesh for regIOobject
|
||||
virtual labelList renumber
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const pointField& cc
|
||||
) const;
|
||||
|
||||
//- Return the order in which cells need to be visited
|
||||
//- (ie. from ordered back to original cell label).
|
||||
virtual labelList renumber
|
||||
(
|
||||
const CompactListList<label>& cellCells,
|
||||
const pointField& cellCentres
|
||||
) const
|
||||
{
|
||||
NotImplemented;
|
||||
return labelList();
|
||||
}
|
||||
|
||||
//- Return the order in which cells need to be visited
|
||||
//- (ie. from ordered back to original cell label).
|
||||
virtual labelList renumber
|
||||
(
|
||||
const labelListList& cellCells,
|
||||
const pointField& cellCentres
|
||||
) const
|
||||
{
|
||||
NotImplemented;
|
||||
return labelList();
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -119,8 +119,8 @@ bool Foam::viscosityModels::Casson::read
|
||||
|
||||
CassonCoeffs_.readEntry("m", m_);
|
||||
CassonCoeffs_.readEntry("tau0", tau0_);
|
||||
CassonCoeffs_.readEntry("nuMin", nuMin_);
|
||||
CassonCoeffs_.readEntry("nuMax", nuMax_);
|
||||
CassonCoeffs_.readEntry("nuMin_", nuMin_);
|
||||
CassonCoeffs_.readEntry("nuMax_", nuMax_);
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
@ -1,26 +1,21 @@
|
||||
# ----------------------------------------------------------------------------
|
||||
# CGAL definitions - several possibilities
|
||||
#
|
||||
# - missing
|
||||
# - header-only
|
||||
# - header-only, no mpfr
|
||||
# - library, no mpfr
|
||||
# - library, with mpfr (default for older CGAL)
|
||||
# 0. missing
|
||||
# 1. header-only
|
||||
# 2. library, no mpfr
|
||||
# 3. library, with mpfr (a likely default)
|
||||
#
|
||||
# Dispatch according to the defined 'CGAL_FLAVOUR'
|
||||
# - names may change [see wmake/scripts/have_cgal]
|
||||
# (no-cgal | cgal-header | cgal-header-no-mpfr | cgal-no-mpfr | cgal-mpfr)
|
||||
# (no-cgal | cgal-header | cgal-no-mpfr | cgal-mpfr)
|
||||
|
||||
cgal_subrule := cgal-mpfr
|
||||
ifneq (,$(findstring header,$(CGAL_FLAVOUR)))
|
||||
cgal_subrule := cgal-header-only
|
||||
endif
|
||||
ifneq (,$(findstring no-mpfr,$(CGAL_FLAVOUR)))
|
||||
cgal_subrule := cgal-no-mpfr
|
||||
ifneq (,$(findstring header,$(CGAL_FLAVOUR)))
|
||||
cgal_subrule := cgal-header-no-mpfr
|
||||
endif
|
||||
else
|
||||
ifneq (,$(findstring header,$(CGAL_FLAVOUR)))
|
||||
cgal_subrule := cgal-header-only
|
||||
endif
|
||||
endif
|
||||
|
||||
# ----------------------------------------------------------------------------
|
||||
|
||||
@ -1,18 +0,0 @@
|
||||
# -----------------------------------------------------------------------------
|
||||
# CGAL (header-only version) without mpfr
|
||||
|
||||
CGAL_INC = -DCGAL_HEADER_ONLY
|
||||
CGAL_LIBS =
|
||||
|
||||
CGAL_INC += \
|
||||
$(foreach dir,$(BOOST_INC_DIR),-I$(dir)) \
|
||||
$(foreach dir,$(CGAL_INC_DIR),-I$(dir))
|
||||
|
||||
CGAL_LIBS += \
|
||||
$(foreach dir,$(BOOST_LIB_DIR),-L$(dir))
|
||||
|
||||
# ----
|
||||
# Extra failsafe - still needed? (2020-05-15)
|
||||
## CGAL_INC += -I/usr/local/include -I/usr/include
|
||||
|
||||
# -----------------------------------------------------------------------------
|
||||
14
wmake/rules/darwin64Clang/cgal
Normal file
14
wmake/rules/darwin64Clang/cgal
Normal file
@ -0,0 +1,14 @@
|
||||
# ----------------------------------------------------------------------------
|
||||
# CGAL on Darwin
|
||||
# CGAL (library version) without mpfr
|
||||
|
||||
CGAL_INC = \
|
||||
$(foreach dir,$(BOOST_INC_DIR),-I$(dir)) \
|
||||
$(foreach dir,$(CGAL_INC_DIR),-I$(dir))
|
||||
|
||||
CGAL_LIBS = \
|
||||
$(foreach dir,$(BOOST_LIB_DIR),-L$(dir)) \
|
||||
$(foreach dir,$(CGAL_LIB_DIR),-L$(dir)) \
|
||||
-lCGAL
|
||||
|
||||
# ----------------------------------------------------------------------------
|
||||
@ -35,7 +35,6 @@ LINKLIBSO = $(CC) $(c++FLAGS) -shared \
|
||||
-Wl,--strip-all
|
||||
|
||||
LINKEXE = $(CC) $(c++FLAGS) \
|
||||
-static-libgcc -static-libstdc++ \
|
||||
-Wl,--enable-auto-import \
|
||||
-Wl,--strip-all \
|
||||
-Wl,--force-exe-suffix
|
||||
|
||||
Reference in New Issue
Block a user