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