src/meshTools/indexedOctree/treeDataTriSurface.C

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