Commit 7544110b authored by henry's avatar henry
Browse files

Improved extendedStencil handling.

parent 4e2027c1
......@@ -123,9 +123,9 @@ public:
List<List<T> >& stencilFld
) const;
//- Given weights interpolate vol field
//- Calculate weighted sum of vol field
template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > interpolate
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > weightedSum
(
const GeometricField<Type, fvPatchField, volMesh>& fld,
const List<List<scalar> >& stencilWeights
......
......@@ -77,7 +77,7 @@ void Foam::extendedStencil::collectData
template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
Foam::extendedStencil::interpolate
Foam::extendedStencil::weightedSum
(
const GeometricField<Type, fvPatchField, volMesh>& fld,
const List<List<scalar> >& stencilWeights
......@@ -120,7 +120,34 @@ Foam::extendedStencil::interpolate
sf[faceI] += stField[i]*stWeight[i];
}
}
// And what for boundaries?
// Coupled boundaries
/*
typename GeometricField<Type, fvsPatchField, surfaceMesh>::
GeometricBoundaryField& bSfCorr = sf.boundaryField();
forAll(bSfCorr, patchi)
{
fvsPatchField<Type>& pSfCorr = bSfCorr[patchi];
if (pSfCorr.coupled())
{
label faceI = pSfCorr.patch().patch().start();
forAll(pSfCorr, i)
{
const List<Type>& stField = stencilFld[faceI];
const List<scalar>& stWeight = stencilWeights[faceI];
forAll(stField, j)
{
pSfCorr[i] += stField[j]*stWeight[j];
}
faceI++;
}
}
}
*/
return tsfCorr;
}
......
......@@ -122,7 +122,7 @@ public:
const extendedStencil& stencil = cfd.stencil();
const List<scalarList>& f = cfd.fit();
return stencil.interpolate(vf, f);
return stencil.weightedSum(vf, f);
}
};
......
......@@ -41,6 +41,8 @@ namespace Foam
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
static int count = 0;
Foam::quadraticFitData::quadraticFitData
(
const fvMesh& mesh,
......@@ -114,13 +116,15 @@ Foam::quadraticFitData::quadraticFitData
interpPolySize[faci] = calcFit(stencilPoints[faci], faci);
}
interpPolySize.write();
Pout<< "count = " << count << endl;
if (debug)
{
Info<< "quadraticFitData::quadraticFitData() :"
<< "Finished constructing polynomialFit data"
<< endl;
interpPolySize.write();
}
}
......@@ -270,8 +274,8 @@ Foam::label Foam::quadraticFitData::calcFit
//goodFit = (fit0 > 0 && fit1 > 0);
goodFit =
(mag(fit0 - w[faci])/w[faci] < 0.5)
&& (mag(fit1 - (1 - w[faci]))/(1 - w[faci]) < 0.5);
(mag(fit0 - w[faci])/w[faci] < 0.15)
&& (mag(fit1 - (1 - w[faci]))/(1 - w[faci]) < 0.15);
//scalar w0Err = fit0/w[faci];
//scalar w1Err = fit1/(1 - w[faci]);
......@@ -312,11 +316,21 @@ Foam::label Foam::quadraticFitData::calcFit
}
}
//static const scalar alpha = 1.5;
//static const scalar beta = alpha/0.5;
// 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)
{
if ((mag(fit_[faci][0] - w[faci])/w[faci] < 0.15)
&& (mag(fit_[faci][1] - (1 - w[faci]))/(1 - w[faci]) < 0.15))
{
count++;
//Pout<< "fit " << mag(fit_[faci][0] - w[faci])/w[faci] << " " << mag(fit_[faci][1] - (1 - w[faci]))/(1 - w[faci]) << endl;
}
// scalar limiter =
// max
// (
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment