diff --git a/src/mesh/blockMesh/Make/files b/src/mesh/blockMesh/Make/files index cc7a1c1bb5..d68caa1c24 100644 --- a/src/mesh/blockMesh/Make/files +++ b/src/mesh/blockMesh/Make/files @@ -1,13 +1,17 @@ +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 diff --git a/src/mesh/blockMesh/curvedEdges/CatmullRomSpline.C b/src/mesh/blockMesh/curvedEdges/CatmullRomSpline.C new file mode 100644 index 0000000000..f82d3f0d58 --- /dev/null +++ b/src/mesh/blockMesh/curvedEdges/CatmullRomSpline.C @@ -0,0 +1,131 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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; +} + + +// ************************************************************************* // diff --git a/src/mesh/blockMesh/curvedEdges/CatmullRomSpline.H b/src/mesh/blockMesh/curvedEdges/CatmullRomSpline.H new file mode 100644 index 0000000000..6f1851a7a6 --- /dev/null +++ b/src/mesh/blockMesh/curvedEdges/CatmullRomSpline.H @@ -0,0 +1,128 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + +// ************************************************************************* // diff --git a/src/mesh/blockMesh/curvedEdges/arcEdge.C b/src/mesh/blockMesh/curvedEdges/arcEdge.C index 5718b731ad..21143ac43a 100644 --- a/src/mesh/blockMesh/curvedEdges/arcEdge.C +++ b/src/mesh/blockMesh/curvedEdges/arcEdge.C @@ -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); } diff --git a/src/mesh/blockMesh/curvedEdges/arcEdge.H b/src/mesh/blockMesh/curvedEdges/arcEdge.H index 1539916aaa..194e5c634b 100644 --- a/src/mesh/blockMesh/curvedEdges/arcEdge.H +++ b/src/mesh/blockMesh/curvedEdges/arcEdge.H @@ -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&); diff --git a/src/mesh/blockMesh/curvedEdges/curvedEdge.C b/src/mesh/blockMesh/curvedEdges/curvedEdge.C index 595a923475..35e58b7692 100644 --- a/src/mesh/blockMesh/curvedEdges/curvedEdge.C +++ b/src/mesh/blockMesh/curvedEdges/curvedEdge.C @@ -108,7 +108,7 @@ Foam::autoPtr Foam::curvedEdge::New // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::pointField Foam::curvedEdge::fullKnotList +Foam::pointField Foam::curvedEdge::appendEndPoints ( const pointField& points, const label start, diff --git a/src/mesh/blockMesh/curvedEdges/curvedEdge.H b/src/mesh/blockMesh/curvedEdges/curvedEdge.H index 629b833cb9..bd1ad3246a 100644 --- a/src/mesh/blockMesh/curvedEdges/curvedEdge.H +++ b/src/mesh/blockMesh/curvedEdges/curvedEdge.H @@ -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, diff --git a/src/mesh/blockMesh/curvedEdges/BSpline.C b/src/mesh/blockMesh/curvedEdges/legacy/BSpline.C similarity index 100% rename from src/mesh/blockMesh/curvedEdges/BSpline.C rename to src/mesh/blockMesh/curvedEdges/legacy/BSpline.C diff --git a/src/mesh/blockMesh/curvedEdges/BSpline.H b/src/mesh/blockMesh/curvedEdges/legacy/BSpline.H similarity index 100% rename from src/mesh/blockMesh/curvedEdges/BSpline.H rename to src/mesh/blockMesh/curvedEdges/legacy/BSpline.H diff --git a/src/mesh/blockMesh/curvedEdges/polySplineEdge.C b/src/mesh/blockMesh/curvedEdges/legacy/polySplineEdge.C similarity index 93% rename from src/mesh/blockMesh/curvedEdges/polySplineEdge.C rename to src/mesh/blockMesh/curvedEdges/legacy/polySplineEdge.C index 94ef3feb70..bf370cd7d3 100644 --- a/src/mesh/blockMesh/curvedEdges/polySplineEdge.C +++ b/src/mesh/blockMesh/curvedEdges/legacy/polySplineEdge.C @@ -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 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; + while (param_[segmentI] < lambda) { - i++; + ++segmentI; } - i--; // we now want the lower bound + --segmentI; // we now want the lower bound // linear interpolation return ( - controlPoints_[i] - + ( controlPoints_[i+1] - controlPoints_[i] ) - * ( lambda - distances_[i] ) / ( distances_[i+1] - distances_[i] ) + points_[segmentI] + + ( points_[segmentI+1] - points_[segmentI] ) + * ( lambda - param_[segmentI] ) + / ( param_[segmentI+1] - param_[segmentI] ) ); } diff --git a/src/mesh/blockMesh/curvedEdges/polyLine.H b/src/mesh/blockMesh/curvedEdges/polyLine.H index e549b8fb85..ea8c59e389 100644 --- a/src/mesh/blockMesh/curvedEdges/polyLine.H +++ b/src/mesh/blockMesh/curvedEdges/polyLine.H @@ -26,11 +26,8 @@ Class Foam::polyLine Description - Define a series of control points, which can also be interpreted as a - series of straight line segments. - - This is the basic polyLine class which implements just the line - (no topology - it is not derived from curvedEdge) + A series of straight line segments, which can also be interpreted as + a series of control points for splines, etc. SourceFiles polyLine.C @@ -67,19 +64,25 @@ protected: // Protected data - //- The control points or ends of each segmen - pointField controlPoints_; - - //- The rational (0-1) cumulative distance for each control-point - scalarList distances_; + //- The control points or ends of each segments + pointField points_; //- The real line length scalar lineLength_; + //- The rational (0-1) cumulative parameter value for each point + scalarList param_; + // Protected member functions - //- Precalculate the rational cumulative distances and the line-length - void calcDistances(); + //- Precalculate the rational cumulative parameter value + // and the line-length + void calcParam(); + + + //- Return the line segment and the local parameter [0..1] + // corresponding to the global lambda [0..1] + label localParameter(scalar& lambda) const; public: @@ -92,7 +95,10 @@ public: // Member Functions //- Return const-access to the control-points - const pointField& controlPoints() const; + const pointField& points() const; + + //- Return the number of line segments + label nSegments() const; //- Return the point position corresponding to the curve parameter // 0 <= lambda <= 1 diff --git a/src/mesh/blockMesh/curvedEdges/polyLineEdge.C b/src/mesh/blockMesh/curvedEdges/polyLineEdge.C index 522c6abe48..1f038a78a8 100644 --- a/src/mesh/blockMesh/curvedEdges/polyLineEdge.C +++ b/src/mesh/blockMesh/curvedEdges/polyLineEdge.C @@ -48,14 +48,20 @@ Foam::polyLineEdge::polyLineEdge ) : curvedEdge(ps, start, end), - polyLine(fullKnotList(ps, start_, end_, otherPoints)) + polyLine(appendEndPoints(ps, start_, end_, otherPoints)) {} Foam::polyLineEdge::polyLineEdge(const pointField& ps, Istream& is) : curvedEdge(ps, is), - polyLine(fullKnotList(ps, start_, end_, pointField(is))) + polyLine(appendEndPoints(ps, start_, end_, pointField(is))) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::polyLineEdge::~polyLineEdge() {} diff --git a/src/mesh/blockMesh/curvedEdges/polyLineEdge.H b/src/mesh/blockMesh/curvedEdges/polyLineEdge.H index 2c79765a98..e941bba716 100644 --- a/src/mesh/blockMesh/curvedEdges/polyLineEdge.H +++ b/src/mesh/blockMesh/curvedEdges/polyLineEdge.H @@ -83,9 +83,8 @@ public: polyLineEdge(const pointField&, Istream&); - // Destructor - - virtual ~polyLineEdge(){} + //- Destructor + virtual ~polyLineEdge(); // Member Functions diff --git a/src/mesh/blockMesh/curvedEdges/splineEdge.C b/src/mesh/blockMesh/curvedEdges/splineEdge.C new file mode 100644 index 0000000000..6c5df771a9 --- /dev/null +++ b/src/mesh/blockMesh/curvedEdges/splineEdge.C @@ -0,0 +1,86 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "splineEdge.H" +#include "addToRunTimeSelectionTable.H" + + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(splineEdge, 0); + addToRunTimeSelectionTable(curvedEdge, splineEdge, Istream); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::splineEdge::splineEdge +( + const pointField& points, + const label start, + const label end, + const pointField& otherknots +) +: + curvedEdge(points, start, end), + CatmullRomSpline(appendEndPoints(points, start, end, otherknots)) +{} + + +Foam::splineEdge::splineEdge(const pointField& points, Istream& is) +: + curvedEdge(points, is), + CatmullRomSpline(appendEndPoints(points, start_, end_, pointField(is))) +{ + token t(is); + is.putBack(t); + + // might have start/end tangents that we currently ignore + if (t == token::BEGIN_LIST) + { + vector fstend(is); + vector sndend(is); + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::point Foam::splineEdge::position(const scalar mu) const +{ + return CatmullRomSpline::position(mu); +} + + +Foam::scalar Foam::splineEdge::length() const +{ + return CatmullRomSpline::length(); +} + + +// ************************************************************************* // diff --git a/src/mesh/blockMesh/curvedEdges/splineEdge.H b/src/mesh/blockMesh/curvedEdges/splineEdge.H new file mode 100644 index 0000000000..c1ad9f5f76 --- /dev/null +++ b/src/mesh/blockMesh/curvedEdges/splineEdge.H @@ -0,0 +1,111 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + +Class + Foam::splineEdge + +Description + A curvedEdge interface for Catmull-Rom splines. + +SourceFiles + splineEdge.C + +\*---------------------------------------------------------------------------*/ + +#ifndef splineEdge_H +#define splineEdge_H + +#include "curvedEdge.H" +#include "CatmullRomSpline.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class splineEdge Declaration +\*---------------------------------------------------------------------------*/ + +class splineEdge +: + public curvedEdge, + public CatmullRomSpline +{ + // Private Member Functions + + //- Disallow default bitwise copy construct + splineEdge(const splineEdge&); + + //- Disallow default bitwise assignment + void operator=(const splineEdge&); + + +public: + + //- Runtime type information + TypeName("spline"); + + + // Constructors + + //- Construct from components + splineEdge + ( + const pointField&, + const label start, + const label end, + const pointField& otherKnots + ); + + //- Construct from Istream setting pointsList + splineEdge(const pointField&, Istream&); + + + // Destructor + + virtual ~splineEdge() + {} + + + // Member Functions + + //- Return the point position corresponding to the curve parameter + // 0 <= lambda <= 1 + virtual point position(const scalar) const; + + //- Return the length of the simple spline curve + virtual scalar length() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //