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