src/dynamicMesh/boundaryMesh/octreeDataFaceList.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 Description
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "octreeDataFaceList.H"
  30 #include "octree.H"
  31
  32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  33
  34 defineTypeNameAndDebug(Foam::octreeDataFaceList, 0);
  35
  36 Foam::scalar Foam::octreeDataFaceList::tol = 1E-6;
  37
  38
  39 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  40
  41 void Foam::octreeDataFaceList::calcBb()
  42 {
  43     allBb_.setSize(faceLabels_.size());
  44     allBb_ = treeBoundBox
  45     (
  46         vector(GREAT, GREAT, GREAT),
  47         vector(-GREAT, -GREAT, -GREAT)
  48     );
  49
  50     forAll (faceLabels_, faceLabelI)
  51     {
  52         label faceI = faceLabels_[faceLabelI];
  53
  54         // Update bb of face
  55         treeBoundBox& myBb = allBb_[faceLabelI];
  56
  57         const face& f = mesh_.localFaces()[faceI];
  58
  59         forAll(f, fp)
  60         {
  61             const point& coord = mesh_.localPoints()[f[fp]];
  62
  63             myBb.min() = min(myBb.min(), coord);
  64             myBb.max() = max(myBb.max(), coord);
  65         }
  66     }
  67 }
  68
  69 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  70
  71 // Construct from all faces in bMesh
  72 Foam::octreeDataFaceList::octreeDataFaceList(const bMesh& mesh)
  73 :
  74     mesh_(mesh),
  75     faceLabels_(mesh_.size()),
  76     allBb_(mesh_.size())
  77 {
  78     forAll(faceLabels_, faceI)
  79     {
  80         faceLabels_[faceI] = faceI;
  81     }
  82
  83     // Generate tight fitting bounding box
  84     calcBb();
  85 }
  86
  87
  88 // Construct from selected faces in bMesh
  89 Foam::octreeDataFaceList::octreeDataFaceList
  90 (
  91     const bMesh& mesh,
  92     const labelList& faceLabels
  93 )
  94 :
  95     mesh_(mesh),
  96     faceLabels_(faceLabels),
  97     allBb_(faceLabels.size())
  98 {
  99     // Generate tight fitting bounding box
 100     calcBb();
 101 }
 102
 103
 104
 105 // Construct as copy
 106 Foam::octreeDataFaceList::octreeDataFaceList(const octreeDataFaceList& shapes)
 107 :
 108     mesh_(shapes.mesh()),
 109     faceLabels_(shapes.faceLabels()),
 110     allBb_(shapes.allBb())
 111 {}
 112
 113
 114 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 115
 116 Foam::octreeDataFaceList::~octreeDataFaceList()
 117 {}
 118
 119
 120 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 121
 122
 123 Foam::label Foam::octreeDataFaceList::getSampleType
 124 (
 125     const octree<octreeDataFaceList>& oc,
 126     const point& sample
 127 ) const
 128 {
 129     // Need to determine whether sample is 'inside' or 'outside'
 130     // Done by finding nearest face. This gives back a face which is
 131     // guaranteed to contain nearest point. This point can be
 132     // - in interior of face: compare to face normal
 133     // - on edge of face: compare to edge normal
 134     // - on point of face: compare to point normal
 135     // Unfortunately the octree does not give us back the intersection point
 136     // or where on the face it has hit so we have to recreate all that
 137     // information.
 138
 139
 140     // Find nearest face to sample
 141     treeBoundBox tightest(treeBoundBox::greatBox);
 142
 143     scalar tightestDist = GREAT;
 144
 145     label index = oc.findNearest(sample, tightest, tightestDist);
 146
 147     if (index == -1)
 148     {
 149         FatalErrorIn
 150         (
 151             "octreeDataFaceList::getSampleType"
 152             "(octree<octreeDataFaceList>&, const point&)"
 153         )   << "Could not find " << sample << " in octree."
 154             << abort(FatalError);
 155     }
 156
 157     label faceI = faceLabels_[index];
 158
 159     // Get actual intersection point on face
 160
 161     if (debug & 2)
 162     {
 163         Pout<< "getSampleType : sample:" << sample
 164             << " nearest face:" << faceI;
 165     }
 166
 167     const face& f = mesh_.localFaces()[faceI];
 168
 169     const pointField& points = mesh_.localPoints();
 170
 171     pointHit curHit = f.nearestPoint(sample, points);
 172
 173     //
 174     // 1] Check whether sample is above face
 175     //
 176
 177     if (curHit.hit())
 178     {
 179         // Simple case. Compare to face normal.
 180
 181         if (debug & 2)
 182         {
 183             Pout<< " -> face hit:" << curHit.hitPoint()
 184                 << " comparing to face normal " << mesh_.faceNormals()[faceI]
 185                 << endl;
 186         }
 187         return octree<octreeDataFaceList>::getVolType
 188         (
 189             mesh_.faceNormals()[faceI],
 190             sample - curHit.hitPoint()
 191         );
 192     }
 193
 194     if (debug & 2)
 195     {
 196         Pout<< " -> face miss:" << curHit.missPoint();
 197     }
 198
 199     //
 200     // 2] Check whether intersection is on one of the face vertices or
 201     //    face centre
 202     //
 203
 204     // typical dimension as sqrt of face area.
 205     scalar typDim = sqrt(mag(f.normal(points))) + VSMALL;
 206
 207     forAll(f, fp)
 208     {
 209         if ((mag(points[f[fp]] - curHit.missPoint())/typDim) < tol)
 210         {
 211             // Face intersection point equals face vertex fp
 212
 213             if (debug & 2)
 214             {
 215                     Pout<< " -> face point hit :" << points[f[fp]]
 216                         << " point normal:" << mesh_.pointNormals()[f[fp]]
 217                         << " distance:"
 218                         << mag(points[f[fp]] - curHit.missPoint())/typDim
 219                         << endl;
 220             }
 221             return octree<octreeDataFaceList>::getVolType
 222             (
 223                 mesh_.pointNormals()[f[fp]],
 224                 sample - curHit.missPoint()
 225             );
 226         }
 227     }
 228
 229     // Get face centre
 230     point ctr(f.centre(points));
 231
 232     if ((mag(ctr - curHit.missPoint())/typDim) < tol)
 233     {
 234         // Face intersection point equals face centre. Normal at face centre
 235         // is already average of face normals
 236
 237         if (debug & 2)
 238         {
 239             Pout<< " -> centre hit:" << ctr
 240                 << " distance:"
 241                 << mag(ctr - curHit.missPoint())/typDim
 242                 << endl;
 243         }
 244
 245         return octree<octreeDataFaceList>::getVolType
 246         (
 247             mesh_.faceNormals()[faceI],
 248             sample - curHit.missPoint()
 249         );
 250     }
 251
 252
 253     //
 254     // 3] Get the 'real' edge the face intersection is on
 255     //
 256
 257     const labelList& myEdges = mesh_.faceEdges()[faceI];
 258
 259     forAll(myEdges, myEdgeI)
 260     {
 261         const edge& e = mesh_.edges()[myEdges[myEdgeI]];
 262
 263         pointHit edgeHit =
 264             line<point, const point&>
 265             (
 266                 points[e.start()],
 267                 points[e.end()]
 268             ).nearestDist(sample);
 269
 270         point edgePoint;
 271         if (edgeHit.hit())
 272         {
 273             edgePoint = edgeHit.hitPoint();
 274         }
 275         else
 276         {
 277             edgePoint = edgeHit.missPoint();
 278         }
 279
 280
 281         if ((mag(edgePoint - curHit.missPoint())/typDim) < tol)
 282         {
 283             // Face intersection point lies on edge e
 284
 285             // Calculate edge normal (wrong: uses face normals instead of
 286             // triangle normals)
 287             const labelList& myFaces = mesh_.edgeFaces()[myEdges[myEdgeI]];
 288
 289             vector edgeNormal(vector::zero);
 290
 291             forAll(myFaces, myFaceI)
 292             {
 293                 edgeNormal += mesh_.faceNormals()[myFaces[myFaceI]];
 294             }
 295
 296             if (debug & 2)
 297             {
 298                 Pout<< " -> real edge hit point:" << edgePoint
 299                     << " comparing to edge normal:" << edgeNormal
 300                     << endl;
 301             }
 302
 303             // Found face intersection point on this edge. Compare to edge
 304             // normal
 305             return octree<octreeDataFaceList>::getVolType
 306             (
 307                 edgeNormal,
 308                 sample - curHit.missPoint()
 309             );
 310         }
 311     }
 312
 313
 314     //
 315     // 4] Get the internal edge (vertex - faceCentre) the face intersection
 316     //    is on
 317     //
 318
 319     forAll(f, fp)
 320     {
 321         pointHit edgeHit =
 322             line<point, const point&>
 323             (
 324                 points[f[fp]],
 325                 ctr
 326             ).nearestDist(sample);
 327
 328         point edgePoint;
 329         if (edgeHit.hit())
 330         {
 331             edgePoint = edgeHit.hitPoint();
 332         }
 333         else
 334         {
 335             edgePoint = edgeHit.missPoint();
 336         }
 337
 338         if ((mag(edgePoint - curHit.missPoint())/typDim) < tol)
 339         {
 340             // Face intersection point lies on edge between two face triangles
 341
 342             // Calculate edge normal as average of the two triangle normals
 343             label fpPrev = f.rcIndex(fp);
 344             label fpNext = f.fcIndex(fp);
 345
 346             vector e = points[f[fp]] - ctr;
 347             vector ePrev = points[f[fpPrev]] - ctr;
 348             vector eNext = points[f[fpNext]] - ctr;
 349
 350             vector nLeft = ePrev ^ e;
 351             nLeft /= mag(nLeft) + VSMALL;
 352
 353             vector nRight = e ^ eNext;
 354             nRight /= mag(nRight) + VSMALL;
 355
 356             if (debug & 2)
 357             {
 358                 Pout<< " -> internal edge hit point:" << edgePoint
 359                     << " comparing to edge normal "
 360                     << 0.5*(nLeft + nRight)
 361                     << endl;
 362             }
 363
 364             // Found face intersection point on this edge. Compare to edge
 365             // normal
 366             return octree<octreeDataFaceList>::getVolType
 367             (
 368                 0.5*(nLeft + nRight),
 369                 sample - curHit.missPoint()
 370             );
 371         }
 372     }
 373
 374     if (debug & 2)
 375     {
 376         Pout<< "Did not find sample " << sample
 377             << " anywhere related to nearest face " << faceI << endl
 378             << "Face:";
 379
 380         forAll(f, fp)
 381         {
 382             Pout<< "    vertex:" << f[fp] << "  coord:" << points[f[fp]]
 383                 << endl;
 384         }
 385     }
 386
 387     // Can't determine status of sample with respect to nearest face.
 388     // Either
 389     // - tolerances are wrong. (if e.g. face has zero area)
 390     // - or (more likely) surface is not closed.
 391
 392     return octree<octreeDataFaceList>::UNKNOWN;
 393 }
 394
 395
 396 bool Foam::octreeDataFaceList::overlaps
 397 (
 398     const label index,
 399     const treeBoundBox& sampleBb
 400 ) const
 401 {
 402     return sampleBb.overlaps(allBb_[index]);
 403 }
 404
 405
 406 bool Foam::octreeDataFaceList::contains
 407 (
 408     const label,
 409     const point&
 410 ) const
 411 {
 412     notImplemented
 413     (
 414         "octreeDataFaceList::contains(const label, const point&)"
 415     );
 416     return false;
 417 }
 418
 419
 420 bool Foam::octreeDataFaceList::intersects
 421 (
 422     const label index,
 423     const point& start,
 424     const point& end,
 425     point& intersectionPoint
 426 ) const
 427 {
 428     label faceI = faceLabels_[index];
 429
 430     const face& f = mesh_.localFaces()[faceI];
 431
 432     const vector dir(end - start);
 433
 434     // Disable picking up intersections behind us.
 435     scalar oldTol = intersection::setPlanarTol(0.0);
 436
 437     pointHit inter =
 438         f.ray
 439         (
 440             start,
 441             dir,
 442             mesh_.localPoints(),
 443             intersection::HALF_RAY,
 444             intersection::VECTOR
 445         );
 446
 447     intersection::setPlanarTol(oldTol);
 448
 449     if (inter.hit() && inter.distance() <= mag(dir))
 450     {
 451         intersectionPoint = inter.hitPoint();
 452
 453         return true;
 454     }
 455     else
 456     {
 457         return false;
 458     }
 459 }
 460
 461
 462 bool Foam::octreeDataFaceList::findTightest
 463 (
 464     const label index,
 465     const point& sample,
 466     treeBoundBox& tightest
 467 ) const
 468 {
 469     // Get nearest and furthest away vertex
 470     point myNear, myFar;
 471     allBb_[index].calcExtremities(sample, myNear, myFar);
 472
 473     const point dist = myFar - sample;
 474     scalar myFarDist = mag(dist);
 475
 476     point tightestNear, tightestFar;
 477     tightest.calcExtremities(sample, tightestNear, tightestFar);
 478
 479     scalar tightestFarDist = mag(tightestFar - sample);
 480
 481     if (tightestFarDist < myFarDist)
 482     {
 483         // Keep current tightest.
 484         return false;
 485     }
 486     else
 487     {
 488         // Construct bb around sample and myFar
 489         const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
 490
 491         tightest.min() = sample - dist2;
 492         tightest.max() = sample + dist2;
 493
 494         return true;
 495     }
 496 }
 497
 498
 499 // Determine numerical value of sign of sample compared to shape at index
 500 Foam::scalar Foam::octreeDataFaceList::calcSign
 501 (
 502     const label index,
 503     const point& sample,
 504     vector&
 505 ) const
 506 {
 507     label faceI = faceLabels_[index];
 508
 509     const face& f = mesh_.localFaces()[faceI];
 510
 511     point ctr = f.centre(mesh_.localPoints());
 512
 513     vector vec = sample - ctr;
 514
 515     vec /= mag(vec) + VSMALL;
 516
 517     return mesh_.faceNormals()[faceI] & vec;
 518 }
 519
 520
 521 // Calculate nearest point on/in shapei
 522 Foam::scalar Foam::octreeDataFaceList::calcNearest
 523 (
 524     const label index,
 525     const point& sample,
 526     point& nearest
 527 ) const
 528 {
 529     label faceI = faceLabels_[index];
 530
 531     const face& f = mesh_.localFaces()[faceI];
 532
 533     pointHit nearHit = f.nearestPoint(sample, mesh_.localPoints());
 534
 535     if (nearHit.hit())
 536     {
 537         nearest = nearHit.hitPoint();
 538     }
 539     else
 540     {
 541         nearest = nearHit.missPoint();
 542     }
 543
 544     if (debug & 1)
 545     {
 546         point ctr = f.centre(mesh_.localPoints());
 547
 548         scalar sign = mesh_.faceNormals()[faceI] & (sample - nearest);
 549
 550         Pout<< "octreeDataFaceList::calcNearest : "
 551             << "sample:" << sample
 552             << "  faceI:" << faceI
 553             << "  ctr:" << ctr
 554             << "  sign:" << sign
 555             << "  nearest point:" << nearest
 556             << "  distance to face:" << nearHit.distance()
 557             << endl;
 558     }
 559     return nearHit.distance();
 560 }
 561
 562
 563 void Foam::octreeDataFaceList::write
 564 (
 565     Ostream& os,
 566     const label index
 567 ) const
 568 {
 569     os << faceLabels_[index] << " " << allBb_[index];
 570 }
 571
 572
 573 // ************************************************************************* //