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

searchableSurfaceWithGaps.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 "searchableSurfaceWithGaps.H"
00027 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00028 #include <OpenFOAM/Time.H>
00029 #include <OpenFOAM/ListOps.H>
00030 
00031 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00032 
00033 namespace Foam
00034 {
00035 
00036 defineTypeNameAndDebug(searchableSurfaceWithGaps, 0);
00037 addToRunTimeSelectionTable(searchableSurface, searchableSurfaceWithGaps, dict);
00038 
00039 }
00040 
00041 
00042 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00043 
00044 Foam::Pair<Foam::vector> Foam::searchableSurfaceWithGaps::offsetVecs
00045 (
00046     const point& start,
00047     const point& end
00048 ) const
00049 {
00050     Pair<vector> offsets(vector::zero, vector::zero);
00051 
00052     vector n(end-start);
00053 
00054     scalar magN = mag(n);
00055 
00056     if (magN > SMALL)
00057     {
00058         n /= magN;
00059 
00060         // Do first offset vector. Is the coordinate axes with the smallest
00061         // component along the vector n.
00062         scalar minMag = GREAT;
00063         direction minCmpt = 0;
00064 
00065         for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
00066         {
00067             if (mag(n[cmpt]) < minMag)
00068             {
00069                 minMag = mag(n[cmpt]);
00070                 minCmpt = cmpt;
00071             }
00072         }
00073 
00074         offsets[0][minCmpt] = 1.0;
00075         // Orthonormalise
00076         offsets[0] -= n[minCmpt]*n;
00077         offsets[0] /= mag(offsets[0]);
00078         // Do second offset vector perp to original edge and first offset vector
00079         offsets[1] = n ^ offsets[0];
00080 
00081         // Scale
00082         offsets[0] *= gap_;
00083         offsets[1] *= gap_;
00084     }
00085 
00086     return offsets;
00087 }
00088 
00089 
00090 void Foam::searchableSurfaceWithGaps::offsetVecs
00091 (
00092     const pointField& start,
00093     const pointField& end,
00094     pointField& offset0,
00095     pointField& offset1
00096 ) const
00097 {
00098     offset0.setSize(start.size());
00099     offset1.setSize(start.size());
00100 
00101     forAll(start, i)
00102     {
00103         const Pair<vector> offsets(offsetVecs(start[i], end[i]));
00104         offset0[i] = offsets[0];
00105         offset1[i] = offsets[1];
00106     }
00107 }
00108 
00109 
00110 Foam::label Foam::searchableSurfaceWithGaps::countMisses
00111 (
00112     const List<pointIndexHit>& info,
00113     labelList& missMap
00114 )
00115 {
00116     label nMiss = 0;
00117     forAll(info, i)
00118     {
00119         if (!info[i].hit())
00120         {
00121             nMiss++;
00122         }
00123     }
00124 
00125     missMap.setSize(nMiss);
00126     nMiss = 0;
00127 
00128     forAll(info, i)
00129     {
00130         if (!info[i].hit())
00131         {
00132             missMap[nMiss++] = i;
00133         }
00134     }
00135 
00136     return nMiss;
00137 }
00138 
00139 
00140 // Anything not a hit in both counts as a hit
00141 Foam::label Foam::searchableSurfaceWithGaps::countMisses
00142 (
00143     const List<pointIndexHit>& plusInfo,
00144     const List<pointIndexHit>& minInfo,
00145     labelList& missMap
00146 )
00147 {
00148     label nMiss = 0;
00149     forAll(plusInfo, i)
00150     {
00151         if (!plusInfo[i].hit() || !minInfo[i].hit())
00152         {
00153             nMiss++;
00154         }
00155     }
00156 
00157     missMap.setSize(nMiss);
00158     nMiss = 0;
00159 
00160     forAll(plusInfo, i)
00161     {
00162         if (!plusInfo[i].hit() || !minInfo[i].hit())
00163         {
00164             missMap[nMiss++] = i;
00165         }
00166     }
00167 
00168     return nMiss;
00169 }
00170 
00171 
00172 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00173 
00174 Foam::searchableSurfaceWithGaps::searchableSurfaceWithGaps
00175 (
00176     const IOobject& io,
00177     const dictionary& dict
00178 )
00179 :
00180     searchableSurface(io),
00181     gap_(readScalar(dict.lookup("gap"))),
00182     subGeom_(1)
00183 {
00184     const word subGeomName(dict.lookup("surface"));
00185 
00186     const searchableSurface& s =
00187         io.db().lookupObject<searchableSurface>(subGeomName);
00188 
00189     subGeom_.set(0, &const_cast<searchableSurface&>(s));
00190 }
00191 
00192 
00193 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
00194 
00195 Foam::searchableSurfaceWithGaps::~searchableSurfaceWithGaps()
00196 {}
00197 
00198 
00199 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00200 
00201 void Foam::searchableSurfaceWithGaps::findLine
00202 (
00203     const pointField& start,
00204     const pointField& end,
00205     List<pointIndexHit>& info
00206 ) const
00207 {
00208 
00209     // Test with unperturbed vectors
00210     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00211 
00212     surface().findLine(start, end, info);
00213 
00214     // Count number of misses. Determine map
00215     labelList compactMap;
00216     label nMiss = countMisses(info, compactMap);
00217 
00218     if (returnReduce(nMiss, sumOp<label>()) > 0)
00219     {
00220         //Pout<< "** retesting with offset0 " << nMiss << " misses out of "
00221         //    << start.size() << endl;
00222 
00223         // extract segments according to map
00224         pointField compactStart(start, compactMap);
00225         pointField compactEnd(end, compactMap);
00226 
00227         // Calculate offset vector
00228         pointField offset0, offset1;
00229         offsetVecs
00230         (
00231             compactStart,
00232             compactEnd,
00233             offset0,
00234             offset1
00235         );
00236 
00237         // Test with offset0 perturbed vectors
00238         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00239 
00240         // test in pairs: only if both perturbations hit something
00241         // do we accept the hit.
00242 
00243         const vectorField smallVec(SMALL*(compactEnd-compactStart));
00244 
00245         List<pointIndexHit> plusInfo;
00246         surface().findLine
00247         (
00248             compactStart+offset0-smallVec,
00249             compactEnd+offset0+smallVec,
00250             plusInfo
00251         );
00252         List<pointIndexHit> minInfo;
00253         surface().findLine
00254         (
00255             compactStart-offset0-smallVec,
00256             compactEnd-offset0+smallVec,
00257             minInfo
00258         );
00259 
00260         // Extract any hits
00261         forAll(plusInfo, i)
00262         {
00263             if (plusInfo[i].hit() && minInfo[i].hit())
00264             {
00265                 info[compactMap[i]] = plusInfo[i];
00266                 info[compactMap[i]].rawPoint() -= offset0[i];
00267             }
00268         }
00269 
00270         labelList plusMissMap;
00271         nMiss = countMisses(plusInfo, minInfo, plusMissMap);
00272 
00273         if (returnReduce(nMiss, sumOp<label>()) > 0)
00274         {
00275             //Pout<< "** retesting with offset1 " << nMiss << " misses out of "
00276             //    << start.size() << endl;
00277 
00278             // Test with offset1 perturbed vectors
00279             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00280 
00281             // Extract (inplace possible because of order)
00282             forAll(plusMissMap, i)
00283             {
00284                 label mapI = plusMissMap[i];
00285                 compactStart[i] = compactStart[mapI];
00286                 compactEnd[i] = compactEnd[mapI];
00287                 compactMap[i] = compactMap[mapI];
00288                 offset0[i] = offset0[mapI];
00289                 offset1[i] = offset1[mapI];
00290             }
00291             compactStart.setSize(plusMissMap.size());
00292             compactEnd.setSize(plusMissMap.size());
00293             compactMap.setSize(plusMissMap.size());
00294             offset0.setSize(plusMissMap.size());
00295             offset1.setSize(plusMissMap.size());
00296 
00297             const vectorField smallVec(SMALL*(compactEnd-compactStart));
00298 
00299             surface().findLine
00300             (
00301                 compactStart+offset1-smallVec,
00302                 compactEnd+offset1+smallVec,
00303                 plusInfo
00304             );
00305             surface().findLine
00306             (
00307                 compactStart-offset1-smallVec,
00308                 compactEnd-offset1+smallVec,
00309                 minInfo
00310             );
00311 
00312             // Extract any hits
00313             forAll(plusInfo, i)
00314             {
00315                 if (plusInfo[i].hit() && minInfo[i].hit())
00316                 {
00317                     info[compactMap[i]] = plusInfo[i];
00318                     info[compactMap[i]].rawPoint() -= offset1[i];
00319                 }
00320             }
00321         }
00322     }
00323 }
00324 
00325 
00326 void Foam::searchableSurfaceWithGaps::findLineAny
00327 (
00328     const pointField& start,
00329     const pointField& end,
00330     List<pointIndexHit>& info
00331 ) const
00332 {
00333     // To be done ...
00334     findLine(start, end, info);
00335 }
00336 
00337 
00338 void Foam::searchableSurfaceWithGaps::findLineAll
00339 (
00340     const pointField& start,
00341     const pointField& end,
00342     List<List<pointIndexHit> >& info
00343 ) const
00344 {
00345     // To be done. Assume for now only one intersection.
00346     List<pointIndexHit> nearestInfo;
00347     findLine(start, end, nearestInfo);
00348 
00349     info.setSize(start.size());
00350     forAll(info, pointI)
00351     {
00352         if (nearestInfo[pointI].hit())
00353         {
00354             info[pointI].setSize(1);
00355             info[pointI][0] = nearestInfo[pointI];
00356         }
00357         else
00358         {
00359             info[pointI].clear();
00360         }
00361     }
00362 }
00363 
00364 
00365 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines