applications/utilities/surface/surfaceToPatch/treeDataPrimitivePatch.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "treeDataPrimitivePatch.H"
  27 #include "indexedOctree.H"
  28 #include "triangleFuncs.H"
  29
  30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  31
  32 template
  33 <
  34     class Face,
  35     template<class> class FaceList,
  36     class PointField,
  37     class PointType
  38 >
  39 Foam::scalar
  40 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
  41 tolSqr = sqr(1E-6);
  42
  43
  44 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  45
  46 template
  47 <
  48     class Face,
  49     template<class> class FaceList,
  50     class PointField,
  51     class PointType
  52 >
  53 Foam::treeBoundBox
  54 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
  55 calcBb
  56 (
  57     const pointField& points,
  58     const face& f
  59 )
  60 {
  61     treeBoundBox bb(points[f[0]], points[f[0]]);
  62
  63     for (label fp = 1; fp < f.size(); fp++)
  64     {
  65         const point& p = points[f[fp]];
  66
  67         bb.min() = min(bb.min(), p);
  68         bb.max() = max(bb.max(), p);
  69     }
  70     return bb;
  71 }
  72
  73
  74 template
  75 <
  76     class Face,
  77     template<class> class FaceList,
  78     class PointField,
  79     class PointType
  80 >
  81 void Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
  82 update()
  83 {
  84     if (cacheBb_)
  85     {
  86         bbs_.setSize(patch_.size());
  87
  88         forAll(patch_, i)
  89         {
  90             bbs_[i] = calcBb(patch_.points(), patch_[i]);
  91         }
  92     }
  93 }
  94
  95
  96 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  97
  98 // Construct from components
  99 template
 100 <
 101     class Face,
 102     template<class> class FaceList,
 103     class PointField,
 104     class PointType
 105 >
 106 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
 107 treeDataPrimitivePatch
 108 (
 109     const bool cacheBb,
 110     const PrimitivePatch<Face, FaceList, PointField, PointType>& patch
 111 )
 112 :
 113     patch_(patch),
 114     cacheBb_(cacheBb)
 115 {
 116     update();
 117 }
 118
 119
 120 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 121
 122 template
 123 <
 124     class Face,
 125     template<class> class FaceList,
 126     class PointField,
 127     class PointType
 128 >
 129 Foam::pointField
 130 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
 131 points() const
 132 {
 133     pointField cc(patch_.size());
 134
 135     forAll(patch_, i)
 136     {
 137         cc[i] = patch_[i].centre(patch_.points());
 138     }
 139
 140     return cc;
 141 }
 142
 143
 144 //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
 145 //  Only makes sense for closed surfaces.
 146 template
 147 <
 148     class Face,
 149     template<class> class FaceList,
 150     class PointField,
 151     class PointType
 152 >
 153 Foam::label
 154 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
 155 getVolumeType
 156 (
 157     const indexedOctree
 158     <
 159         treeDataPrimitivePatch
 160         <
 161             Face,
 162             FaceList,
 163             PointField,
 164             PointType
 165         >
 166     >& oc,
 167     const point& sample
 168 ) const
 169 {
 170     // Need to determine whether sample is 'inside' or 'outside'
 171     // Done by finding nearest face. This gives back a face which is
 172     // guaranteed to contain nearest point. This point can be
 173     // - in interior of face: compare to face normal
 174     // - on edge of face: compare to edge normal
 175     // - on point of face: compare to point normal
 176     // Unfortunately the octree does not give us back the intersection point
 177     // or where on the face it has hit so we have to recreate all that
 178     // information.
 179
 180
 181     // Find nearest face to sample
 182     pointIndexHit info = oc.findNearest(sample, sqr(GREAT));
 183
 184     if (info.index() == -1)
 185     {
 186         FatalErrorIn
 187         (
 188             "treeDataPrimitivePatch::getSampleType"
 189             "(indexedOctree<treeDataPrimitivePatch>&, const point&)"
 190         )   << "Could not find " << sample << " in octree."
 191             << abort(FatalError);
 192     }
 193
 194
 195     // Get actual intersection point on face
 196     label faceI = info.index();
 197
 198     if (debug & 2)
 199     {
 200         Pout<< "getSampleType : sample:" << sample
 201             << " nearest face:" << faceI;
 202     }
 203
 204     const pointField& points = patch_.localPoints();
 205     const face& f = patch_.localFaces()[faceI];
 206
 207     // Retest to classify where on face info is. Note: could be improved. We
 208     // already have point.
 209
 210     pointHit curHit = f.nearestPoint(sample, points);
 211     const vector area = f.normal(points);
 212     const point& curPt = curHit.rawPoint();
 213
 214     //
 215     // 1] Check whether sample is above face
 216     //
 217
 218     if (curHit.hit())
 219     {
 220         // Nearest point inside face. Compare to face normal.
 221
 222         if (debug & 2)
 223         {
 224             Pout<< " -> face hit:" << curPt
 225                 << " comparing to face normal " << area << endl;
 226         }
 227         return indexedOctree<treeDataPrimitivePatch>::getSide
 228         (
 229             area,
 230             sample - curPt
 231         );
 232     }
 233
 234     if (debug & 2)
 235     {
 236         Pout<< " -> face miss:" << curPt;
 237     }
 238
 239     //
 240     // 2] Check whether intersection is on one of the face vertices or
 241     //    face centre
 242     //
 243
 244     const scalar typDimSqr = mag(area) + VSMALL;
 245
 246     forAll(f, fp)
 247     {
 248         if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < tolSqr)
 249         {
 250             // Face intersection point equals face vertex fp
 251
 252             // Calculate point normal (wrong: uses face normals instead of
 253             // triangle normals)
 254
 255             return indexedOctree<treeDataPrimitivePatch>::getSide
 256             (
 257                 patch_.pointNormals()[f[fp]],
 258                 sample - curPt
 259             );
 260         }
 261     }
 262
 263     const point fc(f.centre(points));
 264
 265     if ((magSqr(fc - curPt)/typDimSqr) < tolSqr)
 266     {
 267         // Face intersection point equals face centre. Normal at face centre
 268         // is already average of face normals
 269
 270         if (debug & 2)
 271         {
 272             Pout<< " -> centre hit:" << fc
 273                 << " distance:" << magSqr(fc - curPt)/typDimSqr << endl;
 274         }
 275
 276         return indexedOctree<treeDataPrimitivePatch>::getSide
 277         (
 278             area,
 279             sample - curPt
 280         );
 281     }
 282
 283
 284
 285     //
 286     // 3] Get the 'real' edge the face intersection is on
 287     //
 288
 289     const labelList& fEdges = patch_.faceEdges()[faceI];
 290
 291     forAll(fEdges, fEdgeI)
 292     {
 293         label edgeI = fEdges[fEdgeI];
 294         const edge& e = patch_.edges()[edgeI];
 295
 296         pointHit edgeHit = e.line(points).nearestDist(sample);
 297
 298         if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
 299         {
 300             // Face intersection point lies on edge e
 301
 302             // Calculate edge normal (wrong: uses face normals instead of
 303             // triangle normals)
 304             const labelList& eFaces = patch_.edgeFaces()[edgeI];
 305
 306             vector edgeNormal(vector::zero);
 307
 308             forAll(eFaces, i)
 309             {
 310                 edgeNormal += patch_.faceNormal()[eFaces[i]];
 311             }
 312
 313             if (debug & 2)
 314             {
 315                 Pout<< " -> real edge hit point:" << edgeHit.rawPoint()
 316                     << " comparing to edge normal:" << edgeNormal
 317                     << endl;
 318             }
 319
 320             // Found face intersection point on this edge. Compare to edge
 321             // normal
 322             return indexedOctree<treeDataPrimitivePatch>::getSide
 323             (
 324                 edgeNormal,
 325                 sample - curPt
 326             );
 327         }
 328     }
 329
 330
 331     //
 332     // 4] Get the internal edge the face intersection is on
 333     //
 334
 335     forAll(f, fp)
 336     {
 337         pointHit edgeHit = linePointRef
 338         (
 339             points[f[fp]],
 340             fc
 341         ).nearestDist(sample);
 342
 343         if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
 344         {
 345             // Face intersection point lies on edge between two face triangles
 346
 347             // Calculate edge normal as average of the two triangle normals
 348             vector e = points[f[fp]] - fc;
 349             vector ePrev = points[f[f.rcIndex(fp)]] - fc;
 350             vector eNext = points[f[f.fcIndex(fp)]] - fc;
 351
 352             vector nLeft = ePrev ^ e;
 353             nLeft /= mag(nLeft) + VSMALL;
 354
 355             vector nRight = e ^ eNext;
 356             nRight /= mag(nRight) + VSMALL;
 357
 358             if (debug & 2)
 359             {
 360                 Pout<< " -> internal edge hit point:" << edgeHit.rawPoint()
 361                     << " comparing to edge normal "
 362                     << 0.5*(nLeft + nRight)
 363                     << endl;
 364             }
 365
 366             // Found face intersection point on this edge. Compare to edge
 367             // normal
 368             return indexedOctree<treeDataPrimitivePatch>::getSide
 369             (
 370                 0.5*(nLeft + nRight),
 371                 sample - curPt
 372             );
 373         }
 374     }
 375
 376     if (debug & 2)
 377     {
 378         Pout<< "Did not find sample " << sample
 379             << " anywhere related to nearest face " << faceI << endl
 380             << "Face:";
 381
 382         forAll(f, fp)
 383         {
 384             Pout<< "    vertex:" << f[fp] << "  coord:" << points[f[fp]]
 385                 << endl;
 386         }
 387     }
 388
 389     // Can't determine status of sample with respect to nearest face.
 390     // Either
 391     // - tolerances are wrong. (if e.g. face has zero area)
 392     // - or (more likely) surface is not closed.
 393
 394     return indexedOctree<treeDataPrimitivePatch>::UNKNOWN;
 395 }
 396
 397
 398 // Check if any point on shape is inside cubeBb.
 399 template
 400 <
 401     class Face,
 402     template<class> class FaceList,
 403     class PointField,
 404     class PointType
 405 >
 406 bool
 407 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
 408 overlaps
 409 (
 410     const label index,
 411     const treeBoundBox& cubeBb
 412 ) const
 413 {
 414     // 1. Quick rejection: bb does not intersect face bb at all
 415     if (cacheBb_)
 416     {
 417         if (!cubeBb.overlaps(bbs_[index]))
 418         {
 419             return false;
 420         }
 421     }
 422     else
 423     {
 424         if (!cubeBb.overlaps(calcBb(patch_.points(), patch_[index])))
 425         {
 426             return false;
 427         }
 428     }
 429
 430
 431     // 2. Check if one or more face points inside
 432
 433     const pointField& points = patch_.points();
 434     const face& f = patch_[index];
 435
 436     forAll(f, fp)
 437     {
 438         if (cubeBb.contains(points[f[fp]]))
 439         {
 440             return true;
 441         }
 442     }
 443
 444     // 3. Difficult case: all points are outside but connecting edges might
 445     // go through cube. Use triangle-bounding box intersection.
 446     const point fc = f.centre(points);
 447
 448     forAll(f, fp)
 449     {
 450         bool triIntersects = triangleFuncs::intersectBb
 451         (
 452             points[f[fp]],
 453             points[f[f.fcIndex(fp)]],
 454             fc,
 455             cubeBb
 456         );
 457
 458         if (triIntersects)
 459         {
 460             return true;
 461         }
 462     }
 463     return false;
 464 }
 465
 466
 467 // Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
 468 // nearestPoint.
 469 template
 470 <
 471     class Face,
 472     template<class> class FaceList,
 473     class PointField,
 474     class PointType
 475 >
 476 void
 477 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
 478 findNearest
 479 (
 480     const labelList& indices,
 481     const point& sample,
 482
 483     scalar& nearestDistSqr,
 484     label& minIndex,
 485     point& nearestPoint
 486 ) const
 487 {
 488     const pointField& points = patch_.points();
 489
 490     forAll(indices, i)
 491     {
 492         label index = indices[i];
 493
 494         const face& f = patch_[index];
 495
 496         pointHit nearHit = f.nearestPoint(sample, points);
 497         scalar distSqr = sqr(nearHit.distance());
 498
 499         if (distSqr < nearestDistSqr)
 500         {
 501             nearestDistSqr = distSqr;
 502             minIndex = index;
 503             nearestPoint = nearHit.rawPoint();
 504         }
 505     }
 506 }
 507
 508
 509 template
 510 <
 511     class Face,
 512     template<class> class FaceList,
 513     class PointField,
 514     class PointType
 515 >
 516 bool
 517 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
 518 intersects
 519 (
 520     const label index,
 521     const point& start,
 522     const point& end,
 523     point& intersectionPoint
 524 ) const
 525 {
 526     // Do quick rejection test
 527     if (cacheBb_)
 528     {
 529         const treeBoundBox& faceBb = bbs_[index];
 530
 531         if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
 532         {
 533             // start and end in same block outside of faceBb.
 534             return false;
 535         }
 536     }
 537
 538     const pointField& points = patch_.points();
 539     const face& f = patch_[index];
 540     const point fc = f.centre(points);
 541     const vector dir(end - start);
 542
 543     pointHit inter = patch_[index].intersection
 544     (
 545         start,
 546         dir,
 547         fc,
 548         points,
 549         intersection::HALF_RAY
 550     );
 551
 552     if (inter.hit() && inter.distance() <= 1)
 553     {
 554         // Note: no extra test on whether intersection is in front of us
 555         // since using half_ray
 556         intersectionPoint = inter.hitPoint();
 557         return true;
 558     }
 559     else
 560     {
 561         return false;
 562     }
 563 }
 564
 565
 566 // ************************************************************************* //