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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26Application
27 selectCells
28
29Group
30 grpMeshAdvancedUtilities
31
32Description
33 Select cells in relation to surface.
34
35 Divides cells into three sets:
36 - cutCells : cells cut by surface or close to surface.
37 - outside : cells not in cutCells and reachable from set of
38 user-defined points (outsidePoints)
39 - inside : same but not reachable.
40
41 Finally the wanted sets are combined into a cellSet 'selected'. Apart
42 from straightforward adding the contents there are a few extra rules to
43 make sure that the surface of the 'outside' of the mesh is singly
44 connected.
45
46\*---------------------------------------------------------------------------*/
47
48#include "argList.H"
49#include "Time.H"
50#include "polyMesh.H"
51#include "IOdictionary.H"
52#include "twoDPointCorrector.H"
53#include "OFstream.H"
54#include "meshTools.H"
55
56#include "triSurface.H"
57#include "triSurfaceSearch.H"
58#include "meshSearch.H"
59#include "cellClassification.H"
60#include "cellSet.H"
61#include "cellInfo.H"
62#include "MeshWave.H"
63#include "edgeStats.H"
64#include "treeDataTriSurface.H"
65#include "indexedOctree.H"
66#include "globalMeshData.H"
67
68using namespace Foam;
69
70// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
71
72// cellType for cells included/not included in mesh.
73static const label MESH = cellClassification::INSIDE;
74static const label NONMESH = cellClassification::OUTSIDE;
75
76
77void writeSet(const cellSet& cells, const string& msg)
78{
79 Info<< "Writing " << msg << " (" << cells.size() << ") to cellSet "
80 << cells.instance()/cells.local()/cells.name()
81 << endl << endl;
82 cells.write();
83}
84
85
86void getType(const labelList& elems, const label type, labelHashSet& set)
87{
88 forAll(elems, i)
89 {
90 if (elems[i] == type)
91 {
92 set.insert(i);
93 }
94 }
95}
96
97
98void cutBySurface
99(
100 const polyMesh& mesh,
101 const meshSearch& queryMesh,
102 const triSurfaceSearch& querySurf,
103
104 const pointField& outsidePts,
105 const bool selectCut,
106 const bool selectInside,
107 const bool selectOutside,
108 const scalar nearDist,
109
110 cellClassification& cellType
111)
112{
113 // Cut with surface and classify as inside/outside/cut
114 cellType =
116 (
117 mesh,
118 queryMesh,
119 querySurf,
120 outsidePts
121 );
122
123 // Get inside/outside/cutCells cellSets.
124 cellSet inside(mesh, "inside", mesh.nCells()/10);
125 getType(cellType, cellClassification::INSIDE, inside);
126 writeSet(inside, "inside cells");
127
128 cellSet outside(mesh, "outside", mesh.nCells()/10);
129 getType(cellType, cellClassification::OUTSIDE, outside);
130 writeSet(outside, "outside cells");
131
132 cellSet cutCells(mesh, "cutCells", mesh.nCells()/10);
133 getType(cellType, cellClassification::CUT, cutCells);
134 writeSet(cutCells, "cells cut by surface");
135
136
137 // Change cellType to reflect selected part of mesh. Use
138 // MESH to denote selected part, NONMESH for all
139 // other cells.
140 // Is a bit of a hack but allows us to reuse all the functionality
141 // in cellClassification.
142
143 forAll(cellType, celli)
144 {
145 label cType = cellType[celli];
146
147 if (cType == cellClassification::CUT)
148 {
149 if (selectCut)
150 {
151 cellType[celli] = MESH;
152 }
153 else
154 {
155 cellType[celli] = NONMESH;
156 }
157 }
158 else if (cType == cellClassification::INSIDE)
159 {
160 if (selectInside)
161 {
162 cellType[celli] = MESH;
163 }
164 else
165 {
166 cellType[celli] = NONMESH;
167 }
168 }
169 else if (cType == cellClassification::OUTSIDE)
170 {
171 if (selectOutside)
172 {
173 cellType[celli] = MESH;
174 }
175 else
176 {
177 cellType[celli] = NONMESH;
178 }
179 }
180 else
181 {
183 << "Multiple mesh regions in original mesh" << endl
184 << "Please use splitMeshRegions to separate these"
185 << exit(FatalError);
186 }
187 }
188
189
190 if (nearDist > 0)
191 {
192 Info<< "Removing cells with points closer than " << nearDist
193 << " to the surface ..." << nl << endl;
194
195 const pointField& pts = mesh.points();
196 const indexedOctree<treeDataTriSurface>& tree = querySurf.tree();
197
198 label nRemoved = 0;
199
200 forAll(pts, pointi)
201 {
202 const point& pt = pts[pointi];
203
204 pointIndexHit hitInfo = tree.findNearest(pt, sqr(nearDist));
205
206 if (hitInfo.hit())
207 {
208 const labelList& pCells = mesh.pointCells()[pointi];
209
210 forAll(pCells, i)
211 {
212 if (cellType[pCells[i]] != NONMESH)
213 {
214 cellType[pCells[i]] = NONMESH;
215 nRemoved++;
216 }
217 }
218 }
219 }
220
221// tmp<pointField> tnearest = querySurf.calcNearest(pts);
222// const pointField& nearest = tnearest();
223//
224// label nRemoved = 0;
225//
226// forAll(nearest, pointi)
227// {
228// if (mag(nearest[pointi] - pts[pointi]) < nearDist)
229// {
230// const labelList& pCells = mesh.pointCells()[pointi];
231//
232// forAll(pCells, i)
233// {
234// if (cellType[pCells[i]] != NONMESH)
235// {
236// cellType[pCells[i]] = NONMESH;
237// nRemoved++;
238// }
239// }
240// }
241// }
242
243 Info<< "Removed " << nRemoved << " cells since too close to surface"
244 << nl << endl;
245 }
246}
247
248
249
250// We're meshing the outside. Subset the currently selected mesh cells with the
251// ones reachable from the outsidepoints.
252label selectOutsideCells
253(
254 const polyMesh& mesh,
255 const meshSearch& queryMesh,
256 const pointField& outsidePts,
257 cellClassification& cellType
258)
259{
260 //
261 // Check all outsidePts and for all of them inside a mesh cell
262 // collect the faces to start walking from
263 //
264
265 // Outside faces
266 labelHashSet outsideFacesMap(outsidePts.size() * 6 * 2);
267 DynamicList<label> outsideFaces(outsideFacesMap.size());
268 // CellInfo on outside faces
269 DynamicList<cellInfo> outsideFacesInfo(outsideFacesMap.size());
270
271 // cellInfo for mesh cell
272 const cellInfo meshInfo(MESH);
273
274 forAll(outsidePts, outsidePtI)
275 {
276 // Find cell containing point. Linear search.
277 label celli = queryMesh.findCell(outsidePts[outsidePtI], -1, false);
278
279 if (celli != -1 && cellType[celli] == MESH)
280 {
281 Info<< "Marking cell " << celli << " containing outside point "
282 << outsidePts[outsidePtI] << " with type " << cellType[celli]
283 << " ..." << endl;
284
285 //
286 // Mark this cell and its faces to start walking from
287 //
288
289 // Mark faces of celli
290 const labelList& cFaces = mesh.cells()[celli];
291 forAll(cFaces, i)
292 {
293 label facei = cFaces[i];
294
295 if (outsideFacesMap.insert(facei))
296 {
297 outsideFaces.append(facei);
298 outsideFacesInfo.append(meshInfo);
299 }
300 }
301 }
302 }
303
304 // Floodfill starting from outsideFaces (of type meshInfo)
305 MeshWave<cellInfo> regionCalc
306 (
307 mesh,
308 outsideFaces.shrink(),
309 outsideFacesInfo.shrink(),
310 mesh.globalData().nTotalCells()+1 // max iterations
311 );
312
313 // Now regionCalc should hold info on cells that are reachable from
314 // changedFaces. Use these to subset cellType
315 const List<cellInfo>& allCellInfo = regionCalc.allCellInfo();
316
317 label nChanged = 0;
318
319 forAll(allCellInfo, celli)
320 {
321 if (cellType[celli] == MESH)
322 {
323 // Original cell was selected for meshing. Check if cell was
324 // reached from outsidePoints
325 if (allCellInfo[celli].type() != MESH)
326 {
327 cellType[celli] = NONMESH;
328 nChanged++;
329 }
330 }
331 }
332
333 return nChanged;
334}
335
336
337
338int main(int argc, char *argv[])
339{
340 argList::addNote
341 (
342 "Select cells in relation to surface"
343 );
344
345 argList::noParallel();
346
347 #include "setRootCase.H"
348 #include "createTime.H"
349 #include "createPolyMesh.H"
350
351 // Mesh edge statistics calculator
352 edgeStats edgeCalc(mesh);
353
354
355 IOdictionary refineDict
356 (
358 (
359 "selectCellsDict",
360 runTime.system(),
361 mesh,
362 IOobject::MUST_READ_IF_MODIFIED,
363 IOobject::NO_WRITE
364 )
365 );
366
367 fileName surfName(refineDict.get<fileName>("surface"));
368 pointField outsidePts(refineDict.lookup("outsidePoints"));
369 const bool useSurface(refineDict.get<bool>("useSurface"));
370 const bool selectCut(refineDict.get<bool>("selectCut"));
371 const bool selectInside(refineDict.get<bool>("selectInside"));
372 const bool selectOutside(refineDict.get<bool>("selectOutside"));
373 const scalar nearDist(refineDict.get<scalar>("nearDistance"));
374
375
376 if (useSurface)
377 {
378 Info<< "Cells to be used for meshing (0=false, 1=true):" << nl
379 << " cells cut by surface : " << selectCut << nl
380 << " cells inside of surface : " << selectInside << nl
381 << " cells outside of surface : " << selectOutside << nl
382 << " cells with points further than : " << nearDist << nl
383 << endl;
384 }
385 else
386 {
387 Info<< "Cells to be used for meshing (0=false, 1=true):" << nl
388 << " cells reachable from outsidePoints:" << selectOutside << nl
389 << endl;
390 }
391
392 // Print edge stats on original mesh.
393 (void)edgeCalc.minLen(Info);
394
395 // Search engine on mesh. Face decomposition since faces might be warped.
396 meshSearch queryMesh(mesh);
397
398 // Check all 'outside' points
399 forAll(outsidePts, outsideI)
400 {
401 const point& outsidePoint = outsidePts[outsideI];
402
403 label celli = queryMesh.findCell(outsidePoint, -1, false);
404 if (returnReduce(celli, maxOp<label>()) == -1)
405 {
407 << "outsidePoint " << outsidePoint
408 << " is not inside any cell"
409 << exit(FatalError);
410 }
411 }
412
413 // Cell status (compared to surface if provided): inside/outside/cut.
414 // Start off from everything selected and cut later.
416 (
417 mesh,
419 (
420 mesh.nCells(),
421 cellClassification::MESH
422 )
423 );
424
425
426 if (useSurface)
427 {
428 triSurface surf(surfName);
429
430 // Dump some stats
431 surf.writeStats(Info);
432
433 triSurfaceSearch querySurf(surf);
434
435 // Set cellType[celli] according to relation to surface
436 cutBySurface
437 (
438 mesh,
439 queryMesh,
440 querySurf,
441
442 outsidePts,
443 selectCut,
444 selectInside,
445 selectOutside,
446 nearDist,
447
448 cellType
449 );
450 }
451
452
453 // Now 'trim' all the corners from the mesh so meshing/surface extraction
454 // becomes easier.
455
456 label nHanging, nRegionEdges, nRegionPoints, nOutside;
457
458 do
459 {
460 Info<< "Removing cells which after subsetting would have all points"
461 << " on outside ..." << nl << endl;
462
463 nHanging = cellType.fillHangingCells
464 (
465 MESH, // meshType
466 NONMESH, // fill type
467 mesh.nCells()
468 );
469
470
471 Info<< "Removing edges connecting cells unconnected by faces ..."
472 << nl << endl;
473
474 nRegionEdges = cellType.fillRegionEdges
475 (
476 MESH, // meshType
477 NONMESH, // fill type
478 mesh.nCells()
479 );
480
481
482 Info<< "Removing points connecting cells unconnected by faces ..."
483 << nl << endl;
484
485 nRegionPoints = cellType.fillRegionPoints
486 (
487 MESH, // meshType
488 NONMESH, // fill type
489 mesh.nCells()
490 );
491
492 nOutside = 0;
493 if (selectOutside)
494 {
495 // Since we're selecting the cells reachable from outsidePoints
496 // and the set might have changed, redo the outsideCells
497 // calculation
498 nOutside = selectOutsideCells
499 (
500 mesh,
501 queryMesh,
502 outsidePts,
503 cellType
504 );
505 }
506 } while
507 (
508 nHanging != 0
509 || nRegionEdges != 0
510 || nRegionPoints != 0
511 || nOutside != 0
512 );
513
514 cellSet selectedCells(mesh, "selected", mesh.nCells()/10);
515 getType(cellType, MESH, selectedCells);
516
517 writeSet(selectedCells, "cells selected for meshing");
518
519
520 Info<< "End\n" << endl;
521
522 return 0;
523}
524
525
526// ************************************************************************* //
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
FaceCellWave plus data.
Definition: MeshWave.H:62
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:66
bool hit() const noexcept
Is there a hit?
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
'Cuts' a mesh with a surface.
Holds information regarding type of cell. Used in inside/outside determination in cellClassification.
Definition: cellInfo.H:65
A collection of cell labels.
Definition: cellSet.H:54
Helper class to calculate minimum edge length on mesh.
Definition: edgeStats.H:57
A class for handling file names.
Definition: fileName.H:76
Non-pointer based hierarchical recursive searching.
Definition: indexedOctree.H:74
Various (local, not parallel) searches on polyMesh; uses (demand driven) octree to search.
Definition: meshSearch.H:61
label findCell(const point &location, const label seedCelli=-1, const bool useTreeSearch=true) const
Find cell containing location.
Definition: meshSearch.C:780
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
Helper class to search on triSurface.
const indexedOctree< treeDataTriSurface > & tree() const
Demand driven construction of the octree.
Triangulated surface description with patch information.
Definition: triSurface.H:79
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const cellShapeList & cells
void set(List< bool > &bools, const labelUList &locations)
Set the listed locations (assign 'true').
Definition: BitOps.C:38
type
Types of root.
Definition: Roots.H:55
cellType
Equivalent to enumeration in "vtkCellType.h" (should be uint8_t)
Definition: foamVtkCore.H:90
Namespace for OpenFOAM.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333