Skip to content
Snippets Groups Projects
Commit 79efd0a5 authored by Mark Olesen's avatar Mark Olesen
Browse files

BSpline code cleanup

parent 51491420
Branches
Tags
No related merge requests found
...@@ -38,75 +38,70 @@ Foam::pointField Foam::BSpline::findKnots ...@@ -38,75 +38,70 @@ Foam::pointField Foam::BSpline::findKnots
const vector& sndend const vector& sndend
) )
{ {
label newnKnots(allknots.size() + 2); const label NKnots = allknots.size();
label NKnots(allknots.size());
pointField newknots(newnKnots);
// set up 1/6 and 2/3 which are the matrix elements throughout most // set up 1/6 and 2/3 which are the matrix elements throughout most
// of the matrix // of the matrix
register scalar oneSixth = 1.0/6.0; register const scalar oneSixth = 1.0/6.0;
register scalar twoThird = 2.0/3.0; register const scalar twoThird = 2.0/3.0;
simpleMatrix<vector> M(newnKnots, pTraits<scalar>::zero);
simpleMatrix<vector> M(NKnots+2, pTraits<scalar>::zero);
// set up the matrix // set up the matrix
M[0][0] = -0.5*scalar(NKnots - 1); M[0][0] = -0.5*scalar(NKnots - 1);
M[0][2] = 0.5*scalar(NKnots - 1); M[0][2] = 0.5*scalar(NKnots - 1);
for (register label i=1; i<newnKnots-1; i++) for (register label i = 1; i <= NKnots; i++)
{ {
M[i][i-1] = oneSixth; M[i][i-1] = oneSixth;
M[i][i] = twoThird; M[i][i] = twoThird;
M[i][i+1] = oneSixth; M[i][i+1] = oneSixth;
} }
M[newnKnots - 1][newnKnots - 3] = -0.5*scalar(NKnots - 1); M[NKnots+1][NKnots-1] = -0.5*scalar(NKnots - 1);
M[newnKnots - 1][newnKnots - 1] = 0.5*scalar(NKnots - 1); M[NKnots+1][NKnots+1] = 0.5*scalar(NKnots - 1);
// set up the vector // set up the vector
for (register label i = 1; i <= NKnots; i++)
for (label i=1; i<=NKnots; i++)
{ {
M.source()[i] = allknots[i-1]; M.source()[i] = allknots[i-1];
} }
// set the gradients at the two ends
// set the gradients at the two ends
if (mag(fstend)<1e-8) if (mag(fstend)<1e-8)
{ {
// set to the default : forward differences on the end knots // default : forward differences on the end knots
M.source()[0] = allknots[1] - allknots[0]; M.source()[0] = allknots[1] - allknots[0];
M.source()[0] /= mag(M.source()[0]); M.source()[0] /= mag(M.source()[0]);
}
else
{
// set to the gradient vector provided
M.source()[0] = fstend/mag(fstend);
}
if (mag(sndend)<1e-8)
{
// default : forward differences on the end knots
M.source()[NKnots+1] = M.source()[NKnots-1] - M.source()[NKnots]; M.source()[NKnots+1] = M.source()[NKnots-1] - M.source()[NKnots];
M.source()[NKnots+1] /= mag(M.source()[NKnots+1]); M.source()[NKnots+1] /= mag(M.source()[NKnots+1]);
} }
else else
{ {
// set to the gradient vectors provided // set to the gradient vector provided
M.source()[0] = fstend/mag(fstend);
M.source()[NKnots+1] = sndend/mag(sndend); M.source()[NKnots+1] = sndend/mag(sndend);
} }
// invert the equation to find the control knots
newknots = M.solve(); // invert the equation to find the control knots
return M.solve();
return newknots;
} }
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::BSpline::BSpline(const pointField& Knots)
:
spline(findKnots(Knots))
{}
Foam::BSpline::BSpline Foam::BSpline::BSpline
( (
const pointField& Knots, const pointField& Knots,
......
...@@ -56,8 +56,8 @@ class BSpline ...@@ -56,8 +56,8 @@ class BSpline
pointField findKnots pointField findKnots
( (
const pointField&, const pointField&,
const vector& fstend = vector::zero, const vector& fstend,
const vector& sndend = vector::zero const vector& sndend
); );
//- Disallow default bitwise copy construct //- Disallow default bitwise copy construct
...@@ -71,15 +71,12 @@ public: ...@@ -71,15 +71,12 @@ public:
// Constructors // Constructors
//- Construct from components
BSpline(const pointField& knots);
//- Construct from components //- Construct from components
BSpline BSpline
( (
const pointField& knots, const pointField& knots,
const vector& fstend, const vector& fstend = vector::zero,
const vector& sndend const vector& sndend = vector::zero
); );
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment