Commit 1d08ed9b authored by Mark Olesen's avatar Mark Olesen
Browse files

ENH: support arc edge specification with origin point

- The arc will frequently enclose an angle less than 180 degrees.
  For the case, it is possible to define the arc by its endpoints
  and its centre (origin) point. For example,

      arc 0 1 origin (0 0 0);

  When defined in the way, any discrepancy in the arc radius for the
  endpoints is resolved by adjusting the origin to ensure that the
  average radius is satisfied.

  It is also possible to specify a \em flatness factor as a multiplier
  of the radius. For example,

      arc 0 1 origin 1.1 (0 0 0);

ENH: minor code cleanup for block edges

ENH: expose point appending as polyList::concat
parent 8939a556
......@@ -134,7 +134,7 @@ Foam::point Foam::BSpline::position
Foam::scalar Foam::BSpline::length() const
{
NotImplemented;
return 1.0;
return 1;
}
......
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -77,21 +78,12 @@ class BSpline
:
public polyLine
{
// Private Member Functions
//- No copy construct
BSpline(const BSpline&) = delete;
//- No copy assignment
void operator=(const BSpline&) = delete;
public:
// Constructors
//- Construct from components
BSpline
explicit BSpline
(
const pointField& knots,
const bool notImplementedClosed = false
......@@ -100,15 +92,16 @@ public:
// Member Functions
//- Return the point position corresponding to the curve parameter
//- The point position corresponding to the global 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
//- The point position corresponding to the local lambda (0-1)
//- on the given segment
point position(const label segment, const scalar lambda) const;
//- Return the length of the curve
//- The length of the curve
// \note NotImplemented
scalar length() const;
};
......
......@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2014-2016 OpenFOAM Foundation
Copyright (C) 2019 OpenCFD Ltd.
Copyright (C) 2019-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -27,9 +27,9 @@ License
\*---------------------------------------------------------------------------*/
#include "BSplineEdge.H"
#include "polyLine.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
......@@ -59,7 +59,10 @@ Foam::blockEdges::BSplineEdge::BSplineEdge
)
:
blockEdge(points, start, end),
BSpline(appendEndPoints(points, start, end, internalPoints))
BSpline
(
polyLine::concat(points[start_], internalPoints, points[end_])
)
{}
......@@ -73,13 +76,16 @@ Foam::blockEdges::BSplineEdge::BSplineEdge
)
:
blockEdge(dict, index, points, is),
BSpline(appendEndPoints(points, start_, end_, pointField(is)))
BSpline
(
polyLine::concat(points[start_], pointField(is), points[end_])
)
{
token t(is);
is.putBack(t);
token tok(is);
is.putBack(tok);
// discard unused start/end tangents
if (t == token::BEGIN_LIST)
// Discard unused start/end tangents
if (tok == token::BEGIN_LIST)
{
vector tangent0Ignored(is);
vector tangent1Ignored(is);
......
......@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2014-2016 OpenFOAM Foundation
Copyright (C) 2019 OpenCFD Ltd.
Copyright (C) 2019-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -77,20 +77,20 @@ public:
//- Construct from components
BSplineEdge
(
const pointField&,
const label start,
const label end,
const pointField& points, //!< Referenced point field
const label start, //!< Start point in referenced point field
const label end, //!< End point in referenced point field
const pointField& internalPoints
);
//- Construct from Istream, setting pointsList
//- Construct from Istream and point field.
BSplineEdge
(
const dictionary& dict,
const label index,
const searchableSurfaces& geometry,
const pointField&,
Istream&
const pointField& points, //!< Referenced point field
Istream& is
);
......@@ -100,11 +100,12 @@ public:
// Member Functions
//- Return the point position corresponding to the curve parameter
//- The point position corresponding to the curve parameter
// 0 <= lambda <= 1
virtual point position(const scalar) const;
//- Return the length of the spline curve (not implemented)
//- The length of the spline curve
// \note NotImplemented
virtual scalar length() const;
};
......
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -43,10 +44,15 @@ namespace blockEdges
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::coordSystem::cylindrical Foam::blockEdges::arcEdge::calcAngle()
void Foam::blockEdges::arcEdge::calcFromMidPoint
(
const point& p1,
const point& p3,
const point& p2
)
{
const vector a = p2_ - p1_;
const vector b = p3_ - p1_;
const vector a = p2 - p1;
const vector b = p3 - p1;
// Find centre of arcEdge
const scalar asqr = a & a;
......@@ -55,7 +61,7 @@ Foam::coordSystem::cylindrical Foam::blockEdges::arcEdge::calcAngle()
const scalar denom = asqr*bsqr - adotb*adotb;
if (mag(denom) < VSMALL)
if (mag(denom) < ROOTVSMALL)
{
FatalErrorInFunction
<< denom
......@@ -64,63 +70,157 @@ Foam::coordSystem::cylindrical Foam::blockEdges::arcEdge::calcAngle()
const scalar fact = 0.5*(bsqr - adotb)/denom;
point centre = 0.5*a + fact*((a ^ b) ^ a);
const point centre = p1 + 0.5*a + fact*((a ^ b) ^ a);
// Position vectors from centre
const vector r1(p1 - centre);
const vector r2(p2 - centre);
const vector r3(p3 - centre);
centre += p1_;
const scalar mag1(mag(r1));
const scalar mag3(mag(r3));
// Find position vectors w.r.t. the arcEdge centre
const vector r1(p1_ - centre);
const vector r2(p2_ - centre);
const vector r3(p3_ - centre);
vector arcAxis(r1 ^ r3);
// The radius from r1 and from r3 will be identical
radius_ = mag(r3);
// Find angle (in degrees)
angle_ = radToDeg(acos((r3 & r1)/(mag(r3) * mag(r1))));
// Determine the angle
angle_ = acos((r1 & r3)/(mag1*mag3));
// Check if the vectors define an exterior or an interior arcEdge
if (((r1 ^ r2) & (r1 ^ r3)) < 0.0)
{
angle_ = 360.0 - angle_;
angle_ = constant::mathematical::twoPi - angle_;
}
vector arcAxis;
if (angle_ <= 180.0)
if (angle_ <= constant::mathematical::pi)
{
arcAxis = r1 ^ r3;
if (mag(arcAxis)/(mag(r1)*mag(r3)) < 0.001)
if (mag(arcAxis)/(mag1*mag3) < 0.001)
{
arcAxis = r1 ^ r2;
}
}
else
{
arcAxis = r3 ^ r1;
arcAxis = -arcAxis;
}
radius_ = mag(r3);
// Corresponding local cylindrical coordinate system
cs_ = coordSystem::cylindrical(centre, arcAxis, r1);
}
void Foam::blockEdges::arcEdge::calcFromCentre
(
const point& p1,
const point& p3,
const point& centre,
bool adjustCentre,
scalar rMultiplier
)
{
// Position vectors from centre
const vector r1(p1 - centre);
const vector r3(p3 - centre);
const scalar mag1(mag(r1));
const scalar mag3(mag(r3));
const vector chord(p3 - p1);
// The corresponding local cylindrical coordinate system (radians)
return coordSystem::cylindrical("arc", centre, arcAxis, r1);
const vector arcAxis(r1 ^ r3);
// The average radius
radius_ = 0.5*(mag1 + mag3);
// The included angle
angle_ = acos((r1 & r3)/(mag1*mag3));
// TODO? check for 180 degrees (co-linear points)?
bool needsAdjust = false;
if (adjustCentre)
{
needsAdjust = !equal(mag1, mag3);
if (!equal(rMultiplier, 1))
{
// The min radius is constrained by the chord,
// otherwise bad things will happen.
needsAdjust = true;
radius_ *= rMultiplier;
radius_ = max(radius_, (1.001*0.5*mag(chord)));
}
}
if (needsAdjust)
{
// The centre is not equidistant to p1 and p3.
// Use the chord and the arcAxis to determine the vector to
// the midpoint of the chord and adjust the centre along this
// line.
const point newCentre =
(
(0.5 * (p3 + p1)) // mid-chord point
+ sqrt(sqr(radius_) - 0.25 * magSqr(chord))
* normalised(arcAxis ^ chord) // mid-chord -> centre
);
//// Info<< nl << "Adjust centre. r1=" << mag1 << " r3=" << mag3
//// << " radius=" << radius_ << nl
//// << "angle=" << radToDeg(angle_) << ' '
//// << coordSystem::cylindrical(centre, arcAxis, r1) << nl;
// Recalculate - do attempt to readjust
calcFromCentre(p1, p3, newCentre, false);
}
else
{
// Corresponding local cylindrical coordinate system
cs_ = coordSystem::cylindrical(centre, arcAxis, r1);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::blockEdges::arcEdge::arcEdge
(
const pointField& points,
const point& origin,
const label start,
const label end
)
:
blockEdge(points, start, end),
radius_(0),
angle_(0),
cs_()
{
calcFromCentre(points[start_], points[end_], origin);
}
Foam::blockEdges::arcEdge::arcEdge
(
const pointField& points,
const label start,
const label end,
const point& pMid
const point& midPoint
)
:
blockEdge(points, start, end),
p1_(points_[start_]),
p2_(pMid),
p3_(points_[end_]),
cs_(calcAngle())
{}
radius_(0),
angle_(0),
cs_()
{
calcFromMidPoint(points[start_], points[end_], midPoint);
}
Foam::blockEdges::arcEdge::arcEdge
......@@ -133,41 +233,82 @@ Foam::blockEdges::arcEdge::arcEdge
)
:
blockEdge(dict, index, points, is),
p1_(points_[start_]),
p2_(is),
p3_(points_[end_]),
cs_(calcAngle())
{}
radius_(0),
angle_(0),
cs_()
{
point p;
token tok(is);
if (tok.isWord())
{
// Can be
// - origin (0 0 0)
// - origin 1.2 (0 0 0)
scalar rMultiplier = 1;
is >> tok;
if (tok.isNumber())
{
rMultiplier = tok.number();
}
else
{
is.putBack(tok);
}
is >> p; // The origin (centre)
calcFromCentre(points_[start_], points_[end_], p, true, rMultiplier);
}
else
{
is.putBack(tok);
is >> p; // A mid-point
calcFromMidPoint(points_[start_], points_[end_], p);
}
// Debug information
#if 0
Info<< "arc " << start_ << ' ' << end_
<< ' ' << position(0.5) << ' ' << cs_
// << " radius=" << radius_ << " angle=" << radToDeg(angle_)
<< nl;
#endif
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::point Foam::blockEdges::arcEdge::position(const scalar lambda) const
{
#ifdef FULLDEBUG
if (lambda < -SMALL || lambda > 1 + SMALL)
{
FatalErrorInFunction
<< "Parameter out of range, lambda = " << lambda
<< abort(FatalError);
WarningInFunction
<< "Parameter out of range, lambda = " << lambda << nl;
}
#endif
if (lambda < SMALL)
{
return p1_;
return points_[start_];
}
else if (lambda > 1 - SMALL)
else if (lambda >= 1 - SMALL)
{
return p3_;
return points_[end_];
}
// The angle is degrees, the coordinate system in radians
return cs_.globalPosition(vector(radius_, degToRad(lambda*angle_), 0));
return cs_.globalPosition(vector(radius_, (lambda*angle_), 0));
}
Foam::scalar Foam::blockEdges::arcEdge::length() const
Foam::scalar Foam::blockEdges::arcEdge::length() const noexcept
{
return degToRad(radius_*angle_);
return (radius_*angle_);
}
......
......@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2017 OpenCFD Ltd.
Copyright (C) 2017-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -28,7 +28,32 @@ Class
Foam::blockEdges::arcEdge
Description
Defines the arcEdge of a circle in terms of 3 points on its circumference
A blockEdge defined as an arc of a circle.
The arc is normally defined by its endpoints and a point on
its circumference, typically a midpoint. For example,
\verbatim
points ((1 0 0) (0 1 0));
arc 0 1 (0.707107 0.707107 0);
\endverbatim
The arc can enclose an angle greater than 0 and less than 360 degrees.
The arc will frequently enclose an angle less than 180 degrees.
For the case, it is possible to define the arc by its endpoints and its
centre (origin) point. For example,
\verbatim
arc 0 1 origin (0 0 0);
\endverbatim
When defined in the way, any discrepancy in the arc radius for the
endpoints is resolved by adjusting the origin to ensure that the average
radius is satisfied.
It is also possible to define a \em flatness factor as a multiplier of
the calculated radius. For example,
\verbatim
arc 0 1 origin 1.1 (0 0 0);
\endverbatim
SourceFiles
arcEdge.C
......@@ -56,26 +81,41 @@ class arcEdge
:
public blockEdge
{
// Private data
// Begin, mid, end points
point p1_, p2_, p3_;
//- The arc angle (in degrees)
scalar angle_;
// Private Data
//- The arc radius
scalar radius_;
//- The arc angle (radians)
scalar angle_;
//- The local cylindrical coordinate system
coordSystem::cylindrical cs_;
// Private Member Functions
//- Calculate the angle, radius and axis
// \return the cylindrical coordinate system
coordSystem::cylindrical calcAngle();
//- Calculate angle, radius, cylindrical coordinate system
//- from end points and the given point on the circumference
void calcFromMidPoint
(
const point& p1, //!< Start point
const point& p3, //!< End point
const point& p2 //!< Point on circumference
);
//- Calculate angle, radius, cylindrical coordinate system
//- from end points and the given origin.
// Optionally adjust centre to accommodate deviations in the
// effective radius to the end-points
void calcFromCentre
(
const point& p1, //!< Start point
const point& p3, //!< End point
const point& centre, //!< Centre
bool adjustCentre = false, //!< Adjust centre
scalar rMultiplier = 1 //!< Adjust radius by this factor
);
//- No copy construct
arcEdge(const arcEdge&) = delete;
......@@ -92,23 +132,34 @@ public:
// Constructors
//- Construct from components
//- Construct from components, given the origin of the circle
arcEdge
(
const pointField& points, //!< Referenced point field
const point& origin, //!< The origin of the circle
const label start, //!< Start point in referenced point field
const label end //!< End point in referenced point field
);
//- Construct from components, using a point on the circumference
arcEdge
(
const pointField& points,
const label start,
const label end,
const point& pMid
const pointField& points, //!< Referenced point field
const label start, //!< Start point in referenced point field
const label end, //!< End point in referenced point field
const point& midPoint //!< A point on the circumference
);
//- Construct from Istream setting pointsList
//- Construct from Istream and point field.
// The Istream can specify either a point on the circumference,
// or with a tag to specify the origin.
arcEdge
(
const dictionary& dict,
const label index,
const searchableSurfaces& geometry,
const pointField& points,
Istream&
const searchableSurfaces& geometry, // unsed