src/meshTools/searchableSurface/searchableBox.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 "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::pointIndexHit Foam::searchableBox::findNearest
 233 (
 234     const point& sample,
 235     const scalar nearestDistSqr
 236 ) const
 237 {
 238     return findNearest(midpoint(), sample, nearestDistSqr);
 239 }
 240
 241
 242 Foam::pointIndexHit Foam::searchableBox::findNearestOnEdge
 243 (
 244     const point& sample,
 245     const scalar nearestDistSqr
 246 ) const
 247 {
 248     const point bbMid(midpoint());
 249
 250     // Outside point projected onto cube. Assume faces 0..5.
 251     pointIndexHit info(true, sample, -1);
 252     bool outside = false;
 253
 254     // (for internal points) per direction what nearest cube side is
 255     point near;
 256
 257     for (direction dir = 0; dir < vector::nComponents; dir++)
 258     {
 259         if (info.rawPoint()[dir] < min()[dir])
 260         {
 261             projectOntoCoordPlane(dir, min(), info);
 262             outside = true;
 263         }
 264         else if (info.rawPoint()[dir] > max()[dir])
 265         {
 266             projectOntoCoordPlane(dir, max(), info);
 267             outside = true;
 268         }
 269         else if (info.rawPoint()[dir] > bbMid[dir])
 270         {
 271             near[dir] = max()[dir];
 272         }
 273         else
 274         {
 275             near[dir] = min()[dir];
 276         }
 277     }
 278
 279
 280     // For outside points the info will be correct now. Handle inside points
 281     // using the three near distances. Project onto the nearest two planes.
 282     if (!outside)
 283     {
 284         // Get the per-component distance to nearest wall
 285         vector dist(cmptMag(info.rawPoint() - near));
 286
 287         SortableList<scalar> sortedDist(3);
 288         sortedDist[0] = dist[0];
 289         sortedDist[1] = dist[1];
 290         sortedDist[2] = dist[2];
 291         sortedDist.sort();
 292
 293         // Project onto nearest
 294         projectOntoCoordPlane(sortedDist.indices()[0], near, info);
 295         // Project onto second nearest
 296         projectOntoCoordPlane(sortedDist.indices()[1], near, info);
 297     }
 298
 299
 300     // Check if outside. Optimisation: could do some checks on distance already
 301     // on components above
 302     if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
 303     {
 304         info.setMiss();
 305         info.setIndex(-1);
 306     }
 307
 308     return info;
 309 }
 310
 311
 312 Foam::pointIndexHit Foam::searchableBox::findNearest
 313 (
 314     const linePointRef& ln,
 315     treeBoundBox& tightest,
 316     point& linePoint
 317 ) const
 318 {
 319     notImplemented
 320     (
 321         "searchableBox::findNearest"
 322         "(const linePointRef&, treeBoundBox&, point&)"
 323     );
 324     return pointIndexHit();
 325 }
 326
 327
 328 Foam::pointIndexHit Foam::searchableBox::findLine
 329 (
 330     const point& start,
 331     const point& end
 332 ) const
 333 {
 334     pointIndexHit info(false, start, -1);
 335
 336     bool foundInter;
 337
 338     if (posBits(start) == 0)
 339     {
 340         if (posBits(end) == 0)
 341         {
 342             // Both start and end inside.
 343             foundInter = false;
 344         }
 345         else
 346         {
 347             // end is outside. Clip to bounding box.
 348             foundInter = intersects(end, start, info.rawPoint());
 349         }
 350     }
 351     else
 352     {
 353         // start is outside. Clip to bounding box.
 354         foundInter = intersects(start, end, info.rawPoint());
 355     }
 356
 357
 358     // Classify point
 359     if (foundInter)
 360     {
 361         info.setHit();
 362
 363         for (direction dir = 0; dir < vector::nComponents; dir++)
 364         {
 365             if (info.rawPoint()[dir] == min()[dir])
 366             {
 367                 info.setIndex(2*dir);
 368                 break;
 369             }
 370             else if (info.rawPoint()[dir] == max()[dir])
 371             {
 372                 info.setIndex(2*dir+1);
 373                 break;
 374             }
 375         }
 376
 377         if (info.index() == -1)
 378         {
 379             FatalErrorIn("searchableBox::findLine(const point&, const point&)")
 380                 << "point " << info.rawPoint()
 381                 << " on segment " << start << end
 382                 << " should be on face of " << *this
 383                 << " but it isn't." << abort(FatalError);
 384         }
 385     }
 386
 387     return info;
 388 }
 389
 390
 391 Foam::pointIndexHit Foam::searchableBox::findLineAny
 392 (
 393     const point& start,
 394     const point& end
 395 ) const
 396 {
 397     return findLine(start, end);
 398 }
 399
 400
 401 void Foam::searchableBox::findNearest
 402 (
 403     const pointField& samples,
 404     const scalarField& nearestDistSqr,
 405     List<pointIndexHit>& info
 406 ) const
 407 {
 408     info.setSize(samples.size());
 409
 410     const point bbMid(midpoint());
 411
 412     forAll(samples, i)
 413     {
 414         info[i] = findNearest(bbMid, samples[i], nearestDistSqr[i]);
 415     }
 416 }
 417
 418
 419 void Foam::searchableBox::findLine
 420 (
 421     const pointField& start,
 422     const pointField& end,
 423     List<pointIndexHit>& info
 424 ) const
 425 {
 426     info.setSize(start.size());
 427
 428     forAll(start, i)
 429     {
 430         info[i] = findLine(start[i], end[i]);
 431     }
 432 }
 433
 434
 435 void Foam::searchableBox::findLineAny
 436 (
 437     const pointField& start,
 438     const pointField& end,
 439     List<pointIndexHit>& info
 440 ) const
 441 {
 442     info.setSize(start.size());
 443
 444     forAll(start, i)
 445     {
 446         info[i] = findLineAny(start[i], end[i]);
 447     }
 448 }
 449
 450
 451 void Foam::searchableBox::findLineAll
 452 (
 453     const pointField& start,
 454     const pointField& end,
 455     List<List<pointIndexHit> >& info
 456 ) const
 457 {
 458     info.setSize(start.size());
 459
 460     // Work array
 461     DynamicList<pointIndexHit, 1, 1> hits;
 462
 463     // Tolerances:
 464     // To find all intersections we add a small vector to the last intersection
 465     // This is chosen such that
 466     // - it is significant (SMALL is smallest representative relative tolerance;
 467     //   we need something bigger since we're doing calculations)
 468     // - if the start-end vector is zero we still progress
 469     const vectorField dirVec(end-start);
 470     const scalarField magSqrDirVec(magSqr(dirVec));
 471     const vectorField smallVec
 472     (
 473         Foam::sqrt(SMALL)*dirVec
 474       + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
 475     );
 476
 477     forAll(start, pointI)
 478     {
 479         // See if any intersection between pt and end
 480         pointIndexHit inter = findLine(start[pointI], end[pointI]);
 481
 482         if (inter.hit())
 483         {
 484             hits.clear();
 485             hits.append(inter);
 486
 487             point pt = inter.hitPoint() + smallVec[pointI];
 488
 489             while (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
 490             {
 491                 // See if any intersection between pt and end
 492                 pointIndexHit inter = findLine(pt, end[pointI]);
 493
 494                 // Check for not hit or hit same face as before (can happen
 495                 // if vector along surface of face)
 496                 if
 497                 (
 498                     !inter.hit()
 499                  || (inter.index() == hits[hits.size()-1].index())
 500                 )
 501                 {
 502                     break;
 503                 }
 504                 hits.append(inter);
 505
 506                 pt = inter.hitPoint() + smallVec[pointI];
 507             }
 508
 509             info[pointI].transfer(hits);
 510         }
 511         else
 512         {
 513             info[pointI].clear();
 514         }
 515     }
 516 }
 517
 518
 519 void Foam::searchableBox::getRegion
 520 (
 521     const List<pointIndexHit>& info,
 522     labelList& region
 523 ) const
 524 {
 525     region.setSize(info.size());
 526     region = 0;
 527 }
 528
 529
 530 void Foam::searchableBox::getNormal
 531 (
 532     const List<pointIndexHit>& info,
 533     vectorField& normal
 534 ) const
 535 {
 536     normal.setSize(info.size());
 537     normal = vector::zero;
 538
 539     forAll(info, i)
 540     {
 541         if (info[i].hit())
 542         {
 543             normal[i] = treeBoundBox::faceNormals[info[i].index()];
 544         }
 545         else
 546         {
 547             // Set to what?
 548         }
 549     }
 550 }
 551
 552
 553 void Foam::searchableBox::getVolumeType
 554 (
 555     const pointField& points,
 556     List<volumeType>& volType
 557 ) const
 558 {
 559     volType.setSize(points.size());
 560     volType = INSIDE;
 561
 562     forAll(points, pointI)
 563     {
 564         const point& pt = points[pointI];
 565
 566         for (direction dir = 0; dir < vector::nComponents; dir++)
 567         {
 568             if (pt[dir] < min()[dir] || pt[dir] > max()[dir])
 569             {
 570                 volType[pointI] = OUTSIDE;
 571                 break;
 572             }
 573         }
 574     }
 575 }
 576
 577
 578 // ************************************************************************* //