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