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

searchableBox.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-2010 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 "searchableBox.H"
00027 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00028 #include <OpenFOAM/SortableList.H>
00029 
00030 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00031 
00032 namespace Foam
00033 {
00034     defineTypeNameAndDebug(searchableBox, 0);
00035     addToRunTimeSelectionTable(searchableSurface, searchableBox, dict);
00036 }
00037 
00038 
00039 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00040 
00041 void Foam::searchableBox::projectOntoCoordPlane
00042 (
00043     const direction dir,
00044     const point& planePt,
00045     pointIndexHit& info
00046 ) const
00047 {
00048     // Set point
00049     info.rawPoint()[dir] = planePt[dir];
00050     // Set face
00051     if (planePt[dir] == min()[dir])
00052     {
00053         info.setIndex(dir*2);
00054     }
00055     else if (planePt[dir] == max()[dir])
00056     {
00057         info.setIndex(dir*2+1);
00058     }
00059     else
00060     {
00061         FatalErrorIn("searchableBox::projectOntoCoordPlane(..)")
00062             << "Point on plane " << planePt
00063             << " is not on coordinate " << min()[dir]
00064             << " nor " << max()[dir] << abort(FatalError);
00065     }
00066 }
00067 
00068 
00069 // Returns miss or hit with face (0..5) and region(always 0)
00070 Foam::pointIndexHit Foam::searchableBox::findNearest
00071 (
00072     const point& bbMid,
00073     const point& sample,
00074     const scalar nearestDistSqr
00075 ) const
00076 {
00077     // Point can be inside or outside. For every component direction can be
00078     // left of min, right of max or inbetween.
00079     // - outside points: project first one x plane (either min().x()
00080     // or max().x()), then onto y plane and finally z. You should be left
00081     // with intersection point
00082     // - inside point: find nearest side (compare to mid point). Project onto
00083     //   that.
00084 
00085     // The face is set to the last projected face.
00086 
00087 
00088     // Outside point projected onto cube. Assume faces 0..5.
00089     pointIndexHit info(true, sample, -1);
00090     bool outside = false;
00091 
00092     // (for internal points) per direction what nearest cube side is
00093     point near;
00094 
00095     for (direction dir = 0; dir < vector::nComponents; dir++)
00096     {
00097         if (info.rawPoint()[dir] < min()[dir])
00098         {
00099             projectOntoCoordPlane(dir, min(), info);
00100             outside = true;
00101         }
00102         else if (info.rawPoint()[dir] > max()[dir])
00103         {
00104             projectOntoCoordPlane(dir, max(), info);
00105             outside = true;
00106         }
00107         else if (info.rawPoint()[dir] > bbMid[dir])
00108         {
00109             near[dir] = max()[dir];
00110         }
00111         else
00112         {
00113             near[dir] = min()[dir];
00114         }
00115     }
00116 
00117 
00118     // For outside points the info will be correct now. Handle inside points
00119     // using the three near distances. Project onto the nearest plane.
00120     if (!outside)
00121     {
00122         vector dist(cmptMag(info.rawPoint() - near));
00123 
00124         if (dist.x() < dist.y())
00125         {
00126             if (dist.x() < dist.z())
00127             {
00128                 // Project onto x plane
00129                 projectOntoCoordPlane(vector::X, near, info);
00130             }
00131             else
00132             {
00133                 projectOntoCoordPlane(vector::Z, near, info);
00134             }
00135         }
00136         else
00137         {
00138             if (dist.y() < dist.z())
00139             {
00140                 projectOntoCoordPlane(vector::Y, near, info);
00141             }
00142             else
00143             {
00144                 projectOntoCoordPlane(vector::Z, near, info);
00145             }
00146         }
00147     }
00148 
00149 
00150     // Check if outside. Optimisation: could do some checks on distance already
00151     // on components above
00152     if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
00153     {
00154         info.setMiss();
00155         info.setIndex(-1);
00156     }
00157 
00158     return info;
00159 }
00160 
00161 
00162 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00163 
00164 Foam::searchableBox::searchableBox
00165 (
00166     const IOobject& io,
00167     const treeBoundBox& bb
00168 )
00169 :
00170     searchableSurface(io),
00171     treeBoundBox(bb)
00172 {
00173     if (!contains(midpoint()))
00174     {
00175         FatalErrorIn
00176         (
00177             "Foam::searchableBox::searchableBox\n"
00178             "(\n"
00179             "    const IOobject& io,\n"
00180             "    const treeBoundBox& bb\n"
00181             ")\n"
00182         )   << "Illegal bounding box specification : "
00183             << static_cast<const treeBoundBox>(*this) << exit(FatalError);
00184     }
00185 }
00186 
00187 
00188 Foam::searchableBox::searchableBox
00189 (
00190     const IOobject& io,
00191     const dictionary& dict
00192 )
00193 :
00194     searchableSurface(io),
00195     treeBoundBox(dict.lookup("min"), dict.lookup("max"))
00196 {
00197     if (!contains(midpoint()))
00198     {
00199         FatalErrorIn
00200         (
00201             "Foam::searchableBox::searchableBox\n"
00202             "(\n"
00203             "    const IOobject& io,\n"
00204             "    const treeBoundBox& bb\n"
00205             ")\n"
00206         )   << "Illegal bounding box specification : "
00207             << static_cast<const treeBoundBox>(*this) << exit(FatalError);
00208     }
00209 }
00210 
00211 
00212 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
00213 
00214 Foam::searchableBox::~searchableBox()
00215 {}
00216 
00217 
00218 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00219 
00220 const Foam::wordList& Foam::searchableBox::regions() const
00221 {
00222     if (regions_.empty())
00223     {
00224         regions_.setSize(1);
00225         regions_[0] = "region0";
00226     }
00227     return regions_;
00228 }
00229 
00230 
00231 Foam::pointField Foam::searchableBox::coordinates() const
00232 {
00233     pointField ctrs(6);
00234 
00235     const pointField pts = treeBoundBox::points();
00236     const faceList& fcs = treeBoundBox::faces;
00237 
00238     forAll(fcs, i)
00239     {
00240         ctrs[i] = fcs[i].centre(pts);
00241     }
00242     return ctrs;
00243 }
00244 
00245 
00246 Foam::pointIndexHit Foam::searchableBox::findNearest
00247 (
00248     const point& sample,
00249     const scalar nearestDistSqr
00250 ) const
00251 {
00252     return findNearest(midpoint(), sample, nearestDistSqr);
00253 }
00254 
00255 
00256 Foam::pointIndexHit Foam::searchableBox::findNearestOnEdge
00257 (
00258     const point& sample,
00259     const scalar nearestDistSqr
00260 ) const
00261 {
00262     const point bbMid(midpoint());
00263 
00264     // Outside point projected onto cube. Assume faces 0..5.
00265     pointIndexHit info(true, sample, -1);
00266     bool outside = false;
00267 
00268     // (for internal points) per direction what nearest cube side is
00269     point near;
00270 
00271     for (direction dir = 0; dir < vector::nComponents; dir++)
00272     {
00273         if (info.rawPoint()[dir] < min()[dir])
00274         {
00275             projectOntoCoordPlane(dir, min(), info);
00276             outside = true;
00277         }
00278         else if (info.rawPoint()[dir] > max()[dir])
00279         {
00280             projectOntoCoordPlane(dir, max(), info);
00281             outside = true;
00282         }
00283         else if (info.rawPoint()[dir] > bbMid[dir])
00284         {
00285             near[dir] = max()[dir];
00286         }
00287         else
00288         {
00289             near[dir] = min()[dir];
00290         }
00291     }
00292 
00293 
00294     // For outside points the info will be correct now. Handle inside points
00295     // using the three near distances. Project onto the nearest two planes.
00296     if (!outside)
00297     {
00298         // Get the per-component distance to nearest wall
00299         vector dist(cmptMag(info.rawPoint() - near));
00300 
00301         SortableList<scalar> sortedDist(3);
00302         sortedDist[0] = dist[0];
00303         sortedDist[1] = dist[1];
00304         sortedDist[2] = dist[2];
00305         sortedDist.sort();
00306 
00307         // Project onto nearest
00308         projectOntoCoordPlane(sortedDist.indices()[0], near, info);
00309         // Project onto second nearest
00310         projectOntoCoordPlane(sortedDist.indices()[1], near, info);
00311     }
00312 
00313 
00314     // Check if outside. Optimisation: could do some checks on distance already
00315     // on components above
00316     if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
00317     {
00318         info.setMiss();
00319         info.setIndex(-1);
00320     }
00321 
00322     return info;
00323 }
00324 
00325 
00326 Foam::pointIndexHit Foam::searchableBox::findNearest
00327 (
00328     const linePointRef& ln,
00329     treeBoundBox& tightest,
00330     point& linePoint
00331 ) const
00332 {
00333     notImplemented
00334     (
00335         "searchableBox::findNearest"
00336         "(const linePointRef&, treeBoundBox&, point&)"
00337     );
00338     return pointIndexHit();
00339 }
00340 
00341 
00342 Foam::pointIndexHit Foam::searchableBox::findLine
00343 (
00344     const point& start,
00345     const point& end
00346 ) const
00347 {
00348     pointIndexHit info(false, start, -1);
00349 
00350     bool foundInter;
00351 
00352     if (posBits(start) == 0)
00353     {
00354         if (posBits(end) == 0)
00355         {
00356             // Both start and end inside.
00357             foundInter = false;
00358         }
00359         else
00360         {
00361             // end is outside. Clip to bounding box.
00362             foundInter = intersects(end, start, info.rawPoint());
00363         }
00364     }
00365     else
00366     {
00367         // start is outside. Clip to bounding box.
00368         foundInter = intersects(start, end, info.rawPoint());
00369     }
00370 
00371 
00372     // Classify point
00373     if (foundInter)
00374     {
00375         info.setHit();
00376 
00377         for (direction dir = 0; dir < vector::nComponents; dir++)
00378         {
00379             if (info.rawPoint()[dir] == min()[dir])
00380             {
00381                 info.setIndex(2*dir);
00382                 break;
00383             }
00384             else if (info.rawPoint()[dir] == max()[dir])
00385             {
00386                 info.setIndex(2*dir+1);
00387                 break;
00388             }
00389         }
00390 
00391         if (info.index() == -1)
00392         {
00393             FatalErrorIn("searchableBox::findLine(const point&, const point&)")
00394                 << "point " << info.rawPoint()
00395                 << " on segment " << start << end
00396                 << " should be on face of " << *this
00397                 << " but it isn't." << abort(FatalError);
00398         }
00399     }
00400 
00401     return info;
00402 }
00403 
00404 
00405 Foam::pointIndexHit Foam::searchableBox::findLineAny
00406 (
00407     const point& start,
00408     const point& end
00409 ) const
00410 {
00411     return findLine(start, end);
00412 }
00413 
00414 
00415 void Foam::searchableBox::findNearest
00416 (
00417     const pointField& samples,
00418     const scalarField& nearestDistSqr,
00419     List<pointIndexHit>& info
00420 ) const
00421 {
00422     info.setSize(samples.size());
00423 
00424     const point bbMid(midpoint());
00425 
00426     forAll(samples, i)
00427     {
00428         info[i] = findNearest(bbMid, samples[i], nearestDistSqr[i]);
00429     }
00430 }
00431 
00432 
00433 void Foam::searchableBox::findLine
00434 (
00435     const pointField& start,
00436     const pointField& end,
00437     List<pointIndexHit>& info
00438 ) const
00439 {
00440     info.setSize(start.size());
00441 
00442     forAll(start, i)
00443     {
00444         info[i] = findLine(start[i], end[i]);
00445     }
00446 }
00447 
00448 
00449 void Foam::searchableBox::findLineAny
00450 (
00451     const pointField& start,
00452     const pointField& end,
00453     List<pointIndexHit>& info
00454 ) const
00455 {
00456     info.setSize(start.size());
00457 
00458     forAll(start, i)
00459     {
00460         info[i] = findLineAny(start[i], end[i]);
00461     }
00462 }
00463 
00464 
00465 void Foam::searchableBox::findLineAll
00466 (
00467     const pointField& start,
00468     const pointField& end,
00469     List<List<pointIndexHit> >& info
00470 ) const
00471 {
00472     info.setSize(start.size());
00473 
00474     // Work array
00475     DynamicList<pointIndexHit, 1, 1> hits;
00476 
00477     // Tolerances:
00478     // To find all intersections we add a small vector to the last intersection
00479     // This is chosen such that
00480     // - it is significant (SMALL is smallest representative relative tolerance;
00481     //   we need something bigger since we're doing calculations)
00482     // - if the start-end vector is zero we still progress
00483     const vectorField dirVec(end-start);
00484     const scalarField magSqrDirVec(magSqr(dirVec));
00485     const vectorField smallVec
00486     (
00487         Foam::sqrt(SMALL)*dirVec
00488       + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
00489     );
00490 
00491     forAll(start, pointI)
00492     {
00493         // See if any intersection between pt and end
00494         pointIndexHit inter = findLine(start[pointI], end[pointI]);
00495 
00496         if (inter.hit())
00497         {
00498             hits.clear();
00499             hits.append(inter);
00500 
00501             point pt = inter.hitPoint() + smallVec[pointI];
00502 
00503             while (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
00504             {
00505                 // See if any intersection between pt and end
00506                 pointIndexHit inter = findLine(pt, end[pointI]);
00507 
00508                 // Check for not hit or hit same face as before (can happen
00509                 // if vector along surface of face)
00510                 if
00511                 (
00512                     !inter.hit()
00513                  || (inter.index() == hits[hits.size()-1].index())
00514                 )
00515                 {
00516                     break;
00517                 }
00518                 hits.append(inter);
00519 
00520                 pt = inter.hitPoint() + smallVec[pointI];
00521             }
00522 
00523             info[pointI].transfer(hits);
00524         }
00525         else
00526         {
00527             info[pointI].clear();
00528         }
00529     }
00530 }
00531 
00532 
00533 void Foam::searchableBox::getRegion
00534 (
00535     const List<pointIndexHit>& info,
00536     labelList& region
00537 ) const
00538 {
00539     region.setSize(info.size());
00540     region = 0;
00541 }
00542 
00543 
00544 void Foam::searchableBox::getNormal
00545 (
00546     const List<pointIndexHit>& info,
00547     vectorField& normal
00548 ) const
00549 {
00550     normal.setSize(info.size());
00551     normal = vector::zero;
00552 
00553     forAll(info, i)
00554     {
00555         if (info[i].hit())
00556         {
00557             normal[i] = treeBoundBox::faceNormals[info[i].index()];
00558         }
00559         else
00560         {
00561             // Set to what?
00562         }
00563     }
00564 }
00565 
00566 
00567 void Foam::searchableBox::getVolumeType
00568 (
00569     const pointField& points,
00570     List<volumeType>& volType
00571 ) const
00572 {
00573     volType.setSize(points.size());
00574     volType = INSIDE;
00575 
00576     forAll(points, pointI)
00577     {
00578         const point& pt = points[pointI];
00579 
00580         for (direction dir = 0; dir < vector::nComponents; dir++)
00581         {
00582             if (pt[dir] < min()[dir] || pt[dir] > max()[dir])
00583             {
00584                 volType[pointI] = OUTSIDE;
00585                 break;
00586             }
00587         }
00588     }
00589 }
00590 
00591 
00592 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines