Skip to content
GitLab
Commit e30c626d authored by mattijs's avatar mattijs
Browse files

Merge branch 'master' of /home/hunt2/OpenFOAM/OpenFOAM-dev

parents 3c2f9b33 809bbd5a
Hide whitespace changes
Inline Side-by-side
src/finiteVolume/Make/files
View file @ e30c626d
......@@ -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
......
src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C
View file @ e30c626d
......@@ -128,6 +128,7 @@ Foam::extendedStencil::weightedSum
// Boundaries. Either constrained or calculated so assign value
// directly (instead of nicely using operator==)
/*
typename GeometricField<Type, fvsPatchField, surfaceMesh>::
GeometricBoundaryField& bSfCorr = sf.boundaryField();
......@@ -150,6 +151,7 @@ Foam::extendedStencil::weightedSum
faceI++;
}
}
*/
return tsfCorr;
}
......
src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C → src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C
View file @ e30c626d
......@@ -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<class Polynomial>
Foam::CentredFitData<Polynomial>::CentredFitData
(
const fvMesh& mesh,
const extendedCentredStencil& stencil,
const scalar linearLimitFactor,
const scalar cWeight
const scalar centralWeight
)
:
MeshObject<fvMesh, quadraticFitData>(mesh),
MeshObject<fvMesh, CentredFitData<Polynomial> >(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<Polynomial>" << endl;
}
// check input
if (centralWeight_ < 1 - SMALL)
{
FatalErrorIn("quadraticFitData::quadraticFitData")
FatalErrorIn("CentredFitData<Polynomial>::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<List<point> > 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<Polynomial>::CentredFitData() :"
<< "Finished constructing polynomialFit data"
<< endl;
......@@ -129,7 +113,8 @@ Foam::quadraticFitData::quadraticFitData
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::quadraticFitData::findFaceDirs
template<class Polynomial>
void Foam::CentredFitData<Polynomial>::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<class Polynomial>
Foam::label Foam::CentredFitData<Polynomial>::calcFit
(
const List<point>& 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];
// p0 -> p vector in the face-local coordinate system
vector d;
// calculate the matrix of the polynomial components
// 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;
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);
// Scale the radius vector
d /= scale;
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<scalar, fvsPatchField, surfaceMesh>& 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<stencilSize; i++)
{
fit_[faci][i] = wts[0]*wts[i]*svd.VSinvUt()[0][i];
coeffs_[faci][i] = wts[i]*svd.VSinvUt()[0][i];
}
singVals = svd.S();
......@@ -300,12 +279,6 @@ Foam::label Foam::quadraticFitData::calcFit
wts[0] *= 10;
wts[1] *= 10;
for(label i = 0; i < B.n(); i++)
{
B[i][0] *= 10;
B[i][1] *= 10;
}
for(label j = 0; j < B.m(); j++)
{
B[0][j] *= 10;
......@@ -327,39 +300,52 @@ Foam::label Foam::quadraticFitData::calcFit
// (
// min
// (
// min(alpha - beta*mag(fit_[faci][0] - w[faci])/w[faci], 1),
// min(alpha - beta*mag(fit_[faci][1] - (1 - w[faci]))/(1 - w[faci]), 1)
// min(alpha - beta*mag(coeffs_[faci][0] - w[faci])/w[faci], 1),
// min(alpha - beta*mag(coeffs_[faci][1] - (1 - w[faci]))/(1 - w[faci]), 1)
// ), 0
// );
//Info<< w[faci] << " " << coeffs_[faci][0] << " " << (1 - w[faci]) << " " << coeffs_[faci][1] << endl;
// Remove the uncorrected linear coefficients
fit_[faci][0] -= w[faci];
fit_[faci][1] -= 1 - w[faci];
coeffs_[faci][0] -= w[faci];
coeffs_[faci][1] -= 1 - w[faci];
// if (limiter < 0.99)
// {
// for(label i = 0; i < stencilSize; i++)
// {
// fit_[faci][i] *= limiter;
// coeffs_[faci][i] *= limiter;
// }
// }
}
else
{
Pout<< "Could not fit face " << faci
<< " " << fit_[faci][0] << " " << w[faci]
<< " " << fit_[faci][1] << " " << 1 - w[faci]<< endl;
if (debug)
{
WarningIn
(
"CentredFitData<Polynomial>::calcFit"
"(const List<point>& 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<class Polynomial>
bool Foam::CentredFitData<Polynomial>::movePoints()
{
notImplemented("quadraticFitData::movePoints()");
notImplemented("CentredFitData<Polynomial>::movePoints()");
return true;
}
......
src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H → src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.H
View file @ e30c626d
......@@ -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 Polynomial>
class CentredFitData
:
public MeshObject<fvMesh, quadraticFitData>
public MeshObject<fvMesh, CentredFitData<Polynomial> >
{
// 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<scalarList> fit_;
List<scalarList> 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<scalarList>& fit() const
const List<scalarList>& 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
// ************************************************************************* //
src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H → src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitScheme.H
View file @ e30c626d
......@@ -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 Type>
class quadraticFit
template<class Type, class Polynomial, class Stencil>
class CentredFitScheme
:
public linear<Type>
{
......@@ -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<Type>(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<Type>(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<Polynomial>& cfd =
CentredFitData<Polynomial>::New
(
mesh,
stencil,
......@@ -137,7 +133,7 @@ public:
centralWeight_
);
const List<scalarList>& f = cfd.fit();
const List<scalarList>& f = cfd.coeffs();