src/meshTools/triSurface/octreeData/octreeDataTriSurface.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2010 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
  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 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "octreeDataTriSurface.H"
  29
  30 #include "labelList.H"
  31 #include "treeBoundBox.H"
  32 #include "faceList.H"
  33 #include "triPointRef.H"
  34 #include "octree.H"
  35 #include "triSurfaceTools.H"
  36 #include "triangleFuncs.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 defineTypeNameAndDebug(Foam::octreeDataTriSurface, 0);
  41
  42 Foam::scalar Foam::octreeDataTriSurface::tol(1E-6);
  43
  44
  45 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  46
  47 // Fast distance to triangle calculation. From
  48 // "Distance Between Point and Trangle in 3D"
  49 // David Eberly, Magic Software Inc. Aug. 2003.
  50 // Works on function Q giving distance to point and tries to minimize this.
  51 void Foam::octreeDataTriSurface::nearestCoords
  52 (
  53     const point& base,
  54     const point& E0,
  55     const point& E1,
  56     const scalar a,
  57     const scalar b,
  58     const scalar c,
  59     const point& P,
  60     scalar& s,
  61     scalar& t
  62 )
  63 {
  64     // distance vector
  65     const vector D(base - P);
  66
  67     // Precalculate distance factors.
  68     const scalar d = E0 & D;
  69     const scalar e = E1 & D;
  70
  71     // Do classification
  72     const scalar det = a*c - b*b;
  73
  74     s = b*e - c*d;
  75     t = b*d - a*e;
  76
  77     if (s+t < det)
  78     {
  79         if (s < 0)
  80         {
  81             if (t < 0)
  82             {
  83                 //region 4
  84                 if (e > 0)
  85                 {
  86                     //min on edge t = 0
  87                     t = 0;
  88                     s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
  89                 }
  90                 else
  91                 {
  92                     //min on edge s=0
  93                     s = 0;
  94                     t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
  95                 }
  96             }
  97             else
  98             {
  99                 //region 3. Min on edge s = 0
 100                 s = 0;
 101                 t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
 102             }
 103         }
 104         else if (t < 0)
 105         {
 106             //region 5
 107             t = 0;
 108             s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
 109         }
 110         else
 111         {
 112             //region 0
 113             const scalar invDet = 1/det;
 114             s *= invDet;
 115             t *= invDet;
 116         }
 117     }
 118     else
 119     {
 120         if (s < 0)
 121         {
 122             //region 2
 123             const scalar tmp0 = b + d;
 124             const scalar tmp1 = c + e;
 125             if (tmp1 > tmp0)
 126             {
 127                 //min on edge s+t=1
 128                 const scalar numer = tmp1 - tmp0;
 129                 const scalar denom = a-2*b+c;
 130                 s = (numer >= denom ? 1 : numer/denom);
 131                 t = 1 - s;
 132             }
 133             else
 134             {
 135                 //min on edge s=0
 136                 s = 0;
 137                 t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c));
 138             }
 139         }
 140         else if (t < 0)
 141         {
 142             //region 6
 143             const scalar tmp0 = b + d;
 144             const scalar tmp1 = c + e;
 145             if (tmp1 > tmp0)
 146             {
 147                 //min on edge s+t=1
 148                 const scalar numer = tmp1 - tmp0;
 149                 const scalar denom = a-2*b+c;
 150                 s = (numer >= denom ? 1 : numer/denom);
 151                 t = 1 - s;
 152             }
 153             else
 154             {
 155                 //min on edge t=0
 156                 t = 0;
 157                 s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a));
 158             }
 159         }
 160         else
 161         {
 162             //region 1
 163             const scalar numer = c+e-(b+d);
 164             if (numer <= 0)
 165             {
 166                 s = 0;
 167             }
 168             else
 169             {
 170                 const scalar denom = a-2*b+c;
 171                 s = (numer >= denom ? 1 : numer/denom);
 172             }
 173         }
 174         t = 1 - s;
 175     }
 176
 177
 178     // Calculate distance.
 179     // Note: abs should not be needed but truncation error causes problems
 180     // with points very close to one of the triangle vertices.
 181     // (seen up to -9e-15). Alternatively add some small value.
 182
 183     // const scalar f = D & D;
 184     // return a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f + SMALL;
 185     // return Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f);
 186 }
 187
 188
 189 Foam::point Foam::octreeDataTriSurface::nearestPoint
 190 (
 191     const label index,
 192     const point& p
 193 ) const
 194 {
 195     scalar s;
 196     scalar t;
 197
 198     nearestCoords
 199     (
 200         base_[index],
 201         E0_[index],
 202         E1_[index],
 203         a_[index],
 204         b_[index],
 205         c_[index],
 206         p,
 207         s,
 208         t
 209     );
 210
 211     return base_[index] + s*E0_[index] + t*E1_[index];
 212 }
 213
 214
 215 // Helper function to calculate tight fitting bounding boxes.
 216 Foam::treeBoundBoxList Foam::octreeDataTriSurface::calcBb
 217 (
 218     const triSurface& surf
 219 )
 220 {
 221     treeBoundBoxList allBb(surf.size(), treeBoundBox::invertedBox);
 222
 223     const labelListList& pointFcs = surf.pointFaces();
 224     const pointField& localPts = surf.localPoints();
 225
 226     forAll(pointFcs, pointI)
 227     {
 228         const labelList& myFaces = pointFcs[pointI];
 229         const point& vertCoord = localPts[pointI];
 230
 231         forAll(myFaces, myFaceI)
 232         {
 233             // Update bb
 234             label faceI = myFaces[myFaceI];
 235
 236             treeBoundBox& bb = allBb[faceI];
 237
 238             bb.min() = min(bb.min(), vertCoord);
 239             bb.max() = max(bb.max(), vertCoord);
 240         }
 241     }
 242
 243     return allBb;
 244 }
 245
 246
 247 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 248
 249 // Construct from components
 250 Foam::octreeDataTriSurface::octreeDataTriSurface(const triSurface& surface)
 251 :
 252     surface_(surface),
 253     allBb_(calcBb(surface_)),
 254     base_(surface_.size()),
 255     E0_(surface_.size()),
 256     E1_(surface_.size()),
 257     a_(surface_.size()),
 258     b_(surface_.size()),
 259     c_(surface_.size())
 260 {
 261     // Precalculate factors for distance calculation
 262     const pointField& points = surface_.points();
 263
 264     forAll(surface_, faceI)
 265     {
 266         const labelledTri& f = surface_[faceI];
 267
 268         // Calculate base and spanning vectors of triangles
 269         base_[faceI] = points[f[1]];
 270         E0_[faceI] = points[f[0]] - points[f[1]];
 271         E1_[faceI] = points[f[2]] - points[f[1]];
 272
 273         a_[faceI] = E0_[faceI] & E0_[faceI];
 274         b_[faceI] = E0_[faceI] & E1_[faceI];
 275         c_[faceI] = E1_[faceI] & E1_[faceI];
 276     }
 277 }
 278
 279
 280 // Construct from components
 281 Foam::octreeDataTriSurface::octreeDataTriSurface
 282 (
 283     const triSurface& surface,
 284     const treeBoundBoxList& allBb
 285 )
 286 :
 287     surface_(surface),
 288     allBb_(allBb),
 289     base_(surface_.size()),
 290     E0_(surface_.size()),
 291     E1_(surface_.size()),
 292     a_(surface_.size()),
 293     b_(surface_.size()),
 294     c_(surface_.size())
 295 {
 296     const pointField& points = surface_.points();
 297
 298     forAll(surface_, faceI)
 299     {
 300         const labelledTri& f = surface_[faceI];
 301
 302         // Calculate base and spanning vectors of triangles
 303         base_[faceI] = points[f[1]];
 304         E0_[faceI] = points[f[0]] - points[f[1]];
 305         E1_[faceI] = points[f[2]] - points[f[1]];
 306
 307         a_[faceI] = E0_[faceI] & E0_[faceI];
 308         b_[faceI] = E0_[faceI] & E1_[faceI];
 309         c_[faceI] = E1_[faceI] & E1_[faceI];
 310     }
 311 }
 312
 313
 314 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 315
 316 Foam::label Foam::octreeDataTriSurface::getSampleType
 317 (
 318     const octree<octreeDataTriSurface>& oc,
 319     const point& sample
 320 ) const
 321
 322 {
 323     treeBoundBox tightest(treeBoundBox::greatBox);
 324     scalar tightestDist(treeBoundBox::great);
 325
 326     // Find nearest face to sample
 327     label faceI = oc.findNearest(sample, tightest, tightestDist);
 328
 329     if (debug & 2)
 330     {
 331         Pout<< "getSampleType : sample:" << sample
 332             << " nearest face:" << faceI;
 333     }
 334
 335     if (faceI == -1)
 336     {
 337         FatalErrorIn
 338         (
 339             "octreeDataTriSurface::getSampleType"
 340             "(octree<octreeDataTriSurface>&, const point&)"
 341         )   << "Could not find " << sample << " in octree."
 342             << abort(FatalError);
 343     }
 344
 345     const pointField& pts = surface_.points();
 346     const labelledTri& f = surface_[faceI];
 347
 348     pointHit curHit = triPointRef
 349     (
 350         pts[f[0]],
 351         pts[f[1]],
 352         pts[f[2]]
 353     ).nearestPoint(sample);
 354
 355     // Get normal according to position on face. On point -> pointNormal,
 356     // on edge-> edge normal, face normal on interior.
 357     vector n
 358     (
 359         triSurfaceTools::surfaceNormal
 360         (
 361             surface_,
 362             faceI,
 363             curHit.rawPoint()
 364         )
 365     );
 366
 367     return
 368         octree<octreeDataTriSurface>::getVolType(n, sample - curHit.rawPoint());
 369 }
 370
 371
 372 // Check if any point on triangle is inside cubeBb.
 373 bool Foam::octreeDataTriSurface::overlaps
 374 (
 375     const label index,
 376     const treeBoundBox& cubeBb
 377 ) const
 378 {
 379     //return cubeBb.overlaps(allBb_[index]);
 380
 381     //- Exact test of triangle intersecting bb
 382
 383     // Quick rejection.
 384     if (!cubeBb.overlaps(allBb_[index]))
 385     {
 386         return false;
 387     }
 388
 389     // Triangle points
 390     const pointField& points = surface_.points();
 391     const labelledTri& f = surface_[index];
 392     const point& p0 = points[f[0]];
 393     const point& p1 = points[f[1]];
 394     const point& p2 = points[f[2]];
 395
 396     // Check if one or more triangle point inside
 397     if (cubeBb.contains(p0) || cubeBb.contains(p1) || cubeBb.contains(p2))
 398     {
 399         // One or more points inside
 400         return true;
 401     }
 402
 403     // Now we have the difficult case: all points are outside but connecting
 404     // edges might go through cube. Use fast intersection of bounding box.
 405
 406     return triangleFuncs::intersectBb(p0, p1, p2, cubeBb);
 407 }
 408
 409
 410 bool Foam::octreeDataTriSurface::contains
 411 (
 412     const label,
 413     const point&
 414 ) const
 415 {
 416     notImplemented
 417     (
 418         "octreeDataTriSurface::contains(const label, const point&)"
 419     );
 420
 421     return false;
 422 }
 423
 424
 425 bool Foam::octreeDataTriSurface::intersects
 426 (
 427     const label index,
 428     const point& start,
 429     const point& end,
 430     point& intersectionPoint
 431 ) const
 432 {
 433     if (mag(surface_.faceNormals()[index]) < VSMALL)
 434     {
 435         return false;
 436     }
 437
 438     const pointField& points = surface_.points();
 439
 440     const labelledTri& f = surface_[index];
 441
 442     triPointRef tri(points[f[0]], points[f[1]], points[f[2]]);
 443
 444     const vector dir(end - start);
 445
 446     // Disable picking up intersections behind us.
 447     scalar oldTol = intersection::setPlanarTol(0.0);
 448
 449     pointHit inter = tri.ray(start, dir, intersection::HALF_RAY);
 450
 451     intersection::setPlanarTol(oldTol);
 452
 453     if (inter.hit() && inter.distance() <= mag(dir))
 454     {
 455         // Note: no extra test on whether intersection is in front of us
 456         // since using half_ray AND zero tolerance. (note that tolerance
 457         // is used to look behind us)
 458         intersectionPoint = inter.hitPoint();
 459
 460         return true;
 461     }
 462     else
 463     {
 464         return false;
 465     }
 466 }
 467
 468
 469 bool Foam::octreeDataTriSurface::findTightest
 470 (
 471     const label index,
 472     const point& sample,
 473     treeBoundBox& tightest
 474 ) const
 475 {
 476
 477     // get nearest and furthest away vertex
 478     point myNear, myFar;
 479     allBb_[index].calcExtremities(sample, myNear, myFar);
 480
 481     const point dist = myFar - sample;
 482     scalar myFarDist = mag(dist);
 483
 484     point tightestNear, tightestFar;
 485     tightest.calcExtremities(sample, tightestNear, tightestFar);
 486
 487     scalar tightestFarDist = mag(tightestFar - sample);
 488
 489     if (tightestFarDist < myFarDist)
 490     {
 491         // Keep current tightest.
 492         return false;
 493     }
 494     else
 495     {
 496         // Construct bb around sample and myFar
 497         const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
 498
 499         tightest.min() = sample - dist2;
 500         tightest.max() = sample + dist2;
 501
 502         return true;
 503     }
 504 }
 505
 506
 507 // Determine numerical value of sign of sample compared to shape at index
 508 Foam::scalar Foam::octreeDataTriSurface::calcSign
 509 (
 510     const label index,
 511     const point& sample,
 512     vector& n
 513 ) const
 514 {
 515     n = surface_.faceNormals()[index];
 516
 517     const labelledTri& tri = surface_[index];
 518
 519     // take vector from sample to any point on triangle (we use vertex 0)
 520     vector vec = sample - surface_.points()[tri[0]];
 521
 522     vec /= mag(vec) + VSMALL;
 523
 524     return n & vec;
 525 }
 526
 527
 528 // Calculate nearest point to sample on/in shapei. !Does not set nearest
 529 Foam::scalar Foam::octreeDataTriSurface::calcNearest
 530 (
 531     const label index,
 532     const point& sample,
 533     point&
 534 ) const
 535 {
 536     return mag(nearestPoint(index, sample) - sample);
 537 }
 538
 539
 540 // Calculate nearest point on/in shapei
 541 Foam::scalar Foam::octreeDataTriSurface::calcNearest
 542 (
 543     const label index,
 544     const linePointRef& ln,
 545     point& linePt,
 546     point& shapePt
 547 ) const
 548 {
 549     notImplemented
 550     (
 551         "octreeDataTriSurface::calcNearest"
 552         "(const label, const linePointRef&, point& linePt, point&)"
 553     );
 554     return GREAT;
 555 }
 556
 557
 558 void Foam::octreeDataTriSurface::write
 559 (
 560     Ostream& os,
 561     const label index
 562 ) const
 563 {
 564     os << surface_[index] << token::SPACE << allBb_[index];
 565 }
 566
 567
 568 // ************************************************************************* //