FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2011 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 "surfaceToCell.H"
27 #include <OpenFOAM/polyMesh.H>
28 #include <meshTools/meshSearch.H>
29 #include <triSurface/triSurface.H>
32 
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39 
40 defineTypeNameAndDebug(surfaceToCell, 0);
41 
42 addToRunTimeSelectionTable(topoSetSource, surfaceToCell, word);
43 
44 addToRunTimeSelectionTable(topoSetSource, surfaceToCell, istream);
45 
46 }
47 
48 
49 Foam::topoSetSource::addToUsageTable Foam::surfaceToCell::usage_
50 (
51  surfaceToCell::typeName,
52  "\n Usage: surfaceToCell"
53  "<surface> <outsidePoints> <cut> <inside> <outside> <near> <curvature>\n\n"
54  " <surface> name of triSurface\n"
55  " <outsidePoints> list of points that define outside\n"
56  " <cut> boolean whether to include cells cut by surface\n"
57  " <inside> ,, ,, inside surface\n"
58  " <outside> ,, ,, outside surface\n"
59  " <near> scalar; include cells with centre <= near to surface\n"
60  " <curvature> scalar; include cells close to strong curvature"
61  " on surface\n"
62  " (curvature defined as difference in surface normal at nearest"
63  " point on surface for each vertex of cell)\n\n"
64 );
65 
66 
67 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
68 
69 Foam::label Foam::surfaceToCell::getNearest
70 (
71  const triSurfaceSearch& querySurf,
72  const label pointI,
73  const point& pt,
74  const vector& span,
75  Map<label>& cache
76 )
77 {
78  Map<label>::const_iterator iter = cache.find(pointI);
79 
80  if (iter != cache.end())
81  {
82  // Found cached answer
83  return iter();
84  }
85  else
86  {
87  pointIndexHit inter = querySurf.nearest(pt, span);
88 
89  // Triangle label (can be -1)
90  label triI = inter.index();
91 
92  // Store triangle on point
93  cache.insert(pointI, triI);
94 
95  return triI;
96  }
97 }
98 
99 
100 // Return true if nearest surface to points on cell makes largish angle
101 // with nearest surface to cell centre. Returns false otherwise. Points visited
102 // are cached in pointToNearest
103 bool Foam::surfaceToCell::differingPointNormals
104 (
105  const triSurfaceSearch& querySurf,
106 
107  const vector& span, // current search span
108  const label cellI,
109  const label cellTriI, // nearest (to cell centre) surface triangle
110 
111  Map<label>& pointToNearest // cache for nearest triangle to point
112 ) const
113 {
114  const triSurface& surf = querySurf.surface();
115  const vectorField& normals = surf.faceNormals();
116 
117  const faceList& faces = mesh().faces();
118  const pointField& points = mesh().points();
119 
120  const labelList& cFaces = mesh().cells()[cellI];
121 
122  forAll(cFaces, cFaceI)
123  {
124  const face& f = faces[cFaces[cFaceI]];
125 
126  forAll(f, fp)
127  {
128  label pointI = f[fp];
129 
130  label pointTriI =
131  getNearest
132  (
133  querySurf,
134  pointI,
135  points[pointI],
136  span,
137  pointToNearest
138  );
139 
140  if (pointTriI != -1 && pointTriI != cellTriI)
141  {
142  scalar cosAngle = normals[pointTriI] & normals[cellTriI];
143 
144  if (cosAngle < 0.9)
145  {
146  return true;
147  }
148  }
149  }
150  }
151  return false;
152 }
153 
154 
155 void Foam::surfaceToCell::combine(topoSet& set, const bool add) const
156 {
157  cpuTime timer;
158 
159  if (includeCut_ || includeInside_ || includeOutside_)
160  {
161  //
162  // Cut cells with surface and classify cells
163  //
164 
165 
166  // Construct search engine on mesh
167 
168  meshSearch queryMesh(mesh_, true);
169 
170 
171  // Check all 'outside' points
172  forAll(outsidePoints_, outsideI)
173  {
174  const point& outsidePoint = outsidePoints_[outsideI];
175 
176  // Find cell point is in. Linear search.
177  label cellI = queryMesh.findCell(outsidePoint, -1, false);
178  if (returnReduce(cellI, maxOp<label>()) == -1)
179  {
180  FatalErrorIn("surfaceToCell::combine(topoSet&, const bool)")
181  << "outsidePoint " << outsidePoint
182  << " is not inside any cell"
183  << exit(FatalError);
184  }
185  }
186 
187  // Cut faces with surface and classify cells
188 
189  cellClassification cellType
190  (
191  mesh_,
192  queryMesh,
193  querySurf(),
194  outsidePoints_
195  );
196 
197 
198  Info<< " Marked inside/outside in = "
199  << timer.cpuTimeIncrement() << " s" << endl << endl;
200 
201 
202  forAll(cellType, cellI)
203  {
204  label cType = cellType[cellI];
205 
206  if
207  (
208  (
209  includeCut_
210  && (cType == cellClassification::CUT)
211  )
212  || (
213  includeInside_
214  && (cType == cellClassification::INSIDE)
215  )
216  || (
217  includeOutside_
218  && (cType == cellClassification::OUTSIDE)
219  )
220  )
221  {
222  addOrDelete(set, cellI, add);
223  }
224  }
225  }
226 
227 
228  if (nearDist_ > 0)
229  {
230  //
231  // Determine distance to surface
232  //
233 
234  const pointField& ctrs = mesh_.cellCentres();
235 
236  // Box dimensions to search in octree.
237  const vector span(nearDist_, nearDist_, nearDist_);
238 
239 
240  if (curvature_ < -1)
241  {
242  Info<< " Selecting cells with cellCentre closer than "
243  << nearDist_ << " to surface" << endl;
244 
245  // No need to test curvature. Insert near cells into set.
246 
247  forAll(ctrs, cellI)
248  {
249  const point& c = ctrs[cellI];
250 
251  pointIndexHit inter = querySurf().nearest(c, span);
252 
253  if (inter.hit() && (mag(inter.hitPoint() - c) < nearDist_))
254  {
255  addOrDelete(set, cellI, add);
256  }
257  }
258 
259  Info<< " Determined nearest surface point in = "
260  << timer.cpuTimeIncrement() << " s" << endl << endl;
261 
262  }
263  else
264  {
265  // Test near cells for curvature
266 
267  Info<< " Selecting cells with cellCentre closer than "
268  << nearDist_ << " to surface and curvature factor"
269  << " less than " << curvature_ << endl;
270 
271  // Cache for nearest surface triangle for a point
272  Map<label> pointToNearest(mesh_.nCells()/10);
273 
274  forAll(ctrs, cellI)
275  {
276  const point& c = ctrs[cellI];
277 
278  pointIndexHit inter = querySurf().nearest(c, span);
279 
280  if (inter.hit() && (mag(inter.hitPoint() - c) < nearDist_))
281  {
282  if
283  (
284  differingPointNormals
285  (
286  querySurf(),
287  span,
288  cellI,
289  inter.index(), // nearest surface triangle
290  pointToNearest
291  )
292  )
293  {
294  addOrDelete(set, cellI, add);
295  }
296  }
297  }
298 
299  Info<< " Determined nearest surface point in = "
300  << timer.cpuTimeIncrement() << " s" << endl << endl;
301  }
302  }
303 }
304 
305 
306 void Foam::surfaceToCell::checkSettings() const
307 {
308  if
309  (
310  (nearDist_ < 0)
311  && (curvature_ < -1)
312  && (
313  (includeCut_ && includeInside_ && includeOutside_)
314  || (!includeCut_ && !includeInside_ && !includeOutside_)
315  )
316  )
317  {
319  (
320  "surfaceToCell:checkSettings()"
321  ) << "Illegal include cell specification."
322  << " Result would be either all or no cells." << endl
323  << "Please set one of includeCut, includeInside, includeOutside"
324  << " to true, set nearDistance to a value > 0"
325  << " or set curvature to a value -1 .. 1."
326  << exit(FatalError);
327  }
328 }
329 
330 
331 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
332 
333 // Construct from components
335 (
336  const polyMesh& mesh,
337  const fileName& surfName,
338  const pointField& outsidePoints,
339  const bool includeCut,
340  const bool includeInside,
341  const bool includeOutside,
342  const scalar nearDist,
343  const scalar curvature
344 )
345 :
346  topoSetSource(mesh),
347  surfName_(surfName),
348  outsidePoints_(outsidePoints),
349  includeCut_(includeCut),
350  includeInside_(includeInside),
351  includeOutside_(includeOutside),
352  nearDist_(nearDist),
353  curvature_(curvature),
354  surfPtr_(new triSurface(surfName_)),
355  querySurfPtr_(new triSurfaceSearch(*surfPtr_)),
356  IOwnPtrs_(true)
357 {
358  checkSettings();
359 }
360 
361 
362 // Construct from components. Externally supplied surface.
364 (
365  const polyMesh& mesh,
366  const fileName& surfName,
367  const triSurface& surf,
368  const triSurfaceSearch& querySurf,
369  const pointField& outsidePoints,
370  const bool includeCut,
371  const bool includeInside,
372  const bool includeOutside,
373  const scalar nearDist,
374  const scalar curvature
375 )
376 :
377  topoSetSource(mesh),
378  surfName_(surfName),
379  outsidePoints_(outsidePoints),
380  includeCut_(includeCut),
381  includeInside_(includeInside),
382  includeOutside_(includeOutside),
383  nearDist_(nearDist),
384  curvature_(curvature),
385  surfPtr_(&surf),
386  querySurfPtr_(&querySurf),
387  IOwnPtrs_(false)
388 {
389  checkSettings();
390 }
391 
392 
393 // Construct from dictionary
395 (
396  const polyMesh& mesh,
397  const dictionary& dict
398 )
399 :
400  topoSetSource(mesh),
401  surfName_(dict.lookup("file")),
402  outsidePoints_(dict.lookup("outsidePoints")),
403  includeCut_(readBool(dict.lookup("includeCut"))),
404  includeInside_(readBool(dict.lookup("includeInside"))),
405  includeOutside_(readBool(dict.lookup("includeOutside"))),
406  nearDist_(readScalar(dict.lookup("nearDistance"))),
407  curvature_(readScalar(dict.lookup("curvature"))),
408  surfPtr_(new triSurface(surfName_)),
409  querySurfPtr_(new triSurfaceSearch(*surfPtr_)),
410  IOwnPtrs_(true)
411 {
412  checkSettings();
413 }
414 
415 
416 // Construct from Istream
418 (
419  const polyMesh& mesh,
420  Istream& is
421 )
422 :
423  topoSetSource(mesh),
424  surfName_(checkIs(is)),
425  outsidePoints_(checkIs(is)),
426  includeCut_(readBool(checkIs(is))),
427  includeInside_(readBool(checkIs(is))),
428  includeOutside_(readBool(checkIs(is))),
429  nearDist_(readScalar(checkIs(is))),
430  curvature_(readScalar(checkIs(is))),
431  surfPtr_(new triSurface(surfName_)),
432  querySurfPtr_(new triSurfaceSearch(*surfPtr_)),
433  IOwnPtrs_(true)
434 {
435  checkSettings();
436 }
437 
438 
439 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
440 
442 {
443  if (IOwnPtrs_)
444  {
445  if (surfPtr_)
446  {
447  delete surfPtr_;
448  }
449  if (querySurfPtr_)
450  {
451  delete querySurfPtr_;
452  }
453  }
454 }
455 
456 
457 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
458 
460 (
461  const topoSetSource::setAction action,
462  topoSet& set
463 ) const
464 {
465  if ( (action == topoSetSource::NEW) || (action == topoSetSource::ADD))
466  {
467  Info<< " Adding cells in relation to surface " << surfName_
468  << " ..." << endl;
469 
470  combine(set, true);
471  }
472  else if (action == topoSetSource::DELETE)
473  {
474  Info<< " Removing cells in relation to surface " << surfName_
475  << " ..." << endl;
476 
477  combine(set, false);
478  }
479 }
480 
481 
482 // ************************ vim: set sw=4 sts=4 et: ************************ //