FreeFOAM The Cross-Platform CFD Toolkit
uniformSet.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) 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 "uniformSet.H"
27 #include <meshTools/meshSearch.H>
28 #include <OpenFOAM/DynamicList.H>
29 #include <OpenFOAM/polyMesh.H>
30 
31 #include <lagrangian/Cloud.H>
33 #include <OpenFOAM/IDLList.H>
34 
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41  defineTypeNameAndDebug(uniformSet, 0);
42  addToRunTimeSelectionTable(sampledSet, uniformSet, word);
43 }
44 
45 
46 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
47 
48 // Finds along line (samplePt + t * offset) next sample beyond or equal to
49 // currentPt.
50 // Updates samplePt, sampleI
51 bool Foam::uniformSet::nextSample
52 (
53  const point& currentPt,
54  const vector& offset,
55  const scalar smallDist,
56  point& samplePt,
57  label& sampleI
58 ) const
59 {
60  bool pointFound = false;
61 
62  const vector normOffset = offset/mag(offset);
63 
64  samplePt += offset;
65  sampleI++;
66 
67  for(; sampleI < nPoints_; sampleI++)
68  {
69  scalar s = (samplePt - currentPt) & normOffset;
70 
71  if (s > -smallDist)
72  {
73  // samplePt is close to or beyond currentPt -> use it
74  pointFound = true;
75 
76  break;
77  }
78  samplePt += offset;
79  }
80 
81  return pointFound;
82 }
83 
84 
85 // Sample singly connected segment. Returns false if end_ reached.
86 bool Foam::uniformSet::trackToBoundary
87 (
88  Particle<passiveParticle>& singleParticle,
89  point& samplePt,
90  label& sampleI,
91  DynamicList<point>& samplingPts,
92  DynamicList<label>& samplingCells,
93  DynamicList<label>& samplingFaces,
94  DynamicList<scalar>& samplingCurveDist
95 ) const
96 {
97  // distance vector between sampling points
98  const vector offset = (end_ - start_)/(nPoints_ - 1);
99  const vector smallVec = tol*offset;
100  const scalar smallDist = mag(smallVec);
101 
102  // Alias
103  const point& trackPt = singleParticle.position();
104 
105  while(true)
106  {
107  // Find next samplePt on/after trackPt. Update samplePt, sampleI
108  if (!nextSample(trackPt, offset, smallDist, samplePt, sampleI))
109  {
110  // no more samples.
111  if (debug)
112  {
113  Info<< "trackToBoundary : Reached end : samplePt now:"
114  << samplePt << " sampleI now:" << sampleI << endl;
115  }
116  return false;
117  }
118 
119  if (mag(samplePt - trackPt) < smallDist)
120  {
121  // trackPt corresponds with samplePt. Store and use next samplePt
122  if (debug)
123  {
124  Info<< "trackToBoundary : samplePt corresponds to trackPt : "
125  << " trackPt:" << trackPt << " samplePt:" << samplePt
126  << endl;
127  }
128 
129  samplingPts.append(trackPt);
130  samplingCells.append(singleParticle.cell());
131  samplingFaces.append(-1);
132  samplingCurveDist.append(mag(trackPt - start_));
133 
134  // go to next samplePt
135  if (!nextSample(trackPt, offset, smallDist, samplePt, sampleI))
136  {
137  // no more samples.
138  if (debug)
139  {
140  Info<< "trackToBoundary : Reached end : "
141  << " samplePt now:" << samplePt
142  << " sampleI now:" << sampleI
143  << endl;
144  }
145 
146  return false;
147  }
148  }
149 
150 
151  if (debug)
152  {
153  Info<< "Searching along trajectory from "
154  << " trackPt:" << trackPt
155  << " trackCellI:" << singleParticle.cell()
156  << " to:" << samplePt << endl;
157  }
158 
159  point oldPos = trackPt;
160  label facei = -1;
161  do
162  {
163  singleParticle.stepFraction() = 0;
164  singleParticle.track(samplePt);
165 
166  if (debug)
167  {
168  Info<< "Result of tracking "
169  << " trackPt:" << trackPt
170  << " trackCellI:" << singleParticle.cell()
171  << " trackFaceI:" << singleParticle.face()
172  << " onBoundary:" << singleParticle.onBoundary()
173  << " samplePt:" << samplePt
174  << " smallDist:" << smallDist
175  << endl;
176  }
177  }
178  while
179  (
180  !singleParticle.onBoundary()
181  && (mag(trackPt - oldPos) < smallDist)
182  );
183 
184  if (singleParticle.onBoundary())
185  {
186  //Info<< "trackToBoundary : reached boundary" << endl;
187  if (mag(trackPt - samplePt) < smallDist)
188  {
189  //Info<< "trackToBoundary : boundary is also sampling point"
190  // << endl;
191  // Reached samplePt on boundary
192  samplingPts.append(trackPt);
193  samplingCells.append(singleParticle.cell());
194  samplingFaces.append(facei);
195  samplingCurveDist.append(mag(trackPt - start_));
196  }
197 
198  return true;
199  }
200 
201  //Info<< "trackToBoundary : reached internal sampling point" << endl;
202  // Reached samplePt in cell or on internal face
203  samplingPts.append(trackPt);
204  samplingCells.append(singleParticle.cell());
205  samplingFaces.append(-1);
206  samplingCurveDist.append(mag(trackPt - start_));
207 
208  // go to next samplePt
209  }
210 }
211 
212 
213 void Foam::uniformSet::calcSamples
214 (
215  DynamicList<point>& samplingPts,
216  DynamicList<label>& samplingCells,
217  DynamicList<label>& samplingFaces,
218  DynamicList<label>& samplingSegments,
219  DynamicList<scalar>& samplingCurveDist
220 ) const
221 {
222  // distance vector between sampling points
223  if ((nPoints_ < 2) || (mag(end_ - start_) < SMALL))
224  {
225  FatalErrorIn("uniformSet::calcSamples()")
226  << "Incorrect sample specification. Either too few points or"
227  << " start equals end point." << endl
228  << "nPoints:" << nPoints_
229  << " start:" << start_
230  << " end:" << end_
231  << exit(FatalError);
232  }
233 
234  const vector offset = (end_ - start_)/(nPoints_ - 1);
235  const vector normOffset = offset/mag(offset);
236  const vector smallVec = tol*offset;
237  const scalar smallDist = mag(smallVec);
238 
239  // Get all boundary intersections
240  List<pointIndexHit> bHits = searchEngine().intersections
241  (
242  start_ - smallVec,
243  end_ + smallVec
244  );
245 
246  point bPoint(GREAT, GREAT, GREAT);
247  label bFaceI = -1;
248 
249  if (bHits.size())
250  {
251  bPoint = bHits[0].hitPoint();
252  bFaceI = bHits[0].index();
253  }
254 
255  // Get first tracking point. Use bPoint, bFaceI if provided.
256 
257  point trackPt;
258  label trackCellI = -1;
259  label trackFaceI = -1;
260 
261  bool isSample =
262  getTrackingPoint
263  (
264  offset,
265  start_,
266  bPoint,
267  bFaceI,
268 
269  trackPt,
270  trackCellI,
271  trackFaceI
272  );
273 
274  if (trackCellI == -1)
275  {
276  // Line start_ - end_ does not intersect domain at all.
277  // (or is along edge)
278  // Set points and cell/face labels to empty lists
279 
280  return;
281  }
282 
283  if (isSample)
284  {
285  samplingPts.append(start_);
286  samplingCells.append(trackCellI);
287  samplingFaces.append(trackFaceI);
288  samplingCurveDist.append(0.0);
289  }
290 
291  //
292  // Track until hit end of all boundary intersections
293  //
294 
295  // current segment number
296  label segmentI = 0;
297 
298  // starting index of current segment in samplePts
299  label startSegmentI = 0;
300 
301  label sampleI = 0;
302  point samplePt = start_;
303 
304  // index in bHits; current boundary intersection
305  label bHitI = 1;
306 
307  while(true)
308  {
309  // Initialize tracking starting from trackPt
310  Cloud<passiveParticle> particles(mesh(), IDLList<passiveParticle>());
311 
312  passiveParticle singleParticle
313  (
314  particles,
315  trackPt,
316  trackCellI
317  );
318 
319  bool reachedBoundary = trackToBoundary
320  (
321  singleParticle,
322  samplePt,
323  sampleI,
324  samplingPts,
325  samplingCells,
326  samplingFaces,
327  samplingCurveDist
328  );
329 
330  // fill sampleSegments
331  for(label i = samplingPts.size() - 1; i >= startSegmentI; --i)
332  {
333  samplingSegments.append(segmentI);
334  }
335 
336 
337  if (!reachedBoundary)
338  {
339  if (debug)
340  {
341  Info<< "calcSamples : Reached end of samples: "
342  << " samplePt now:" << samplePt
343  << " sampleI now:" << sampleI
344  << endl;
345  }
346  break;
347  }
348 
349 
350  bool foundValidB = false;
351 
352  while (bHitI < bHits.size())
353  {
354  scalar dist =
355  (bHits[bHitI].hitPoint() - singleParticle.position())
356  & normOffset;
357 
358  if (debug)
359  {
360  Info<< "Finding next boundary : "
361  << "bPoint:" << bHits[bHitI].hitPoint()
362  << " tracking:" << singleParticle.position()
363  << " dist:" << dist
364  << endl;
365  }
366 
367  if (dist > smallDist)
368  {
369  // hitpoint is past tracking position
370  foundValidB = true;
371  break;
372  }
373  else
374  {
375  bHitI++;
376  }
377  }
378 
379  if (!foundValidB)
380  {
381  // No valid boundary intersection found beyond tracking position
382  break;
383  }
384 
385  // Update starting point for tracking
386  trackFaceI = bFaceI;
387  trackPt = pushIn(bPoint, trackFaceI);
388  trackCellI = getBoundaryCell(trackFaceI);
389 
390  segmentI++;
391 
392  startSegmentI = samplingPts.size();
393  }
394 }
395 
396 
397 void Foam::uniformSet::genSamples()
398 {
399  // Storage for sample points
400  DynamicList<point> samplingPts;
401  DynamicList<label> samplingCells;
402  DynamicList<label> samplingFaces;
403  DynamicList<label> samplingSegments;
404  DynamicList<scalar> samplingCurveDist;
405 
406  calcSamples
407  (
408  samplingPts,
409  samplingCells,
410  samplingFaces,
411  samplingSegments,
412  samplingCurveDist
413  );
414 
415  samplingPts.shrink();
416  samplingCells.shrink();
417  samplingFaces.shrink();
418  samplingSegments.shrink();
419  samplingCurveDist.shrink();
420 
421  setSamples
422  (
423  samplingPts,
424  samplingCells,
425  samplingFaces,
426  samplingSegments,
427  samplingCurveDist
428  );
429 }
430 
431 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
432 
434 (
435  const word& name,
436  const polyMesh& mesh,
437  meshSearch& searchEngine,
438  const word& axis,
439  const point& start,
440  const point& end,
441  const label nPoints
442 )
443 :
444  sampledSet(name, mesh, searchEngine, axis),
445  start_(start),
446  end_(end),
447  nPoints_(nPoints)
448 {
449  genSamples();
450 
451  if (debug)
452  {
453  write(Info);
454  }
455 }
456 
457 
459 (
460  const word& name,
461  const polyMesh& mesh,
462  meshSearch& searchEngine,
463  const dictionary& dict
464 )
465 :
466  sampledSet(name, mesh, searchEngine, dict),
467  start_(dict.lookup("start")),
468  end_(dict.lookup("end")),
469  nPoints_(readLabel(dict.lookup("nPoints")))
470 {
471  genSamples();
472 
473  if (debug)
474  {
475  write(Info);
476  }
477 }
478 
479 
480 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
481 
483 {}
484 
485 
486 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
487 
488 
490 {
491  // Use start point as reference for 'distance'
492  return start_;
493 }
494 
495 
496 // ************************ vim: set sw=4 sts=4 et: ************************ //