ENH: refactored mesh movePoints to enable meshPhi update in the fvGeometryScheme

This commit is contained in:
Andrew Heather
2022-03-16 20:03:53 +00:00
committed by Andrew Heather
parent 5a1307f41a
commit 823641ab9b
13 changed files with 176 additions and 122 deletions

View File

@ -1156,10 +1156,7 @@ const Foam::pointField& Foam::polyMesh::oldCellCentres() const
}
Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
(
const pointField& newPoints
)
void Foam::polyMesh::movePoints(const pointField& newPoints)
{
DebugInFunction
<< "Moving points for time " << time().value()
@ -1227,11 +1224,17 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
tetBasePtIsPtr_->eventNo() = getEvent();
}
tmp<scalarField> sweptVols = primitiveMesh::movePoints
(
points_,
oldPoints()
);
// Currently a no-op; earlier versions set meshPhi and call
// primitiveMesh::clearGeom
(void)primitiveMesh::movePoints(points_, oldPoints());
// Update the mesh geometry (via fvGeometryScheme)
// - updateGeom is virtual -> calls fvMesh::updateGeom (or higher)
// - fvMesh::updateGeom defers to surfaceInterpolation::updateGeom(),
// which defers to fvGeometryScheme::movePoints()
// - set the mesh flux
// - clear out/recalculate stale geometry
updateGeom();
// Adjust parallel shared points
if (globalMeshDataPtr_)
@ -1279,8 +1282,6 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
// e.g. fvMesh::write since meshPhi not yet complete.
polyMesh::write();
}
return sweptVols;
}

View File

@ -556,8 +556,8 @@ public:
return (moving() || topoChanging());
}
//- Move points, returns volumes swept by faces in motion
virtual tmp<scalarField> movePoints(const pointField&);
//- Move points
virtual void movePoints(const pointField&);
//- Reset motion
void resetMotion() const;

