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