src/sampling/sampledSet/polyLine/polyLineSet.C

   1 /*---------------------------------------------------------------------------* \
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2011 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
  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 "polyLineSet.H"
  27 #include "meshSearch.H"
  28 #include "DynamicList.H"
  29 #include "polyMesh.H"
  30
  31 #include "addToRunTimeSelectionTable.H"
  32
  33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37     defineTypeNameAndDebug(polyLineSet, 0);
  38     addToRunTimeSelectionTable(sampledSet, polyLineSet, word);
  39 }
  40
  41
  42 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  43
  44 bool Foam::polyLineSet::trackToBoundary
  45 (
  46     passiveParticle& singleParticle,
  47     label& sampleI,
  48     DynamicList<point>& samplingPts,
  49     DynamicList<label>& samplingCells,
  50     DynamicList<label>& samplingFaces,
  51     DynamicList<scalar>& samplingCurveDist
  52 ) const
  53 {
  54     passiveParticleCloud particleCloud(mesh());
  55     particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
  56
  57     // Alias
  58     const point& trackPt = singleParticle.position();
  59
  60     while (true)
  61     {
  62         // Local geometry info
  63         const vector offset = sampleCoords_[sampleI+1] - sampleCoords_[sampleI];
  64         const scalar smallDist = mag(tol*offset);
  65
  66         point oldPos = trackPt;
  67         label facei = -1;
  68         do
  69         {
  70             singleParticle.stepFraction() = 0;
  71             singleParticle.track(sampleCoords_[sampleI+1], trackData);
  72         }
  73         while
  74         (
  75             !singleParticle.onBoundary()
  76          && (mag(trackPt - oldPos) < smallDist)
  77         );
  78
  79         if (singleParticle.onBoundary())
  80         {
  81             //Info<< "trackToBoundary : reached boundary"
  82             //    << "  trackPt:" << trackPt << endl;
  83             if
  84             (
  85                 mag(trackPt - sampleCoords_[sampleI+1])
  86               < smallDist
  87             )
  88             {
  89                 // Reached samplePt on boundary
  90                 //Info<< "trackToBoundary : boundary. also sampling."
  91                 //    << "  trackPt:" << trackPt << " sampleI+1:" << sampleI+1
  92                 //    << endl;
  93                 samplingPts.append(trackPt);
  94                 samplingCells.append(singleParticle.cell());
  95                 samplingFaces.append(facei);
  96
  97                 // trackPt is at sampleI+1
  98                 samplingCurveDist.append(1.0*(sampleI+1));
  99             }
 100             return true;
 101         }
 102
 103         // Reached samplePt in cell
 104         samplingPts.append(trackPt);
 105         samplingCells.append(singleParticle.cell());
 106         samplingFaces.append(-1);
 107
 108         // Convert trackPt to fraction inbetween sampleI and sampleI+1
 109         scalar dist =
 110             mag(trackPt - sampleCoords_[sampleI])
 111           / mag(sampleCoords_[sampleI+1] - sampleCoords_[sampleI]);
 112         samplingCurveDist.append(sampleI + dist);
 113
 114         // go to next samplePt
 115         sampleI++;
 116
 117         if (sampleI == sampleCoords_.size() - 1)
 118         {
 119             // no more samples.
 120             //Info<< "trackToBoundary : Reached end : sampleI now:" << sampleI
 121             //    << endl;
 122             return false;
 123         }
 124     }
 125 }
 126
 127
 128 void Foam::polyLineSet::calcSamples
 129 (
 130     DynamicList<point>& samplingPts,
 131     DynamicList<label>& samplingCells,
 132     DynamicList<label>& samplingFaces,
 133     DynamicList<label>& samplingSegments,
 134     DynamicList<scalar>& samplingCurveDist
 135 ) const
 136 {
 137     // Check sampling points
 138     if (sampleCoords_.size() < 2)
 139     {
 140         FatalErrorIn("polyLineSet::calcSamples()")
 141             << "Incorrect sample specification. Too few points:"
 142             << sampleCoords_ << exit(FatalError);
 143     }
 144     point oldPoint = sampleCoords_[0];
 145     for (label sampleI = 1; sampleI < sampleCoords_.size(); sampleI++)
 146     {
 147         if (mag(sampleCoords_[sampleI] - oldPoint) < SMALL)
 148         {
 149             FatalErrorIn("polyLineSet::calcSamples()")
 150                 << "Incorrect sample specification."
 151                 << " Point " << sampleCoords_[sampleI-1]
 152                 << " at position " << sampleI-1
 153                 << " and point " << sampleCoords_[sampleI]
 154                 << " at position " << sampleI
 155                 << " are too close" << exit(FatalError);
 156         }
 157         oldPoint = sampleCoords_[sampleI];
 158     }
 159
 160     // Force calculation of minimum-tet decomposition.
 161     (void) mesh().tetBasePtIs();
 162
 163     // current segment number
 164     label segmentI = 0;
 165
 166     // starting index of current segment in samplePts
 167     label startSegmentI = 0;
 168
 169     label sampleI = 0;
 170
 171     point lastSample(GREAT, GREAT, GREAT);
 172     while (true)
 173     {
 174         // Get boundary intersection
 175         point trackPt;
 176         label trackCellI = -1;
 177         label trackFaceI = -1;
 178
 179         do
 180         {
 181             const vector offset =
 182                 sampleCoords_[sampleI+1] - sampleCoords_[sampleI];
 183             const scalar smallDist = mag(tol*offset);
 184
 185
 186             // Get all boundary intersections
 187             List<pointIndexHit> bHits = searchEngine().intersections
 188             (
 189                 sampleCoords_[sampleI],
 190                 sampleCoords_[sampleI+1]
 191             );
 192
 193             point bPoint(GREAT, GREAT, GREAT);
 194             label bFaceI = -1;
 195
 196             if (bHits.size())
 197             {
 198                 bPoint = bHits[0].hitPoint();
 199                 bFaceI = bHits[0].index();
 200             }
 201
 202             // Get tracking point
 203
 204             bool isSample =
 205                 getTrackingPoint
 206                 (
 207                     sampleCoords_[sampleI+1] - sampleCoords_[sampleI],
 208                     sampleCoords_[sampleI],
 209                     bPoint,
 210                     bFaceI,
 211
 212                     trackPt,
 213                     trackCellI,
 214                     trackFaceI
 215                 );
 216
 217             if (isSample && (mag(lastSample - trackPt) > smallDist))
 218             {
 219                 //Info<< "calcSamples : getTrackingPoint returned valid sample "
 220                 //    << "  trackPt:" << trackPt
 221                 //    << "  trackFaceI:" << trackFaceI
 222                 //    << "  trackCellI:" << trackCellI
 223                 //    << "  sampleI:" << sampleI
 224                 //    << "  dist:" << dist
 225                 //    << endl;
 226
 227                 samplingPts.append(trackPt);
 228                 samplingCells.append(trackCellI);
 229                 samplingFaces.append(trackFaceI);
 230
 231                 // Convert sampling position to unique curve parameter. Get
 232                 // fraction of distance between sampleI and sampleI+1.
 233                 scalar dist =
 234                     mag(trackPt - sampleCoords_[sampleI])
 235                   / mag(sampleCoords_[sampleI+1] - sampleCoords_[sampleI]);
 236                 samplingCurveDist.append(sampleI + dist);
 237
 238                 lastSample = trackPt;
 239             }
 240
 241             if (trackCellI == -1)
 242             {
 243                 // No intersection found. Go to next point
 244                 sampleI++;
 245             }
 246         } while ((trackCellI == -1) && (sampleI < sampleCoords_.size() - 1));
 247
 248         if (sampleI == sampleCoords_.size() - 1)
 249         {
 250             //Info<< "calcSamples : Reached end of samples: "
 251             //    << "  sampleI now:" << sampleI
 252             //    << endl;
 253             break;
 254         }
 255
 256         //
 257         // Segment sampleI .. sampleI+1 intersected by domain
 258         //
 259
 260         // Initialize tracking starting from sampleI
 261         passiveParticle singleParticle
 262         (
 263             mesh(),
 264             trackPt,
 265             trackCellI
 266         );
 267
 268         bool bReached = trackToBoundary
 269         (
 270             singleParticle,
 271             sampleI,
 272             samplingPts,
 273             samplingCells,
 274             samplingFaces,
 275             samplingCurveDist
 276         );
 277
 278         // fill sampleSegments
 279         for (label i = samplingPts.size() - 1; i >= startSegmentI; --i)
 280         {
 281             samplingSegments.append(segmentI);
 282         }
 283
 284         if (!bReached)
 285         {
 286             //Info<< "calcSamples : Reached end of samples: "
 287             //    << "  sampleI now:" << sampleI
 288             //    << endl;
 289             break;
 290         }
 291         lastSample = singleParticle.position();
 292
 293
 294         // Find next boundary.
 295         sampleI++;
 296
 297         if (sampleI == sampleCoords_.size() - 1)
 298         {
 299             //Info<< "calcSamples : Reached end of samples: "
 300             //    << "  sampleI now:" << sampleI
 301             //    << endl;
 302             break;
 303         }
 304
 305         segmentI++;
 306
 307         startSegmentI = samplingPts.size();
 308     }
 309 }
 310
 311
 312 void Foam::polyLineSet::genSamples()
 313 {
 314     // Storage for sample points
 315     DynamicList<point> samplingPts;
 316     DynamicList<label> samplingCells;
 317     DynamicList<label> samplingFaces;
 318     DynamicList<label> samplingSegments;
 319     DynamicList<scalar> samplingCurveDist;
 320
 321     calcSamples
 322     (
 323         samplingPts,
 324         samplingCells,
 325         samplingFaces,
 326         samplingSegments,
 327         samplingCurveDist
 328     );
 329
 330     samplingPts.shrink();
 331     samplingCells.shrink();
 332     samplingFaces.shrink();
 333     samplingSegments.shrink();
 334     samplingCurveDist.shrink();
 335
 336     setSamples
 337     (
 338         samplingPts,
 339         samplingCells,
 340         samplingFaces,
 341         samplingSegments,
 342         samplingCurveDist
 343     );
 344 }
 345
 346
 347 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 348
 349 Foam::polyLineSet::polyLineSet
 350 (
 351     const word& name,
 352     const polyMesh& mesh,
 353     meshSearch& searchEngine,
 354     const word& axis,
 355     const List<point>& sampleCoords
 356 )
 357 :
 358     sampledSet(name, mesh, searchEngine, axis),
 359     sampleCoords_(sampleCoords)
 360 {
 361     genSamples();
 362
 363     if (debug)
 364     {
 365         write(Info);
 366     }
 367 }
 368
 369
 370 Foam::polyLineSet::polyLineSet
 371 (
 372     const word& name,
 373     const polyMesh& mesh,
 374     meshSearch& searchEngine,
 375     const dictionary& dict
 376 )
 377 :
 378     sampledSet(name, mesh, searchEngine, dict),
 379     sampleCoords_(dict.lookup("points"))
 380 {
 381     genSamples();
 382
 383     if (debug)
 384     {
 385         write(Info);
 386     }
 387 }
 388
 389
 390 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 391
 392 Foam::polyLineSet::~polyLineSet()
 393 {}
 394
 395
 396 // ************************************************************************* //