Commit 4354a9af authored by Henry's avatar Henry
Browse files

finiteVolume: non-orthogonal correction is now selected by the user and not overridden

Previously non-orthogonal correction was not applied if the mesh appeared to be
orthogonal based on a hard-coded test.  Now the non-orthogonal correction is
applied as selected by the user irrespective of the mesh quality.  This avoids
problems with the built-in mesh test being inappropriate for some cases.
parent d242f68b
......@@ -334,6 +334,7 @@ $(snGradSchemes)/snGradScheme/snGradSchemes.C
$(snGradSchemes)/correctedSnGrad/correctedSnGrads.C
$(snGradSchemes)/limitedSnGrad/limitedSnGrads.C
$(snGradSchemes)/uncorrectedSnGrad/uncorrectedSnGrads.C
$(snGradSchemes)/orthogonalSnGrad/orthogonalSnGrads.C
/*
$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGradData.C
$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGrads.C
......
......@@ -110,38 +110,34 @@ void Foam::extendedLeastSquaresVectors::makeLeastSquaresVectors() const
const labelUList& owner = mesh_.owner();
const labelUList& neighbour = mesh_.neighbour();
// Build the d-vectors
surfaceVectorField d
(
mesh_.Sf()/(mesh_.magSf()*mesh_.deltaCoeffs())
);
if (!mesh_.orthogonal())
{
d -= mesh_.correctionVectors()/mesh_.deltaCoeffs();
}
const volVectorField& C = mesh.C();
// Set up temporary storage for the dd tensor (before inversion)
symmTensorField dd(mesh_.nCells(), symmTensor::zero);
forAll(owner, faceI)
forAll(owner, facei)
{
const symmTensor wdd(1.0/magSqr(d[faceI])*sqr(d[faceI]));
label own = owner[facei];
label nei = neighbour[facei];
vector d = C[nei] - C[own];
dd[owner[faceI]] += wdd;
dd[neighbour[faceI]] += wdd;
const symmTensor wdd(1.0/magSqr(d[facei])*sqr(d[facei]));
dd[own] += wdd;
dd[nei] += wdd;
}
// Visit the boundaries. Coupled boundaries are taken into account
// in the construction of d vectors.
forAll(d.boundaryField(), patchI)
forAll(lsP.boundaryField(), patchi)
{
const fvsPatchVectorField& pd = d.boundaryField()[patchI];
const fvPatch& p = pd.patch();
const fvPatch& p = lsP.boundaryField()[patchi].patch();
const labelUList& faceCells = p.faceCells();
// Build the d-vectors
vectorField pd = p.delta();
forAll(pd, patchFaceI)
{
dd[faceCells[patchFaceI]] +=
......@@ -232,24 +228,30 @@ void Foam::extendedLeastSquaresVectors::makeLeastSquaresVectors() const
// Revisit all faces and calculate the lsP and lsN vectors
forAll(owner, faceI)
forAll(owner, facei)
{
lsP[faceI] =
(1.0/magSqr(d[faceI]))*(invDd[owner[faceI]] & d[faceI]);
label own = owner[facei];
label nei = neighbour[facei];
vector d = C[nei] - C[own];
lsN[faceI] =
((-1.0)/magSqr(d[faceI]))*(invDd[neighbour[faceI]] & d[faceI]);
lsP[facei] =
(1.0/magSqr(d[facei]))*(invDd[owner[facei]] & d);
lsN[facei] =
((-1.0)/magSqr(d[facei]))*(invDd[neighbour[facei]] & d);
}
forAll(lsP.boundaryField(), patchI)
{
const fvsPatchVectorField& pd = d.boundaryField()[patchI];
fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchI];
const fvPatch& p = patchLsP.patch();
const labelUList& faceCells = p.faceCells();
// Build the d-vectors
vectorField pd = p.delta();
forAll(p, patchFaceI)
{
patchLsP[patchFaceI] =
......
......@@ -121,24 +121,12 @@ Foam::fv::fourthGrad<Type>::calcGrad
const scalarField& lambdap = lambda.boundaryField()[patchi];
// Build the d-vectors
vectorField pd
(
mesh.Sf().boundaryField()[patchi]
/ (
mesh.magSf().boundaryField()[patchi]
* mesh.deltaCoeffs().boundaryField()[patchi]
)
);
const fvPatch& p = fGrad.boundaryField()[patchi].patch();
if (!mesh.orthogonal())
{
pd -= mesh.correctionVectors().boundaryField()[patchi]
/mesh.deltaCoeffs().boundaryField()[patchi];
}
const labelUList& faceCells = p.faceCells();
const labelUList& faceCells =
fGrad.boundaryField()[patchi].patch().faceCells();
// Build the d-vectors
vectorField pd = p.delta();
const Field<GradType> neighbourSecondfGrad
(
......
......@@ -130,20 +130,7 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
const labelUList& faceCells = p.patch().faceCells();
// Build the d-vectors
vectorField pd
(
mesh.Sf().boundaryField()[patchi]
/ (
mesh.magSf().boundaryField()[patchi]
* mesh.deltaCoeffs().boundaryField()[patchi]
)
);
if (!mesh.orthogonal())
{
pd -= mesh.correctionVectors().boundaryField()[patchi]
/mesh.deltaCoeffs().boundaryField()[patchi];
}
vectorField pd = p.delta();
if (p.coupled())
{
......@@ -196,21 +183,7 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
const labelUList& faceCells = p.faceCells();
// Build the d-vectors
vectorField pd
(
mesh.Sf().boundaryField()[patchi]
/(
mesh.magSf().boundaryField()[patchi]
*mesh.deltaCoeffs().boundaryField()[patchi]
)
);
if (!mesh.orthogonal())
{
pd -= mesh.correctionVectors().boundaryField()[patchi]
/mesh.deltaCoeffs().boundaryField()[patchi];
}
vectorField pd = p.delta();
if (p.coupled())
{
......
......@@ -211,7 +211,7 @@ public:
// Add the patch constructor functions to the hash tables
#define makeFvLaplacianTypeScheme(SS, Type, GType) \
#define makeFvLaplacianTypeScheme(SS, GType, Type) \
\
typedef SS<Type, GType> SS##Type##GType; \
defineNamedTemplateTypeNameAndDebug(SS##Type##GType, 0); \
......@@ -224,13 +224,20 @@ public:
#define makeFvLaplacianScheme(SS) \
\
makeFvLaplacianTypeScheme(SS, scalar, scalar) \
makeFvLaplacianTypeScheme(SS, scalar, symmTensor) \
makeFvLaplacianTypeScheme(SS, scalar, tensor) \
makeFvLaplacianTypeScheme(SS, vector, scalar) \
makeFvLaplacianTypeScheme(SS, sphericalTensor, scalar) \
makeFvLaplacianTypeScheme(SS, symmTensor, scalar) \
makeFvLaplacianTypeScheme(SS, tensor, scalar) \
makeFvLaplacianTypeScheme(SS, scalar, vector) \
makeFvLaplacianTypeScheme(SS, symmTensor, vector) \
makeFvLaplacianTypeScheme(SS, tensor, vector) \
makeFvLaplacianTypeScheme(SS, scalar, sphericalTensor) \
makeFvLaplacianTypeScheme(SS, symmTensor, sphericalTensor) \
makeFvLaplacianTypeScheme(SS, tensor, sphericalTensor) \
makeFvLaplacianTypeScheme(SS, scalar, symmTensor) \
makeFvLaplacianTypeScheme(SS, symmTensor, symmTensor) \
makeFvLaplacianTypeScheme(SS, tensor, symmTensor) \
makeFvLaplacianTypeScheme(SS, scalar, tensor) \
makeFvLaplacianTypeScheme(SS, symmTensor, tensor) \
makeFvLaplacianTypeScheme(SS, tensor, tensor)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -50,7 +50,7 @@ Foam::fv::correctedSnGrad<Type>::fullGradCorrection
// construct GeometricField<Type, fvsPatchField, surfaceMesh>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tssf =
mesh.correctionVectors()
mesh.nonOrthCorrectionVectors()
& linear<typename outerProduct<vector, Type>::type>(mesh).interpolate
(
gradScheme<Type>::New
......@@ -88,7 +88,7 @@ Foam::fv::correctedSnGrad<Type>::correction
IOobject::NO_WRITE
),
mesh,
vf.dimensions()*mesh.deltaCoeffs().dimensions()
vf.dimensions()*mesh.nonOrthDeltaCoeffs().dimensions()
)
);
GeometricField<Type, fvsPatchField, surfaceMesh>& ssf = tssf();
......@@ -98,7 +98,7 @@ Foam::fv::correctedSnGrad<Type>::correction
ssf.replace
(
cmpt,
mesh.correctionVectors()
mesh.nonOrthCorrectionVectors()
& linear
<
typename
......
......@@ -96,13 +96,13 @@ public:
const GeometricField<Type, fvPatchField, volMesh>&
) const
{
return this->mesh().deltaCoeffs();
return this->mesh().nonOrthDeltaCoeffs();
}
//- Return true if this scheme uses an explicit correction
virtual bool corrected() const
{
return !this->mesh().orthogonal();
return true;
}
//- Return the explicit correction to the correctedSnGrad
......
......@@ -120,13 +120,13 @@ public:
const GeometricField<Type, fvPatchField, volMesh>&
) const
{
return this->mesh().deltaCoeffs();
return this->mesh().nonOrthDeltaCoeffs();
}
//- Return true if this scheme uses an explicit correction
virtual bool corrected() const
{
return !this->mesh().orthogonal();
return true;
}
//- Return the explicit correction to the limitedSnGrad
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Description
Simple central-difference snGrad scheme without non-orthogonal correction.
\*---------------------------------------------------------------------------*/
#include "orthogonalSnGrad.H"
#include "volFields.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace fv
{
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Type>
orthogonalSnGrad<Type>::~orthogonalSnGrad()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
orthogonalSnGrad<Type>::correction
(
const GeometricField<Type, fvPatchField, volMesh>&
) const
{
notImplemented
(
"orthogonalSnGrad<Type>::correction"
"(const GeometricField<Type, fvPatchField, volMesh>&)"
);
return tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >(NULL);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::fv::orthogonalSnGrad
Description
Simple central-difference snGrad scheme without non-orthogonal correction.
SourceFiles
orthogonalSnGrad.C
\*---------------------------------------------------------------------------*/
#ifndef orthogonalSnGrad_H
#define orthogonalSnGrad_H
#include "snGradScheme.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace fv
{
/*---------------------------------------------------------------------------*\
Class orthogonalSnGrad Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class orthogonalSnGrad
:
public snGradScheme<Type>
{
// Private Member Functions
//- Disallow default bitwise assignment
void operator=(const orthogonalSnGrad&);
public:
//- Runtime type information
TypeName("uncorrected");
// Constructors
//- Construct from mesh
orthogonalSnGrad(const fvMesh& mesh)
:
snGradScheme<Type>(mesh)
{}
//- Construct from mesh and data stream
orthogonalSnGrad(const fvMesh& mesh, Istream&)
:
snGradScheme<Type>(mesh)
{}
//- Destructor
virtual ~orthogonalSnGrad();
// Member Functions
//- Return the interpolation weighting factors for the given field
virtual tmp<surfaceScalarField> deltaCoeffs
(
const GeometricField<Type, fvPatchField, volMesh>&
) const
{
return this->mesh().deltaCoeffs();
}
//- Return true if this scheme uses an explicit correction
virtual bool corrected() const
{
return false;
}
//- Return the explicit correction to the orthogonalSnGrad
// for the given field
virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
correction(const GeometricField<Type, fvPatchField, volMesh>&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "orthogonalSnGrad.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Description
Simple central-difference snGrad scheme without non-orthogonal correction.
\*---------------------------------------------------------------------------*/
#include "orthogonalSnGrad.H"
#include "fvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace fv
{
makeSnGradScheme(orthogonalSnGrad)
}
}
// ************************************************************************* //
......@@ -109,7 +109,7 @@ public:
const GeometricField<Type, fvPatchField, volMesh>&
) const
{
return this->mesh().deltaCoeffs();
return this->mesh().nonOrthDeltaCoeffs();
}
//- Return true if this scheme uses an explicit correction
......
......@@ -247,7 +247,7 @@ Foam::label Foam::quadraticFitSnGradData::calcFit
scalarList singVals(minSize_);
label nSVDzeros = 0;
const scalar deltaCoeff = mesh().deltaCoeffs()[faci];
const scalar deltaCoeff = deltaCoeffs()[faci];
bool goodFit = false;
for (int iIt = 0; iIt < 10 && !goodFit; iIt++)
......
......@@ -96,7 +96,7 @@ public:
const GeometricField<Type, fvPatchField, volMesh>&
) const
{
return this->mesh().deltaCoeffs();
return this->mesh().nonOrthDeltaCoeffs();
}
//- Return true if this scheme uses an explicit correction
......
......@@ -66,9 +66,6 @@ protected:
//- Make patch weighting factors
virtual void makeWeights(scalarField&) const = 0;
//- Make patch face - neighbour cell distances
virtual void makeDeltaCoeffs(scalarField&) const = 0;
public:
......
......@@ -57,25 +57,6 @@ void Foam::cyclicFvPatch::makeWeights(scalarField& w) const
}
// Make patch face - neighbour cell distances
void Foam::cyclicFvPatch::makeDeltaCoeffs(scalarField& dc) const
{
//const cyclicPolyPatch& nbrPatch = cyclicPolyPatch_.neighbPatch();
const cyclicFvPatch& nbrPatch = neighbFvPatch();
const scalarField deltas(nf() & fvPatch::delta());
const scalarField nbrDeltas(nbrPatch.nf() & nbrPatch.fvPatch::delta());
forAll(deltas, facei)