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