songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Mapping injected faces through interpolation

Browse Source
This commit is contained in:
Mattijs Janssens 2018-11-15 14:58:56 +00:00
parent edf5cf4588
commit 8ceec6456a
14 changed files with 1906 additions and 127 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
README.txt Normal file
View File

@ -0,0 +1,47 @@
Notes from merging Integration-TUD-addOns and plus.
0)
- Changed
dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
to created internal faces out-of-nothing.
- Added the mapNewInternalFaces to dynamicRefineFvMesh
- unset FOAM_SETNAN, FOAM_SIGFPE and run
testCases/testAMRandLoadBalancing/damBreakWithObstacle
- However when writing V0 we noticed a nan in the V0 field.
(this doesn't get written anymore but it might indicate a bug)
1) Tried moving additional mapping of surface fields to
virtual dynamicRefineFvMesh::mapFields(const mapPolyMesh&);
so it would get called from the fvMesh::updateMesh. However
this mapping (explicitly) gets done using unadapted addressing
and it causes the cell addressing to be wrong:
[1] --> FOAM FATAL ERROR:
[1] index 826015320 out of range 0 ... 71079
[1]
[1] From function void Foam::UList<T>::checkIndex(Foam::label) const [with T = Foam::cell; Foam::label = int]
[1] in file /home/preston2/mattijs/OpenFOAM/work/OpenFOAM-plus.integration-TUD/src/OpenFOAM/lnInclude/UListI.H at line 106.
From dynamicRefineFvMesh::mapNewInternalFaces : (I think)
cell faceOwner = this->cells()[owner[facei]];
So this adaptation should be done in a separate pass
2) Mapping new internal faces:
- currently done by averaging the values of 'properly' mapped
faces from owner and neighbour.
- this would not work if a cell gets split into 3x3x3
- instead this should be done by some geometric interpolation
from point values?
- have selection mechanism instead or use oriented flag on surfaceFields
3) Current averaging is valid only for internal faces & only scalar/vector

1
src/dynamicFvMesh/Make/files
View File

@ -5,6 +5,7 @@ dynamicMotionSolverFvMesh/dynamicMotionSolverFvMesh.C
dynamicMultiMotionSolverFvMesh/dynamicMultiMotionSolverFvMesh.C
dynamicInkJetFvMesh/dynamicInkJetFvMesh.C
dynamicRefineFvMesh/dynamicRefineFvMesh.C
dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.C
dynamicMotionSolverListFvMesh/dynamicMotionSolverListFvMesh.C
simplifiedDynamicFvMesh/simplifiedDynamicFvMeshes.C

3
src/dynamicFvMesh/Make/options
View File

@ -1,10 +1,13 @@
EXE_INC = \
-g -DFULLDEBUG -O0 -g \
-I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
LIB_LIBS = \
-ldecompositionMethods \
-lmeshTools \
-ldynamicMesh \
-lfiniteVolume

29
src/dynamicFvMesh/dynamicRefineBalancedFvMesh/balanceParDict Normal file
View File

@ -0,0 +1,29 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object balanceParDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
numberOfSubdomains 2;
method ptscotch;
constraints
{
refinementHistory
{
type refinementHistory;
}
}
// ************************************************************************* //

125
src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicMeshDict Normal file
View File

@ -0,0 +1,125 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object dynamicMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dynamicFvMesh dynamicRefineBalancedFvMesh;
refinementControls
{
enableRefinementControl true;
fields // must be scalarFields
(
//alpha (min max refineLevel)
alpha (0.01 0.99 2) // refine cells where alpha in [0.01:0.99] with maximal 2 refinement layers
);
interface // must be a scalarField (only one dictionary!)
(
alpha // refine interface (found based on snGrad of alpha > 0.1)
{
innerRefLayers 2; // describes how many cell layers inside phase alpha are to be refined
outerRefLayers 5; // describes how many cell layers outside phase alpha are to be refined
// optional settings:
maxRefineLevel 4; // max refinement layers; Default: maxRefinement from dynamicRefineFvMeshCoeffs is used
// to get slower than 2:1 refinement; add #nAddLayers between each refinement level at the interface
nAddLayers 1; //Default: 0
}
);
gradients // must be scalars
(
// arguments as in 'fields'
// min/max values are based on mag(fvc::grad(volScalarField)) * cellVolume
T (0.01 10 1)
);
curls // must be vectors
(
// arguments as in 'fields'
// min/max values are based on mag(fvc::curl(volVectorField))
U (0.5 1 2)
);
regions
(
boxToCell
{
minLevel 1;
box (-1 0.001 0.002)(1 0.005 0.003);
}
);
}
dynamicRefineFvMeshCoeffs
{
// Extra entries for balancing
enableBalancing true;
allowableImbalance 0.15;
// How often to refine
refineInterval 10;
// Field to be refinement on (set it to 'internalRefinementField' to use the
// refinementControls dictionary entries above)
field internalRefinementField;
// Refine field inbetween lower..upper
lowerRefineLevel 0.5; // do not change
upperRefineLevel 3.5; // maxRefinement+0.5
// If value < unrefineLevel unrefine
unrefineLevel -0.5; // do not change
// Have slower than 2:1 refinement
nBufferLayers 4;
// Refine cells only up to maxRefinement levels
maxRefinement 3;
// Stop refinement if maxCells reached
maxCells 200000;
// Flux field and corresponding velocity field. Fluxes on changed
// faces get recalculated by interpolating the velocity. Use 'none'
// on surfaceScalarFields that do not need to be reinterpolated.
correctFluxes
(
(phi Urel)
(phiAbs U)
(phiAbs_0 U_0)
(nHatf none)
(rho*phi none)
(ghf none)
);
// List of non-flux surface<Type>Fields to be mapped
// only for new internal faces (AMR refine)
mapSurfaceFields
(
Uf
Uf_0
);
// Write the refinement level as a volScalarField
dumpLevel true;
}
// ************************************************************************* //

843
src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.C Normal file
View File

