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