FreeFOAM The Cross-Platform CFD Toolkit
selectCells.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 Application
25  selectCells
26 
27 Description
28  Select cells in relation to surface.
29 
30  Divides cells into three sets:
31  - cutCells : cells cut by surface or close to surface.
32  - outside : cells not in cutCells and reachable from set of
33  user-defined points (outsidePoints)
34  - inside : same but not reachable.
35 
36  Finally the wanted sets are combined into a cellSet 'selected'. Apart
37  from straightforward adding the contents there are a few extra rules to
38  make sure that the surface of the 'outside' of the mesh is singly
39  connected.
40 
41 Usage
42 
43  - selectCells [OPTIONS]
44 
45  @param -case <dir>\n
46  Case directory.
47 
48  @param -help \n
49  Display help message.
50 
51  @param -doc \n
52  Display Doxygen API documentation page for this application.
53 
54  @param -srcDoc \n
55  Display Doxygen source documentation page for this application.
56 
57 \*---------------------------------------------------------------------------*/
58 
59 #include <OpenFOAM/argList.H>
60 #include <OpenFOAM/Time.H>
61 #include <OpenFOAM/polyMesh.H>
62 #include <OpenFOAM/IOdictionary.H>
64 #include <OpenFOAM/OFstream.H>
65 #include <meshTools/meshTools.H>
66 
67 #include <triSurface/triSurface.H>
69 #include <meshTools/meshSearch.H>
71 #include <meshTools/cellSet.H>
72 #include <meshTools/cellInfo.H>
73 #include <OpenFOAM/MeshWave.H>
74 #include "edgeStats.H"
77 
78 using namespace Foam;
79 
80 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
81 
82 // cellType for cells included/not included in mesh.
83 static const label MESH = cellClassification::INSIDE;
84 static const label NONMESH = cellClassification::OUTSIDE;
85 
86 
87 void writeSet(const cellSet& cells, const string& msg)
88 {
89  Info<< "Writing " << msg << " (" << cells.size() << ") to cellSet "
90  << cells.instance()/cells.local()/cells.name()
91  << endl << endl;
92  cells.write();
93 }
94 
95 
96 void getType(const labelList& elems, const label type, labelHashSet& set)
97 {
98  forAll(elems, i)
99  {
100  if (elems[i] == type)
101  {
102  set.insert(i);
103  }
104  }
105 }
106 
107 
108 void cutBySurface
109 (
110  const polyMesh& mesh,
111  const meshSearch& queryMesh,
112  const triSurfaceSearch& querySurf,
113 
114  const pointField& outsidePts,
115  const bool selectCut,
116  const bool selectInside,
117  const bool selectOutside,
118  const scalar nearDist,
119 
120  cellClassification& cellType
121 )
122 {
123  // Cut with surface and classify as inside/outside/cut
124  cellType =
126  (
127  mesh,
128  queryMesh,
129  querySurf,
130  outsidePts
131  );
132 
133  // Get inside/outside/cutCells cellSets.
134  cellSet inside(mesh, "inside", mesh.nCells()/10);
135  getType(cellType, cellClassification::INSIDE, inside);
136  writeSet(inside, "inside cells");
137 
138  cellSet outside(mesh, "outside", mesh.nCells()/10);
139  getType(cellType, cellClassification::OUTSIDE, outside);
140  writeSet(outside, "outside cells");
141 
142  cellSet cutCells(mesh, "cutCells", mesh.nCells()/10);
143  getType(cellType, cellClassification::CUT, cutCells);
144  writeSet(cutCells, "cells cut by surface");
145 
146 
147  // Change cellType to reflect selected part of mesh. Use
148  // MESH to denote selected part, NONMESH for all
149  // other cells.
150  // Is a bit of a hack but allows us to reuse all the functionality
151  // in cellClassification.
152 
153  forAll(cellType, cellI)
154  {
155  label cType = cellType[cellI];
156 
157  if (cType == cellClassification::CUT)
158  {
159  if (selectCut)
160  {
161  cellType[cellI] = MESH;
162  }
163  else
164  {
165  cellType[cellI] = NONMESH;
166  }
167  }
168  else if (cType == cellClassification::INSIDE)
169  {
170  if (selectInside)
171  {
172  cellType[cellI] = MESH;
173  }
174  else
175  {
176  cellType[cellI] = NONMESH;
177  }
178  }
179  else if (cType == cellClassification::OUTSIDE)
180  {
181  if (selectOutside)
182  {
183  cellType[cellI] = MESH;
184  }
185  else
186  {
187  cellType[cellI] = NONMESH;
188  }
189  }
190  else
191  {
192  FatalErrorIn("cutBySurface")
193  << "Multiple mesh regions in original mesh" << endl
194  << "Please use splitMeshRegions to separate these"
195  << exit(FatalError);
196  }
197  }
198 
199 
200  if (nearDist > 0)
201  {
202  Info<< "Removing cells with points closer than " << nearDist
203  << " to the surface ..." << nl << endl;
204 
205  const pointField& pts = mesh.points();
206  const indexedOctree<treeDataTriSurface>& tree = querySurf.tree();
207 
208  label nRemoved = 0;
209 
210  forAll(pts, pointI)
211  {
212  const point& pt = pts[pointI];
213 
214  pointIndexHit hitInfo = tree.findNearest(pt, sqr(nearDist));
215 
216  if (hitInfo.hit())
217  {
218  const labelList& pCells = mesh.pointCells()[pointI];
219 
220  forAll(pCells, i)
221  {
222  if (cellType[pCells[i]] != NONMESH)
223  {
224  cellType[pCells[i]] = NONMESH;
225  nRemoved++;
226  }
227  }
228  }
229  }
230 
231 // tmp<pointField> tnearest = querySurf.calcNearest(pts);
232 // const pointField& nearest = tnearest();
233 //
234 // label nRemoved = 0;
235 //
236 // forAll(nearest, pointI)
237 // {
238 // if (mag(nearest[pointI] - pts[pointI]) < nearDist)
239 // {
240 // const labelList& pCells = mesh.pointCells()[pointI];
241 //
242 // forAll(pCells, i)
243 // {
244 // if (cellType[pCells[i]] != NONMESH)
245 // {
246 // cellType[pCells[i]] = NONMESH;
247 // nRemoved++;
248 // }
249 // }
250 // }
251 // }
252 
253  Info<< "Removed " << nRemoved << " cells since too close to surface"
254  << nl << endl;
255  }
256 }
257 
258 
259 
260 // We're meshing the outside. Subset the currently selected mesh cells with the
261 // ones reachable from the outsidepoints.
262 label selectOutsideCells
263 (
264  const polyMesh& mesh,
265  const meshSearch& queryMesh,
266  const pointField& outsidePts,
267  cellClassification& cellType
268 )
269 {
270  //
271  // Check all outsidePts and for all of them inside a mesh cell
272  // collect the faces to start walking from
273  //
274 
275  // Outside faces
276  labelHashSet outsideFacesMap(outsidePts.size() * 6 * 2);
277  DynamicList<label> outsideFaces(outsideFacesMap.size());
278  // CellInfo on outside faces
279  DynamicList<cellInfo> outsideFacesInfo(outsideFacesMap.size());
280 
281  // cellInfo for mesh cell
282  const cellInfo meshInfo(MESH);
283 
284  forAll(outsidePts, outsidePtI)
285  {
286  // Find cell containing point. Linear search.
287  label cellI = queryMesh.findCell(outsidePts[outsidePtI], -1, false);
288 
289  if (cellI != -1 && cellType[cellI] == MESH)
290  {
291  Info<< "Marking cell " << cellI << " containing outside point "
292  << outsidePts[outsidePtI] << " with type " << cellType[cellI]
293  << " ..." << endl;
294 
295  //
296  // Mark this cell and its faces to start walking from
297  //
298 
299  // Mark faces of cellI
300  const labelList& cFaces = mesh.cells()[cellI];
301  forAll(cFaces, i)
302  {
303  label faceI = cFaces[i];
304 
305  if (outsideFacesMap.insert(faceI))
306  {
307  outsideFaces.append(faceI);
308  outsideFacesInfo.append(meshInfo);
309  }
310  }
311  }
312  }
313 
314  // Floodfill starting from outsideFaces (of type meshInfo)
315  MeshWave<cellInfo> regionCalc
316  (
317  mesh,
318  outsideFaces.shrink(),
319  outsideFacesInfo.shrink(),
320  mesh.nCells() // max iterations
321  );
322 
323  // Now regionCalc should hold info on cells that are reachable from
324  // changedFaces. Use these to subset cellType
325  const List<cellInfo>& allCellInfo = regionCalc.allCellInfo();
326 
327  label nChanged = 0;
328 
329  forAll(allCellInfo, cellI)
330  {
331  if (cellType[cellI] == MESH)
332  {
333  // Original cell was selected for meshing. Check if cell was
334  // reached from outsidePoints
335  if (allCellInfo[cellI].type() != MESH)
336  {
337  cellType[cellI] = NONMESH;
338  nChanged++;
339  }
340  }
341  }
342 
343  return nChanged;
344 }
345 
346 
347 // Main program:
348 
349 int main(int argc, char *argv[])
350 {
352 
353 # include <OpenFOAM/setRootCase.H>
354 # include <OpenFOAM/createTime.H>
355 # include <OpenFOAM/createPolyMesh.H>
356 
357  // Mesh edge statistics calculator
358  edgeStats edgeCalc(mesh);
359 
360 
361  IOdictionary refineDict
362  (
363  IOobject
364  (
365  "selectCellsDict",
366  runTime.system(),
367  mesh,
370  )
371  );
372 
373  fileName surfName(refineDict.lookup("surface"));
374  pointField outsidePts(refineDict.lookup("outsidePoints"));
375  bool useSurface(readBool(refineDict.lookup("useSurface")));
376  bool selectCut(readBool(refineDict.lookup("selectCut")));
377  bool selectInside(readBool(refineDict.lookup("selectInside")));
378  bool selectOutside(readBool(refineDict.lookup("selectOutside")));
379  scalar nearDist(readScalar(refineDict.lookup("nearDistance")));
380 
381 
382  if (useSurface)
383  {
384  Info<< "Cells to be used for meshing (0=false, 1=true):" << nl
385  << " cells cut by surface : " << selectCut << nl
386  << " cells inside of surface : " << selectInside << nl
387  << " cells outside of surface : " << selectOutside << nl
388  << " cells with points further than : " << nearDist << nl
389  << endl;
390  }
391  else
392  {
393  Info<< "Cells to be used for meshing (0=false, 1=true):" << nl
394  << " cells reachable from outsidePoints:" << selectOutside << nl
395  << endl;
396  }
397 
398  // Print edge stats on original mesh.
399  (void)edgeCalc.minLen(Info);
400 
401  // Search engine on mesh. Face decomposition since faces might be warped.
402  meshSearch queryMesh(mesh, true);
403 
404  // Check all 'outside' points
405  forAll(outsidePts, outsideI)
406  {
407  const point& outsidePoint = outsidePts[outsideI];
408 
409  label cellI = queryMesh.findCell(outsidePoint, -1, false);
410  if (returnReduce(cellI, maxOp<label>()) == -1)
411  {
413  << "outsidePoint " << outsidePoint
414  << " is not inside any cell"
415  << exit(FatalError);
416  }
417  }
418 
419  // Cell status (compared to surface if provided): inside/outside/cut.
420  // Start off from everything selected and cut later.
421  cellClassification cellType
422  (
423  mesh,
424  labelList
425  (
426  mesh.nCells(),
428  )
429  );
430 
431 
432  // Surface
433  autoPtr<triSurface> surf(NULL);
434  // Search engine on surface.
435  autoPtr<triSurfaceSearch> querySurf(NULL);
436 
437  if (useSurface)
438  {
439  surf.reset(new triSurface(surfName));
440 
441  // Dump some stats
442  surf().writeStats(Info);
443 
444  // Search engine on surface.
445  querySurf.reset(new triSurfaceSearch(surf));
446 
447  // Set cellType[cellI] according to relation to surface
448  cutBySurface
449  (
450  mesh,
451  queryMesh,
452  querySurf,
453 
454  outsidePts,
455  selectCut,
456  selectInside,
457  selectOutside,
458  nearDist,
459 
460  cellType
461  );
462  }
463 
464 
465  // Now 'trim' all the corners from the mesh so meshing/surface extraction
466  // becomes easier.
467 
468  label nHanging, nRegionEdges, nRegionPoints, nOutside;
469 
470  do
471  {
472  Info<< "Removing cells which after subsetting would have all points"
473  << " on outside ..." << nl << endl;
474 
475  nHanging = cellType.fillHangingCells
476  (
477  MESH, // meshType
478  NONMESH, // fill type
479  mesh.nCells()
480  );
481 
482 
483  Info<< "Removing edges connecting cells unconnected by faces ..."
484  << nl << endl;
485 
486  nRegionEdges = cellType.fillRegionEdges
487  (
488  MESH, // meshType
489  NONMESH, // fill type
490  mesh.nCells()
491  );
492 
493 
494  Info<< "Removing points connecting cells unconnected by faces ..."
495  << nl << endl;
496 
497  nRegionPoints = cellType.fillRegionPoints
498  (
499  MESH, // meshType
500  NONMESH, // fill type
501  mesh.nCells()
502  );
503 
504  nOutside = 0;
505  if (selectOutside)
506  {
507  // Since we're selecting the cells reachable from outsidePoints
508  // and the set might have changed, redo the outsideCells
509  // calculation
510  nOutside = selectOutsideCells
511  (
512  mesh,
513  queryMesh,
514  outsidePts,
515  cellType
516  );
517  }
518  } while
519  (
520  nHanging != 0
521  || nRegionEdges != 0
522  || nRegionPoints != 0
523  || nOutside != 0
524  );
525 
526  cellSet selectedCells(mesh, "selected", mesh.nCells()/10);
527  getType(cellType, MESH, selectedCells);
528 
529  writeSet(selectedCells, "cells selected for meshing");
530 
531 
532  Info << "End\n" << endl;
533 
534  return 0;
535 }
536 
537 
538 // ************************ vim: set sw=4 sts=4 et: ************************ //