src/meshTools/indexedOctree/treeDataTriSurface.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 "treeDataTriSurface.H"
  27 #include "triSurfaceTools.H"
  28 #include "triangleFuncs.H"
  29 #include "indexedOctree.H"
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 defineTypeNameAndDebug(Foam::treeDataTriSurface, 0);
  34
  35
  36 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  37
  38 // // Fast distance to triangle calculation. From
  39 // // "Distance Between Point and Triangle in 3D"
  40 // // David Eberly, Magic Software Inc. Aug. 2003.
  41 // // Works on function Q giving distance to point and tries to minimize this.
  42 // Foam::scalar Foam::treeDataTriSurface::nearestCoords
  43 // (
  44 //     const point& base,
  45 //     const point& E0,
  46 //     const point& E1,
  47 //     const scalar a,
  48 //     const scalar b,
  49 //     const scalar c,
  50 //     const point& P,
  51 //     scalar& s,
  52 //     scalar& t
  53 // )
  54 // {
  55 //     // distance vector
  56 //     const vector D(base - P);
  57
  58 //     // Precalculate distance factors.
  59 //     const scalar d = E0 & D;
  60 //     const scalar e = E1 & D;
  61
  62 //     // Do classification
  63 //     const scalar det = a*c - b*b;
  64
  65 //     s = b*e - c*d;
  66 //     t = b*d - a*e;
  67
  68 //     if (s+t < det)
  69 //     {
  70 //         if (s < 0)
  71 //         {
  72 //             if (t < 0)
  73 //             {
  74 //                 //region 4
  75 //                 if (e > 0)
  76 //                 {
  77 //                     //min on edge t = 0
  78 //                     t = 0;
  79 //                     s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
  80 //                 }
  81 //                 else
  82 //                 {
  83 //                     //min on edge s=0
  84 //                     s = 0;
  85 //                     t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
  86 //                 }
  87 //             }
  88 //             else
  89 //             {
  90 //                 //region 3. Min on edge s = 0
  91 //                 s = 0;
  92 //                 t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
  93 //             }
  94 //         }
  95 //         else if (t < 0)
  96 //         {
  97 //             //region 5
  98 //             t = 0;
  99 //             s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
 100 //         }
 101 //         else
 102 //         {
 103 //             //region 0
 104 //             const scalar invDet = 1/det;
 105 //             s *= invDet;
 106 //             t *= invDet;
 107 //         }
 108 //     }
 109 //     else
 110 //     {
 111 //         if (s < 0)
 112 //         {
 113 //             //region 2
 114 //             const scalar tmp0 = b + d;
 115 //             const scalar tmp1 = c + e;
 116 //             if (tmp1 > tmp0)
 117 //             {
 118 //                 //min on edge s+t=1
 119 //                 const scalar numer = tmp1 - tmp0;
 120 //                 const scalar denom = a-2*b+c;
 121 //                 s = (numer >= denom ? 1 : numer/denom);
 122 //                 t = 1 - s;
 123 //             }
 124 //             else
 125 //             {
 126 //                 //min on edge s=0
 127 //                 s = 0;
 128 //                 t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c));
 129 //             }
 130 //         }
 131 //         else if (t < 0)
 132 //         {
 133 //             //region 6
 134 //             const scalar tmp0 = b + d;
 135 //             const scalar tmp1 = c + e;
 136 //             if (tmp1 > tmp0)
 137 //             {
 138 //                 //min on edge s+t=1
 139 //                 const scalar numer = tmp1 - tmp0;
 140 //                 const scalar denom = a-2*b+c;
 141 //                 s = (numer >= denom ? 1 : numer/denom);
 142 //                 t = 1 - s;
 143 //             }
 144 //             else
 145 //             {
 146 //                 //min on edge t=0
 147 //                 t = 0;
 148 //                 s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a));
 149 //             }
 150 //         }
 151 //         else
 152 //         {
 153 //             //region 1
 154 //             const scalar numer = c+e-(b+d);
 155 //             if (numer <= 0)
 156 //             {
 157 //                 s = 0;
 158 //             }
 159 //             else
 160 //             {
 161 //                 const scalar denom = a-2*b+c;
 162 //                 s = (numer >= denom ? 1 : numer/denom);
 163 //             }
 164 //         }
 165 //         t = 1 - s;
 166 //     }
 167
 168
 169 //     // Calculate distance.
 170 //     // Note: abs should not be needed but truncation error causes problems
 171 //     // with points very close to one of the triangle vertices.
 172 //     // (seen up to -9e-15). Alternatively add some small value.
 173
 174 //     const scalar f = D & D;
 175 //     return Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f);
 176 // }
 177
 178
 179 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 180
 181 // Construct from components
 182 Foam::treeDataTriSurface::treeDataTriSurface
 183 (
 184     const triSurface& surface,
 185     const scalar planarTol
 186 )
 187 :
 188     surface_(surface),
 189     planarTol_(planarTol)
 190 {}
 191
 192
 193 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 194
 195 Foam::pointField Foam::treeDataTriSurface::shapePoints() const
 196 {
 197     const pointField& points = surface_.points();
 198
 199     pointField centres(surface_.size());
 200
 201     forAll(surface_, triI)
 202     {
 203         centres[triI] = surface_[triI].centre(points);
 204     }
 205     return centres;
 206 }
 207
 208
 209 //- Get type of sample (inside/outside/mixed) w.r.t. surface.
 210 Foam::label Foam::treeDataTriSurface::getVolumeType
 211 (
 212     const indexedOctree<treeDataTriSurface>& tree,
 213     const point& sample
 214 ) const
 215 {
 216     // Find nearest point
 217     const treeBoundBox& treeBb = tree.bb();
 218
 219     pointIndexHit pHit = tree.findNearest
 220     (
 221         sample,
 222         max
 223         (
 224             Foam::sqr(GREAT),
 225             Foam::magSqr(treeBb.span())
 226         )
 227     );
 228
 229     if (!pHit.hit())
 230     {
 231         FatalErrorIn("treeDataTriSurface::getVolumeType(..)")
 232             << "treeBb:" << treeBb
 233             << " sample:" << sample
 234             << " pHit:" << pHit
 235             << abort(FatalError);
 236     }
 237
 238     triSurfaceTools::sideType t = triSurfaceTools::surfaceSide
 239     (
 240         surface_,
 241         sample,
 242         pHit.index()
 243     );
 244
 245     if (t == triSurfaceTools::UNKNOWN)
 246     {
 247         return indexedOctree<treeDataTriSurface>::UNKNOWN;
 248     }
 249     else if (t == triSurfaceTools::INSIDE)
 250     {
 251         return indexedOctree<treeDataTriSurface>::INSIDE;
 252     }
 253     else if (t == triSurfaceTools::OUTSIDE)
 254     {
 255         return indexedOctree<treeDataTriSurface>::OUTSIDE;
 256     }
 257     else
 258     {
 259         FatalErrorIn("treeDataTriSurface::getVolumeType(..)")
 260             << "problem" << abort(FatalError);
 261         return indexedOctree<treeDataTriSurface>::UNKNOWN;
 262     }
 263 }
 264
 265
 266 // Check if any point on triangle is inside cubeBb.
 267 bool Foam::treeDataTriSurface::overlaps
 268 (
 269     const label index,
 270     const treeBoundBox& cubeBb
 271 ) const
 272 {
 273     const pointField& points = surface_.points();
 274     const labelledTri& f = surface_[index];
 275
 276     treeBoundBox triBb(points, surface_[index]);
 277
 278     //- For testing: robust one
 279     //return cubeBb.overlaps(triBb);
 280
 281     //- Exact test of triangle intersecting bb
 282
 283     // Quick rejection. If whole bounding box of tri is outside cubeBb then
 284     // there will be no intersection.
 285     if (!cubeBb.overlaps(triBb))
 286     {
 287         return false;
 288     }
 289
 290     // Check if one or more triangle point inside
 291     if (cubeBb.containsAny(points, f))
 292     {
 293         return true;
 294     }
 295
 296     // Triangle points
 297     const point& p0 = points[f[0]];
 298     const point& p1 = points[f[1]];
 299     const point& p2 = points[f[2]];
 300
 301
 302     // Now we have the difficult case: all points are outside but connecting
 303     // edges might go through cube. Use fast intersection of bounding box.
 304
 305     //return triangleFuncs::intersectBbExact(p0, p1, p2, cubeBb);
 306     return triangleFuncs::intersectBb(p0, p1, p2, cubeBb);
 307 }
 308
 309
 310 // Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
 311 // nearestPoint.
 312 void Foam::treeDataTriSurface::findNearest
 313 (
 314     const labelUList& indices,
 315     const point& sample,
 316
 317     scalar& nearestDistSqr,
 318     label& minIndex,
 319     point& nearestPoint
 320 ) const
 321 {
 322     const pointField& points = surface_.points();
 323
 324     forAll(indices, i)
 325     {
 326         label index = indices[i];
 327         const triSurface::FaceType& f = surface_[index];
 328
 329         ////- Possible optimization: do quick rejection of triangle if bounding
 330         ////  sphere does not intersect triangle bounding box. From simplistic
 331         ////  test was not found to speed up things.
 332         //
 333         //// Triangle bounding box.
 334         //point triBbMin = min(p0, min(p1, p2));
 335         //point triBbMax = max(p0, max(p1, p2));
 336         //
 337         //if
 338         //(
 339         //    indexedOctree<treeDataTriSurface>::intersects
 340         //    (
 341         //        triBbMin,
 342         //        triBbMax,
 343         //        nearestDistSqr,
 344         //        sample
 345         //    )
 346         //)
 347         {
 348             // // Get spanning vectors of triangle
 349             // vector base(p1);
 350             // vector E0(p0 - p1);
 351             // vector E1(p2 - p1);
 352
 353             // scalar a(E0& E0);
 354             // scalar b(E0& E1);
 355             // scalar c(E1& E1);
 356
 357             // // Get nearest point in s,t coordinates (s is along E0, t
 358             // // is along E1)
 359             // scalar s;
 360             // scalar t;
 361
 362             // scalar distSqr = nearestCoords
 363             // (
 364             //     base,
 365             //     E0,
 366             //     E1,
 367             //     a,
 368             //     b,
 369             //     c,
 370             //     sample,
 371
 372             //     s,
 373             //     t
 374             // );
 375
 376             pointHit pHit = f.nearestPoint(sample, points);
 377
 378             scalar distSqr = sqr(pHit.distance());
 379
 380             if (distSqr < nearestDistSqr)
 381             {
 382                 nearestDistSqr = distSqr;
 383                 minIndex = index;
 384                 nearestPoint = pHit.rawPoint();
 385             }
 386         }
 387     }
 388 }
 389
 390
 391 // Calculate nearest point to line. Updates (if any) nearestDistSqr, minIndex,
 392 // nearestPoint.
 393 void Foam::treeDataTriSurface::findNearest
 394 (
 395     const labelUList& indices,
 396     const linePointRef& ln,
 397
 398     treeBoundBox& tightest,
 399     label& minIndex,
 400     point& linePoint,
 401     point& nearestPoint
 402 ) const
 403 {
 404     notImplemented
 405     (
 406         "treeDataTriSurface::findNearest(const labelUList&"
 407         ", const linePointRef&, treeBoundBox&, label&, point&, point&) const"
 408     );
 409 }
 410
 411
 412 bool Foam::treeDataTriSurface::intersects
 413 (
 414     const label index,
 415     const point& start,
 416     const point& end,
 417     point& intersectionPoint
 418 ) const
 419 {
 420     const pointField& points = surface_.points();
 421
 422     const triSurface::FaceType& f = surface_[index];
 423
 424     // Do quick rejection test
 425     treeBoundBox triBb(points, f);
 426
 427     const direction startBits(triBb.posBits(start));
 428     const direction endBits(triBb.posBits(end));
 429
 430     if ((startBits & endBits) != 0)
 431     {
 432         // start and end in same block outside of triBb.
 433         return false;
 434     }
 435
 436     const vector dir(end - start);
 437
 438     // Use relative tolerance (from octree) to determine intersection.
 439     pointHit inter = f.intersection
 440     (
 441         start,
 442         dir,
 443         points,
 444         intersection::HALF_RAY,
 445         planarTol_
 446     );
 447
 448     if (inter.hit() && inter.distance() <= 1)
 449     {
 450         // Note: no extra test on whether intersection is in front of us
 451         // since using half_ray.
 452         intersectionPoint = inter.hitPoint();
 453
 454         return true;
 455     }
 456     else
 457     {
 458         return false;
 459     }
 460
 461
 462     //- Using exact intersections
 463     //pointHit info = f.tri(points).intersectionExact(start, end);
 464     //
 465     //if (info.hit())
 466     //{
 467     //    intersectionPoint = info.hitPoint();
 468     //}
 469     //return info.hit();
 470 }
 471
 472
 473 // ************************************************************************* //