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