blockMesh: Angle-and-axis specification for arc edges

Arc-edges can now be specified with a sector angle (in degrees) and an
axis of the circle of which the arc forms a part. The new syntax is as
follows:

   edges
   (
       arc <vertex-0> <vertex-1> <angle> (<axis-x> <axis-y> <axis-z>)
   );

This is often more convenient than the alternative specification where a
third third point somewhere in the arc is given; it usually does not
require any additional calculation on the part of the user, and multiple
entries are very likely to be identical.

Which specification is used depends on the form of the entry that comes
after the two vertices. If the entry is a vector then it is assumed to
be a point in the arc; if it is scalar then is is taken to be the angle
and the axis is assumed to follow.

For example, to put a 90 degree arc between the vertices 12 and 13, at
(1 0 0) and (0 1 0) respectively, the following specification can now be
used:

   edges
   (
       arc 12 13 90.0 (0 0 1)
   );

This is equivalent to the existing point-in-arc speficiation below:

   edges
   (
       arc 12 13 (0.707107 0.707107 0)
   );

An edge's points are ordered on the perimeter of the circle according to
a right-hand screw rule on the given axis. Changing the the side of the
edge on which the arc is defined can therefore be achieved by reversing
either the edge or the direction of the axis.

If the given axis is not perpendicular to the line between the vertices,
then the arc gains some axial length and becomes a helix.

