FreeFOAM The Cross-Platform CFD Toolkit
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) 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 "sampledSet.H"
27 #include <OpenFOAM/polyMesh.H>
28 #include <OpenFOAM/primitiveMesh.H>
29 #include <meshTools/meshSearch.H>
30 #include <sampling/writer.H>
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36  const scalar sampledSet::tol = 1e-6;
37 
38  defineTypeNameAndDebug(sampledSet, 0);
39  defineRunTimeSelectionTable(sampledSet, word);
40 }
41 
42 
43 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44 
45 Foam::label Foam::sampledSet::getBoundaryCell(const label faceI) const
46 {
47  return mesh().faceOwner()[faceI];
48 }
49 
50 
51 Foam::label Foam::sampledSet::getCell
52 (
53  const label faceI,
54  const point& sample
55 ) const
56 {
57  if (faceI == -1)
58  {
60  (
61  "sampledSet::getCell(const label, const point&)"
62  ) << "Illegal face label " << faceI
63  << abort(FatalError);
64  }
65 
66  if (faceI >= mesh().nInternalFaces())
67  {
68  label cellI = getBoundaryCell(faceI);
69 
70  if (!mesh().pointInCell(sample, cellI))
71  {
73  (
74  "sampledSet::getCell(const label, const point&)"
75  ) << "Found cell " << cellI << " using face " << faceI
76  << ". But cell does not contain point " << sample
77  << abort(FatalError);
78  }
79  return cellI;
80  }
81  else
82  {
83  // Try owner and neighbour to see which one contains sample
84 
85  label cellI = mesh().faceOwner()[faceI];
86 
87  if (mesh().pointInCell(sample, cellI))
88  {
89  return cellI;
90  }
91  else
92  {
93  cellI = mesh().faceNeighbour()[faceI];
94 
95  if (mesh().pointInCell(sample, cellI))
96  {
97  return cellI;
98  }
99  else
100  {
102  (
103  "sampledSet::getCell(const label, const point&)"
104  ) << "None of the neighbours of face "
105  << faceI << " contains point " << sample
106  << abort(FatalError);
107 
108  return -1;
109  }
110  }
111  }
112 }
113 
114 
115 Foam::scalar Foam::sampledSet::calcSign
116 (
117  const label faceI,
118  const point& sample
119 ) const
120 {
121  vector vec = sample - mesh().faceCentres()[faceI];
122 
123  scalar magVec = mag(vec);
124 
125  if (magVec < VSMALL)
126  {
127  // sample on face centre. Regard as inside
128  return -1;
129  }
130 
131  vec /= magVec;
132 
133  vector n = mesh().faceAreas()[faceI];
134 
135  n /= mag(n) + VSMALL;
136 
137  return n & vec;
138 }
139 
140 
141 // Return face (or -1) of face which is within smallDist of sample
143 (
144  const label cellI,
145  const point& sample,
146  const scalar smallDist
147 ) const
148 {
149  const cell& myFaces = mesh().cells()[cellI];
150 
151  forAll(myFaces, myFaceI)
152  {
153  const face& f = mesh().faces()[myFaces[myFaceI]];
154 
155  pointHit inter = f.nearestPoint(sample, mesh().points());
156 
157  scalar dist;
158 
159  if (inter.hit())
160  {
161  dist = mag(inter.hitPoint() - sample);
162  }
163  else
164  {
165  dist = mag(inter.missPoint() - sample);
166  }
167 
168  if (dist < smallDist)
169  {
170  return myFaces[myFaceI];
171  }
172  }
173  return -1;
174 }
175 
176 
177 // 'Pushes' point facePt (which is almost on face) in direction of cell centre
178 // so it is clearly inside.
180 (
181  const point& facePt,
182  const label faceI
183 ) const
184 {
185  label cellI = mesh().faceOwner()[faceI];
186 
187  const point& cellCtr = mesh().cellCentres()[cellI];
188 
189  point newSample =
190  facePt + tol*(cellCtr - facePt);
191 
192  if (!searchEngine().pointInCell(newSample, cellI))
193  {
195  (
196  "sampledSet::pushIn(const point&, const label)"
197  ) << "After pushing " << facePt << " to " << newSample
198  << " it is still outside faceI " << faceI << endl
199  << "Please change your starting point"
200  << abort(FatalError);
201  }
202  //Info<< "pushIn : moved " << facePt << " to " << newSample
203  // << endl;
204 
205  return newSample;
206 }
207 
208 
209 // Calculates start of tracking given samplePt and first boundary intersection
210 // (bPoint, bFaceI). bFaceI == -1 if no boundary intersection.
211 // Returns true if trackPt is sampling point
213 (
214  const vector& offset,
215  const point& samplePt,
216  const point& bPoint,
217  const label bFaceI,
218 
219  point& trackPt,
220  label& trackCellI,
221  label& trackFaceI
222 ) const
223 {
224  const scalar smallDist = mag(tol*offset);
225 
226  bool isGoodSample = false;
227 
228  if (bFaceI == -1)
229  {
230  // No boundary intersection. Try and find cell samplePt is in
231  trackCellI = mesh().findCell(samplePt);
232 
233  if
234  (
235  (trackCellI == -1)
236  || !mesh().pointInCell(samplePt, trackCellI)
237  )
238  {
239  // Line samplePt - end_ does not intersect domain at all.
240  // (or is along edge)
241  //Info<< "getTrackingPoint : samplePt outside domain : "
242  // << " samplePt:" << samplePt
243  // << endl;
244 
245  trackCellI = -1;
246  trackFaceI = -1;
247 
248  isGoodSample = false;
249  }
250  else
251  {
252  // start is inside. Use it as tracking point
253  //Info<< "getTrackingPoint : samplePt inside :"
254  // << " samplePt:" << samplePt
255  // << " trackCellI:" << trackCellI
256  // << endl;
257 
258  trackPt = samplePt;
259  trackFaceI = -1;
260 
261  isGoodSample = true;
262  }
263  }
264  else if (mag(samplePt - bPoint) < smallDist)
265  {
266  //Info<< "getTrackingPoint : samplePt:" << samplePt
267  // << " close to bPoint:"
268  // << bPoint << endl;
269 
270  // samplePt close to bPoint. Snap to it
271  trackPt = pushIn(bPoint, bFaceI);
272  trackFaceI = bFaceI;
273  trackCellI = getBoundaryCell(trackFaceI);
274 
275  isGoodSample = true;
276  }
277  else
278  {
279  scalar sign = calcSign(bFaceI, samplePt);
280 
281  if (sign < 0)
282  {
283  // samplePt inside or marginally outside.
284  trackPt = samplePt;
285  trackFaceI = -1;
286  trackCellI = mesh().findCell(trackPt);
287 
288  isGoodSample = true;
289  }
290  else
291  {
292  // samplePt outside. use bPoint
293  trackPt = pushIn(bPoint, bFaceI);
294  trackFaceI = bFaceI;
295  trackCellI = getBoundaryCell(trackFaceI);
296 
297  isGoodSample = false;
298  }
299  }
300 
301  if (debug)
302  {
303  Info<< "sampledSet::getTrackingPoint :"
304  << " offset:" << offset
305  << " samplePt:" << samplePt
306  << " bPoint:" << bPoint
307  << " bFaceI:" << bFaceI
308  << endl << " Calculated first tracking point :"
309  << " trackPt:" << trackPt
310  << " trackCellI:" << trackCellI
311  << " trackFaceI:" << trackFaceI
312  << " isGoodSample:" << isGoodSample
313  << endl;
314  }
315 
316  return isGoodSample;
317 }
318 
319 
321 (
322  const List<point>& samplingPts,
323  const labelList& samplingCells,
324  const labelList& samplingFaces,
325  const labelList& samplingSegments,
326  const scalarList& samplingCurveDist
327 )
328 {
329  setSize(samplingPts.size());
330  cells_.setSize(samplingCells.size());
331  faces_.setSize(samplingFaces.size());
332  segments_.setSize(samplingSegments.size());
333  curveDist_.setSize(samplingCurveDist.size());
334 
335  if
336  (
337  (cells_.size() != size())
338  || (faces_.size() != size())
339  || (segments_.size() != size())
340  || (curveDist_.size() != size())
341  )
342  {
343  FatalErrorIn("sampledSet::setSamples()")
344  << "sizes not equal : "
345  << " points:" << size()
346  << " cells:" << cells_.size()
347  << " faces:" << faces_.size()
348  << " segments:" << segments_.size()
349  << " curveDist:" << curveDist_.size()
350  << abort(FatalError);
351  }
352 
353  forAll(samplingPts, sampleI)
354  {
355  operator[](sampleI) = samplingPts[sampleI];
356  }
357  cells_ = samplingCells;
358  faces_ = samplingFaces;
359  segments_ = samplingSegments;
360  curveDist_ = samplingCurveDist;
361 }
362 
363 
364 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
365 
367 (
368  const word& name,
369  const polyMesh& mesh,
370  meshSearch& searchEngine,
371  const word& axis
372 )
373 :
374  coordSet(name, axis),
375  mesh_(mesh),
376  searchEngine_(searchEngine),
377  segments_(0),
378  curveDist_(0),
379  cells_(0),
380  faces_(0)
381 {}
382 
383 
385 (
386  const word& name,
387  const polyMesh& mesh,
388  meshSearch& searchEngine,
389  const dictionary& dict
390 )
391 :
392  coordSet(name, dict.lookup("axis")),
393  mesh_(mesh),
394  searchEngine_(searchEngine),
395  segments_(0),
396  curveDist_(0),
397  cells_(0),
398  faces_(0)
399 {}
400 
401 
402 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
403 
405 {}
406 
407 
408 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
409 
411 (
412  const word& name,
413  const polyMesh& mesh,
414  meshSearch& searchEngine,
415  const dictionary& dict
416 )
417 {
418  word sampleType(dict.lookup("type"));
419 
420  wordConstructorTable::iterator cstrIter =
421  wordConstructorTablePtr_->find(sampleType);
422 
423  if (cstrIter == wordConstructorTablePtr_->end())
424  {
426  (
427  "sampledSet::New(const word&, "
428  "const polyMesh&, meshSearch&, const dictionary&)"
429  ) << "Unknown sample type " << sampleType
430  << nl << nl
431  << "Valid sample types : " << nl
432  << wordConstructorTablePtr_->sortedToc()
433  << exit(FatalError);
434  }
435 
436  return autoPtr<sampledSet>
437  (
438  cstrIter()
439  (
440  name,
441  mesh,
442  searchEngine,
443  dict
444  )
445  );
446 }
447 
448 
450 {
451  coordSet::write(os);
452 
453  os << endl << "\t(cellI)\t(faceI)" << endl;
454 
455  forAll(*this, sampleI)
456  {
457  os << '\t' << cells_[sampleI]
458  << '\t' << faces_[sampleI]
459  << endl;
460  }
461 
462  return os;
463 }
464 
465 
466 // ************************ vim: set sw=4 sts=4 et: ************************ //