From 4e4f9e0a0a1518fa4124ea8955fe0201f6b8e7ab Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Sun, 12 Oct 2008 11:42:20 +0100
Subject: [PATCH] Added quadraticFit snGrad.
---
.../quadraticFitSnGrad/quadraticFitSnGrad.H | 161 +++++++++
.../quadraticFitSnGradData.C | 325 ++++++++++++++++++
.../quadraticFitSnGradData.H | 131 +++++++
.../quadraticFitSnGrad/quadraticFitSnGrads.C | 44 +++
4 files changed, 661 insertions(+)
create mode 100644 src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrad.H
create mode 100644 src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C
create mode 100644 src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.H
create mode 100644 src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrads.C
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrad.H b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrad.H
new file mode 100644
index 00000000000..b74c3d2f205
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrad.H
@@ -0,0 +1,161 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 1991-2007 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
+ quadraticFitSnGrad
+
+Description
+ Simple central-difference snGrad scheme with quadratic fit correction from
+ a larger stencil.
+
+SourceFiles
+ quadraticFitSnGrad.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef quadraticFitSnGrad_H
+#define quadraticFitSnGrad_H
+
+#include "snGradScheme.H"
+#include "quadraticFitSnGradData.H"
+#include "extendedStencil.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+ Class quadraticFitSnGrad Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class quadraticFitSnGrad
+:
+ public snGradScheme<Type>
+{
+ // Private Data
+ //- weights for central stencil
+ const scalar centralWeight_;
+
+ // Private Member Functions
+
+ //- Disallow default bitwise assignment
+ void operator=(const quadraticFitSnGrad&);
+
+
+public:
+
+ //- Runtime type information
+ TypeName("quadraticFit");
+
+
+ // Constructors
+
+ //- Construct from mesh and scheme data
+ quadraticFitSnGrad
+ (
+ const fvMesh& mesh,
+ const scalar centralWeight
+ )
+ :
+ snGradScheme<Type>(mesh),
+ centralWeight_(centralWeight)
+ {}
+
+
+ //- Construct from mesh and data stream
+ quadraticFitSnGrad(const fvMesh& mesh, Istream& is)
+ :
+ snGradScheme<Type>(mesh),
+ centralWeight_(readScalar(is))
+ {}
+
+
+ // Destructor
+
+ virtual ~quadraticFitSnGrad() {}
+
+
+ // Member Functions
+
+ //- Return the interpolation weighting factors for the given field
+ virtual tmp<surfaceScalarField> deltaCoeffs
+ (
+ const GeometricField<Type, fvPatchField, volMesh>&
+ ) const
+ {
+ return this->mesh().deltaCoeffs();
+ }
+
+ //- Return true if this scheme uses an explicit correction
+ virtual bool corrected() const
+ {
+ return true;
+ }
+
+ //- Return the explicit correction to the quadraticFitSnGrad
+ // for the given field
+ virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+ correction(const GeometricField<Type, fvPatchField, volMesh>& vf) const
+ {
+ const fvMesh& mesh = this->mesh();
+
+ const quadraticFitSnGradData& qfd = quadraticFitSnGradData::New
+ (
+ mesh,
+ centralWeight_
+ );
+
+ const extendedStencil& stencil = qfd.stencil();
+ const List<scalarList>& f = qfd.fit();
+
+ tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > sft
+ = stencil.weightedSum(vf, f);
+
+ sft().dimensions() /= dimLength;
+
+ return sft;
+ }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C
new file mode 100644
index 00000000000..133f1e79332
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C
@@ -0,0 +1,325 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 1991-2007 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 "quadraticFitSnGradData.H"
+#include "surfaceFields.H"
+#include "volFields.H"
+#include "SVD.H"
+#include "syncTools.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(quadraticFitSnGradData, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::quadraticFitSnGradData::quadraticFitSnGradData
+(
+ const fvMesh& mesh,
+ const scalar cWeight
+)
+:
+ MeshObject<fvMesh, quadraticFitSnGradData>(mesh),
+ centralWeight_(cWeight),
+ #ifdef SPHERICAL_GEOMETRY
+ dim_(2),
+ #else
+ dim_(mesh.nGeometricD()),
+ #endif
+ minSize_
+ (
+ dim_ == 1 ? 3 :
+ dim_ == 2 ? 6 :
+ dim_ == 3 ? 9 : 0
+ ),
+ stencil_(mesh),
+ fit_(mesh.nInternalFaces())
+{
+ if (debug)
+ {
+ Info << "Contructing quadraticFitSnGradData" << endl;
+ }
+
+ // check input
+ if (centralWeight_ < 1 - SMALL)
+ {
+ FatalErrorIn("quadraticFitSnGradData::quadraticFitSnGradData")
+ << "centralWeight requested = " << centralWeight_
+ << " should not be less than one"
+ << exit(FatalError);
+ }
+
+ if (minSize_ == 0)
+ {
+ FatalErrorIn("quadraticFitSnGradData")
+ << " dimension must be 1,2 or 3, not" << dim_ << exit(FatalError);
+ }
+
+ // store the polynomial size for each face to write out
+ surfaceScalarField snGradPolySize
+ (
+ IOobject
+ (
+ "quadraticFitSnGradPolySize",
+ "constant",
+ mesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("quadraticFitSnGradPolySize", dimless, scalar(0))
+ );
+
+ // Get the cell/face centres in stencil order.
+ // Centred face stencils no good for triangles of tets. Need bigger stencils
+ List<List<point> > stencilPoints(stencil_.stencil().size());
+ stencil_.collectData
+ (
+ mesh.C(),
+ stencilPoints
+ );
+
+ // find the fit coefficients for every face in the mesh
+
+ for(label faci = 0; faci < mesh.nInternalFaces(); faci++)
+ {
+ snGradPolySize[faci] = calcFit(stencilPoints[faci], faci);
+ }
+
+ if (debug)
+ {
+ snGradPolySize.write();
+ Info<< "quadraticFitSnGradData::quadraticFitSnGradData() :"
+ << "Finished constructing polynomialFit data"
+ << endl;
+ }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::quadraticFitSnGradData::findFaceDirs
+(
+ vector& idir, // value changed in return
+ vector& jdir, // value changed in return
+ vector& kdir, // value changed in return
+ const fvMesh& mesh,
+ const label faci
+)
+{
+ idir = mesh.Sf()[faci];
+ idir /= mag(idir);
+
+ #ifndef SPHERICAL_GEOMETRY
+ if (mesh.nGeometricD() <= 2) // find the normal direcion
+ {
+ if (mesh.directions()[0] == -1)
+ {
+ kdir = vector(1, 0, 0);
+ }
+ else if (mesh.directions()[1] == -1)
+ {
+ kdir = vector(0, 1, 0);
+ }
+ else
+ {
+ kdir = vector(0, 0, 1);
+ }
+ }
+ else // 3D so find a direction in the place of the face
+ {
+ const face& f = mesh.faces()[faci];
+ kdir = mesh.points()[f[0]] - mesh.points()[f[1]];
+ }
+ #else
+ // Spherical geometry so kdir is the radial direction
+ kdir = mesh.Cf()[faci];
+ #endif
+
+ if (mesh.nGeometricD() == 3)
+ {
+ // Remove the idir component from kdir and normalise
+ kdir -= (idir & kdir)*idir;
+
+ scalar magk = mag(kdir);
+
+ if (magk < SMALL)
+ {
+ FatalErrorIn("findFaceDirs") << " calculated kdir = zero"
+ << exit(FatalError);
+ }
+ else
+ {
+ kdir /= magk;
+ }
+ }
+
+ jdir = kdir ^ idir;
+}
+
+
+Foam::label Foam::quadraticFitSnGradData::calcFit
+(
+ const List<point>& C,
+ const label faci
+)
+{
+ vector idir(1,0,0);
+ vector jdir(0,1,0);
+ vector kdir(0,0,1);
+ findFaceDirs(idir, jdir, kdir, mesh(), faci);
+
+ scalarList wts(C.size(), scalar(1));
+ wts[0] = centralWeight_;
+ wts[1] = centralWeight_;
+
+ point p0 = mesh().faceCentres()[faci];
+ scalar scale = 0;
+
+ // calculate the matrix of the polynomial components
+ scalarRectangularMatrix B(C.size(), minSize_, scalar(0));
+
+ for(label ip = 0; ip < C.size(); ip++)
+ {
+ const point& p = C[ip];
+
+ scalar px = (p - p0)&idir;
+ scalar py = (p - p0)&jdir;
+ #ifdef SPHERICAL_GEOMETRY
+ scalar pz = mag(p) - mag(p0);
+ #else
+ scalar pz = (p - p0)&kdir;
+ #endif
+
+ if (ip == 0) scale = max(max(mag(px), mag(py)), mag(pz));
+
+ px /= scale;
+ py /= scale;
+ pz /= scale;
+
+ label is = 0;
+
+ B[ip][is++] = wts[0]*wts[ip];
+ B[ip][is++] = wts[0]*wts[ip]*px;
+ B[ip][is++] = wts[ip]*sqr(px);
+
+ if (dim_ >= 2)
+ {
+ B[ip][is++] = wts[ip]*py;
+ B[ip][is++] = wts[ip]*px*py;
+ B[ip][is++] = wts[ip]*sqr(py);
+ }
+ if (dim_ == 3)
+ {
+ B[ip][is++] = wts[ip]*pz;
+ B[ip][is++] = wts[ip]*px*pz;
+ //B[ip][is++] = wts[ip]*py*pz;
+ B[ip][is++] = wts[ip]*sqr(pz);
+ }
+ }
+
+ // Set the fit
+ label stencilSize = C.size();
+ fit_[faci].setSize(stencilSize);
+ scalarList singVals(minSize_);
+ label nSVDzeros = 0;
+
+ const scalar& deltaCoeff = mesh().deltaCoeffs()[faci];
+
+ bool goodFit = false;
+ for(int iIt = 0; iIt < 10 && !goodFit; iIt++)
+ {
+ SVD svd(B, SMALL);
+
+ scalar fit0 = wts[0]*wts[0]*svd.VSinvUt()[1][0]/scale;
+ scalar fit1 = wts[0]*wts[1]*svd.VSinvUt()[1][1]/scale;
+
+ goodFit =
+ fit0 < 0 && fit1 > 0
+ && mag(fit0 + deltaCoeff) < 0.5*deltaCoeff
+ && mag(fit1 - deltaCoeff) < 0.5*deltaCoeff;
+
+ if (goodFit)
+ {
+ fit_[faci][0] = fit0;
+ fit_[faci][1] = fit1;
+ for(label i = 2; i < stencilSize; i++)
+ {
+ fit_[faci][i] = wts[0]*wts[i]*svd.VSinvUt()[1][i]/scale;
+ }
+ singVals = svd.S();
+ nSVDzeros = svd.nZeros();
+ }
+ else // (not good fit so increase weight in the centre and for linear)
+ {
+ wts[0] *= 10;
+ wts[1] *= 10;
+
+ for(label i = 0; i < B.n(); i++)
+ {
+ B[i][0] *= 10;
+ B[i][1] *= 10;
+ }
+
+ for(label j = 0; j < B.m(); j++)
+ {
+ B[0][j] *= 10;
+ B[1][j] *= 10;
+ }
+ }
+ }
+
+ if (goodFit)
+ {
+ // remove the uncorrected snGradScheme coefficients
+ fit_[faci][0] += deltaCoeff;
+ fit_[faci][1] -= deltaCoeff;
+ }
+ else
+ {
+ Pout<< "quadratifFitSnGradData could not fit face " << faci
+ << " fit_[faci][0] = " << fit_[faci][0]
+ << " fit_[faci][1] = " << fit_[faci][1]
+ << " deltaCoeff = " << deltaCoeff << endl;
+ fit_[faci] = 0;
+ }
+
+ return minSize_ - nSVDzeros;
+}
+
+bool Foam::quadraticFitSnGradData::movePoints()
+{
+ notImplemented("quadraticFitSnGradData::movePoints()");
+
+ return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.H b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.H
new file mode 100644
index 00000000000..0f5fda5ce6e
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.H
@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 1991-2007 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
+ quadraticFitSnGradData
+
+Description
+ Data for the quadratic fit correction snGrad scheme
+
+SourceFiles
+ quadraticFitSnGradData.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef quadraticFitSnGradData_H
+#define quadraticFitSnGradData_H
+
+#include "MeshObject.H"
+#include "fvMesh.H"
+#include "extendedStencil.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class quadraticFitSnGradData Declaration
+\*---------------------------------------------------------------------------*/
+
+class quadraticFitSnGradData
+:
+ public MeshObject<fvMesh, quadraticFitSnGradData>
+{
+ // Private data
+
+ const scalar centralWeight_;
+ const label dim_;
+
+ //- minimum stencil size
+ const label minSize_;
+
+ //- Extended stencil addressing
+ extendedStencil stencil_;
+
+ //- For each cell in the mesh store the values which multiply the
+ // values of the stencil to obtain the gradient for each direction
+ List<scalarList> fit_;
+
+
+ // Private member functions
+
+ //- Find the normal direction and i, j and k directions for face faci
+ static void findFaceDirs
+ (
+ vector& idir, // value changed in return
+ vector& jdir, // value changed in return
+ vector& kdir, // value changed in return
+ const fvMesh& mesh,
+ const label faci
+ );
+
+ label calcFit(const List<point>&, const label faci);
+
+
+public:
+
+ TypeName("quadraticFitSnGradData");
+
+
+ // Constructors
+
+ explicit quadraticFitSnGradData
+ (
+ const fvMesh& mesh,
+ const scalar cWeight
+ );
+
+
+ // Destructor
+
+ virtual ~quadraticFitSnGradData()
+ {};
+
+
+ // Member functions
+
+ //- Return reference to the stencil
+ const extendedStencil& stencil() const
+ {
+ return stencil_;
+ }
+
+ //- Return reference to fit coefficients
+ const List<scalarList>& fit() const { return fit_; }
+
+ //- Delete the data when the mesh moves not implemented
+ virtual bool movePoints();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrads.C b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrads.C
new file mode 100644
index 00000000000..aee4af29e3b
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrads.C
@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 1991-2007 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
+
+Description
+ Simple central-difference snGrad scheme with quadratic fit correction from
+ a larger stencil.
+
+\*---------------------------------------------------------------------------*/
+
+#include "quadraticFitSnGrad.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+ makeSnGradScheme(quadraticFitSnGrad)
+}
+}
+
+// ************************************************************************* //
--
GitLab