diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 7bf27115a5..2f301c0406 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -186,10 +186,8 @@ $(schemes)/harmonic/harmonic.C $(schemes)/localBlended/localBlended.C $(schemes)/localMax/localMax.C $(schemes)/localMin/localMin.C -//$(schemes)/linearFit/linearFit.C -//$(schemes)/linearFit/linearFitData.C -$(schemes)/quadraticFit/quadraticFit.C -$(schemes)/quadraticFit/quadraticFitData.C +$(schemes)/linearFit/linearFit.C +$(schemes)/quadraticLinearFit/quadraticLinearFit.C limitedSchemes = $(surfaceInterpolation)/limitedSchemes $(limitedSchemes)/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationSchemes.C diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C b/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C index 480cf15a2a..4aac55df75 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C +++ b/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C @@ -128,6 +128,7 @@ Foam::extendedStencil::weightedSum // Boundaries. Either constrained or calculated so assign value // directly (instead of nicely using operator==) + /* typename GeometricField:: GeometricBoundaryField& bSfCorr = sf.boundaryField(); @@ -150,6 +151,7 @@ Foam::extendedStencil::weightedSum faceI++; } } + */ return tsfCorr; } diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C similarity index 67% rename from src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C rename to src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C index 868d33b6f3..dc8720ac89 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C @@ -24,56 +24,44 @@ License \*---------------------------------------------------------------------------*/ -#include "quadraticFitData.H" +#include "CentredFitData.H" #include "surfaceFields.H" #include "volFields.H" #include "SVD.H" #include "syncTools.H" #include "extendedCentredStencil.H" -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -namespace Foam -{ - defineTypeNameAndDebug(quadraticFitData, 0); -} - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::quadraticFitData::quadraticFitData +template +Foam::CentredFitData::CentredFitData ( const fvMesh& mesh, const extendedCentredStencil& stencil, const scalar linearLimitFactor, - const scalar cWeight + const scalar centralWeight ) : - MeshObject(mesh), + MeshObject >(mesh), linearLimitFactor_(linearLimitFactor), - centralWeight_(cWeight), + centralWeight_(centralWeight), # ifdef SPHERICAL_GEOMETRY dim_(2), # else dim_(mesh.nGeometricD()), # endif - minSize_ - ( - dim_ == 1 ? 3 : - dim_ == 2 ? 6 : - dim_ == 3 ? 7 : 0 - ), - fit_(mesh.nInternalFaces()) + minSize_(Polynomial::nTerms(dim_)), + coeffs_(mesh.nInternalFaces()) { if (debug) { - Info << "Contructing quadraticFitData" << endl; + Info<< "Contructing CentredFitData" << endl; } // check input if (centralWeight_ < 1 - SMALL) { - FatalErrorIn("quadraticFitData::quadraticFitData") + FatalErrorIn("CentredFitData::CentredFitData") << "centralWeight requested = " << centralWeight_ << " should not be less than one" << exit(FatalError); @@ -81,7 +69,7 @@ Foam::quadraticFitData::quadraticFitData if (minSize_ == 0) { - FatalErrorIn("quadraticFitSnGradData") + FatalErrorIn("CentredFitSnGradData") << " dimension must be 1,2 or 3, not" << dim_ << exit(FatalError); } @@ -90,24 +78,20 @@ Foam::quadraticFitData::quadraticFitData ( IOobject ( - "quadraticFitInterpPolySize", + "CentredFitInterpPolySize", "constant", mesh, IOobject::NO_READ, IOobject::NO_WRITE ), mesh, - dimensionedScalar("quadraticFitInterpPolySize", dimless, scalar(0)) + dimensionedScalar("CentredFitInterpPolySize", dimless, scalar(0)) ); // Get the cell/face centres in stencil order. // Centred face stencils no good for triangles of tets. Need bigger stencils List > stencilPoints(mesh.nFaces()); - stencil.collectData - ( - mesh.C(), - stencilPoints - ); + stencil.collectData(mesh.C(), stencilPoints); // find the fit coefficients for every face in the mesh @@ -118,7 +102,7 @@ Foam::quadraticFitData::quadraticFitData if (debug) { - Info<< "quadraticFitData::quadraticFitData() :" + Info<< "CentredFitData::CentredFitData() :" << "Finished constructing polynomialFit data" << endl; @@ -129,7 +113,8 @@ Foam::quadraticFitData::quadraticFitData // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void Foam::quadraticFitData::findFaceDirs +template +void Foam::CentredFitData::findFaceDirs ( vector& idir, // value changed in return vector& jdir, // value changed in return @@ -139,6 +124,7 @@ void Foam::quadraticFitData::findFaceDirs ) { idir = mesh.Sf()[faci]; + //idir = mesh.C()[mesh.neighbour()[faci]] - mesh.C()[mesh.owner()[faci]]; idir /= mag(idir); # ifndef SPHERICAL_GEOMETRY @@ -189,7 +175,8 @@ void Foam::quadraticFitData::findFaceDirs } -Foam::label Foam::quadraticFitData::calcFit +template +Foam::label Foam::CentredFitData::calcFit ( const List& C, const label faci @@ -198,82 +185,74 @@ Foam::label Foam::quadraticFitData::calcFit vector idir(1,0,0); vector jdir(0,1,0); vector kdir(0,0,1); - findFaceDirs(idir, jdir, kdir, mesh(), faci); + findFaceDirs(idir, jdir, kdir, this->mesh(), faci); + // Setup the point weights scalarList wts(C.size(), scalar(1)); wts[0] = centralWeight_; wts[1] = centralWeight_; - point p0 = mesh().faceCentres()[faci]; - scalar scale = 0; + // Reference point + point p0 = this->mesh().faceCentres()[faci]; - // calculate the matrix of the polynomial components + // p0 -> p vector in the face-local coordinate system + vector d; + + // Local coordinate scaling + scalar scale = 1; + + // 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; + d.x() = (p - p0)&idir; + d.y() = (p - p0)&jdir; # ifndef SPHERICAL_GEOMETRY - scalar pz = (p - p0)&kdir; + d.z() = (p - p0)&kdir; # else - scalar pz = mag(p) - mag(p0); + d.z() = mag(p) - mag(p0); # endif if (ip == 0) { - scale = max(max(mag(px), mag(py)), mag(pz)); + scale = cmptMax(cmptMag((d))); } - px /= scale; - py /= scale; - pz /= scale; + // Scale the radius vector + d /= 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); - } + Polynomial::addCoeffs + ( + B[ip], + d, + wts[ip], + dim_ + ); } // Set the fit label stencilSize = C.size(); - fit_[faci].setSize(stencilSize); + coeffs_[faci].setSize(stencilSize); scalarList singVals(minSize_); label nSVDzeros = 0; const GeometricField& w = - mesh().surfaceInterpolation::weights(); + this->mesh().surfaceInterpolation::weights(); bool goodFit = false; for(int iIt = 0; iIt < 10 && !goodFit; iIt++) { SVD svd(B, SMALL); - scalar fit0 = wts[0]*wts[0]*svd.VSinvUt()[0][0]; - scalar fit1 = wts[0]*wts[1]*svd.VSinvUt()[0][1]; - - //goodFit = (fit0 > 0 && fit1 > 0); + scalar fit0 = wts[0]*svd.VSinvUt()[0][0]; + scalar fit1 = wts[1]*svd.VSinvUt()[0][1]; goodFit = - (mag(fit0 - w[faci])/w[faci] < linearLimitFactor_) - && (mag(fit1 - (1 - w[faci]))/(1 - w[faci]) < linearLimitFactor_); + (mag(fit0 - w[faci]) < linearLimitFactor_*w[faci]) + && (mag(fit1 - (1 - w[faci])) < linearLimitFactor_*(1 - w[faci])); //scalar w0Err = fit0/w[faci]; //scalar w1Err = fit1/(1 - w[faci]); @@ -284,12 +263,12 @@ Foam::label Foam::quadraticFitData::calcFit if (goodFit) { - fit_[faci][0] = fit0; - fit_[faci][1] = fit1; + coeffs_[faci][0] = fit0; + coeffs_[faci][1] = fit1; for(label i=2; i::calcFit" + "(const List& C, const label faci" + ) << "Could not fit face " << faci + << ", reverting to linear." << nl + << " Weights " + << coeffs_[faci][0] << " " << w[faci] << nl + << " Linear weights " + << coeffs_[faci][1] << " " << 1 - w[faci] << endl; + } - fit_[faci] = 0; + coeffs_[faci] = 0; } return minSize_ - nSVDzeros; } -bool Foam::quadraticFitData::movePoints() +template +bool Foam::CentredFitData::movePoints() { - notImplemented("quadraticFitData::movePoints()"); + notImplemented("CentredFitData::movePoints()"); return true; } diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.H similarity index 84% rename from src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H rename to src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.H index f0da60901d..c225642964 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.H @@ -23,18 +23,18 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Class - Foam::quadraticFitData + Foam::CentredFitData Description Data for the quadratic fit correction interpolation scheme SourceFiles - quadraticFitData.C + CentredFitData.C \*---------------------------------------------------------------------------*/ -#ifndef quadraticFitData_H -#define quadraticFitData_H +#ifndef CentredFitData_H +#define CentredFitData_H #include "MeshObject.H" #include "fvMesh.H" @@ -47,12 +47,13 @@ namespace Foam class extendedCentredStencil; /*---------------------------------------------------------------------------*\ - Class quadraticFitData Declaration + Class CentredFitData Declaration \*---------------------------------------------------------------------------*/ -class quadraticFitData +template +class CentredFitData : - public MeshObject + public MeshObject > { // Private data @@ -72,7 +73,7 @@ class quadraticFitData //- For each cell in the mesh store the values which multiply the // values of the stencil to obtain the gradient for each direction - List fit_; + List coeffs_; // Private member functions @@ -92,12 +93,13 @@ class quadraticFitData public: - TypeName("quadraticFitData"); + TypeName("CentredFitData"); // Constructors - explicit quadraticFitData + //- Construct from components + CentredFitData ( const fvMesh& mesh, const extendedCentredStencil& stencil, @@ -107,16 +109,16 @@ public: //- Destructor - virtual ~quadraticFitData() + virtual ~CentredFitData() {} // Member functions //- Return reference to fit coefficients - const List& fit() const + const List& coeffs() const { - return fit_; + return coeffs_; } //- Delete the data when the mesh moves not implemented @@ -130,6 +132,12 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +#ifdef NoRepository +# include "CentredFitData.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + #endif // ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitScheme.H similarity index 56% rename from src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H rename to src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitScheme.H index c1308c6ba3..6c5ad7aa02 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitScheme.H @@ -23,23 +23,19 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Class - Foam::quadraticFit + Foam::CentredFitScheme Description - Quadratic fit interpolation scheme which applies an explicit correction to - linear. - -SourceFiles - quadraticFit.C + Centred fit surface interpolation scheme which applies an explicit + correction to linear. \*---------------------------------------------------------------------------*/ -#ifndef quadraticFit_H -#define quadraticFit_H +#ifndef CentredFitScheme_H +#define CentredFitScheme_H -#include "quadraticFitData.H" +#include "CentredFitData.H" #include "linear.H" -#include "centredCFCStencilObject.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -47,11 +43,11 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class quadraticFit Declaration + Class CentredFitScheme Declaration \*---------------------------------------------------------------------------*/ -template -class quadraticFit +template +class CentredFitScheme : public linear { @@ -69,31 +65,31 @@ class quadraticFit // Private Member Functions //- Disallow default bitwise copy construct - quadraticFit(const quadraticFit&); + CentredFitScheme(const CentredFitScheme&); //- Disallow default bitwise assignment - void operator=(const quadraticFit&); + void operator=(const CentredFitScheme&); public: //- Runtime type information - TypeName("quadraticFit"); + TypeName("CentredFitScheme"); // Constructors //- Construct from mesh and Istream - quadraticFit(const fvMesh& mesh, Istream& is) + CentredFitScheme(const fvMesh& mesh, Istream& is) : linear(mesh), linearLimitFactor_(readScalar(is)), - centralWeight_(readScalar(is)) + centralWeight_(1000) {} //- Construct from mesh, faceFlux and Istream - quadraticFit + CentredFitScheme ( const fvMesh& mesh, const surfaceScalarField& faceFlux, @@ -102,7 +98,7 @@ public: : linear(mesh), linearLimitFactor_(readScalar(is)), - centralWeight_(readScalar(is)) + centralWeight_(1000) {} @@ -123,13 +119,13 @@ public: { const fvMesh& mesh = this->mesh(); - const extendedCentredStencil& stencil = - centredCFCStencilObject::New + const extendedCentredStencil& stencil = Stencil::New ( mesh ); - const quadraticFitData& cfd = quadraticFitData::New + const CentredFitData& cfd = + CentredFitData::New ( mesh, stencil, @@ -137,7 +133,7 @@ public: centralWeight_ ); - const List& f = cfd.fit(); + const List& f = cfd.coeffs(); return stencil.weightedSum(vf, f); } @@ -148,6 +144,34 @@ public: } // End namespace Foam +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// Add the patch constructor functions to the hash tables + +#define makeCentredFitSurfaceInterpolationTypeScheme(SS, POLYNOMIAL, STENCIL, TYPE) \ + \ +typedef CentredFitScheme \ + CentredFitScheme##TYPE##POLYNOMIAL##STENCIL##_; \ +defineTemplateTypeNameAndDebugWithName \ + (CentredFitScheme##TYPE##POLYNOMIAL##STENCIL##_, #SS, 0); \ + \ +surfaceInterpolationScheme::addMeshConstructorToTable \ + > \ + add##SS##STENCIL##TYPE##MeshConstructorToTable_; \ + \ +surfaceInterpolationScheme::addMeshFluxConstructorToTable \ + > \ + add##SS##STENCIL##TYPE##MeshFluxConstructorToTable_; + +#define makeCentredFitSurfaceInterpolationScheme(SS, POLYNOMIAL, STENCIL) \ + \ +makeCentredFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,scalar) \ +makeCentredFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,vector) \ +makeCentredFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,sphericalTensor) \ +makeCentredFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,symmTensor)\ +makeCentredFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,tensor) + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #endif diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.C index a7dead621f..beff7881b2 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.C @@ -24,13 +24,26 @@ License \*---------------------------------------------------------------------------*/ -#include "linearFit.H" +#include "CentredFitScheme.H" +#include "linearFitPolynomial.H" +#include "centredCFCStencilObject.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { - makeSurfaceInterpolationScheme(linearFit); + defineTemplateTypeNameAndDebug + ( + CentredFitData, + 0 + ); + + makeCentredFitSurfaceInterpolationScheme + ( + linearFit, + linearFitPolynomial, + centredCFCStencilObject + ); } // ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.H deleted file mode 100644 index 882f3e94bc..0000000000 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.H +++ /dev/null @@ -1,138 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / 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 - linearFit - -Description - Linear fit interpolation scheme which applies an explicit correction to - linear. - -SourceFiles - linearFit.C - -\*---------------------------------------------------------------------------*/ - -#ifndef linearFit_H -#define linearFit_H - -#include "linear.H" -#include "linearFitData.H" -#include "extendedStencil.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -/*---------------------------------------------------------------------------*\ - Class linearFit Declaration -\*---------------------------------------------------------------------------*/ - -template -class linearFit -: - public linear -{ - // Private Data - const scalar centralWeight_; - - // Private Member Functions - - //- Disallow default bitwise copy construct - linearFit(const linearFit&); - - //- Disallow default bitwise assignment - void operator=(const linearFit&); - - -public: - - //- Runtime type information - TypeName("linearFit"); - - - // Constructors - - //- Construct from mesh and Istream - linearFit(const fvMesh& mesh, Istream& is) - : - linear(mesh), - centralWeight_(readScalar(is)) - {} - - - //- Construct from mesh, faceFlux and Istream - linearFit - ( - const fvMesh& mesh, - const surfaceScalarField& faceFlux, - Istream& is - ) - : - linear(mesh), - centralWeight_(readScalar(is)) - {} - - - // 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 > - correction - ( - const GeometricField& vf - ) const - { - const fvMesh& mesh = this->mesh(); - - const linearFitData& cfd = linearFitData::New - ( - mesh, - centralWeight_ - ); - - const extendedStencil& stencil = cfd.stencil(); - const List& f = cfd.fit(); - - return stencil.weightedSum(vf, f); - } -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitData.C deleted file mode 100644 index 996d915e1c..0000000000 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitData.C +++ /dev/null @@ -1,371 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / 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 "linearFitData.H" -#include "surfaceFields.H" -#include "volFields.H" -#include "SVD.H" -#include "syncTools.H" - - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -namespace Foam -{ - defineTypeNameAndDebug(linearFitData, 0); -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -static int count = 0; - -Foam::linearFitData::linearFitData -( - const fvMesh& mesh, - const scalar cWeight -) -: - MeshObject(mesh), - centralWeight_(cWeight), -# ifdef SPHERICAL_GEOMETRY - dim_(2), -# else - dim_(mesh.nGeometricD()), -# endif - minSize_ - ( - dim_ == 1 ? 2 : - dim_ == 2 ? 3 : - dim_ == 3 ? 4 : 0 - ), - stencil_(mesh), - fit_(mesh.nInternalFaces()) -{ - if (debug) - { - Info << "Contructing linearFitData" << endl; - } - - // check input - if (centralWeight_ < 1 - SMALL) - { - FatalErrorIn("linearFitData::linearFitData") - << "centralWeight requested = " << centralWeight_ - << " should not be less than one" - << exit(FatalError); - } - - if (minSize_ == 0) - { - FatalErrorIn("linearFitSnGradData") - << " dimension must be 1,2 or 3, not" << dim_ << exit(FatalError); - } - - // store the polynomial size for each cell to write out - surfaceScalarField interpPolySize - ( - IOobject - ( - "linearFitInterpPolySize", - "constant", - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh, - dimensionedScalar("linearFitInterpPolySize", dimless, scalar(0)) - ); - - // Get the cell/face centres in stencil order. - // Centred face stencils no good for triangles of tets. Need bigger stencils - List > 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++) - { - interpPolySize[faci] = calcFit(stencilPoints[faci], faci); - } - - Pout<< "count = " << count << endl; - - if (debug) - { - Info<< "linearFitData::linearFitData() :" - << "Finished constructing polynomialFit data" - << endl; - - interpPolySize.write(); - } -} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -void Foam::linearFitData::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 = mesh.C()[mesh.neighbour()[faci]] - mesh.C()[mesh.owner()[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::linearFitData::calcFit -( - const List& 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; -# ifndef SPHERICAL_GEOMETRY - scalar pz = (p - p0)&kdir; -# else - scalar pz = mag(p) - mag(p0); -# 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; - - if (dim_ >= 2) - { - B[ip][is++] = wts[ip]*py; - } - if (dim_ == 3) - { - B[ip][is++] = wts[ip]*pz; - } - } - - // Set the fit - label stencilSize = C.size(); - fit_[faci].setSize(stencilSize); - scalarList singVals(minSize_); - label nSVDzeros = 0; - - const GeometricField& w = - mesh().surfaceInterpolation::weights(); - - bool goodFit = false; - for(int iIt = 0; iIt < 10 && !goodFit; iIt++) - { - SVD svd(B, SMALL); - - scalar fit0 = wts[0]*wts[0]*svd.VSinvUt()[0][0]; - scalar fit1 = wts[0]*wts[1]*svd.VSinvUt()[0][1]; - - //goodFit = (fit0 > 0 && fit1 > 0); - - goodFit = - (mag(fit0 - w[faci])/w[faci] < 0.5) - && (mag(fit1 - (1 - w[faci]))/(1 - w[faci]) < 0.5); - - //scalar w0Err = fit0/w[faci]; - //scalar w1Err = fit1/(1 - w[faci]); - - //goodFit = - // (w0Err > 0.5 && w0Err < 1.5) - // && (w1Err > 0.5 && w1Err < 1.5); - - if (goodFit) - { - fit_[faci][0] = fit0; - fit_[faci][1] = fit1; - - for(label i=2; i -{ - // Private data - - //- weights for central stencil - const scalar centralWeight_; - - //- dimensionality of the geometry - 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 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&, const label faci); - - -public: - - TypeName("linearFitData"); - - - // Constructors - - explicit linearFitData - ( - const fvMesh& mesh, - scalar cWeightDim - ); - - - // Destructor - - virtual ~linearFitData() - {} - - - // Member functions - - - //- Return reference to the stencil - const extendedStencil& stencil() const - { - return stencil_; - } - - //- Return reference to fit coefficients - const List& 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/interpolation/surfaceInterpolation/schemes/linearFit/linearFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitPolynomial.H new file mode 100644 index 0000000000..70275f9ee1 --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitPolynomial.H @@ -0,0 +1,98 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 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::linearFitPolynomial + +Description + Linear polynomial for interpolation fitting. + Can be used with the CentredFit scheme to crate a linear surface + interpolation scheme + +\*---------------------------------------------------------------------------*/ + +#ifndef linearFitPolynomial_H +#define linearFitPolynomial_H + +#include "vector.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class linearFitPolynomial Declaration +\*---------------------------------------------------------------------------*/ + +class linearFitPolynomial +{ +public: + + // Member functions + + static label nTerms(const direction dim) + { + return + ( + dim == 1 ? 2 : + dim == 2 ? 3 : + dim == 3 ? 4 : 0 + ); + } + + static void addCoeffs + ( + scalar* coeffs, + const vector& d, + const scalar weight, + const direction dim + ) + { + register label curIdx = 0; + + coeffs[curIdx++] = weight; + coeffs[curIdx++] = weight*d.x(); + + if (dim >= 2) + { + coeffs[curIdx++] = weight*d.y(); + } + if (dim == 3) + { + coeffs[curIdx++] = weight*d.z(); + } + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFit.C similarity index 78% rename from src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C rename to src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFit.C index a215576330..0ecd5ece93 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFit.C @@ -24,13 +24,26 @@ License \*---------------------------------------------------------------------------*/ -#include "quadraticFit.H" +#include "CentredFitScheme.H" +#include "quadraticLinearFitPolynomial.H" +#include "centredCFCStencilObject.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { - makeSurfaceInterpolationScheme(quadraticFit); + defineTemplateTypeNameAndDebug + ( + CentredFitData, + 0 + ); + + makeCentredFitSurfaceInterpolationScheme + ( + quadraticLinearFit, + quadraticLinearFitPolynomial, + centredCFCStencilObject + ); } // ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFitPolynomial.H new file mode 100644 index 0000000000..4250e71c1e --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFitPolynomial.H @@ -0,0 +1,104 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 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::quadraticLinearFitPolynomial + +Description + Quadratic/linear polynomial for interpolation fitting: + quadratic normal to the face, + linear in the plane of the face for consistency with 2nd-order Gauss. + + Can be used with the CentredFit scheme to crate a quadratic surface + interpolation scheme + +\*---------------------------------------------------------------------------*/ + +#ifndef quadraticLinearFitPolynomial_H +#define quadraticLinearFitPolynomial_H + +#include "vector.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class quadraticLinearFitPolynomial Declaration +\*---------------------------------------------------------------------------*/ + +class quadraticLinearFitPolynomial +{ +public: + + // Member functions + + static label nTerms(const direction dim) + { + return + ( + dim == 1 ? 3 : + dim == 2 ? 5 : + dim == 3 ? 7 : 0 + ); + } + + static void addCoeffs + ( + scalar* coeffs, + const vector& d, + const scalar weight, + const direction dim + ) + { + register label curIdx = 0; + + coeffs[curIdx++] = weight; + coeffs[curIdx++] = weight*d.x(); + coeffs[curIdx++] = weight*sqr(d.x()); + + if (dim >= 2) + { + coeffs[curIdx++] = weight*d.y(); + coeffs[curIdx++] = weight*d.x()*d.y(); + } + if (dim == 3) + { + coeffs[curIdx++] = weight*d.z(); + coeffs[curIdx++] = weight*d.x()*d.z(); + } + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //