src/meshTools/searchableSurface/searchableSurfaceWithGaps.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 "searchableSurfaceWithGaps.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "Time.H"
  29 #include "ListOps.H"
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35
  36 defineTypeNameAndDebug(searchableSurfaceWithGaps, 0);
  37 addToRunTimeSelectionTable(searchableSurface, searchableSurfaceWithGaps, dict);
  38
  39 }
  40
  41
  42 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  43
  44 Foam::Pair<Foam::vector> Foam::searchableSurfaceWithGaps::offsetVecs
  45 (
  46     const point& start,
  47     const point& end
  48 ) const
  49 {
  50     Pair<vector> offsets(vector::zero, vector::zero);
  51
  52     vector n(end-start);
  53
  54     scalar magN = mag(n);
  55
  56     if (magN > SMALL)
  57     {
  58         n /= magN;
  59
  60         // Do first offset vector. Is the coordinate axes with the smallest
  61         // component along the vector n.
  62         scalar minMag = GREAT;
  63         direction minCmpt = 0;
  64
  65         for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
  66         {
  67             if (mag(n[cmpt]) < minMag)
  68             {
  69                 minMag = mag(n[cmpt]);
  70                 minCmpt = cmpt;
  71             }
  72         }
  73
  74         offsets[0][minCmpt] = 1.0;
  75         // Orthonormalise
  76         offsets[0] -= n[minCmpt]*n;
  77         offsets[0] /= mag(offsets[0]);
  78         // Do second offset vector perp to original edge and first offset vector
  79         offsets[1] = n ^ offsets[0];
  80
  81         // Scale
  82         offsets[0] *= gap_;
  83         offsets[1] *= gap_;
  84     }
  85
  86     return offsets;
  87 }
  88
  89
  90 void Foam::searchableSurfaceWithGaps::offsetVecs
  91 (
  92     const pointField& start,
  93     const pointField& end,
  94     pointField& offset0,
  95     pointField& offset1
  96 ) const
  97 {
  98     offset0.setSize(start.size());
  99     offset1.setSize(start.size());
 100
 101     forAll(start, i)
 102     {
 103         const Pair<vector> offsets(offsetVecs(start[i], end[i]));
 104         offset0[i] = offsets[0];
 105         offset1[i] = offsets[1];
 106     }
 107 }
 108
 109
 110 Foam::label Foam::searchableSurfaceWithGaps::countMisses
 111 (
 112     const List<pointIndexHit>& info,
 113     labelList& missMap
 114 )
 115 {
 116     label nMiss = 0;
 117     forAll(info, i)
 118     {
 119         if (!info[i].hit())
 120         {
 121             nMiss++;
 122         }
 123     }
 124
 125     missMap.setSize(nMiss);
 126     nMiss = 0;
 127
 128     forAll(info, i)
 129     {
 130         if (!info[i].hit())
 131         {
 132             missMap[nMiss++] = i;
 133         }
 134     }
 135
 136     return nMiss;
 137 }
 138
 139
 140 // Anything not a hit in both counts as a hit
 141 Foam::label Foam::searchableSurfaceWithGaps::countMisses
 142 (
 143     const List<pointIndexHit>& plusInfo,
 144     const List<pointIndexHit>& minInfo,
 145     labelList& missMap
 146 )
 147 {
 148     label nMiss = 0;
 149     forAll(plusInfo, i)
 150     {
 151         if (!plusInfo[i].hit() || !minInfo[i].hit())
 152         {
 153             nMiss++;
 154         }
 155     }
 156
 157     missMap.setSize(nMiss);
 158     nMiss = 0;
 159
 160     forAll(plusInfo, i)
 161     {
 162         if (!plusInfo[i].hit() || !minInfo[i].hit())
 163         {
 164             missMap[nMiss++] = i;
 165         }
 166     }
 167
 168     return nMiss;
 169 }
 170
 171
 172 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 173
 174 Foam::searchableSurfaceWithGaps::searchableSurfaceWithGaps
 175 (
 176     const IOobject& io,
 177     const dictionary& dict
 178 )
 179 :
 180     searchableSurface(io),
 181     gap_(readScalar(dict.lookup("gap"))),
 182     subGeom_(1)
 183 {
 184     const word subGeomName(dict.lookup("surface"));
 185
 186     const searchableSurface& s =
 187         io.db().lookupObject<searchableSurface>(subGeomName);
 188
 189     subGeom_.set(0, &const_cast<searchableSurface&>(s));
 190 }
 191
 192
 193 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 194
 195 Foam::searchableSurfaceWithGaps::~searchableSurfaceWithGaps()
 196 {}
 197
 198
 199 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 200
 201 void Foam::searchableSurfaceWithGaps::findLine
 202 (
 203     const pointField& start,
 204     const pointField& end,
 205     List<pointIndexHit>& info
 206 ) const
 207 {
 208
 209     // Test with unperturbed vectors
 210     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 211
 212     surface().findLine(start, end, info);
 213
 214     // Count number of misses. Determine map
 215     labelList compactMap;
 216     label nMiss = countMisses(info, compactMap);
 217
 218     if (returnReduce(nMiss, sumOp<label>()) > 0)
 219     {
 220         //Pout<< "** retesting with offset0 " << nMiss << " misses out of "
 221         //    << start.size() << endl;
 222
 223         // extract segments according to map
 224         pointField compactStart(start, compactMap);
 225         pointField compactEnd(end, compactMap);
 226
 227         // Calculate offset vector
 228         pointField offset0, offset1;
 229         offsetVecs
 230         (
 231             compactStart,
 232             compactEnd,
 233             offset0,
 234             offset1
 235         );
 236
 237         // Test with offset0 perturbed vectors
 238         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 239
 240         // test in pairs: only if both perturbations hit something
 241         // do we accept the hit.
 242
 243         const vectorField smallVec(1E-6*(compactEnd-compactStart));
 244
 245         List<pointIndexHit> plusInfo;
 246         surface().findLine
 247         (
 248             compactStart+offset0-smallVec,
 249             compactEnd+offset0+smallVec,
 250             plusInfo
 251         );
 252         List<pointIndexHit> minInfo;
 253         surface().findLine
 254         (
 255             compactStart-offset0-smallVec,
 256             compactEnd-offset0+smallVec,
 257             minInfo
 258         );
 259
 260         // Extract any hits
 261         forAll(plusInfo, i)
 262         {
 263             if (plusInfo[i].hit() && minInfo[i].hit())
 264             {
 265                 info[compactMap[i]] = plusInfo[i];
 266                 info[compactMap[i]].rawPoint() -= offset0[i];
 267             }
 268         }
 269
 270         labelList plusMissMap;
 271         nMiss = countMisses(plusInfo, minInfo, plusMissMap);
 272
 273         if (returnReduce(nMiss, sumOp<label>()) > 0)
 274         {
 275             //Pout<< "** retesting with offset1 " << nMiss << " misses out of "
 276             //    << start.size() << endl;
 277
 278             // Test with offset1 perturbed vectors
 279             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 280
 281             // Extract (inplace possible because of order)
 282             forAll(plusMissMap, i)
 283             {
 284                 label mapI = plusMissMap[i];
 285                 compactStart[i] = compactStart[mapI];
 286                 compactEnd[i] = compactEnd[mapI];
 287                 compactMap[i] = compactMap[mapI];
 288                 offset0[i] = offset0[mapI];
 289                 offset1[i] = offset1[mapI];
 290             }
 291             compactStart.setSize(plusMissMap.size());
 292             compactEnd.setSize(plusMissMap.size());
 293             compactMap.setSize(plusMissMap.size());
 294             offset0.setSize(plusMissMap.size());
 295             offset1.setSize(plusMissMap.size());
 296
 297             const vectorField smallVec(1E-6*(compactEnd-compactStart));
 298
 299             surface().findLine
 300             (
 301                 compactStart+offset1-smallVec,
 302                 compactEnd+offset1+smallVec,
 303                 plusInfo
 304             );
 305             surface().findLine
 306             (
 307                 compactStart-offset1-smallVec,
 308                 compactEnd-offset1+smallVec,
 309                 minInfo
 310             );
 311
 312             // Extract any hits
 313             forAll(plusInfo, i)
 314             {
 315                 if (plusInfo[i].hit() && minInfo[i].hit())
 316                 {
 317                     info[compactMap[i]] = plusInfo[i];
 318                     info[compactMap[i]].rawPoint() -= offset1[i];
 319                 }
 320             }
 321         }
 322     }
 323 }
 324
 325
 326 void Foam::searchableSurfaceWithGaps::findLineAny
 327 (
 328     const pointField& start,
 329     const pointField& end,
 330     List<pointIndexHit>& info
 331 ) const
 332 {
 333     // To be done ...
 334     findLine(start, end, info);
 335 }
 336
 337
 338 void Foam::searchableSurfaceWithGaps::findLineAll
 339 (
 340     const pointField& start,
 341     const pointField& end,
 342     List<List<pointIndexHit> >& info
 343 ) const
 344 {
 345     // To be done. Assume for now only one intersection.
 346     List<pointIndexHit> nearestInfo;
 347     findLine(start, end, nearestInfo);
 348
 349     info.setSize(start.size());
 350     forAll(info, pointI)
 351     {
 352         if (nearestInfo[pointI].hit())
 353         {
 354             info[pointI].setSize(1);
 355             info[pointI][0] = nearestInfo[pointI];
 356         }
 357         else
 358         {
 359             info[pointI].clear();
 360         }
 361     }
 362 }
 363
 364
 365 // ************************************************************************* //