src/meshTools/searchableSurface/searchableBox.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 "searchableBox.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "SortableList.H"
  29
  30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34     defineTypeNameAndDebug(searchableBox, 0);
  35     addToRunTimeSelectionTable(searchableSurface, searchableBox, dict);
  36 }
  37
  38
  39 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  40
  41 void Foam::searchableBox::projectOntoCoordPlane
  42 (
  43     const direction dir,
  44     const point& planePt,
  45     pointIndexHit& info
  46 ) const
  47 {
  48     // Set point
  49     info.rawPoint()[dir] = planePt[dir];
  50     // Set face
  51     if (planePt[dir] == min()[dir])
  52     {
  53         info.setIndex(dir*2);
  54     }
  55     else if (planePt[dir] == max()[dir])
  56     {
  57         info.setIndex(dir*2+1);
  58     }
  59     else
  60     {
  61         FatalErrorIn("searchableBox::projectOntoCoordPlane(..)")
  62             << "Point on plane " << planePt
  63             << " is not on coordinate " << min()[dir]
  64             << " nor " << max()[dir] << abort(FatalError);
  65     }
  66 }
  67
  68
  69 // Returns miss or hit with face (0..5) and region(always 0)
  70 Foam::pointIndexHit Foam::searchableBox::findNearest
  71 (
  72     const point& bbMid,
  73     const point& sample,
  74     const scalar nearestDistSqr
  75 ) const
  76 {
  77     // Point can be inside or outside. For every component direction can be
  78     // left of min, right of max or inbetween.
  79     // - outside points: project first one x plane (either min().x()
  80     // or max().x()), then onto y plane and finally z. You should be left
  81     // with intersection point
  82     // - inside point: find nearest side (compare to mid point). Project onto
  83     //   that.
  84
  85     // The face is set to the last projected face.
  86
  87
  88     // Outside point projected onto cube. Assume faces 0..5.
  89     pointIndexHit info(true, sample, -1);
  90     bool outside = false;
  91
  92     // (for internal points) per direction what nearest cube side is
  93     point near;
  94
  95     for (direction dir = 0; dir < vector::nComponents; dir++)
  96     {
  97         if (info.rawPoint()[dir] < min()[dir])
  98         {
  99             projectOntoCoordPlane(dir, min(), info);
 100             outside = true;
 101         }
 102         else if (info.rawPoint()[dir] > max()[dir])
 103         {
 104             projectOntoCoordPlane(dir, max(), info);
 105             outside = true;
 106         }
 107         else if (info.rawPoint()[dir] > bbMid[dir])
 108         {
 109             near[dir] = max()[dir];
 110         }
 111         else
 112         {
 113             near[dir] = min()[dir];
 114         }
 115     }
 116
 117
 118     // For outside points the info will be correct now. Handle inside points
 119     // using the three near distances. Project onto the nearest plane.
 120     if (!outside)
 121     {
 122         vector dist(cmptMag(info.rawPoint() - near));
 123
 124         if (dist.x() < dist.y())
 125         {
 126             if (dist.x() < dist.z())
 127             {
 128                 // Project onto x plane
 129                 projectOntoCoordPlane(vector::X, near, info);
 130             }
 131             else
 132             {
 133                 projectOntoCoordPlane(vector::Z, near, info);
 134             }
 135         }
 136         else
 137         {
 138             if (dist.y() < dist.z())
 139             {
 140                 projectOntoCoordPlane(vector::Y, near, info);
 141             }
 142             else
 143             {
 144                 projectOntoCoordPlane(vector::Z, near, info);
 145             }
 146         }
 147     }
 148
 149
 150     // Check if outside. Optimisation: could do some checks on distance already
 151     // on components above
 152     if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
 153     {
 154         info.setMiss();
 155         info.setIndex(-1);
 156     }
 157
 158     return info;
 159 }
 160
 161
 162 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 163
 164 Foam::searchableBox::searchableBox
 165 (
 166     const IOobject& io,
 167     const treeBoundBox& bb
 168 )
 169 :
 170     searchableSurface(io),
 171     treeBoundBox(bb)
 172 {
 173     if (!contains(midpoint()))
 174     {
 175         FatalErrorIn
 176         (
 177             "Foam::searchableBox::searchableBox\n"
 178             "(\n"
 179             "    const IOobject& io,\n"
 180             "    const treeBoundBox& bb\n"
 181             ")\n"
 182         )   << "Illegal bounding box specification : "
 183             << static_cast<const treeBoundBox>(*this) << exit(FatalError);
 184     }
 185 }
 186
 187
 188 Foam::searchableBox::searchableBox
 189 (
 190     const IOobject& io,
 191     const dictionary& dict
 192 )
 193 :
 194     searchableSurface(io),
 195     treeBoundBox(dict.lookup("min"), dict.lookup("max"))
 196 {
 197     if (!contains(midpoint()))
 198     {
 199         FatalErrorIn
 200         (
 201             "Foam::searchableBox::searchableBox\n"
 202             "(\n"
 203             "    const IOobject& io,\n"
 204             "    const treeBoundBox& bb\n"
 205             ")\n"
 206         )   << "Illegal bounding box specification : "
 207             << static_cast<const treeBoundBox>(*this) << exit(FatalError);
 208     }
 209 }
 210
 211
 212 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 213
 214 Foam::searchableBox::~searchableBox()
 215 {}
 216
 217
 218 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 219
 220 const Foam::wordList& Foam::searchableBox::regions() const
 221 {
 222     if (regions_.empty())
 223     {
 224         regions_.setSize(1);
 225         regions_[0] = "region0";
 226     }
 227     return regions_;
 228 }
 229
 230
 231 Foam::pointField Foam::searchableBox::coordinates() const
 232 {
 233     pointField ctrs(6);
 234
 235     const pointField pts(treeBoundBox::points());
 236     const faceList& fcs = treeBoundBox::faces;
 237
 238     forAll(fcs, i)
 239     {
 240         ctrs[i] = fcs[i].centre(pts);
 241     }
 242     return ctrs;
 243 }
 244
 245
 246 Foam::pointIndexHit Foam::searchableBox::findNearest
 247 (
 248     const point& sample,
 249     const scalar nearestDistSqr
 250 ) const
 251 {
 252     return findNearest(midpoint(), sample, nearestDistSqr);
 253 }
 254
 255
 256 Foam::pointIndexHit Foam::searchableBox::findNearestOnEdge
 257 (
 258     const point& sample,
 259     const scalar nearestDistSqr
 260 ) const
 261 {
 262     const point bbMid(midpoint());
 263
 264     // Outside point projected onto cube. Assume faces 0..5.
 265     pointIndexHit info(true, sample, -1);
 266     bool outside = false;
 267
 268     // (for internal points) per direction what nearest cube side is
 269     point near;
 270
 271     for (direction dir = 0; dir < vector::nComponents; dir++)
 272     {
 273         if (info.rawPoint()[dir] < min()[dir])
 274         {
 275             projectOntoCoordPlane(dir, min(), info);
 276             outside = true;
 277         }
 278         else if (info.rawPoint()[dir] > max()[dir])
 279         {
 280             projectOntoCoordPlane(dir, max(), info);
 281             outside = true;
 282         }
 283         else if (info.rawPoint()[dir] > bbMid[dir])
 284         {
 285             near[dir] = max()[dir];
 286         }
 287         else
 288         {
 289             near[dir] = min()[dir];
 290         }
 291     }
 292
 293
 294     // For outside points the info will be correct now. Handle inside points
 295     // using the three near distances. Project onto the nearest two planes.
 296     if (!outside)
 297     {
 298         // Get the per-component distance to nearest wall
 299         vector dist(cmptMag(info.rawPoint() - near));
 300
 301         SortableList<scalar> sortedDist(3);
 302         sortedDist[0] = dist[0];
 303         sortedDist[1] = dist[1];
 304         sortedDist[2] = dist[2];
 305         sortedDist.sort();
 306
 307         // Project onto nearest
 308         projectOntoCoordPlane(sortedDist.indices()[0], near, info);
 309         // Project onto second nearest
 310         projectOntoCoordPlane(sortedDist.indices()[1], near, info);
 311     }
 312
 313
 314     // Check if outside. Optimisation: could do some checks on distance already
 315     // on components above
 316     if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
 317     {
 318         info.setMiss();
 319         info.setIndex(-1);
 320     }
 321
 322     return info;
 323 }
 324
 325
 326 Foam::pointIndexHit Foam::searchableBox::findNearest
 327 (
 328     const linePointRef& ln,
 329     treeBoundBox& tightest,
 330     point& linePoint
 331 ) const
 332 {
 333     notImplemented
 334     (
 335         "searchableBox::findNearest"
 336         "(const linePointRef&, treeBoundBox&, point&)"
 337     );
 338     return pointIndexHit();
 339 }
 340
 341
 342 Foam::pointIndexHit Foam::searchableBox::findLine
 343 (
 344     const point& start,
 345     const point& end
 346 ) const
 347 {
 348     pointIndexHit info(false, start, -1);
 349
 350     bool foundInter;
 351
 352     if (posBits(start) == 0)
 353     {
 354         if (posBits(end) == 0)
 355         {
 356             // Both start and end inside.
 357             foundInter = false;
 358         }
 359         else
 360         {
 361             // end is outside. Clip to bounding box.
 362             foundInter = intersects(end, start, info.rawPoint());
 363         }
 364     }
 365     else
 366     {
 367         // start is outside. Clip to bounding box.
 368         foundInter = intersects(start, end, info.rawPoint());
 369     }
 370
 371
 372     // Classify point
 373     if (foundInter)
 374     {
 375         info.setHit();
 376
 377         for (direction dir = 0; dir < vector::nComponents; dir++)
 378         {
 379             if (info.rawPoint()[dir] == min()[dir])
 380             {
 381                 info.setIndex(2*dir);
 382                 break;
 383             }
 384             else if (info.rawPoint()[dir] == max()[dir])
 385             {
 386                 info.setIndex(2*dir+1);
 387                 break;
 388             }
 389         }
 390
 391         if (info.index() == -1)
 392         {
 393             FatalErrorIn("searchableBox::findLine(const point&, const point&)")
 394                 << "point " << info.rawPoint()
 395                 << " on segment " << start << end
 396                 << " should be on face of " << *this
 397                 << " but it isn't." << abort(FatalError);
 398         }
 399     }
 400
 401     return info;
 402 }
 403
 404
 405 Foam::pointIndexHit Foam::searchableBox::findLineAny
 406 (
 407     const point& start,
 408     const point& end
 409 ) const
 410 {
 411     return findLine(start, end);
 412 }
 413
 414
 415 void Foam::searchableBox::findNearest
 416 (
 417     const pointField& samples,
 418     const scalarField& nearestDistSqr,
 419     List<pointIndexHit>& info
 420 ) const
 421 {
 422     info.setSize(samples.size());
 423
 424     const point bbMid(midpoint());
 425
 426     forAll(samples, i)
 427     {
 428         info[i] = findNearest(bbMid, samples[i], nearestDistSqr[i]);
 429     }
 430 }
 431
 432
 433 void Foam::searchableBox::findLine
 434 (
 435     const pointField& start,
 436     const pointField& end,
 437     List<pointIndexHit>& info
 438 ) const
 439 {
 440     info.setSize(start.size());
 441
 442     forAll(start, i)
 443     {
 444         info[i] = findLine(start[i], end[i]);
 445     }
 446 }
 447
 448
 449 void Foam::searchableBox::findLineAny
 450 (
 451     const pointField& start,
 452     const pointField& end,
 453     List<pointIndexHit>& info
 454 ) const
 455 {
 456     info.setSize(start.size());
 457
 458     forAll(start, i)
 459     {
 460         info[i] = findLineAny(start[i], end[i]);
 461     }
 462 }
 463
 464
 465 void Foam::searchableBox::findLineAll
 466 (
 467     const pointField& start,
 468     const pointField& end,
 469     List<List<pointIndexHit> >& info
 470 ) const
 471 {
 472     info.setSize(start.size());
 473
 474     // Work array
 475     DynamicList<pointIndexHit, 1, 1> hits;
 476
 477     // Tolerances:
 478     // To find all intersections we add a small vector to the last intersection
 479     // This is chosen such that
 480     // - it is significant (SMALL is smallest representative relative tolerance;
 481     //   we need something bigger since we're doing calculations)
 482     // - if the start-end vector is zero we still progress
 483     const vectorField dirVec(end-start);
 484     const scalarField magSqrDirVec(magSqr(dirVec));
 485     const vectorField smallVec
 486     (
 487         Foam::sqrt(SMALL)*dirVec
 488       + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
 489     );
 490
 491     forAll(start, pointI)
 492     {
 493         // See if any intersection between pt and end
 494         pointIndexHit inter = findLine(start[pointI], end[pointI]);
 495
 496         if (inter.hit())
 497         {
 498             hits.clear();
 499             hits.append(inter);
 500
 501             point pt = inter.hitPoint() + smallVec[pointI];
 502
 503             while (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
 504             {
 505                 // See if any intersection between pt and end
 506                 pointIndexHit inter = findLine(pt, end[pointI]);
 507
 508                 // Check for not hit or hit same face as before (can happen
 509                 // if vector along surface of face)
 510                 if
 511                 (
 512                     !inter.hit()
 513                  || (inter.index() == hits.last().index())
 514                 )
 515                 {
 516                     break;
 517                 }
 518                 hits.append(inter);
 519
 520                 pt = inter.hitPoint() + smallVec[pointI];
 521             }
 522
 523             info[pointI].transfer(hits);
 524         }
 525         else
 526         {
 527             info[pointI].clear();
 528         }
 529     }
 530 }
 531
 532
 533 void Foam::searchableBox::getRegion
 534 (
 535     const List<pointIndexHit>& info,
 536     labelList& region
 537 ) const
 538 {
 539     region.setSize(info.size());
 540     region = 0;
 541 }
 542
 543
 544 void Foam::searchableBox::getNormal
 545 (
 546     const List<pointIndexHit>& info,
 547     vectorField& normal
 548 ) const
 549 {
 550     normal.setSize(info.size());
 551     normal = vector::zero;
 552
 553     forAll(info, i)
 554     {
 555         if (info[i].hit())
 556         {
 557             normal[i] = treeBoundBox::faceNormals[info[i].index()];
 558         }
 559         else
 560         {
 561             // Set to what?
 562         }
 563     }
 564 }
 565
 566
 567 void Foam::searchableBox::getVolumeType
 568 (
 569     const pointField& points,
 570     List<volumeType>& volType
 571 ) const
 572 {
 573     volType.setSize(points.size());
 574     volType = INSIDE;
 575
 576     forAll(points, pointI)
 577     {
 578         const point& pt = points[pointI];
 579
 580         for (direction dir = 0; dir < vector::nComponents; dir++)
 581         {
 582             if (pt[dir] < min()[dir] || pt[dir] > max()[dir])
 583             {
 584                 volType[pointI] = OUTSIDE;
 585                 break;
 586             }
 587         }
 588     }
 589 }
 590
 591
 592 // ************************************************************************* //