From 6c36054c43fcbf44ce5e1957fb9785e83a7eb1d8 Mon Sep 17 00:00:00 2001 From: henry <Henry Weller h.weller@opencfd.co.uk> Date: Fri, 30 Jan 2009 15:05:33 +0000 Subject: [PATCH] Added new centred and upwinded polynomial-fit schemes. --- src/finiteVolume/Make/files | 7 + .../extendedUpwindStencilTemplates.C | 2 +- .../schemes/CentredFitScheme/CentredFitData.C | 273 +---------------- .../schemes/CentredFitScheme/CentredFitData.H | 43 +-- .../schemes/FitData/FitData.C | 287 ++++++++++++++++++ .../schemes/FitData/FitData.H | 154 ++++++++++ .../schemes/UpwindFitScheme/UpwindFitData.C | 140 +++++++++ .../schemes/UpwindFitScheme/UpwindFitData.H | 128 ++++++++ .../schemes/UpwindFitScheme/UpwindFitScheme.H | 188 ++++++++++++ .../schemes/biLinearFit/biLinearFit.C | 49 +++ .../biLinearFit/biLinearFitPolynomial.H | 100 ++++++ .../schemes/cubicUpwindFit/cubicUpwindFit.C | 49 +++ .../cubicUpwindFit/cubicUpwindFitPolynomial.H | 109 +++++++ .../schemes/quadraticFit/quadraticFit.C | 49 +++ .../quadraticFit/quadraticFitPolynomial.H | 104 +++++++ .../quadraticLinearUpwindFit.C | 49 +++ .../quadraticLinearUpwindFitPolynomial.H | 102 +++++++ .../quadraticUpwindFit/quadraticUpwindFit.C | 49 +++ .../quadraticUpwindFitPolynomial.H | 105 +++++++ 19 files changed, 1685 insertions(+), 302 deletions(-) create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitData.C create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitData.H create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitScheme.H create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/biLinearFit/biLinearFit.C create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/biLinearFit/biLinearFitPolynomial.H create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubicUpwindFit/cubicUpwindFit.C create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubicUpwindFit/cubicUpwindFitPolynomial.H create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitPolynomial.H create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFit.C create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFitPolynomial.H create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticUpwindFit/quadraticUpwindFit.C create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticUpwindFit/quadraticUpwindFitPolynomial.H diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 85471d1ce9c..c6a95a976bd 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -187,8 +187,15 @@ $(schemes)/harmonic/harmonic.C $(schemes)/localBlended/localBlended.C $(schemes)/localMax/localMax.C $(schemes)/localMin/localMin.C + $(schemes)/linearFit/linearFit.C +$(schemes)/biLinearFit/biLinearFit.C $(schemes)/quadraticLinearFit/quadraticLinearFit.C +$(schemes)/quadraticFit/quadraticFit.C + +$(schemes)/quadraticLinearUpwindFit/quadraticLinearUpwindFit.C +$(schemes)/quadraticUpwindFit/quadraticUpwindFit.C +$(schemes)/cubicUpwindFit/cubicUpwindFit.C limitedSchemes = $(surfaceInterpolation)/limitedSchemes $(limitedSchemes)/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationSchemes.C diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencilTemplates.C b/src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencilTemplates.C index 00cd58e17d7..5c183290444 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencilTemplates.C +++ b/src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencilTemplates.C @@ -108,7 +108,7 @@ Foam::extendedUpwindStencil::weightedSum forAll(pSfCorr, i) { - if (phi[faceI] > 0) + if (phi.boundaryField()[patchi][i] > 0) { // Flux out of owner. Use upwind (= owner side) stencil. const List<Type>& stField = ownFld[faceI]; diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C index dea1b411437..7f684998447 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C @@ -42,16 +42,10 @@ Foam::CentredFitData<Polynomial>::CentredFitData const scalar centralWeight ) : - MeshObject<fvMesh, CentredFitData<Polynomial> >(mesh), - stencil_(stencil), - linearLimitFactor_(linearLimitFactor), - centralWeight_(centralWeight), -# ifdef SPHERICAL_GEOMETRY - dim_(2), -# else - dim_(mesh.nGeometricD()), -# endif - minSize_(Polynomial::nTerms(dim_)), + FitData<CentredFitData<Polynomial>, extendedCentredStencil, Polynomial> + ( + mesh, stencil, linearLimitFactor, centralWeight + ), coeffs_(mesh.nFaces()) { if (debug) @@ -59,15 +53,6 @@ Foam::CentredFitData<Polynomial>::CentredFitData Info<< "Contructing CentredFitData<Polynomial>" << endl; } - // Check input - if (linearLimitFactor > 1) - { - FatalErrorIn("CentredFitData<Polynomial>::CentredFitData") - << "linearLimitFactor requested = " << linearLimitFactor - << " should not be less than one" - << exit(FatalError); - } - calcFit(); if (debug) @@ -81,85 +66,25 @@ Foam::CentredFitData<Polynomial>::CentredFitData // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template<class Polynomial> -void Foam::CentredFitData<Polynomial>::findFaceDirs -( - vector& idir, // value changed in return - vector& jdir, // value changed in return - vector& kdir, // value changed in return - const fvMesh& mesh, - const label facei -) -{ - idir = mesh.faceAreas()[facei]; - idir /= mag(idir); - -# ifndef SPHERICAL_GEOMETRY - if (mesh.nGeometricD() <= 2) // find the normal direction - { - 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 plane of the face - { - const face& f = mesh.faces()[facei]; - kdir = mesh.points()[f[0]] - mesh.faceCentres()[facei]; - } -# else - // Spherical geometry so kdir is the radial direction - kdir = mesh.faceCentres()[facei]; -# 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; -} - - template<class Polynomial> void Foam::CentredFitData<Polynomial>::calcFit() { const fvMesh& mesh = this->mesh(); // Get the cell/face centres in stencil order. - // Centred face stencils no good for triangles of tets. + // Centred face stencils no good for triangles or tets. // Need bigger stencils List<List<point> > stencilPoints(mesh.nFaces()); - stencil_.collectData(mesh.C(), stencilPoints); + this->stencil().collectData(mesh.C(), stencilPoints); // find the fit coefficients for every face in the mesh - const surfaceScalarField& w = this->mesh().surfaceInterpolation::weights(); + const surfaceScalarField& w = mesh.surfaceInterpolation::weights(); for(label facei = 0; facei < mesh.nInternalFaces(); facei++) { - calcFit(stencilPoints[facei], w[facei], facei); + FitData<CentredFitData<Polynomial>, extendedCentredStencil, Polynomial>:: + calcFit(coeffs_[facei], stencilPoints[facei], w[facei], facei); } const surfaceScalarField::GeometricBoundaryField& bw = w.boundaryField(); @@ -174,7 +99,10 @@ void Foam::CentredFitData<Polynomial>::calcFit() forAll(pw, i) { - calcFit(stencilPoints[facei], pw[i], facei); + FitData + < + CentredFitData<Polynomial>, extendedCentredStencil, Polynomial + >::calcFit(coeffs_[facei], stencilPoints[facei], pw[i], facei); facei++; } } @@ -182,179 +110,4 @@ void Foam::CentredFitData<Polynomial>::calcFit() } -template<class Polynomial> -Foam::label Foam::CentredFitData<Polynomial>::calcFit -( - const List<point>& C, - const scalar wLin, - const label facei -) -{ - vector idir(1,0,0); - vector jdir(0,1,0); - vector kdir(0,0,1); - findFaceDirs(idir, jdir, kdir, this->mesh(), facei); - - // Setup the point weights - scalarList wts(C.size(), scalar(1)); - wts[0] = centralWeight_; - wts[1] = centralWeight_; - - // Reference point - point p0 = this->mesh().faceCentres()[facei]; - - // 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]; - - d.x() = (p - p0)&idir; - d.y() = (p - p0)&jdir; -# ifndef SPHERICAL_GEOMETRY - d.z() = (p - p0)&kdir; -# else - d.z() = mag(p) - mag(p0); -# endif - - if (ip == 0) - { - scale = cmptMax(cmptMag((d))); - } - - // Scale the radius vector - d /= scale; - - Polynomial::addCoeffs - ( - B[ip], - d, - wts[ip], - dim_ - ); - } - - // Set the fit - label stencilSize = C.size(); - coeffs_[facei].setSize(stencilSize); - scalarList singVals(minSize_); - label nSVDzeros = 0; - - bool goodFit = false; - for(int iIt = 0; iIt < 10 && !goodFit; iIt++) - { - SVD svd(B, SMALL); - - scalar fit0 = wts[0]*svd.VSinvUt()[0][0]; - scalar fit1 = wts[1]*svd.VSinvUt()[0][1]; - - goodFit = - (mag(fit0 - wLin) < linearLimitFactor_*wLin) - && (mag(fit1 - (1 - wLin)) < linearLimitFactor_*(1 - wLin)); - - //scalar w0Err = fit0/wLin; - //scalar w1Err = fit1/(1 - wLin); - - //goodFit = - // (w0Err > linearLimitFactor_ && w0Err < (1 + linearLimitFactor_)) - // && (w1Err > linearLimitFactor_ && w1Err < (1 + linearLimitFactor_)); - - if (goodFit) - { - coeffs_[facei][0] = fit0; - coeffs_[facei][1] = fit1; - - for(label i=2; i<stencilSize; i++) - { - coeffs_[facei][i] = wts[i]*svd.VSinvUt()[0][i]; - } - - 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 j = 0; j < B.m(); j++) - { - B[0][j] *= 10; - B[1][j] *= 10; - } - } - } - - // static const scalar L = 0.1; - // static const scalar R = 0.2; - - // static const scalar beta = 1.0/(R - L); - // static const scalar alpha = R*beta; - - if (goodFit) - { - // scalar limiter = - // max - // ( - // min - // ( - // min(alpha - beta*mag(coeffs_[facei][0] - wLin)/wLin, 1), - // min(alpha - beta*mag(coeffs_[facei][1] - (1 - wLin)) - // /(1 - wLin), 1) - // ), 0 - // ); - - //Info<< wLin << " " << coeffs_[facei][0] - // << " " << (1 - wLin) << " " << coeffs_[facei][1] << endl; - - // Remove the uncorrected linear ocefficients - coeffs_[facei][0] -= wLin; - coeffs_[facei][1] -= 1 - wLin; - - // if (limiter < 0.99) - // { - // for(label i = 0; i < stencilSize; i++) - // { - // coeffs_[facei][i] *= limiter; - // } - // } - } - else - { - if (debug) - { - WarningIn - ( - "CentredFitData<Polynomial>::calcFit" - "(const List<point>& C, const label facei" - ) << "Could not fit face " << facei - << ", reverting to linear." << nl - << " Weights " - << coeffs_[facei][0] << " " << wLin << nl - << " Linear weights " - << coeffs_[facei][1] << " " << 1 - wLin << endl; - } - - coeffs_[facei] = 0; - } - - return minSize_ - nSVDzeros; -} - - -template<class Polynomial> -bool Foam::CentredFitData<Polynomial>::movePoints() -{ - calcFit(); - return true; -} - - // ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.H index 1630e06a927..cb01840260d 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.H @@ -36,8 +36,7 @@ SourceFiles #ifndef CentredFitData_H #define CentredFitData_H -#include "MeshObject.H" -#include "fvMesh.H" +#include "FitData.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -53,27 +52,10 @@ class extendedCentredStencil; template<class Polynomial> class CentredFitData : - public MeshObject<fvMesh, CentredFitData<Polynomial> > + public FitData<CentredFitData<Polynomial>, extendedCentredStencil, Polynomial> { // Private data - //- The stencil the fit is based on - const extendedCentredStencil& stencil_; - - //- 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_; - - //- Dimensionality of the geometry - const label dim_; - - //- Minimum stencil size - const label minSize_; - //- 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> coeffs_; @@ -81,28 +63,10 @@ class CentredFitData // 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 - ); - //- Calculate the fit for the all the mesh faces // and set the coefficients void calcFit(); - //- Calculate the fit for the specified face and set the coefficients - label calcFit - ( - const List<point>&, // Stencil points - const scalar wLin, // Linear weight - const label faci // Current face index - ); - public: @@ -133,9 +97,6 @@ public: { return coeffs_; } - - //- Delete the data when the mesh moves not implemented - virtual bool movePoints(); }; diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C new file mode 100644 index 00000000000..0872cb81a9c --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C @@ -0,0 +1,287 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "FitData.H" +#include "surfaceFields.H" +#include "volFields.H" +#include "SVD.H" +#include "syncTools.H" +#include "extendedStencil.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template<class Form, class extendedStencil, class Polynomial> +Foam::FitData<Form, extendedStencil, Polynomial>::FitData +( + const fvMesh& mesh, + const extendedStencil& stencil, + const scalar linearLimitFactor, + const scalar centralWeight +) +: + MeshObject<fvMesh, Form>(mesh), + stencil_(stencil), + linearLimitFactor_(linearLimitFactor), + centralWeight_(centralWeight), +# ifdef SPHERICAL_GEOMETRY + dim_(2), +# else + dim_(mesh.nGeometricD()), +# endif + minSize_(Polynomial::nTerms(dim_)) +{ + // Check input + if (linearLimitFactor <= SMALL || linearLimitFactor > 3) + { + FatalErrorIn("FitData<Polynomial>::FitData") + << "linearLimitFactor requested = " << linearLimitFactor + << " should be between zero and 3" + << exit(FatalError); + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class FitDataType, class ExtendedStencil, class Polynomial> +void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::findFaceDirs +( + vector& idir, // value changed in return + vector& jdir, // value changed in return + vector& kdir, // value changed in return + const label facei +) +{ + const fvMesh& mesh = this->mesh(); + + idir = mesh.faceAreas()[facei]; + idir /= mag(idir); + +# ifndef SPHERICAL_GEOMETRY + if (mesh.nGeometricD() <= 2) // find the normal direction + { + 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 plane of the face + { + const face& f = mesh.faces()[facei]; + kdir = mesh.points()[f[0]] - mesh.faceCentres()[facei]; + } +# else + // Spherical geometry so kdir is the radial direction + kdir = mesh.faceCentres()[facei]; +# 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; +} + + +template<class FitDataType, class ExtendedStencil, class Polynomial> +void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit +( + scalarList& coeffsi, + const List<point>& C, + const scalar wLin, + const label facei +) +{ + vector idir(1,0,0); + vector jdir(0,1,0); + vector kdir(0,0,1); + findFaceDirs(idir, jdir, kdir, facei); + + // Setup the point weights + scalarList wts(C.size(), scalar(1)); + wts[0] = centralWeight_; + wts[1] = centralWeight_; + + // Reference point + point p0 = this->mesh().faceCentres()[facei]; + + // Info << "Face " << facei << " at " << p0 << " stencil points at:\n" + // << C - p0 << endl; + + // 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]; + + d.x() = (p - p0)&idir; + d.y() = (p - p0)&jdir; +# ifndef SPHERICAL_GEOMETRY + d.z() = (p - p0)&kdir; +# else + d.z() = mag(p) - mag(p0); +# endif + + if (ip == 0) + { + scale = cmptMax(cmptMag((d))); + } + + // Scale the radius vector + d /= scale; + + Polynomial::addCoeffs + ( + B[ip], + d, + wts[ip], + dim_ + ); + } + + // Set the fit + label stencilSize = C.size(); + coeffsi.setSize(stencilSize); + + bool goodFit = false; + for(int iIt = 0; iIt < 8 && !goodFit; iIt++) + { + SVD svd(B, SMALL); + + scalar maxCoeff = 0; + label maxCoeffi = 0; + + for(label i=0; i<stencilSize; i++) + { + coeffsi[i] = wts[i]*svd.VSinvUt()[0][i]; + if (mag(coeffsi[i]) > maxCoeff) + { + maxCoeff = mag(coeffsi[i]); + maxCoeffi = i; + } + } + + goodFit = + (mag(coeffsi[0] - wLin) < linearLimitFactor_*wLin) + && (mag(coeffsi[1] - (1 - wLin)) < linearLimitFactor_*(1 - wLin)) + && maxCoeffi <= 1; + + // if (goodFit && iIt > 0) + // { + // Info << "FitData<Polynomial>::calcFit" + // << "(const List<point>& C, const label facei" << nl + // << "Can now fit face " << facei << " iteration " << iIt + // << " with sum of weights " << sum(coeffsi) << nl + // << " Weights " << coeffsi << nl + // << " Linear weights " << wLin << " " << 1 - wLin << nl + // << " sing vals " << svd.S() << endl; + // } + + if (!goodFit) // (not good fit so increase weight in the centre) + { + // if (iIt == 7) + // { + // WarningIn + // ( + // "FitData<Polynomial>::calcFit" + // "(const List<point>& C, const label facei" + // ) << "Cannot fit face " << facei + // << " sing vals " << svd.S() << endl; + // } + + wts[0] *= 10; + wts[1] *= 10; + + for(label j = 0; j < B.m(); j++) + { + B[0][j] *= 10; + B[1][j] *= 10; + } + } + } + + if (goodFit) + { + // Remove the uncorrected linear ocefficients + coeffsi[0] -= wLin; + coeffsi[1] -= 1 - wLin; + } + else + { + // if (debug) + // { + WarningIn + ( + "FitData<Polynomial>::calcFit" + "(const List<point>& C, const label facei" + ) << "Could not fit face " << facei + << " Weights = " << coeffsi + << ", reverting to linear." << nl + << " Linear weights " << wLin << " " << 1 - wLin << endl; + // } + + coeffsi = 0; + } +} + + +template<class FitDataType, class ExtendedStencil, class Polynomial> +bool Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::movePoints() +{ + calcFit(); + return true; +} + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H new file mode 100644 index 00000000000..810c8d7e4c3 --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H @@ -0,0 +1,154 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + Foam::FitData + +Description + Data for the upwinded and centred polynomial fit interpolation schemes + +SourceFiles + FitData.C + +\*---------------------------------------------------------------------------*/ + +#ifndef FitData_H +#define FitData_H + +#include "MeshObject.H" +#include "fvMesh.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class FitData Declaration +\*---------------------------------------------------------------------------*/ + +template<class FitDataType, class ExtendedStencil, class Polynomial> +class FitData +: + public MeshObject<fvMesh, FitDataType> +{ + // Private data + + //- The stencil the fit is based on + const ExtendedStencil& stencil_; + + //- 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_; + + //- Dimensionality of the geometry + const label dim_; + + //- Minimum stencil size + const label minSize_; + + + // Private member functions + + //- Find the normal direction (i) and j and k directions for face faci + void findFaceDirs + ( + vector& idir, // value changed in return + vector& jdir, // value changed in return + vector& kdir, // value changed in return + const label faci + ); + + //- Calculate the fit for the all the mesh faces + // and set the coefficients + // virtual void calcFit(); + + +public: + + //TypeName("FitData"); + + + // Constructors + + //- Construct from components + FitData + ( + const fvMesh& mesh, + const ExtendedStencil& stencil, + const scalar linearLimitFactor, + const scalar centralWeight + ); + + + //- Destructor + virtual ~FitData() + {} + + + // Member functions + + //- Return reference to the stencil + const ExtendedStencil& stencil() const + { + return stencil_; + } + + //- 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 label faci // Current face index + ); + + //- Calculate the fit for all the faces + virtual void calcFit() = 0; + + + //- Delete the data when the mesh moves not implemented + bool movePoints(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "FitData.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitData.C new file mode 100644 index 00000000000..fd218100cce --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitData.C @@ -0,0 +1,140 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "UpwindFitData.H" +#include "surfaceFields.H" +#include "volFields.H" +#include "SVD.H" +#include "syncTools.H" +#include "extendedUpwindStencil.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template<class Polynomial> +Foam::UpwindFitData<Polynomial>::UpwindFitData +( + const fvMesh& mesh, + const extendedUpwindStencil& stencil, + const scalar linearLimitFactor, + const scalar centralWeight +) +: + FitData<UpwindFitData<Polynomial>, extendedUpwindStencil, Polynomial> + ( + mesh, stencil, linearLimitFactor, centralWeight + ), + owncoeffs_(mesh.nFaces()), + neicoeffs_(mesh.nFaces()) +{ + if (debug) + { + Info<< "Contructing UpwindFitData<Polynomial>" << endl; + } + + calcFit(); + + if (debug) + { + Info<< "UpwindFitData<Polynomial>::UpwindFitData() :" + << "Finished constructing polynomialFit data" + << endl; + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class Polynomial> +void Foam::UpwindFitData<Polynomial>::calcFit() +{ + const fvMesh& mesh = this->mesh(); + + // 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()); + 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 + ); + + // find the fit coefficients for every owner and neighbour of ever face + + const surfaceScalarField& w = mesh.surfaceInterpolation::weights(); + + for(label facei = 0; facei < mesh.nInternalFaces(); facei++) + { + FitData<UpwindFitData<Polynomial>, extendedUpwindStencil, Polynomial>:: + calcFit(owncoeffs_[facei], ownStencilPoints[facei], w[facei], facei); + FitData<UpwindFitData<Polynomial>, extendedUpwindStencil, Polynomial>:: + calcFit(neicoeffs_[facei], neiStencilPoints[facei], w[facei], facei); + } + + const surfaceScalarField::GeometricBoundaryField& bw = w.boundaryField(); + + forAll(bw, patchi) + { + const fvsPatchScalarField& pw = bw[patchi]; + + if (pw.coupled()) + { + label facei = pw.patch().patch().start(); + + forAll(pw, i) + { + FitData + < + UpwindFitData<Polynomial>, extendedUpwindStencil, Polynomial + >::calcFit + ( + owncoeffs_[facei], ownStencilPoints[facei], pw[i], facei + ); + FitData + < + UpwindFitData<Polynomial>, extendedUpwindStencil, Polynomial + >::calcFit + ( + neicoeffs_[facei], neiStencilPoints[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 new file mode 100644 index 00000000000..7f01135b9e6 --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitData.H @@ -0,0 +1,128 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + Foam::UpwindFitData + +Description + Data for the quadratic fit correction interpolation scheme to be used with + upwind biased stencil + +SourceFiles + UpwindFitData.C + +\*---------------------------------------------------------------------------*/ + +#ifndef UpwindFitData_H +#define UpwindFitData_H + +#include "FitData.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class extendedUpwindStencil; + +/*---------------------------------------------------------------------------*\ + Class UpwindFitData Declaration +\*---------------------------------------------------------------------------*/ + +template<class Polynomial> +class UpwindFitData +: + public FitData<UpwindFitData<Polynomial>, extendedUpwindStencil, Polynomial> +{ + // Private data + + //- For each face of the mesh store the coefficients to multiply the + // stencil cell values by if the flow is from the owner + List<scalarList> owncoeffs_; + + //- For each face of the mesh store the coefficients to multiply the + // stencil cell values by if the flow is from the neighbour + List<scalarList> neicoeffs_; + + + // Private member functions + + //- Calculate the fit for the all the mesh faces + // and set the coefficients + void calcFit(); + + +public: + + TypeName("UpwindFitData"); + + + // Constructors + + //- Construct from components + UpwindFitData + ( + const fvMesh& mesh, + const extendedUpwindStencil& stencil, + const scalar linearLimitFactor, + const scalar centralWeight + ); + + + //- Destructor + virtual ~UpwindFitData() + {} + + + // Member functions + + //- Return reference to owner fit coefficients + const List<scalarList>& owncoeffs() const + { + return owncoeffs_; + } + + //- Return reference to neighbour fit coefficients + const List<scalarList>& neicoeffs() const + { + return neicoeffs_; + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "UpwindFitData.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitScheme.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitScheme.H new file mode 100644 index 00000000000..0921dcfea96 --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/UpwindFitScheme/UpwindFitScheme.H @@ -0,0 +1,188 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + Foam::UpwindFitScheme + +Description + Upwind biased fit surface interpolation scheme which applies an explicit + correction to linear. + +\*---------------------------------------------------------------------------*/ + +#ifndef UpwindFitScheme_H +#define UpwindFitScheme_H + +#include "UpwindFitData.H" +#include "linear.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class UpwindFitScheme Declaration +\*---------------------------------------------------------------------------*/ + +template<class Type, class Polynomial, class Stencil> +class UpwindFitScheme +: + public linear<Type> +{ + // Private Data + + //- Reference to the surface flux used to choose upwind direction + const surfaceScalarField& faceFlux_; + + //- 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 + UpwindFitScheme(const UpwindFitScheme&); + + //- Disallow default bitwise assignment + void operator=(const UpwindFitScheme&); + + +public: + + //- Runtime type information + TypeName("UpwindFitScheme"); + + + // Constructors + + //- Construct from mesh and Istream + // The name of the flux field is read from the Istream and looked-up + // from the mesh objectRegistry + UpwindFitScheme(const fvMesh& mesh, Istream& is) + : + linear<Type>(mesh), + faceFlux_(mesh.lookupObject<surfaceScalarField>(word(is))), + linearLimitFactor_(readScalar(is)), + centralWeight_(1000) + {} + + + //- Construct from mesh, faceFlux and Istream + UpwindFitScheme + ( + const fvMesh& mesh, + const surfaceScalarField& faceFlux, + Istream& is + ) + : + linear<Type>(mesh), + faceFlux_(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(); + + const extendedUpwindStencil& stencil = Stencil::New + ( + mesh, + scalar(0.5) + ); + + const UpwindFitData<Polynomial>& ufd = + UpwindFitData<Polynomial>::New + ( + mesh, + stencil, + linearLimitFactor_, + centralWeight_ + ); + + const List<scalarList>& fo = ufd.owncoeffs(); + const List<scalarList>& fn = ufd.neicoeffs(); + + return stencil.weightedSum(faceFlux_, vf, fo, fn); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// Add the patch constructor functions to the hash tables + +#define makeUpwindFitSurfaceInterpolationTypeScheme(SS, POLYNOMIAL, STENCIL, TYPE) \ + \ +typedef UpwindFitScheme<TYPE, POLYNOMIAL, STENCIL> \ + UpwindFitScheme##TYPE##POLYNOMIAL##STENCIL##_; \ +defineTemplateTypeNameAndDebugWithName \ + (UpwindFitScheme##TYPE##POLYNOMIAL##STENCIL##_, #SS, 0); \ + \ +surfaceInterpolationScheme<TYPE>::addMeshConstructorToTable \ +<UpwindFitScheme<TYPE, POLYNOMIAL, STENCIL> > \ + add##SS##STENCIL##TYPE##MeshConstructorToTable_; \ + \ +surfaceInterpolationScheme<TYPE>::addMeshFluxConstructorToTable \ +<UpwindFitScheme<TYPE, POLYNOMIAL, STENCIL> > \ + add##SS##STENCIL##TYPE##MeshFluxConstructorToTable_; + +#define makeUpwindFitSurfaceInterpolationScheme(SS, POLYNOMIAL, STENCIL) \ + \ +makeUpwindFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,scalar) \ +makeUpwindFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,vector) \ +makeUpwindFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,sphericalTensor) \ +makeUpwindFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,symmTensor)\ +makeUpwindFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,tensor) + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/biLinearFit/biLinearFit.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/biLinearFit/biLinearFit.C new file mode 100644 index 00000000000..8d49862a971 --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/biLinearFit/biLinearFit.C @@ -0,0 +1,49 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "CentredFitScheme.H" +#include "biLinearFitPolynomial.H" +#include "centredCFCStencilObject.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + defineTemplateTypeNameAndDebug + ( + CentredFitData<biLinearFitPolynomial>, + 0 + ); + + makeCentredFitSurfaceInterpolationScheme + ( + biLinearFit, + biLinearFitPolynomial, + centredCFCStencilObject + ); +} + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/biLinearFit/biLinearFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/biLinearFit/biLinearFitPolynomial.H new file mode 100644 index 00000000000..5960f54ee2b --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/biLinearFit/biLinearFitPolynomial.H @@ -0,0 +1,100 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::biLinearFitPolynomial + +Description + BiLinear polynomial for interpolation fitting. + Can be used with the CentredFit scheme to crate a biLinear surface + interpolation scheme + +\*---------------------------------------------------------------------------*/ + +#ifndef biLinearFitPolynomial_H +#define biLinearFitPolynomial_H + +#include "vector.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class biLinearFitPolynomial Declaration +\*---------------------------------------------------------------------------*/ + +class biLinearFitPolynomial +{ +public: + + // Member functions + + static label nTerms(const direction dim) + { + return + ( + dim == 1 ? 2 : + dim == 2 ? 4 : + dim == 3 ? 6 : 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(); + 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 + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubicUpwindFit/cubicUpwindFit.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubicUpwindFit/cubicUpwindFit.C new file mode 100644 index 00000000000..63f239f5295 --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubicUpwindFit/cubicUpwindFit.C @@ -0,0 +1,49 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "UpwindFitScheme.H" +#include "cubicUpwindFitPolynomial.H" +#include "upwindCFCStencilObject.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + defineTemplateTypeNameAndDebug + ( + UpwindFitData<cubicUpwindFitPolynomial>, + 0 + ); + + makeUpwindFitSurfaceInterpolationScheme + ( + cubicUpwindFit, + cubicUpwindFitPolynomial, + upwindCFCStencilObject + ); +} + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubicUpwindFit/cubicUpwindFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubicUpwindFit/cubicUpwindFitPolynomial.H new file mode 100644 index 00000000000..dab0871b8e4 --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubicUpwindFit/cubicUpwindFitPolynomial.H @@ -0,0 +1,109 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::cubicUpwindFitPolynomial + +Description + Cubic polynomial for upwind biased interpolation fitting. + + Can be used with the UpwindFit scheme to crate a cubic surface + interpolation scheme + +\*---------------------------------------------------------------------------*/ + +#ifndef cubicUpwindFitPolynomial_H +#define cubicUpwindFitPolynomial_H + +#include "vector.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class cubicUpwindFitPolynomial Declaration +\*---------------------------------------------------------------------------*/ + +class cubicUpwindFitPolynomial +{ +public: + + // Member functions + + static label nTerms(const direction dim) + { + return + ( + dim == 1 ? 4 : + dim == 2 ? 8 : + dim == 3 ? 14 : 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()); + coeffs[curIdx++] = weight*pow(d.x(),3); + + if (dim >= 2) + { + coeffs[curIdx++] = weight*d.y(); + coeffs[curIdx++] = weight*d.x()*d.y(); + coeffs[curIdx++] = weight*sqr(d.y()); + coeffs[curIdx++] = weight*d.x()*sqr(d.y()); + } + if (dim == 3) + { + coeffs[curIdx++] = weight*d.z(); + coeffs[curIdx++] = weight*d.x()*d.z(); + coeffs[curIdx++] = weight*d.y()*d.z(); + coeffs[curIdx++] = weight*sqr(d.z()); + coeffs[curIdx++] = weight*d.x()*d.y()*d.z(); + coeffs[curIdx++] = weight*d.x()*sqr(d.z()); + } + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C new file mode 100644 index 00000000000..fd0bd3f8bdf --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C @@ -0,0 +1,49 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "CentredFitScheme.H" +#include "quadraticFitPolynomial.H" +#include "centredCFCStencilObject.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + defineTemplateTypeNameAndDebug + ( + CentredFitData<quadraticFitPolynomial>, + 0 + ); + + makeCentredFitSurfaceInterpolationScheme + ( + quadraticFit, + quadraticFitPolynomial, + centredCFCStencilObject + ); +} + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitPolynomial.H new file mode 100644 index 00000000000..5ba7520df14 --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitPolynomial.H @@ -0,0 +1,104 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::quadraticFitPolynomial + +Description + Quadratic polynomial for centred interpolation fitting. + + Can be used with the CentredFit scheme to crate a quadratic surface + interpolation scheme + +\*---------------------------------------------------------------------------*/ + +#ifndef quadraticFitPolynomial_H +#define quadraticFitPolynomial_H + +#include "vector.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class quadraticFitPolynomial Declaration +\*---------------------------------------------------------------------------*/ + +class quadraticFitPolynomial +{ +public: + + // Member functions + + static label nTerms(const direction dim) + { + return + ( + dim == 1 ? 3 : + dim == 2 ? 6 : + dim == 3 ? 9 : 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(); + coeffs[curIdx++] = weight*sqr(d.y()); + } + if (dim == 3) + { + coeffs[curIdx++] = weight*d.z(); + coeffs[curIdx++] = weight*d.x()*d.z(); + coeffs[curIdx++] = weight*sqr(d.z()); + } + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFit.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFit.C new file mode 100644 index 00000000000..fb344602dca --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFit.C @@ -0,0 +1,49 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "UpwindFitScheme.H" +#include "quadraticLinearUpwindFitPolynomial.H" +#include "upwindCFCStencilObject.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + defineTemplateTypeNameAndDebug + ( + UpwindFitData<quadraticLinearUpwindFitPolynomial>, + 0 + ); + + makeUpwindFitSurfaceInterpolationScheme + ( + quadraticLinearUpwindFit, + quadraticLinearUpwindFitPolynomial, + upwindCFCStencilObject + ); +} + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFitPolynomial.H new file mode 100644 index 00000000000..b904d6baa82 --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFitPolynomial.H @@ -0,0 +1,102 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::quadraticLinearUpwindFitPolynomial + +Description + Quadratic polynomial for upwind biased interpolation fitting. + + Can be used with the UpwindFit scheme to crate a quadratic surface + interpolation scheme + +\*---------------------------------------------------------------------------*/ + +#ifndef quadraticLinearUpwindFitPolynomial_H +#define quadraticLinearUpwindFitPolynomial_H + +#include "vector.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class quadraticLinearUpwindFitPolynomial Declaration +\*---------------------------------------------------------------------------*/ + +class quadraticLinearUpwindFitPolynomial +{ +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 + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticUpwindFit/quadraticUpwindFit.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticUpwindFit/quadraticUpwindFit.C new file mode 100644 index 00000000000..66aadef7dbb --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticUpwindFit/quadraticUpwindFit.C @@ -0,0 +1,49 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "UpwindFitScheme.H" +#include "quadraticUpwindFitPolynomial.H" +#include "upwindFECStencilObject.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + defineTemplateTypeNameAndDebug + ( + UpwindFitData<quadraticUpwindFitPolynomial>, + 0 + ); + + makeUpwindFitSurfaceInterpolationScheme + ( + quadraticUpwindFit, + quadraticUpwindFitPolynomial, + upwindFECStencilObject + ); +} + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticUpwindFit/quadraticUpwindFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticUpwindFit/quadraticUpwindFitPolynomial.H new file mode 100644 index 00000000000..7256d601a38 --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticUpwindFit/quadraticUpwindFitPolynomial.H @@ -0,0 +1,105 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::quadraticUpwindFitPolynomial + +Description + Quadratic polynomial for upwind biased interpolation fitting. + + Can be used with the UpwindFit scheme to crate a quadratic surface + interpolation scheme + +\*---------------------------------------------------------------------------*/ + +#ifndef quadraticUpwindFitPolynomial_H +#define quadraticUpwindFitPolynomial_H + +#include "vector.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class quadraticUpwindFitPolynomial Declaration +\*---------------------------------------------------------------------------*/ + +class quadraticUpwindFitPolynomial +{ +public: + + // Member functions + + static label nTerms(const direction dim) + { + return + ( + dim == 1 ? 3 : + dim == 2 ? 6 : + dim == 3 ? 9 : 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(); + //coeffs[curIdx++] = weight*d.x()*sqr(d.y()); + coeffs[curIdx++] = weight*sqr(d.y()); + } + if (dim == 3) + { + coeffs[curIdx++] = weight*d.z(); + coeffs[curIdx++] = weight*d.x()*d.z(); + coeffs[curIdx++] = weight*sqr(d.z()); + } + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // -- GitLab