View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -136,7 +137,7 @@ bool Foam::primitiveMesh::calcPointOrder
// Internal points are points that are not used by a boundary face.
// Map from old to new position
oldToNew.setSize(nPoints);
oldToNew.resize_nocopy(nPoints);
oldToNew = -1;
@ -144,14 +145,12 @@ bool Foam::primitiveMesh::calcPointOrder
// from 0 inside oldToNew. (shifted up later on)
label nBoundaryPoints = 0;
for (label facei = nInternalFaces; facei < faces.size(); facei++)
for (label facei = nInternalFaces; facei < faces.size(); ++facei)
{
const face& f = faces[facei];
forAll(f, fp)
for (label pointi : f)
{
label pointi = f[fp];
if (oldToNew[pointi] == -1)
{
oldToNew[pointi] = nBoundaryPoints++;
@ -184,10 +183,8 @@ bool Foam::primitiveMesh::calcPointOrder
{
const face& f = faces[facei];
forAll(f, fp)
for (label pointi : f)
{
label pointi = f[fp];
if (oldToNew[pointi] == -1)
{
if (pointi >= nInternalPoints)
@ -318,35 +315,17 @@ void Foam::primitiveMesh::resetGeometry
}
Foam::tmp<Foam::scalarField> Foam::primitiveMesh::movePoints
void Foam::primitiveMesh::movePoints
(
const pointField& newPoints,
const pointField& oldPoints
)
{
if (newPoints.size() < nPoints() || oldPoints.size() < nPoints())
{
FatalErrorInFunction
<< "Cannot move points: size of given point list smaller "
<< "than the number of active points"
<< abort(FatalError);
}
// Create swept volumes
const faceList& fcs = faces();
auto tsweptVols = tmp<scalarField>::New(fcs.size());
auto& sweptVols = tsweptVols.ref();
forAll(fcs, facei)
{
sweptVols[facei] = fcs[facei].sweptVol(oldPoints, newPoints);
}
// Note: the following clearout is now handled by the fvGeometryScheme
// triggered by the call to updateGeom() in polyMesh::movePoints
// Force recalculation of all geometric data with new points
clearGeom();
return tsweptVols;
//clearGeom();
}
@ -361,7 +340,6 @@ const Foam::cellShapeList& Foam::primitiveMesh::cellShapes() const
}
void Foam::primitiveMesh::updateGeom()
{
if (!faceCentresPtr_ || !faceAreasPtr_)

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2018 OpenCFD Ltd.
Copyright (C) 2018-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -598,8 +598,8 @@ public:
// Mesh motion
//- Move points, returns volumes swept by faces in motion
tmp<scalarField> movePoints
//- Move points
void movePoints
(
const pointField& p,
const pointField& oldP

View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -52,11 +53,11 @@ void Foam::primitiveMesh::calcCellCentresAndVols() const
}
// set the accumulated cell centre to zero vector
cellCentresPtr_ = new vectorField(nCells());
cellCentresPtr_ = new vectorField(nCells(), Zero);
vectorField& cellCtrs = *cellCentresPtr_;
// Initialise cell volumes to 0
cellVolumesPtr_ = new scalarField(nCells());
cellVolumesPtr_ = new scalarField(nCells(), Zero);
scalarField& cellVols = *cellVolumesPtr_;
// Make centres and volumes

View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.

View File

@ -55,11 +55,14 @@ Foam::basicFvGeometryScheme::basicFvGeometryScheme
void Foam::basicFvGeometryScheme::movePoints()
{
fvGeometryScheme::movePoints();
if (debug)
{
Pout<< "basicFvGeometryScheme::movePoints() : "
<< "recalculating primitiveMesh centres" << endl;
}
// Use lower level to calculate the geometry
const_cast<fvMesh&>(mesh_).primitiveMesh::updateGeom();
}
@ -74,9 +77,8 @@ Foam::tmp<Foam::surfaceScalarField> Foam::basicFvGeometryScheme::weights() const
<< endl;
}
tmp<surfaceScalarField> tweights
(
new surfaceScalarField
auto tweights =
tmp<surfaceScalarField>::New
(
IOobject
(
@ -89,9 +91,9 @@ Foam::tmp<Foam::surfaceScalarField> Foam::basicFvGeometryScheme::weights() const
),
mesh_,
dimless
)
);
surfaceScalarField& weights = tweights.ref();
);
auto& weights = tweights.ref();
weights.setOriented();
// Set local references to mesh data
@ -128,7 +130,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::basicFvGeometryScheme::weights() const
}
}
surfaceScalarField::Boundary& wBf = weights.boundaryFieldRef();
auto& wBf = weights.boundaryFieldRef();
forAll(mesh_.boundary(), patchi)
{
@ -159,9 +161,8 @@ Foam::basicFvGeometryScheme::deltaCoeffs() const
// needed to make sure deltaCoeffs are calculated for parallel runs.
(void)mesh_.weights();
tmp<surfaceScalarField> tdeltaCoeffs
(
new surfaceScalarField
auto tdeltaCoeffs =
tmp<surfaceScalarField>::New
(
IOobject
(
@ -174,14 +175,14 @@ Foam::basicFvGeometryScheme::deltaCoeffs() const
),
mesh_,
dimless/dimLength
)
);
surfaceScalarField& deltaCoeffs = tdeltaCoeffs.ref();
);
auto& deltaCoeffs = tdeltaCoeffs.ref();
deltaCoeffs.setOriented();
// Set local references to mesh data
const volVectorField& C = mesh_.C();
const vectorField& C = mesh_.cellCentres();
const labelUList& owner = mesh_.owner();
const labelUList& neighbour = mesh_.neighbour();
@ -190,8 +191,7 @@ Foam::basicFvGeometryScheme::deltaCoeffs() const
deltaCoeffs[facei] = 1.0/mag(C[neighbour[facei]] - C[owner[facei]]);
}
surfaceScalarField::Boundary& deltaCoeffsBf =
deltaCoeffs.boundaryFieldRef();
auto& deltaCoeffsBf = deltaCoeffs.boundaryFieldRef();
forAll(deltaCoeffsBf, patchi)
{
@ -220,9 +220,8 @@ Foam::basicFvGeometryScheme::nonOrthDeltaCoeffs() const
// needed to make sure deltaCoeffs are calculated for parallel runs.
weights();
tmp<surfaceScalarField> tnonOrthDeltaCoeffs
(
new surfaceScalarField
auto tnonOrthDeltaCoeffs =
tmp<surfaceScalarField>::New
(
IOobject
(
@ -235,9 +234,9 @@ Foam::basicFvGeometryScheme::nonOrthDeltaCoeffs() const
),
mesh_,
dimless/dimLength
)
);
surfaceScalarField& nonOrthDeltaCoeffs = tnonOrthDeltaCoeffs.ref();
);
auto& nonOrthDeltaCoeffs = tnonOrthDeltaCoeffs.ref();
nonOrthDeltaCoeffs.setOriented();
@ -266,8 +265,7 @@ Foam::basicFvGeometryScheme::nonOrthDeltaCoeffs() const
nonOrthDeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta));
}
surfaceScalarField::Boundary& nonOrthDeltaCoeffsBf =
nonOrthDeltaCoeffs.boundaryFieldRef();
auto& nonOrthDeltaCoeffsBf = nonOrthDeltaCoeffs.boundaryFieldRef();
forAll(nonOrthDeltaCoeffsBf, patchi)
{
@ -306,9 +304,8 @@ Foam::basicFvGeometryScheme::nonOrthCorrectionVectors() const
<< endl;
}
tmp<surfaceVectorField> tnonOrthCorrectionVectors
(
new surfaceVectorField
auto tnonOrthCorrectionVectors =
tmp<surfaceVectorField>::New
(
IOobject
(
@ -321,9 +318,9 @@ Foam::basicFvGeometryScheme::nonOrthCorrectionVectors() const
),
mesh_,
dimless
)
);
surfaceVectorField& corrVecs = tnonOrthCorrectionVectors.ref();
);
auto& corrVecs = tnonOrthCorrectionVectors.ref();
corrVecs.setOriented();
// Set local references to mesh data
@ -337,8 +334,8 @@ Foam::basicFvGeometryScheme::nonOrthCorrectionVectors() const
forAll(owner, facei)
{
vector unitArea = Sf[facei]/magSf[facei];
vector delta = C[neighbour[facei]] - C[owner[facei]];
vector unitArea(Sf[facei]/magSf[facei]);
vector delta(C[neighbour[facei]] - C[owner[facei]]);
corrVecs[facei] = unitArea - delta*NonOrthDeltaCoeffs[facei];
}
@ -347,7 +344,7 @@ Foam::basicFvGeometryScheme::nonOrthCorrectionVectors() const
// and calculated consistently with internal corrections for
// coupled patches
surfaceVectorField::Boundary& corrVecsBf = corrVecs.boundaryFieldRef();
auto& corrVecsBf = corrVecs.boundaryFieldRef();
forAll(corrVecsBf, patchi)
{
@ -361,7 +358,7 @@ Foam::basicFvGeometryScheme::nonOrthCorrectionVectors() const
}
else
{
const fvsPatchScalarField& patchNonOrthDeltaCoeffs =
const auto& patchNonOrthDeltaCoeffs =
NonOrthDeltaCoeffs.boundaryField()[patchi];
const vectorField patchDeltas(mesh_.boundary()[patchi].delta());

View File

@ -26,6 +26,8 @@ License
\*---------------------------------------------------------------------------*/
#include "fvGeometryScheme.H"
#include "fvMesh.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -37,10 +39,63 @@ namespace Foam
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
bool Foam::fvGeometryScheme::setMeshPhi() const
{
if (!mesh_.moving())
{
return false;
}
const pointField& oldPoints = mesh_.oldPoints();
const pointField& currPoints = mesh_.points();
if (oldPoints.size() != currPoints.size())
{
FatalErrorInFunction
<< "Old and current points sizes must be the same. "
<< "Old points:" << oldPoints.size()
<< " Current points:" << currPoints.size()
<< abort(FatalError);
}
const faceList& faces = mesh_.faces();
const scalar rdt = 1.0/mesh_.time().deltaTValue();
auto& meshPhi = const_cast<fvMesh&>(mesh_).setPhi();
auto& meshPhii = meshPhi.primitiveFieldRef();
forAll(meshPhii, facei)
{
const face& f = faces[facei];
meshPhii[facei] = f.sweptVol(oldPoints, currPoints)*rdt;
}
auto& meshPhiBf = meshPhi.boundaryFieldRef();
for (auto& meshPhip : meshPhiBf)
{
if (!meshPhip.size())
{
// Empty patches
continue;
}
const auto& pp = meshPhip.patch().patch();
forAll(pp, facei)
{
const face& f = pp[facei];
meshPhip[facei] = f.sweptVol(oldPoints, currPoints)*rdt;
}
}
return true;
}
// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
Foam::tmp<Foam::fvGeometryScheme>
Foam::fvGeometryScheme::New
Foam::tmp<Foam::fvGeometryScheme> Foam::fvGeometryScheme::New
(
const fvMesh& mesh,
const dictionary& dict,
@ -55,10 +110,7 @@ Foam::fvGeometryScheme::New
: dict.getOrDefault<word>("type", defaultScheme)
);
if (debug)
{
InfoInFunction << "Geometry scheme = " << schemeName << endl;
}
DebugInFunction << "Geometry scheme = " << schemeName << endl;
auto* ctorPtr = dictConstructorTable(schemeName);
@ -77,4 +129,20 @@ Foam::fvGeometryScheme::New
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::fvGeometryScheme::movePoints()
{
if (mesh_.moving())
{
// Set the mesh motion fluxes
setMeshPhi();
// Clear out old geometry
// Note: this recreates the old primitiveMesh::movePoints behaviour
const_cast<fvMesh&>(mesh_).primitiveMesh::clearGeom();
}
}
// ************************************************************************* //

View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -69,10 +69,18 @@ class fvGeometryScheme
protected:
// Protected Data
//- Hold reference to mesh
const fvMesh& mesh_;
// Protected Member Functions
//- Set the mesh motion flux
bool setMeshPhi() const;
public:
//- Runtime type information
@ -127,8 +135,7 @@ public:
}
//- Update basic geometric properties from provided points
virtual void movePoints()
{}
virtual void movePoints();
//- Return linear difference weighting factors
virtual tmp<surfaceScalarField> weights() const = 0;

