src/meshTools/searchableSurface/searchableBox.C

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