src/meshTools/indexedOctree/treeDataFace.C

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