sampledSet.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 "sampledSet.H"
27 #include "polyMesh.H"
28 #include "primitiveMesh.H"
29 #include "meshSearch.H"
30 #include "writer.H"
31 #include "particle.H"
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37  defineTypeNameAndDebug(sampledSet, 0);
38  defineRunTimeSelectionTable(sampledSet, word);
39 }
40 
41 
42 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43 
45 {
46  return mesh().faceOwner()[facei];
47 }
48 
49 
51 {
52  if (facei >= mesh().nInternalFaces())
53  {
54  return mesh().faceOwner()[facei];
55  }
56  else
57  {
58  return mesh().faceNeighbour()[facei];
59  }
60 }
61 
62 
64 (
65  const point& p,
66  const label samplei
67 ) const
68 {
69  // Collect the face owner and neighbour cells of the sample into an array
70  // for convenience
71  label cells[4] =
72  {
73  mesh().faceOwner()[faces_[samplei]],
74  getNeighbourCell(faces_[samplei]),
75  mesh().faceOwner()[faces_[samplei+1]],
76  getNeighbourCell(faces_[samplei+1])
77  };
78 
79  // Find the sampled cell by checking the owners and neighbours of the
80  // sampled faces
81  label cellm =
82  (cells[0] == cells[2] || cells[0] == cells[3]) ? cells[0]
83  : (cells[1] == cells[2] || cells[1] == cells[3]) ? cells[1]
84  : -1;
85 
86  if (cellm != -1)
87  {
88  // If found the sampled cell check the point is in the cell
89  // otherwise ignore
90  if (!mesh().pointInCell(p, cellm, searchEngine_.decompMode()))
91  {
92  cellm = -1;
93 
94  if (debug)
95  {
97  << "Could not find mid-point " << p
98  << " cell " << cellm << endl;
99  }
100  }
101  }
102  else
103  {
104  // If the sample does not pass through a single cell check if the point
105  // is in any of the owners or neighbours otherwise ignore
106  for (label i=0; i<4; i++)
107  {
108  if (mesh().pointInCell(p, cells[i], searchEngine_.decompMode()))
109  {
110  return cells[i];
111  }
112  }
113 
114  if (debug)
115  {
117  << "Could not find cell for mid-point" << nl
118  << " samplei: " << samplei
119  << " pts[samplei]: " << operator[](samplei)
120  << " face[samplei]: " << faces_[samplei]
121  << " pts[samplei+1]: " << operator[](samplei+1)
122  << " face[samplei+1]: " << faces_[samplei+1]
123  << " cellio: " << cells[0]
124  << " cellin: " << cells[1]
125  << " celljo: " << cells[2]
126  << " celljn: " << cells[3]
127  << endl;
128  }
129  }
130 
131  return cellm;
132 }
133 
134 
135 Foam::scalar Foam::sampledSet::calcSign
136 (
137  const label facei,
138  const point& sample
139 ) const
140 {
141  vector vec = sample - mesh().faceCentres()[facei];
142 
143  scalar magVec = mag(vec);
144 
145  if (magVec < VSMALL)
146  {
147  // Sample on face centre. Regard as inside
148  return -1;
149  }
150 
151  vec /= magVec;
152 
153  vector n = mesh().faceAreas()[facei];
154 
155  n /= mag(n) + VSMALL;
156 
157  return n & vec;
158 }
159 
160 
162 (
163  const label celli,
164  const point& sample,
165  const scalar smallDist
166 ) const
167 {
168  const cell& myFaces = mesh().cells()[celli];
169 
170  forAll(myFaces, myFacei)
171  {
172  const face& f = mesh().faces()[myFaces[myFacei]];
173 
174  pointHit inter = f.nearestPoint(sample, mesh().points());
175 
176  scalar dist;
177 
178  if (inter.hit())
179  {
180  dist = mag(inter.hitPoint() - sample);
181  }
182  else
183  {
184  dist = mag(inter.missPoint() - sample);
185  }
186 
187  if (dist < smallDist)
188  {
189  return myFaces[myFacei];
190  }
191  }
192  return -1;
193 }
194 
195 
197 (
198  const point& facePt,
199  const label facei
200 ) const
201 {
202  label celli = mesh().faceOwner()[facei];
203  const point& cC = mesh().cellCentres()[celli];
204 
205  point newPosition = facePt;
206 
207  // Taken from particle::initCellFacePt()
208  label tetFacei;
209  label tetPtI;
210  mesh().findTetFacePt(celli, facePt, tetFacei, tetPtI);
211 
212  if (tetFacei == -1 || tetPtI == -1)
213  {
214  newPosition = facePt;
215 
217 
218  label iterNo = 0;
219 
220  do
221  {
222  newPosition += particle::trackingCorrectionTol*(cC - facePt);
223 
225  (
226  celli,
227  newPosition,
228  tetFacei,
229  tetPtI
230  );
231 
232  iterNo++;
233 
234  } while (tetFacei < 0 && iterNo <= trap);
235  }
236 
237  if (tetFacei == -1)
238  {
240  << "After pushing " << facePt << " to " << newPosition
241  << " it is still outside face " << facei
242  << " at " << mesh().faceCentres()[facei]
243  << " of cell " << celli
244  << " at " << cC << endl
245  << "Please change your starting point"
246  << abort(FatalError);
247  }
248 
249  return newPosition;
250 }
251 
252 
254 (
255  const point& samplePt,
256  const point& bPoint,
257  const label bFacei,
258  const scalar smallDist,
259 
260  point& trackPt,
261  label& trackCelli,
262  label& trackFacei
263 ) const
264 {
265  bool isGoodSample = false;
266 
267  if (bFacei == -1)
268  {
269  // No boundary intersection. Try and find cell samplePt is in
270  trackCelli = mesh().findCell(samplePt, searchEngine_.decompMode());
271 
272  if
273  (
274  (trackCelli == -1)
275  || !mesh().pointInCell
276  (
277  samplePt,
278  trackCelli,
279  searchEngine_.decompMode()
280  )
281  )
282  {
283  // Line samplePt - end_ does not intersect domain at all.
284  // (or is along edge)
285 
286  trackCelli = -1;
287  trackFacei = -1;
288 
289  isGoodSample = false;
290  }
291  else
292  {
293  // Start is inside. Use it as tracking point
294 
295  trackPt = samplePt;
296  trackFacei = -1;
297 
298  isGoodSample = true;
299  }
300  }
301  else if (mag(samplePt - bPoint) < smallDist)
302  {
303  // samplePt close to bPoint. Snap to it
304  trackPt = pushIn(bPoint, bFacei);
305  trackFacei = bFacei;
306  trackCelli = getBoundaryCell(trackFacei);
307 
308  isGoodSample = true;
309  }
310  else
311  {
312  scalar sign = calcSign(bFacei, samplePt);
313 
314  if (sign < 0)
315  {
316  // samplePt inside or marginally outside.
317  trackPt = samplePt;
318  trackFacei = -1;
319  trackCelli = mesh().findCell(trackPt, searchEngine_.decompMode());
320 
321  isGoodSample = true;
322  }
323  else
324  {
325  // samplePt outside. use bPoint
326  trackPt = pushIn(bPoint, bFacei);
327  trackFacei = bFacei;
328  trackCelli = getBoundaryCell(trackFacei);
329 
330  isGoodSample = false;
331  }
332  }
333 
334  if (debug)
335  {
337  << " samplePt:" << samplePt
338  << " bPoint:" << bPoint
339  << " bFacei:" << bFacei
340  << endl << " Calculated first tracking point :"
341  << " trackPt:" << trackPt
342  << " trackCelli:" << trackCelli
343  << " trackFacei:" << trackFacei
344  << " isGoodSample:" << isGoodSample
345  << endl;
346  }
347 
348  return isGoodSample;
349 }
350 
351 
353 (
354  const List<point>& samplingPts,
355  const labelList& samplingCells,
356  const labelList& samplingFaces,
357  const labelList& samplingSegments,
358  const scalarList& samplingCurveDist
359 )
360 {
361  setSize(samplingPts.size());
362  cells_.setSize(samplingCells.size());
363  faces_.setSize(samplingFaces.size());
364  segments_.setSize(samplingSegments.size());
365  curveDist_.setSize(samplingCurveDist.size());
366 
367  if
368  (
369  (cells_.size() != size())
370  || (faces_.size() != size())
371  || (segments_.size() != size())
372  || (curveDist_.size() != size())
373  )
374  {
376  << "sizes not equal : "
377  << " points:" << size()
378  << " cells:" << cells_.size()
379  << " faces:" << faces_.size()
380  << " segments:" << segments_.size()
381  << " curveDist:" << curveDist_.size()
382  << abort(FatalError);
383  }
384 
385  forAll(samplingPts, sampleI)
386  {
387  operator[](sampleI) = samplingPts[sampleI];
388  }
389  curveDist_ = samplingCurveDist;
390 
391  cells_ = samplingCells;
392  faces_ = samplingFaces;
393  segments_ = samplingSegments;
394 }
395 
396 
397 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
398 
400 (
401  const word& name,
402  const polyMesh& mesh,
403  const meshSearch& searchEngine,
404  const word& axis
405 )
406 :
407  coordSet(name, axis),
408  mesh_(mesh),
409  searchEngine_(searchEngine),
410  segments_(0),
411  cells_(0),
412  faces_(0)
413 {}
414 
415 
417 (
418  const word& name,
419  const polyMesh& mesh,
420  const meshSearch& searchEngine,
421  const dictionary& dict
422 )
423 :
424  coordSet(name, dict.lookup("axis")),
425  mesh_(mesh),
426  searchEngine_(searchEngine),
427  segments_(0),
428  cells_(0),
429  faces_(0)
430 {}
431 
432 
433 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
434 
436 {}
437 
438 
439 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
440 
442 (
443  const word& name,
444  const polyMesh& mesh,
445  const meshSearch& searchEngine,
446  const dictionary& dict
447 )
448 {
449  const word sampleType(dict.lookup("type"));
450 
451  wordConstructorTable::iterator cstrIter =
452  wordConstructorTablePtr_->find(sampleType);
453 
454  if (cstrIter == wordConstructorTablePtr_->end())
455  {
457  << "Unknown sample type "
458  << sampleType << nl << nl
459  << "Valid sample types : " << endl
460  << wordConstructorTablePtr_->sortedToc()
461  << exit(FatalError);
462  }
463 
464  return autoPtr<sampledSet>
465  (
466  cstrIter()
467  (
468  name,
469  mesh,
470  searchEngine,
471  dict
472  )
473  );
474 }
475 
476 
478 {
479  coordSet::write(os);
480 
481  os << endl << "\t(celli)\t(facei)" << endl;
482 
483  forAll(*this, sampleI)
484  {
485  os << '\t' << cells_[sampleI]
486  << '\t' << faces_[sampleI]
487  << endl;
488  }
489 
490  return os;
491 }
492 
493 
494 // ************************************************************************* //
label getBoundaryCell(const label) const
Returns cell next to boundary face.
Definition: sampledSet.C:44
dimensionedScalar sign(const dimensionedScalar &ds)
label findNearFace(const label celli, const point &sample, const scalar smallDist) const
Returns face label (or -1) of face which is close to sample.
Definition: sampledSet.C:162
Various (local, not parallel) searches on polyMesh; uses (demand driven) octree to search...
Definition: meshSearch.H:57
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
static const scalar trackingCorrectionTol
Fraction of distance to tet centre to move a particle to.
Definition: particle.H:326
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
const Point & missPoint() const
Return miss point.
Definition: PointHit.H:145
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
static autoPtr< sampledSet > New(const word &name, const polyMesh &mesh, const meshSearch &searchEngine, const dictionary &dict)
Return a reference to the selected sampledSet.
Definition: sampledSet.C:442
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
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
const vectorField & faceAreas() const
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
virtual ~sampledSet()
Destructor.
Definition: sampledSet.C:435
Ostream & write(Ostream &) const
Output for debugging.
Definition: sampledSet.C:477
const vectorField & faceCentres() const
const cellList & cells() const
void setSamples(const List< point > &samplingPts, const labelList &samplingCells, const labelList &samplingFaces, const labelList &samplingSegments, const scalarList &samplingCurveDist)
Sets sample data.
Definition: sampledSet.C:353
bool getTrackingPoint(const point &samplePt, const point &bPoint, const label bFacei, const scalar smallDist, point &trackPt, label &trackCelli, label &trackFacei) const
Calculates start of tracking given samplePt and first boundary.
Definition: sampledSet.C:254
label pointInCell(const point &p, const label samplei) const
Return the cell in which the point on the sample line.
Definition: sampledSet.C:64
points setSize(newPointi)
dynamicFvMesh & mesh
Holds list of sampling positions.
Definition: coordSet.H:49
const cellShapeList & cells
const pointField & points
Ostream & write(Ostream &os) const
Definition: coordSet.C:134
A class for handling words, derived from string.
Definition: word.H:59
point pushIn(const point &sample, const label facei) const
Moves sample in direction of -n to it is &#39;inside&#39; of facei.
Definition: sampledSet.C:197
const vectorField & cellCentres() const
void findTetFacePt(const label celli, const point &p, label &tetFacei, label &tetPti) const
Find the tetFacei and tetPti for point p in celli.
Definition: polyMesh.C:1273
errorManip< error > abort(error &err)
Definition: errorManip.H:131
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
static const char nl
Definition: Ostream.H:262
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
defineTypeNameAndDebug(combustionModel, 0)
labelList f(nPoints)
pointHit nearestPoint(const point &p, const pointField &) const
Return nearest point to face.
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
sampledSet(const word &name, const polyMesh &mesh, const meshSearch &searchEngine, const word &axis)
Construct from components.
Definition: sampledSet.C:400
bool pointInCell(const point &p, label celli, const cellDecomposition=CELL_TETS) const
Test if point p is in the celli.
Definition: polyMesh.C:1289
scalar calcSign(const label facei, const point &sample) const
Calculates inproduct of face normal and vector sample-face centre.
Definition: sampledSet.C:136
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1023
#define WarningInFunction
Report a warning using Foam::Warning.
const Point & hitPoint() const
Return hit point.
Definition: PointHit.H:126
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
label getNeighbourCell(const label) const
Returns the neigbour cell or the owner if face in on the boundary.
Definition: sampledSet.C:50
label findCell(const point &p, const cellDecomposition=CELL_TETS) const
Find cell enclosing this location and return index.
Definition: polyMesh.C:1408
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
bool hit() const
Is there a hit.
Definition: PointHit.H:120
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1017
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1004
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
#define InfoInFunction
Report an information message using Foam::Info.