From 73d253c34b3e184802efb316f996f244cc795ec6 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Fri, 28 Jun 2019 12:17:13 +0100 Subject: [PATCH] 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 ( ) ); 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. --- .../blockMesh/blockEdges/arcEdge/arcEdge.C | 161 +++++++++--------- .../blockMesh/blockEdges/arcEdge/arcEdge.H | 49 ++++-- .../cylinder/system/blockMeshDict | 32 ++-- 3 files changed, 125 insertions(+), 117 deletions(-) diff --git a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C index b2f5c57b66..88bb8ebd7d 100644 --- a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C +++ b/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_; - vector b = p3_ - p1_; + const tensor dps(pA - p0_, pA - p1_, (pA - p0_)^(pA - p1_)); + const tensor pMs((p0_ + pA)/2, (p1_ + pA)/2, pA); - // find centre of arcEdge - scalar asqr = a & a; - scalar bsqr = b & b; - scalar adotb = a & b; - - scalar denom = asqr*bsqr - adotb*adotb; - - if (mag(denom) < vSmall) + if (det(dps) < 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 - vector r1(p1_ - centre); - vector r2(p2_ - centre); - vector r3(p3_ - centre); + const vector chord = dp - l*axis; + const scalar magChord = mag(chord); - // find angles - const scalar cosAngle = (r3 & r1)/(mag(r3) * mag(r1)); - angle_ = radToDeg(acos(max(-1, min(cosAngle, 1)))); - - // 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); + centre_ = pM - l*axis/2 - rM*magChord/2/tan(degToRad(theta)/2); + axis_ = axis; + theta_ = degToRad(theta); + length_ = l; } // * * * * * * * * * * * * * * * * 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_]), - p2_(is), - p3_(points_[end_]), - cs_(calcAngle()) -{} + p0_(points_[start_]), + p1_(points_[end_]), + centre_(), + 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) - { - return p1_; - } - else if (lambda > 1 - small) - { - return p3_; - } - else - { - return cs_.globalPosition(vector(radius_, lambda*angle_, 0.0)); - } + const scalar t = theta_*lambda; + const vector r1 = p0_ - centre_, r2 = axis_ ^ r1; + + return centre_ + r1*cos(t) + r2*sin(t) + axis_*length_*lambda; } 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_))); } diff --git a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H index 9cfdaed3fc..a9886bb573 100644 --- a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H +++ b/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_; - cylindricalCS cs_; - scalar angle_; - scalar radius_; + //- The start point of the arc + const point p0_; + + //- 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 - cylindricalCS calcAngle(); + //- Calculate the coordinate system, angle and radius from a mid point + // 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 ( diff --git a/tutorials/basic/potentialFoam/cylinder/system/blockMeshDict b/tutorials/basic/potentialFoam/cylinder/system/blockMeshDict index 8aace2d3c7..7766ff8e15 100644 --- a/tutorials/basic/potentialFoam/cylinder/system/blockMeshDict +++ b/tutorials/basic/potentialFoam/cylinder/system/blockMeshDict @@ -76,22 +76,22 @@ blocks edges ( - arc 0 5 (0.469846 0.17101 -0.5) - arc 5 10 (0.17101 0.469846 -0.5) - arc 1 4 (0.939693 0.34202 -0.5) - arc 4 9 (0.34202 0.939693 -0.5) - arc 19 24 (0.469846 0.17101 0.5) - arc 24 29 (0.17101 0.469846 0.5) - arc 20 23 (0.939693 0.34202 0.5) - arc 23 28 (0.34202 0.939693 0.5) - arc 11 16 (-0.469846 0.17101 -0.5) - arc 16 10 (-0.17101 0.469846 -0.5) - arc 12 15 (-0.939693 0.34202 -0.5) - arc 15 9 (-0.34202 0.939693 -0.5) - arc 30 35 (-0.469846 0.17101 0.5) - arc 35 29 (-0.17101 0.469846 0.5) - arc 31 34 (-0.939693 0.34202 0.5) - arc 34 28 (-0.34202 0.939693 0.5) + arc 0 5 45.0 (0 0 1) + arc 5 10 45.0 (0 0 1) + arc 1 4 45.0 (0 0 1) + arc 4 9 45.0 (0 0 1) + arc 19 24 45.0 (0 0 1) + arc 24 29 45.0 (0 0 1) + arc 20 23 45.0 (0 0 1) + arc 23 28 45.0 (0 0 1) + arc 11 16 45.0 (0 0 -1) + arc 16 10 45.0 (0 0 -1) + arc 12 15 45.0 (0 0 -1) + arc 15 9 45.0 (0 0 -1) + arc 30 35 45.0 (0 0 -1) + arc 35 29 45.0 (0 0 -1) + arc 31 34 45.0 (0 0 -1) + arc 34 28 45.0 (0 0 -1) ); boundary