From 3efaa302a96808128a4dc5d96f50baccdc366524 Mon Sep 17 00:00:00 2001 From: mattijs <mattijs@hunt.opencfd.co.uk> Date: Fri, 12 Jun 2009 16:41:13 +0100 Subject: [PATCH] upwind versions of stencils --- etc/controlDict | 1 + .../mapPolyMesh/mapDistribute/mapDistribute.C | 29 +++ .../mapPolyMesh/mapDistribute/mapDistribute.H | 18 +- src/finiteVolume/Make/files | 5 + .../fullStencils/CECCellToCellStencil.C | 3 - .../centredCECCellToFaceStencilObject.H | 9 +- .../centredCFCCellToFaceStencilObject.H | 9 +- .../centredCPCCellToFaceStencilObject.H | 9 +- .../centredFECCellToFaceStencilObject.H | 9 +- .../pureUpwindCFCCellToFaceStencilObject.C | 38 ++++ .../pureUpwindCFCCellToFaceStencilObject.H | 95 +++++++++ .../upwindCECCellToFaceStencilObject.H | 11 +- .../upwindCFCCellToFaceStencilObject.H | 15 +- .../upwindCPCCellToFaceStencilObject.H | 11 +- .../upwindFECCellToFaceStencilObject.H | 11 +- .../cellToFace/extendedCellToFaceStencil.C | 67 ++++++- .../cellToFace/extendedCellToFaceStencil.H | 14 ++ .../extendedCellToFaceStencilTemplates.C | 5 +- .../extendedUpwindCellToFaceStencil.C | 132 ++++++++++++- .../extendedUpwindCellToFaceStencil.H | 29 ++- ...extendedUpwindCellToFaceStencilTemplates.C | 5 +- .../faceToCell/extendedFaceToCellStencil.H | 2 + .../schemes/CentredFitScheme/CentredFitData.C | 2 +- .../schemes/FitData/FitData.C | 44 ++++- .../schemes/FitData/FitData.H | 28 ++- .../PureUpwindFitScheme/PureUpwindFitScheme.H | 186 ++++++++++++++++++ .../schemes/UpwindFitScheme/UpwindFitData.C | 101 +++++++--- .../schemes/UpwindFitScheme/UpwindFitData.H | 5 +- .../schemes/UpwindFitScheme/UpwindFitScheme.H | 2 + .../linearPureUpwindFit/linearPureUpwindFit.C | 49 +++++ .../quadraticLinearPureUpwindFit.C | 45 +++++ .../quadraticLinearUpwindFit.C | 4 +- 32 files changed, 915 insertions(+), 78 deletions(-) create mode 100644 src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/pureUpwindCFCCellToFaceStencilObject.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/pureUpwindCFCCellToFaceStencilObject.H create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/PureUpwindFitScheme/PureUpwindFitScheme.H create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearPureUpwindFit/linearPureUpwindFit.C create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearPureUpwindFit/quadraticLinearPureUpwindFit.C diff --git a/etc/controlDict b/etc/controlDict index c5f3831a30c..d7824f102ed 100644 --- a/etc/controlDict +++ b/etc/controlDict @@ -95,6 +95,7 @@ DebugSwitches Euler 0; EulerImplicit 0; EulerRotation 0; + extendedCellToFaceStencil 0; FDIC 0; FaceCellWave 0; GAMG 0; diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C index 4bd093e6412..cc4278e2218 100644 --- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C @@ -273,6 +273,15 @@ Foam::mapDistribute::mapDistribute } +Foam::mapDistribute::mapDistribute(const mapDistribute& map) +: + constructSize_(map.constructSize_), + subMap_(map.subMap_), + constructMap_(map.constructMap_), + schedulePtr_() +{} + + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // void Foam::mapDistribute::compact(const boolList& elemIsUsed) @@ -413,4 +422,24 @@ void Foam::mapDistribute::compact(const boolList& elemIsUsed) } +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +void Foam::mapDistribute::operator=(const mapDistribute& rhs) +{ + // Check for assignment to self + if (this == &rhs) + { + FatalErrorIn + ( + "Foam::mapDistribute::operator=(const Foam::mapDistribute&)" + ) << "Attempted assignment to self" + << abort(FatalError); + } + constructSize_ = rhs.constructSize_; + subMap_ = rhs.subMap_; + constructMap_ = rhs.constructMap_; + schedulePtr_.clear(); +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H index 3ccf8104769..d2cfe64ca53 100644 --- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H @@ -81,16 +81,6 @@ class mapDistribute mutable autoPtr<List<labelPair> > schedulePtr_; - // Private Member Functions - - //- Disallow default bitwise copy construct - mapDistribute(const mapDistribute&); - - //- Disallow default bitwise assignment - void operator=(const mapDistribute&); - - - public: // Constructors @@ -120,6 +110,9 @@ public: const labelList& recvProcs ); + //- Construct copy + mapDistribute(const mapDistribute&); + // Member Functions @@ -262,6 +255,11 @@ public: "mapDistribute::updateMesh(const mapPolyMesh&)" ); } + + // Member Operators + + void operator=(const mapDistribute&); + }; diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index a438aa25b0a..0367b8acf63 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -65,6 +65,7 @@ $(cellToFace)/MeshObjects/upwindCECCellToFaceStencilObject.C $(cellToFace)/MeshObjects/upwindCFCCellToFaceStencilObject.C $(cellToFace)/MeshObjects/upwindCPCCellToFaceStencilObject.C $(cellToFace)/MeshObjects/upwindFECCellToFaceStencilObject.C +$(cellToFace)/MeshObjects/pureUpwindCFCCellToFaceStencilObject.C faceToCell = $(extendedStencil)/faceToCell $(faceToCell)/fullStencils/faceToCellStencil.C @@ -216,6 +217,10 @@ $(schemes)/quadraticFit/quadraticFit.C $(schemes)/quadraticLinearUpwindFit/quadraticLinearUpwindFit.C $(schemes)/quadraticUpwindFit/quadraticUpwindFit.C $(schemes)/cubicUpwindFit/cubicUpwindFit.C +/* +$(schemes)/quadraticLinearPureUpwindFit/quadraticLinearPureUpwindFit.C +*/ +$(schemes)/linearPureUpwindFit/linearPureUpwindFit.C limitedSchemes = $(surfaceInterpolation)/limitedSchemes $(limitedSchemes)/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationSchemes.C diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/CECCellToCellStencil.C b/src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/CECCellToCellStencil.C index fc167a24391..b6d2dfa91ef 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/CECCellToCellStencil.C +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/CECCellToCellStencil.C @@ -26,9 +26,6 @@ License #include "CECCellToCellStencil.H" #include "syncTools.H" -//#include "meshTools.H" -//#include "OFstream.H" -//#include "Time.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCECCellToFaceStencilObject.H b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCECCellToFaceStencilObject.H index 7f267d643cc..b639ee910aa 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCECCellToFaceStencilObject.H +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCECCellToFaceStencilObject.H @@ -67,7 +67,14 @@ public: : MeshObject<fvMesh, centredCECCellToFaceStencilObject>(mesh), extendedCentredCellToFaceStencil(CECCellToFaceStencil(mesh)) - {} + { + if (extendedCellToFaceStencil::debug) + { + Info<< "Generated centred stencil " << type() + << nl << endl; + writeStencilStats(Info, stencil(), map()); + } + } // Destructor diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCFCCellToFaceStencilObject.H b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCFCCellToFaceStencilObject.H index 6ae971e2a79..d6961220b5f 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCFCCellToFaceStencilObject.H +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCFCCellToFaceStencilObject.H @@ -67,7 +67,14 @@ public: : MeshObject<fvMesh, centredCFCCellToFaceStencilObject>(mesh), extendedCentredCellToFaceStencil(CFCCellToFaceStencil(mesh)) - {} + { + if (extendedCellToFaceStencil::debug) + { + Info<< "Generated centred stencil " << type() + << nl << endl; + writeStencilStats(Info, stencil(), map()); + } + } //- Destructor diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCPCCellToFaceStencilObject.H b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCPCCellToFaceStencilObject.H index 392b544f409..68608f210bc 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCPCCellToFaceStencilObject.H +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCPCCellToFaceStencilObject.H @@ -67,7 +67,14 @@ public: : MeshObject<fvMesh, centredCPCCellToFaceStencilObject>(mesh), extendedCentredCellToFaceStencil(CPCCellToFaceStencil(mesh)) - {} + { + if (extendedCellToFaceStencil::debug) + { + Info<< "Generated centred stencil " << type() + << nl << endl; + writeStencilStats(Info, stencil(), map()); + } + } // Destructor diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredFECCellToFaceStencilObject.H b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredFECCellToFaceStencilObject.H index d8efe704a7c..3a802374f7a 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredFECCellToFaceStencilObject.H +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredFECCellToFaceStencilObject.H @@ -67,7 +67,14 @@ public: : MeshObject<fvMesh, centredFECCellToFaceStencilObject>(mesh), extendedCentredCellToFaceStencil(FECCellToFaceStencil(mesh)) - {} + { + if (extendedCellToFaceStencil::debug) + { + Info<< "Generated centred stencil " << type() + << nl << endl; + writeStencilStats(Info, stencil(), map()); + } + } // Destructor diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/pureUpwindCFCCellToFaceStencilObject.C b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/pureUpwindCFCCellToFaceStencilObject.C new file mode 100644 index 00000000000..541e7f40d0c --- /dev/null +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/pureUpwindCFCCellToFaceStencilObject.C @@ -0,0 +1,38 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "pureUpwindCFCCellToFaceStencilObject.H" + + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(pureUpwindCFCCellToFaceStencilObject, 0); +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/pureUpwindCFCCellToFaceStencilObject.H b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/pureUpwindCFCCellToFaceStencilObject.H new file mode 100644 index 00000000000..c3be0efe7e0 --- /dev/null +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/pureUpwindCFCCellToFaceStencilObject.H @@ -0,0 +1,95 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::upwindCFCCellToFaceStencilObject + +Description + +SourceFiles + +\*---------------------------------------------------------------------------*/ + +#ifndef pureUpwindCFCCellToFaceStencilObject_H +#define pureUpwindCFCCellToFaceStencilObject_H + +#include "extendedUpwindCellToFaceStencil.H" +#include "CFCCellToFaceStencil.H" +#include "MeshObject.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class pureUpwindCFCCellToFaceStencilObject Declaration +\*---------------------------------------------------------------------------*/ + +class pureUpwindCFCCellToFaceStencilObject +: + public MeshObject<fvMesh, pureUpwindCFCCellToFaceStencilObject>, + public extendedUpwindCellToFaceStencil +{ + +public: + + TypeName("pureUpwindCFCCellToFaceStencil"); + + // Constructors + + //- Construct from uncompacted face stencil + explicit pureUpwindCFCCellToFaceStencilObject + ( + const fvMesh& mesh + ) + : + MeshObject<fvMesh, pureUpwindCFCCellToFaceStencilObject>(mesh), + extendedUpwindCellToFaceStencil(CFCCellToFaceStencil(mesh)) + { + if (extendedCellToFaceStencil::debug) + { + Info<< "Generated pure upwind stencil " << type() + << nl << endl; + writeStencilStats(Info, ownStencil(), ownMap()); + } + } + + + // Destructor + + virtual ~pureUpwindCFCCellToFaceStencilObject() + {} +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindCECCellToFaceStencilObject.H b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindCECCellToFaceStencilObject.H index 05c0c397795..d23cdacfef1 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindCECCellToFaceStencilObject.H +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindCECCellToFaceStencilObject.H @@ -63,6 +63,7 @@ public: explicit upwindCECCellToFaceStencilObject ( const fvMesh& mesh, + const bool pureUpwind, const scalar minOpposedness ) : @@ -70,9 +71,17 @@ public: extendedUpwindCellToFaceStencil ( CECCellToFaceStencil(mesh), + pureUpwind, minOpposedness ) - {} + { + if (extendedCellToFaceStencil::debug) + { + Info<< "Generated off-centred stencil " << type() + << nl << endl; + writeStencilStats(Info, ownStencil(), ownMap()); + } + } // Destructor diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindCFCCellToFaceStencilObject.H b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindCFCCellToFaceStencilObject.H index b335d724d0c..f543f6d8e63 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindCFCCellToFaceStencilObject.H +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindCFCCellToFaceStencilObject.H @@ -35,7 +35,7 @@ SourceFiles #define upwindCFCCellToFaceStencilObject_H #include "extendedUpwindCellToFaceStencil.H" -#include "FECCellToFaceStencil.H" +#include "CFCCellToFaceStencil.H" #include "MeshObject.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -63,16 +63,25 @@ public: explicit upwindCFCCellToFaceStencilObject ( const fvMesh& mesh, + const bool pureUpwind, const scalar minOpposedness ) : MeshObject<fvMesh, upwindCFCCellToFaceStencilObject>(mesh), extendedUpwindCellToFaceStencil ( - FECCellToFaceStencil(mesh), + CFCCellToFaceStencil(mesh), + pureUpwind, minOpposedness ) - {} + { + if (extendedCellToFaceStencil::debug) + { + Info<< "Generated off-centred stencil " << type() + << nl << endl; + writeStencilStats(Info, ownStencil(), ownMap()); + } + } // Destructor diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindCPCCellToFaceStencilObject.H b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindCPCCellToFaceStencilObject.H index a0677af97d7..b86e0b00b5b 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindCPCCellToFaceStencilObject.H +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindCPCCellToFaceStencilObject.H @@ -63,6 +63,7 @@ public: explicit upwindCPCCellToFaceStencilObject ( const fvMesh& mesh, + const bool pureUpwind, const scalar minOpposedness ) : @@ -70,9 +71,17 @@ public: extendedUpwindCellToFaceStencil ( CPCCellToFaceStencil(mesh), + pureUpwind, minOpposedness ) - {} + { + if (extendedCellToFaceStencil::debug) + { + Info<< "Generated off-centred stencil " << type() + << nl << endl; + writeStencilStats(Info, ownStencil(), ownMap()); + } + } // Destructor diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindFECCellToFaceStencilObject.H b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindFECCellToFaceStencilObject.H index 84167c9be84..9bd07a16073 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindFECCellToFaceStencilObject.H +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindFECCellToFaceStencilObject.H @@ -63,6 +63,7 @@ public: explicit upwindFECCellToFaceStencilObject ( const fvMesh& mesh, + const bool pureUpwind, const scalar minOpposedness ) : @@ -70,9 +71,17 @@ public: extendedUpwindCellToFaceStencil ( FECCellToFaceStencil(mesh), + pureUpwind, minOpposedness ) - {} + { + if (extendedCellToFaceStencil::debug) + { + Info<< "Generated off-centred stencil " << type() + << nl << endl; + writeStencilStats(Info, ownStencil(), ownMap()); + } + } // Destructor diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencil.C b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencil.C index 968941177b7..5b1cb12e342 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencil.C +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencil.C @@ -29,8 +29,70 @@ License #include "syncTools.H" #include "SortableList.H" +/* * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * */ + +defineTypeNameAndDebug(Foam::extendedCellToFaceStencil, 0); + + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +void Foam::extendedCellToFaceStencil::writeStencilStats +( + Ostream& os, + const labelListList& stencil, + const mapDistribute& map +) +{ + label sumSize = 0; + label nSum = 0; + label minSize = labelMax; + label maxSize = labelMin; + + forAll(stencil, i) + { + const labelList& sCells = stencil[i]; + + if (sCells.size() > 0) + { + sumSize += sCells.size(); + nSum++; + minSize = min(minSize, sCells.size()); + maxSize = max(maxSize, sCells.size()); + } + } + reduce(sumSize, sumOp<label>()); + reduce(nSum, sumOp<label>()); + + reduce(minSize, minOp<label>()); + reduce(maxSize, maxOp<label>()); + + os << "Stencil size :" << nl + << " average : " << scalar(sumSize)/nSum << nl + << " min : " << minSize << nl + << " max : " << maxSize << nl + << endl; + + // Sum all sent data + label nSent = 0; + label nLocal = 0; + forAll(map.subMap(), procI) + { + if (procI != Pstream::myProcNo()) + { + nSent += map.subMap()[procI].size(); + } + else + { + nLocal += map.subMap()[procI].size(); + } + } + + os << "Local data size : " << returnReduce(nLocal, sumOp<label>()) << nl + << "Sent data size : " << returnReduce(nSent, sumOp<label>()) << nl + << endl; +} + + Foam::autoPtr<Foam::mapDistribute> Foam::extendedCellToFaceStencil::calcDistributeMap ( @@ -211,8 +273,9 @@ Foam::extendedCellToFaceStencil::calcDistributeMap } } + // Constuct map for distribution of compact data. - return autoPtr<mapDistribute> + autoPtr<mapDistribute> mapPtr ( new mapDistribute ( @@ -222,6 +285,8 @@ Foam::extendedCellToFaceStencil::calcDistributeMap true // reuse send/recv maps. ) ); + + return mapPtr; } diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencil.H b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencil.H index 269eb391ab3..bb0e92d11fc 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencil.H +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencil.H @@ -88,8 +88,22 @@ private: void operator=(const extendedCellToFaceStencil&); +protected: + + //- Write some statistics about stencil + static void writeStencilStats + ( + Ostream& os, + const labelListList& stencil, + const mapDistribute& map + ); + public: + // Declare name of the class and its debug switch + ClassName("extendedCellToFaceStencil"); + + // Constructors //- Construct from mesh diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencilTemplates.C b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencilTemplates.C index 8e55a5266b6..9ee125de9b5 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencilTemplates.C +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencilTemplates.C @@ -101,7 +101,10 @@ Foam::extendedCellToFaceStencil::weightedSum ( fld.name(), mesh.time().timeName(), - mesh + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false ), mesh, dimensioned<Type> diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencil.C b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencil.C index 11f094415da..80e47855a8c 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencil.C +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencil.C @@ -375,10 +375,12 @@ void Foam::extendedUpwindCellToFaceStencil::transportStencils Foam::extendedUpwindCellToFaceStencil::extendedUpwindCellToFaceStencil ( const cellToFaceStencil& stencil, + const bool pureUpwind, const scalar minOpposedness ) : - extendedCellToFaceStencil(stencil.mesh()) + extendedCellToFaceStencil(stencil.mesh()), + pureUpwind_(pureUpwind) { //forAll(stencil, faceI) //{ @@ -430,6 +432,134 @@ Foam::extendedUpwindCellToFaceStencil::extendedUpwindCellToFaceStencil stencil.globalNumbering(), neiStencil_ ); + + // stencil now in compact form + if (pureUpwind_) + { + const fvMesh& mesh = dynamic_cast<const fvMesh&>(stencil.mesh()); + + List<List<point> > stencilPoints(ownStencil_.size()); + + // Owner stencil + // ~~~~~~~~~~~~~ + + collectData(ownMapPtr_(), ownStencil_, mesh.C(), stencilPoints); + + // Mask off all stencil points on wrong side of face + forAll(stencilPoints, faceI) + { + const point& fc = mesh.faceCentres()[faceI]; + const vector& fArea = mesh.faceAreas()[faceI]; + + const List<point>& points = stencilPoints[faceI]; + const labelList& stencil = ownStencil_[faceI]; + + DynamicList<label> newStencil(stencil.size()); + forAll(points, i) + { + if (((points[i]-fc) & fArea) < 0) + { + newStencil.append(stencil[i]); + } + } + if (newStencil.size() != stencil.size()) + { + ownStencil_[faceI].transfer(newStencil); + } + } + + + // Neighbour stencil + // ~~~~~~~~~~~~~~~~~ + + collectData(neiMapPtr_(), neiStencil_, mesh.C(), stencilPoints); + + // Mask off all stencil points on wrong side of face + forAll(stencilPoints, faceI) + { + const point& fc = mesh.faceCentres()[faceI]; + const vector& fArea = mesh.faceAreas()[faceI]; + + const List<point>& points = stencilPoints[faceI]; + const labelList& stencil = neiStencil_[faceI]; + + DynamicList<label> newStencil(stencil.size()); + forAll(points, i) + { + if (((points[i]-fc) & fArea) > 0) + { + newStencil.append(stencil[i]); + } + } + if (newStencil.size() != stencil.size()) + { + neiStencil_[faceI].transfer(newStencil); + } + } + + // Note: could compact schedule as well. for if cells are not needed + // across any boundary anymore. However relatively rare. + } +} + + +Foam::extendedUpwindCellToFaceStencil::extendedUpwindCellToFaceStencil +( + const cellToFaceStencil& stencil +) +: + extendedCellToFaceStencil(stencil.mesh()), + pureUpwind_(true) +{ + // Calculate stencil points with full stencil + + ownStencil_ = stencil; + + ownMapPtr_ = calcDistributeMap + ( + stencil.mesh(), + stencil.globalNumbering(), + ownStencil_ + ); + + const fvMesh& mesh = dynamic_cast<const fvMesh&>(stencil.mesh()); + + List<List<point> > stencilPoints(ownStencil_.size()); + collectData(ownMapPtr_(), ownStencil_, mesh.C(), stencilPoints); + + // Split stencil into owner and neighbour + neiStencil_.setSize(ownStencil_.size()); + + forAll(stencilPoints, faceI) + { + const point& fc = mesh.faceCentres()[faceI]; + const vector& fArea = mesh.faceAreas()[faceI]; + + const List<point>& points = stencilPoints[faceI]; + const labelList& stencil = ownStencil_[faceI]; + + DynamicList<label> newOwnStencil(stencil.size()); + DynamicList<label> newNeiStencil(stencil.size()); + forAll(points, i) + { + if (((points[i]-fc) & fArea) > 0) + { + newNeiStencil.append(stencil[i]); + } + else + { + newOwnStencil.append(stencil[i]); + } + } + if (newNeiStencil.size() > 0) + { + ownStencil_[faceI].transfer(newOwnStencil); + neiStencil_[faceI].transfer(newNeiStencil); + } + } + + // Should compact schedule. Or have both return the same schedule. + neiMapPtr_.reset(new mapDistribute(ownMapPtr_())); } diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencil.H b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencil.H index 9b9ef1fd577..408d28ad891 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencil.H +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencil.H @@ -27,7 +27,14 @@ Class Description Creates upwind stencil by shifting a centred stencil to upwind and downwind - faces. + faces and optionally removing all non-(up/down)wind faces ('pureUpwind'). + + Note: the minOpposedness parameter is to decide which upwind and + downwind faces to combine the stencils from. If myArea is the + local area and upwindArea + the area of the possible upwind candidate it will be included if + (upwindArea & myArea)/magSqr(myArea) > minOpposedness + so this includes both cosine and area. WIP. SourceFiles extendedUpwindCellToFaceStencil.C @@ -57,6 +64,9 @@ class extendedUpwindCellToFaceStencil { // Private data + //- Does stencil contain upwind points only + const bool pureUpwind_; + //- Swap map for getting neigbouring data autoPtr<mapDistribute> ownMapPtr_; autoPtr<mapDistribute> neiMapPtr_; @@ -115,16 +125,31 @@ public: // Constructors - //- Construct from mesh and uncompacted face stencil + //- Construct from mesh and uncompacted centred face stencil. + // Transports facestencil to create owner and neighbour versions. + // pureUpwind to remove any remaining downwind cells. extendedUpwindCellToFaceStencil ( const cellToFaceStencil&, + const bool pureUpwind, const scalar minOpposedness ); + //- Construct from mesh and uncompacted centred face stencil. Splits + // stencil into owner and neighbour (so always pure upwind) + extendedUpwindCellToFaceStencil + ( + const cellToFaceStencil& + ); + // Member Functions + bool pureUpwind() const + { + return pureUpwind_; + } + //- Return reference to the parallel distribution map const mapDistribute& ownMap() const { diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencilTemplates.C b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencilTemplates.C index be9409c9f9b..ee545754fa2 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencilTemplates.C +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencilTemplates.C @@ -54,7 +54,10 @@ Foam::extendedUpwindCellToFaceStencil::weightedSum ( fld.name(), mesh.time().timeName(), - mesh + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false ), mesh, dimensioned<Type> diff --git a/src/finiteVolume/fvMesh/extendedStencil/faceToCell/extendedFaceToCellStencil.H b/src/finiteVolume/fvMesh/extendedStencil/faceToCell/extendedFaceToCellStencil.H index 38681d15f56..5c334786c60 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/faceToCell/extendedFaceToCellStencil.H +++ b/src/finiteVolume/fvMesh/extendedStencil/faceToCell/extendedFaceToCellStencil.H @@ -26,6 +26,8 @@ Class Foam::extendedFaceToCellStencil Description + Note: transformations on coupled patches not supported. Problem is the + positions of cells reachable through these patches. SourceFiles extendedFaceToCellStencil.C diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C index e795d965048..fb9611e9c8a 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C @@ -49,7 +49,7 @@ Foam::CentredFitData<Polynomial>::CentredFitData Polynomial > ( - mesh, stencil, linearLimitFactor, centralWeight + mesh, stencil, true, linearLimitFactor, centralWeight ), coeffs_(mesh.nFaces()) { diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C index c3a5bf77757..dcea0902668 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C @@ -36,12 +36,14 @@ Foam::FitData<Form, ExtendedStencil, Polynomial>::FitData ( const fvMesh& mesh, const ExtendedStencil& stencil, + const bool linearCorrection, const scalar linearLimitFactor, const scalar centralWeight ) : MeshObject<fvMesh, Form>(mesh), stencil_(stencil), + linearCorrection_(linearCorrection), linearLimitFactor_(linearLimitFactor), centralWeight_(centralWeight), # ifdef SPHERICAL_GEOMETRY @@ -143,7 +145,10 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit // Setup the point weights scalarList wts(C.size(), scalar(1)); wts[0] = centralWeight_; - wts[1] = centralWeight_; + if (linearCorrection_) + { + wts[1] = centralWeight_; + } // Reference point point p0 = this->mesh().faceCentres()[facei]; @@ -218,10 +223,20 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit } } - goodFit = - (mag(coeffsi[0] - wLin) < linearLimitFactor_*wLin) - && (mag(coeffsi[1] - (1 - wLin)) < linearLimitFactor_*(1 - wLin)) - && maxCoeffi <= 1; + if (linearCorrection_) + { + goodFit = + (mag(coeffsi[0] - wLin) < linearLimitFactor_*wLin) + && (mag(coeffsi[1] - (1 - wLin)) < linearLimitFactor_*(1 - wLin)) + && maxCoeffi <= 1; + } + else + { + // Upwind: weight on face is 1. + goodFit = + (mag(coeffsi[0] - 1.0) < linearLimitFactor_*1.0) + && maxCoeffi <= 1; + } // if (goodFit && iIt > 0) // { @@ -251,7 +266,10 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit // } wts[0] *= 10; - wts[1] *= 10; + if (linearCorrection_) + { + wts[1] *= 10; + } for(label j = 0; j < B.m(); j++) { @@ -269,9 +287,17 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit if (goodFit) { - // Remove the uncorrected linear ocefficients - coeffsi[0] -= wLin; - coeffsi[1] -= 1 - wLin; + if (linearCorrection_) + { + // Remove the uncorrected linear coefficients + coeffsi[0] -= wLin; + coeffsi[1] -= 1 - wLin; + } + else + { + // Remove the uncorrected upwind coefficients + coeffsi[0] -= 1.0; + } } else { diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H index b151a752041..cb7bd9ceae5 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H @@ -26,7 +26,11 @@ Class Foam::FitData Description - Data for the upwinded and centred polynomial fit interpolation schemes + Data for the upwinded and centred polynomial fit interpolation schemes. + The linearCorrection_ determines whether the fit is for a corrected + linear scheme (first two coefficients are corrections for owner and + neighbour) or a pure upwind scheme (first coefficient is correction for + owner ; weight on face taken as 1). SourceFiles FitData.C @@ -58,7 +62,11 @@ class FitData //- The stencil the fit is based on const ExtendedStencil& stencil_; - //- Factor the fit is allowed to deviate from linear. + //- Is scheme correction on linear (true) or on upwind (false) + const bool linearCorrection_; + + //- Factor the fit is allowed to deviate from the base scheme + // (linear or pure upwind) // This limits the amount of high-order correction and increases // stability on bad meshes const scalar linearLimitFactor_; @@ -84,11 +92,6 @@ class FitData const label faci ); - //- Calculate the fit for the all the mesh faces - // and set the coefficients - // virtual void calcFit(); - - public: //TypeName("FitData"); @@ -101,6 +104,7 @@ public: ( const fvMesh& mesh, const ExtendedStencil& stencil, + const bool linearCorrection, const scalar linearLimitFactor, const scalar centralWeight ); @@ -119,12 +123,18 @@ public: return stencil_; } + bool linearCorrection() const + { + return linearCorrection_; + } + //- Calculate the fit for the specified face and set the coefficients void calcFit ( scalarList& coeffsi, // coefficients to be set const List<point>&, // Stencil points - const scalar wLin, // Linear weight + const scalar wLin, // Weight for linear approximation (weights + // nearest neighbours) const label faci // Current face index ); @@ -132,7 +142,7 @@ public: virtual void calcFit() = 0; - //- Delete the data when the mesh moves not implemented + //- Recalculate weights (but not stencil) when the mesh moves bool movePoints(); }; diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/PureUpwindFitScheme/PureUpwindFitScheme.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/PureUpwindFitScheme/PureUpwindFitScheme.H new file mode 100644 index 00000000000..6b0b792f9a2 --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/PureUpwindFitScheme/PureUpwindFitScheme.H @@ -0,0 +1,186 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::PureUpwindFitScheme + +Description + Upwind biased fit surface interpolation scheme that applies an explicit + correction to upwind. + +\*---------------------------------------------------------------------------*/ + +#ifndef PureUpwindFitScheme_H +#define PureUpwindFitScheme_H + +#include "UpwindFitData.H" +#include "upwind.H" +#include "Switch.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class PureUpwindFitScheme Declaration +\*---------------------------------------------------------------------------*/ + +template<class Type, class Polynomial, class Stencil> +class PureUpwindFitScheme +: + public upwind<Type> +{ + // Private Data + + //- Factor the fit is allowed to deviate from linear. + // This limits the amount of high-order correction and increases + // stability on bad meshes + const scalar linearLimitFactor_; + + //- Weights for central stencil + const scalar centralWeight_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + PureUpwindFitScheme(const PureUpwindFitScheme&); + + //- Disallow default bitwise assignment + void operator=(const PureUpwindFitScheme&); + + +public: + + //- Runtime type information + TypeName("PureUpwindFitScheme"); + + + // Constructors + + //- Construct from mesh and Istream + // The name of the flux field is read from the Istream and looked-up + // from the mesh objectRegistry + PureUpwindFitScheme(const fvMesh& mesh, Istream& is) + : + upwind<Type> + ( + mesh, + mesh.lookupObject<surfaceScalarField>(word(is)) + ), + linearLimitFactor_(readScalar(is)), + centralWeight_(1000) + {} + + + //- Construct from mesh, faceFlux and Istream + PureUpwindFitScheme + ( + const fvMesh& mesh, + const surfaceScalarField& faceFlux, + Istream& is + ) + : + upwind<Type>(mesh, faceFlux), + linearLimitFactor_(readScalar(is)), + centralWeight_(1000) + {} + + + // Member Functions + + //- Return true if this scheme uses an explicit correction + virtual bool corrected() const + { + return true; + } + + //- Return the explicit correction to the face-interpolate + virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > + correction + ( + const GeometricField<Type, fvPatchField, volMesh>& vf + ) const + { + const fvMesh& mesh = this->mesh(); + + // Use the owner/neighbour splitting constructor + const extendedUpwindCellToFaceStencil& stencil = Stencil::New(mesh); + + const UpwindFitData<Polynomial>& ufd = + UpwindFitData<Polynomial>::New + ( + mesh, + stencil, + false, //offset to upwind + linearLimitFactor_, + centralWeight_ + ); + + const List<scalarList>& fo = ufd.owncoeffs(); + const List<scalarList>& fn = ufd.neicoeffs(); + + return stencil.weightedSum(this->faceFlux_, vf, fo, fn); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// Add the patch constructor functions to the hash tables + +#define makePureUpwindFitSurfaceInterpolationTypeScheme(SS, POLYNOMIAL, STENCIL, TYPE) \ + \ +typedef PureUpwindFitScheme<TYPE, POLYNOMIAL, STENCIL> \ + PureUpwindFitScheme##TYPE##POLYNOMIAL##STENCIL##_; \ +defineTemplateTypeNameAndDebugWithName \ + (PureUpwindFitScheme##TYPE##POLYNOMIAL##STENCIL##_, #SS, 0); \ + \ +surfaceInterpolationScheme<TYPE>::addMeshConstructorToTable \ +<PureUpwindFitScheme<TYPE, POLYNOMIAL, STENCIL> > \ + add##SS##STENCIL##TYPE##MeshConstructorToTable_; \ + \ +surfaceInterpolationScheme<TYPE>::addMeshFluxConstructorToTable \ +<PureUpwindFitScheme<TYPE, POLYNOMIAL, STENCIL> > \ + add##SS##STENCIL##TYPE##MeshFluxConstructorToTable_; + +#define makePureUpwindFitSurfaceInterpolationScheme(SS, POLYNOMIAL, STENCIL) \ + \ +makePureUpwindFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,scalar) \ +makePureUpwindFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,vector) \ +makePureUpwindFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,sphericalTensor) \ +makePureUpwindFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,symmTensor) \ +makePureUpwindFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,tensor) + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitData.C index 9bd5451cf06..00c490560be 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitData.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitData.C @@ -38,6 +38,7 @@ Foam::UpwindFitData<Polynomial>::UpwindFitData ( const fvMesh& mesh, const extendedUpwindCellToFaceStencil& stencil, + const bool linearCorrection, const scalar linearLimitFactor, const scalar centralWeight ) @@ -49,7 +50,7 @@ Foam::UpwindFitData<Polynomial>::UpwindFitData Polynomial > ( - mesh, stencil, linearLimitFactor, centralWeight + mesh, stencil, linearCorrection, linearLimitFactor, centralWeight ), owncoeffs_(mesh.nFaces()), neicoeffs_(mesh.nFaces()) @@ -77,30 +78,25 @@ void Foam::UpwindFitData<Polynomial>::calcFit() { const fvMesh& mesh = this->mesh(); + const surfaceScalarField& w = mesh.surfaceInterpolation::weights(); + const surfaceScalarField::GeometricBoundaryField& bw = w.boundaryField(); + + // Owner stencil weights + // ~~~~~~~~~~~~~~~~~~~~~ + // Get the cell/face centres in stencil order. - // Upwind face stencils no good for triangles or tets. - // Need bigger stencils - List<List<point> > ownStencilPoints(mesh.nFaces()); + List<List<point> > stencilPoints(mesh.nFaces()); this->stencil().collectData ( this->stencil().ownMap(), this->stencil().ownStencil(), mesh.C(), - ownStencilPoints - ); - List<List<point> > neiStencilPoints(mesh.nFaces()); - this->stencil().collectData - ( - this->stencil().neiMap(), - this->stencil().neiStencil(), - mesh.C(), - neiStencilPoints + stencilPoints ); - // find the fit coefficients for every owner and neighbour of ever face - - const surfaceScalarField& w = mesh.surfaceInterpolation::weights(); + // find the fit coefficients for every owner + Pout<< "-- Owner --" << endl; for(label facei = 0; facei < mesh.nInternalFaces(); facei++) { FitData @@ -108,16 +104,17 @@ void Foam::UpwindFitData<Polynomial>::calcFit() UpwindFitData<Polynomial>, extendedUpwindCellToFaceStencil, Polynomial - >::calcFit(owncoeffs_[facei], ownStencilPoints[facei], w[facei], facei); - FitData - < - UpwindFitData<Polynomial>, - extendedUpwindCellToFaceStencil, - Polynomial - >::calcFit(neicoeffs_[facei], neiStencilPoints[facei], w[facei], facei); - } + >::calcFit(owncoeffs_[facei], stencilPoints[facei], w[facei], facei); - const surfaceScalarField::GeometricBoundaryField& bw = w.boundaryField(); + Pout<< " facei:" << facei + << " at:" << mesh.faceCentres()[facei] << endl; + forAll(owncoeffs_[facei], i) + { + Pout<< " point:" << stencilPoints[facei][i] + << "\tweight:" << owncoeffs_[facei][i] + << endl; + } + } forAll(bw, patchi) { @@ -136,8 +133,58 @@ void Foam::UpwindFitData<Polynomial>::calcFit() Polynomial >::calcFit ( - owncoeffs_[facei], ownStencilPoints[facei], pw[i], facei + owncoeffs_[facei], stencilPoints[facei], pw[i], facei ); + facei++; + } + } + } + + + // Neighbour stencil weights + // ~~~~~~~~~~~~~~~~~~~~~~~~~ + + // Note:reuse stencilPoints since is major storage + this->stencil().collectData + ( + this->stencil().neiMap(), + this->stencil().neiStencil(), + mesh.C(), + stencilPoints + ); + + // find the fit coefficients for every neighbour + + Pout<< "-- Neighbour --" << endl; + for(label facei = 0; facei < mesh.nInternalFaces(); facei++) + { + FitData + < + UpwindFitData<Polynomial>, + extendedUpwindCellToFaceStencil, + Polynomial + >::calcFit(neicoeffs_[facei], stencilPoints[facei], w[facei], facei); + + Pout<< " facei:" << facei + << " at:" << mesh.faceCentres()[facei] << endl; + forAll(neicoeffs_[facei], i) + { + Pout<< " point:" << stencilPoints[facei][i] + << "\tweight:" << neicoeffs_[facei][i] + << endl; + } + } + + forAll(bw, patchi) + { + const fvsPatchScalarField& pw = bw[patchi]; + + if (pw.coupled()) + { + label facei = pw.patch().patch().start(); + + forAll(pw, i) + { FitData < UpwindFitData<Polynomial>, @@ -145,7 +192,7 @@ void Foam::UpwindFitData<Polynomial>::calcFit() Polynomial >::calcFit ( - neicoeffs_[facei], neiStencilPoints[facei], pw[i], facei + neicoeffs_[facei], stencilPoints[facei], pw[i], facei ); facei++; } diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitData.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitData.H index f5367aed2d4..a62a3002c32 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitData.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitData.H @@ -27,7 +27,9 @@ Class Description Data for the quadratic fit correction interpolation scheme to be used with - upwind biased stencil + upwind biased stencil. + - linearCorrection = true : fit calculated for corrected linear scheme + - linearCorrection = false : fit calculated for corrected upwind scheme SourceFiles UpwindFitData.C @@ -90,6 +92,7 @@ public: ( const fvMesh& mesh, const extendedUpwindCellToFaceStencil& stencil, + const bool linearCorrection, const scalar linearLimitFactor, const scalar centralWeight ); diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitScheme.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitScheme.H index 465ce019a46..46fe3fe1ba0 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitScheme.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitScheme.H @@ -129,6 +129,7 @@ public: const extendedUpwindCellToFaceStencil& stencil = Stencil::New ( mesh, + false, //pureUpwind scalar(0.5) ); @@ -137,6 +138,7 @@ public: ( mesh, stencil, + true, //calculate as offset to linear linearLimitFactor_, centralWeight_ ); diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearPureUpwindFit/linearPureUpwindFit.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearPureUpwindFit/linearPureUpwindFit.C new file mode 100644 index 00000000000..28aa8e7434b --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearPureUpwindFit/linearPureUpwindFit.C @@ -0,0 +1,49 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "PureUpwindFitScheme.H" +#include "linearFitPolynomial.H" +#include "pureUpwindCFCCellToFaceStencilObject.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + defineTemplateTypeNameAndDebug + ( + UpwindFitData<linearFitPolynomial>, + 0 + ); + + makePureUpwindFitSurfaceInterpolationScheme + ( + linearPureUpwindFit, + linearFitPolynomial, + pureUpwindCFCCellToFaceStencilObject + ); +} + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearPureUpwindFit/quadraticLinearPureUpwindFit.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearPureUpwindFit/quadraticLinearPureUpwindFit.C new file mode 100644 index 00000000000..513de347d75 --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearPureUpwindFit/quadraticLinearPureUpwindFit.C @@ -0,0 +1,45 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "PureUpwindFitScheme.H" +#include "quadraticLinearUpwindFitPolynomial.H" +#include "upwindCFCCellToFaceStencilObject.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + // Use stencil with three upwind cells: + // upwindCFCCellToFaceStencilObject + pureUpwind + makePureUpwindFitSurfaceInterpolationScheme + ( + quadraticLinearPureUpwindFit, + quadraticLinearUpwindFitPolynomial, + upwindCFCCellToFaceStencilObject + ); +} + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFit.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFit.C index 75be9da01e8..d89107ab71e 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFit.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFit.C @@ -26,7 +26,7 @@ License #include "UpwindFitScheme.H" #include "quadraticLinearUpwindFitPolynomial.H" -#include "upwindCFCCellToFaceStencilObject.H" +#include "upwindFECCellToFaceStencilObject.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -42,7 +42,7 @@ namespace Foam ( quadraticLinearUpwindFit, quadraticLinearUpwindFitPolynomial, - upwindCFCCellToFaceStencilObject + upwindFECCellToFaceStencilObject ); } -- GitLab