src/sampling/sampledSet/face/faceOnlySet.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 "faceOnlySet.H"
  28 #include "meshSearch.H"
  29 #include "DynamicList.H"
  30 #include "polyMesh.H"
  31
  32 #include "Cloud.H"
  33 #include "passiveParticle.H"
  34 #include "IDLList.H"
  35
  36 #include "addToRunTimeSelectionTable.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42     defineTypeNameAndDebug(faceOnlySet, 0);
  43     addToRunTimeSelectionTable(sampledSet, faceOnlySet, word);
  44 }
  45
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49
  50 // Sample singly connected segment. Returns false if end_ reached.
  51 bool Foam::faceOnlySet::trackToBoundary
  52 (
  53     Particle<passiveParticle>& singleParticle,
  54     DynamicList<point>& samplingPts,
  55     DynamicList<label>& samplingCells,
  56     DynamicList<label>& samplingFaces,
  57     DynamicList<scalar>& samplingCurveDist
  58 ) const
  59 {
  60     // distance vector between sampling points
  61     const vector offset = end_ - start_;
  62     const vector smallVec = tol*offset;
  63     const scalar smallDist = mag(smallVec);
  64
  65     // Alias
  66     const point& trackPt = singleParticle.position();
  67
  68     while(true)
  69     {
  70         point oldPoint = trackPt;
  71
  72         singleParticle.trackToFace(end_);
  73
  74         if (singleParticle.face() != -1 && mag(oldPoint - trackPt) > smallDist)
  75         {
  76             // Reached face. Sample.
  77             samplingPts.append(trackPt);
  78             samplingCells.append(singleParticle.cell());
  79             samplingFaces.append(singleParticle.face());
  80             samplingCurveDist.append(mag(trackPt - start_));
  81         }
  82
  83         if (mag(trackPt - end_) < smallDist)
  84         {
  85             // end reached
  86             return false;
  87         }
  88         else if (singleParticle.onBoundary())
  89         {
  90             // Boundary reached.
  91             return true;
  92         }
  93     }
  94 }
  95
  96
  97 void Foam::faceOnlySet::calcSamples
  98 (
  99     DynamicList<point>& samplingPts,
 100     DynamicList<label>& samplingCells,
 101     DynamicList<label>& samplingFaces,
 102     DynamicList<label>& samplingSegments,
 103     DynamicList<scalar>& samplingCurveDist
 104 ) const
 105 {
 106     // distance vector between sampling points
 107     if (mag(end_ - start_) < SMALL)
 108     {
 109         FatalErrorIn("faceOnlySet::calcSamples()")
 110             << "Incorrect sample specification :"
 111             << " start equals end point." << endl
 112             << "  start:" << start_
 113             << "  end:" << end_
 114             << exit(FatalError);
 115     }
 116
 117     const vector offset = (end_ - start_);
 118     const vector normOffset = offset/mag(offset);
 119     const vector smallVec = tol*offset;
 120     const scalar smallDist = mag(smallVec);
 121
 122
 123     // Get all boundary intersections
 124     List<pointIndexHit> bHits = searchEngine().intersections
 125     (
 126         start_ - smallVec,
 127         end_ + smallVec
 128     );
 129
 130     point bPoint(GREAT, GREAT, GREAT);
 131     label bFaceI = -1;
 132
 133     if (bHits.size())
 134     {
 135         bPoint = bHits[0].hitPoint();
 136         bFaceI = bHits[0].index();
 137     }
 138
 139     // Get first tracking point. Use bPoint, bFaceI if provided.
 140
 141     point trackPt;
 142     label trackCellI = -1;
 143     label trackFaceI = -1;
 144
 145     //Info<< "before getTrackingPoint : bPoint:" << bPoint
 146     //    << " bFaceI:" << bFaceI << endl;
 147
 148     getTrackingPoint
 149     (
 150         offset,
 151         start_,
 152         bPoint,
 153         bFaceI,
 154
 155         trackPt,
 156         trackCellI,
 157         trackFaceI
 158     );
 159
 160     //Info<< "after getTrackingPoint : "
 161     //    << " trackPt:" << trackPt
 162     //    << " trackCellI:" << trackCellI
 163     //    << " trackFaceI:" << trackFaceI
 164     //    << endl;
 165
 166     if (trackCellI == -1)
 167     {
 168         // Line start_ - end_ does not intersect domain at all.
 169         // (or is along edge)
 170         // Set points and cell/face labels to empty lists
 171         //Info<< "calcSamples : Both start_ and end_ outside domain"
 172         //    << endl;
 173
 174         return;
 175     }
 176
 177     if (trackFaceI == -1)
 178     {
 179         // No boundary face. Check for nearish internal face
 180         trackFaceI = findNearFace(trackCellI, trackPt, smallDist);
 181     }
 182
 183     //Info<< "calcSamples : got first point to track from :"
 184     //    << "  trackPt:" << trackPt
 185     //    << "  trackCell:" << trackCellI
 186     //    << "  trackFace:" << trackFaceI
 187     //    << endl;
 188
 189     //
 190     // Track until hit end of all boundary intersections
 191     //
 192
 193     // current segment number
 194     label segmentI = 0;
 195
 196     // starting index of current segment in samplePts
 197     label startSegmentI = 0;
 198
 199     // index in bHits; current boundary intersection
 200     label bHitI = 1;
 201
 202     while(true)
 203     {
 204         if (trackFaceI != -1)
 205         {
 206             //Info<< "trackPt:" << trackPt << " on face so use." << endl;
 207             samplingPts.append(trackPt);
 208             samplingCells.append(trackCellI);
 209             samplingFaces.append(trackFaceI);
 210             samplingCurveDist.append(mag(trackPt - start_));
 211         }
 212
 213         // Initialize tracking starting from trackPt
 214         Cloud<passiveParticle> particles(mesh(), IDLList<passiveParticle>());
 215
 216         passiveParticle singleParticle
 217         (
 218             particles,
 219             trackPt,
 220             trackCellI
 221         );
 222
 223         bool reachedBoundary = trackToBoundary
 224         (
 225             singleParticle,
 226             samplingPts,
 227             samplingCells,
 228             samplingFaces,
 229             samplingCurveDist
 230         );
 231
 232         // fill sampleSegments
 233         for(label i = samplingPts.size() - 1; i >= startSegmentI; --i)
 234         {
 235             samplingSegments.append(segmentI);
 236         }
 237
 238
 239         if (!reachedBoundary)
 240         {
 241             //Info<< "calcSamples : Reached end of samples: "
 242             //    << "  samplePt now:" << singleParticle.position()
 243             //    << endl;
 244             break;
 245         }
 246
 247
 248         // Go past boundary intersection where tracking stopped
 249         // Use coordinate comparison instead of face comparison for
 250         // accuracy reasons
 251
 252         bool foundValidB = false;
 253
 254         while (bHitI < bHits.size())
 255         {
 256             scalar dist =
 257                 (bHits[bHitI].hitPoint() - singleParticle.position())
 258               & normOffset;
 259
 260             //Info<< "Finding next boundary : "
 261             //    << "bPoint:" << bHits[bHitI].hitPoint()
 262             //    << "  tracking:" << singleParticle.position()
 263             //    << "  dist:" << dist
 264             //    << endl;
 265
 266             if (dist > smallDist)
 267             {
 268                 // hitpoint is past tracking position
 269                 foundValidB = true;
 270                 break;
 271             }
 272             else
 273             {
 274                 bHitI++;
 275             }
 276         }
 277
 278         if (!foundValidB)
 279         {
 280             // No valid boundary intersection found beyond tracking position
 281             break;
 282         }
 283
 284         // Update starting point for tracking
 285         trackFaceI = bHits[bHitI].index();
 286         trackPt = pushIn(bHits[bHitI].hitPoint(), trackFaceI);
 287         trackCellI = getBoundaryCell(trackFaceI);
 288
 289         segmentI++;
 290
 291         startSegmentI = samplingPts.size();
 292     }
 293 }
 294
 295
 296 void Foam::faceOnlySet::genSamples()
 297 {
 298     // Storage for sample points
 299     DynamicList<point> samplingPts;
 300     DynamicList<label> samplingCells;
 301     DynamicList<label> samplingFaces;
 302     DynamicList<label> samplingSegments;
 303     DynamicList<scalar> samplingCurveDist;
 304
 305     calcSamples
 306     (
 307         samplingPts,
 308         samplingCells,
 309         samplingFaces,
 310         samplingSegments,
 311         samplingCurveDist
 312     );
 313
 314     samplingPts.shrink();
 315     samplingCells.shrink();
 316     samplingFaces.shrink();
 317     samplingSegments.shrink();
 318     samplingCurveDist.shrink();
 319
 320     // Copy into *this
 321     setSamples
 322     (
 323         samplingPts,
 324         samplingCells,
 325         samplingFaces,
 326         samplingSegments,
 327         samplingCurveDist
 328     );
 329 }
 330
 331
 332 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 333
 334 Foam::faceOnlySet::faceOnlySet
 335 (
 336     const word& name,
 337     const polyMesh& mesh,
 338     meshSearch& searchEngine,
 339     const word& axis,
 340     const point& start,
 341     const point& end
 342 )
 343 :
 344     sampledSet(name, mesh, searchEngine, axis),
 345     start_(start),
 346     end_(end)
 347 {
 348     genSamples();
 349
 350     if (debug)
 351     {
 352         write(Info);
 353     }
 354 }
 355
 356
 357 Foam::faceOnlySet::faceOnlySet
 358 (
 359     const word& name,
 360     const polyMesh& mesh,
 361     meshSearch& searchEngine,
 362     const dictionary& dict
 363 )
 364 :
 365     sampledSet(name, mesh, searchEngine, dict),
 366     start_(dict.lookup("start")),
 367     end_(dict.lookup("end"))
 368 {
 369     genSamples();
 370
 371     if (debug)
 372     {
 373         write(Info);
 374     }
 375 }
 376
 377
 378 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 379
 380 Foam::faceOnlySet::~faceOnlySet()
 381 {}
 382
 383
 384 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 385
 386 Foam::point Foam::faceOnlySet::getRefPoint(const List<point>& pts) const
 387 {
 388     return start_;
 389 }
 390
 391
 392 // ************************************************************************* //