FreeFOAM The Cross-Platform CFD Toolkit
arcEdge.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
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"
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_ = (acos((r3 & r1)/(mag(r3) * mag(r1))))
75 
76  // check if the vectors define an exterior or an interior arcEdge
77  if (((r1 ^ r2) & (r1 ^ r3)) < 0.0)
78  {
79  angle_ = 360.0 - angle_;
80  }
81 
82  vector tempAxis;
83 
84  if (angle_ <= 180.0)
85  {
86  tempAxis = r1 ^ r3;
87 
88  if (mag(tempAxis)/(mag(r1)*mag(r3)) < 0.001)
89  {
90  tempAxis = r1 ^ r2;
91  }
92  }
93  else
94  {
95  tempAxis = r3 ^ r1;
96  }
97 
98  radius_ = mag(r3);
99 
100  // set up and return the local coordinate system
101  return cylindricalCS("arcEdgeCS", centre, tempAxis, r1);
102 }
103 
104 
105 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
106 
107 Foam::arcEdge::arcEdge
108 (
109  const pointField& points,
110  const label start,
111  const label end,
112  const point& pMid
113 )
114 :
115  curvedEdge(points, start, end),
116  p1_(points_[start_]),
117  p2_(pMid),
118  p3_(points_[end_]),
119  cs_(calcAngle())
120 {}
121 
122 
123 Foam::arcEdge::arcEdge(const pointField& points, Istream& is)
124 :
125  curvedEdge(points, is),
126  p1_(points_[start_]),
127  p2_(is),
128  p3_(points_[end_]),
129  cs_(calcAngle())
130 {}
131 
132 
133 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
134 
135 Foam::point Foam::arcEdge::position(const scalar lambda) const
136 {
137  if (lambda < 0 || lambda > 1)
138  {
139  FatalErrorIn("arcEdge::position(const scalar lambda) const")
140  << "Parameter out of range, lambda = " << lambda
141  << abort(FatalError);
142  }
143 
144  if (lambda < SMALL)
145  {
146  return p1_;
147  }
148  else if (lambda > 1 - SMALL)
149  {
150  return p3_;
151  }
152  else
153  {
154  return cs_.globalPosition(vector(radius_, lambda*angle_, 0.0));
155  }
156 }
157 
158 
159 Foam::scalar Foam::arcEdge::length() const
160 {
161  return (angle_ * radius_) * mathematicalConstant::pi/180.0;
162 }
163 
164 
165 // ************************ vim: set sw=4 sts=4 et: ************************ //