BUG: UListIO: byteSize overflowing on really big faceLists
[OpenFOAM-2.0.x.git] / src / sampling / sampledSet / uniform / uniformSet.C
blob84015c7f2fae205d1dc26cf0105d703ade0d9ead
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.
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.
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.
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/>.
24 \*---------------------------------------------------------------------------*/
26 #include "uniformSet.H"
27 #include "meshSearch.H"
28 #include "DynamicList.H"
29 #include "polyMesh.H"
31 #include "addToRunTimeSelectionTable.H"
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 namespace Foam
37     defineTypeNameAndDebug(uniformSet, 0);
38     addToRunTimeSelectionTable(sampledSet, uniformSet, word);
42 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
44 bool Foam::uniformSet::nextSample
46     const point& currentPt,
47     const vector& offset,
48     const scalar smallDist,
49     point& samplePt,
50     label& sampleI
51 ) const
53     bool pointFound = false;
55     const vector normOffset = offset/mag(offset);
57     samplePt += offset;
58     sampleI++;
60     for (; sampleI < nPoints_; sampleI++)
61     {
62         scalar s = (samplePt - currentPt) & normOffset;
64         if (s > -smallDist)
65         {
66             // samplePt is close to or beyond currentPt -> use it
67             pointFound = true;
69             break;
70         }
71         samplePt += offset;
72     }
74     return pointFound;
78 bool Foam::uniformSet::trackToBoundary
80     passiveParticle& singleParticle,
81     point& samplePt,
82     label& sampleI,
83     DynamicList<point>& samplingPts,
84     DynamicList<label>& samplingCells,
85     DynamicList<label>& samplingFaces,
86     DynamicList<scalar>& samplingCurveDist
87 ) const
89     // distance vector between sampling points
90     const vector offset = (end_ - start_)/(nPoints_ - 1);
91     const vector smallVec = tol*offset;
92     const scalar smallDist = mag(smallVec);
94     // Alias
95     const point& trackPt = singleParticle.position();
97     passiveParticleCloud particleCloud(mesh());
98     particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
100     while(true)
101     {
102         // Find next samplePt on/after trackPt. Update samplePt, sampleI
103         if (!nextSample(trackPt, offset, smallDist, samplePt, sampleI))
104         {
105             // no more samples.
106             if (debug)
107             {
108                 Info<< "trackToBoundary : Reached end : samplePt now:"
109                     << samplePt << "  sampleI now:" << sampleI << endl;
110             }
111             return false;
112         }
114         if (mag(samplePt - trackPt) < smallDist)
115         {
116             // trackPt corresponds with samplePt. Store and use next samplePt
117             if (debug)
118             {
119                 Info<< "trackToBoundary : samplePt corresponds to trackPt : "
120                     << "  trackPt:" << trackPt << "  samplePt:" << samplePt
121                     << endl;
122             }
124             samplingPts.append(trackPt);
125             samplingCells.append(singleParticle.cell());
126             samplingFaces.append(-1);
127             samplingCurveDist.append(mag(trackPt - start_));
129             // go to next samplePt
130             if (!nextSample(trackPt, offset, smallDist, samplePt, sampleI))
131             {
132                 // no more samples.
133                 if (debug)
134                 {
135                     Info<< "trackToBoundary : Reached end : "
136                         << "  samplePt now:" << samplePt
137                         << "  sampleI now:" << sampleI
138                         << endl;
139                 }
141                 return false;
142             }
143         }
146         if (debug)
147         {
148             Info<< "Searching along trajectory from "
149                 << "  trackPt:" << trackPt
150                 << "  trackCellI:" << singleParticle.cell()
151                 << "  to:" << samplePt << endl;
152         }
154         point oldPos = trackPt;
155         label facei = -1;
156         do
157         {
158             singleParticle.stepFraction() = 0;
159             singleParticle.track(samplePt, trackData);
161             if (debug)
162             {
163                 Info<< "Result of tracking "
164                     << "  trackPt:" << trackPt
165                     << "  trackCellI:" << singleParticle.cell()
166                     << "  trackFaceI:" << singleParticle.face()
167                     << "  onBoundary:" << singleParticle.onBoundary()
168                     << "  samplePt:" << samplePt
169                     << "  smallDist:" << smallDist
170                     << endl;
171             }
172         }
173         while
174         (
175             !singleParticle.onBoundary()
176          && (mag(trackPt - oldPos) < smallDist)
177         );
179         if (singleParticle.onBoundary())
180         {
181             //Info<< "trackToBoundary : reached boundary" << endl;
182             if (mag(trackPt - samplePt) < smallDist)
183             {
184                 //Info<< "trackToBoundary : boundary is also sampling point"
185                 //    << endl;
186                 // Reached samplePt on boundary
187                 samplingPts.append(trackPt);
188                 samplingCells.append(singleParticle.cell());
189                 samplingFaces.append(facei);
190                 samplingCurveDist.append(mag(trackPt - start_));
191             }
193             return true;
194         }
196         //Info<< "trackToBoundary : reached internal sampling point" << endl;
197         // Reached samplePt in cell or on internal face
198         samplingPts.append(trackPt);
199         samplingCells.append(singleParticle.cell());
200         samplingFaces.append(-1);
201         samplingCurveDist.append(mag(trackPt - start_));
203         // go to next samplePt
204     }
208 void Foam::uniformSet::calcSamples
210     DynamicList<point>& samplingPts,
211     DynamicList<label>& samplingCells,
212     DynamicList<label>& samplingFaces,
213     DynamicList<label>& samplingSegments,
214     DynamicList<scalar>& samplingCurveDist
215 ) const
217     // distance vector between sampling points
218     if ((nPoints_ < 2) || (mag(end_ - start_) < SMALL))
219     {
220         FatalErrorIn("uniformSet::calcSamples()")
221             << "Incorrect sample specification. Either too few points or"
222             << " start equals end point." << endl
223             << "nPoints:" << nPoints_
224             << "  start:" << start_
225             << "  end:" << end_
226             << exit(FatalError);
227     }
229     const vector offset = (end_ - start_)/(nPoints_ - 1);
230     const vector normOffset = offset/mag(offset);
231     const vector smallVec = tol*offset;
232     const scalar smallDist = mag(smallVec);
234     // Force calculation of minimum-tet decomposition.
235     (void) mesh().tetBasePtIs();
237     // Get all boundary intersections
238     List<pointIndexHit> bHits = searchEngine().intersections
239     (
240         start_ - smallVec,
241         end_ + smallVec
242     );
244     point bPoint(GREAT, GREAT, GREAT);
245     label bFaceI = -1;
247     if (bHits.size())
248     {
249         bPoint = bHits[0].hitPoint();
250         bFaceI = bHits[0].index();
251     }
253     // Get first tracking point. Use bPoint, bFaceI if provided.
255     point trackPt;
256     label trackCellI = -1;
257     label trackFaceI = -1;
259     bool isSample =
260         getTrackingPoint
261         (
262             offset,
263             start_,
264             bPoint,
265             bFaceI,
267             trackPt,
268             trackCellI,
269             trackFaceI
270         );
272     if (trackCellI == -1)
273     {
274         // Line start_ - end_ does not intersect domain at all.
275         // (or is along edge)
276         // Set points and cell/face labels to empty lists
278         return;
279     }
281     if (isSample)
282     {
283         samplingPts.append(start_);
284         samplingCells.append(trackCellI);
285         samplingFaces.append(trackFaceI);
286         samplingCurveDist.append(0.0);
287     }
289     //
290     // Track until hit end of all boundary intersections
291     //
293     // current segment number
294     label segmentI = 0;
296     // starting index of current segment in samplePts
297     label startSegmentI = 0;
299     label sampleI = 0;
300     point samplePt = start_;
302     // index in bHits; current boundary intersection
303     label bHitI = 1;
305     while(true)
306     {
307         // Initialize tracking starting from trackPt
308         passiveParticle singleParticle(mesh(), trackPt, trackCellI);
310         bool reachedBoundary = trackToBoundary
311         (
312             singleParticle,
313             samplePt,
314             sampleI,
315             samplingPts,
316             samplingCells,
317             samplingFaces,
318             samplingCurveDist
319         );
321         // fill sampleSegments
322         for (label i = samplingPts.size() - 1; i >= startSegmentI; --i)
323         {
324             samplingSegments.append(segmentI);
325         }
328         if (!reachedBoundary)
329         {
330             if (debug)
331             {
332                 Info<< "calcSamples : Reached end of samples: "
333                     << "  samplePt now:" << samplePt
334                     << "  sampleI now:" << sampleI
335                     << endl;
336             }
337             break;
338         }
341         bool foundValidB = false;
343         while (bHitI < bHits.size())
344         {
345             scalar dist =
346                 (bHits[bHitI].hitPoint() - singleParticle.position())
347               & normOffset;
349             if (debug)
350             {
351                 Info<< "Finding next boundary : "
352                     << "bPoint:" << bHits[bHitI].hitPoint()
353                     << "  tracking:" << singleParticle.position()
354                     << "  dist:" << dist
355                     << endl;
356             }
358             if (dist > smallDist)
359             {
360                 // hitpoint is past tracking position
361                 foundValidB = true;
362                 break;
363             }
364             else
365             {
366                 bHitI++;
367             }
368         }
370         if (!foundValidB)
371         {
372             // No valid boundary intersection found beyond tracking position
373             break;
374         }
376         // Update starting point for tracking
377         trackFaceI = bFaceI;
378         trackPt = pushIn(bPoint, trackFaceI);
379         trackCellI = getBoundaryCell(trackFaceI);
381         segmentI++;
383         startSegmentI = samplingPts.size();
384     }
388 void Foam::uniformSet::genSamples()
390     // Storage for sample points
391     DynamicList<point> samplingPts;
392     DynamicList<label> samplingCells;
393     DynamicList<label> samplingFaces;
394     DynamicList<label> samplingSegments;
395     DynamicList<scalar> samplingCurveDist;
397     calcSamples
398     (
399         samplingPts,
400         samplingCells,
401         samplingFaces,
402         samplingSegments,
403         samplingCurveDist
404     );
406     samplingPts.shrink();
407     samplingCells.shrink();
408     samplingFaces.shrink();
409     samplingSegments.shrink();
410     samplingCurveDist.shrink();
412     setSamples
413     (
414         samplingPts,
415         samplingCells,
416         samplingFaces,
417         samplingSegments,
418         samplingCurveDist
419     );
422 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
424 Foam::uniformSet::uniformSet
426     const word& name,
427     const polyMesh& mesh,
428     meshSearch& searchEngine,
429     const word& axis,
430     const point& start,
431     const point& end,
432     const label nPoints
435     sampledSet(name, mesh, searchEngine, axis),
436     start_(start),
437     end_(end),
438     nPoints_(nPoints)
440     genSamples();
442     if (debug)
443     {
444         write(Info);
445     }
449 Foam::uniformSet::uniformSet
451     const word& name,
452     const polyMesh& mesh,
453     meshSearch& searchEngine,
454     const dictionary& dict
457     sampledSet(name, mesh, searchEngine, dict),
458     start_(dict.lookup("start")),
459     end_(dict.lookup("end")),
460     nPoints_(readLabel(dict.lookup("nPoints")))
462     genSamples();
464     if (debug)
465     {
466         write(Info);
467     }
471 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
473 Foam::uniformSet::~uniformSet()
477 // ************************************************************************* //