FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

searchableCylinder.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2011 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "searchableCylinder.H"
00027 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00028 
00029 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00030 
00031 namespace Foam
00032 {
00033 
00034 defineTypeNameAndDebug(searchableCylinder, 0);
00035 addToRunTimeSelectionTable(searchableSurface, searchableCylinder, dict);
00036 
00037 }
00038 
00039 
00040 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00041 
00042 Foam::pointField Foam::searchableCylinder::coordinates() const
00043 {
00044     pointField ctrs(1, 0.5*(point1_ + point2_));
00045 
00046     return ctrs;
00047 }
00048 
00049 
00050 Foam::pointIndexHit Foam::searchableCylinder::findNearest
00051 (
00052     const point& sample,
00053     const scalar nearestDistSqr
00054 ) const
00055 {
00056     pointIndexHit info(false, sample, -1);
00057 
00058     vector v(sample - point1_);
00059 
00060     // Decompose sample-point1 into normal and parallel component
00061     scalar parallel = (v & unitDir_);
00062 
00063     // Remove the parallel component and normalise
00064     v -= parallel*unitDir_;
00065     scalar magV = mag(v);
00066 
00067     if (magV < ROOTVSMALL)
00068     {
00069         v = vector::zero;
00070     }
00071     else
00072     {
00073         v /= magV;
00074     }
00075 
00076     if (parallel <= 0)
00077     {
00078         // nearest is at point1 end cap. Clip to radius.
00079         info.setPoint(point1_ + min(magV, radius_)*v);
00080     }
00081     else if (parallel >= magDir_)
00082     {
00083         // nearest is at point2 end cap. Clip to radius.
00084         info.setPoint(point2_ + min(magV, radius_)*v);
00085     }
00086     else
00087     {
00088         // inbetween endcaps. Might either be nearer endcaps or cylinder wall
00089 
00090         // distance to endpoint: parallel or parallel-magDir
00091         // distance to cylinder wall: magV-radius_
00092 
00093         // Nearest cylinder point
00094         point cylPt;
00095         if (magV < ROOTVSMALL)
00096         {
00097             // Point exactly on centre line. Take any point on wall.
00098             vector e1 = point(1,0,0) ^ unitDir_;
00099             scalar magE1 = mag(e1);
00100             if (magE1 < SMALL)
00101             {
00102                 e1 = point(0,1,0) ^ unitDir_;
00103                 magE1 = mag(e1);
00104             }
00105             e1 /= magE1;
00106             cylPt = sample + radius_*e1;
00107         }
00108         else
00109         {
00110             cylPt = sample + (radius_-magV)*v;
00111         }
00112 
00113         if (parallel < 0.5*magDir_)
00114         {
00115             // Project onto p1 endcap
00116             point end1Pt = point1_ + min(magV, radius_)*v;
00117 
00118             if (magSqr(sample-cylPt) < magSqr(sample-end1Pt))
00119             {
00120                 info.setPoint(cylPt);
00121             }
00122             else
00123             {
00124                 info.setPoint(end1Pt);
00125             }
00126         }
00127         else
00128         {
00129             // Project onto p2 endcap
00130             point end2Pt = point2_ + min(magV, radius_)*v;
00131 
00132             if (magSqr(sample-cylPt) < magSqr(sample-end2Pt))
00133             {
00134                 info.setPoint(cylPt);
00135             }
00136             else
00137             {
00138                 info.setPoint(end2Pt);
00139             }
00140         }
00141     }
00142 
00143     if (magSqr(sample - info.rawPoint()) < nearestDistSqr)
00144     {
00145         info.setHit();
00146         info.setIndex(0);
00147     }
00148 
00149     return info;
00150 }
00151 
00152 
00153 Foam::scalar Foam::searchableCylinder::radius2(const point& pt) const
00154 {
00155     const vector x = (pt-point1_) ^ unitDir_;
00156     return x&x;
00157 }
00158 
00159 
00160 // From http://www.gamedev.net/community/forums/topic.asp?topic_id=467789 -
00161 // intersection of cylinder with ray
00162 void Foam::searchableCylinder::findLineAll
00163 (
00164     const point& start,
00165     const point& end,
00166     pointIndexHit& near,
00167     pointIndexHit& far
00168 ) const
00169 {
00170     near.setMiss();
00171     far.setMiss();
00172 
00173     vector point1Start(start-point1_);
00174     vector point2Start(start-point2_);
00175     vector point1End(end-point1_);
00176 
00177     // Quick rejection of complete vector outside endcaps
00178     scalar s1 = point1Start&unitDir_;
00179     scalar s2 = point1End&unitDir_;
00180 
00181     if ((s1 < 0 && s2 < 0) || (s1 > magDir_ && s2 > magDir_))
00182     {
00183         return;
00184     }
00185 
00186     // Line as P = start+t*V  where V is unit vector and t=[0..mag(end-start)]
00187     vector V(end-start);
00188     scalar magV = mag(V);
00189     if (magV < ROOTVSMALL)
00190     {
00191         return;
00192     }
00193     V /= magV;
00194 
00195 
00196     // We now get the nearest intersections to start. This can either be
00197     // the intersection with the end plane or with the cylinder side.
00198 
00199     // Get the two points (expressed in t) on the end planes. This is to
00200     // clip any cylinder intersection against.
00201     scalar tPoint1;
00202     scalar tPoint2;
00203 
00204     // Maintain the two intersections with the endcaps
00205     scalar tNear = VGREAT;
00206     scalar tFar = VGREAT;
00207 
00208     {
00209         scalar s = (V&unitDir_);
00210         if (mag(s) > VSMALL)
00211         {
00212             tPoint1 = -s1/s;
00213             tPoint2 = -(point2Start&unitDir_)/s;
00214             if (tPoint2 < tPoint1)
00215             {
00216                 Swap(tPoint1, tPoint2);
00217             }
00218             if (tPoint1 > magV || tPoint2 < 0)
00219             {
00220                 return;
00221             }
00222 
00223             // See if the points on the endcaps are actually inside the cylinder
00224             if (tPoint1 >= 0 && tPoint1 <= magV)
00225             {
00226                 if (radius2(start+tPoint1*V) <= sqr(radius_))
00227                 {
00228                     tNear = tPoint1;
00229                 }
00230             }
00231             if (tPoint2 >= 0 && tPoint2 <= magV)
00232             {
00233                 if (radius2(start+tPoint2*V) <= sqr(radius_))
00234                 {
00235                     // Check if already have a near hit from point1
00236                     if (tNear <= 1)
00237                     {
00238                         tFar = tPoint2;
00239                     }
00240                     else
00241                     {
00242                         tNear = tPoint2;
00243                     }
00244                 }
00245             }
00246         }
00247         else
00248         {
00249             // Vector perpendicular to cylinder. Check for outside already done
00250             // above so just set tpoint to allow all.
00251             tPoint1 = -VGREAT;
00252             tPoint2 = VGREAT;
00253         }
00254     }
00255 
00256 
00257     const vector x = point1Start ^ unitDir_;
00258     const vector y = V ^ unitDir_;
00259     const scalar d = sqr(radius_);
00260 
00261     // Second order equation of the form a*t^2 + b*t + c
00262     const scalar a = (y&y);
00263     const scalar b = 2*(x&y);
00264     const scalar c = (x&x)-d;
00265 
00266     const scalar disc = b*b-4*a*c;
00267 
00268     scalar t1 = -VGREAT;
00269     scalar t2 = VGREAT;
00270 
00271     if (disc < 0)
00272     {
00273         // Fully outside
00274         return;
00275     }
00276     else if (disc < ROOTVSMALL)
00277     {
00278         // Single solution
00279         if (mag(a) > ROOTVSMALL)
00280         {
00281             t1 = -b/(2*a);
00282 
00283             //Pout<< "single solution t:" << t1
00284             //    << " for start:" << start << " end:" << end
00285             //    << " c:" << c << endl;
00286 
00287             if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
00288             {
00289                 // valid. Insert sorted.
00290                 if (t1 < tNear)
00291                 {
00292                     tFar = tNear;
00293                     tNear = t1;
00294                 }
00295                 else if (t1 < tFar)
00296                 {
00297                     tFar = t1;
00298                 }
00299             }
00300             else
00301             {
00302                 return;
00303             }
00304         }
00305         else
00306         {
00307             // Aligned with axis. Check if outside radius
00308             //Pout<< "small discriminant:" << disc
00309             //    << " for start:" << start << " end:" << end
00310             //    << " magV:" << magV
00311             //    << " c:" << c << endl;
00312             if (c > 0)
00313             {
00314                 return;
00315             }
00316         }
00317     }
00318     else
00319     {
00320         if (mag(a) > ROOTVSMALL)
00321         {
00322             scalar sqrtDisc = sqrt(disc);
00323 
00324             t1 = (-b - sqrtDisc)/(2*a);
00325             t2 = (-b + sqrtDisc)/(2*a);
00326             if (t2 < t1)
00327             {
00328                 Swap(t1, t2);
00329             }
00330 
00331             if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
00332             {
00333                 // valid. Insert sorted.
00334                 if (t1 < tNear)
00335                 {
00336                     tFar = tNear;
00337                     tNear = t1;
00338                 }
00339                 else if (t1 < tFar)
00340                 {
00341                     tFar = t1;
00342                 }
00343             }
00344             if (t2 >= 0 && t2 <= magV && t2 >= tPoint1 && t2 <= tPoint2)
00345             {
00346                 // valid. Insert sorted.
00347                 if (t2 < tNear)
00348                 {
00349                     tFar = tNear;
00350                     tNear = t2;
00351                 }
00352                 else if (t2 < tFar)
00353                 {
00354                     tFar = t2;
00355                 }
00356             }
00357             //Pout<< "two solutions t1:" << t1 << " t2:" << t2
00358             //    << " for start:" << start << " end:" << end
00359             //    << " magV:" << magV
00360             //    << " c:" << c << endl;
00361         }
00362         else
00363         {
00364             // Aligned with axis. Check if outside radius
00365             //Pout<< "large discriminant:" << disc
00366             //    << " small a:" << a
00367             //    << " for start:" << start << " end:" << end
00368             //    << " magV:" << magV
00369             //    << " c:" << c << endl;
00370             if (c > 0)
00371             {
00372                 return;
00373             }
00374         }
00375     }
00376 
00377     // Check tNear, tFar
00378     if (tNear >= 0 && tNear <= magV)
00379     {
00380         near.setPoint(start+tNear*V);
00381         near.setHit();
00382         near.setIndex(0);
00383 
00384         if (tFar <= magV)
00385         {
00386             far.setPoint(start+tFar*V);
00387             far.setHit();
00388             far.setIndex(0);
00389         }
00390     }
00391     else if (tFar >= 0 && tFar <= magV)
00392     {
00393         near.setPoint(start+tFar*V);
00394         near.setHit();
00395         near.setIndex(0);
00396     }
00397 }
00398 
00399 
00400 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00401 
00402 Foam::searchableCylinder::searchableCylinder
00403 (
00404     const IOobject& io,
00405     const point& point1,
00406     const point& point2,
00407     const scalar radius
00408 )
00409 :
00410     searchableSurface(io),
00411     point1_(point1),
00412     point2_(point2),
00413     magDir_(mag(point2_-point1_)),
00414     unitDir_((point2_-point1_)/magDir_),
00415     radius_(radius)
00416 {}
00417 
00418 
00419 Foam::searchableCylinder::searchableCylinder
00420 (
00421     const IOobject& io,
00422     const dictionary& dict
00423 )
00424 :
00425     searchableSurface(io),
00426     point1_(dict.lookup("point1")),
00427     point2_(dict.lookup("point2")),
00428     magDir_(mag(point2_-point1_)),
00429     unitDir_((point2_-point1_)/magDir_),
00430     radius_(readScalar(dict.lookup("radius")))
00431 {}
00432 
00433 
00434 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
00435 
00436 Foam::searchableCylinder::~searchableCylinder()
00437 {}
00438 
00439 
00440 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00441 
00442 const Foam::wordList& Foam::searchableCylinder::regions() const
00443 {
00444     if (regions_.empty())
00445     {
00446         regions_.setSize(1);
00447         regions_[0] = "region0";
00448     }
00449     return regions_;
00450 }
00451 
00452 
00453 void Foam::searchableCylinder::findNearest
00454 (
00455     const pointField& samples,
00456     const scalarField& nearestDistSqr,
00457     List<pointIndexHit>& info
00458 ) const
00459 {
00460     info.setSize(samples.size());
00461 
00462     forAll(samples, i)
00463     {
00464         info[i] = findNearest(samples[i], nearestDistSqr[i]);
00465     }
00466 }
00467 
00468 
00469 void Foam::searchableCylinder::findLine
00470 (
00471     const pointField& start,
00472     const pointField& end,
00473     List<pointIndexHit>& info
00474 ) const
00475 {
00476     info.setSize(start.size());
00477 
00478     forAll(start, i)
00479     {
00480         // Pick nearest intersection. If none intersected take second one.
00481         pointIndexHit b;
00482         findLineAll(start[i], end[i], info[i], b);
00483         if (!info[i].hit() && b.hit())
00484         {
00485             info[i] = b;
00486         }
00487     }
00488 }
00489 
00490 
00491 void Foam::searchableCylinder::findLineAny
00492 (
00493     const pointField& start,
00494     const pointField& end,
00495     List<pointIndexHit>& info
00496 ) const
00497 {
00498     info.setSize(start.size());
00499 
00500     forAll(start, i)
00501     {
00502         // Discard far intersection
00503         pointIndexHit b;
00504         findLineAll(start[i], end[i], info[i], b);
00505         if (!info[i].hit() && b.hit())
00506         {
00507             info[i] = b;
00508         }
00509     }
00510 }
00511 
00512 
00513 void Foam::searchableCylinder::findLineAll
00514 (
00515     const pointField& start,
00516     const pointField& end,
00517     List<List<pointIndexHit> >& info
00518 ) const
00519 {
00520     info.setSize(start.size());
00521 
00522     forAll(start, i)
00523     {
00524         pointIndexHit near, far;
00525         findLineAll(start[i], end[i], near, far);
00526 
00527         if (near.hit())
00528         {
00529             if (far.hit())
00530             {
00531                 info[i].setSize(2);
00532                 info[i][0] = near;
00533                 info[i][1] = far;
00534             }
00535             else
00536             {
00537                 info[i].setSize(1);
00538                 info[i][0] = near;
00539             }
00540         }
00541         else
00542         {
00543             if (far.hit())
00544             {
00545                 info[i].setSize(1);
00546                 info[i][0] = far;
00547             }
00548             else
00549             {
00550                 info[i].clear();
00551             }
00552         }
00553     }
00554 }
00555 
00556 
00557 void Foam::searchableCylinder::getRegion
00558 (
00559     const List<pointIndexHit>& info,
00560     labelList& region
00561 ) const
00562 {
00563     region.setSize(info.size());
00564     region = 0;
00565 }
00566 
00567 
00568 void Foam::searchableCylinder::getNormal
00569 (
00570     const List<pointIndexHit>& info,
00571     vectorField& normal
00572 ) const
00573 {
00574     normal.setSize(info.size());
00575     normal = vector::zero;
00576 
00577     forAll(info, i)
00578     {
00579         if (info[i].hit())
00580         {
00581             vector v(info[i].hitPoint() - point1_);
00582 
00583             // Decompose sample-point1 into normal and parallel component
00584             scalar parallel = v & unitDir_;
00585 
00586             if (parallel < 0)
00587             {
00588                 normal[i] = -unitDir_;
00589             }
00590             else if (parallel > magDir_)
00591             {
00592                 normal[i] = -unitDir_;
00593             }
00594             else
00595             {
00596                 // Remove the parallel component
00597                 v -= parallel*unitDir_;
00598                 normal[i] = v/mag(v);
00599             }
00600         }
00601     }
00602 }
00603 
00604 
00605 void Foam::searchableCylinder::getVolumeType
00606 (
00607     const pointField& points,
00608     List<volumeType>& volType
00609 ) const
00610 {
00611     volType.setSize(points.size());
00612     volType = INSIDE;
00613 
00614     forAll(points, pointI)
00615     {
00616         const point& pt = points[pointI];
00617 
00618         vector v(pt - point1_);
00619 
00620         // Decompose sample-point1 into normal and parallel component
00621         scalar parallel = v & unitDir_;
00622 
00623         if (parallel < 0)
00624         {
00625             // left of point1 endcap
00626             volType[pointI] = OUTSIDE;
00627         }
00628         else if (parallel > magDir_)
00629         {
00630             // right of point2 endcap
00631             volType[pointI] = OUTSIDE;
00632         }
00633         else
00634         {
00635             // Remove the parallel component
00636             v -= parallel*unitDir_;
00637 
00638             if (mag(v) > radius_)
00639             {
00640                 volType[pointI] = OUTSIDE;
00641             }
00642             else
00643             {
00644                 volType[pointI] = INSIDE;
00645             }
00646         }
00647     }
00648 }
00649 
00650 
00651 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines