From 46dbfabd9dbd3b89796b880c399536129b9b38c7 Mon Sep 17 00:00:00 2001 From: Mattijs Janssens Date: Fri, 11 Dec 2020 10:31:34 +0000 Subject: [PATCH] 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. --- .../extrude/extrudeMesh/extrudeMesh.C | 2 +- .../mesh/manipulation/checkMesh/checkMesh.C | 1 + .../mesh/manipulation/checkMesh/writeFields.C | 27 + .../splitMeshRegions/splitMeshRegions.C | 2 +- .../redistributePar/loadOrCreateMesh.C | 14 + .../redistributePar/redistributePar.C | 31 + src/OpenFOAM/db/IOobject/IOobject.C | 14 + src/OpenFOAM/db/IOobject/IOobject.H | 8 + src/OpenFOAM/include/createMesh.H | 4 +- src/OpenFOAM/include/createNamedMesh.H | 4 +- src/OpenFOAM/matrices/solution/solution.C | 44 + src/OpenFOAM/matrices/solution/solution.H | 12 +- src/OpenFOAM/meshes/data/data.C | 19 + src/OpenFOAM/meshes/data/data.H | 4 + src/OpenFOAM/meshes/polyMesh/polyMesh.C | 77 +- src/OpenFOAM/meshes/polyMesh/polyMesh.H | 7 +- .../polyMesh/polyMeshCheck/polyMeshCheck.C | 17 +- .../meshes/primitiveMesh/primitiveMesh.C | 60 + .../meshes/primitiveMesh/primitiveMesh.H | 30 +- .../primitiveMeshCellCentresAndVols.C | 115 +- .../primitiveMeshCheck/primitiveMeshTools.C | 171 + .../primitiveMeshCheck/primitiveMeshTools.H | 27 +- .../primitiveMeshFaceCentresAndAreas.C | 77 +- .../dynamicFvMesh/dynamicFvMesh.H | 14 +- .../include/createDynamicFvMesh.H | 5 +- src/dynamicMesh/fvMeshSubset/fvMeshSubset.C | 5 +- .../polyMeshFilter/polyMeshFilter.C | 2 +- src/finiteVolume/Make/files | 14 +- .../finiteVolume/fvSchemes/fvSchemes.C | 136 +- .../finiteVolume/fvSchemes/fvSchemes.H | 6 + .../finiteVolume/fvSolution/fvSolution.H | 7 + .../averageNeighbourFvGeometryScheme.C | 880 ++++ .../averageNeighbourFvGeometryScheme.H | 178 + .../basic/basicFvGeometryScheme.C | 388 ++ .../basic/basicFvGeometryScheme.H | 107 + .../fvGeometryScheme/fvGeometryScheme.C | 80 + .../fvGeometryScheme/fvGeometryScheme.H | 164 + .../highAspectRatio/cellAspectRatio.C | 123 + .../highAspectRatio/cellAspectRatio.H | 97 + .../highAspectRatioFvGeometryScheme.C | 466 +++ .../highAspectRatioFvGeometryScheme.H | 137 + .../stabilised/stabilisedFvGeometryScheme.C | 210 + .../stabilised/stabilisedFvGeometryScheme.H | 126 + src/finiteVolume/fvMesh/fvMesh.C | 141 +- src/finiteVolume/fvMesh/fvMesh.H | 39 +- .../fvMesh/fvPatches/fvPatch/fvPatch.H | 16 +- .../surfaceInterpolation.C | 376 +- .../surfaceInterpolation.H | 54 +- .../layerParameters/layerParameters.C | 808 +++- .../layerParameters/layerParameters.H | 113 +- .../snappyHexMeshDriver/snappyLayerDriver.C | 325 +- .../snappyHexMeshDriver/snappyLayerDriver.H | 3 + .../snappyHexMeshDriver/snappyRefineDriver.C | 267 ++ .../snappyHexMeshDriver/snappyRefineDriver.H | 15 + src/meshTools/polyTopoChange/polyTopoChange.H | 3 +- .../polyTopoChange/polyTopoChangeTemplates.C | 8 +- .../snappyHexMesh/airfoilWithLayers/Allclean | 9 + .../snappyHexMesh/airfoilWithLayers/Allrun | 22 + .../constant/transportProperties | 32 + .../constant/triSurface/aerofoil.curvature | Bin 0 -> 5122 bytes .../constant/triSurface/aerofoil.eMesh | Bin 0 -> 17438 bytes .../constant/triSurface/aerofoil.stl | 3586 +++++++++++++++++ .../constant/triSurface/internalFace.stl | 72 + .../airfoilWithLayers/system/blockMeshDict | 56 + .../airfoilWithLayers/system/controlDict | 55 + .../airfoilWithLayers/system/fvSchemes | 73 + .../airfoilWithLayers/system/fvSolution | 71 + .../airfoilWithLayers/system/meshQualityDict | 78 + .../system/snappyHexMeshDict | 753 ++++ .../system/surfaceFeatureExtractDict | 47 + 70 files changed, 10039 insertions(+), 895 deletions(-) create mode 100644 src/finiteVolume/fvMesh/fvGeometryScheme/averageNeighbour/averageNeighbourFvGeometryScheme.C create mode 100644 src/finiteVolume/fvMesh/fvGeometryScheme/averageNeighbour/averageNeighbourFvGeometryScheme.H create mode 100644 src/finiteVolume/fvMesh/fvGeometryScheme/basic/basicFvGeometryScheme.C create mode 100644 src/finiteVolume/fvMesh/fvGeometryScheme/basic/basicFvGeometryScheme.H create mode 100644 src/finiteVolume/fvMesh/fvGeometryScheme/fvGeometryScheme/fvGeometryScheme.C create mode 100644 src/finiteVolume/fvMesh/fvGeometryScheme/fvGeometryScheme/fvGeometryScheme.H create mode 100644 src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/cellAspectRatio.C create mode 100644 src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/cellAspectRatio.H create mode 100644 src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.C create mode 100644 src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.H create mode 100644 src/finiteVolume/fvMesh/fvGeometryScheme/stabilised/stabilisedFvGeometryScheme.C create mode 100644 src/finiteVolume/fvMesh/fvGeometryScheme/stabilised/stabilisedFvGeometryScheme.H create mode 100755 tutorials/mesh/snappyHexMesh/airfoilWithLayers/Allclean create mode 100755 tutorials/mesh/snappyHexMesh/airfoilWithLayers/Allrun create mode 100644 tutorials/mesh/snappyHexMesh/airfoilWithLayers/constant/transportProperties create mode 100644 tutorials/mesh/snappyHexMesh/airfoilWithLayers/constant/triSurface/aerofoil.curvature create mode 100644 tutorials/mesh/snappyHexMesh/airfoilWithLayers/constant/triSurface/aerofoil.eMesh create mode 100644 tutorials/mesh/snappyHexMesh/airfoilWithLayers/constant/triSurface/aerofoil.stl create mode 100644 tutorials/mesh/snappyHexMesh/airfoilWithLayers/constant/triSurface/internalFace.stl create mode 100644 tutorials/mesh/snappyHexMesh/airfoilWithLayers/system/blockMeshDict create mode 100644 tutorials/mesh/snappyHexMesh/airfoilWithLayers/system/controlDict create mode 100644 tutorials/mesh/snappyHexMesh/airfoilWithLayers/system/fvSchemes create mode 100644 tutorials/mesh/snappyHexMesh/airfoilWithLayers/system/fvSolution create mode 100644 tutorials/mesh/snappyHexMesh/airfoilWithLayers/system/meshQualityDict create mode 100644 tutorials/mesh/snappyHexMesh/airfoilWithLayers/system/snappyHexMeshDict create mode 100644 tutorials/mesh/snappyHexMesh/airfoilWithLayers/system/surfaceFeatureExtractDict diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C index 7d0909ba5b..8f9c0b9d7b 100644 --- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C +++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C @@ -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 ), diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C b/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C index 836b2515e2..ec1a2314b1 100644 --- a/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C +++ b/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C @@ -154,6 +154,7 @@ int main(int argc, char *argv[]) "cellShapes", "cellVolume", "cellVolumeRatio", + "cellAspectRatio", "minTetVolume", "minPyrVolume", "cellRegion", diff --git a/applications/utilities/mesh/manipulation/checkMesh/writeFields.C b/applications/utilities/mesh/manipulation/checkMesh/writeFields.C index dfc825a817..92e6e73f47 100644 --- a/applications/utilities/mesh/manipulation/checkMesh/writeFields.C +++ b/applications/utilities/mesh/manipulation/checkMesh/writeFields.C @@ -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 // ~~~~~~~~~ diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C index 22060652a6..aa84e3c65a 100644 --- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C +++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C @@ -680,7 +680,7 @@ autoPtr createRegionMesh regionName, mesh.time().timeName(), mesh.time(), - IOobject::NO_READ, + IOobject::READ_IF_PRESENT, // read fv* if present IOobject::AUTO_WRITE ), mesh diff --git a/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.C b/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.C index dbcec2a27a..5983fc8538 100644 --- a/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.C +++ b/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.C @@ -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::loadOrCreateMesh auto meshPtr = autoPtr::New(io); fvMesh& mesh = *meshPtr; + // Make sure to use a non-parallel geometry calculation method + { + tmp basicGeometry + ( + fvGeometryScheme::New + ( + mesh, + dictionary(), + basicFvGeometryScheme::typeName + ) + ); + mesh.geometry(basicGeometry); + } // Sync patches diff --git a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C index 1d2f5a5a11..32e8ea53e4 100644 --- a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C +++ b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C @@ -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 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 meshPtr = loadOrCreateMesh diff --git a/src/OpenFOAM/db/IOobject/IOobject.C b/src/OpenFOAM/db/IOobject/IOobject.C index f3062e9be0..e25418fef3 100644 --- a/src/OpenFOAM/db/IOobject/IOobject.C +++ b/src/OpenFOAM/db/IOobject/IOobject.C @@ -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 diff --git a/src/OpenFOAM/db/IOobject/IOobject.H b/src/OpenFOAM/db/IOobject/IOobject.H index d6cfba3ba7..a4ed735066 100644 --- a/src/OpenFOAM/db/IOobject/IOobject.H +++ b/src/OpenFOAM/db/IOobject/IOobject.H @@ -340,6 +340,14 @@ public: const word& name ); + //- Copy construct, resetting io options + IOobject + ( + const IOobject& io, + readOption, + writeOption + ); + //- Clone autoPtr clone() const { diff --git a/src/OpenFOAM/include/createMesh.H b/src/OpenFOAM/include/createMesh.H index 8d53556a9c..a4452382f5 100644 --- a/src/OpenFOAM/include/createMesh.H +++ b/src/OpenFOAM/include/createMesh.H @@ -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(); diff --git a/src/OpenFOAM/include/createNamedMesh.H b/src/OpenFOAM/include/createNamedMesh.H index 707effa653..a66327b470 100644 --- a/src/OpenFOAM/include/createNamedMesh.H +++ b/src/OpenFOAM/include/createNamedMesh.H @@ -21,5 +21,7 @@ Foam::fvMesh mesh runTime.timeName(), runTime, Foam::IOobject::MUST_READ - ) + ), + false ); +mesh.init(true); // initialise all (lower levels and current) diff --git a/src/OpenFOAM/matrices/solution/solution.C b/src/OpenFOAM/matrices/solution/solution.C index b38406b32d..498503ff47 100644 --- a/src/OpenFOAM/matrices/solution/solution.C +++ b/src/OpenFOAM/matrices/solution/solution.C @@ -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 diff --git a/src/OpenFOAM/matrices/solution/solution.H b/src/OpenFOAM/matrices/solution/solution.H index 88da0ece22..793cd3520e 100644 --- a/src/OpenFOAM/matrices/solution/solution.H +++ b/src/OpenFOAM/matrices/solution/solution.H @@ -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 diff --git a/src/OpenFOAM/meshes/data/data.C b/src/OpenFOAM/meshes/data/data.C index 95712bb7d0..15e419193f 100644 --- a/src/OpenFOAM/meshes/data/data.C +++ b/src/OpenFOAM/meshes/data/data.C @@ -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 diff --git a/src/OpenFOAM/meshes/data/data.H b/src/OpenFOAM/meshes/data/data.H index 7c5c95951f..91da9a93d4 100644 --- a/src/OpenFOAM/meshes/data/data.H +++ b/src/OpenFOAM/meshes/data/data.H @@ -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 diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C index 81fcd21f7f..232cb3de79 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C @@ -70,20 +70,29 @@ void Foam::polyMesh::calcDirections() const forAll(boundaryMesh(), patchi) { - if (boundaryMesh()[patchi].size()) + const polyPatch& pp = boundaryMesh()[patchi]; + if (isA(pp)) { - if (isA(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(boundaryMesh()[patchi])) - { - const wedgePolyPatch& wpp = refCast - ( - boundaryMesh()[patchi] - ); + } + else if (isA(pp)) + { + const wedgePolyPatch& wpp = refCast(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::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())) { @@ -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_ ( diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.H b/src/OpenFOAM/meshes/polyMesh/polyMesh.H index 5441e564ca..6851fb68eb 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMesh.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.H @@ -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& 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(); diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C index 34d9a2f517..04fc7a63b7 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C @@ -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 diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C index 027d3d6825..2f3cb10efc 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C @@ -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) + { + Pout<< "primitiveMesh::resetGeometry : geometry reset for" + << " nFaces:" << faceCentresPtr_->size() + << " nCells:" << cellCentresPtr_->size() << endl; + } +} + + Foam::tmp Foam::primitiveMesh::movePoints ( const pointField& newPoints, @@ -324,4 +362,26 @@ const Foam::cellShapeList& Foam::primitiveMesh::cellShapes() const } + +void Foam::primitiveMesh::updateGeom() +{ + if (!faceCentresPtr_) + { + calcFaceCentresAndAreas(); + } + if (!faceAreasPtr_) + { + calcFaceCentresAndAreas(); + } + if (!cellCentresPtr_) + { + calcCellCentresAndVols(); + } + if (!cellVolumesPtr_) + { + calcCellCentresAndVols(); + } +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H index d7cc8bf323..3ef4a22261 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H @@ -255,22 +255,9 @@ protected: //- Calculate face centres and areas void calcFaceCentresAndAreas() const; - void makeFaceCentresAndAreas - ( - const pointField& p, - vectorField& fCtrs, - vectorField& fAreas - ) const; //- Calculate cell centres and volumes void calcCellCentresAndVols() const; - void makeCellCentresAndVols - ( - const vectorField& fCtrs, - const vectorField& fAreas, - vectorField& cellCtrs, - scalarField& cellVols - ) const; //- Calculate edge vectors void calcEdgeVectors() const; @@ -476,6 +463,21 @@ public: cellList& cells ); + //- Reset the local geometry + void resetGeometry + ( + pointField&& faceCentres, + pointField&& faceAreas, + pointField&& cellCentres, + scalarField&& cellVolumes + ); + + //- Initialise all non-demand-driven data + virtual bool init(const bool doInit) + { + return false; + } + // Access @@ -891,6 +893,8 @@ public: const labelList& cellEdges(const label celli) const; + //- Update all geometric data + virtual void updateGeom(); //- Clear geometry void clearGeom(); diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C index 5f0ff8bcef..f8a1e32340 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C @@ -30,7 +30,7 @@ Description \*---------------------------------------------------------------------------*/ #include "primitiveMesh.H" -#include "PrecisionAdaptor.H" +#include "primitiveMeshTools.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // @@ -61,7 +61,14 @@ void Foam::primitiveMesh::calcCellCentresAndVols() const scalarField& cellVols = *cellVolumesPtr_; // Make centres and volumes - makeCellCentresAndVols(faceCentres(), faceAreas(), cellCtrs, cellVols); + primitiveMeshTools::makeCellCentresAndVols + ( + *this, + faceCentres(), + faceAreas(), + cellCtrs, + cellVols + ); if (debug) { @@ -72,111 +79,14 @@ void Foam::primitiveMesh::calcCellCentresAndVols() const } -void Foam::primitiveMesh::makeCellCentresAndVols -( - const vectorField& fCtrs, - const vectorField& fAreas, - vectorField& cellCtrs_s, - scalarField& cellVols_s -) const -{ - typedef Vector solveVector; - - PrecisionAdaptor tcellCtrs(cellCtrs_s, false); - Field& cellCtrs = tcellCtrs.ref(); - PrecisionAdaptor tcellVols(cellVols_s, false); - Field& cellVols = tcellVols.ref(); - - cellCtrs = Zero; - cellVols = 0.0; - - const labelList& own = faceOwner(); - const labelList& nei = faceNeighbour(); - - // first estimate the approximate cell centre as the average of - // face centres - - Field cEst(nCells(), Zero); - labelField nCellFaces(nCells(), Zero); - - forAll(own, facei) - { - cEst[own[facei]] += solveVector(fCtrs[facei]); - ++nCellFaces[own[facei]]; - } - - forAll(nei, facei) - { - cEst[nei[facei]] += solveVector(fCtrs[facei]); - ++nCellFaces[nei[facei]]; - } - - forAll(cEst, celli) - { - cEst[celli] /= nCellFaces[celli]; - } - - forAll(own, facei) - { - const solveVector fc(fCtrs[facei]); - const solveVector fA(fAreas[facei]); - - // Calculate 3*face-pyramid volume - solveScalar pyr3Vol = - fA & (fc - cEst[own[facei]]); - - // Calculate face-pyramid centre - solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[own[facei]]; - - // Accumulate volume-weighted face-pyramid centre - cellCtrs[own[facei]] += pyr3Vol*pc; - - // Accumulate face-pyramid volume - cellVols[own[facei]] += pyr3Vol; - } - - forAll(nei, facei) - { - const solveVector fc(fCtrs[facei]); - const solveVector fA(fAreas[facei]); - - // Calculate 3*face-pyramid volume - solveScalar pyr3Vol = - fA & (cEst[nei[facei]] - fc); - - // Calculate face-pyramid centre - solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[nei[facei]]; - - // Accumulate volume-weighted face-pyramid centre - cellCtrs[nei[facei]] += pyr3Vol*pc; - - // Accumulate face-pyramid volume - cellVols[nei[facei]] += pyr3Vol; - } - - forAll(cellCtrs, celli) - { - if (mag(cellVols[celli]) > VSMALL) - { - cellCtrs[celli] /= cellVols[celli]; - } - else - { - cellCtrs[celli] = cEst[celli]; - } - } - - cellVols *= (1.0/3.0); -} - - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // const Foam::vectorField& Foam::primitiveMesh::cellCentres() const { if (!cellCentresPtr_) { - calcCellCentresAndVols(); + //calcCellCentresAndVols(); + const_cast(*this).updateGeom(); } return *cellCentresPtr_; @@ -187,7 +97,8 @@ const Foam::scalarField& Foam::primitiveMesh::cellVolumes() const { if (!cellVolumesPtr_) { - calcCellCentresAndVols(); + //calcCellCentresAndVols(); + const_cast(*this).updateGeom(); } return *cellVolumesPtr_; diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C index 98a8b032b7..46b62d3466 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C @@ -27,11 +27,182 @@ License \*---------------------------------------------------------------------------*/ #include "primitiveMeshTools.H" +#include "primitiveMesh.H" #include "syncTools.H" #include "pyramidPointFaceRef.H" +#include "PrecisionAdaptor.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +void Foam::primitiveMeshTools::makeFaceCentresAndAreas +( + const primitiveMesh& mesh, + const pointField& p, + vectorField& fCtrs, + vectorField& fAreas +) +{ + const faceList& fs = mesh.faces(); + + forAll(fs, facei) + { + const labelList& f = fs[facei]; + const label nPoints = f.size(); + + // If the face is a triangle, do a direct calculation for efficiency + // and to avoid round-off error-related problems + if (nPoints == 3) + { + fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]); + fAreas[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]])); + } + else + { + typedef Vector solveVector; + + solveVector sumN = Zero; + solveScalar sumA = 0.0; + solveVector sumAc = Zero; + + solveVector fCentre = p[f[0]]; + for (label pi = 1; pi < nPoints; pi++) + { + fCentre += solveVector(p[f[pi]]); + } + + fCentre /= nPoints; + + for (label pi = 0; pi < nPoints; pi++) + { + const label nextPi(pi == nPoints-1 ? 0 : pi+1); + const solveVector nextPoint(p[f[nextPi]]); + const solveVector thisPoint(p[f[pi]]); + + solveVector c = thisPoint + nextPoint + fCentre; + solveVector n = (nextPoint - thisPoint)^(fCentre - thisPoint); + solveScalar a = mag(n); + sumN += n; + sumA += a; + sumAc += a*c; + } + + // This is to deal with zero-area faces. Mark very small faces + // to be detected in e.g., processorPolyPatch. + if (sumA < ROOTVSMALL) + { + fCtrs[facei] = fCentre; + fAreas[facei] = Zero; + } + else + { + fCtrs[facei] = (1.0/3.0)*sumAc/sumA; + fAreas[facei] = 0.5*sumN; + } + } + } +} + + +void Foam::primitiveMeshTools::makeCellCentresAndVols +( + const primitiveMesh& mesh, + const vectorField& fCtrs, + const vectorField& fAreas, + vectorField& cellCtrs_s, + scalarField& cellVols_s +) +{ + typedef Vector solveVector; + + PrecisionAdaptor tcellCtrs(cellCtrs_s); + Field& cellCtrs = tcellCtrs.ref(); + PrecisionAdaptor tcellVols(cellVols_s); + Field& cellVols = tcellVols.ref(); + + // Clear the fields for accumulation + cellCtrs = Zero; + cellVols = 0.0; + + const labelList& own = mesh.faceOwner(); + const labelList& nei = mesh.faceNeighbour(); + + // first estimate the approximate cell centre as the average of + // face centres + + Field cEst(mesh.nCells(), Zero); + labelField nCellFaces(mesh.nCells(), Zero); + + forAll(own, facei) + { + cEst[own[facei]] += solveVector(fCtrs[facei]); + ++nCellFaces[own[facei]]; + } + + forAll(nei, facei) + { + cEst[nei[facei]] += solveVector(fCtrs[facei]); + ++nCellFaces[nei[facei]]; + } + + forAll(cEst, celli) + { + cEst[celli] /= nCellFaces[celli]; + } + + forAll(own, facei) + { + const solveVector fc(fCtrs[facei]); + const solveVector fA(fAreas[facei]); + + // Calculate 3*face-pyramid volume + solveScalar pyr3Vol = + fA & (fc - cEst[own[facei]]); + + // Calculate face-pyramid centre + solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[own[facei]]; + + // Accumulate volume-weighted face-pyramid centre + cellCtrs[own[facei]] += pyr3Vol*pc; + + // Accumulate face-pyramid volume + cellVols[own[facei]] += pyr3Vol; + } + + forAll(nei, facei) + { + const solveVector fc(fCtrs[facei]); + const solveVector fA(fAreas[facei]); + + // Calculate 3*face-pyramid volume + solveScalar pyr3Vol = + fA & (cEst[nei[facei]] - fc); + + // Calculate face-pyramid centre + solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[nei[facei]]; + + // Accumulate volume-weighted face-pyramid centre + cellCtrs[nei[facei]] += pyr3Vol*pc; + + // Accumulate face-pyramid volume + cellVols[nei[facei]] += pyr3Vol; + } + + forAll(cellCtrs, celli) + { + if (mag(cellVols[celli]) > VSMALL) + { + cellCtrs[celli] /= cellVols[celli]; + } + else + { + cellCtrs[celli] = cEst[celli]; + } + } + + cellVols *= (1.0/3.0); +} + + Foam::scalar Foam::primitiveMeshTools::faceSkewness ( const primitiveMesh& mesh, diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H index 9803acac36..3bfef13593 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H @@ -38,22 +38,45 @@ SourceFiles #ifndef primitiveMeshTools_H #define primitiveMeshTools_H -#include "primitiveMesh.H" #include "bitSet.H" +#include "primitiveFields.H" +#include "pointField.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { +// Forward declarations +class primitiveMesh; + /*---------------------------------------------------------------------------*\ - Namespace primitiveMeshTools Declaration + Namespace primitiveMeshTools Declaration \*---------------------------------------------------------------------------*/ class primitiveMeshTools { public: + //- Calculate face centres and areas + static void makeFaceCentresAndAreas + ( + const primitiveMesh& mesh, + const pointField& p, + vectorField& fCtrs, + vectorField& fAreas + ); + + //- Calculate cell centres and volumes from face properties + static void makeCellCentresAndVols + ( + const primitiveMesh& mesh, + const vectorField& fCtrs, + const vectorField& fAreas, + vectorField& cellCtrs, + scalarField& cellVols + ); + //- Generate non-orthogonality field (internal faces only) static tmp faceOrthogonality ( diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C index defcfae773..dda3e3634c 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C @@ -33,6 +33,7 @@ Description \*---------------------------------------------------------------------------*/ #include "primitiveMesh.H" +#include "primitiveMeshTools.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // @@ -60,7 +61,7 @@ void Foam::primitiveMesh::calcFaceCentresAndAreas() const faceAreasPtr_ = new vectorField(nFaces()); vectorField& fAreas = *faceAreasPtr_; - makeFaceCentresAndAreas(points(), fCtrs, fAreas); + primitiveMeshTools::makeFaceCentresAndAreas(*this, points(), fCtrs, fAreas); if (debug) { @@ -71,81 +72,14 @@ void Foam::primitiveMesh::calcFaceCentresAndAreas() const } -void Foam::primitiveMesh::makeFaceCentresAndAreas -( - const pointField& p, - vectorField& fCtrs, - vectorField& fAreas -) const -{ - const faceList& fs = faces(); - - forAll(fs, facei) - { - const labelList& f = fs[facei]; - const label nPoints = f.size(); - - // If the face is a triangle, do a direct calculation for efficiency - // and to avoid round-off error-related problems - if (nPoints == 3) - { - fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]); - fAreas[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]])); - } - else - { - typedef Vector solveVector; - - solveVector sumN = Zero; - solveScalar sumA = 0.0; - solveVector sumAc = Zero; - - solveVector fCentre = p[f[0]]; - for (label pi = 1; pi < nPoints; pi++) - { - fCentre += solveVector(p[f[pi]]); - } - - fCentre /= nPoints; - - for (label pi = 0; pi < nPoints; pi++) - { - const label nextPi(pi == nPoints-1 ? 0 : pi+1); - const solveVector nextPoint(p[f[nextPi]]); - const solveVector thisPoint(p[f[pi]]); - - solveVector c = thisPoint + nextPoint + fCentre; - solveVector n = (nextPoint - thisPoint)^(fCentre - thisPoint); - solveScalar a = mag(n); - sumN += n; - sumA += a; - sumAc += a*c; - } - - // This is to deal with zero-area faces. Mark very small faces - // to be detected in e.g., processorPolyPatch. - if (sumA < ROOTVSMALL) - { - fCtrs[facei] = fCentre; - fAreas[facei] = Zero; - } - else - { - fCtrs[facei] = (1.0/3.0)*sumAc/sumA; - fAreas[facei] = 0.5*sumN; - } - } - } -} - - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // const Foam::vectorField& Foam::primitiveMesh::faceCentres() const { if (!faceCentresPtr_) { - calcFaceCentresAndAreas(); + //calcFaceCentresAndAreas(); + const_cast(*this).updateGeom(); } return *faceCentresPtr_; @@ -156,7 +90,8 @@ const Foam::vectorField& Foam::primitiveMesh::faceAreas() const { if (!faceAreasPtr_) { - calcFaceCentresAndAreas(); + //calcFaceCentresAndAreas(); + const_cast(*this).updateGeom(); } return *faceAreasPtr_; diff --git a/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMesh.H b/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMesh.H index 8a0a7e3563..c59bdbd23c 100644 --- a/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMesh.H +++ b/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMesh.H @@ -104,7 +104,7 @@ public: // Constructors //- Construct from an IOobject - explicit dynamicFvMesh(const IOobject& io); + explicit dynamicFvMesh(const IOobject& io); //, const bool doInit=true); //- Construct from components without boundary. // Boundary is added using addFvPatches() member function @@ -139,7 +139,11 @@ public: //- Select, construct and return the dynamicFvMesh // If the constant/dynamicMeshDict does not exist // a staticFvMesh is returned - static autoPtr New(const IOobject& io); + static autoPtr New + ( + const IOobject& io //, + //const bool doInit=true + ); //- Select, construct and return the dynamicFvMesh @@ -148,7 +152,8 @@ public: static autoPtr New ( const argList& args, - const Time& runTime + const Time& runTime //, + //const bool doInit=true ); @@ -158,6 +163,9 @@ public: // Member Functions + //- Initialise all non-demand-driven data + //virtual bool init(const bool doInit); + //- Is mesh dynamic virtual bool dynamic() const { diff --git a/src/dynamicFvMesh/include/createDynamicFvMesh.H b/src/dynamicFvMesh/include/createDynamicFvMesh.H index 543bcbfe34..d447d12dd6 100644 --- a/src/dynamicFvMesh/include/createDynamicFvMesh.H +++ b/src/dynamicFvMesh/include/createDynamicFvMesh.H @@ -1,6 +1,9 @@ Info<< "Create mesh for time = " << runTime.timeName() << nl << endl; -autoPtr meshPtr(dynamicFvMesh::New(args, runTime)); +autoPtr meshPtr(dynamicFvMesh::New(args, runTime));//, false)); dynamicFvMesh& mesh = meshPtr(); + +// initialise all (lower levels and current) +mesh.init(true); diff --git a/src/dynamicMesh/fvMeshSubset/fvMeshSubset.C b/src/dynamicMesh/fvMeshSubset/fvMeshSubset.C index 98a5fac060..416442f4fb 100644 --- a/src/dynamicMesh/fvMeshSubset/fvMeshSubset.C +++ b/src/dynamicMesh/fvMeshSubset/fvMeshSubset.C @@ -268,7 +268,7 @@ void Foam::fvMeshSubset::removeCellsImpl baseMesh().name(), baseMesh().time().timeName(), baseMesh().time(), - IOobject::NO_READ, + IOobject::READ_IF_PRESENT, // read fv* if present IOobject::NO_WRITE ), baseMesh(), @@ -980,9 +980,10 @@ void Foam::fvMeshSubset::setCellSubset baseMesh().name(), baseMesh().time().timeName(), baseMesh().time(), - IOobject::NO_READ, + IOobject::NO_READ, // do not read any dictionaries IOobject::NO_WRITE ), + baseMesh(), // get dictionaries from base mesh std::move(newPoints), std::move(newFaces), std::move(newCells), diff --git a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C index d563d0e5e9..3f2e76dc4f 100644 --- a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C +++ b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C @@ -79,7 +79,7 @@ Foam::autoPtr Foam::polyMeshFilter::copyMesh(const fvMesh& mesh) mesh.name(), mesh.polyMesh::instance(), mesh.time(), - IOobject::NO_READ, + IOobject::READ_IF_PRESENT, // read fv* if present IOobject::NO_WRITE, false ), diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 1f175aafaa..e2f15cf369 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -1,6 +1,17 @@ fvMesh/fvMeshGeometry.C fvMesh/fvMesh.C +fvGeometryScheme = fvMesh/fvGeometryScheme +$(fvGeometryScheme)/fvGeometryScheme/fvGeometryScheme.C +$(fvGeometryScheme)/basic/basicFvGeometryScheme.C +$(fvGeometryScheme)/highAspectRatio/highAspectRatioFvGeometryScheme.C +$(fvGeometryScheme)/highAspectRatio/cellAspectRatio.C +$(fvGeometryScheme)/averageNeighbour/averageNeighbourFvGeometryScheme.C +$(fvGeometryScheme)/stabilised/stabilisedFvGeometryScheme.C + +surfaceInterpolation = interpolation/surfaceInterpolation +$(surfaceInterpolation)/surfaceInterpolation/surfaceInterpolation.C + fvMesh/singleCellFvMesh/singleCellFvMesh.C fvMesh/simplifiedFvMesh/simplifiedFvMesh/simplifiedFvMesh.C @@ -53,7 +64,6 @@ fvMeshMapper = fvMesh/fvMeshMapper $(fvMeshMapper)/fvPatchMapper.C $(fvMeshMapper)/fvSurfaceMapper.C - extendedStencil = fvMesh/extendedStencil cellToCell = $(extendedStencil)/cellToCell @@ -316,8 +326,6 @@ volPointInterpolation = interpolation/volPointInterpolation $(volPointInterpolation)/volPointInterpolation.C $(volPointInterpolation)/pointConstraints.C -surfaceInterpolation = interpolation/surfaceInterpolation -$(surfaceInterpolation)/surfaceInterpolation/surfaceInterpolation.C $(surfaceInterpolation)/surfaceInterpolationScheme/surfaceInterpolationSchemes.C $(surfaceInterpolation)/blendedSchemeBase/blendedSchemeBaseName.C diff --git a/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C b/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C index f65ae4c1b1..b67b65de6b 100644 --- a/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C +++ b/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2015 OpenFOAM Foundation - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -327,6 +327,140 @@ Foam::fvSchemes::fvSchemes(const objectRegistry& obr) } +Foam::fvSchemes::fvSchemes(const objectRegistry& obr, const dictionary& dict) +: + IOdictionary + ( + IOobject + ( + "fvSchemes", + 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 + ), + ddtSchemes_ + ( + ITstream + ( + objectPath() + ".ddtSchemes", + tokenList() + )() + ), + defaultDdtScheme_ + ( + ddtSchemes_.name() + ".default", + tokenList() + ), + d2dt2Schemes_ + ( + ITstream + ( + objectPath() + ".d2dt2Schemes", + tokenList() + )() + ), + defaultD2dt2Scheme_ + ( + d2dt2Schemes_.name() + ".default", + tokenList() + ), + interpolationSchemes_ + ( + ITstream + ( + objectPath() + ".interpolationSchemes", + tokenList() + )() + ), + defaultInterpolationScheme_ + ( + interpolationSchemes_.name() + ".default", + tokenList() + ), + divSchemes_ + ( + ITstream + ( + objectPath() + ".divSchemes", + tokenList() + )() + ), + defaultDivScheme_ + ( + divSchemes_.name() + ".default", + tokenList() + ), + gradSchemes_ + ( + ITstream + ( + objectPath() + ".gradSchemes", + tokenList() + )() + ), + defaultGradScheme_ + ( + gradSchemes_.name() + ".default", + tokenList() + ), + snGradSchemes_ + ( + ITstream + ( + objectPath() + ".snGradSchemes", + tokenList() + )() + ), + defaultSnGradScheme_ + ( + snGradSchemes_.name() + ".default", + tokenList() + ), + laplacianSchemes_ + ( + ITstream + ( + objectPath() + ".laplacianSchemes", + tokenList() + )() + ), + defaultLaplacianScheme_ + ( + laplacianSchemes_.name() + ".default", + tokenList() + ), + fluxRequired_ + ( + ITstream + ( + objectPath() + ".fluxRequired", + tokenList() + )() + ), + defaultFluxRequired_(false), + steady_(false) +{ + if + ( + readOpt() == IOobject::MUST_READ + || readOpt() == IOobject::MUST_READ_IF_MODIFIED + || (readOpt() == IOobject::READ_IF_PRESENT && headerOk()) + ) + { + read(schemesDict()); + } +} + + + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // bool Foam::fvSchemes::read() diff --git a/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.H b/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.H index 81f1edecaf..f129c23208 100644 --- a/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.H +++ b/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.H @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2015 OpenFOAM Foundation + Copyright (C) 2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -111,6 +112,11 @@ public: //- Construct for objectRegistry fvSchemes(const objectRegistry& obr); + //- Construct from objectRegistry and supplied (optional) content + // (gets used in case of NO_READ or fvSchemes dictionary cannot be + // read) + fvSchemes(const objectRegistry& obr, const dictionary& dict); + // Member Functions diff --git a/src/finiteVolume/finiteVolume/fvSolution/fvSolution.H b/src/finiteVolume/finiteVolume/fvSolution/fvSolution.H index 7cf1aed1f3..bdafb4c86b 100644 --- a/src/finiteVolume/finiteVolume/fvSolution/fvSolution.H +++ b/src/finiteVolume/finiteVolume/fvSolution/fvSolution.H @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011 OpenFOAM Foundation + Copyright (C) 2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -69,6 +70,12 @@ public: : solution(obr, "fvSolution") {} + + //- Construct for objectRegistry and optional contents + fvSolution(const objectRegistry& obr, const dictionary& dict) + : + solution(obr, "fvSolution", dict) + {} }; diff --git a/src/finiteVolume/fvMesh/fvGeometryScheme/averageNeighbour/averageNeighbourFvGeometryScheme.C b/src/finiteVolume/fvMesh/fvGeometryScheme/averageNeighbour/averageNeighbourFvGeometryScheme.C new file mode 100644 index 0000000000..507f3121d3 --- /dev/null +++ b/src/finiteVolume/fvMesh/fvGeometryScheme/averageNeighbour/averageNeighbourFvGeometryScheme.C @@ -0,0 +1,880 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2020 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "averageNeighbourFvGeometryScheme.H" +#include "addToRunTimeSelectionTable.H" +#include "fvMesh.H" +#include "cellAspectRatio.H" +#include "syncTools.H" +#include "polyMeshTools.H" +#include "unitConversion.H" +#include "OBJstream.H" +#include "surfaceWriter.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(averageNeighbourFvGeometryScheme, 0); + addToRunTimeSelectionTable + ( + fvGeometryScheme, + averageNeighbourFvGeometryScheme, + dict + ); +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::label Foam::averageNeighbourFvGeometryScheme::clipFaceTet +( + const scalar minRatio, + const vectorField& faceCentres, + const vectorField& faceNormals, + + vectorField& faceCorrection +) const +{ + // Clip correction vector if any triangle becomes too small. Return number + // of correction vectors clipped + + typedef Vector solveVector; + + const pointField& p = mesh_.points(); + + label nClipped = 0; + for (label facei = 0; facei < mesh_.nFaces(); facei++) + { + const vector& fcCorr = faceCorrection[facei]; + if (fcCorr != vector::zero) + { + const vector& fn = faceNormals[facei]; + const point& fc = faceCentres[facei]; + const face& f = mesh_.faces()[facei]; + + forAll(f, fp) + { + const solveVector thisPt(p[f[fp]]); + const solveVector nextPt(p[f.fcValue(fp)]); + const solveVector d(nextPt-thisPt); + + // Calculate triangle area with correction + const solveVector nCorr(d^(fc+fcCorr - thisPt)); + + if ((nCorr & fn) < 0) + { + // Triangle points wrong way + faceCorrection[facei] = vector::zero; + nClipped++; + break; + } + else + { + // Calculate triangle area without correction + const solveVector n(d^(fc - thisPt)); + if ((n & fn) < 0) + { + // Original triangle points the wrong way, new one is ok + } + else + { + // Both point correctly. Make sure triangle doesn't get + // too small + if (mag(nCorr) < minRatio*mag(n)) + { + faceCorrection[facei] = vector::zero; + nClipped++; + break; + } + } + } + } + } + } + return returnReduce(nClipped, sumOp