surfaceToCell.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 "surfaceToCell.H"
27 #include "polyMesh.H"
28 #include "meshSearch.H"
29 #include "triSurface.H"
30 #include "triSurfaceSearch.H"
31 #include "cellClassification.H"
32 #include "cpuTime.H"
33 #include "demandDrivenData.H"
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
38 namespace Foam
39 {
40  defineTypeNameAndDebug(surfaceToCell, 0);
41  addToRunTimeSelectionTable(topoSetSource, surfaceToCell, word);
42  addToRunTimeSelectionTable(topoSetSource, surfaceToCell, istream);
43 }
44 
45 
46 Foam::topoSetSource::addToUsageTable Foam::surfaceToCell::usage_
47 (
48  surfaceToCell::typeName,
49  "\n Usage: surfaceToCell"
50  "<surface> <outsidePoints> <cut> <inside> <outside> <near> <curvature>\n\n"
51  " <surface> name of triSurface\n"
52  " <outsidePoints> list of points that define outside\n"
53  " <cut> boolean whether to include cells cut by surface\n"
54  " <inside> ,, ,, inside surface\n"
55  " <outside> ,, ,, outside surface\n"
56  " <near> scalar; include cells with centre <= near to surface\n"
57  " <curvature> scalar; include cells close to strong curvature"
58  " on surface\n"
59  " (curvature defined as difference in surface normal at nearest"
60  " point on surface for each vertex of cell)\n\n"
61 );
62 
63 
64 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
65 
66 Foam::label Foam::surfaceToCell::getNearest
67 (
68  const triSurfaceSearch& querySurf,
69  const label pointi,
70  const point& pt,
71  const vector& span,
72  Map<label>& cache
73 )
74 {
75  Map<label>::const_iterator iter = cache.find(pointi);
76 
77  if (iter != cache.end())
78  {
79  // Found cached answer
80  return iter();
81  }
82  else
83  {
84  pointIndexHit inter = querySurf.nearest(pt, span);
85 
86  // Triangle label (can be -1)
87  label triI = inter.index();
88 
89  // Store triangle on point
90  cache.insert(pointi, triI);
91 
92  return triI;
93  }
94 }
95 
96 
97 bool Foam::surfaceToCell::differingPointNormals
98 (
99  const triSurfaceSearch& querySurf,
100 
101  const vector& span, // Current search span
102  const label celli,
103  const label cellTriI, // Nearest (to cell centre) surface triangle
104 
105  Map<label>& pointToNearest // Cache for nearest triangle to point
106 ) const
107 {
108  const triSurface& surf = querySurf.surface();
109  const vectorField& normals = surf.faceNormals();
110 
111  const faceList& faces = mesh().faces();
112  const pointField& points = mesh().points();
113 
114  const labelList& cFaces = mesh().cells()[celli];
115 
116  forAll(cFaces, cFacei)
117  {
118  const face& f = faces[cFaces[cFacei]];
119 
120  forAll(f, fp)
121  {
122  label pointi = f[fp];
123 
124  label pointTriI =
125  getNearest
126  (
127  querySurf,
128  pointi,
129  points[pointi],
130  span,
131  pointToNearest
132  );
133 
134  if (pointTriI != -1 && pointTriI != cellTriI)
135  {
136  scalar cosAngle = normals[pointTriI] & normals[cellTriI];
137 
138  if (cosAngle < 0.9)
139  {
140  return true;
141  }
142  }
143  }
144  }
145  return false;
146 }
147 
148 
149 void Foam::surfaceToCell::combine(topoSet& set, const bool add) const
150 {
151  cpuTime timer;
152 
153 
154  if (useSurfaceOrientation_ && (includeInside_ || includeOutside_))
155  {
156  const meshSearch queryMesh(mesh_);
157 
158  //- Calculate for each searchPoint inside/outside status.
159  boolList isInside(querySurf().calcInside(mesh_.cellCentres()));
160 
161  Info<< " Marked inside/outside using surface orientation in = "
162  << timer.cpuTimeIncrement() << " s" << endl << endl;
163 
164  forAll(isInside, celli)
165  {
166  if (isInside[celli] && includeInside_)
167  {
168  addOrDelete(set, celli, add);
169  }
170  else if (!isInside[celli] && includeOutside_)
171  {
172  addOrDelete(set, celli, add);
173  }
174  }
175  }
176  else if (includeCut_ || includeInside_ || includeOutside_)
177  {
178  //
179  // Cut cells with surface and classify cells
180  //
181 
182 
183  // Construct search engine on mesh
184 
185  const meshSearch queryMesh(mesh_);
186 
187 
188  // Check all 'outside' points
189  forAll(outsidePoints_, outsideI)
190  {
191  const point& outsidePoint = outsidePoints_[outsideI];
192 
193  // Find cell point is in. Linear search.
194  label celli = queryMesh.findCell(outsidePoint, -1, false);
195  if (returnReduce(celli, maxOp<label>()) == -1)
196  {
198  << "outsidePoint " << outsidePoint
199  << " is not inside any cell"
200  << exit(FatalError);
201  }
202  }
203 
204  // Cut faces with surface and classify cells
205 
206  cellClassification cellType
207  (
208  mesh_,
209  queryMesh,
210  querySurf(),
211  outsidePoints_
212  );
213 
214 
215  Info<< " Marked inside/outside using surface intersection in = "
216  << timer.cpuTimeIncrement() << " s" << endl << endl;
217 
218  //- Add/remove cells using set
219  forAll(cellType, celli)
220  {
221  label cType = cellType[celli];
222 
223  if
224  (
225  (
226  includeCut_
227  && (cType == cellClassification::CUT)
228  )
229  || (
230  includeInside_
231  && (cType == cellClassification::INSIDE)
232  )
233  || (
234  includeOutside_
235  && (cType == cellClassification::OUTSIDE)
236  )
237  )
238  {
239  addOrDelete(set, celli, add);
240  }
241  }
242  }
243 
244 
245  if (nearDist_ > 0)
246  {
247  //
248  // Determine distance to surface
249  //
250 
251  const pointField& ctrs = mesh_.cellCentres();
252 
253  // Box dimensions to search in octree.
254  const vector span(nearDist_, nearDist_, nearDist_);
255 
256 
257  if (curvature_ < -1)
258  {
259  Info<< " Selecting cells with cellCentre closer than "
260  << nearDist_ << " to surface" << endl;
261 
262  // No need to test curvature. Insert near cells into set.
263 
264  forAll(ctrs, celli)
265  {
266  const point& c = ctrs[celli];
267 
268  pointIndexHit inter = querySurf().nearest(c, span);
269 
270  if (inter.hit() && (mag(inter.hitPoint() - c) < nearDist_))
271  {
272  addOrDelete(set, celli, add);
273  }
274  }
275 
276  Info<< " Determined nearest surface point in = "
277  << timer.cpuTimeIncrement() << " s" << endl << endl;
278 
279  }
280  else
281  {
282  // Test near cells for curvature
283 
284  Info<< " Selecting cells with cellCentre closer than "
285  << nearDist_ << " to surface and curvature factor"
286  << " less than " << curvature_ << endl;
287 
288  // Cache for nearest surface triangle for a point
289  Map<label> pointToNearest(mesh_.nCells()/10);
290 
291  forAll(ctrs, celli)
292  {
293  const point& c = ctrs[celli];
294 
295  pointIndexHit inter = querySurf().nearest(c, span);
296 
297  if (inter.hit() && (mag(inter.hitPoint() - c) < nearDist_))
298  {
299  if
300  (
301  differingPointNormals
302  (
303  querySurf(),
304  span,
305  celli,
306  inter.index(), // nearest surface triangle
307  pointToNearest
308  )
309  )
310  {
311  addOrDelete(set, celli, add);
312  }
313  }
314  }
315 
316  Info<< " Determined nearest surface point in = "
317  << timer.cpuTimeIncrement() << " s" << endl << endl;
318  }
319  }
320 }
321 
322 
323 void Foam::surfaceToCell::checkSettings() const
324 {
325  if
326  (
327  (nearDist_ < 0)
328  && (curvature_ < -1)
329  && (
330  (includeCut_ && includeInside_ && includeOutside_)
331  || (!includeCut_ && !includeInside_ && !includeOutside_)
332  )
333  )
334  {
336  << "Illegal include cell specification."
337  << " Result would be either all or no cells." << endl
338  << "Please set one of includeCut, includeInside, includeOutside"
339  << " to true, set nearDistance to a value > 0"
340  << " or set curvature to a value -1 .. 1."
341  << exit(FatalError);
342  }
343 
344  if (useSurfaceOrientation_ && includeCut_)
345  {
347  << "Illegal include cell specification."
348  << " You cannot specify both 'useSurfaceOrientation'"
349  << " and 'includeCut'"
350  << " since 'includeCut' specifies a topological split"
351  << exit(FatalError);
352  }
353 }
354 
355 
356 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
357 
359 (
360  const polyMesh& mesh,
361  const fileName& surfName,
362  const pointField& outsidePoints,
363  const bool includeCut,
364  const bool includeInside,
365  const bool includeOutside,
366  const bool useSurfaceOrientation,
367  const scalar nearDist,
368  const scalar curvature
369 )
370 :
371  topoSetSource(mesh),
372  surfName_(surfName),
373  outsidePoints_(outsidePoints),
374  includeCut_(includeCut),
375  includeInside_(includeInside),
376  includeOutside_(includeOutside),
377  useSurfaceOrientation_(useSurfaceOrientation),
378  nearDist_(nearDist),
379  curvature_(curvature),
380  surfPtr_(new triSurface(surfName_)),
381  querySurfPtr_(new triSurfaceSearch(*surfPtr_)),
382  IOwnPtrs_(true)
383 {
384  checkSettings();
385 }
386 
387 
389 (
390  const polyMesh& mesh,
391  const fileName& surfName,
392  const triSurface& surf,
393  const triSurfaceSearch& querySurf,
394  const pointField& outsidePoints,
395  const bool includeCut,
396  const bool includeInside,
397  const bool includeOutside,
398  const bool useSurfaceOrientation,
399  const scalar nearDist,
400  const scalar curvature
401 )
402 :
403  topoSetSource(mesh),
404  surfName_(surfName),
405  outsidePoints_(outsidePoints),
406  includeCut_(includeCut),
407  includeInside_(includeInside),
408  includeOutside_(includeOutside),
409  useSurfaceOrientation_(useSurfaceOrientation),
410  nearDist_(nearDist),
411  curvature_(curvature),
412  surfPtr_(&surf),
413  querySurfPtr_(&querySurf),
414  IOwnPtrs_(false)
415 {
416  checkSettings();
417 }
418 
419 
421 (
422  const polyMesh& mesh,
423  const dictionary& dict
424 )
425 :
426  topoSetSource(mesh),
427  surfName_(fileName(dict.lookup("file")).expand()),
428  outsidePoints_(dict.lookup("outsidePoints")),
429  includeCut_(readBool(dict.lookup("includeCut"))),
430  includeInside_(readBool(dict.lookup("includeInside"))),
431  includeOutside_(readBool(dict.lookup("includeOutside"))),
432  useSurfaceOrientation_
433  (
434  dict.lookupOrDefault<bool>("useSurfaceOrientation", false)
435  ),
436  nearDist_(readScalar(dict.lookup("nearDistance"))),
437  curvature_(readScalar(dict.lookup("curvature"))),
438  surfPtr_(new triSurface(surfName_)),
439  querySurfPtr_(new triSurfaceSearch(*surfPtr_)),
440  IOwnPtrs_(true)
441 {
442  checkSettings();
443 }
444 
445 
447 (
448  const polyMesh& mesh,
449  Istream& is
450 )
451 :
452  topoSetSource(mesh),
453  surfName_(checkIs(is)),
454  outsidePoints_(checkIs(is)),
455  includeCut_(readBool(checkIs(is))),
456  includeInside_(readBool(checkIs(is))),
457  includeOutside_(readBool(checkIs(is))),
458  useSurfaceOrientation_(false),
459  nearDist_(readScalar(checkIs(is))),
460  curvature_(readScalar(checkIs(is))),
461  surfPtr_(new triSurface(surfName_)),
462  querySurfPtr_(new triSurfaceSearch(*surfPtr_)),
463  IOwnPtrs_(true)
464 {
465  checkSettings();
466 }
467 
468 
469 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
470 
472 {
473  if (IOwnPtrs_)
474  {
475  deleteDemandDrivenData(surfPtr_);
476  deleteDemandDrivenData(querySurfPtr_);
477  }
478 }
479 
480 
481 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
482 
484 (
485  const topoSetSource::setAction action,
486  topoSet& set
487 ) const
488 {
489  if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
490  {
491  Info<< " Adding cells in relation to surface " << surfName_
492  << " ..." << endl;
493 
494  combine(set, true);
495  }
496  else if (action == topoSetSource::DELETE)
497  {
498  Info<< " Removing cells in relation to surface " << surfName_
499  << " ..." << endl;
500 
501  combine(set, false);
502  }
503 }
504 
505 
506 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
string expand(const string &, const HashTable< string, word, string::hash > &mapping, const char sigil='$')
Expand occurences of variables according to the mapping.
Definition: stringOps.C:75
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
A class for handling file names.
Definition: fileName.H:69
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
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
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
virtual void applyToSet(const topoSetSource::setAction action, topoSet &) const
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
HashTable< label, label, Hash< label > >::const_iterator const_iterator
Definition: Map.H:59
const cellList & cells() const
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
bool readBool(Istream &)
Definition: boolIO.C:60
AccessType combine(const List< T > &, AccessOp aop=accessOp< T >())
Combines sublists into one big list.
Definition: ListListOps.C:34
virtual ~surfaceToCell()
Destructor.
Macros for easy insertion into run-time selection tables.
Base class of a source for a topoSet.
Definition: topoSetSource.H:63
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1011
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
Helper class to search on triSurface.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
dynamicFvMesh & mesh
setAction
Enumeration defining the valid actions.
Definition: topoSetSource.H:82
List< label > labelList
A List of labels.
Definition: labelList.H:56
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1036
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
General set of labels of mesh quantity (points, cells, faces).
Definition: topoSet.H:61
surfaceToCell(const polyMesh &mesh, const fileName &surfName, const pointField &outsidePoints, const bool includeCut, const bool includeInside, const bool includeOutside, const bool useSurfaceOrientation, const scalar nearDist, const scalar curvature)
Construct from components.
Template functions to aid in the implementation of demand driven data.
vector point
Point is a vector.
Definition: point.H:41
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
Class with constructor to add usage string to table.
const dimensionedScalar c
Speed of light in a vacuum.
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Triangulated surface description with patch information.
Definition: triSurface.H:65
void deleteDemandDrivenData(DataPtr &dataPtr)
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:576