src/meshTools/searchableSurface/searchableCylinder.C

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