src/meshTools/searchableSurface/searchableCylinder.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "searchableCylinder.H"
  27 #include "addToRunTimeSelectionTable.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33
  34 defineTypeNameAndDebug(searchableCylinder, 0);
  35 addToRunTimeSelectionTable(searchableSurface, searchableCylinder, dict);
  36
  37 }
  38
  39
  40 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  41
  42 Foam::pointField Foam::searchableCylinder::coordinates() const
  43 {
  44     pointField ctrs(1, 0.5*(point1_ + point2_));
  45
  46     return ctrs;
  47 }
  48
  49
  50 Foam::pointIndexHit Foam::searchableCylinder::findNearest
  51 (
  52     const point& sample,
  53     const scalar nearestDistSqr
  54 ) const
  55 {
  56     pointIndexHit info(false, sample, -1);
  57
  58     vector v(sample - point1_);
  59
  60     // Decompose sample-point1 into normal and parallel component
  61     scalar parallel = (v & unitDir_);
  62
  63     // Remove the parallel component
  64     v -= parallel*unitDir_;
  65     scalar magV = mag(v);
  66
  67     if (parallel <= 0)
  68     {
  69         // nearest is at point1 end cap. Clip to radius.
  70         if (magV < ROOTVSMALL)
  71         {
  72             info.setPoint(point1_);
  73         }
  74         else
  75         {
  76             //info.setPoint(point1_ + min(magV, radius_)*v/magV);
  77             info.setPoint(point1_ + radius_*(v/magV));
  78         }
  79     }
  80     else if (parallel >= magDir_)
  81     {
  82         // nearest is at point2
  83         if (magV < ROOTVSMALL)
  84         {
  85             info.setPoint(point2_);
  86         }
  87         else
  88         {
  89             info.setPoint(point2_ + min(magV, radius_)*v/magV);
  90         }
  91     }
  92     else
  93     {
  94         if (magV < ROOTVSMALL)
  95         {
  96             info.setPoint(point1_ + parallel*unitDir_);
  97         }
  98         else
  99         {
 100             // Project onto radius
 101             info.setPoint(point1_ + parallel*unitDir_ + radius_*v/magV);
 102         }
 103     }
 104
 105     if (magSqr(sample - info.rawPoint()) < nearestDistSqr)
 106     {
 107         info.setHit();
 108         info.setIndex(0);
 109     }
 110
 111     return info;
 112 }
 113
 114
 115 Foam::scalar Foam::searchableCylinder::radius2(const point& pt) const
 116 {
 117     const vector x = (pt-point1_) ^ unitDir_;
 118     return x&x;
 119 }
 120
 121
 122 // From http://www.gamedev.net/community/forums/topic.asp?topic_id=467789 -
 123 // intersection of cylinder with ray
 124 void Foam::searchableCylinder::findLineAll
 125 (
 126     const point& start,
 127     const point& end,
 128     pointIndexHit& near,
 129     pointIndexHit& far
 130 ) const
 131 {
 132     near.setMiss();
 133     far.setMiss();
 134
 135     vector point1Start(start-point1_);
 136     vector point2Start(start-point2_);
 137     vector point1End(end-point1_);
 138
 139     // Quick rejection of complete vector outside endcaps
 140     scalar s1 = point1Start&unitDir_;
 141     scalar s2 = point1End&unitDir_;
 142
 143     if ((s1 < 0 && s2 < 0) || (s1 > magDir_ && s2 > magDir_))
 144     {
 145         return;
 146     }
 147
 148     // Line as P = start+t*V  where V is unit vector and t=[0..mag(end-start)]
 149     vector V(end-start);
 150     scalar magV = mag(V);
 151     if (magV < ROOTVSMALL)
 152     {
 153         return;
 154     }
 155     V /= magV;
 156
 157
 158     // We now get the nearest intersections to start. This can either be
 159     // the intersection with the end plane or with the cylinder side.
 160
 161     // Get the two points (expressed in t) on the end planes. This is to
 162     // clip any cylinder intersection against.
 163     scalar tPoint1;
 164     scalar tPoint2;
 165
 166     // Maintain the two intersections with the endcaps
 167     scalar tNear = VGREAT;
 168     scalar tFar = VGREAT;
 169
 170     {
 171         scalar s = (V&unitDir_);
 172         if (mag(s) > VSMALL)
 173         {
 174             tPoint1 = -s1/s;
 175             tPoint2 = -(point2Start&unitDir_)/s;
 176             if (tPoint2 < tPoint1)
 177             {
 178                 Swap(tPoint1, tPoint2);
 179             }
 180             if (tPoint1 > magV || tPoint2 < 0)
 181             {
 182                 return;
 183             }
 184
 185             // See if the points on the endcaps are actually inside the cylinder
 186             if (tPoint1 >= 0 && tPoint1 <= magV)
 187             {
 188                 if (radius2(start+tPoint1*V) <= sqr(radius_))
 189                 {
 190                     tNear = tPoint1;
 191                 }
 192             }
 193             if (tPoint2 >= 0 && tPoint2 <= magV)
 194             {
 195                 if (radius2(start+tPoint2*V) <= sqr(radius_))
 196                 {
 197                     // Check if already have a near hit from point1
 198                     if (tNear <= 1)
 199                     {
 200                         tFar = tPoint2;
 201                     }
 202                     else
 203                     {
 204                         tNear = tPoint2;
 205                     }
 206                 }
 207             }
 208         }
 209         else
 210         {
 211             // Vector perpendicular to cylinder. Check for outside already done
 212             // above so just set tpoint to allow all.
 213             tPoint1 = -VGREAT;
 214             tPoint2 = VGREAT;
 215         }
 216     }
 217
 218
 219     const vector x = point1Start ^ unitDir_;
 220     const vector y = V ^ unitDir_;
 221     const scalar d = sqr(radius_);
 222
 223     // Second order equation of the form a*t^2 + b*t + c
 224     const scalar a = (y&y);
 225     const scalar b = 2*(x&y);
 226     const scalar c = (x&x)-d;
 227
 228     const scalar disc = b*b-4*a*c;
 229
 230     scalar t1 = -VGREAT;
 231     scalar t2 = VGREAT;
 232
 233     if (disc < 0)
 234     {
 235         // Fully outside
 236         return;
 237     }
 238     else if (disc < ROOTVSMALL)
 239     {
 240         // Single solution
 241         if (mag(a) > ROOTVSMALL)
 242         {
 243             t1 = -b/(2*a);
 244
 245             //Pout<< "single solution t:" << t1
 246             //    << " for start:" << start << " end:" << end
 247             //    << " c:" << c << endl;
 248
 249             if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
 250             {
 251                 // valid. Insert sorted.
 252                 if (t1 < tNear)
 253                 {
 254                     tFar = tNear;
 255                     tNear = t1;
 256                 }
 257                 else if (t1 < tFar)
 258                 {
 259                     tFar = t1;
 260                 }
 261             }
 262             else
 263             {
 264                 return;
 265             }
 266         }
 267         else
 268         {
 269             // Aligned with axis. Check if outside radius
 270             //Pout<< "small discriminant:" << disc
 271             //    << " for start:" << start << " end:" << end
 272             //    << " magV:" << magV
 273             //    << " c:" << c << endl;
 274             if (c > 0)
 275             {
 276                 return;
 277             }
 278         }
 279     }
 280     else
 281     {
 282         if (mag(a) > ROOTVSMALL)
 283         {
 284             scalar sqrtDisc = sqrt(disc);
 285
 286             t1 = (-b - sqrtDisc)/(2*a);
 287             t2 = (-b + sqrtDisc)/(2*a);
 288             if (t2 < t1)
 289             {
 290                 Swap(t1, t2);
 291             }
 292
 293             if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
 294             {
 295                 // valid. Insert sorted.
 296                 if (t1 < tNear)
 297                 {
 298                     tFar = tNear;
 299                     tNear = t1;
 300                 }
 301                 else if (t1 < tFar)
 302                 {
 303                     tFar = t1;
 304                 }
 305             }
 306             if (t2 >= 0 && t2 <= magV && t2 >= tPoint1 && t2 <= tPoint2)
 307             {
 308                 // valid. Insert sorted.
 309                 if (t2 < tNear)
 310                 {
 311                     tFar = tNear;
 312                     tNear = t2;
 313                 }
 314                 else if (t2 < tFar)
 315                 {
 316                     tFar = t2;
 317                 }
 318             }
 319             //Pout<< "two solutions t1:" << t1 << " t2:" << t2
 320             //    << " for start:" << start << " end:" << end
 321             //    << " magV:" << magV
 322             //    << " c:" << c << endl;
 323         }
 324         else
 325         {
 326             // Aligned with axis. Check if outside radius
 327             //Pout<< "large discriminant:" << disc
 328             //    << " small a:" << a
 329             //    << " for start:" << start << " end:" << end
 330             //    << " magV:" << magV
 331             //    << " c:" << c << endl;
 332             if (c > 0)
 333             {
 334                 return;
 335             }
 336         }
 337     }
 338
 339     // Check tNear, tFar
 340     if (tNear >= 0 && tNear <= magV)
 341     {
 342         near.setPoint(start+tNear*V);
 343         near.setHit();
 344         near.setIndex(0);
 345
 346         if (tFar <= magV)
 347         {
 348             far.setPoint(start+tFar*V);
 349             far.setHit();
 350             far.setIndex(0);
 351         }
 352     }
 353     else if (tFar >= 0 && tFar <= magV)
 354     {
 355         near.setPoint(start+tFar*V);
 356         near.setHit();
 357         near.setIndex(0);
 358     }
 359 }
 360
 361
 362 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 363
 364 Foam::searchableCylinder::searchableCylinder
 365 (
 366     const IOobject& io,
 367     const point& point1,
 368     const point& point2,
 369     const scalar radius
 370 )
 371 :
 372     searchableSurface(io),
 373     point1_(point1),
 374     point2_(point2),
 375     magDir_(mag(point2_-point1_)),
 376     unitDir_((point2_-point1_)/magDir_),
 377     radius_(radius)
 378 {}
 379
 380
 381 Foam::searchableCylinder::searchableCylinder
 382 (
 383     const IOobject& io,
 384     const dictionary& dict
 385 )
 386 :
 387     searchableSurface(io),
 388     point1_(dict.lookup("point1")),
 389     point2_(dict.lookup("point2")),
 390     magDir_(mag(point2_-point1_)),
 391     unitDir_((point2_-point1_)/magDir_),
 392     radius_(readScalar(dict.lookup("radius")))
 393 {}
 394
 395
 396 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 397
 398 Foam::searchableCylinder::~searchableCylinder()
 399 {}
 400
 401
 402 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 403
 404 const Foam::wordList& Foam::searchableCylinder::regions() const
 405 {
 406     if (regions_.empty())
 407     {
 408         regions_.setSize(1);
 409         regions_[0] = "region0";
 410     }
 411     return regions_;
 412 }
 413
 414
 415 void Foam::searchableCylinder::findNearest
 416 (
 417     const pointField& samples,
 418     const scalarField& nearestDistSqr,
 419     List<pointIndexHit>& info
 420 ) const
 421 {
 422     info.setSize(samples.size());
 423
 424     forAll(samples, i)
 425     {
 426         info[i] = findNearest(samples[i], nearestDistSqr[i]);
 427     }
 428 }
 429
 430
 431 void Foam::searchableCylinder::findLine
 432 (
 433     const pointField& start,
 434     const pointField& end,
 435     List<pointIndexHit>& info
 436 ) const
 437 {
 438     info.setSize(start.size());
 439
 440     forAll(start, i)
 441     {
 442         // Pick nearest intersection. If none intersected take second one.
 443         pointIndexHit b;
 444         findLineAll(start[i], end[i], info[i], b);
 445         if (!info[i].hit() && b.hit())
 446         {
 447             info[i] = b;
 448         }
 449     }
 450 }
 451
 452
 453 void Foam::searchableCylinder::findLineAny
 454 (
 455     const pointField& start,
 456     const pointField& end,
 457     List<pointIndexHit>& info
 458 ) const
 459 {
 460     info.setSize(start.size());
 461
 462     forAll(start, i)
 463     {
 464         // Discard far intersection
 465         pointIndexHit b;
 466         findLineAll(start[i], end[i], info[i], b);
 467         if (!info[i].hit() && b.hit())
 468         {
 469             info[i] = b;
 470         }
 471     }
 472 }
 473
 474
 475 void Foam::searchableCylinder::findLineAll
 476 (
 477     const pointField& start,
 478     const pointField& end,
 479     List<List<pointIndexHit> >& info
 480 ) const
 481 {
 482     info.setSize(start.size());
 483
 484     forAll(start, i)
 485     {
 486         pointIndexHit near, far;
 487         findLineAll(start[i], end[i], near, far);
 488
 489         if (near.hit())
 490         {
 491             if (far.hit())
 492             {
 493                 info[i].setSize(2);
 494                 info[i][0] = near;
 495                 info[i][1] = far;
 496             }
 497             else
 498             {
 499                 info[i].setSize(1);
 500                 info[i][0] = near;
 501             }
 502         }
 503         else
 504         {
 505             if (far.hit())
 506             {
 507                 info[i].setSize(1);
 508                 info[i][0] = far;
 509             }
 510             else
 511             {
 512                 info[i].clear();
 513             }
 514         }
 515     }
 516 }
 517
 518
 519 void Foam::searchableCylinder::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::searchableCylinder::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             vector v(info[i].hitPoint() - point1_);
 544
 545             // Decompose sample-point1 into normal and parallel component
 546             scalar parallel = v & unitDir_;
 547
 548             if (parallel < 0)
 549             {
 550                 normal[i] = -unitDir_;
 551             }
 552             else if (parallel > magDir_)
 553             {
 554                 normal[i] = -unitDir_;
 555             }
 556             else
 557             {
 558                 // Remove the parallel component
 559                 v -= parallel*unitDir_;
 560                 normal[i] = v/mag(v);
 561             }
 562         }
 563     }
 564 }
 565
 566
 567 void Foam::searchableCylinder::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         vector v(pt - point1_);
 581
 582         // Decompose sample-point1 into normal and parallel component
 583         scalar parallel = v & unitDir_;
 584
 585         if (parallel < 0)
 586         {
 587             // left of point1 endcap
 588             volType[pointI] = OUTSIDE;
 589         }
 590         else if (parallel > magDir_)
 591         {
 592             // right of point2 endcap
 593             volType[pointI] = OUTSIDE;
 594         }
 595         else
 596         {
 597             // Remove the parallel component
 598             v -= parallel*unitDir_;
 599
 600             if (mag(v) > radius_)
 601             {
 602                 volType[pointI] = OUTSIDE;
 603             }
 604             else
 605             {
 606                 volType[pointI] = INSIDE;
 607             }
 608         }
 609     }
 610 }
 611
 612
 613 // ************************************************************************* //