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