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