Commit 73f9f7f7 authored by Mark Olesen's avatar Mark Olesen
Browse files

Remove legacy splines code and use CatmullRomSpline as 'spline'

- compatibility:
  * 'polySpline' and 'simpleSpline' accepted
  * detect and discard end tangent specifications

- a BSpline is also included (eg, for future support of NURBS), but is
  not selectable from blockMesh since it really isn't as nice as the
  Catmull-Rom (ie, doesn't *exactly* go through each point).
parent 1b0cf102
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2009-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -29,11 +29,17 @@ License
#include "IFstream.H"
#include "BSpline.H"
#include "BSpline2.H"
#include "CatmullRomSpline.H"
using namespace Foam;
inline Ostream& printPoint(Ostream& os, const point& p)
{
os << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
return os;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
......@@ -41,8 +47,7 @@ int main(int argc, char *argv[])
{
argList::noParallel();
argList::validArgs.insert("file .. fileN");
argList::addBoolOption("Bold", "B-Spline implementation OLD");
argList::addBoolOption("Bnew", "B-Spline implementation NEW");
argList::addBoolOption("B", "B-Spline implementation");
argList::addBoolOption("CMR", "catmull-rom spline (default)");
argList::addOption
(
......@@ -58,12 +63,11 @@ int main(int argc, char *argv[])
args.printUsage();
}
bool useBSplineOld = args.optionFound("Bold");
bool useBSplineNew = args.optionFound("Bnew");
bool useBSpline = args.optionFound("B");
bool useCatmullRom = args.optionFound("CMR");
label nSeg = args.optionLookupOrDefault<label>("n", 20);
if (!useCatmullRom && !useBSplineOld && !useBSplineNew)
if (!useCatmullRom && !useBSpline)
{
Info<<"defaulting to Catmull-Rom spline" << endl;
useCatmullRom = true;
......@@ -77,55 +81,50 @@ int main(int argc, char *argv[])
List<pointField> pointFields(ifs);
forAll(pointFields, splineI)
{
Info<<"\noriginal points: " << pointFields[splineI] << nl;
Info<<"\n# points:" << endl;
forAll(pointFields[splineI], ptI)
{
printPoint(Info, pointFields[splineI][ptI]);
}
if (useBSplineOld)
if (useBSpline)
{
BSpline spl(pointFields[splineI]);
Info<< nl
<< "B-Spline interpolation: OLD" << nl
<< "----------------------" << endl;
Info<< nl << "# B-Spline" << endl;
for (label segI = 0; segI <= nSeg; ++segI)
{
scalar lambda = scalar(segI)/scalar(nSeg);
Info<< spl.position(lambda) << " // " << lambda << endl;
printPoint(Info, spl.position(lambda));
}
}
if (useBSplineNew)
if (useCatmullRom)
{
BSpline2 spl(pointFields[splineI]);
CatmullRomSpline spl(pointFields[splineI]);
Info<< nl
<< "B-Spline interpolation: NEW" << nl
<< "----------------------" << endl;
Info<< nl <<"# Catmull-Rom" << endl;
for (label segI = 0; segI <= nSeg; ++segI)
{
scalar lambda = scalar(segI)/scalar(nSeg);
Info<< spl.position(lambda) << " // " << lambda << endl;
printPoint(Info, spl.position(lambda));
}
}
if (useCatmullRom)
{
CatmullRomSpline spl
(
pointFields[splineI]
);
polyLine pl(pointFields[splineI]);
Info<< nl
<<"Catmull-Rom interpolation:" << nl
<< "-------------------------" << endl;
Info<< nl <<"# polyList" << endl;
for (label segI = 0; segI <= nSeg; ++segI)
{
scalar lambda = scalar(segI)/scalar(nSeg);
Info<< spl.position(lambda) << " // " << lambda << endl;
printPoint(Info, pl.position(lambda));
}
}
}
......
......@@ -6,64 +6,41 @@
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
(
// Upper body longitudinal splines.
// Upper body longitudinal splines
(
(-0.22685 -0.01125166 0) // 7
(-0.21685 -0.01340204 0)
(-0.20685 -0.01529684 0)
(-0.19685 -0.01694748 0)
(-0.18685 -0.01836538 0)
(-0.17685 -0.01956197 0)
(-0.16685 -0.02054868 0)
(-0.15685 -0.02133693 0)
(-0.14685 -0.02193816 0)
(-0.13685 -0.02236377 0)
(-0.12685 -0.02262521 0)
(-0.11685 -0.02273389 0) // 2
(0.22685 0.01125166 0) // 7
(0.21685 0.01340204 0)
(0.20685 0.01529684 0)
(0.19685 0.01694748 0)
(0.18685 0.01836538 0)
(0.17685 0.01956197 0)
(0.16685 0.02054868 0)
(0.15685 0.02133693 0)
(0.14685 0.02193816 0)
(0.13685 0.02236377 0)
(0.12685 0.02262521 0)
(0.11685 0.02273389 0) // 2
)
// sine function
(
(-0.22685 0 0.01125166) // 8
(-0.21685 0 0.01340204)
(-0.20685 0 0.01529684)
(-0.19685 0 0.01694748)
(-0.18685 0 0.01836538)
(-0.17685 0 0.01956197)
(-0.16685 0 0.02054868)
(-0.15685 0 0.02133693)
(-0.14685 0 0.02193816)
(-0.13685 0 0.02236377)
(-0.12685 0 0.02262521)
(-0.11685 0 0.02273389) // 3
(0 0 0)
(45 0.70707 0)
(90 1 0)
(135 0.70707 0)
(180 0 0)
(225 -0.70707 0)
(270 -1 0)
(315 -0.70707 0)
(360 0 0)
)
// cosine function, but with extremely few points
(
(-0.22685 0.01125166 0) // 9
(-0.21685 0.01340204 0)
(-0.20685 0.01529684 0)
(-0.19685 0.01694748 0)
(-0.18685 0.01836538 0)
(-0.17685 0.01956197 0)
(-0.16685 0.02054868 0)
(-0.15685 0.02133693 0)
(-0.14685 0.02193816 0)
(-0.13685 0.02236377 0)
(-0.12685 0.02262521 0)
(-0.11685 0.02273389 0) // 4
(0 1 0)
(180 -1 0)
(360 1 0)
)
(
(-0.22685 0 -0.01125166) // 6
(-0.21685 0 -0.01340204)
(-0.20685 0 -0.01529684)
(-0.19685 0 -0.01694748)
(-0.18685 0 -0.01836538)
(-0.17685 0 -0.01956197)
(-0.16685 0 -0.02054868)
(-0.15685 0 -0.02133693)
(-0.14685 0 -0.02193816)
(-0.13685 0 -0.02236377)
(-0.12685 0 -0.02262521)
(-0.11685 0 -0.02273389) // 1
)
);
curvedEdges/BSpline2.C
curvedEdges/BSpline.C
curvedEdges/CatmullRomSpline.C
curvedEdges/polyLine.C
......@@ -9,11 +9,6 @@ curvedEdges/polyLineEdge.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
......
......@@ -25,24 +25,23 @@ License
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "BSpline2.H"
#include "BSpline.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::BSpline2::BSpline2
Foam::BSpline::BSpline
(
const pointField& Knots,
const vector&,
const vector&
const pointField& knots,
const bool closed
)
:
polyLine(Knots)
polyLine(knots, closed)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::point Foam::BSpline2::position(const scalar mu) const
Foam::point Foam::BSpline::position(const scalar mu) const
{
// endpoints
if (mu < SMALL)
......@@ -60,25 +59,36 @@ Foam::point Foam::BSpline2::position(const scalar mu) const
}
Foam::point Foam::BSpline2::position
Foam::point Foam::BSpline::position
(
const label segment,
const scalar mu
) const
{
// out-of-bounds
if (segment < 0)
{
return points().first();
}
else if (segment > nSegments())
{
return points().last();
}
const point& p0 = points()[segment];
const point& p1 = points()[segment+1];
// special cases - no calculation needed
if (segment < 0 || mu < 0.0)
if (mu <= 0.0)
{
return p0;
}
else if (segment > nSegments() || mu >= 1.0)
else if (mu >= 1.0)
{
return p1;
}
// determine the end points
point e0;
point e1;
......@@ -86,7 +96,7 @@ Foam::point Foam::BSpline2::position
if (segment == 0)
{
// end: simple reflection
e0 = 2.0 * p0 - p1;
e0 = 2*p0 - p1;
}
else
{
......@@ -96,7 +106,7 @@ Foam::point Foam::BSpline2::position
if (segment+1 == nSegments())
{
// end: simple reflection
e1 = 2.0 * p1 - p0;
e1 = 2*p1 - p0;
}
else
{
......@@ -121,9 +131,9 @@ Foam::point Foam::BSpline2::position
}
Foam::scalar Foam::BSpline2::length() const
Foam::scalar Foam::BSpline::length() const
{
notImplemented("BSpline2::length() const");
notImplemented("BSpline::length() const");
return 1.0;
}
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2009-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -26,7 +26,33 @@ Class
Foam::BSpline
Description
A cubic spline going through all the knots
An implementation of B-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) = 1/6 * [ t^3 t^2 t 1 ] * [ -1 3 -3 1 ] * [ P-1 ]
[ 3 -6 3 0 ] [ P0 ]
[ -3 0 3 0 ] [ P1 ]
[ 1 4 1 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
CatmullRomSpline
ToDo
A future implementation could also handle closed splines.
SourceFiles
BSpline.C
......@@ -36,7 +62,7 @@ SourceFiles
#ifndef BSpline_H
#define BSpline_H
#include "spline.H"
#include "polyLine.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -49,17 +75,10 @@ namespace Foam
class BSpline
:
public spline
public polyLine
{
// Private Member Functions
pointField findKnots
(
const pointField&,
const vector& fstend,
const vector& sndend
);
//- Disallow default bitwise copy construct
BSpline(const BSpline&);
......@@ -75,21 +94,20 @@ public:
BSpline
(
const pointField& knots,
const vector& fstend = vector::zero,
const vector& sndend = vector::zero
const bool notImplementedClosed = false
);
// Member Functions
//- Return the real point position corresponding to the curve parameter
// 0 <= lambda <= 1
point realPosition(const scalar lambda) const;
//- 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;
};
......
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::BSpline2
Description
An implementation of B-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) = 1/6 * [ t^3 t^2 t 1 ] * [ -1 3 -3 1 ] * [ P-1 ]
[ 3 -6 3 0 ] [ P0 ]
[ -3 0 3 0 ] [ P1 ]
[ 1 4 1 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
CatmullRomSpline
ToDo
A future implementation could also handle closed splines - either
when the start/end points are identically or when specified.
SourceFiles
BSpline2.C
\*---------------------------------------------------------------------------*/
#ifndef BSpline2_H
#define BSpline2_H
#include "polyLine.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class BSpline2 Declaration
\*---------------------------------------------------------------------------*/
class BSpline2
:
public polyLine
{
// Private Member Functions
//- Disallow default bitwise copy construct
BSpline2(const BSpline2&);
//- Disallow default bitwise assignment
void operator=(const BSpline2&);
public:
// Constructors
//- Construct from components
BSpline2
(
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
// ************************************************************************* //
......@@ -31,12 +31,11 @@ License
Foam::CatmullRomSpline::CatmullRomSpline
(
const pointField& Knots,
const vector&,
const vector&
const pointField& knots,
const bool closed
)
:
polyLine(Knots)
polyLine(knots, closed)
{}
......@@ -66,19 +65,30 @@ Foam::point Foam::CatmullRomSpline::position
const scalar mu
) const
{
// out-of-bounds
if (segment < 0)
{
return points().first();
}
else if (segment > nSegments())
{
return points().last();
}
const point& p0 = points()[segment];
const point& p1 = points()[segment+1];
// special cases - no calculation needed
if (segment < 0 || mu < 0.0)
if (mu <= 0.0)
{
return p0;
}
else if (segment > nSegments() || mu >= 1.0)
else if (mu >= 1.0)
{
return p1;
}
// determine the end points
point e0;
point e1;
......@@ -86,7 +96,7 @@ Foam::point Foam::CatmullRomSpline::position
if (segment == 0)
{
// end: simple reflection
e0 = 2.0 * p0 - p1;
e0 = 2*p0 - p1;
}
else
{
......@@ -96,7 +106,7 @@ Foam::point Foam::CatmullRomSpline::position
if (segment+1 == nSegments())
{
// end: simple reflection
e1 = 2.0 * p1 - p0;
e1 = 2*p1 - p0;
}
else
{
......