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