src/mesh/blockMesh/curvedEdges/arcEdge.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "arcEdge.H"
  27 #include "unitConversion.H"
  28 #include "addToRunTimeSelectionTable.H"
  29
  30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34     defineTypeNameAndDebug(arcEdge, 0);
  35     addToRunTimeSelectionTable(curvedEdge, arcEdge, Istream);
  36 }
  37
  38
  39 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  40
  41 Foam::cylindricalCS Foam::arcEdge::calcAngle()
  42 {
  43     vector a = p2_ - p1_;
  44     vector b = p3_ - p1_;
  45
  46     // find centre of arcEdge
  47     scalar asqr = a & a;
  48     scalar bsqr = b & b;
  49     scalar adotb = a & b;
  50
  51     scalar denom = asqr*bsqr - adotb*adotb;
  52
  53     if (mag(denom) < VSMALL)
  54     {
  55         FatalErrorIn("cylindricalCS arcEdge::calcAngle()")
  56             << "Invalid arc definition - are the points co-linear?  Denom ="
  57             << denom
  58             << abort(FatalError);
  59     }
  60
  61     scalar fact = 0.5*(bsqr - adotb)/denom;
  62
  63     point centre = 0.5*a + fact*((a ^ b) ^ a);
  64
  65     centre += p1_;
  66
  67     // find position vectors w.r.t. the arcEdge centre
  68     vector r1(p1_ - centre);
  69     vector r2(p2_ - centre);
  70     vector r3(p3_ - centre);
  71
  72     // find angles
  73     angle_ = radToDeg(acos((r3 & r1)/(mag(r3) * mag(r1))));
  74
  75     // check if the vectors define an exterior or an interior arcEdge
  76     if (((r1 ^ r2) & (r1 ^ r3)) < 0.0)
  77     {
  78         angle_ = 360.0 - angle_;
  79     }
  80
  81     vector tempAxis;
  82
  83     if (angle_ <= 180.0)
  84     {
  85         tempAxis = r1 ^ r3;
  86
  87         if (mag(tempAxis)/(mag(r1)*mag(r3)) < 0.001)
  88         {
  89             tempAxis = r1 ^ r2;
  90         }
  91     }
  92     else
  93     {
  94         tempAxis = r3 ^ r1;
  95     }
  96
  97     radius_ = mag(r3);
  98
  99     // set up and return the local coordinate system
 100     return cylindricalCS("arcEdgeCS", centre, tempAxis, r1);
 101 }
 102
 103
 104 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 105
 106 Foam::arcEdge::arcEdge
 107 (
 108     const pointField& points,
 109     const label start,
 110     const label end,
 111     const point& pMid
 112 )
 113 :
 114     curvedEdge(points, start, end),
 115     p1_(points_[start_]),
 116     p2_(pMid),
 117     p3_(points_[end_]),
 118     cs_(calcAngle())
 119 {}
 120
 121
 122 Foam::arcEdge::arcEdge(const pointField& points, Istream& is)
 123 :
 124     curvedEdge(points, is),
 125     p1_(points_[start_]),
 126     p2_(is),
 127     p3_(points_[end_]),
 128     cs_(calcAngle())
 129 {}
 130
 131
 132 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 133
 134 Foam::point Foam::arcEdge::position(const scalar lambda) const
 135 {
 136     if (lambda < 0 || lambda > 1)
 137     {
 138         FatalErrorIn("arcEdge::position(const scalar lambda) const")
 139             << "Parameter out of range, lambda = " << lambda
 140             << abort(FatalError);
 141     }
 142
 143     if (lambda < SMALL)
 144     {
 145         return p1_;
 146     }
 147     else if (lambda > 1 - SMALL)
 148     {
 149         return p3_;
 150     }
 151     else
 152     {
 153         return cs_.globalPosition(vector(radius_, lambda*angle_, 0.0));
 154     }
 155 }
 156
 157
 158 Foam::scalar Foam::arcEdge::length() const
 159 {
 160     return degToRad(angle_*radius_);
 161 }
 162
 163
 164 // ************************************************************************* //