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