@ -0,0 +1,843 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2014 Tyler Voskuilen
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is a derivative work 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 "dynamicRefineBalancedFvMesh.H"
#include "addToRunTimeSelectionTable.H"
#include "surfaceInterpolate.H"
#include "volFields.H"
#include "polyTopoChange.H"
#include "surfaceFields.H"
#include "syncTools.H"
#include "pointFields.H"
#include "fvCFD.H"
#include "volPointInterpolation.H"
#include "pointMesh.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(dynamicRefineBalancedFvMesh, 0);
addToRunTimeSelectionTable(dynamicFvMesh, dynamicRefineBalancedFvMesh, IOobject);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::label Foam::dynamicRefineBalancedFvMesh::topParentID(label p)
{
label nextP = meshCutter().history().splitCells()[p].parent_;
if( nextP < 0 )
{
return p;
}
else
{
return topParentID(nextP);
}
}
Foam::List<Foam::scalar> Foam::dynamicRefineBalancedFvMesh::readRefinementPoints()
{
dictionary refineDict
(
IOdictionary
(
IOobject
(
"dynamicMeshDict",
time().constant(),
*this,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
).subDict("dynamicRefineFvMeshCoeffs")
);
List<scalar> refData(4, scalar(0));
refData[0] = readScalar(refineDict.lookup("unrefineLevel"));
refData[1] = readScalar(refineDict.lookup("lowerRefineLevel"));
refData[2] = readScalar(refineDict.lookup("refineInterval"));
refData[3] = readScalar(refineDict.lookup("maxRefinement"));
return refData;
}
void Foam::dynamicRefineBalancedFvMesh::updateRefinementField()
{
Info<< "Calculating internal refinement field" << endl;
volScalarField& intRefFld = *internalRefinementFieldPtr_;
volScalarField& targetFld = *targetLevelPtr_;
volScalarField& currFld = *isLevelPtr_;
// Set the internal refinement field to zero to start with
intRefFld = dimensionedScalar("zero",dimless,0.0);
// Get the cell level field from dynamicRefineFvMesh
const labelList& cellLevel = meshCutter().cellLevel();
// Read the points at which refinement and unrefinement occur from the
// dynamicMeshDict entries
List<scalar> refinePoints = readRefinementPoints();
// init list for refine/unrefine field (-1=unrefine, 1=refine, 0=do nothing)
labelList markRefineCells (this->nCells(), 0);
// init list for target refinement level per cell
labelList targetLevel (this->nCells(), 0);
// First fields
List<word> fieldNames = fields_.toc();
Field<scalar> refFld(nCells(),0.0);
forAll(fieldNames, i)
{
word fldName = fieldNames[i];
scalar minValue = fields_[fldName][0];
scalar maxValue = fields_[fldName][1];
label refineLevel = static_cast<label>(fields_[fldName][2]);
const volScalarField& fld = this->lookupObject<volScalarField>(fldName);
// Limit the value of refFld based on its max level
forAll(fld, cellI)
{
if ((fld[cellI] >= minValue) && (fld[cellI] <= maxValue))
{
// increase targetLevel up to refineLevel
// BUT: do not decrease if cell already marked for higher refinement level by previous criterion
targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
}
}
}
// Then gradients
List<word> gradFieldNames = gradFields_.toc();
Field<scalar> cubeRtV = Foam::pow(this->V(),1.0/3.0);
forAll(gradFieldNames, i)
{
word fldName = gradFieldNames[i];
scalar minValue = gradFields_[fldName][0];
scalar maxValue = gradFields_[fldName][1];
label refineLevel = static_cast<label>(gradFields_[fldName][2]);
const volScalarField& fld = this->lookupObject<volScalarField>(fldName);
refFld = mag(fvc::grad(fld)) * cubeRtV;
// Limit the value of refFld based on its max level
forAll(refFld, cellI)
{
if ((refFld[cellI] >= minValue) && (refFld[cellI] <= maxValue))
{
targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
}
}
}
// Then curls
List<word> curlFieldNames = curlFields_.toc();
forAll(curlFieldNames, i)
{
word fldName = curlFieldNames[i];
scalar minValue = curlFields_[fldName][0];
scalar maxValue = curlFields_[fldName][1];
label refineLevel = static_cast<label>(curlFields_[fldName][2]);
const volVectorField& fld = this->lookupObject<volVectorField>(fldName);
refFld = mag(fvc::curl(fld));
// Limit the value of refFld based on its max level
forAll(refFld, cellI)
{
if ((refFld[cellI] >= minValue) && (refFld[cellI] <= maxValue))
{
targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
}
}
}
// At the interface, assumed to be always the maximum refinement level
List<word> interfaceRefineField = interface_.toc();
forAll(interfaceRefineField, i)
{
word fldName = interfaceRefineField[i];
// read region of maximum refinement levels inside and outside of interface indicator field
// (does not need to be alpha, can also be a concentration field)
scalar innerRefLayers = readScalar(interface_[fldName].lookup("innerRefLayers"));
scalar outerRefLayers = readScalar(interface_[fldName].lookup("outerRefLayers"));
scalar nAddLayers(0);
if (interface_[fldName].found("nAddLayers"))
{
nAddLayers = readScalar(interface_[fldName].lookup("nAddLayers"));
}
label refineLevel = refinePoints[3];
if (interface_[fldName].found("maxRefineLevel"))
{
refineLevel = readScalar(interface_[fldName].lookup("maxRefineLevel"));
// to avoid user input mistakes, limit the value with the maximum allowed
refineLevel = min(refinePoints[3], refineLevel);
}
const volScalarField& fld = this->lookupObject<volScalarField>(fldName);
volScalarField isInterface(intRefFld * 0.0);
isInterface = dimensionedScalar("isInterface",dimless,0.0);
surfaceScalarField deltaAlpha = mag(fvc::snGrad(fld) / this->deltaCoeffs());
const labelUList& owner = this->owner();
const labelUList& neighbour = this->neighbour();
forAll(deltaAlpha, faceI)
{
//TODO: add a min max value of the specified field to be marked for refinement
// NOW: hard-coded method snGradAlpha: checks mag(alpha_N - alpha_P) > 0.1 ? 1:0
if (deltaAlpha[faceI] > 0.1) // currently a fixed prescribed value; should be read from díctionary
{
label own = owner[faceI];
label nei = neighbour[faceI];
// set isInterface field to one
isInterface[own] = 1.0;
isInterface[nei] = 1.0;
}
}
// assumed fld=0.5*(fldMax+fldMin) defines the interface
dimensionedScalar fldInterfaceValue(0.5*(max(fld)+min(fld)));
//-DD: old implementation based on face interpolation
// which results in slower transport in diagonal direction
// add inner refinement layers
//for(label i=0; i < innerRefLayers; i++)
//{
// isInterface += neg(- fvc::average(fvc::interpolate(isInterface)) * pos(fld - fldInterfaceValue));
// isInterface = neg(- isInterface);
//}
//
// add outer refinement layers
//for(label i=0; i < outerRefLayers; i++)
//{
// isInterface += neg(- fvc::average(fvc::interpolate(isInterface)) * pos(fldInterfaceValue - fld));
// isInterface = neg(- isInterface);
//}
//
//forAll(isInterface, cellI)
//{
// if (isInterface[cellI] > 0.5)
// {
// targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
// }
//}
//-DD: new version using volPointInterpolation (direction independent buffer layer)
const volPointInterpolation& pInterp = volPointInterpolation::New(*this);
const fvMesh& mesh = fld.mesh();
// add inner refinement layers
for(label i=0; i < innerRefLayers; i++)
{
volScalarField markInner(isInterface*pos(fld - fldInterfaceValue));
pointScalarField markLayerP(pInterp.interpolate(markInner));
forAll(mesh.C(), cellI)
{
scalar sum = 0.;
label nPoints = 0;
forAll(mesh.cellPoints()[cellI], pointI)
{
sum += markLayerP[mesh.cellPoints()[cellI][pointI]];
nPoints++;
}
if (nPoints > 0)
{
sum /= nPoints;
}
isInterface[cellI] += sum;
}
}
isInterface = pos(isInterface - SMALL);
// add outer refinement layers
for(label i=0; i < outerRefLayers; i++)
{
volScalarField markOuter(isInterface*pos(fldInterfaceValue - fld));
pointScalarField markLayerP(pInterp.interpolate(markOuter));
forAll(mesh.C(), cellI)
{
scalar sum = 0.;
label nPoints = 0;
forAll(mesh.cellPoints()[cellI], pointI)
{
sum += markLayerP[mesh.cellPoints()[cellI][pointI]];
nPoints++;
}
if (nPoints > 0)
{
sum /= nPoints;
}
isInterface[cellI] += sum;
}
}
isInterface = pos(isInterface - SMALL);
forAll(isInterface, cellI)
{
if (isInterface[cellI] > 0.5)
{
targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
}
}
//-DD: old implementation based on face interpolation
// which results in slower transport in diagonal direction
// expand additional layers if specified:
if (nAddLayers > 0)
{
// loop over additional layers
for(label i=1; i < refineLevel; i++)
{
// #nAddLayers cells per additional level
for(label j=0; j < nAddLayers; j++)
{
isInterface += neg(- fvc::average(fvc::interpolate(isInterface)));
isInterface = neg(- isInterface);
}
forAll(isInterface, cellI)
{
if (isInterface[cellI] == 1.0)
{
targetLevel[cellI] = max(targetLevel[cellI], (refineLevel - i));
}
}
}
}
}
// The set refinement physical regions (force the mesh to stay refined
// near key features)
forAll(refinedRegions_, regionI)
{
const entry& region = refinedRegions_[regionI];
autoPtr<topoSetSource> source =
topoSetSource::New(region.keyword(), *this, region.dict());
cellSet selectedCellSet
(
*this,
"cellSet",
nCells()/10+1 //estimate
);
source->applyToSet
(
topoSetSource::NEW,
selectedCellSet
);
const labelList cells = selectedCellSet.toc();
label minLevel = readLabel(region.dict().lookup("minLevel"));
forAll(cells, i)
{
const label& cellI = cells[i];
targetLevel[cellI] = max(targetLevel[cellI], minLevel);
}
}
//-DD: add buffer layer based on targetLevel field to prevent 2-to-1 refinement
volScalarField blockedLevel = targetFld * 0.;
for (label currLayer=refinePoints[3]; currLayer>=1; currLayer--)
{
forAll (targetLevel, cellI)
{
if (targetLevel[cellI] >= currLayer)
{
blockedLevel[cellI] = targetLevel[cellI];
}
}
for (label i=0; i<nBufferLayers_; i++)
{
blockedLevel = max(blockedLevel, pos(fvc::average(fvc::interpolate(blockedLevel)) - SMALL)*(currLayer-1));
}
labelList blockLev(blockedLevel.internalField().size(), 0);
forAll (blockLev, cellI)
{
blockLev[cellI] = blockedLevel[cellI];
}
targetLevel = max(targetLevel, blockLev);
}
// mark cells to be refined/unrefined based on above criteria:
// simply check, if targetLevel lower or higher cellLevel
forAll (intRefFld.internalField(), cellI)
{
intRefFld.primitiveFieldRef()[cellI] = targetLevel[cellI] - cellLevel[cellI];
targetFld.primitiveFieldRef()[cellI] = targetLevel[cellI];
currFld.primitiveFieldRef()[cellI] = cellLevel[cellI];
}
intRefFld.correctBoundaryConditions();
Info<<"Min,max refinement field = " << Foam::min(intRefFld).value() << ", "
<< Foam::max(intRefFld).value() << endl;
}
void Foam::dynamicRefineBalancedFvMesh::readRefinementDict()
{
dictionary dynamicMeshDict
(
IOdictionary
(
IOobject
(
"dynamicMeshDict",
time().constant(),
*this,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
)
);
const dictionary refineDict
(
IOdictionary
(
IOobject
(
"dynamicMeshDict",
time().constant(),
*this,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
).subDict("dynamicRefineFvMeshCoeffs")
);
{
wordList surfFlds;
if (refineDict.readIfPresent("mapSurfaceFields", surfFlds))
{
// Rework into hashtable.
mapSurfaceFields_.resize(surfFlds.size());
forAll(surfFlds, i)
{
mapSurfaceFields_.insert(surfFlds[i], surfFlds[i]);
}
}
}
nBufferLayers_ = readLabel(refineDict.lookup("nBufferLayers"));
if( dynamicMeshDict.isDict("refinementControls") )
{
dictionary refineControlDict =
dynamicMeshDict.subDict("refinementControls");
enableRefinementControl_ =
Switch(refineControlDict.lookup("enableRefinementControl"));
if( enableRefinementControl_ )
{
// Overwrite field name entry in dynamicRefineFvMeshCoeffs?
// For now you just have to be smart and enter
// 'internalRefinementField' for the field name manually
// Read HashTable of field-refinement scalars
if( refineControlDict.found("fields") )
{
fields_ = HashTable< List<scalar> >
(
refineControlDict.lookup("fields")
);
}
// Read HashTable of gradient-refinement scalars
if( refineControlDict.found("gradients") )
{
gradFields_ = HashTable< List<scalar> >
(
refineControlDict.lookup("gradients")
);
}
// Read HashTable of curl-refinement vectors
if( refineControlDict.found("curls") )
{
curlFields_ = HashTable< List<scalar> >
(
refineControlDict.lookup("curls")
);
}
// Read interface refinement data
if( refineControlDict.found("interface") )
{
interface_ = HashTable< dictionary >
(
refineControlDict.lookup("interface")
);
}
// Read refinement regions
if( refineControlDict.found("regions") )
{
refinedRegions_ = PtrList<entry>
(
refineControlDict.lookup("regions")
);
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::dynamicRefineBalancedFvMesh::dynamicRefineBalancedFvMesh
(
const IOobject& io
)
:
dynamicRefineFvMesh(io),
internalRefinementFieldPtr_(NULL),
targetLevelPtr_(NULL),
isLevelPtr_(NULL),
fields_(),
gradFields_(),
curlFields_(),
interface_(),
refinedRegions_(),
enableRefinementControl_(false),
nBufferLayers_(0),
rebalance_(false)
{
readRefinementDict();
if( enableRefinementControl_ )
{
internalRefinementFieldPtr_ = new volScalarField
(
IOobject
(
"internalRefinementField",
time().timeName(),
*this,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
*this,
dimensionedScalar("zero", dimless, 0.0)
);
targetLevelPtr_ = new volScalarField
(
IOobject
(
"targetLevel",
time().timeName(),
*this,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
*this,
dimensionedScalar("zero", dimless, 0.0)
);
isLevelPtr_ = new volScalarField
(
IOobject
(
"isLevel",
time().timeName(),
*this,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
*this,
dimensionedScalar("zero", dimless, 0.0)
);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::dynamicRefineBalancedFvMesh::~dynamicRefineBalancedFvMesh()
{
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::dynamicRefineBalancedFvMesh::mapFields(const mapPolyMesh& mpm)
{
DebugVar(mpm.nOldCells());
dynamicRefineFvMesh::mapFields(mpm);
// Correct surface fields on introduced internal faces. These get
// created out-of-nothing so get an interpolated value.
mapNewInternalFaces<scalar>(mpm.faceMap());
mapNewInternalFaces<vector>(mpm.faceMap());
mapNewInternalFaces<sphericalTensor>(mpm.faceMap());
mapNewInternalFaces<symmTensor>(mpm.faceMap());
mapNewInternalFaces<tensor>(mpm.faceMap());
}
bool Foam::dynamicRefineBalancedFvMesh::update()
{
//Part 0 - Update internally calculated refinement field
readRefinementDict();
List<scalar> refinePoints = readRefinementPoints();
label refineInterval = refinePoints[2];
if( enableRefinementControl_ && time().timeIndex() > 0 && time().timeIndex() % refineInterval == 0 )
{
updateRefinementField();
}
//Part 1 - Call normal update from dynamicRefineFvMesh
bool hasChanged = dynamicRefineFvMesh::update();
if( Pstream::parRun() && hasChanged)
{
//Correct values on all coupled patches
correctBoundaries<scalar>();
correctBoundaries<vector>();
correctBoundaries<sphericalTensor>();
correctBoundaries<symmTensor>();
correctBoundaries<tensor>();
}
dictionary decomposeParDict
(
IOdictionary
(
IOobject
(
"decomposeParDict",
time().system(),
*this,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
)
);
rebalance_ = false;
// Part 2 - Load Balancing
{
dictionary refineDict
(
IOdictionary
(
IOobject
(
"dynamicMeshDict",
time().constant(),
*this,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
).subDict("dynamicRefineFvMeshCoeffs")
);
Switch enableBalancing = refineDict.lookup("enableBalancing");
if ( Pstream::parRun() && hasChanged )
{
const scalar allowableImbalance =
readScalar(refineDict.lookup("allowableImbalance"));
//First determine current level of imbalance - do this for all
// parallel runs with a changing mesh, even if balancing is disabled
label nGlobalCells = globalData().nTotalCells();
scalar idealNCells = scalar(nGlobalCells)/scalar(Pstream::nProcs());
scalar localImbalance = mag(scalar(nCells()) - idealNCells);
Foam::reduce(localImbalance, maxOp<scalar>());
scalar maxImbalance = localImbalance/idealNCells;
Info<<"Maximum imbalance = " << 100*maxImbalance << " %" << endl;
//If imbalanced, construct weighted coarse graph (level 0) with node
// weights equal to their number of subcells. This partitioning works
// as long as the number of level 0 cells is several times greater than
// the number of processors.
if( maxImbalance > allowableImbalance && enableBalancing)
{
Info << "\n**Solver hold for redistribution at time = " << time().timeName() << " s" << endl;
rebalance_ = true;
const labelIOList& cellLevel = meshCutter().cellLevel();
Map<label> coarseIDmap(100);
labelList uniqueIndex(nCells(),0);
label nCoarse = 0;
forAll(cells(), cellI)
{
if( cellLevel[cellI] > 0 )
{
uniqueIndex[cellI] = nCells() + topParentID
(
meshCutter().history().parentIndex(cellI)
);
}
else
{
uniqueIndex[cellI] = cellI;
}
if( coarseIDmap.insert(uniqueIndex[cellI], nCoarse) )
{
++nCoarse;
}
}
// Convert to local sequential indexing and calculate coarse
// points and weights
labelList localIndex(nCells(),0);
pointField coarsePoints(nCoarse,vector::zero);
scalarField coarseWeights(nCoarse,0.0);
forAll(uniqueIndex, cellI)
{
localIndex[cellI] = coarseIDmap[uniqueIndex[cellI]];
// If 2D refinement (quadtree) is ever implemented, this '3'
// should be set in general as the number of refinement
// dimensions.
label w = (1 << (3*cellLevel[cellI]));
coarseWeights[localIndex[cellI]] += 1.0;
coarsePoints[localIndex[cellI]] += C()[cellI]/w;
}
// Set up decomposer - a separate dictionary is used here so
// you can use a simple partitioning for decomposePar and
// ptscotch for the rebalancing (or any chosen algorithms)
autoPtr<decompositionMethod> decomposer
(
decompositionMethod::New
(
IOdictionary
(
IOobject
(
"balanceParDict",
time().system(),
*this,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
)
)
);
labelList finalDecomp = decomposer().decompose
(
*this,
localIndex,
coarsePoints,
coarseWeights
);
scalar tolDim = globalMeshData::matchTol_ * bounds().mag();
Info<< "Distributing the mesh ..." << endl;
fvMeshDistribute distributor(*this, tolDim);
Info<< "Mapping the fields ..." << endl;
autoPtr<mapDistributePolyMesh> map =
distributor.distribute(finalDecomp);
Info<< "Distribute the map ..." << endl;
meshCutter_.distribute(map);
Info << "Successfully distributed mesh" << endl;
scalarList procLoadNew (Pstream::nProcs(), 0.0);
procLoadNew[Pstream::myProcNo()] = this->nCells();
reduce(procLoadNew, sumOp<List<scalar> >());
scalar overallLoadNew = sum(procLoadNew);
scalar averageLoadNew = overallLoadNew/double(Pstream::nProcs());
Info << "New distribution: " << procLoadNew << endl;
Info << "Max deviation: " << max(Foam::mag(procLoadNew-averageLoadNew)/averageLoadNew)*100.0 << " %" << endl;
}
}
}
if( Pstream::parRun() && rebalance_)
{
//Correct values on all coupled patches
correctBoundaries<scalar>();
correctBoundaries<vector>();
correctBoundaries<sphericalTensor>();
correctBoundaries<symmTensor>();
correctBoundaries<tensor>();
}
return hasChanged;
}
// ************************************************************************* //

206
src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.H Normal file
View File

@ -0,0 +1,206 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2014 Tyler Voskuilen
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is a derivative work 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::dynamicRefineBalancedFvMesh
SourceFiles
dynamicRefineBalancedFvMesh.C
dynamicRefineBalancedFvMeshTemplates.C
Authors
T.G. Voskuilen ( https://github.com/tgvoskuilen/meshBalancing )
Daniel Deising <deising@mma.tu-darmstadt.de>
Daniel Rettenmaier <rettenmaier@gsc.tu-darmstadt.de>
All rights reserved.
Description
A fvMesh with built-in multi-criterion refinement
and run-time load balancing.
You may refer to this software as :
//- full bibliographic data to be provided
This code has been developed by :
Daniel Rettenmaier (main developer).
Method Development and Intellectual Property :
T.G. Voskuilen (Purdue University)
Timothée Pourpoint <timothee@purdue.edu> (Purdue University
Daniel Rettenmaier <rettenmaier@gsc.tu-darmstadt.de>
Daniel Deising <deising@mma.tu-darmstadt.de>
Holger Marschall <marschall@csi.tu-darmstadt.de>
Dieter Bothe <bothe@csi.tu-darmstadt.de>
Cameron Tropea <ctropea@sla.tu-darmstadt.de>
School of Aeronautics and Astronautics
Purdue University
Mathematical Modeling and Analysis
Institute for Fluid Mechanics and Aerodynamics
Center of Smart Interfaces
Technische Universitaet Darmstadt
If you use this software for your scientific work or your publications,
please don't forget to acknowledge explicitly the use of it.
\*---------------------------------------------------------------------------*/
#ifndef dynamicRefineBalancedFvMesh_H
#define dynamicRefineBalancedFvMesh_H
#include "dynamicFvMesh.H"
#include "dynamicRefineFvMesh.H"
#include "hexRef8.H"
#include "PackedBoolList.H"
#include "Switch.H"
#include "decompositionMethod.H"
#include "fvMeshDistribute.H"
#include "mapDistributePolyMesh.H"
#include "HashTable.H"
#include "topoSetSource.H"
#include "cellSet.H"
#include "PtrDictionary.H"
#include "dictionaryEntry.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class dynamicRefineBalancedFvMesh Declaration
\*---------------------------------------------------------------------------*/
class dynamicRefineBalancedFvMesh
:
public dynamicRefineFvMesh
{
private:
//- (non-flux)SurfaceFields to map
HashTable<word> mapSurfaceFields_;
//-
volScalarField* internalRefinementFieldPtr_;
volScalarField* targetLevelPtr_;
volScalarField* isLevelPtr_;
//-
HashTable< List<scalar> > fields_;
//-
HashTable< List<scalar> > gradFields_;
//-
HashTable< List<scalar> > curlFields_;
//-
HashTable< dictionary > interface_;
//-
PtrList<entry> refinedRegions_;
//-
Switch enableRefinementControl_;
//-
void updateRefinementField();
//-
void readRefinementDict();
//-
List<scalar> readRefinementPoints();
//- Disallow default bitwise copy construct
dynamicRefineBalancedFvMesh(const dynamicRefineBalancedFvMesh&);
//- Disallow default bitwise assignment
void operator=(const dynamicRefineBalancedFvMesh&);
//-
label topParentID(label p);
//-
label nBufferLayers_;
//-
bool rebalance_;
//- Map non-flux surface<Type>Fields for new internal faces
// (from cell splitting)
template <class T>
void mapNewInternalFaces(const labelList& faceMap);
public:
//- Runtime type information
TypeName("dynamicRefineBalancedFvMesh");
// Constructors
//- Construct from IOobject
explicit dynamicRefineBalancedFvMesh(const IOobject& io);
//- Destructor
virtual ~dynamicRefineBalancedFvMesh();
// Member Functions
//- Update the mesh for both mesh motion and topology change
virtual bool update();
//- Map all fields in time using given map
virtual void mapFields(const mapPolyMesh& mpm);
//- Template to update all volField boundaries
template<class Type> void correctBoundaries();
//-
bool rebalance()
{
return rebalance_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "dynamicRefineBalancedFvMeshTemplates.C"
#endif
#endif
// ************************************************************************* //

217
src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMeshTemplates.C Normal file
View File

@ -0,0 +1,217 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2014 Tyler Voskuilen
\\/ M anipulation |
-------------------------------------------------------------------------------
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 "volMesh.H"
#include "fvPatchField.H"
#include "surfaceFields.H"
template <class T>
void Foam::dynamicRefineBalancedFvMesh::mapNewInternalFaces
(
const labelList& faceMap
)
{
typedef GeometricField<T, fvsPatchField, surfaceMesh> GeoField;
HashTable<GeoField*> sFlds(this->objectRegistry::lookupClass<GeoField>());
const labelUList& owner = this->faceOwner();
const labelUList& neighbour = this->faceNeighbour();
const dimensionedScalar deltaN = 1e-8 / pow(average(this->V()), 1.0 / 3.0);
forAllIter(typename HashTable<GeoField*>, sFlds, iter)
{
GeoField& sFld = *iter();
if (mapSurfaceFields_.found(iter.key()))
{
if (debug)
{
InfoInFunction << iter.key()<< endl;
}
// Create flat version of field
Field<T> tsFld(this->nFaces(), pTraits<T>::zero);
{
forAll(sFld.internalField(), iFace)
{
tsFld[iFace] = sFld.internalField()[iFace];
}
forAll(sFld.boundaryField(), iPatch)
{
const fvsPatchField<T>& pfld = sFld.boundaryField()[iPatch];
label start = pfld.patch().start();
forAll(pfld, faceI)
{
tsFld[faceI+start] = pfld[faceI];
}
}
}
// Loop over all faces
for (label facei = 0; facei < nInternalFaces(); facei++)
{
label oldFacei = faceMap[facei];
//- map surface field on newly generated faces
if (oldFacei == -1)
{
//- find owner and neighbour cell
const cell& ownFaces = this->cells()[owner[facei]];
const cell& neiFaces = this->cells()[neighbour[facei]];
//- loop over all owner/neighbour cell faces
//- and find already mapped ones (master-faces):
T tmpValue = pTraits<T>::zero;
scalar magFld = 0;
label counter = 0;
//- simple averaging of all neighbour master-faces
forAll(ownFaces, iFace)
{
if (faceMap[ownFaces[iFace]] != -1)
{
tmpValue += tsFld[ownFaces[iFace]];
magFld += mag(tsFld[ownFaces[iFace]]);
counter++;
}
}
forAll(neiFaces, iFace)
{
if (faceMap[neiFaces[iFace]] != -1)
{
tmpValue = tmpValue + tsFld[neiFaces[iFace]];
magFld += mag(tsFld[neiFaces[iFace]]);
counter++;
}
}
if(counter > 0)
{
if
(
GeometricField<T, fvsPatchField, surfaceMesh>::typeName
== "surfaceScalarField"
)
{
tmpValue /= counter;
}
else if
(
GeometricField<T, fvsPatchField, surfaceMesh>::typeName
== "surfaceVectorField"
)
{
magFld /= counter;
tmpValue *= magFld/(mag(tmpValue)+deltaN.value());
}
else
{
FatalErrorInFunction
<< "mapping implementation only valid for"
<< " scalar and vector fields! \n Field "
<< sFld.name() << " is of type: "
<< GeometricField<T, fvsPatchField, surfaceMesh>::typeName
<< abort(FatalError);
}
}
sFld[facei] = tmpValue;
}
}
}
}
}
template<class Type>
void Foam::dynamicRefineBalancedFvMesh::correctBoundaries()
{
typedef GeometricField<Type, fvPatchField, volMesh> GeoField;
HashTable<GeoField*> flds(this->objectRegistry::lookupClass<GeoField>());
forAllIter(typename HashTable<GeoField*>, flds, iter)
{
GeoField& fld = *iter();
//mimic "evaluate" but only for coupled patches (processor or cyclic)
// and only for blocking or nonBlocking comms (no scheduled comms)
if
(
Pstream::defaultCommsType == Pstream::commsTypes::blocking
|| Pstream::defaultCommsType == Pstream::commsTypes::nonBlocking
)
{
label nReq = Pstream::nRequests();
forAll(fld.boundaryField(), patchi)
{
if(fld.boundaryField()[patchi].coupled())
{
fld.boundaryFieldRef()[patchi].initEvaluate
(
Pstream::defaultCommsType
);
}
}
// Block for any outstanding requests
if
(
Pstream::parRun()
&& Pstream::defaultCommsType == Pstream::commsTypes::nonBlocking
)
{
Pstream::waitRequests(nReq);
}
forAll(fld.boundaryField(), patchi)
{
if(fld.boundaryField()[patchi].coupled())
{
fld.boundaryFieldRef()[patchi].evaluate
(
Pstream::defaultCommsType
);
}
}
}
else
{
//Scheduled patch updates not supported
FatalErrorIn
(
"dynamicRefineBalancedFvMeshTemplates::correctBoundaries"
) << "Unsuported communications type "
<< Pstream::commsTypeNames[Pstream::defaultCommsType]
<< exit(FatalError);
}
}
}
// ************************************************************************* //

309
src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
View File

@ -46,7 +46,7 @@ namespace Foam
Foam::label Foam::dynamicRefineFvMesh::count
(
const bitSet& l,
const PackedBoolList& l,
const unsigned int val
)
{
@ -58,9 +58,9 @@ Foam::label Foam::dynamicRefineFvMesh::count
n++;
}
// debug also serves to get-around Clang compiler trying to optimise
// debug also serves to get-around Clang compiler trying to optimsie
// out this forAll loop under O3 optimisation
if (debug)
if (debug & 128)
{
Info<< "n=" << n << endl;
}
@ -72,7 +72,7 @@ Foam::label Foam::dynamicRefineFvMesh::count
void Foam::dynamicRefineFvMesh::calculateProtectedCells
(
bitSet& unrefineableCell
PackedBoolList& unrefineableCell
) const
{
if (protectedCell_.empty())
@ -86,7 +86,7 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
unrefineableCell = protectedCell_;
// Get neighbouring cell level
labelList neiLevel(nBoundaryFaces());
labelList neiLevel(nFaces()-nInternalFaces());
for (label facei = nInternalFaces(); facei < nFaces(); facei++)
{
@ -103,10 +103,9 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
forAll(faceNeighbour(), facei)
{
label own = faceOwner()[facei];
bool ownProtected = unrefineableCell.get(own);
label nei = faceNeighbour()[facei];
bool ownProtected = unrefineableCell.test(own);
bool neiProtected = unrefineableCell.test(nei);
bool neiProtected = unrefineableCell.get(nei);
if (ownProtected && (cellLevel[nei] > cellLevel[own]))
{
@ -120,7 +119,7 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
for (label facei = nInternalFaces(); facei < nFaces(); facei++)
{
label own = faceOwner()[facei];
bool ownProtected = unrefineableCell.test(own);
bool ownProtected = unrefineableCell.get(own);
if
(
ownProtected
@ -142,14 +141,16 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
if (seedFace[facei])
{
label own = faceOwner()[facei];
if (unrefineableCell.set(own))
if (unrefineableCell.get(own) == 0)
{
unrefineableCell.set(own, 1);
hasExtended = true;
}
label nei = faceNeighbour()[facei];
if (unrefineableCell.set(nei))
if (unrefineableCell.get(nei) == 0)
{
unrefineableCell.set(nei, 1);
hasExtended = true;
}
}
@ -159,8 +160,9 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
if (seedFace[facei])
{
label own = faceOwner()[facei];
if (unrefineableCell.set(own))
if (unrefineableCell.get(own) == 0)
{
unrefineableCell.set(own, 1);
hasExtended = true;
}
}
@ -205,63 +207,27 @@ void Foam::dynamicRefineFvMesh::readDict()
}
// Refines cells, maps fields and recalculates (an approximate) flux
Foam::autoPtr<Foam::mapPolyMesh>
Foam::dynamicRefineFvMesh::refine
(
const labelList& cellsToRefine
)
void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
{
// Mesh changing engine.
polyTopoChange meshMod(*this);
// Play refinement commands into mesh changer.
meshCutter_.setRefinement(cellsToRefine, meshMod);
// Create mesh (with inflation), return map from old to new mesh.
autoPtr<mapPolyMesh> mapPtr = meshMod.changeMesh(*this, false);
mapPolyMesh& map = *mapPtr;
Info<< "Refined from "
<< returnReduce(map.nOldCells(), sumOp<label>())
<< " to " << globalData().nTotalCells() << " cells." << endl;
if (debug)
{
// Check map.
for (label facei = 0; facei < nInternalFaces(); facei++)
{
label oldFacei = map.faceMap()[facei];
if (oldFacei >= nInternalFaces())
{
FatalErrorInFunction
<< "New internal face:" << facei
<< " fc:" << faceCentres()[facei]
<< " originates from boundary oldFace:" << oldFacei
<< abort(FatalError);
}
}
}
// // Remove the stored tet base points
// tetBasePtIsPtr_.clear();
// // Remove the cell tree
// cellTreePtr_.clear();
// Update fields
updateMesh(map);
dynamicFvMesh::mapFields(mpm);
// Correct the flux for modified/added faces. All the faces which only
// have been renumbered will already have been handled by the mapping.
{
const labelList& faceMap = map.faceMap();
const labelList& reverseFaceMap = map.reverseFaceMap();
const labelList& faceMap = mpm.faceMap();
const labelList& reverseFaceMap = mpm.reverseFaceMap();
// Storage for any master faces. These will be the original faces
// on the coarse cell that get split into four (or rather the
// master face gets modified and three faces get added from the master)
labelHashSet masterFaces(4*cellsToRefine.size());
// Estimate number of faces created
labelHashSet masterFaces
(
max
(
mag(nFaces()-mpm.nOldFaces())/4,
nFaces()/100
)
);
forAll(faceMap, facei)
{
@ -393,8 +359,10 @@ Foam::dynamicRefineFvMesh::refine
}
// Update master faces
for (const label facei : masterFaces)
forAllConstIter(labelHashSet, masterFaces, iter)
{
label facei = iter.key();
if (isInternalFace(facei))
{
phi[facei] = phiU[facei];
@ -415,6 +383,86 @@ Foam::dynamicRefineFvMesh::refine
}
}
// Correct the flux for injected faces - these are the faces which have
// no correspondence to the old mesh (i.e. added without a masterFace, edge
// or point). An example is the internal faces from hexRef8.
{
const labelList& faceMap = mpm.faceMap();
mapNewInternalFaces<scalar>(this->Sf(), this->magSf(), faceMap);
mapNewInternalFaces<vector>(this->Sf(), this->magSf(), faceMap);
// No oriented fields of more complex type
mapNewInternalFaces<sphericalTensor>(faceMap);
mapNewInternalFaces<symmTensor>(faceMap);
mapNewInternalFaces<tensor>(faceMap);
}
}
// Refines cells, maps fields and recalculates (an approximate) flux
Foam::autoPtr<Foam::mapPolyMesh>
Foam::dynamicRefineFvMesh::refine
(
const labelList& cellsToRefine
)
{
// Mesh changing engine.
polyTopoChange meshMod(*this);
// Play refinement commands into mesh changer.
meshCutter_.setRefinement(cellsToRefine, meshMod);
// Create mesh (with inflation), return map from old to new mesh.
//autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
Info<< "Refined from "
<< returnReduce(map().nOldCells(), sumOp<label>())
<< " to " << globalData().nTotalCells() << " cells." << endl;
if (debug)
{
// Check map.
for (label facei = 0; facei < nInternalFaces(); facei++)
{
label oldFacei = map().faceMap()[facei];
if (oldFacei >= nInternalFaces())
{
FatalErrorInFunction
<< "New internal face:" << facei
<< " fc:" << faceCentres()[facei]
<< " originates from boundary oldFace:" << oldFacei
<< abort(FatalError);
}
}
}
// // Remove the stored tet base points
// tetBasePtIsPtr_.clear();
// // Remove the cell tree
// cellTreePtr_.clear();
// Update fields
updateMesh(map);
// Move mesh
/*
pointField newPoints;
if (map().hasMotionPoints())
{
newPoints = map().preMotionPoints();
}
else
{
newPoints = points();
}
movePoints(newPoints);
*/
// Update numbering of cells/vertices.
meshCutter_.updateMesh(map);
@ -422,12 +470,12 @@ Foam::dynamicRefineFvMesh::refine
// Update numbering of protectedCell_
if (protectedCell_.size())
{
bitSet newProtectedCell(nCells());
PackedBoolList newProtectedCell(nCells());
forAll(newProtectedCell, celli)
{
const label oldCelli = map.cellMap()[celli];
newProtectedCell.set(celli, protectedCell_.test(oldCelli));
label oldCelli = map().cellMap()[celli];
newProtectedCell.set(celli, protectedCell_.get(oldCelli));
}
protectedCell_.transfer(newProtectedCell);
}
@ -435,7 +483,7 @@ Foam::dynamicRefineFvMesh::refine
// Debug: Check refinement levels (across faces only)
meshCutter_.checkRefinementLevels(-1, labelList(0));
return mapPtr;
return map;
}
@ -482,21 +530,36 @@ Foam::dynamicRefineFvMesh::unrefine
// Change mesh and generate map.
autoPtr<mapPolyMesh> mapPtr = meshMod.changeMesh(*this, false);
mapPolyMesh& map = *mapPtr;
//autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
Info<< "Unrefined from "
<< returnReduce(map.nOldCells(), sumOp<label>())
<< returnReduce(map().nOldCells(), sumOp<label>())
<< " to " << globalData().nTotalCells() << " cells."
<< endl;
// Update fields
updateMesh(map);
// Move mesh
/*
pointField newPoints;
if (map().hasMotionPoints())
{
newPoints = map().preMotionPoints();
}
else
{
newPoints = points();
}
movePoints(newPoints);
*/
// Correct the flux for modified faces.
{
const labelList& reversePointMap = map.reversePointMap();
const labelList& reverseFaceMap = map.reverseFaceMap();
const labelList& reversePointMap = map().reversePointMap();
const labelList& reverseFaceMap = map().reverseFaceMap();
HashTable<surfaceScalarField*> fluxes
(
@ -585,14 +648,14 @@ Foam::dynamicRefineFvMesh::unrefine
// Update numbering of protectedCell_
if (protectedCell_.size())
{
bitSet newProtectedCell(nCells());
PackedBoolList newProtectedCell(nCells());
forAll(newProtectedCell, celli)
{
label oldCelli = map.cellMap()[celli];
label oldCelli = map().cellMap()[celli];
if (oldCelli >= 0)
{
newProtectedCell.set(celli, protectedCell_.test(oldCelli));
newProtectedCell.set(celli, protectedCell_.get(oldCelli));
}
}
protectedCell_.transfer(newProtectedCell);
@ -601,7 +664,7 @@ Foam::dynamicRefineFvMesh::unrefine
// Debug: Check refinement levels (across faces only)
meshCutter_.checkRefinementLevels(-1, labelList(0));
return mapPtr;
return map;
}
@ -688,7 +751,7 @@ void Foam::dynamicRefineFvMesh::selectRefineCandidates
const scalar lowerRefineLevel,
const scalar upperRefineLevel,
const scalarField& vFld,
bitSet& candidateCell
PackedBoolList& candidateCell
) const
{
// Get error per cell. Is -1 (not to be refined) to >0 (to be refined,
@ -711,7 +774,7 @@ void Foam::dynamicRefineFvMesh::selectRefineCandidates
{
if (cellError[celli] > 0)
{
candidateCell.set(celli);
candidateCell.set(celli, 1);
}
}
}
@ -721,7 +784,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
(
const label maxCells,
const label maxRefinement,
const bitSet& candidateCell
const PackedBoolList& candidateCell
) const
{
// Every refined cell causes 7 extra cells
@ -731,7 +794,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
// Mark cells that cannot be refined since they would trigger refinement
// of protected cells (since 2:1 cascade)
bitSet unrefineableCell;
PackedBoolList unrefineableCell;
calculateProtectedCells(unrefineableCell);
// Count current selection
@ -748,8 +811,11 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
if
(
cellLevel[celli] < maxRefinement
&& candidateCell.test(celli)
&& !unrefineableCell.test(celli)
&& candidateCell.get(celli)
&& (
unrefineableCell.empty()
|| !unrefineableCell.get(celli)
)
)
{
candidates.append(celli);
@ -766,8 +832,11 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
if
(
cellLevel[celli] == level
&& candidateCell.test(celli)
&& !unrefineableCell.test(celli)
&& candidateCell.get(celli)
&& (
unrefineableCell.empty()
|| !unrefineableCell.get(celli)
)
)
{
candidates.append(celli);
@ -802,7 +871,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
(
const scalar unrefineLevel,
const bitSet& markedCell,
const PackedBoolList& markedCell,
const scalarField& pFld
) const
{
@ -815,7 +884,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
// If we have any protected cells make sure they also are not being
// unrefined
bitSet protectedPoint(nPoints());
PackedBoolList protectedPoint(nPoints());
if (protectedCell_.size())
{
@ -828,9 +897,9 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
{
label cellI = pCells[pCellI];
if (protectedCell_.test(cellI))
if (protectedCell_[cellI])
{
protectedPoint.set(pointI);
protectedPoint[pointI] = true;
break;
}
}
@ -870,7 +939,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
forAll(pCells, pCelli)
{
if (markedCell.test(pCells[pCelli]))
if (markedCell.get(pCells[pCelli]))
{
hasMarked = true;
break;
@ -907,7 +976,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
void Foam::dynamicRefineFvMesh::extendMarkedCells
(
bitSet& markedCell
PackedBoolList& markedCell
) const
{
// Mark faces using any marked cell
@ -915,7 +984,7 @@ void Foam::dynamicRefineFvMesh::extendMarkedCells
forAll(markedCell, celli)
{
if (markedCell.test(celli))
if (markedCell.get(celli))
{
const cell& cFaces = cells()[celli];
@ -933,15 +1002,15 @@ void Foam::dynamicRefineFvMesh::extendMarkedCells
{
if (markedFace[facei])
{
markedCell.set(faceOwner()[facei]);
markedCell.set(faceNeighbour()[facei]);
markedCell.set(faceOwner()[facei], 1);
markedCell.set(faceNeighbour()[facei], 1);
}
}
for (label facei = nInternalFaces(); facei < nFaces(); facei++)
{
if (markedFace[facei])
{
markedCell.set(faceOwner()[facei]);
markedCell.set(faceOwner()[facei], 1);
}
}
}
@ -949,7 +1018,7 @@ void Foam::dynamicRefineFvMesh::extendMarkedCells
void Foam::dynamicRefineFvMesh::checkEightAnchorPoints
(
bitSet& protectedCell,
PackedBoolList& protectedCell,
label& nProtected
) const
{
@ -971,7 +1040,7 @@ void Foam::dynamicRefineFvMesh::checkEightAnchorPoints
// Check if cell has already 8 anchor points -> protect cell
if (nAnchorPoints[celli] == 8)
{
if (protectedCell.set(celli))
if (protectedCell.set(celli, true))
{
nProtected++;
}
@ -988,9 +1057,9 @@ void Foam::dynamicRefineFvMesh::checkEightAnchorPoints
forAll(protectedCell, celli)
{
if (!protectedCell.test(celli) && nAnchorPoints[celli] != 8)
if (!protectedCell[celli] && nAnchorPoints[celli] != 8)
{
protectedCell.set(celli);
protectedCell.set(celli, true);
nProtected++;
}
}
@ -1005,7 +1074,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
meshCutter_(*this),
dumpLevel_(false),
nRefinementIterations_(0),
protectedCell_(nCells(), false)
protectedCell_(nCells(), 0)
{
// Read static part of dictionary
readDict();
@ -1034,7 +1103,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
{
label celli = pCells[i];
if (!protectedCell_.test(celli))
if (!protectedCell_.get(celli))
{
if (pointLevel[pointi] <= cellLevel[celli])
{
@ -1042,7 +1111,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
if (nAnchors[celli] > 8)
{
protectedCell_.set(celli);
protectedCell_.set(celli, 1);
nProtected++;
}
}
@ -1105,9 +1174,9 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
{
if (protectedFace[facei])
{
protectedCell_.set(faceOwner()[facei]);
protectedCell_.set(faceOwner()[facei], 1);
nProtected++;
protectedCell_.set(faceNeighbour()[facei]);
protectedCell_.set(faceNeighbour()[facei], 1);
nProtected++;
}
}
@ -1115,7 +1184,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
{
if (protectedFace[facei])
{
protectedCell_.set(faceOwner()[facei]);
protectedCell_.set(faceOwner()[facei], 1);
nProtected++;
}
}
@ -1127,7 +1196,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
if (cFaces.size() < 6)
{
if (protectedCell_.set(celli))
if (protectedCell_.set(celli, 1))
{
nProtected++;
}
@ -1138,7 +1207,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
{
if (faces()[cFaces[cFacei]].size() < 4)
{
if (protectedCell_.set(celli))
if (protectedCell_.set(celli, 1))
{
nProtected++;
}
@ -1158,6 +1227,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
}
else
{
cellSet protectedCells(*this, "protectedCells", nProtected);
forAll(protectedCell_, celli)
{
@ -1188,9 +1258,9 @@ Foam::dynamicRefineFvMesh::~dynamicRefineFvMesh()
bool Foam::dynamicRefineFvMesh::update()
{
// Re-read dictionary. Usually small so takes trivial amount of time
// compared to actual refinement. Also very useful to be able to modify
// on-the-fly.
// Re-read dictionary. Choosen since usually -small so trivial amount
// of time compared to actual refinement. Also very useful to be able
// to modify on-the-fly.
dictionary refineDict
(
IOdictionary
@ -1264,18 +1334,16 @@ bool Foam::dynamicRefineFvMesh::update()
refineDict.get<scalar>("lowerRefineLevel");
const scalar upperRefineLevel =
refineDict.get<scalar>("upperRefineLevel");
const scalar unrefineLevel = refineDict.lookupOrDefault<scalar>
(
"unrefineLevel",
GREAT
);
const label nBufferLayers =
refineDict.get<label>("nBufferLayers");
// Cells marked for refinement or otherwise protected from unrefinement.
bitSet refineCell(nCells());
PackedBoolList refineCell(nCells());
// Determine candidates for refinement (looking at field only)
selectRefineCandidates
@ -1316,7 +1384,7 @@ bool Foam::dynamicRefineFvMesh::update()
const labelList& cellMap = map().cellMap();
const labelList& reverseCellMap = map().reverseCellMap();
bitSet newRefineCell(cellMap.size());
PackedBoolList newRefineCell(cellMap.size());
forAll(cellMap, celli)
{
@ -1324,15 +1392,15 @@ bool Foam::dynamicRefineFvMesh::update()
if (oldCelli < 0)
{
newRefineCell.set(celli);
newRefineCell.set(celli, 1);
}
else if (reverseCellMap[oldCelli] != celli)
{
newRefineCell.set(celli);
newRefineCell.set(celli, 1);
}
else
{
newRefineCell.set(celli, refineCell.test(oldCelli));
newRefineCell.set(celli, refineCell.get(oldCelli));
}
}
refineCell.transfer(newRefineCell);
@ -1380,7 +1448,7 @@ bool Foam::dynamicRefineFvMesh::update()
if ((nRefinementIterations_ % 10) == 0)
{
// Compact refinement history occasionally (how often?).
// Compact refinement history occassionally (how often?).
// Unrefinement causes holes in the refinementHistory.
const_cast<refinementHistory&>(meshCutter().history()).compact();
}
@ -1403,8 +1471,7 @@ bool Foam::dynamicRefineFvMesh::writeObject
(
IOstream::streamFormat fmt,
IOstream::versionNumber ver,
IOstream::compressionType cmp,
const bool valid
IOstream::compressionType cmp
) const
{
// Force refinement data to go to the current time directory.
@ -1412,8 +1479,8 @@ bool Foam::dynamicRefineFvMesh::writeObject
bool writeOk =
(
dynamicFvMesh::writeObject(fmt, ver, cmp, valid)
&& meshCutter_.write(valid)
dynamicFvMesh::writeObject(fmt, ver, cmp)
&& meshCutter_.write()
);
if (dumpLevel_)
@ -1430,7 +1497,7 @@ bool Foam::dynamicRefineFvMesh::writeObject
false
),
*this,
dimensionedScalar(dimless, Zero)
dimensionedScalar("level", dimless, 0)
);
const labelList& cellLevel = meshCutter_.cellLevel();

44
src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H
View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2017-2018 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -56,6 +56,7 @@ Description
(rho*phi none) //none)
(ghf none) //NaN) //none)
);
// Write the refinement level as a volScalarField
dumpLevel true;
@ -117,10 +118,10 @@ protected:
//- Refine cells. Update mesh and fields.
autoPtr<mapPolyMesh> refine(const labelList&);
virtual autoPtr<mapPolyMesh> refine(const labelList&);
//- Unrefine cells. Gets passed in centre points of cells to combine.
autoPtr<mapPolyMesh> unrefine(const labelList&);
virtual autoPtr<mapPolyMesh> unrefine(const labelList&);
// Selection of cells to un/refine
@ -185,6 +186,32 @@ protected:
label& nProtected
) const;
//- Map single non-flux surface<Type>Field
// for new internal faces (e.g. AMR refine). This currently
// interpolates values from surrounding faces (faces on
// neighbouring cells) that do have a value.
template <class T>
void mapNewInternalFaces
(
const labelList& faceMap,
GeometricField<T, fvsPatchField, surfaceMesh>&
);
//- Correct surface fields for new faces
template <class T>
void mapNewInternalFaces(const labelList& faceMap);
//- Correct surface fields for new faces. Converts any oriented
// fields into non-oriented (i.e. phi -> Uf) before mapping
template <class T>
void mapNewInternalFaces
(
const surfaceVectorField& Sf,
const surfaceScalarField& magSf,
const labelList& faceMap
);
private:
//- No copy construct
@ -232,6 +259,9 @@ public:
//- Update the mesh for both mesh motion and topology change
virtual bool update();
//- Map all fields in time using given map.
virtual void mapFields(const mapPolyMesh& mpm);
// Writing
@ -243,16 +273,20 @@ public:
IOstream::compressionType cmp,
const bool valid
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "dynamicRefineFvMeshTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

196
src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMeshTemplates.C Normal file
View File

@ -0,0 +1,196 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2018 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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 "surfaceFields.H"
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
template<class T>
void Foam::dynamicRefineFvMesh::mapNewInternalFaces
(
const labelList& faceMap,
GeometricField<T, fvsPatchField, surfaceMesh>& sFld
)
{
typedef GeometricField<T, fvsPatchField, surfaceMesh> GeoField;
//- Make flat field for ease of looping
Field<T> tsFld(this->nFaces(), pTraits<T>::zero);
SubField<T>(tsFld, this->nInternalFaces()) = sFld.internalField();
const typename GeoField::Boundary& bFld = sFld.boundaryField();
forAll(bFld, patchi)
{
label facei = this->boundaryMesh()[patchi].start();
for (const T& val : bFld[patchi])
{
tsFld[facei++] = val;
}
}
const labelUList& owner = this->faceOwner();
const labelUList& neighbour = this->faceNeighbour();
const cellList& cells = this->cells();
for (label facei = 0; facei < nInternalFaces(); facei++)
{
label oldFacei = faceMap[facei];
// Map surface field on newly generated faces by obtaining the
// hull of the outside faces
if (oldFacei == -1)
{
// Loop over all owner/neighbour cell faces
// and find already mapped ones (master-faces):
T tmpValue = pTraits<T>::zero;
label counter = 0;
const cell& ownFaces = cells[owner[facei]];
for (auto ownFacei : ownFaces)
{
if (faceMap[ownFacei] != -1)
{
tmpValue += tsFld[ownFacei];
counter++;
}
}
const cell& neiFaces = cells[neighbour[facei]];
for (auto neiFacei : neiFaces)
{
if (faceMap[neiFacei] != -1)
{
tmpValue += tsFld[neiFacei];
counter++;
}
}
if (counter > 0)
{
sFld[facei] = tmpValue/counter;
}
}
}
}
template<class T>
void Foam::dynamicRefineFvMesh::mapNewInternalFaces
(
const labelList& faceMap
)
{
typedef GeometricField<T, fvsPatchField, surfaceMesh> GeoField;
HashTable<GeoField*> sFlds(this->objectRegistry::lookupClass<GeoField>());
forAllIter(typename HashTable<GeoField*>, sFlds, iter)
{
//if (mapSurfaceFields_.found(iter.key()))
{
if (debug)
{
Info<< "dynamicRefineFvMesh::mapNewInternalFaces():"
<< " Mapping new internal faces by interpolation on "
<< iter.key()<< endl;
}
GeoField& sFld = *iter();
if (sFld.oriented()())
{
WarningInFunction << "Ignoring mapping oriented field "
<< sFld.name() << " since of type " << sFld.type()
<< endl;
}
else
{
mapNewInternalFaces(faceMap, sFld);
}
}
}
}
template<class T>
void Foam::dynamicRefineFvMesh::mapNewInternalFaces
(
const surfaceVectorField& Sf,
const surfaceScalarField& magSf,
const labelList& faceMap
)
{
typedef GeometricField<T, fvsPatchField, surfaceMesh> GeoField;
HashTable<GeoField*> sFlds(this->objectRegistry::lookupClass<GeoField>());
forAllIter(typename HashTable<GeoField*>, sFlds, iter)
{
//if (mapSurfaceFields_.found(iter.key()))
{
if (debug)
{
Info<< "dynamicRefineFvMesh::mapNewInternalFaces():"
<< " Mapping new internal faces by interpolation on "
<< iter.key()<< endl;
}
GeoField& sFld = *iter();
if (sFld.oriented()())
{
if (debug)
{
Info<< "dynamicRefineFvMesh::mapNewInternalFaces(): "
<< "Converting oriented field " << iter.key()
<< " to intensive field and mapping" << endl;
}
// Assume any oriented field is face area weighted (i.e. a flux)
// Convert to intensive (& oriented) before mapping. Untested.
typedef GeometricField
<
typename outerProduct<vector, T>::type,
fvsPatchField,
surfaceMesh
> NormalGeoField;
// Convert to intensive and non oriented
NormalGeoField fFld(sFld*Sf/Foam::sqr(magSf));
// Interpolate
mapNewInternalFaces(faceMap, fFld);
// Convert back to extensive and oriented
sFld = (fFld & Sf);
}
else
{
mapNewInternalFaces(faceMap, sFld);
}
}
}
}
// ************************************************************************* //

2
src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
View File

@ -210,7 +210,7 @@ Foam::label Foam::hexRef8::addInternalFace
nei, // neighbour
-1, // master point
-1, // master edge
meshFacei, // master face for addition
-1, // master face for addition
false, // flux flip
-1, // patch for face
-1, // zone for face

10
src/finiteVolume/fvMesh/fvMesh.C
View File

@ -571,6 +571,13 @@ const Foam::lduAddressing& Foam::fvMesh::lduAddr() const
{
if (!lduPtr_)
{
if (debug)
{
InfoInFunction
<< " calculating fvMeshLduAddressing from nFaces:"
<< nFaces() << endl;
}
lduPtr_ = new fvMeshLduAddressing(*this);
}
@ -821,6 +828,9 @@ void Foam::fvMesh::updateMesh(const mapPolyMesh& mpm)
// Update polyMesh. This needs to keep volume existent!
polyMesh::updateMesh(mpm);
// Our slice of the addressing is no longer valid
deleteDemandDrivenData(lduPtr_);
if (VPtr_)
{
// Grab old time volumes if the time has been incremented

1
tutorials/multiphase/interFoam/laminar/damBreakWithObstacle/constant/dynamicMeshDict
View File

@ -49,6 +49,7 @@ correctFluxes
(rhoPhi none)
(alphaPhi0.water none)
(ghf none)
(alphaPhiUn none)
);
// Write the refinement level as a volScalarField