FreeFOAM The Cross-Platform CFD Toolkit
curveSet.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 "curveSet.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(curveSet, 0);
42  addToRunTimeSelectionTable(sampledSet, curveSet, word);
43 }
44 
45 
46 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
47 
48 // Sample till hits boundary.
49 bool Foam::curveSet::trackToBoundary
50 (
51  Particle<passiveParticle>& singleParticle,
52  label& sampleI,
53  DynamicList<point>& samplingPts,
54  DynamicList<label>& samplingCells,
55  DynamicList<label>& samplingFaces,
56  DynamicList<scalar>& samplingCurveDist
57 ) const
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]);
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::curveSet::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  {
142  FatalErrorIn("curveSet::calcSamples()")
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  {
151  FatalErrorIn("curveSet::calcSamples()")
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  // current segment number
163  label segmentI = 0;
164 
165  // starting index of current segment in samplePts
166  label startSegmentI = 0;
167 
168  label sampleI = 0;
169 
170  point lastSample(GREAT, GREAT, GREAT);
171  while(true)
172  {
173  // Get boundary intersection
174  point trackPt;
175  label trackCellI = -1;
176  label trackFaceI = -1;
177 
178  do
179  {
180  const vector offset =
181  sampleCoords_[sampleI+1] - sampleCoords_[sampleI];
182  const scalar smallDist = mag(tol*offset);
183 
184 
185  // Get all boundary intersections
186  List<pointIndexHit> bHits = searchEngine().intersections
187  (
188  sampleCoords_[sampleI],
189  sampleCoords_[sampleI+1]
190  );
191 
192  point bPoint(GREAT, GREAT, GREAT);
193  label bFaceI = -1;
194 
195  if (bHits.size())
196  {
197  bPoint = bHits[0].hitPoint();
198  bFaceI = bHits[0].index();
199  }
200 
201  // Get tracking point
202 
203  bool isSample =
204  getTrackingPoint
205  (
206  sampleCoords_[sampleI+1] - sampleCoords_[sampleI],
207  sampleCoords_[sampleI],
208  bPoint,
209  bFaceI,
210 
211  trackPt,
212  trackCellI,
213  trackFaceI
214  );
215 
216  if (isSample && (mag(lastSample - trackPt) > smallDist))
217  {
218  //Info<< "calcSamples : getTrackingPoint returned valid sample "
219  // << " trackPt:" << trackPt
220  // << " trackFaceI:" << trackFaceI
221  // << " trackCellI:" << trackCellI
222  // << " sampleI:" << sampleI
223  // << " dist:" << dist
224  // << endl;
225 
226  samplingPts.append(trackPt);
227  samplingCells.append(trackCellI);
228  samplingFaces.append(trackFaceI);
229 
230  // Convert sampling position to unique curve parameter. Get
231  // fraction of distance between sampleI and sampleI+1.
232  scalar dist =
233  mag(trackPt - sampleCoords_[sampleI])
234  / mag(sampleCoords_[sampleI+1] - sampleCoords_[sampleI]);
235  samplingCurveDist.append(sampleI + dist);
236 
237  lastSample = trackPt;
238  }
239 
240  if (trackCellI == -1)
241  {
242  // No intersection found. Go to next point
243  sampleI++;
244  }
245  } while((trackCellI == -1) && (sampleI < sampleCoords_.size() - 1));
246 
247  if (sampleI == sampleCoords_.size() - 1)
248  {
249  //Info<< "calcSamples : Reached end of samples: "
250  // << " sampleI now:" << sampleI
251  // << endl;
252  break;
253  }
254 
255  //
256  // Segment sampleI .. sampleI+1 intersected by domain
257  //
258 
259  // Initialize tracking starting from sampleI
260  Cloud<passiveParticle> particles(mesh(), IDLList<passiveParticle>());
261 
262  passiveParticle singleParticle
263  (
264  particles,
265  trackPt,
266  trackCellI
267  );
268 
269  bool bReached = trackToBoundary
270  (
271  singleParticle,
272  sampleI,
273  samplingPts,
274  samplingCells,
275  samplingFaces,
276  samplingCurveDist
277  );
278 
279  // fill sampleSegments
280  for(label i = samplingPts.size() - 1; i >= startSegmentI; --i)
281  {
282  samplingSegments.append(segmentI);
283  }
284 
285  if (!bReached)
286  {
287  //Info<< "calcSamples : Reached end of samples: "
288  // << " sampleI now:" << sampleI
289  // << endl;
290  break;
291  }
292  lastSample = singleParticle.position();
293 
294 
295  // Find next boundary.
296  sampleI++;
297 
298  if (sampleI == sampleCoords_.size() - 1)
299  {
300  //Info<< "calcSamples : Reached end of samples: "
301  // << " sampleI now:" << sampleI
302  // << endl;
303  break;
304  }
305 
306  segmentI++;
307 
308  startSegmentI = samplingPts.size();
309  }
310 }
311 
312 
313 void Foam::curveSet::genSamples()
314 {
315  // Storage for sample points
316  DynamicList<point> samplingPts;
317  DynamicList<label> samplingCells;
318  DynamicList<label> samplingFaces;
319  DynamicList<label> samplingSegments;
320  DynamicList<scalar> samplingCurveDist;
321 
322  calcSamples
323  (
324  samplingPts,
325  samplingCells,
326  samplingFaces,
327  samplingSegments,
328  samplingCurveDist
329  );
330 
331  samplingPts.shrink();
332  samplingCells.shrink();
333  samplingFaces.shrink();
334  samplingSegments.shrink();
335  samplingCurveDist.shrink();
336 
337  setSamples
338  (
339  samplingPts,
340  samplingCells,
341  samplingFaces,
342  samplingSegments,
343  samplingCurveDist
344  );
345 }
346 
347 
348 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
349 
351 (
352  const word& name,
353  const polyMesh& mesh,
354  meshSearch& searchEngine,
355  const word& axis,
356  const List<point>& sampleCoords
357 )
358 :
359  sampledSet(name, mesh, searchEngine, axis),
360  sampleCoords_(sampleCoords)
361 {
362  genSamples();
363 
364  if (debug)
365  {
366  write(Info);
367  }
368 }
369 
370 
372 (
373  const word& name,
374  const polyMesh& mesh,
375  meshSearch& searchEngine,
376  const dictionary& dict
377 )
378 :
379  sampledSet(name, mesh, searchEngine, dict),
380  sampleCoords_(dict.lookup("points"))
381 {
382  genSamples();
383 
384  if (debug)
385  {
386  write(Info);
387  }
388 }
389 
390 
391 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
392 
394 {}
395 
396 
397 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
398 
400 {
401  if (pts.size())
402  {
403  // Use first samplePt as starting point
404  return pts[0];
405  }
406  else
407  {
408  return vector::zero;
409  }
410 }
411 
412 
413 // ************************ vim: set sw=4 sts=4 et: ************************ //