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