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