src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.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 "treeBoundBox.H"
  27 #include "ListOps.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 const Foam::scalar Foam::treeBoundBox::great(GREAT);
  32
  33 const Foam::treeBoundBox Foam::treeBoundBox::greatBox
  34 (
  35     vector(-GREAT, -GREAT, -GREAT),
  36     vector(GREAT, GREAT, GREAT)
  37 );
  38
  39
  40 const Foam::treeBoundBox Foam::treeBoundBox::invertedBox
  41 (
  42     vector(GREAT, GREAT, GREAT),
  43     vector(-GREAT, -GREAT, -GREAT)
  44 );
  45
  46
  47 //! \cond - skip documentation : local scope only
  48 const Foam::label facesArray[6][4] =
  49 {
  50     {0, 4, 6, 2}, // left
  51     {1, 3, 7, 5}, // right
  52     {0, 1, 5, 4}, // bottom
  53     {2, 6, 7, 3}, // top
  54     {0, 2, 3, 1}, // back
  55     {4, 5, 7, 6}  // front
  56 };
  57 //! \endcond
  58
  59
  60 const Foam::faceList Foam::treeBoundBox::faces
  61 (
  62     initListList<face, label, 6, 4>(facesArray)
  63 );
  64
  65
  66 //! \cond - skip documentation : local scope only
  67 const Foam::label edgesArray[12][2] =
  68 {
  69     {0, 1}, // 0
  70     {1, 3},
  71     {2, 3}, // 2
  72     {0, 2},
  73     {4, 5}, // 4
  74     {5, 7},
  75     {6, 7}, // 6
  76     {4, 6},
  77     {0, 4}, // 8
  78     {1, 5},
  79     {3, 7}, // 10
  80     {2, 6}
  81 };
  82 //! \endcond
  83
  84
  85 const Foam::edgeList Foam::treeBoundBox::edges
  86 (
  87     //initListList<edge, label, 12, 2>(edgesArray)
  88     calcEdges(edgesArray)
  89 );
  90
  91
  92 const Foam::FixedList<Foam::vector, 6> Foam::treeBoundBox::faceNormals
  93 (
  94     calcFaceNormals()
  95 );
  96
  97
  98 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  99
 100 Foam::edgeList Foam::treeBoundBox::calcEdges(const label edgesArray[12][2])
 101 {
 102     edgeList edges(12);
 103     forAll(edges, edgeI)
 104     {
 105         edges[edgeI][0] = edgesArray[edgeI][0];
 106         edges[edgeI][1] = edgesArray[edgeI][1];
 107     }
 108     return edges;
 109 }
 110
 111
 112 Foam::FixedList<Foam::vector, 6> Foam::treeBoundBox::calcFaceNormals()
 113 {
 114     FixedList<vector, 6> normals;
 115     normals[LEFT]   = vector(-1,  0,  0);
 116     normals[RIGHT]  = vector( 1,  0,  0);
 117     normals[BOTTOM] = vector( 0, -1,  0);
 118     normals[TOP]    = vector( 0,  1,  0);
 119     normals[BACK]   = vector( 0,  0, -1);
 120     normals[FRONT]  = vector( 0,  0,  1);
 121     return normals;
 122 }
 123
 124
 125 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 126
 127 Foam::treeBoundBox::treeBoundBox(const UList<point>& points)
 128 :
 129     boundBox(points, false)
 130 {
 131     if (points.empty())
 132     {
 133         WarningIn
 134         (
 135             "treeBoundBox::treeBoundBox(const UList<point>&)"
 136         )   << "cannot find bounding box for zero-sized pointField, "
 137             << "returning zero" << endl;
 138
 139         return;
 140     }
 141 }
 142
 143
 144 Foam::treeBoundBox::treeBoundBox
 145 (
 146     const UList<point>& points,
 147     const labelUList& indices
 148 )
 149 :
 150     boundBox(points, indices, false)
 151 {
 152     if (points.empty() || indices.empty())
 153     {
 154         WarningIn
 155         (
 156             "treeBoundBox::treeBoundBox"
 157             "(const UList<point>&, const labelUList&)"
 158         )   << "cannot find bounding box for zero-sized pointField, "
 159             << "returning zero" << endl;
 160
 161         return;
 162     }
 163 }
 164
 165
 166 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 167
 168 Foam::tmp<Foam::pointField> Foam::treeBoundBox::points() const
 169 {
 170     tmp<pointField> tPts = tmp<pointField>(new pointField(8));
 171
 172     pointField& points = tPts();
 173
 174     forAll(points, octant)
 175     {
 176         points[octant] = corner(octant);
 177
 178     }
 179
 180     return tPts;
 181 }
 182
 183
 184 Foam::treeBoundBox Foam::treeBoundBox::subBbox(const direction octant) const
 185 {
 186     return subBbox(midpoint(), octant);
 187 }
 188
 189
 190 // Octant to bounding box using permutation only.
 191 Foam::treeBoundBox Foam::treeBoundBox::subBbox
 192 (
 193     const point& mid,
 194     const direction octant
 195 ) const
 196 {
 197     if (octant > 7)
 198     {
 199         FatalErrorIn
 200         (
 201             "treeBoundBox::subBbox(const point&, const direction)"
 202         )   << "octant should be [0..7]"
 203             << abort(FatalError);
 204     }
 205
 206     // start with a copy of this bounding box and adjust limits accordingly
 207     treeBoundBox subBb(*this);
 208     point& bbMin = subBb.min();
 209     point& bbMax = subBb.max();
 210
 211     if (octant & treeBoundBox::RIGHTHALF)
 212     {
 213         bbMin.x() = mid.x();    // mid -> max
 214     }
 215     else
 216     {
 217         bbMax.x() = mid.x();    // min -> mid
 218     }
 219
 220     if (octant & treeBoundBox::TOPHALF)
 221     {
 222         bbMin.y() = mid.y();    // mid -> max
 223     }
 224     else
 225     {
 226         bbMax.y() = mid.y();    // min -> mid
 227     }
 228
 229     if (octant & treeBoundBox::FRONTHALF)
 230     {
 231         bbMin.z() = mid.z();    // mid -> max
 232     }
 233     else
 234     {
 235         bbMax.z() = mid.z();    // min -> mid
 236     }
 237
 238     return subBb;
 239 }
 240
 241
 242 bool Foam::treeBoundBox::overlaps
 243 (
 244     const point& centre,
 245     const scalar radiusSqr
 246 ) const
 247 {
 248     // Find out where centre is in relation to bb.
 249     // Find nearest point on bb.
 250     scalar distSqr = 0;
 251
 252     for (direction dir = 0; dir < vector::nComponents; dir++)
 253     {
 254         scalar d0 = min()[dir] - centre[dir];
 255         scalar d1 = max()[dir] - centre[dir];
 256
 257         if ((d0 > 0) != (d1 > 0))
 258         {
 259             // centre inside both extrema. This component does not add any
 260             // distance.
 261         }
 262         else if (Foam::mag(d0) < Foam::mag(d1))
 263         {
 264             distSqr += d0*d0;
 265         }
 266         else
 267         {
 268             distSqr += d1*d1;
 269         }
 270
 271         if (distSqr > radiusSqr)
 272         {
 273             return false;
 274         }
 275     }
 276
 277     return true;
 278 }
 279
 280
 281 // line intersection. Returns true if line (start to end) inside
 282 // bb or intersects bb. Sets pt to intersection.
 283 //
 284 // Sutherlands algorithm:
 285 //   loop
 286 //     - start = intersection of line with one of the planes bounding
 287 //       the bounding box
 288 //     - stop if start inside bb (return true)
 289 //     - stop if start and end in same 'half' (e.g. both above bb)
 290 //       (return false)
 291 //
 292 // Uses posBits to efficiently determine 'half' in which start and end
 293 // point are.
 294 //
 295 // Note:
 296 //   - sets coordinate to exact position: e.g. pt.x() = min().x()
 297 //     since plane intersect routine might have truncation error.
 298 //     This makes sure that posBits tests 'inside'
 299 bool Foam::treeBoundBox::intersects
 300 (
 301     const point& overallStart,
 302     const vector& overallVec,
 303     const point& start,
 304     const point& end,
 305     point& pt,
 306     direction& ptOnFaces
 307 ) const
 308 {
 309     const direction endBits = posBits(end);
 310     pt = start;
 311
 312     while (true)
 313     {
 314         direction ptBits = posBits(pt);
 315
 316         if (ptBits == 0)
 317         {
 318             // pt inside bb
 319             ptOnFaces = faceBits(pt);
 320             return true;
 321         }
 322
 323         if ((ptBits & endBits) != 0)
 324         {
 325             // pt and end in same block outside of bb
 326             ptOnFaces = faceBits(pt);
 327             return false;
 328         }
 329
 330         if (ptBits & LEFTBIT)
 331         {
 332             // Intersect with plane V=min, n=-1,0,0
 333             if (Foam::mag(overallVec.x()) > VSMALL)
 334             {
 335                 scalar s = (min().x() - overallStart.x())/overallVec.x();
 336                 pt.x() = min().x();
 337                 pt.y() = overallStart.y() + overallVec.y()*s;
 338                 pt.z() = overallStart.z() + overallVec.z()*s;
 339             }
 340             else
 341             {
 342                 // Vector not in x-direction. But still intersecting bb planes.
 343                 // So must be close - just snap to bb.
 344                 pt.x() = min().x();
 345             }
 346         }
 347         else if (ptBits & RIGHTBIT)
 348         {
 349             // Intersect with plane V=max, n=1,0,0
 350             if (Foam::mag(overallVec.x()) > VSMALL)
 351             {
 352                 scalar s = (max().x() - overallStart.x())/overallVec.x();
 353                 pt.x() = max().x();
 354                 pt.y() = overallStart.y() + overallVec.y()*s;
 355                 pt.z() = overallStart.z() + overallVec.z()*s;
 356             }
 357             else
 358             {
 359                 pt.x() = max().x();
 360             }
 361         }
 362         else if (ptBits & BOTTOMBIT)
 363         {
 364             // Intersect with plane V=min, n=0,-1,0
 365             if (Foam::mag(overallVec.y()) > VSMALL)
 366             {
 367                 scalar s = (min().y() - overallStart.y())/overallVec.y();
 368                 pt.x() = overallStart.x() + overallVec.x()*s;
 369                 pt.y() = min().y();
 370                 pt.z() = overallStart.z() + overallVec.z()*s;
 371             }
 372             else
 373             {
 374                 pt.x() = min().y();
 375             }
 376         }
 377         else if (ptBits & TOPBIT)
 378         {
 379             // Intersect with plane V=max, n=0,1,0
 380             if (Foam::mag(overallVec.y()) > VSMALL)
 381             {
 382                 scalar s = (max().y() - overallStart.y())/overallVec.y();
 383                 pt.x() = overallStart.x() + overallVec.x()*s;
 384                 pt.y() = max().y();
 385                 pt.z() = overallStart.z() + overallVec.z()*s;
 386             }
 387             else
 388             {
 389                 pt.y() = max().y();
 390             }
 391         }
 392         else if (ptBits & BACKBIT)
 393         {
 394             // Intersect with plane V=min, n=0,0,-1
 395             if (Foam::mag(overallVec.z()) > VSMALL)
 396             {
 397                 scalar s = (min().z() - overallStart.z())/overallVec.z();
 398                 pt.x() = overallStart.x() + overallVec.x()*s;
 399                 pt.y() = overallStart.y() + overallVec.y()*s;
 400                 pt.z() = min().z();
 401             }
 402             else
 403             {
 404                 pt.z() = min().z();
 405             }
 406         }
 407         else if (ptBits & FRONTBIT)
 408         {
 409             // Intersect with plane V=max, n=0,0,1
 410             if (Foam::mag(overallVec.z()) > VSMALL)
 411             {
 412                 scalar s = (max().z() - overallStart.z())/overallVec.z();
 413                 pt.x() = overallStart.x() + overallVec.x()*s;
 414                 pt.y() = overallStart.y() + overallVec.y()*s;
 415                 pt.z() = max().z();
 416             }
 417             else
 418             {
 419                 pt.z() = max().z();
 420             }
 421         }
 422     }
 423 }
 424
 425
 426 bool Foam::treeBoundBox::intersects
 427 (
 428     const point& start,
 429     const point& end,
 430     point& pt
 431 ) const
 432 {
 433     direction ptBits;
 434     return intersects(start, end-start, start, end, pt, ptBits);
 435 }
 436
 437
 438 bool Foam::treeBoundBox::contains(const vector& dir, const point& pt) const
 439 {
 440     //
 441     // Compare all components against min and max of bb
 442     //
 443     for (direction cmpt=0; cmpt<3; cmpt++)
 444     {
 445         if (pt[cmpt] < min()[cmpt])
 446         {
 447             return false;
 448         }
 449         else if (pt[cmpt] == min()[cmpt])
 450         {
 451             // On edge. Outside if direction points outwards.
 452             if (dir[cmpt] < 0)
 453             {
 454                 return false;
 455             }
 456         }
 457
 458         if (pt[cmpt] > max()[cmpt])
 459         {
 460             return false;
 461         }
 462         else if (pt[cmpt] == max()[cmpt])
 463         {
 464             // On edge. Outside if direction points outwards.
 465             if (dir[cmpt] > 0)
 466             {
 467                 return false;
 468             }
 469         }
 470     }
 471
 472     // All components inside bb
 473     return true;
 474 }
 475
 476
 477 // Code position of pt on bounding box faces
 478 Foam::direction Foam::treeBoundBox::faceBits(const point& pt) const
 479 {
 480     direction faceBits = 0;
 481     if (pt.x() == min().x())
 482     {
 483         faceBits |= LEFTBIT;
 484     }
 485     else if (pt.x() == max().x())
 486     {
 487         faceBits |= RIGHTBIT;
 488     }
 489
 490     if (pt.y() == min().y())
 491     {
 492         faceBits |= BOTTOMBIT;
 493     }
 494     else if (pt.y() == max().y())
 495     {
 496         faceBits |= TOPBIT;
 497     }
 498
 499     if (pt.z() == min().z())
 500     {
 501         faceBits |= BACKBIT;
 502     }
 503     else if (pt.z() == max().z())
 504     {
 505         faceBits |= FRONTBIT;
 506     }
 507     return faceBits;
 508 }
 509
 510
 511 // Code position of point relative to box
 512 Foam::direction Foam::treeBoundBox::posBits(const point& pt) const
 513 {
 514     direction posBits = 0;
 515
 516     if (pt.x() < min().x())
 517     {
 518         posBits |= LEFTBIT;
 519     }
 520     else if (pt.x() > max().x())
 521     {
 522         posBits |= RIGHTBIT;
 523     }
 524
 525     if (pt.y() < min().y())
 526     {
 527         posBits |= BOTTOMBIT;
 528     }
 529     else if (pt.y() > max().y())
 530     {
 531         posBits |= TOPBIT;
 532     }
 533
 534     if (pt.z() < min().z())
 535     {
 536         posBits |= BACKBIT;
 537     }
 538     else if (pt.z() > max().z())
 539     {
 540         posBits |= FRONTBIT;
 541     }
 542     return posBits;
 543 }
 544
 545
 546 // nearest and furthest corner coordinate.
 547 // !names of treeBoundBox::min() and treeBoundBox::max() are confusing!
 548 void Foam::treeBoundBox::calcExtremities
 549 (
 550     const point& pt,
 551     point& nearest,
 552     point& furthest
 553 ) const
 554 {
 555     scalar nearX, nearY, nearZ;
 556     scalar farX, farY, farZ;
 557
 558     if (Foam::mag(min().x() - pt.x()) < Foam::mag(max().x() - pt.x()))
 559     {
 560         nearX = min().x();
 561         farX = max().x();
 562     }
 563     else
 564     {
 565         nearX = max().x();
 566         farX = min().x();
 567     }
 568
 569     if (Foam::mag(min().y() - pt.y()) < Foam::mag(max().y() - pt.y()))
 570     {
 571         nearY = min().y();
 572         farY = max().y();
 573     }
 574     else
 575     {
 576         nearY = max().y();
 577         farY = min().y();
 578     }
 579
 580     if (Foam::mag(min().z() - pt.z()) < Foam::mag(max().z() - pt.z()))
 581     {
 582         nearZ = min().z();
 583         farZ = max().z();
 584     }
 585     else
 586     {
 587         nearZ = max().z();
 588         farZ = min().z();
 589     }
 590
 591     nearest = point(nearX, nearY, nearZ);
 592     furthest = point(farX, farY, farZ);
 593 }
 594
 595
 596 Foam::scalar Foam::treeBoundBox::maxDist(const point& pt) const
 597 {
 598     point near, far;
 599     calcExtremities(pt, near, far);
 600
 601     return Foam::mag(far - pt);
 602 }
 603
 604
 605 // Distance comparator
 606 // Compare all vertices of bounding box against all of other bounding
 607 // box to see if all vertices of one are nearer
 608 Foam::label Foam::treeBoundBox::distanceCmp
 609 (
 610     const point& pt,
 611     const treeBoundBox& other
 612 ) const
 613 {
 614     //
 615     // Distance point <-> nearest and furthest away vertex of this
 616     //
 617
 618     point nearThis, farThis;
 619
 620     // get nearest and furthest away vertex
 621     calcExtremities(pt, nearThis, farThis);
 622
 623     const scalar minDistThis =
 624         sqr(nearThis.x() - pt.x())
 625      +  sqr(nearThis.y() - pt.y())
 626      +  sqr(nearThis.z() - pt.z());
 627     const scalar maxDistThis =
 628         sqr(farThis.x() - pt.x())
 629      +  sqr(farThis.y() - pt.y())
 630      +  sqr(farThis.z() - pt.z());
 631
 632     //
 633     // Distance point <-> other
 634     //
 635
 636     point nearOther, farOther;
 637
 638     // get nearest and furthest away vertex
 639     other.calcExtremities(pt, nearOther, farOther);
 640
 641     const scalar minDistOther =
 642         sqr(nearOther.x() - pt.x())
 643      +  sqr(nearOther.y() - pt.y())
 644      +  sqr(nearOther.z() - pt.z());
 645     const scalar maxDistOther =
 646         sqr(farOther.x() - pt.x())
 647      +  sqr(farOther.y() - pt.y())
 648      +  sqr(farOther.z() - pt.z());
 649
 650     //
 651     // Categorize
 652     //
 653     if (maxDistThis < minDistOther)
 654     {
 655         // All vertices of this are nearer to point than any vertex of other
 656         return -1;
 657     }
 658     else if (minDistThis > maxDistOther)
 659     {
 660         // All vertices of this are further from point than any vertex of other
 661         return 1;
 662     }
 663     else
 664     {
 665         // Mixed bag
 666         return 0;
 667     }
 668 }
 669
 670
 671 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 672
 673 bool Foam::operator==(const treeBoundBox& a, const treeBoundBox& b)
 674 {
 675     return operator==
 676     (
 677         static_cast<const boundBox&>(a),
 678         static_cast<const boundBox&>(b)
 679     );
 680 }
 681
 682
 683 bool Foam::operator!=(const treeBoundBox& a, const treeBoundBox& b)
 684 {
 685     return !(a == b);
 686 }
 687
 688
 689 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 690
 691 Foam::Ostream& Foam::operator<<(Ostream& os, const treeBoundBox& bb)
 692 {
 693     return os << static_cast<const boundBox&>(bb);
 694 }
 695
 696
 697 Foam::Istream& Foam::operator>>(Istream& is, treeBoundBox& bb)
 698 {
 699     return is >> static_cast<boundBox&>(bb);
 700 }
 701
 702
 703 // ************************************************************************* //