Go to the documentation of this file.
35 template<
class Form,
class ExtendedStencil,
class Polynomial>
39 const ExtendedStencil& stencil,
40 const bool linearCorrection,
41 const scalar linearLimitFactor,
42 const scalar centralWeight
47 linearCorrection_(linearCorrection),
48 linearLimitFactor_(linearLimitFactor),
49 centralWeight_(centralWeight),
50 #ifdef SPHERICAL_GEOMETRY
53 dim_(
mesh.nGeometricD()),
55 minSize_(Polynomial::nTerms(dim_))
58 if (linearLimitFactor <= SMALL || linearLimitFactor > 3)
61 <<
"linearLimitFactor requested = " << linearLimitFactor
62 <<
" should be between zero and 3"
70 template<
class FitDataType,
class ExtendedStencil,
class Polynomial>
83 #ifndef SPHERICAL_GEOMETRY
84 if (
mesh.nGeometricD() <= 2)
86 if (
mesh.geometricD()[0] == -1)
90 else if (
mesh.geometricD()[1] == -1)
102 kdir =
mesh.points()[
f[0]] -
mesh.faceCentres()[facei];
106 kdir =
mesh.faceCentres()[facei];
109 if (
mesh.nGeometricD() == 3)
112 kdir -= (idir & kdir)*idir;
114 scalar magk =
mag(kdir);
131 template<
class FitDataType,
class ExtendedStencil,
class Polynomial>
143 findFaceDirs(idir, jdir, kdir, facei);
147 wts[0] = centralWeight_;
148 if (linearCorrection_)
150 wts[1] = centralWeight_;
172 d.
x() = (
p -
p0)&idir;
173 d.
y() = (
p -
p0)&jdir;
174 #ifndef SPHERICAL_GEOMETRY
175 d.
z() = (
p -
p0)&kdir;
188 Polynomial::addCoeffs
198 for (label i = 0; i <
B.m(); i++)
205 label stencilSize =
C.size();
208 bool goodFit =
false;
209 for (
int iIt = 0; iIt < 8 && !goodFit; iIt++)
217 for (label i=0; i<stencilSize; i++)
219 coeffsi[i] = wts[0]*wts[i]*invB(0, i);
220 if (
mag(coeffsi[i]) > maxCoeff)
222 maxCoeff =
mag(coeffsi[i]);
227 if (linearCorrection_)
230 (
mag(coeffsi[0] - wLin) < linearLimitFactor_*wLin)
231 && (
mag(coeffsi[1] - (1 - wLin)) < linearLimitFactor_*(1 - wLin))
238 (
mag(coeffsi[0] - 1.0) < linearLimitFactor_*1.0)
267 if (linearCorrection_)
272 for (label j = 0; j <
B.n(); j++)
278 for (label i = 0; i <
B.m(); i++)
288 if (linearCorrection_)
292 coeffsi[1] -= 1 - wLin;
305 <<
"Could not fit face " << facei
306 <<
" Weights = " << coeffsi
307 <<
", reverting to linear." <<
nl
308 <<
" Linear weights " << wLin <<
" " << 1 - wLin <<
endl;
316 template<
class FitDataType,
class ExtendedStencil,
class Polynomial>
const Cmpt & x() const
Access to the vector x component.
static constexpr const zero Zero
Global zero (0)
static const Foam::dimensionedScalar B("", Foam::dimless, 18.678)
Ostream & endl(Ostream &os)
Add newline and flush stream.
void findFaceDirs(vector &idir, vector &jdir, vector &kdir, const label faci)
Find the normal direction (i) and j and k directions for face faci.
const Cmpt & z() const
Access to the vector z component.
#define forAll(list, i)
Loop across all elements in list.
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
scalarRectangularMatrix VSinvUt() const
Return the matrix product V S^(-1) U^T (the pseudo inverse)
bool movePoints()
Recalculate weights (but not stencil) when the mesh moves.
Mesh data needed to do the Finite Volume discretisation.
Vector< scalar > vector
A scalar version of the templated Vector.
errorManipArg< error, int > exit(error &err, const int errNo=1)
VectorSpace< Form, Cmpt, Ncmpts > normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
FitData(const fvMesh &mesh, const ExtendedStencil &stencil, const bool linearCorrection, const scalar linearLimitFactor, const scalar centralWeight)
Construct from components.
const Cmpt & y() const
Access to the vector y component.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
virtual void calcFit()=0
Calculate the fit for all the faces.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Singular value decomposition of a rectangular matrix.
Templated abstract base-class for optional mesh objects used to automate their allocation to the mesh...
A face is a list of labels corresponding to mesh vertices.
const volScalarField & p0
Graphite solid properties.
void setSize(const label newSize)
Alias for resize(const label)
#define WarningInFunction
Report a warning using Foam::Warning.