src/mesh/cfMesh/meshLibrary/utilities/octrees/meshOctree/meshOctreeCube/meshOctreeCubeCoordinatesIntersections.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     cfMesh 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 cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "triSurf.H"
  29 #include "meshOctreeCubeCoordinates.H"
  30 #include "helperFunctions.H"
  31
  32 //#define DEBUGSearch
  33
  34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40 // Static data
  41
  42 const label meshOctreeCubeCoordinates::edgeNodes_[12][2] =
  43     {
  44         //- edges in x-direction
  45         {0, 1},
  46         {2, 3},
  47         {4, 5},
  48         {6, 7},
  49         //- edges in y-direction
  50         {0, 2},
  51         {1, 3},
  52         {4, 6},
  53         {5, 7},
  54         //- edges in z-direction
  55         {0, 4},
  56         {1, 5},
  57         {2, 6},
  58         {3, 7}
  59     };
  60
  61 const label meshOctreeCubeCoordinates::faceNodes_[6][4] =
  62     {
  63         {0, 4, 6, 2},
  64         {1, 3, 7, 5},
  65         {0, 1, 5, 4},
  66         {2, 6, 7, 3},
  67         {0, 2, 3, 1},
  68         {4, 5, 7, 6}
  69     };
  70
  71 const label meshOctreeCubeCoordinates::nodeFaces_[8][3] =
  72     {
  73         {0, 2, 4},
  74         {1, 2, 4},
  75         {0, 3, 4},
  76         {1, 3, 4},
  77         {0, 2, 5},
  78         {1, 2, 5},
  79         {0, 3, 5},
  80         {1, 3, 5}
  81     };
  82
  83 const label meshOctreeCubeCoordinates::faceEdges_[6][4] =
  84     {
  85         {8, 6, 10, 4},
  86         {5, 11, 7, 9},
  87         {0, 9, 2, 8},
  88         {10, 3, 11, 1},
  89         {4, 1, 5, 0},
  90         {2, 7, 3, 6}
  91     };
  92
  93 const label meshOctreeCubeCoordinates::edgeFaces_[12][2] =
  94     {
  95         {2, 4},
  96         {3, 4},
  97         {2, 5},
  98         {3, 5},
  99         {0, 4},
 100         {1, 4},
 101         {0, 5},
 102         {1, 5},
 103         {0, 2},
 104         {1, 2},
 105         {0, 3},
 106         {1, 3}
 107     };
 108
 109 const label meshOctreeCubeCoordinates::oppositeFace_[6] = {1, 0, 3, 2, 5, 4};
 110
 111 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 112
 113 void meshOctreeCubeCoordinates::vertices
 114 (
 115     const boundBox& rootBox,
 116     FixedList<point, 8>& vrt
 117 ) const
 118 {
 119     const vector tol = SMALL * (rootBox.max() - rootBox.min());
 120
 121     point min_, max_;
 122     cubeBox(rootBox, min_, max_);
 123     min_ -= tol;
 124     max_ += tol;
 125
 126     vrt[0] = point(min_.x(), min_.y(), min_.z());
 127     vrt[1] = point(max_.x(), min_.y(), min_.z());
 128     vrt[2] = point(min_.x(), max_.y(), min_.z());
 129     vrt[3] = point(max_.x(), max_.y(), min_.z());
 130     vrt[4] = point(min_.x(), min_.y(), max_.z());
 131     vrt[5] = point(max_.x(), min_.y(), max_.z());
 132     vrt[6] = point(min_.x(), max_.y(), max_.z());
 133     vrt[7] = point(max_.x(), max_.y(), max_.z());
 134 }
 135
 136 void meshOctreeCubeCoordinates::edgeVertices
 137 (
 138     const boundBox& rootBox,
 139     FixedList<FixedList<point, 2>, 12>& e
 140 ) const
 141 {
 142     FixedList<point, 8> vertices;
 143     this->vertices(rootBox, vertices);
 144
 145     for(label i=0;i<12;++i)
 146     {
 147         e[i][0] = vertices[edgeNodes_[i][0]];
 148         e[i][1] = vertices[edgeNodes_[i][1]];
 149     }
 150 }
 151
 152 bool meshOctreeCubeCoordinates::intersectsTriangle
 153 (
 154     const triSurf& surface,
 155     const boundBox& rootBox,
 156     const label tI
 157 ) const
 158 {
 159     const pointField& points = surface.points();
 160     const labelledTri& ltri = surface[tI];
 161
 162     const vector tol = SMALL * (rootBox.max() - rootBox.min());
 163
 164     //- calculate the bound box of the octree cube
 165     boundBox cBox;
 166     cubeBox(rootBox, cBox.min(), cBox.max());
 167     cBox.min() -= tol;
 168     cBox.max() += tol;
 169
 170     //- calculate the bounding box of the triangle
 171     boundBox tBox;
 172     tBox.min() = tBox.max() = points[ltri[0]];
 173
 174     for(label pI=1;pI<3;++pI)
 175     {
 176         tBox.max() = Foam::max(points[ltri[pI]], tBox.max());
 177         tBox.min() = Foam::min(points[ltri[pI]], tBox.min());
 178     }
 179
 180     tBox.min() -= tol;
 181     tBox.max() += tol;
 182
 183     return cBox.overlaps(tBox);
 184 }
 185
 186 bool meshOctreeCubeCoordinates::intersectsTriangleExact
 187 (
 188     const triSurf& surface,
 189     const boundBox& rootBox,
 190     const label tI
 191 ) const
 192 {
 193     if( !intersectsTriangle(surface, rootBox, tI) )
 194         return false;
 195
 196     const vector tol = SMALL * (rootBox.max() - rootBox.min());
 197
 198     const pointField& points = surface.points();
 199     const labelledTri& ltri = surface[tI];
 200
 201     //- check if any of the vertices is in the cube
 202     forAll(ltri, pI)
 203         if( isVertexInside(rootBox, points[ltri[pI]]) )
 204             return true;
 205
 206     //- check if any edges of the triangle intersect the cube
 207     boundBox bb;
 208     cubeBox(rootBox, bb.min(), bb.max());
 209     bb.min() -= tol;
 210     bb.max() += tol;
 211
 212     for(label eI=0;eI<3;++eI)
 213     {
 214         const edge edg(ltri[eI], ltri[(eI+1)%3]);
 215         const point& s = points[edg.start()];
 216         const point& e = points[edg.end()];
 217
 218         if( help::boundBoxLineIntersection(s, e, bb) )
 219             return true;
 220     }
 221
 222     //- check if any cube edges intersects the triangle
 223     FixedList<FixedList<point, 2>, 12> e;
 224     this->edgeVertices(rootBox, e);
 225
 226     point intersection;
 227     forAll(e, eI)
 228         if(
 229             help::triLineIntersection
 230             (
 231                 surface,
 232                 tI,
 233                 e[eI][0],
 234                 e[eI][1],
 235                 intersection
 236             )
 237         )
 238             return true;
 239
 240     return false;
 241 }
 242
 243 bool meshOctreeCubeCoordinates::isVertexInside
 244 (
 245     const boundBox& rootBox,
 246     const point& p
 247 ) const
 248 {
 249     const vector tol = SMALL * (rootBox.max() - rootBox.min());
 250
 251     point min, max;
 252     cubeBox(rootBox, min, max);
 253     max += tol;
 254     min -= tol;
 255
 256     if(
 257         ((p.x() - max.x()) > 0.0) ||
 258         ((p.y() - max.y()) > 0.0) ||
 259         ((p.z() - max.z()) > 0.0) ||
 260         ((p.x() - min.x()) < 0.0) ||
 261         ((p.y() - min.y()) < 0.0) ||
 262         ((p.z() - min.z()) < 0.0)
 263     )
 264         return false;
 265
 266     return true;
 267 }
 268
 269 bool meshOctreeCubeCoordinates::isPositionInside
 270 (
 271     const meshOctreeCubeCoordinates& cc
 272 ) const
 273 {
 274     # ifdef DEBUGSearch
 275     Info << "Checking cube " << *this << endl;
 276     Info << "level " << short(l) << " x: " << px
 277         << " y: " << py << " z:" << pz << endl;
 278     # endif
 279
 280     if( cc.level() >= this->level() )
 281     {
 282         const direction diff = cc.level() - this->level();
 283         meshOctreeCubeCoordinates reducedLevel =
 284             cc.reduceLevelBy(diff);
 285
 286         # ifdef DEBUGSearch
 287         Info << "diff " << label(diff) << endl;
 288         Info << "Divider " << divider << endl;
 289         Info << "Coordinates at level are " << reducedLevel << endl;
 290         # endif
 291
 292         if( reducedLevel == *this )
 293             return true;
 294     }
 295     else
 296     {
 297         FatalErrorIn
 298         (
 299             "bool meshOctreeCubeCoordinates::isPositionInside"
 300             "(const label px,const label py,"
 301             "const label pz,const direction l)"
 302         ) << "Cannot find exact position of finer cube" << exit(FatalError);
 303     }
 304
 305     return false;
 306 }
 307
 308 bool meshOctreeCubeCoordinates::intersectsLine
 309 (
 310     const boundBox& rootBox,
 311     const point& s,
 312     const point& e
 313 ) const
 314 {
 315     const scalar tol = SMALL * (rootBox.max().x() - rootBox.min().x());
 316
 317     point min, max;
 318     cubeBox(rootBox, min, max);
 319
 320     //- check if the cube contains start point or end point
 321     min -= vector(tol,tol,tol);
 322     max += vector(tol,tol,tol);
 323
 324     //- check for intersections of line with the cube faces
 325     const vector v(e - s);
 326     scalar t;
 327     point i;
 328
 329     if( mag(v.x()) > tol )
 330     {
 331         //- x-min face
 332         t = (min.x() - s.x()) / v.x();
 333         i = s + t * v;
 334         if(
 335             (t > -tol) && (t < (1.0+tol)) &&
 336             (i.y() - min.y() > -tol) &&
 337             (i.y() - max.y() < tol) &&
 338             (i.z() - min.z() > -tol) &&
 339             (i.z() - max.z() < tol)
 340         )
 341             return true;
 342
 343         //- x-max face
 344         t = (max.x() - s.x()) / v.x();
 345         i = s + t * v;
 346         if(
 347             (t > -tol) && (t < (1.0+tol)) &&
 348             (i.y() - min.y() > -tol) &&
 349             (i.y() - max.y() < tol) &&
 350             (i.z() - min.z() > -tol) &&
 351             (i.z() - max.z() < tol)
 352         )
 353             return true;
 354     }
 355
 356     if( mag(v.y()) > tol)
 357     {
 358         //- y-min face
 359         t = (min.y() - s.y()) / v.y();
 360         i = s + t * v;
 361         if(
 362             (t > -tol) && (t < (1.0+tol)) &&
 363             (i.x() - min.x() > -tol) &&
 364             (i.x() - max.x() < tol) &&
 365             (i.z() - min.z() > -tol) &&
 366             (i.z() - max.z() < tol)
 367         )
 368             return true;
 369
 370         //- y-max face
 371         t = (max.y() - s.y()) / v.y();
 372         i = s + t * v;
 373         if(
 374             (t > -tol) && (t < (1.0+tol)) &&
 375             (i.x() - min.x() > -tol) &&
 376             (i.x() - max.x() < tol) &&
 377             (i.z() - min.z() > -tol) &&
 378             (i.z() - max.z() < tol)
 379         )
 380             return true;
 381     }
 382
 383     if( mag(v.z()) > tol )
 384     {
 385         //- z-min face
 386         t = (min.z() - s.z()) / v.z();
 387         i = s + t * v;
 388         if(
 389             (t > -tol) && (t < (1.0+tol)) &&
 390             (i.x() - min.x() > -tol) &&
 391             (i.x() - max.x() < tol) &&
 392             (i.y() - min.y() > -tol) &&
 393             (i.y() - max.y() < tol)
 394         )
 395             return true;
 396
 397         //- z-min face
 398         t = (max.z() - s.z()) / v.z();
 399         i = s + t * v;
 400         if(
 401             (t > -tol) && (t < (1.0+tol)) &&
 402             (i.x() - min.x() > -tol) &&
 403             (i.x() - max.x() < tol) &&
 404             (i.y() - min.y() > -tol) &&
 405             (i.y() - max.y() < tol)
 406         )
 407             return true;
 408     }
 409
 410     if( isVertexInside(rootBox, s) )
 411         return true;
 412
 413     return false;
 414 }
 415
 416 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 417
 418 } // End namespace Foam
 419
 420 // ************************************************************************* //