src/meshTools/octree/octreeDataFace.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 "octreeDataFace.H"
  27 #include "labelList.H"
  28 #include "polyMesh.H"
  29 #include "octree.H"
  30 #include "polyPatch.H"
  31 #include "triangleFuncs.H"
  32 #include "linePointRef.H"
  33
  34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  35
  36 defineTypeNameAndDebug(Foam::octreeDataFace, 0);
  37
  38 Foam::scalar Foam::octreeDataFace::tol(1E-6);
  39
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43 void Foam::octreeDataFace::calcBb()
  44 {
  45     allBb_.setSize(meshFaces_.size());
  46     allBb_ = treeBoundBox::invertedBox;
  47
  48     forAll(meshFaces_, i)
  49     {
  50         // Update bb of face
  51         treeBoundBox& myBb = allBb_[i];
  52
  53         const face& f = mesh_.faces()[meshFaces_[i]];
  54
  55         forAll(f, faceVertexI)
  56         {
  57             const point& coord = mesh_.points()[f[faceVertexI]];
  58
  59             myBb.min() = min(myBb.min(), coord);
  60             myBb.max() = max(myBb.max(), coord);
  61         }
  62     }
  63 }
  64
  65
  66 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  67
  68 // Construct from selected mesh faces
  69 Foam::octreeDataFace::octreeDataFace
  70 (
  71     const primitiveMesh& mesh,
  72     const labelList& meshFaces,
  73     const treeBoundBoxList& allBb
  74 )
  75 :
  76     mesh_(mesh),
  77     meshFaces_(meshFaces),
  78     allBb_(allBb)
  79 {}
  80
  81
  82 // Construct from selected mesh faces. Bounding box calculated.
  83 Foam::octreeDataFace::octreeDataFace
  84 (
  85     const primitiveMesh& mesh,
  86     const labelList& meshFaces
  87 )
  88 :
  89     mesh_(mesh),
  90     meshFaces_(meshFaces),
  91     allBb_(meshFaces_.size())
  92 {
  93     // Generate tight fitting bounding box
  94     calcBb();
  95 }
  96
  97
  98 // Construct from selected mesh faces
  99 Foam::octreeDataFace::octreeDataFace
 100 (
 101     const primitiveMesh& mesh,
 102     const UList<const labelList*>& meshFaceListPtrs,
 103     const UList<const treeBoundBoxList*>& bbListPtrs
 104 )
 105 :
 106     mesh_(mesh),
 107     meshFaces_(0),
 108     allBb_(0)
 109 {
 110     label faceI = 0;
 111
 112     forAll(meshFaceListPtrs, listI)
 113     {
 114         faceI += meshFaceListPtrs[listI]->size();
 115     }
 116
 117     meshFaces_.setSize(faceI);
 118     allBb_.setSize(faceI);
 119
 120     faceI = 0;
 121
 122     forAll(meshFaceListPtrs, listI)
 123     {
 124         const labelList& meshFaces = *meshFaceListPtrs[listI];
 125         const treeBoundBoxList& allBb = *bbListPtrs[listI];
 126
 127         forAll(meshFaces, meshFaceI)
 128         {
 129             meshFaces_[faceI] = meshFaces[meshFaceI];
 130             allBb_[faceI] = allBb[meshFaceI];
 131             faceI++;
 132         }
 133     }
 134 }
 135
 136
 137 // Construct from selected mesh faces. Bounding box calculated.
 138 Foam::octreeDataFace::octreeDataFace
 139 (
 140     const primitiveMesh& mesh,
 141     const UList<const labelList*>& meshFaceListPtrs
 142 )
 143 :
 144     mesh_(mesh),
 145     meshFaces_(0)
 146 {
 147     label faceI = 0;
 148
 149     forAll(meshFaceListPtrs, listI)
 150     {
 151         faceI += meshFaceListPtrs[listI]->size();
 152     }
 153
 154     meshFaces_.setSize(faceI);
 155
 156     faceI = 0;
 157
 158     forAll(meshFaceListPtrs, listI)
 159     {
 160         const labelList& meshFaces = *meshFaceListPtrs[listI];
 161
 162         forAll(meshFaces, meshFaceI)
 163         {
 164             meshFaces_[faceI++] = meshFaces[meshFaceI];
 165         }
 166     }
 167
 168     // Generate tight fitting bounding box
 169     calcBb();
 170 }
 171
 172
 173 // Construct from all faces in polyPatch. Bounding box calculated.
 174 Foam::octreeDataFace::octreeDataFace(const polyPatch& patch)
 175 :
 176     mesh_(patch.boundaryMesh().mesh()),
 177     meshFaces_(patch.size())
 178 {
 179     forAll(patch, patchFaceI)
 180     {
 181         meshFaces_[patchFaceI] = patch.start() + patchFaceI;
 182     }
 183
 184     // Generate tight fitting bounding box
 185     calcBb();
 186 }
 187
 188
 189 // Construct from primitiveMesh. Inserts all boundary faces.
 190 Foam::octreeDataFace::octreeDataFace(const primitiveMesh& mesh)
 191 :
 192     mesh_(mesh),
 193     meshFaces_(0),
 194     allBb_(0)
 195 {
 196     // Size storage
 197     meshFaces_.setSize(mesh_.nFaces() - mesh_.nInternalFaces());
 198
 199     // Set info for all boundary faces.
 200     label boundaryFaceI = 0;
 201
 202     for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
 203     {
 204         meshFaces_[boundaryFaceI++] = faceI;
 205     }
 206
 207     // Generate tight fitting bounding box
 208     calcBb();
 209 }
 210
 211
 212 // Construct as copy
 213 Foam::octreeDataFace::octreeDataFace(const octreeDataFace& shapes)
 214 :
 215     mesh_(shapes.mesh()),
 216     meshFaces_(shapes.meshFaces()),
 217     allBb_(shapes.allBb())
 218 {}
 219
 220
 221 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 222
 223 Foam::octreeDataFace::~octreeDataFace()
 224 {}
 225
 226
 227 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 228
 229 Foam::label Foam::octreeDataFace::getSampleType
 230 (
 231     const octree<octreeDataFace>& oc,
 232     const point& sample
 233 ) const
 234 {
 235     // Need to determine whether sample is 'inside' or 'outside'
 236     // Done by finding nearest face. This gives back a face which is
 237     // guaranteed to contain nearest point. This point can be
 238     // - in interior of face: compare to face normal
 239     // - on edge of face: compare to edge normal
 240     // - on point of face: compare to point normal
 241     // Unfortunately the octree does not give us back the intersection point
 242     // or where on the face it has hit so we have to recreate all that
 243     // information.
 244
 245     treeBoundBox tightest(treeBoundBox::greatBox);
 246     scalar tightestDist(treeBoundBox::great);
 247     // Find nearest face to sample
 248     label index = oc.findNearest(sample, tightest, tightestDist);
 249
 250     if (index == -1)
 251     {
 252         FatalErrorIn
 253         (
 254             "octreeDataFace::getSampleType"
 255             "(octree<octreeDataFace>&, const point&)"
 256         )   << "Could not find " << sample << " in octree."
 257             << abort(FatalError);
 258     }
 259
 260
 261     // Get actual intersection point on face
 262     label faceI = meshFaces_[index];
 263
 264     if (debug & 2)
 265     {
 266         Pout<< "getSampleType : sample:" << sample
 267             << " nearest face:" << faceI;
 268     }
 269
 270     const face& f = mesh_.faces()[faceI];
 271
 272     const pointField& points = mesh_.points();
 273
 274     pointHit curHit = f.nearestPoint(sample, points);
 275
 276     //
 277     // 1] Check whether sample is above face
 278     //
 279
 280     if (curHit.hit())
 281     {
 282         // Simple case. Compare to face normal.
 283
 284         if (debug & 2)
 285         {
 286             Pout<< " -> face hit:" << curHit.hitPoint()
 287                 << " comparing to face normal " << mesh_.faceAreas()[faceI]
 288                 << endl;
 289         }
 290         return octree<octreeDataFace>::getVolType
 291         (
 292             mesh_.faceAreas()[faceI],
 293             sample - curHit.hitPoint()
 294         );
 295     }
 296
 297     if (debug & 2)
 298     {
 299         Pout<< " -> face miss:" << curHit.missPoint();
 300     }
 301
 302     //
 303     // 2] Check whether intersection is on one of the face vertices or
 304     //    face centre
 305     //
 306
 307     scalar typDim = sqrt(mag(mesh_.faceAreas()[faceI])) + VSMALL;
 308
 309     forAll(f, fp)
 310     {
 311         if ((mag(points[f[fp]] - curHit.missPoint())/typDim) < tol)
 312         {
 313             // Face intersection point equals face vertex fp
 314
 315             // Calculate point normal (wrong: uses face normals instead of
 316             // triangle normals)
 317             const labelList& myFaces = mesh_.pointFaces()[f[fp]];
 318
 319             vector pointNormal(vector::zero);
 320
 321             forAll(myFaces, myFaceI)
 322             {
 323                 if (myFaces[myFaceI] >= mesh_.nInternalFaces())
 324                 {
 325                     vector n = mesh_.faceAreas()[myFaces[myFaceI]];
 326                     n /= mag(n) + VSMALL;
 327
 328                     pointNormal += n;
 329                 }
 330             }
 331
 332             if (debug & 2)
 333             {
 334                     Pout<< " -> face point hit :" << points[f[fp]]
 335                         << " point normal:" << pointNormal
 336                         << " distance:"
 337                         << mag(points[f[fp]] - curHit.missPoint())/typDim
 338                         << endl;
 339             }
 340             return octree<octreeDataFace>::getVolType
 341             (
 342                 pointNormal,
 343                 sample - curHit.missPoint()
 344             );
 345         }
 346     }
 347     if ((mag(mesh_.faceCentres()[faceI] - curHit.missPoint())/typDim) < tol)
 348     {
 349         // Face intersection point equals face centre. Normal at face centre
 350         // is already average of face normals
 351
 352         if (debug & 2)
 353         {
 354             Pout<< " -> centre hit:" << mesh_.faceCentres()[faceI]
 355                 << " distance:"
 356                 << mag(mesh_.faceCentres()[faceI] - curHit.missPoint())/typDim
 357                 << endl;
 358         }
 359
 360         return octree<octreeDataFace>::getVolType
 361         (
 362             mesh_.faceAreas()[faceI],
 363             sample - curHit.missPoint()
 364         );
 365     }
 366
 367
 368     //
 369     // 3] Get the 'real' edge the face intersection is on
 370     //
 371
 372     const labelList& myEdges = mesh_.faceEdges()[faceI];
 373
 374     forAll(myEdges, myEdgeI)
 375     {
 376         const edge& e = mesh_.edges()[myEdges[myEdgeI]];
 377
 378         pointHit edgeHit = line<point, const point&>
 379         (
 380             points[e.start()],
 381             points[e.end()]
 382         ).nearestDist(sample);
 383
 384
 385         if ((mag(edgeHit.rawPoint() - curHit.missPoint())/typDim) < tol)
 386         {
 387             // Face intersection point lies on edge e
 388
 389             // Calculate edge normal (wrong: uses face normals instead of
 390             // triangle normals)
 391             const labelList& myFaces = mesh_.edgeFaces()[myEdges[myEdgeI]];
 392
 393             vector edgeNormal(vector::zero);
 394
 395             forAll(myFaces, myFaceI)
 396             {
 397                 if (myFaces[myFaceI] >= mesh_.nInternalFaces())
 398                 {
 399                     vector n = mesh_.faceAreas()[myFaces[myFaceI]];
 400                     n /= mag(n) + VSMALL;
 401
 402                     edgeNormal += n;
 403                 }
 404             }
 405
 406             if (debug & 2)
 407             {
 408                 Pout<< " -> real edge hit point:" << edgeHit.rawPoint()
 409                     << " comparing to edge normal:" << edgeNormal
 410                     << endl;
 411             }
 412
 413             // Found face intersection point on this edge. Compare to edge
 414             // normal
 415             return octree<octreeDataFace>::getVolType
 416             (
 417                 edgeNormal,
 418                 sample - curHit.missPoint()
 419             );
 420         }
 421     }
 422
 423
 424     //
 425     // 4] Get the internal edge the face intersection is on
 426     //
 427
 428     forAll(f, fp)
 429     {
 430         pointHit edgeHit =
 431             line<point, const point&>
 432             (
 433                 points[f[fp]],
 434                 mesh_.faceCentres()[faceI]
 435             ).nearestDist(sample);
 436
 437         if ((mag(edgeHit.rawPoint() - curHit.missPoint())/typDim) < tol)
 438         {
 439             // Face intersection point lies on edge between two face triangles
 440
 441             // Calculate edge normal as average of the two triangle normals
 442             const label fpPrev = f.rcIndex(fp);
 443             const label fpNext = f.fcIndex(fp);
 444
 445             vector e = points[f[fp]] - mesh_.faceCentres()[faceI];
 446             vector ePrev = points[f[fpPrev]] - mesh_.faceCentres()[faceI];
 447             vector eNext = points[f[fpNext]] - mesh_.faceCentres()[faceI];
 448
 449             vector nLeft = ePrev ^ e;
 450             nLeft /= mag(nLeft) + VSMALL;
 451
 452             vector nRight = e ^ eNext;
 453             nRight /= mag(nRight) + VSMALL;
 454
 455             if (debug & 2)
 456             {
 457                 Pout<< " -> internal edge hit point:" << edgeHit.rawPoint()
 458                     << " comparing to edge normal "
 459                     << 0.5*(nLeft + nRight)
 460                     << endl;
 461             }
 462
 463             // Found face intersection point on this edge. Compare to edge
 464             // normal
 465             return octree<octreeDataFace>::getVolType
 466             (
 467                 0.5*(nLeft + nRight),
 468                 sample - curHit.missPoint()
 469             );
 470         }
 471     }
 472
 473     if (debug & 2)
 474     {
 475         Pout<< "Did not find sample " << sample
 476             << " anywhere related to nearest face " << faceI << endl
 477             << "Face:";
 478
 479         forAll(f, fp)
 480         {
 481             Pout<< "    vertex:" << f[fp] << "  coord:" << points[f[fp]]
 482                 << endl;
 483         }
 484     }
 485
 486     // Can't determine status of sample with respect to nearest face.
 487     // Either
 488     // - tolerances are wrong. (if e.g. face has zero area)
 489     // - or (more likely) surface is not closed.
 490
 491     return octree<octreeDataFace>::UNKNOWN;
 492 }
 493
 494
 495 bool Foam::octreeDataFace::overlaps
 496 (
 497     const label index,
 498     const treeBoundBox& sampleBb
 499 ) const
 500 {
 501     //return sampleBb.overlaps(allBb_[index]);
 502
 503     //- Exact test of face intersecting bb
 504
 505     // 1. Quick rejection: bb does not intersect face bb at all
 506     if (!sampleBb.overlaps(allBb_[index]))
 507     {
 508         return false;
 509     }
 510
 511     // 2. Check if one or more face points inside
 512     label faceI = meshFaces_[index];
 513
 514     const face& f = mesh_.faces()[faceI];
 515
 516     const pointField& points = mesh_.points();
 517     if (sampleBb.containsAny(points, f))
 518     {
 519         return true;
 520     }
 521
 522     // 3. Difficult case: all points are outside but connecting edges might
 523     // go through cube. Use triangle-bounding box intersection.
 524     const point& fc = mesh_.faceCentres()[faceI];
 525
 526     forAll(f, fp)
 527     {
 528         const label fp1 = f.fcIndex(fp);
 529
 530         bool triIntersects = triangleFuncs::intersectBb
 531         (
 532             points[f[fp]],
 533             points[f[fp1]],
 534             fc,
 535             sampleBb
 536         );
 537
 538         if (triIntersects)
 539         {
 540             return true;
 541         }
 542     }
 543     return false;
 544 }
 545
 546
 547 bool Foam::octreeDataFace::contains(const label, const point&) const
 548 {
 549     notImplemented
 550     (
 551         "octreeDataFace::contains(const label, const point&)"
 552     );
 553     return false;
 554 }
 555
 556
 557 bool Foam::octreeDataFace::intersects
 558 (
 559     const label index,
 560     const point& start,
 561     const point& end,
 562     point& intersectionPoint
 563 ) const
 564 {
 565     label faceI = meshFaces_[index];
 566
 567     const face& f = mesh_.faces()[faceI];
 568
 569     const vector dir(end - start);
 570
 571     // Disable picking up intersections behind us.
 572     scalar oldTol = intersection::setPlanarTol(0.0);
 573
 574     pointHit inter = f.ray
 575     (
 576         start,
 577         dir,
 578         mesh_.points(),
 579         intersection::HALF_RAY,
 580         intersection::VECTOR
 581     );
 582
 583     intersection::setPlanarTol(oldTol);
 584
 585     if (inter.hit() && inter.distance() <= mag(dir))
 586     {
 587         intersectionPoint = inter.hitPoint();
 588
 589         return true;
 590     }
 591     else
 592     {
 593         return false;
 594     }
 595 }
 596
 597
 598 bool Foam::octreeDataFace::findTightest
 599 (
 600     const label index,
 601     const point& sample,
 602     treeBoundBox& tightest
 603 ) const
 604 {
 605     // Get furthest away vertex
 606     point myNear, myFar;
 607     allBb_[index].calcExtremities(sample, myNear, myFar);
 608
 609     const point dist = myFar - sample;
 610     scalar myFarDist = mag(dist);
 611
 612     point tightestNear, tightestFar;
 613     tightest.calcExtremities(sample, tightestNear, tightestFar);
 614
 615     scalar tightestFarDist = mag(tightestFar - sample);
 616
 617     if (tightestFarDist < myFarDist)
 618     {
 619         // Keep current tightest.
 620         return false;
 621     }
 622     else
 623     {
 624         // Construct bb around sample and myFar
 625         const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
 626
 627         tightest.min() = sample - dist2;
 628         tightest.max() = sample + dist2;
 629
 630         return true;
 631     }
 632 }
 633
 634
 635 // Determine numerical value of sign of sample compared to shape at index
 636 Foam::scalar Foam::octreeDataFace::calcSign
 637 (
 638     const label index,
 639     const point& sample,
 640     point& n
 641 ) const
 642 {
 643     label faceI = meshFaces_[index];
 644
 645     n = mesh_.faceAreas()[faceI];
 646
 647     n /= mag(n) + VSMALL;
 648
 649     vector vec = sample - mesh_.faceCentres()[faceI];
 650
 651     vec /= mag(vec) + VSMALL;
 652
 653     return n & vec;
 654 }
 655
 656
 657 // Calculate nearest point on/in shapei
 658 Foam::scalar Foam::octreeDataFace::calcNearest
 659 (
 660     const label index,
 661     const point& sample,
 662     point& nearest
 663 ) const
 664 {
 665     label faceI = meshFaces_[index];
 666
 667     const face& f = mesh_.faces()[faceI];
 668
 669     pointHit nearHit = f.nearestPoint(sample, mesh_.points());
 670
 671     nearest = nearHit.rawPoint();
 672
 673     if (debug & 1)
 674     {
 675         const point& ctr = mesh_.faceCentres()[faceI];
 676
 677         scalar sign = mesh_.faceAreas()[faceI] & (sample - nearest);
 678
 679         Pout<< "octreeDataFace::calcNearest : "
 680             << "sample:" << sample
 681             << "  index:" << index
 682             << "  faceI:" << faceI
 683             << "  ctr:" << ctr
 684             << "  sign:" << sign
 685             << "  nearest point:" << nearest
 686             << "  distance to face:" << nearHit.distance()
 687             << endl;
 688     }
 689     return nearHit.distance();
 690 }
 691
 692
 693 // Calculate nearest point on/in shapei
 694 Foam::scalar Foam::octreeDataFace::calcNearest
 695 (
 696     const label index,
 697     const linePointRef& ln,
 698     point& linePt,
 699     point& shapePt
 700 ) const
 701 {
 702     notImplemented
 703     (
 704         "octreeDataFace::calcNearest(const label, const linePointRef&"
 705         ", point&, point&)"
 706     );
 707     return GREAT;
 708 }
 709
 710
 711 void Foam::octreeDataFace::write(Ostream& os, const label index) const
 712 {
 713     os << meshFaces_[index] << " " << allBb_[index];
 714 }
 715
 716
 717 // ************************************************************************* //