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