Commit 5648be03 authored by Mark Olesen's avatar Mark Olesen
Browse files

Added Catmull-Rom splines to blockMesh.

- the blockMesh interface is splineEdge.H, selectable as "spline"

The first tests look fine - it works as expected for the case with
buggy polySpline reported on the forum. Should of course do some more
extensive testing.

The advantages compared to the current B-Spline implementation:

- Doesn't need a matrix solver.
- The coding resembles something that can be found in the literature.
- In contrast to the B-Spline implementation, it is fairly clear what
  is actually going on. I don't even know if the B-Spline are actually
  B-Spline, Beta-Splines or something else.
- Catmull-Rom splines seem to be what all the graphics people have as
  their stable workhorse.

We now have 3 different names for splines in blockMesh:
- "spline" - *new* Catmull-Rom for arbitrary segments.
- "simpleSpline" - B-Spline for a single segment
- "polySpline" - B-Spline for a multiple segments

Assuming the Catmull-Rom splines continue to behave nicely, there is
no reason to keep the other (broken) B-Splines. This would help clean
up the blockMesh interface too.

Placed the older ones under legacy/ for easier identification in the
future.

TODO:
- currently no handling of non-zero end tangents
- could be extended to handle closed loops, which might be useful
  for feature edges from CAD (eg, for the cvm mesher)
