Commit 46dbfabd authored by Mattijs Janssens's avatar Mattijs Janssens Committed by Andrew Heather
Browse files

ENH: primitiveMesh: make geometry calculation runtime selectable

This adds a 'geometry' scheme section to the system/fvSchemes:

geometry
{
    type            highAspectRatio;
}

These 'fvGeometryMethod's are used to calculate
- deltaCoeffs
- nonOrthoCoeffs
etc and can even modify the basic face/cellCentres calculation.
parent 4a166c6f
......@@ -704,7 +704,7 @@ int main(int argc, char *argv[])
regionName,
runTimeExtruded.constant(),
runTimeExtruded,
IOobject::NO_READ,
IOobject::READ_IF_PRESENT, // Read fv* if present
IOobject::AUTO_WRITE,
false
),
......
......@@ -154,6 +154,7 @@ int main(int argc, char *argv[])
"cellShapes",
"cellVolume",
"cellVolumeRatio",
"cellAspectRatio",
"minTetVolume",
"minPyrVolume",
"cellRegion",
......
......@@ -7,6 +7,7 @@
#include "tetPointRef.H"
#include "regionSplit.H"
#include "wallDist.H"
#include "cellAspectRatio.H"
using namespace Foam;
......@@ -318,6 +319,32 @@ void Foam::writeFields
aspectRatio.write();
}
if (selectedFields.found("cellAspectRatio"))
{
volScalarField aspectRatio
(
IOobject
(
"cellAspectRatio",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
mesh,
dimensionedScalar(dimless, Zero),
zeroGradientFvPatchScalarField::typeName
);
aspectRatio.ref().field() = cellAspectRatio(mesh);
aspectRatio.correctBoundaryConditions();
Info<< " Writing approximate aspect ratio to "
<< aspectRatio.name() << endl;
aspectRatio.write();
}
// cell type
// ~~~~~~~~~
......
......@@ -680,7 +680,7 @@ autoPtr<mapPolyMesh> createRegionMesh
regionName,
mesh.time().timeName(),
mesh.time(),
IOobject::NO_READ,
IOobject::READ_IF_PRESENT, // read fv* if present
IOobject::AUTO_WRITE
),
mesh
......
......@@ -36,6 +36,7 @@ License
#include "pointSet.H"
#include "faceSet.H"
#include "cellSet.H"
#include "basicFvGeometryScheme.H"
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
......@@ -254,6 +255,19 @@ Foam::autoPtr<Foam::fvMesh> Foam::loadOrCreateMesh
auto meshPtr = autoPtr<fvMesh>::New(io);
fvMesh& mesh = *meshPtr;
// Make sure to use a non-parallel geometry calculation method
{
tmp<fvGeometryScheme> basicGeometry
(
fvGeometryScheme::New
(
mesh,
dictionary(),
basicFvGeometryScheme::typeName
)
);
mesh.geometry(basicGeometry);
}
// Sync patches
......
......@@ -92,6 +92,7 @@ Usage
#include "zeroGradientFvPatchFields.H"
#include "topoSet.H"
#include "regionProperties.H"
#include "basicFvGeometryScheme.H"
#include "parFvFieldReconstructor.H"
#include "parLagrangianRedistributor.H"
......@@ -152,6 +153,30 @@ scalar getMergeDistance
}
void setBasicGeometry(fvMesh& mesh)
{
// Set the fvGeometryScheme to basic since it does not require
// any parallel communication to construct the geometry. During
// redistributePar one might temporarily end up with processors
// with zero procBoundaries. Normally procBoundaries trigger geometry
// calculation (e.g. send over cellCentres) so on the processors without
// procBoundaries this will not happen. The call to the geometry calculation
// is not synchronised and this might lead to a hang for geometry schemes
// that do require synchronisation
tmp<fvGeometryScheme> basicGeometry
(
fvGeometryScheme::New
(
mesh,
dictionary(),
basicFvGeometryScheme::typeName
)
);
mesh.geometry(basicGeometry);
}
void printMeshData(const polyMesh& mesh)
{
// Collect all data on master
......@@ -2675,6 +2700,10 @@ int main(int argc, char *argv[])
);
fvMesh& mesh = meshPtr();
// Use basic geometry calculation to avoid synchronisation
// problems. See comment in routine
setBasicGeometry(mesh);
// Global matching tolerance
const scalar tolDim = getMergeDistance
(
......@@ -2736,6 +2765,8 @@ int main(int argc, char *argv[])
),
true // read on master only
);
setBasicGeometry(baseMeshPtr());
Info<< "Reading local, decomposed mesh" << endl;
autoPtr<fvMesh> meshPtr = loadOrCreateMesh
......
......@@ -438,6 +438,20 @@ Foam::IOobject::IOobject
{}
Foam::IOobject::IOobject
(
const IOobject& io,
readOption ro,
writeOption wo
)
:
IOobject(io)
{
rOpt_ = ro;
wOpt_ = wo;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::objectRegistry& Foam::IOobject::db() const
......
......@@ -340,6 +340,14 @@ public:
const word& name
);
//- Copy construct, resetting io options
IOobject
(
const IOobject& io,
readOption,
writeOption
);
//- Clone
autoPtr<IOobject> clone() const
{
......
......@@ -63,9 +63,11 @@ else
runTime.timeName(),
runTime,
Foam::IOobject::MUST_READ
)
),
false
)
);
meshPtr().init(true); // initialise all (lower levels and current)
}
Foam::fvMesh& mesh = meshPtr();
......@@ -21,5 +21,7 @@ Foam::fvMesh mesh
runTime.timeName(),
runTime,
Foam::IOobject::MUST_READ
)
),
false
);
mesh.init(true); // initialise all (lower levels and current)
......@@ -155,6 +155,50 @@ Foam::solution::solution
}
Foam::solution::solution
(
const objectRegistry& obr,
const fileName& dictName,
const dictionary& dict
)
:
IOdictionary
(
IOobject
(
dictName,
obr.time().system(),
obr,
(
obr.readOpt() == IOobject::MUST_READ
|| obr.readOpt() == IOobject::READ_IF_PRESENT
? IOobject::MUST_READ_IF_MODIFIED
: obr.readOpt()
),
IOobject::NO_WRITE
),
dict
),
cache_(dictionary::null),
caching_(false),
fieldRelaxDict_(dictionary::null),
eqnRelaxDict_(dictionary::null),
fieldRelaxDefault_(0),
eqnRelaxDefault_(0),
solvers_(dictionary::null)
{
if
(
readOpt() == IOobject::MUST_READ
|| readOpt() == IOobject::MUST_READ_IF_MODIFIED
|| (readOpt() == IOobject::READ_IF_PRESENT && headerOk())
)
{
read(solutionDict());
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::solution::upgradeSolverDict
......
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -101,13 +102,22 @@ public:
// Constructors
//- Construct for given objectRegistry and dictionary
//- Construct for given objectRegistry and dictionary name
solution
(
const objectRegistry& obr,
const fileName& dictName
);
//- Construct for given objectRegistry, dictionary name and (optional)
// content (gets used in case of NO_READ or dictionary cannot be read)
solution
(
const objectRegistry& obr,
const fileName& dictName,
const dictionary& dict
);
// Member Functions
......
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2015 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -54,6 +55,24 @@ Foam::data::data(const objectRegistry& obr)
}
Foam::data::data(const objectRegistry& obr, const dictionary& dict)
:
IOdictionary
(
IOobject
(
"data",
obr.time().system(),
obr,
IOobject::NO_READ,
IOobject::NO_WRITE
),
dict
),
prevTimeIndex_(0)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::dictionary& Foam::data::solverPerformanceDict() const
......
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -82,6 +83,9 @@ public:
//- Construct for objectRegistry
data(const objectRegistry& obr);
//- Construct for objectRegistry and initial contents
data(const objectRegistry& obr, const dictionary& dict);
// Member Functions
......
......@@ -70,20 +70,29 @@ void Foam::polyMesh::calcDirections() const
forAll(boundaryMesh(), patchi)
{
if (boundaryMesh()[patchi].size())
const polyPatch& pp = boundaryMesh()[patchi];
if (isA<emptyPolyPatch>(pp))
{
if (isA<emptyPolyPatch>(boundaryMesh()[patchi]))
// Force calculation of geometric properties, independent of
// size. This avoids parallel synchronisation problems.
const vectorField::subField fa(pp.faceAreas());
if (pp.size())
{
nEmptyPatches++;
emptyDirVec += sum(cmptMag(boundaryMesh()[patchi].faceAreas()));
emptyDirVec += sum(cmptMag(fa));
}
else if (isA<wedgePolyPatch>(boundaryMesh()[patchi]))
{
const wedgePolyPatch& wpp = refCast<const wedgePolyPatch>
(
boundaryMesh()[patchi]
);
}
else if (isA<wedgePolyPatch>(pp))
{
const wedgePolyPatch& wpp = refCast<const wedgePolyPatch>(pp);
// Force calculation of geometric properties, independent of
// size. This avoids parallel synchronisation problems.
(void)wpp.faceNormals();
if (pp.size())
{
nWedgePatches++;
wedgeDirVec += cmptMag(wpp.centreNormal());
}
......@@ -161,7 +170,7 @@ Foam::autoPtr<Foam::labelIOList> Foam::polyMesh::readTetBasePtIs() const
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::polyMesh::polyMesh(const IOobject& io)
Foam::polyMesh::polyMesh(const IOobject& io, const bool doInit)
:
objectRegistry(io),
primitiveMesh(),
......@@ -328,12 +337,6 @@ Foam::polyMesh::polyMesh(const IOobject& io)
neighbour_.write();
}
// Calculate topology for the patches (processor-processor comms etc.)
boundary_.updateMesh();
// Calculate the geometry for the patches (transformation tensors etc.)
boundary_.calcGeometry();
// Warn if global empty mesh
if (returnReduce(boundary_.empty(), orOp<bool>()))
{
......@@ -352,8 +355,30 @@ Foam::polyMesh::polyMesh(const IOobject& io)
}
}
if (doInit)
{
polyMesh::init(false); // do not init lower levels
}
}
bool Foam::polyMesh::init(const bool doInit)
{
if (doInit)
{
primitiveMesh::init(doInit);
}
// Calculate topology for the patches (processor-processor comms etc.)
boundary_.updateMesh();
// Calculate the geometry for the patches (transformation tensors etc.)
boundary_.calcGeometry();
// Initialise demand-driven data
calcDirections();
return false;
}
......@@ -377,7 +402,7 @@ Foam::polyMesh::polyMesh
instance(),
meshSubDir,
*this,
io.readOpt(),
IOobject::NO_READ, //io.readOpt(),
io.writeOpt()
),
std::move(points)
......@@ -390,7 +415,7 @@ Foam::polyMesh::polyMesh
instance(),
meshSubDir,
*this,
io.readOpt(),
IOobject::NO_READ, //io.readOpt(),
io.writeOpt()
),
std::move(faces)
......@@ -403,7 +428,7 @@ Foam::polyMesh::polyMesh
instance(),
meshSubDir,
*this,
io.readOpt(),
IOobject::NO_READ, //io.readOpt(),
io.writeOpt()
),
std::move(owner)
......@@ -416,7 +441,7 @@ Foam::polyMesh::polyMesh
instance(),
meshSubDir,
*this,
io.readOpt(),
IOobject::NO_READ, //io.readOpt(),
io.writeOpt()
),
std::move(neighbour)
......@@ -430,7 +455,7 @@ Foam::polyMesh::polyMesh
instance(),
meshSubDir,
*this,
io.readOpt(),
IOobject::NO_READ, //io.readOpt(),
io.writeOpt()
),
*this,
......@@ -440,7 +465,7 @@ Foam::polyMesh::polyMesh
comm_(UPstream::worldComm),
geometricD_(Zero),
solutionD_(Zero),
tetBasePtIsPtr_(readTetBasePtIs()),
tetBasePtIsPtr_(nullptr),
cellTreePtr_(nullptr),
pointZones_
(
......@@ -450,7 +475,7 @@ Foam::polyMesh::polyMesh
instance(),
meshSubDir,
*this,
io.readOpt(),
IOobject::NO_READ, //io.readOpt(),
IOobject::NO_WRITE
),
*this,
......@@ -464,7 +489,7 @@ Foam::polyMesh::polyMesh
instance(),
meshSubDir,
*this,
io.readOpt(),
IOobject::NO_READ, //io.readOpt(),
IOobject::NO_WRITE
),
*this,
......@@ -478,7 +503,7 @@ Foam::polyMesh::polyMesh
instance(),
meshSubDir,
*this,
io.readOpt(),
IOobject::NO_READ, //io.readOpt(),
IOobject::NO_WRITE
),
*this,
......@@ -591,7 +616,7 @@ Foam::polyMesh::polyMesh
comm_(UPstream::worldComm),
geometricD_(Zero),
solutionD_(Zero),
tetBasePtIsPtr_(readTetBasePtIs()),
tetBasePtIsPtr_(nullptr),
cellTreePtr_(nullptr),
pointZones_
(
......
......@@ -319,11 +319,11 @@ public:
// Constructors
//- Read construct from IOobject
explicit polyMesh(const IOobject& io);
polyMesh(const IOobject& io, const bool doInit = true);
//- Construct from IOobject or as zero-sized mesh
// Boundary is added using addPatches() member function
polyMesh(const IOobject& io, const zero, bool syncPar=true);
polyMesh(const IOobject& io, const zero, const bool syncPar=true);
//- Construct from IOobject or from components.
// Boundary is added using addPatches() member function
......@@ -593,6 +593,9 @@ public:
const List<cellZone*>& cz
);
//- Initialise all non-demand-driven data
virtual bool init(const bool doInit);
//- Update the mesh based on the mesh files saved in
// time directories
virtual readUpdateState readUpdate();
......
......@@ -753,13 +753,26 @@ bool Foam::polyMesh::checkMeshMotion
vectorField fCtrs(nFaces());
vectorField fAreas(nFaces());
makeFaceCentresAndAreas(newPoints, fCtrs, fAreas);
primitiveMeshTools::makeFaceCentresAndAreas
(
*this,
newPoints,
fCtrs,
fAreas
);
// Check cell volumes and calculate new cell centres
vectorField cellCtrs(nCells());
scalarField cellVols(nCells());
makeCellCentresAndVols(fCtrs, fAreas, cellCtrs, cellVols);
primitiveMeshTools::makeCellCentresAndVols
(
*this,
fCtrs,
fAreas,
cellCtrs,
cellVols
);
// Check cell volumes
bool error = checkCellVolumes
......
......@@ -281,6 +281,44 @@ void Foam::primitiveMesh::reset
}
void Foam::primitiveMesh::resetGeometry
(
pointField&& faceCentres,
pointField&& faceAreas,
pointField&& cellCentres,
scalarField&& cellVolumes
)
{
if
(
faceCentres.size() != nFaces_
|| faceAreas.size() != nFaces_
|| cellCentres.size() != nCells_
|| cellVolumes.size() != nCells_
)
{
FatalErrorInFunction
<< "Wrong sizes of passed in face/cell data"
<< abort(FatalError);
}
// Remove old geometry
clearGeom();
faceCentresPtr_ = new pointField(std::move(faceCentres));
faceAreasPtr_ = new pointField(std::move(faceAreas));
cellCentresPtr_ = new pointField(std::move(cellCentres));
cellVolumesPtr_ = new scalarField(std::move(cellVolumes));
if (debug)
{