src/OpenFOAM/meshes/primitiveShapes/plane/plane.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 "plane.H"
  27 #include "tensor.H"
  28
  29 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  30
  31 // Calculate base point and unit normal vector from plane equation
  32 void Foam::plane::calcPntAndVec(const scalarList& C)
  33 {
  34     if (mag(C[0]) > VSMALL)
  35     {
  36         basePoint_ = vector((-C[3]/C[0]), 0, 0);
  37     }
  38     else
  39     {
  40         if (mag(C[1]) > VSMALL)
  41         {
  42             basePoint_ = vector(0, (-C[3]/C[1]), 0);
  43         }
  44         else
  45         {
  46             if (mag(C[2]) > VSMALL)
  47             {
  48                 basePoint_ = vector(0, 0, (-C[3]/C[2]));
  49             }
  50             else
  51             {
  52                 FatalErrorIn("void plane::calcPntAndVec(const scalarList&)")
  53                     << "At least one plane coefficient must have a value"
  54                     << abort(FatalError);
  55             }
  56         }
  57     }
  58
  59     unitVector_ = vector(C[0], C[1], C[2]);
  60     scalar magUnitVector(mag(unitVector_));
  61
  62     if (magUnitVector < VSMALL)
  63     {
  64         FatalErrorIn("void plane::calcPntAndVec(const scalarList&)")
  65             << "Plane normal defined with zero length"
  66             << abort(FatalError);
  67     }
  68
  69     unitVector_ /= magUnitVector;
  70 }
  71
  72
  73 void Foam::plane::calcPntAndVec
  74 (
  75     const point& point1,
  76     const point& point2,
  77     const point& point3
  78 )
  79 {
  80     basePoint_ = (point1 + point2 + point3)/3;
  81     vector line12 = point1 - point2;
  82     vector line23 = point2 - point3;
  83
  84     if
  85     (
  86         mag(line12) < VSMALL
  87      || mag(line23) < VSMALL
  88      || mag(point3-point1) < VSMALL
  89     )
  90     {
  91         FatalErrorIn
  92         (
  93             "void plane::calcPntAndVec\n"
  94             "(\n"
  95             "    const point&,\n"
  96             "    const point&,\n"
  97             "    const point&\n"
  98             ")\n"
  99         )   << "Bad points:" << point1 << ' ' << point2 << ' ' << point3
 100             << abort(FatalError);
 101     }
 102
 103     unitVector_ = line12 ^ line23;
 104     scalar magUnitVector(mag(unitVector_));
 105
 106     if (magUnitVector < VSMALL)
 107     {
 108         FatalErrorIn
 109         (
 110             "void plane::calcPntAndVec\n"
 111             "(\n"
 112             "    const point&,\n"
 113             "    const point&,\n"
 114             "    const point&\n"
 115             ")\n"
 116         )   << "Plane normal defined with zero length" << nl
 117             << "Bad points:" << point1 << ' ' << point2 << ' ' << point3
 118             << abort(FatalError);
 119     }
 120
 121     unitVector_ /= magUnitVector;
 122 }
 123
 124
 125 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 126
 127 // Construct from normal vector through the origin
 128 Foam::plane::plane(const vector& normalVector)
 129 :
 130     unitVector_(normalVector),
 131     basePoint_(vector::zero)
 132 {
 133     scalar magUnitVector(mag(unitVector_));
 134
 135     if (magUnitVector > VSMALL)
 136     {
 137         unitVector_ /= magUnitVector;
 138     }
 139     else
 140     {
 141         FatalErrorIn("plane::plane(const vector&)")
 142             << "plane normal has zero length. basePoint:" << basePoint_
 143             << abort(FatalError);
 144     }
 145 }
 146
 147
 148 // Construct from point and normal vector
 149 Foam::plane::plane(const point& basePoint, const vector& normalVector)
 150 :
 151     unitVector_(normalVector),
 152     basePoint_(basePoint)
 153 {
 154     scalar magUnitVector(mag(unitVector_));
 155
 156     if (magUnitVector > VSMALL)
 157     {
 158         unitVector_ /= magUnitVector;
 159     }
 160     else
 161     {
 162         FatalErrorIn("plane::plane(const point&, const vector&)")
 163             << "plane normal has zero length. basePoint:" << basePoint_
 164             << abort(FatalError);
 165     }
 166 }
 167
 168
 169 // Construct from plane equation
 170 Foam::plane::plane(const scalarList& C)
 171 {
 172     calcPntAndVec(C);
 173 }
 174
 175
 176 // Construct from three points
 177 Foam::plane::plane
 178 (
 179     const point& a,
 180     const point& b,
 181     const point& c
 182 )
 183 {
 184     calcPntAndVec(a, b, c);
 185 }
 186
 187
 188 // Construct from dictionary
 189 Foam::plane::plane(const dictionary& dict)
 190 :
 191     unitVector_(vector::zero),
 192     basePoint_(point::zero)
 193 {
 194     const word planeType(dict.lookup("planeType"));
 195
 196     if (planeType == "planeEquation")
 197     {
 198         const dictionary& subDict = dict.subDict("planeEquationDict");
 199         scalarList C(4);
 200
 201         C[0] = readScalar(subDict.lookup("a"));
 202         C[1] = readScalar(subDict.lookup("b"));
 203         C[2] = readScalar(subDict.lookup("c"));
 204         C[3] = readScalar(subDict.lookup("d"));
 205
 206         calcPntAndVec(C);
 207
 208     }
 209     else if (planeType == "embeddedPoints")
 210     {
 211         const dictionary& subDict = dict.subDict("embeddedPoints");
 212
 213         point point1(subDict.lookup("point1"));
 214         point point2(subDict.lookup("point2"));
 215         point point3(subDict.lookup("point3"));
 216
 217         calcPntAndVec(point1, point2, point3);
 218     }
 219     else if (planeType == "pointAndNormal")
 220     {
 221         const dictionary& subDict = dict.subDict("pointAndNormalDict");
 222
 223         basePoint_ = subDict.lookup("basePoint");
 224         unitVector_ = subDict.lookup("normalVector");
 225         unitVector_ /= mag(unitVector_);
 226     }
 227     else
 228     {
 229         FatalIOErrorIn
 230         (
 231             "plane::plane(const dictionary&)",
 232             dict
 233         )   << "Invalid plane type: " << planeType
 234             << abort(FatalIOError);
 235     }
 236 }
 237
 238
 239 // Construct from Istream. Assumes point and normal vector.
 240 Foam::plane::plane(Istream& is)
 241 :
 242     unitVector_(is),
 243     basePoint_(is)
 244 {
 245     scalar magUnitVector(mag(unitVector_));
 246
 247     if (magUnitVector > VSMALL)
 248     {
 249         unitVector_ /= magUnitVector;
 250     }
 251     else
 252     {
 253         FatalErrorIn("plane::plane(Istream& is)")
 254             << "plane normal has zero length. basePoint:" << basePoint_
 255             << abort(FatalError);
 256     }
 257 }
 258
 259
 260 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 261
 262 // Return plane normal vector
 263 const Foam::vector& Foam::plane::normal() const
 264 {
 265     return unitVector_;
 266 }
 267
 268
 269 // Return plane base point
 270 const Foam::point& Foam::plane::refPoint() const
 271 {
 272     return basePoint_;
 273 }
 274
 275
 276 // Return coefficients for plane equation: ax + by + cz + d = 0
 277 Foam::FixedList<Foam::scalar, 4> Foam::plane::planeCoeffs() const
 278 {
 279     FixedList<scalar, 4> C(4);
 280
 281     scalar magX = mag(unitVector_.x());
 282     scalar magY = mag(unitVector_.y());
 283     scalar magZ = mag(unitVector_.z());
 284
 285     if (magX > magY)
 286     {
 287         if (magX > magZ)
 288         {
 289             C[0] = 1;
 290             C[1] = unitVector_.y()/unitVector_.x();
 291             C[2] = unitVector_.z()/unitVector_.x();
 292         }
 293         else
 294         {
 295             C[0] = unitVector_.x()/unitVector_.z();
 296             C[1] = unitVector_.y()/unitVector_.z();
 297             C[2] = 1;
 298         }
 299     }
 300     else
 301     {
 302         if (magY > magZ)
 303         {
 304             C[0] = unitVector_.x()/unitVector_.y();
 305             C[1] = 1;
 306             C[2] = unitVector_.z()/unitVector_.y();
 307         }
 308         else
 309         {
 310             C[0] = unitVector_.x()/unitVector_.z();
 311             C[1] = unitVector_.y()/unitVector_.z();
 312             C[2] = 1;
 313         }
 314     }
 315
 316     C[3] = - C[0] * basePoint_.x()
 317            - C[1] * basePoint_.y()
 318            - C[2] * basePoint_.z();
 319
 320     return C;
 321 }
 322
 323
 324 // Return nearest point in the plane for the given point
 325 Foam::point Foam::plane::nearestPoint(const point& p) const
 326 {
 327     return p - unitVector_*((p - basePoint_) & unitVector_);
 328 }
 329
 330
 331 // Return distance from the given point to the plane
 332 Foam::scalar Foam::plane::distance(const point& p) const
 333 {
 334     return mag((p - basePoint_) & unitVector_);
 335 }
 336
 337
 338 // Cutting point for plane and line defined by origin and direction
 339 Foam::scalar Foam::plane::normalIntersect
 340 (
 341     const point& pnt0,
 342     const vector& dir
 343 ) const
 344 {
 345     scalar denom = stabilise((dir & unitVector_), VSMALL);
 346
 347     return ((basePoint_ - pnt0) & unitVector_)/denom;
 348 }
 349
 350
 351 // Cutting line of two planes
 352 Foam::plane::ray Foam::plane::planeIntersect(const plane& plane2) const
 353 {
 354     // Mathworld plane-plane intersection. Assume there is a point on the
 355     // intersection line with z=0 and solve the two plane equations
 356     // for that (now 2x2 equation in x and y)
 357     // Better: use either z=0 or x=0 or y=0.
 358
 359     const vector& n1 = normal();
 360     const vector& n2 = plane2.normal();
 361
 362     const point& p1 = refPoint();
 363     const point& p2 = plane2.refPoint();
 364
 365     scalar n1p1 = n1&p1;
 366     scalar n2p2 = n2&p2;
 367
 368     vector dir = n1 ^ n2;
 369
 370     // Determine zeroed out direction (can be x,y or z) by looking at which
 371     // has the largest component in dir.
 372     scalar magX = mag(dir.x());
 373     scalar magY = mag(dir.y());
 374     scalar magZ = mag(dir.z());
 375
 376     direction iZero, i1, i2;
 377
 378     if (magX > magY)
 379     {
 380         if (magX > magZ)
 381         {
 382             iZero = 0;
 383             i1 = 1;
 384             i2 = 2;
 385         }
 386         else
 387         {
 388             iZero = 2;
 389             i1 = 0;
 390             i2 = 1;
 391         }
 392     }
 393     else
 394     {
 395         if (magY > magZ)
 396         {
 397             iZero = 1;
 398             i1 = 2;
 399             i2 = 0;
 400         }
 401         else
 402         {
 403             iZero = 2;
 404             i1 = 0;
 405             i2 = 1;
 406         }
 407     }
 408
 409     vector pt;
 410
 411     pt[iZero] = 0;
 412     pt[i1] = (n2[i2]*n1p1 - n1[i2]*n2p2) / (n1[i1]*n2[i2] - n2[i1]*n1[i2]);
 413     pt[i2] = (n2[i1]*n1p1 - n1[i1]*n2p2) / (n1[i2]*n2[i1] - n1[i1]*n2[i2]);
 414
 415     return ray(pt, dir);
 416 }
 417
 418
 419 // Cutting point of three planes
 420 Foam::point Foam::plane::planePlaneIntersect
 421 (
 422     const plane& plane2,
 423     const plane& plane3
 424 ) const
 425 {
 426     FixedList<scalar, 4> coeffs1(planeCoeffs());
 427     FixedList<scalar, 4> coeffs2(plane2.planeCoeffs());
 428     FixedList<scalar, 4> coeffs3(plane3.planeCoeffs());
 429
 430     tensor a
 431     (
 432         coeffs1[0],coeffs1[1],coeffs1[2],
 433         coeffs2[0],coeffs2[1],coeffs2[2],
 434         coeffs3[0],coeffs3[1],coeffs3[2]
 435     );
 436
 437     vector b(coeffs1[3],coeffs2[3],coeffs3[3]);
 438
 439     return (inv(a) & (-b));
 440 }
 441
 442
 443 void Foam::plane::writeDict(Ostream& os) const
 444 {
 445     os.writeKeyword("planeType") << "pointAndNormal"
 446         << token::END_STATEMENT << nl;
 447     os  << indent << "pointAndNormalDict" << nl
 448         << indent << token::BEGIN_BLOCK << incrIndent << nl;
 449     os.writeKeyword("basePoint") << basePoint_ << token::END_STATEMENT << nl;
 450     os.writeKeyword("normalVector") << unitVector_ << token::END_STATEMENT
 451         << nl;
 452     os << decrIndent << indent << token::END_BLOCK << endl;
 453 }
 454
 455
 456 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 457
 458 bool Foam::operator==(const plane& a, const plane& b)
 459 {
 460     if (a.basePoint_ == b.basePoint_ && a.unitVector_ == b.unitVector_)
 461     {
 462         return true;
 463     }
 464     else
 465     {
 466         return false;
 467     }
 468 }
 469
 470 bool Foam::operator!=(const plane& a, const plane& b)
 471 {
 472     return !(a == b);
 473 }
 474
 475
 476 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 477
 478 Foam::Ostream& Foam::operator<<(Ostream& os, const plane& a)
 479 {
 480     os  << a.unitVector_ << token::SPACE << a.basePoint_;
 481
 482     return os;
 483 }
 484
 485
 486 // ************************************************************************* //