View File

@ -387,6 +387,9 @@ Foam::highAspectRatioFvGeometryScheme::highAspectRatioFvGeometryScheme
void Foam::highAspectRatioFvGeometryScheme::movePoints()
{
//basicFvGeometryScheme::movePoints();
fvGeometryScheme::movePoints();
if (debug)
{
Pout<< "highAspectRatioFvGeometryScheme::movePoints() : "

View File

@ -157,6 +157,8 @@ Foam::stabilisedFvGeometryScheme::stabilisedFvGeometryScheme
void Foam::stabilisedFvGeometryScheme::movePoints()
{
fvGeometryScheme::movePoints();
if (debug)
{
Pout<< "stabilisedFvGeometryScheme::movePoints() : "

View File

@ -90,18 +90,22 @@ void Foam::fvMesh::updateGeomNotOldVol()
{
(void)V();
}
if (haveSf)
{
(void)Sf();
}
if (haveMagSf)
{
(void)magSf();
}
if (haveCP)
{
(void)C();
}
if (haveCf)
{
(void)Cf();
@ -277,7 +281,7 @@ bool Foam::fvMesh::init(const bool doInit)
// doing anything with primitiveMesh::cellCentres etc.
(void)geometry();
// Intialise my data
// Initialise my data
polyMesh::init(doInit);
}
@ -285,6 +289,10 @@ bool Foam::fvMesh::init(const bool doInit)
// and set the storage of V00
if (fileHandler().isFile(time().timePath()/dbDir()/"V0"))
{
// Set the moving flag early in case the demand-driven geometry
// construction checks for it
moving(true);
V0Ptr_ = new DimensionedField<scalar, volMesh>
(
IOobject
@ -306,6 +314,10 @@ bool Foam::fvMesh::init(const bool doInit)
// mesh to be moving
if (fileHandler().isFile(time().timePath()/dbDir()/"meshPhi"))
{
// Set the moving flag early in case the demand-driven geometry
// construction checks for it
moving(true);
phiPtr_ = new surfaceScalarField
(
IOobject
@ -339,8 +351,6 @@ bool Foam::fvMesh::init(const bool doInit)
V()
);
}
moving(true);
}
// Assume something changed
@ -643,7 +653,7 @@ void Foam::fvMesh::removeFvBoundary()
Foam::polyMesh::readUpdateState Foam::fvMesh::readUpdate()
{
DebugInFunction << "Updating fvMesh. ";
DebugInFunction << "Updating fvMesh";
polyMesh::readUpdateState state = polyMesh::readUpdate();
@ -861,7 +871,8 @@ void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap)
}
Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
void Foam::fvMesh::movePoints(const pointField& p)
{
DebugInFunction << endl;
@ -873,12 +884,8 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
}
// Move the polyMesh and set the mesh motion fluxes to the swept-volumes
scalar rDeltaT = 1.0/time().deltaTValue();
tmp<scalarField> tsweptVols = polyMesh::movePoints(p);
scalarField& sweptVols = tsweptVols.ref();
// Move the polyMesh and initialise the mesh motion fluxes field
// Note: mesh flux updated by the fvGeometryScheme
if (!phiPtr_)
{
@ -895,7 +902,7 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
false
),
*this,
dimVolume/dimTime
dimensionedScalar(dimVolume/dimTime, Zero)
);
}
else
@ -907,19 +914,7 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
}
}
surfaceScalarField& phi = *phiPtr_;
phi.primitiveFieldRef() =
scalarField::subField(sweptVols, nInternalFaces());
phi.primitiveFieldRef() *= rDeltaT;
const fvPatchList& patches = boundary();
surfaceScalarField::Boundary& phibf = phi.boundaryFieldRef();
forAll(patches, patchi)
{
phibf[patchi] = patches[patchi].patchSlice(sweptVols);
phibf[patchi] *= rDeltaT;
}
polyMesh::movePoints(p);
// Update or delete the local geometric properties as early as possible so
// they can be used if necessary. These get recreated here instead of
@ -929,20 +924,21 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
// should use the local geometric properties.
updateGeomNotOldVol();
// Update other local data
boundary_.movePoints();
surfaceInterpolation::movePoints();
// Clear weights, deltaCoeffs, nonOrthoDeltaCoeffs, nonOrthCorrectionVectors
surfaceInterpolation::clearOut();
meshObject::movePoints<fvMesh>(*this);
meshObject::movePoints<lduMesh>(*this);
return tsweptVols;
}
void Foam::fvMesh::updateGeom()
{
DebugInFunction << endl;
// Let surfaceInterpolation handle geometry calculation. Note: this does
// lower levels updateGeom
surfaceInterpolation::updateGeom();

View File

@ -469,7 +469,7 @@ public:
using surfaceInterpolation::movePoints;
//- Move points, returns volumes swept by faces in motion
virtual tmp<scalarField> movePoints(const pointField&);
virtual void movePoints(const pointField&);
//- Update all geometric data. This gets redirected up from
//- primitiveMesh level