src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -41,87 +41,66 @@ namespace blockEdges
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::cylindricalCS Foam::blockEdges::arcEdge::calcAngle()
+void Foam::blockEdges::arcEdge::calc(const point& pA)
 {
-    vector a = p2_ - p1_;
+    const tensor dps(pA - p0_, pA - p1_, (pA - p0_)^(pA - p1_));
-    vector b = p3_ - p1_;
+    const tensor pMs((p0_ + pA)/2, (p1_ + pA)/2, pA);
 
-    // find centre of arcEdge
+    if (det(dps) < vSmall)
-    scalar asqr = a & a;
-    scalar bsqr = b & b;
-    scalar adotb = a & b;
-
-    scalar denom = asqr*bsqr - adotb*adotb;
-
-    if (mag(denom) < vSmall)
     {
         FatalErrorInFunction
-            << denom
+            << "Arc point " << pA << " lies on the line of edge "
+            << edge(start_, end_) << abort(FatalError);
+    }
+
+    const point c =
+        inv(dps)
+      & vector(dps.x() & pMs.x(), dps.y() & pMs.y(), dps.z() & pMs.z());
+
+    const vector r0 = p0_ - c, r1 = p1_ - c;
+    const scalar cosT = (r0 & r1)/(mag(r0)*mag(r1));
+    scalar t = acos(max(-1, min(cosT, 1)));
+    if (((r0 ^ r1) & dps.z()) > 0)
+    {
+        t = 2*constant::mathematical::pi - t;
+    }
+
+    centre_ = c;
+    axis_ = - normalised(dps.z());
+    theta_ = t;
+    length_ = 0;
+}
+
+
+void Foam::blockEdges::arcEdge::calc(const scalar theta, const vector& axis)
+{
+    if (0 >= theta || theta >= 360)
+    {
+        FatalErrorInFunction
+            << "Arc angle for edge " << edge(start_, end_)
+            << " must take a value between 0 and 360 degrees"
             << abort(FatalError);
     }
 
-    scalar fact = 0.5*(bsqr - adotb)/denom;
+    const vector dp = p1_ - p0_;
 
-    point centre = 0.5*a + fact*((a ^ b) ^ a);
+    const vector pM = (p0_ + p1_)/2;
+    const vector rM = normalised(dp ^ axis);
 
-    centre += p1_;
+    const scalar l = dp & axis;
 
-    // find position vectors w.r.t. the arcEdge centre
+    const vector chord = dp - l*axis;
-    vector r1(p1_ - centre);
+    const scalar magChord = mag(chord);
-    vector r2(p2_ - centre);
-    vector r3(p3_ - centre);
 
-    // find angles
+    centre_ = pM - l*axis/2 - rM*magChord/2/tan(degToRad(theta)/2);
-    const scalar cosAngle = (r3 & r1)/(mag(r3) * mag(r1));
+    axis_ = axis;
-    angle_ = radToDeg(acos(max(-1, min(cosAngle, 1))));
+    theta_ = degToRad(theta);
-
+    length_ = l;
-    // check if the vectors define an exterior or an interior arcEdge
-    if (((r1 ^ r2) & (r1 ^ r3)) < 0.0)
-    {
-        angle_ = 360.0 - angle_;
-    }
-
-    vector tempAxis;
-
-    if (angle_ <= 180.0)
-    {
-        tempAxis = r1 ^ r3;
-
-        if (mag(tempAxis)/(mag(r1)*mag(r3)) < 0.001)
-        {
-            tempAxis = r1 ^ r2;
-        }
-    }
-    else
-    {
-        tempAxis = r3 ^ r1;
-    }
-
-    radius_ = mag(r3);
-
-    // set up and return the local coordinate system
-    return cylindricalCS("arcEdgeCS", centre, tempAxis, r1);
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::blockEdges::arcEdge::arcEdge
-(
-    const pointField& points,
-    const label start,
-    const label end,
-    const point& pMid
-)
-:
-    blockEdge(points, start, end),
-    p1_(points_[start_]),
-    p2_(pMid),
-    p3_(points_[end_]),
-    cs_(calcAngle())
-{}
-
-
 Foam::blockEdges::arcEdge::arcEdge
 (
     const dictionary& dict,
@@ -132,42 +111,54 @@ Foam::blockEdges::arcEdge::arcEdge
 )
 :
     blockEdge(dict, index, points, is),
-    p1_(points_[start_]),
+    p0_(points_[start_]),
-    p2_(is),
+    p1_(points_[end_]),
-    p3_(points_[end_]),
+    centre_(),
-    cs_(calcAngle())
+    axis_(),
-{}
+    theta_(),
+    length_()
+{
+    token firstToken(is);
+    is.putBack(firstToken);
+
+    if (firstToken == token::BEGIN_LIST)
+    {
+        const point pA(is);
+        calc(pA);
+    }
+    else
+    {
+        const scalar theta = readScalar(is);
+        const vector axis(is);
+        calc(theta, normalised(axis));
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::point Foam::blockEdges::arcEdge::position(const scalar lambda) const
 {
-    if (lambda < -small || lambda > 1 + small)
+    if (lambda < - small || lambda > 1 + small)
     {
         FatalErrorInFunction
             << "Parameter out of range, lambda = " << lambda
             << abort(FatalError);
     }
 
-    if (lambda < small)
+    const scalar t = theta_*lambda;
-    {
+    const vector r1 = p0_ - centre_, r2 = axis_ ^ r1;
-        return p1_;
+
-    }
+    return centre_ + r1*cos(t) + r2*sin(t) + axis_*length_*lambda;
-    else if (lambda > 1 - small)
-    {
-        return p3_;
-    }
-    else
-    {
-        return cs_.globalPosition(vector(radius_, lambda*angle_, 0.0));
-    }
 }
 
 
 Foam::scalar Foam::blockEdges::arcEdge::length() const
 {
-    return degToRad(angle_*radius_);
+    const vector r1 = p0_ - centre_;
+
+    // Length of a helical segment
+    return degToRad(theta_*sqrt(magSqr(r1) + sqr(length_)));
 }
 
 

src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H

@@ -25,7 +25,14 @@ Class
     Foam::blockEdges::arcEdge
 
 Description
-    Defines the arcEdge of a circle in terms of 3 points on its circumference
+    An arcEdge between two points on a circle. The arc is defined either by a
+    third point that the arc passes through, or by the angle of the sector and
+    the axis of the circle.
+
+    If the angle-axis specification is used, the axis may not be not
+    perpendicular to the line connecting the end points. In that case, the axis
+    is considered to be that of a cylinder, and the arc represents a portion of
+    a helix on the surface of that cylinder.
 
 SourceFiles
     arcEdge.C
@@ -36,7 +43,6 @@ SourceFiles
 #define blockEdges_arcEdge_H
 
 #include "blockEdge.H"
-#include "cylindricalCS.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -55,15 +61,34 @@ class arcEdge
 {
     // Private Data
 
-        point p1_, p2_, p3_;
+        //- The start point of the arc
-        cylindricalCS cs_;
+        const point p0_;
-        scalar angle_;
+
-        scalar radius_;
+        //- The end point of the arc
+        const point p1_;
+
+        //- The centre of the circle of which the arc is a part
+        point centre_;
+
+        //- The axis of the circle of which the arc is a part
+        vector axis_;
+
+        //- The sector angle of the arc
+        scalar theta_;
+
+        //- The axial length of the arc/helix
+        scalar length_;
+
 
     // Private Member Functions
 
-        //- Calculate the coordinate system, angle and radius
+        //- Calculate the coordinate system, angle and radius from a mid point
-        cylindricalCS calcAngle();
+        //  somewhere within the arc
+        void calc(const point& pM);
+
+        //- Calculate the coordinate system, angle and radius from a sector
+        //  angle and a normal to the plane of the arc
+        void calc(const scalar theta, const vector& axis);
 
 
 public:
@@ -74,14 +99,6 @@ public:
 
     // Constructors
 
-        //- Construct from components
-        arcEdge
-        (
-            const pointField& points,
-            const label start, const label end,
-            const point& pMid
-        );
-
         //- Construct from Istream setting pointsList
         arcEdge
         (

tutorials/basic/potentialFoam/cylinder/system/blockMeshDict

@@ -76,22 +76,22 @@ blocks
 
 edges
 (
-    arc 0 5 (0.469846 0.17101 -0.5)
+    arc 0 5 45.0 (0 0 1)
-    arc 5 10 (0.17101 0.469846 -0.5)
+    arc 5 10 45.0 (0 0 1)
-    arc 1 4 (0.939693 0.34202 -0.5)
+    arc 1 4 45.0 (0 0 1)
-    arc 4 9 (0.34202 0.939693 -0.5)
+    arc 4 9 45.0 (0 0 1)
-    arc 19 24 (0.469846 0.17101 0.5)
+    arc 19 24 45.0 (0 0 1)
-    arc 24 29 (0.17101 0.469846 0.5)
+    arc 24 29 45.0 (0 0 1)
-    arc 20 23 (0.939693 0.34202 0.5)
+    arc 20 23 45.0 (0 0 1)
-    arc 23 28 (0.34202 0.939693 0.5)
+    arc 23 28 45.0 (0 0 1)
-    arc 11 16 (-0.469846 0.17101 -0.5)
+    arc 11 16 45.0 (0 0 -1)
-    arc 16 10 (-0.17101 0.469846 -0.5)
+    arc 16 10 45.0 (0 0 -1)
-    arc 12 15 (-0.939693 0.34202 -0.5)
+    arc 12 15 45.0 (0 0 -1)
-    arc 15 9 (-0.34202 0.939693 -0.5)
+    arc 15 9 45.0 (0 0 -1)
-    arc 30 35 (-0.469846 0.17101 0.5)
+    arc 30 35 45.0 (0 0 -1)
-    arc 35 29 (-0.17101 0.469846 0.5)
+    arc 35 29 45.0 (0 0 -1)
-    arc 31 34 (-0.939693 0.34202 0.5)
+    arc 31 34 45.0 (0 0 -1)
-    arc 34 28 (-0.34202 0.939693 0.5)
+    arc 34 28 45.0 (0 0 -1)
 );
 
 boundary