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