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