polyLineSet.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------* \
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 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 "polyLineSet.H"
27 #include "meshSearch.H"
28 #include "DynamicList.H"
29 #include "polyMesh.H"
30 
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37  defineTypeNameAndDebug(polyLineSet, 0);
38  addToRunTimeSelectionTable(sampledSet, polyLineSet, word);
39 
40  const scalar polyLineSet::tol = 1e-6;
41 }
42 
43 
44 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
45 
46 bool Foam::polyLineSet::trackToBoundary
47 (
48  passiveParticleCloud& particleCloud,
49  passiveParticle& singleParticle,
50  label& sampleI,
51  DynamicList<point>& samplingPts,
52  DynamicList<label>& samplingCells,
53  DynamicList<label>& samplingFaces,
54  DynamicList<scalar>& samplingCurveDist
55 ) const
56 {
57  particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
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], trackData);
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::polyLineSet::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  {
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  {
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  // Force calculation of cloud addressing on all processors
163  const bool oldMoving = const_cast<polyMesh&>(mesh()).moving(false);
164  passiveParticleCloud particleCloud(mesh());
165 
166  // current segment number
167  label segmentI = 0;
168 
169  // starting index of current segment in samplePts
170  label startSegmentI = 0;
171 
172  label sampleI = 0;
173 
174  point lastSample(GREAT, GREAT, GREAT);
175  while (true)
176  {
177  // Get boundary intersection
178  point trackPt;
179  label trackCelli = -1;
180  label trackFacei = -1;
181 
182  do
183  {
184  const vector offset =
185  sampleCoords_[sampleI+1] - sampleCoords_[sampleI];
186  const scalar smallDist = mag(tol*offset);
187 
188 
189  // Get all boundary intersections
190  List<pointIndexHit> bHits = searchEngine().intersections
191  (
192  sampleCoords_[sampleI],
193  sampleCoords_[sampleI+1]
194  );
195 
196  point bPoint(GREAT, GREAT, GREAT);
197  label bFacei = -1;
198 
199  if (bHits.size())
200  {
201  bPoint = bHits[0].hitPoint();
202  bFacei = bHits[0].index();
203  }
204 
205  // Get tracking point
206 
207  bool isSample =
208  getTrackingPoint
209  (
210  sampleCoords_[sampleI],
211  bPoint,
212  bFacei,
213  smallDist,
214 
215  trackPt,
216  trackCelli,
217  trackFacei
218  );
219 
220  if (isSample && (mag(lastSample - trackPt) > smallDist))
221  {
222  //Info<< "calcSamples : getTrackingPoint returned valid sample "
223  // << " trackPt:" << trackPt
224  // << " trackFacei:" << trackFacei
225  // << " trackCelli:" << trackCelli
226  // << " sampleI:" << sampleI
227  // << " dist:" << dist
228  // << endl;
229 
230  samplingPts.append(trackPt);
231  samplingCells.append(trackCelli);
232  samplingFaces.append(trackFacei);
233 
234  // Convert sampling position to unique curve parameter. Get
235  // fraction of distance between sampleI and sampleI+1.
236  scalar dist =
237  mag(trackPt - sampleCoords_[sampleI])
238  / mag(sampleCoords_[sampleI+1] - sampleCoords_[sampleI]);
239  samplingCurveDist.append(sampleI + dist);
240 
241  lastSample = trackPt;
242  }
243 
244  if (trackCelli == -1)
245  {
246  // No intersection found. Go to next point
247  sampleI++;
248  }
249  } while ((trackCelli == -1) && (sampleI < sampleCoords_.size() - 1));
250 
251  if (sampleI == sampleCoords_.size() - 1)
252  {
253  //Info<< "calcSamples : Reached end of samples: "
254  // << " sampleI now:" << sampleI
255  // << endl;
256  break;
257  }
258 
259  //
260  // Segment sampleI .. sampleI+1 intersected by domain
261  //
262 
263  // Initialize tracking starting from sampleI
264  passiveParticle singleParticle
265  (
266  mesh(),
267  trackPt,
268  trackCelli
269  );
270 
271  bool bReached = trackToBoundary
272  (
273  particleCloud,
274  singleParticle,
275  sampleI,
276  samplingPts,
277  samplingCells,
278  samplingFaces,
279  samplingCurveDist
280  );
281 
282  // fill sampleSegments
283  for (label i = samplingPts.size() - 1; i >= startSegmentI; --i)
284  {
285  samplingSegments.append(segmentI);
286  }
287 
288  if (!bReached)
289  {
290  //Info<< "calcSamples : Reached end of samples: "
291  // << " sampleI now:" << sampleI
292  // << endl;
293  break;
294  }
295  lastSample = singleParticle.position();
296 
297 
298  // Find next boundary.
299  sampleI++;
300 
301  if (sampleI == sampleCoords_.size() - 1)
302  {
303  //Info<< "calcSamples : Reached end of samples: "
304  // << " sampleI now:" << sampleI
305  // << endl;
306  break;
307  }
308 
309  segmentI++;
310 
311  startSegmentI = samplingPts.size();
312  }
313 
314  const_cast<polyMesh&>(mesh()).moving(oldMoving);
315 }
316 
317 
318 void Foam::polyLineSet::genSamples()
319 {
320  // Storage for sample points
321  DynamicList<point> samplingPts;
322  DynamicList<label> samplingCells;
323  DynamicList<label> samplingFaces;
324  DynamicList<label> samplingSegments;
325  DynamicList<scalar> samplingCurveDist;
326 
327  calcSamples
328  (
329  samplingPts,
330  samplingCells,
331  samplingFaces,
332  samplingSegments,
333  samplingCurveDist
334  );
335 
336  samplingPts.shrink();
337  samplingCells.shrink();
338  samplingFaces.shrink();
339  samplingSegments.shrink();
340  samplingCurveDist.shrink();
341 
342  setSamples
343  (
344  samplingPts,
345  samplingCells,
346  samplingFaces,
347  samplingSegments,
348  samplingCurveDist
349  );
350 }
351 
352 
353 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
354 
356 (
357  const word& name,
358  const polyMesh& mesh,
359  const meshSearch& searchEngine,
360  const word& axis,
361  const List<point>& sampleCoords
362 )
363 :
364  sampledSet(name, mesh, searchEngine, axis),
365  sampleCoords_(sampleCoords)
366 {
367  genSamples();
368 
369  if (debug)
370  {
371  write(Info);
372  }
373 }
374 
375 
377 (
378  const word& name,
379  const polyMesh& mesh,
380  const meshSearch& searchEngine,
381  const dictionary& dict
382 )
383 :
384  sampledSet(name, mesh, searchEngine, dict),
385  sampleCoords_(dict.lookup("points"))
386 {
387  genSamples();
388 
389  if (debug)
390  {
391  write(Info);
392  }
393 }
394 
395 
396 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
397 
399 {}
400 
401 
402 // ************************************************************************* //
Various (local, not parallel) searches on polyMesh; uses (demand driven) octree to search...
Definition: meshSearch.H:57
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
List< pointIndexHit > intersections(const point &pStart, const point &pEnd) const
Find all intersections of boundary within segment pStart .. pEnd.
Definition: meshSearch.C:891
const double e
Elementary charge.
Definition: doubleFloat.H:78
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
virtual ~polyLineSet()
Destructor.
Definition: polyLineSet.C:398
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Macros for easy insertion into run-time selection tables.
dynamicFvMesh & mesh
static const scalar tol
Tolerance when comparing points relative to difference between.
Definition: polyLineSet.H:101
Holds list of sampling points which is filled at construction time. Various implementations of this b...
Definition: sampledSet.H:64
A class for handling words, derived from string.
Definition: word.H:59
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
vector point
Point is a vector.
Definition: point.H:41
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
polyLineSet(const word &name, const polyMesh &mesh, const meshSearch &searchEngine, const word &axis, const List< point > &samplePoints)
Construct from components.
Definition: polyLineSet.C:356
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
runTime write()
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:451