src/sampling/sampledSet/uniform/uniformSet.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 "uniformSet.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(uniformSet, 0);
  38     addToRunTimeSelectionTable(sampledSet, uniformSet, word);
  39 }
  40
  41
  42 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  43
  44 bool Foam::uniformSet::nextSample
  45 (
  46     const point& currentPt,
  47     const vector& offset,
  48     const scalar smallDist,
  49     point& samplePt,
  50     label& sampleI
  51 ) const
  52 {
  53     bool pointFound = false;
  54
  55     const vector normOffset = offset/mag(offset);
  56
  57     samplePt += offset;
  58     sampleI++;
  59
  60     for (; sampleI < nPoints_; sampleI++)
  61     {
  62         scalar s = (samplePt - currentPt) & normOffset;
  63
  64         if (s > -smallDist)
  65         {
  66             // samplePt is close to or beyond currentPt -> use it
  67             pointFound = true;
  68
  69             break;
  70         }
  71         samplePt += offset;
  72     }
  73
  74     return pointFound;
  75 }
  76
  77
  78 bool Foam::uniformSet::trackToBoundary
  79 (
  80     passiveParticle& singleParticle,
  81     point& samplePt,
  82     label& sampleI,
  83     DynamicList<point>& samplingPts,
  84     DynamicList<label>& samplingCells,
  85     DynamicList<label>& samplingFaces,
  86     DynamicList<scalar>& samplingCurveDist
  87 ) const
  88 {
  89     // distance vector between sampling points
  90     const vector offset = (end_ - start_)/(nPoints_ - 1);
  91     const vector smallVec = tol*offset;
  92     const scalar smallDist = mag(smallVec);
  93
  94     // Alias
  95     const point& trackPt = singleParticle.position();
  96
  97     passiveParticleCloud particleCloud(mesh());
  98     particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
  99
 100     while(true)
 101     {
 102         // Find next samplePt on/after trackPt. Update samplePt, sampleI
 103         if (!nextSample(trackPt, offset, smallDist, samplePt, sampleI))
 104         {
 105             // no more samples.
 106             if (debug)
 107             {
 108                 Info<< "trackToBoundary : Reached end : samplePt now:"
 109                     << samplePt << "  sampleI now:" << sampleI << endl;
 110             }
 111             return false;
 112         }
 113
 114         if (mag(samplePt - trackPt) < smallDist)
 115         {
 116             // trackPt corresponds with samplePt. Store and use next samplePt
 117             if (debug)
 118             {
 119                 Info<< "trackToBoundary : samplePt corresponds to trackPt : "
 120                     << "  trackPt:" << trackPt << "  samplePt:" << samplePt
 121                     << endl;
 122             }
 123
 124             samplingPts.append(trackPt);
 125             samplingCells.append(singleParticle.cell());
 126             samplingFaces.append(-1);
 127             samplingCurveDist.append(mag(trackPt - start_));
 128
 129             // go to next samplePt
 130             if (!nextSample(trackPt, offset, smallDist, samplePt, sampleI))
 131             {
 132                 // no more samples.
 133                 if (debug)
 134                 {
 135                     Info<< "trackToBoundary : Reached end : "
 136                         << "  samplePt now:" << samplePt
 137                         << "  sampleI now:" << sampleI
 138                         << endl;
 139                 }
 140
 141                 return false;
 142             }
 143         }
 144
 145
 146         if (debug)
 147         {
 148             Info<< "Searching along trajectory from "
 149                 << "  trackPt:" << trackPt
 150                 << "  trackCellI:" << singleParticle.cell()
 151                 << "  to:" << samplePt << endl;
 152         }
 153
 154         point oldPos = trackPt;
 155         label facei = -1;
 156         do
 157         {
 158             singleParticle.stepFraction() = 0;
 159             singleParticle.track(samplePt, trackData);
 160
 161             if (debug)
 162             {
 163                 Info<< "Result of tracking "
 164                     << "  trackPt:" << trackPt
 165                     << "  trackCellI:" << singleParticle.cell()
 166                     << "  trackFaceI:" << singleParticle.face()
 167                     << "  onBoundary:" << singleParticle.onBoundary()
 168                     << "  samplePt:" << samplePt
 169                     << "  smallDist:" << smallDist
 170                     << endl;
 171             }
 172         }
 173         while
 174         (
 175             !singleParticle.onBoundary()
 176          && (mag(trackPt - oldPos) < smallDist)
 177         );
 178
 179         if (singleParticle.onBoundary())
 180         {
 181             //Info<< "trackToBoundary : reached boundary" << endl;
 182             if (mag(trackPt - samplePt) < smallDist)
 183             {
 184                 //Info<< "trackToBoundary : boundary is also sampling point"
 185                 //    << endl;
 186                 // Reached samplePt on boundary
 187                 samplingPts.append(trackPt);
 188                 samplingCells.append(singleParticle.cell());
 189                 samplingFaces.append(facei);
 190                 samplingCurveDist.append(mag(trackPt - start_));
 191             }
 192
 193             return true;
 194         }
 195
 196         //Info<< "trackToBoundary : reached internal sampling point" << endl;
 197         // Reached samplePt in cell or on internal face
 198         samplingPts.append(trackPt);
 199         samplingCells.append(singleParticle.cell());
 200         samplingFaces.append(-1);
 201         samplingCurveDist.append(mag(trackPt - start_));
 202
 203         // go to next samplePt
 204     }
 205 }
 206
 207
 208 void Foam::uniformSet::calcSamples
 209 (
 210     DynamicList<point>& samplingPts,
 211     DynamicList<label>& samplingCells,
 212     DynamicList<label>& samplingFaces,
 213     DynamicList<label>& samplingSegments,
 214     DynamicList<scalar>& samplingCurveDist
 215 ) const
 216 {
 217     // distance vector between sampling points
 218     if ((nPoints_ < 2) || (mag(end_ - start_) < SMALL))
 219     {
 220         FatalErrorIn("uniformSet::calcSamples()")
 221             << "Incorrect sample specification. Either too few points or"
 222             << " start equals end point." << endl
 223             << "nPoints:" << nPoints_
 224             << "  start:" << start_
 225             << "  end:" << end_
 226             << exit(FatalError);
 227     }
 228
 229     const vector offset = (end_ - start_)/(nPoints_ - 1);
 230     const vector normOffset = offset/mag(offset);
 231     const vector smallVec = tol*offset;
 232     const scalar smallDist = mag(smallVec);
 233
 234     // Force calculation of minimum-tet decomposition.
 235     (void) mesh().tetBasePtIs();
 236
 237     // Get all boundary intersections
 238     List<pointIndexHit> bHits = searchEngine().intersections
 239     (
 240         start_ - smallVec,
 241         end_ + smallVec
 242     );
 243
 244     point bPoint(GREAT, GREAT, GREAT);
 245     label bFaceI = -1;
 246
 247     if (bHits.size())
 248     {
 249         bPoint = bHits[0].hitPoint();
 250         bFaceI = bHits[0].index();
 251     }
 252
 253     // Get first tracking point. Use bPoint, bFaceI if provided.
 254
 255     point trackPt;
 256     label trackCellI = -1;
 257     label trackFaceI = -1;
 258
 259     bool isSample =
 260         getTrackingPoint
 261         (
 262             offset,
 263             start_,
 264             bPoint,
 265             bFaceI,
 266
 267             trackPt,
 268             trackCellI,
 269             trackFaceI
 270         );
 271
 272     if (trackCellI == -1)
 273     {
 274         // Line start_ - end_ does not intersect domain at all.
 275         // (or is along edge)
 276         // Set points and cell/face labels to empty lists
 277
 278         return;
 279     }
 280
 281     if (isSample)
 282     {
 283         samplingPts.append(start_);
 284         samplingCells.append(trackCellI);
 285         samplingFaces.append(trackFaceI);
 286         samplingCurveDist.append(0.0);
 287     }
 288
 289     //
 290     // Track until hit end of all boundary intersections
 291     //
 292
 293     // current segment number
 294     label segmentI = 0;
 295
 296     // starting index of current segment in samplePts
 297     label startSegmentI = 0;
 298
 299     label sampleI = 0;
 300     point samplePt = start_;
 301
 302     // index in bHits; current boundary intersection
 303     label bHitI = 1;
 304
 305     while(true)
 306     {
 307         // Initialize tracking starting from trackPt
 308         passiveParticle singleParticle(mesh(), trackPt, trackCellI);
 309
 310         bool reachedBoundary = trackToBoundary
 311         (
 312             singleParticle,
 313             samplePt,
 314             sampleI,
 315             samplingPts,
 316             samplingCells,
 317             samplingFaces,
 318             samplingCurveDist
 319         );
 320
 321         // fill sampleSegments
 322         for (label i = samplingPts.size() - 1; i >= startSegmentI; --i)
 323         {
 324             samplingSegments.append(segmentI);
 325         }
 326
 327
 328         if (!reachedBoundary)
 329         {
 330             if (debug)
 331             {
 332                 Info<< "calcSamples : Reached end of samples: "
 333                     << "  samplePt now:" << samplePt
 334                     << "  sampleI now:" << sampleI
 335                     << endl;
 336             }
 337             break;
 338         }
 339
 340
 341         bool foundValidB = false;
 342
 343         while (bHitI < bHits.size())
 344         {
 345             scalar dist =
 346                 (bHits[bHitI].hitPoint() - singleParticle.position())
 347               & normOffset;
 348
 349             if (debug)
 350             {
 351                 Info<< "Finding next boundary : "
 352                     << "bPoint:" << bHits[bHitI].hitPoint()
 353                     << "  tracking:" << singleParticle.position()
 354                     << "  dist:" << dist
 355                     << endl;
 356             }
 357
 358             if (dist > smallDist)
 359             {
 360                 // hitpoint is past tracking position
 361                 foundValidB = true;
 362                 break;
 363             }
 364             else
 365             {
 366                 bHitI++;
 367             }
 368         }
 369
 370         if (!foundValidB)
 371         {
 372             // No valid boundary intersection found beyond tracking position
 373             break;
 374         }
 375
 376         // Update starting point for tracking
 377         trackFaceI = bFaceI;
 378         trackPt = pushIn(bPoint, trackFaceI);
 379         trackCellI = getBoundaryCell(trackFaceI);
 380
 381         segmentI++;
 382
 383         startSegmentI = samplingPts.size();
 384     }
 385 }
 386
 387
 388 void Foam::uniformSet::genSamples()
 389 {
 390     // Storage for sample points
 391     DynamicList<point> samplingPts;
 392     DynamicList<label> samplingCells;
 393     DynamicList<label> samplingFaces;
 394     DynamicList<label> samplingSegments;
 395     DynamicList<scalar> samplingCurveDist;
 396
 397     calcSamples
 398     (
 399         samplingPts,
 400         samplingCells,
 401         samplingFaces,
 402         samplingSegments,
 403         samplingCurveDist
 404     );
 405
 406     samplingPts.shrink();
 407     samplingCells.shrink();
 408     samplingFaces.shrink();
 409     samplingSegments.shrink();
 410     samplingCurveDist.shrink();
 411
 412     setSamples
 413     (
 414         samplingPts,
 415         samplingCells,
 416         samplingFaces,
 417         samplingSegments,
 418         samplingCurveDist
 419     );
 420 }
 421
 422 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 423
 424 Foam::uniformSet::uniformSet
 425 (
 426     const word& name,
 427     const polyMesh& mesh,
 428     meshSearch& searchEngine,
 429     const word& axis,
 430     const point& start,
 431     const point& end,
 432     const label nPoints
 433 )
 434 :
 435     sampledSet(name, mesh, searchEngine, axis),
 436     start_(start),
 437     end_(end),
 438     nPoints_(nPoints)
 439 {
 440     genSamples();
 441
 442     if (debug)
 443     {
 444         write(Info);
 445     }
 446 }
 447
 448
 449 Foam::uniformSet::uniformSet
 450 (
 451     const word& name,
 452     const polyMesh& mesh,
 453     meshSearch& searchEngine,
 454     const dictionary& dict
 455 )
 456 :
 457     sampledSet(name, mesh, searchEngine, dict),
 458     start_(dict.lookup("start")),
 459     end_(dict.lookup("end")),
 460     nPoints_(readLabel(dict.lookup("nPoints")))
 461 {
 462     genSamples();
 463
 464     if (debug)
 465     {
 466         write(Info);
 467     }
 468 }
 469
 470
 471 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 472
 473 Foam::uniformSet::~uniformSet()
 474 {}
 475
 476
 477 // ************************************************************************* //