parent 9b92e3c4
curvedEdges/CatmullRomSpline.C
curvedEdges/polyLine.C
curvedEdges/arcEdge.C
curvedEdges/curvedEdge.C
curvedEdges/lineEdge.C
curvedEdges/polyLine.C
curvedEdges/polyLineEdge.C
curvedEdges/arcEdge.C
curvedEdges/spline.C
curvedEdges/BSpline.C
curvedEdges/simpleSplineEdge.C
curvedEdges/polySplineEdge.C
curvedEdges/lineDivide.C
curvedEdges/splineEdge.C
curvedEdges/legacy/spline.C
curvedEdges/legacy/BSpline.C
curvedEdges/legacy/simpleSplineEdge.C
curvedEdges/legacy/polySplineEdge.C
blockDescriptor/blockDescriptor.C
blockDescriptor/blockDescriptorEdges.C
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "CatmullRomSpline.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::CatmullRomSpline::CatmullRomSpline
(
const pointField& Knots,
const vector&,
const vector&
)
:
polyLine(Knots)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::point Foam::CatmullRomSpline::position(const scalar mu) const
{
// endpoints
if (mu < SMALL)
{
return points()[0];
}
else if (mu > 1 - SMALL)
{
return points()[points().size()-1];
}
scalar lambda = mu;
label segment = localParameter(lambda);
return position(segment, lambda);
}
Foam::point Foam::CatmullRomSpline::position
(
const label segment,
const scalar mu
) const
{
const point& p0 = points()[segment];
const point& p1 = points()[segment+1];
// special cases - no calculation needed
if (segment < 0 || mu < 0.0)
{
return p0;
}
else if (segment > nSegments() || mu >= 1.0)
{
return p1;
}
// determine the end points
point e0;
point e1;
if (segment == 0)
{
// end: simple reflection
e0 = 2.0 * p0 - p1;
}
else
{
e0 = points()[segment-1];
}
if (segment+1 == nSegments())
{
// end: simple reflection
e1 = 2.0 * p1 - p0;
}
else
{
e1 = points()[segment+2];
}
return 0.5 *
(
( 2 * p0 )
+ mu *
(
( -e0 + p1 )
+ mu *
(
( 2*e0 - 5*p0 + 4*p1 - e1 )
+ mu *
( -e0 + 3*p0 - 3*p1 + e1 )
)
)
);
}
Foam::scalar Foam::CatmullRomSpline::length() const
{
notImplemented("CatmullRomSpline::length() const");
return 1.0;
}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2009-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::CatmullRomSpline
Description
An implementation of Catmull-Rom splines (sometime as known as
Overhauser splines).
In this implementation, the end tangents are created
automatically by reflection.
In matrix form, the @e local interpolation on the interval t=[0..1] is
described as follows:
@verbatim
P(t) = 0.5 * [ t^3 t^2 t 1 ] * [ -1 3 -3 1 ] * [ P-1 ]
[ 2 -5 4 -1 ] [ P0 ]
[ -1 0 1 0 ] [ P1 ]
[ 0 2 0 0 ] [ P2 ]
@endverbatim
Where P-1 and P2 represent the neighbouring points or the
extrapolated end points. Simple reflection is used to
automatically create the end points.
The spline is discretized based on the chord length of the
individual segments. In rare cases (sections with very high
curvatures), the resulting distribution may be sub-optimal.
SeeAlso
http://www.algorithmist.net/catmullrom.html provides a nice
introduction
ToDo
A future implementation could also handle closed splines - either
when the start/end points are identically or when specified.
SourceFiles
CatmullRomSpline.C
\*---------------------------------------------------------------------------*/
#ifndef CatmullRomSpline_H
#define CatmullRomSpline_H
#include "polyLine.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class CatmullRomSpline Declaration
\*---------------------------------------------------------------------------*/
class CatmullRomSpline
:
public polyLine
{
// Private Member Functions
//- Disallow default bitwise copy construct
CatmullRomSpline(const CatmullRomSpline&);
//- Disallow default bitwise assignment
void operator=(const CatmullRomSpline&);
public:
// Constructors
//- Construct from components
CatmullRomSpline
(
const pointField& knots,
const vector& begTangentNotImplemented = vector::zero,
const vector& endTangentNotImplemented = vector::zero
);
// Member Functions
//- Return the point position corresponding to the curve parameter
// 0 <= lambda <= 1
point position(const scalar lambda) const;
//- Return the point position corresponding to the local parameter
// 0 <= lambda <= 1 on the given segment
point position(const label segment, const scalar lambda) const;
//- Return the length of the curve
scalar length() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
......@@ -71,8 +71,7 @@ Foam::cylindricalCS Foam::arcEdge::calcAngle()
vector r3(p3_ - centre);
// find angles
scalar tmp = (r3&r1)/(mag(r3)*mag(r1));
angle_ = radToDeg(acos(tmp));
angle_ = radToDeg(acos((r3 & r1)/(mag(r3) * mag(r1))));
// check if the vectors define an exterior or an interior arcEdge
if (((r1 ^ r2) & (r1 ^ r3)) < 0.0)
......@@ -99,7 +98,7 @@ Foam::cylindricalCS Foam::arcEdge::calcAngle()
radius_ = mag(r3);
// set up and return the local coordinate system
return cylindricalCS("tmpCS", centre, tempAxis, r1);
return cylindricalCS("arcEdgeCS", centre, tempAxis, r1);
}
......
......@@ -55,14 +55,15 @@ class arcEdge
// Private data
point p1_, p2_, p3_;
cylindricalCS cs_;
scalar angle_;
scalar radius_;
cylindricalCS cs_;
cylindricalCS calcAngle();
// Private Member Functions
//- Calculate the coordinate system, angle and radius
cylindricalCS calcAngle();
//- Disallow default bitwise copy construct
arcEdge(const arcEdge&);
......
......@@ -108,7 +108,7 @@ Foam::autoPtr<Foam::curvedEdge> Foam::curvedEdge::New
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::pointField Foam::curvedEdge::fullKnotList
Foam::pointField Foam::curvedEdge::appendEndPoints
(
const pointField& points,
const label start,
......
......@@ -64,9 +64,9 @@ protected:
// Protected Member Functions
//- Return a complete knotList by adding the start/end points
//- Return a complete point field by appending the start/end points
// to the given list
static pointField fullKnotList
static pointField appendEndPoints
(
const pointField&,
const label start,
......
......@@ -56,7 +56,7 @@ Foam::pointField Foam::polySplineEdge::intervening
{
BSpline spl
(
fullKnotList(points_, start_, end_, otherknots),
appendEndPoints(curvedEdge::points_, start_, end_, otherknots),
fstend,
sndend
);
......@@ -73,7 +73,7 @@ Foam::pointField Foam::polySplineEdge::intervening
interval /= nBetweenKnots + 1;
pointField ans(nSize);
ans[0] = points_[start_];
ans[0] = curvedEdge::points_[start_];
register scalar index(init);
for (register label i=1; i<nSize-1; i++)
......@@ -82,7 +82,7 @@ Foam::pointField Foam::polySplineEdge::intervening
ans[i] = spl.realPosition(index);
}
ans[nSize-1] = points_[end_];
ans[nSize-1] = curvedEdge::points_[end_];
return ans;
}
......@@ -128,8 +128,8 @@ Foam::polySplineEdge::polySplineEdge
vector fstend(is);
vector sndend(is);
controlPoints_ = intervening(otherKnots_, nInterKnots, fstend, sndend);
calcDistances();
polyLine::points_ = intervening(otherKnots_, nInterKnots, fstend, sndend);
calcParam();
}
......
......@@ -26,7 +26,7 @@ Class
Foam::polySplineEdge
Description
A spline representation via a polyLine
A curvedEdge interface for B-splines.
SourceFiles
polySplineEdge.C
......
......@@ -48,7 +48,7 @@ Foam::simpleSplineEdge::simpleSplineEdge
)
:
curvedEdge(points, start, end),
BSpline(fullKnotList(points, start, end, otherknots))
BSpline(appendEndPoints(points, start, end, otherknots))
{}
......@@ -63,14 +63,14 @@ Foam::simpleSplineEdge::simpleSplineEdge
)
:
curvedEdge(points, start, end),
BSpline(fullKnotList(points, start, end, otherknots), fstend, sndend)
BSpline(appendEndPoints(points, start, end, otherknots), fstend, sndend)
{}
Foam::simpleSplineEdge::simpleSplineEdge(const pointField& points, Istream& is)
:
curvedEdge(points, is),
BSpline(fullKnotList(points, start_, end_, pointField(is)))
BSpline(appendEndPoints(points, start_, end_, pointField(is)))
{}
......
......@@ -26,7 +26,7 @@ Class
Foam::simpleSplineEdge
Description
The actual access class for Bspline
A curvedEdge interface for B-splines.
SourceFiles
simpleSplineEdge.C
......
......@@ -55,6 +55,10 @@ Foam::lineEdge::lineEdge(const pointField& points, Istream& is)
curvedEdge(points, is)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * //
Foam::lineEdge::~lineEdge()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
......
......@@ -26,7 +26,7 @@ Class
Foam::lineEdge
Description
Defines a straight line between the start point and the end point.
A straight edge between the start point and the end point.
SourceFiles
lineEdge.C
......@@ -70,13 +70,12 @@ public:
//- Construct from components
lineEdge(const pointField&, const label start, const label end);
//- Construct from Istream setting pointsList
//- Construct from Istream with a pointField
lineEdge(const pointField&, Istream&);
//- Destructor
virtual ~lineEdge()
{}
//- Destructor
virtual ~lineEdge();
// Member Functions
......
......@@ -29,30 +29,27 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// calcDistances generates the distances_ lookup table (cumulative
// distance along the line) from the individual vectors to the points
void Foam::polyLine::calcDistances()
void Foam::polyLine::calcParam()
{
distances_.setSize(controlPoints_.size());
param_.setSize(points_.size());
if (distances_.size())
if (param_.size())
{
distances_[0] = 0.0;
param_[0] = 0.0;
for (label i=1; i < distances_.size(); i++)
for (label i=1; i < param_.size(); i++)
{
distances_[i] = distances_[i-1] +
mag(controlPoints_[i] - controlPoints_[i-1]);
param_[i] = param_[i-1] +
mag(points_[i] - points_[i-1]);
}
// normalize on the interval 0-1
lineLength_ = distances_[distances_.size()-1];
for (label i=1; i < distances_.size() - 1; i++)
lineLength_ = param_[param_.size()-1];
for (label i=1; i < param_.size() - 1; i++)
{
distances_[i] /= lineLength_;
param_[i] /= lineLength_;
}
distances_[distances_.size()-1] = 1.0;
param_[param_.size()-1] = 1.0;
}
else
{
......@@ -66,19 +63,69 @@ void Foam::polyLine::calcDistances()
Foam::polyLine::polyLine(const pointField& ps)
:
controlPoints_(ps),
distances_(0),
lineLength_(0.0)
points_(ps),
lineLength_(0.0),
param_(0)
{
calcDistances();
calcParam();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::pointField& Foam::polyLine::controlPoints() const
const Foam::pointField& Foam::polyLine::points() const
{
return points_;
}
Foam::label Foam::polyLine::nSegments() const
{
return points_.size()-1;
}
Foam::label Foam::polyLine::localParameter(scalar& lambda) const
{
return controlPoints_;
// check range of lambda
if (lambda < 0 || lambda > 1)
{
FatalErrorIn("polyLine::localParameter(scalar&)")
<< "Parameter out-of-range, "
<< "lambda = " << lambda
<< abort(FatalError);
}
// check endpoints
if (lambda < SMALL)
{
lambda = 0;
return 0;
}
else if (lambda > 1 - SMALL)
{
lambda = 1;
return nSegments();
}
// search table of cumulative distances to find which line-segment
// we are on. Check the upper bound.
label segmentI = 1;
while (param_[segmentI] < lambda)
{
segmentI++;
}
segmentI--; // we want the corresponding lower bound
// the local parameter [0-1] on this line segment
lambda =
(
( lambda - param_[segmentI] )
/ ( param_[segmentI+1] - param_[segmentI] )
);
return segmentI;
}
......@@ -96,31 +143,32 @@ Foam::point Foam::polyLine::position(const scalar lambda) const
// check endpoints
if (lambda < SMALL)
{
return controlPoints_[0];
return points_[0];
}
else if (lambda > 1 - SMALL)
{
return controlPoints_[controlPoints_.size()-1];
return points_[points_.size()-1];
}
// search table of cumulative distances to find which line-segment
// we are on. Check the upper bound.
label i = 1;
while (distances_[i] < lambda)
label segmentI = 1;