quickly implemented BSpline2 as an alternative B-Spline implementation

- also looks reasonable and doesn't deviate much from Catmull-Rom

applications/test/spline/splineTest.C

@@ -27,7 +27,9 @@ License
 
 #include "vector.H"
 #include "IFstream.H"
+
 #include "BSpline.H"
+#include "BSpline2.H"
 #include "CatmullRomSpline.H"
 
 using namespace Foam;
@@ -39,8 +41,9 @@ int main(int argc, char *argv[])
 {
     argList::noParallel();
     argList::validArgs.insert("file .. fileN");
-    argList::addBoolOption("B", "B-Spline");
-    argList::addBoolOption("cmr", "catmull-rom spline (default)");
+    argList::addBoolOption("Bold", "B-Spline implementation OLD");
+    argList::addBoolOption("Bnew", "B-Spline implementation NEW");
+    argList::addBoolOption("CMR", "catmull-rom spline (default)");
     argList::addOption
     (
         "n",
@@ -55,11 +58,12 @@ int main(int argc, char *argv[])
         args.printUsage();
     }
 
-    bool useBSpline    = args.optionFound("B");
-    bool useCatmullRom = args.optionFound("cmr");
+    bool useBSplineOld = args.optionFound("Bold");
+    bool useBSplineNew = args.optionFound("Bnew");
+    bool useCatmullRom = args.optionFound("CMR");
     label nSeg = args.optionLookupOrDefault<label>("n", 20);
 
-    if (!useBSpline && !useCatmullRom)
+    if (!useCatmullRom && !useBSplineOld && !useBSplineNew)
     {
         Info<<"defaulting to Catmull-Rom spline" << endl;
         useCatmullRom = true;
@@ -77,12 +81,27 @@ int main(int argc, char *argv[])
         {
             Info<<"\noriginal points: " << pointFields[splineI] << nl;
 
-            if (useBSpline)
+            if (useBSplineOld)
             {
-                BSpline spl(pointFields[splineI], vector::zero, vector::zero);
+                BSpline spl(pointFields[splineI]);
 
                 Info<< nl
-                    << "B-Spline interpolation:" << nl
+                    << "B-Spline interpolation: OLD" << nl
+                    << "----------------------" << endl;
+
+                for (label segI = 0; segI <= nSeg; ++segI)
+                {
+                    scalar lambda = scalar(segI)/scalar(nSeg);
+                    Info<< spl.position(lambda) << "    // " << lambda << endl;
+                }
+            }
+
+            if (useBSplineNew)
+            {
+                BSpline2 spl(pointFields[splineI]);
+
+                Info<< nl
+                    << "B-Spline interpolation: NEW" << nl
                     << "----------------------" << endl;
 
                 for (label segI = 0; segI <= nSeg; ++segI)

src/mesh/blockMesh/Make/files

@@ -1,3 +1,4 @@
+curvedEdges/BSpline2.C
 curvedEdges/CatmullRomSpline.C
 curvedEdges/polyLine.C
 

src/mesh/blockMesh/curvedEdges/BSpline2.C

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "BSpline2.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::BSpline2::BSpline2
+(
+    const pointField& Knots,
+    const vector&,
+    const vector&
+)
+:
+    polyLine(Knots)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::point Foam::BSpline2::position(const scalar mu) const
+{
+    // endpoints
+    if (mu < SMALL)
+    {
+        return points().first();
+    }
+    else if (mu > 1 - SMALL)
+    {
+        return points().last();
+    }
+
+    scalar lambda = mu;
+    label segment = localParameter(lambda);
+    return position(segment, lambda);
+}
+
+
+Foam::point Foam::BSpline2::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 1.0/6.0 *
+    (
+        ( e0 + 4*p0 + p1 )
+      + mu *
+        (
+            ( -3*e0 + 3*p1 )
+          + mu *
+            (
+                ( 3*e0 - 6*p0 + 3*p1 )
+              + mu *
+                ( -e0 + 3*p0 - 3*p1 + e1 )
+            )
+        )
+    );
+}
+
+
+Foam::scalar Foam::BSpline2::length() const
+{
+    notImplemented("BSpline2::length() const");
+    return 1.0;
+}
+
+
+// ************************************************************************* //

src/mesh/blockMesh/curvedEdges/BSpline2.H

@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+// ************************************************************************* //

src/mesh/blockMesh/curvedEdges/CatmullRomSpline.C

@@ -106,7 +106,7 @@ Foam::point Foam::CatmullRomSpline::position
 
     return 0.5 *
     (
-        ( 2 * p0 )
+        ( 2*p0 )
       + mu *
         (
             ( -e0 + p1 )