regionToCell.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 "regionToCell.H"
27 #include "regionSplit.H"
28 #include "emptyPolyPatch.H"
29 #include "cellSet.H"
30 #include "syncTools.H"
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37 
38 defineTypeNameAndDebug(regionToCell, 0);
39 
40 addToRunTimeSelectionTable(topoSetSource, regionToCell, word);
41 
42 addToRunTimeSelectionTable(topoSetSource, regionToCell, istream);
43 
44 }
45 
46 
47 Foam::topoSetSource::addToUsageTable Foam::regionToCell::usage_
48 (
49  regionToCell::typeName,
50  "\n Usage: regionToCell subCellSet (pt0 .. ptn)\n\n"
51  " Select all cells in the connected region containing"
52  " points (pt0..ptn).\n"
53 );
54 
55 
56 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
57 
58 void Foam::regionToCell::markRegionFaces
59 (
60  const boolList& selectedCell,
61  boolList& regionFace
62 ) const
63 {
64  // Internal faces
65  const labelList& faceOwner = mesh_.faceOwner();
66  const labelList& faceNeighbour = mesh_.faceNeighbour();
67  forAll(faceNeighbour, facei)
68  {
69  if
70  (
71  selectedCell[faceOwner[facei]]
72  != selectedCell[faceNeighbour[facei]]
73  )
74  {
75  regionFace[facei] = true;
76  }
77  }
78 
79  // Swap neighbour selectedCell state
80  boolList nbrSelected;
81  syncTools::swapBoundaryCellList(mesh_, selectedCell, nbrSelected);
82 
83  // Boundary faces
84  const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
85  forAll(pbm, patchi)
86  {
87  const polyPatch& pp = pbm[patchi];
88  const labelUList& faceCells = pp.faceCells();
89  forAll(faceCells, i)
90  {
91  label facei = pp.start()+i;
92  label bFacei = facei-mesh_.nInternalFaces();
93  if
94  (
95  selectedCell[faceCells[i]]
96  != selectedCell[nbrSelected[bFacei]]
97  )
98  {
99  regionFace[facei] = true;
100  }
101  }
102  }
103 }
104 
105 
106 Foam::boolList Foam::regionToCell::findRegions
107 (
108  const bool verbose,
109  const regionSplit& cellRegion
110 ) const
111 {
112  boolList keepRegion(cellRegion.nRegions(), false);
113 
114  forAll(insidePoints_, i)
115  {
116  // Find the region containing the insidePoint
117 
118  label celli = mesh_.findCell(insidePoints_[i]);
119 
120  label keepRegionI = -1;
121  label keepProci = -1;
122  if (celli != -1)
123  {
124  keepRegionI = cellRegion[celli];
125  keepProci = Pstream::myProcNo();
126  }
127  reduce(keepRegionI, maxOp<label>());
128  keepRegion[keepRegionI] = true;
129 
130  reduce(keepProci, maxOp<label>());
131 
132  if (keepProci == -1)
133  {
135  << "Did not find " << insidePoints_[i]
136  << " in mesh." << " Mesh bounds are " << mesh_.bounds()
137  << exit(FatalError);
138  }
139 
140  if (verbose)
141  {
142  Info<< " Found location " << insidePoints_[i]
143  << " in cell " << celli << " on processor " << keepProci
144  << " in global region " << keepRegionI
145  << " out of " << cellRegion.nRegions() << " regions." << endl;
146  }
147  }
148 
149  return keepRegion;
150 }
151 
152 
153 void Foam::regionToCell::unselectOutsideRegions
154 (
155  boolList& selectedCell
156 ) const
157 {
158  // Determine faces on the edge of selectedCell
159  boolList blockedFace(mesh_.nFaces(), false);
160  markRegionFaces(selectedCell, blockedFace);
161 
162  // Determine regions
163  regionSplit cellRegion(mesh_, blockedFace);
164 
165  // Determine regions containing insidePoints_
166  boolList keepRegion(findRegions(true, cellRegion));
167 
168  // Go back to bool per cell
169  forAll(cellRegion, celli)
170  {
171  if (!keepRegion[cellRegion[celli]])
172  {
173  selectedCell[celli] = false;
174  }
175  }
176 }
177 
178 
179 void Foam::regionToCell::shrinkRegions
180 (
181  boolList& selectedCell
182 ) const
183 {
184  // Select points on unselected cells and boundary
185  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
186 
187  boolList boundaryPoint(mesh_.nPoints(), false);
188 
189  const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
190 
191  forAll(pbm, patchi)
192  {
193  const polyPatch& pp = pbm[patchi];
194 
195  if (!pp.coupled() && !isA<emptyPolyPatch>(pp))
196  {
197  forAll(pp, i)
198  {
199  const face& f = pp[i];
200  forAll(f, fp)
201  {
202  boundaryPoint[f[fp]] = true;
203  }
204  }
205  }
206  }
207 
208  forAll(selectedCell, celli)
209  {
210  if (!selectedCell[celli])
211  {
212  const labelList& cPoints = mesh_.cellPoints(celli);
213  forAll(cPoints, i)
214  {
215  boundaryPoint[cPoints[i]] = true;
216  }
217  }
218  }
219 
220  syncTools::syncPointList(mesh_, boundaryPoint, orEqOp<bool>(), false);
221 
222 
223  // Select all cells using these points
224  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
225 
226  label nChanged = 0;
227  forAll(boundaryPoint, pointi)
228  {
229  if (boundaryPoint[pointi])
230  {
231  const labelList& pCells = mesh_.pointCells(pointi);
232  forAll(pCells, i)
233  {
234  label celli = pCells[i];
235  if (selectedCell[celli])
236  {
237  selectedCell[celli] = false;
238  nChanged++;
239  }
240  }
241  }
242  }
243  Info<< " Eroded " << returnReduce(nChanged, sumOp<label>())
244  << " cells." << endl;
245 }
246 
247 
248 void Foam::regionToCell::erode
249 (
250  boolList& selectedCell
251 ) const
252 {
253  //Info<< "Entering shrinkRegions:" << count(selectedCell) << endl;
254  //generateField("selectedCell_before", selectedCell)().write();
255 
256  // Now erode and see which regions get disconnected
257  boolList shrunkSelectedCell(selectedCell);
258 
259  for (label iter = 0; iter < nErode_; iter++)
260  {
261  shrinkRegions(shrunkSelectedCell);
262  }
263 
264  //Info<< "After shrinking:" << count(shrunkSelectedCell) << endl;
265  //generateField("shrunkSelectedCell", shrunkSelectedCell)().write();
266 
267 
268 
269  // Determine faces on the edge of shrunkSelectedCell
270  boolList blockedFace(mesh_.nFaces(), false);
271  markRegionFaces(shrunkSelectedCell, blockedFace);
272 
273  // Find disconnected regions
274  regionSplit cellRegion(mesh_, blockedFace);
275 
276  // Determine regions containing insidePoints
277  boolList keepRegion(findRegions(true, cellRegion));
278 
279 
280  // Extract cells in regions that are not to be kept.
281  boolList removeCell(mesh_.nCells(), false);
282  forAll(cellRegion, celli)
283  {
284  if (shrunkSelectedCell[celli] && !keepRegion[cellRegion[celli]])
285  {
286  removeCell[celli] = true;
287  }
288  }
289 
290  //Info<< "removeCell before:" << count(removeCell) << endl;
291  //generateField("removeCell_before", removeCell)().write();
292 
293 
294 
295  // Grow removeCell
296  for (label iter = 0; iter < nErode_; iter++)
297  {
298  // Grow selected cell in regions that are not for keeping
299  boolList boundaryPoint(mesh_.nPoints(), false);
300  forAll(removeCell, celli)
301  {
302  if (removeCell[celli])
303  {
304  const labelList& cPoints = mesh_.cellPoints(celli);
305  forAll(cPoints, i)
306  {
307  boundaryPoint[cPoints[i]] = true;
308  }
309  }
310  }
311  syncTools::syncPointList(mesh_, boundaryPoint, orEqOp<bool>(), false);
312 
313  // Select all cells using these points
314 
315  label nChanged = 0;
316  forAll(boundaryPoint, pointi)
317  {
318  if (boundaryPoint[pointi])
319  {
320  const labelList& pCells = mesh_.pointCells(pointi);
321  forAll(pCells, i)
322  {
323  label celli = pCells[i];
324  if (!removeCell[celli])
325  {
326  removeCell[celli] = true;
327  nChanged++;
328  }
329  }
330  }
331  }
332  }
333 
334  //Info<< "removeCell after:" << count(removeCell) << endl;
335  //generateField("removeCell_after", removeCell)().write();
336 
337 
338  // Unmark removeCell
339  forAll(removeCell, celli)
340  {
341  if (removeCell[celli])
342  {
343  selectedCell[celli] = false;
344  }
345  }
346 }
347 
348 
349 void Foam::regionToCell::combine(topoSet& set, const bool add) const
350 {
351  // Note: wip. Select cells first
352  boolList selectedCell(mesh_.nCells(), true);
353 
354  if (setName_.size() && setName_ != "none")
355  {
356  Info<< " Loading subset " << setName_ << " to delimit search region."
357  << endl;
358  cellSet subSet(mesh_, setName_);
359 
360  selectedCell = false;
361  forAllConstIter(cellSet, subSet, iter)
362  {
363  selectedCell[iter.key()] = true;
364  }
365  }
366 
367 
368  unselectOutsideRegions(selectedCell);
369 
370  if (nErode_ > 0)
371  {
372  erode(selectedCell);
373  }
374 
375 
376  forAll(selectedCell, celli)
377  {
378  if (selectedCell[celli])
379  {
380  addOrDelete(set, celli, add);
381  }
382  }
383 }
384 
385 
386 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
387 
388 // Construct from components
390 (
391  const polyMesh& mesh,
392  const word& setName,
393  const pointField& insidePoints,
394  const label nErode
395 )
396 :
397  topoSetSource(mesh),
398  setName_(setName),
399  insidePoints_(insidePoints),
400  nErode_(nErode)
401 {}
402 
403 
404 // Construct from dictionary
406 (
407  const polyMesh& mesh,
408  const dictionary& dict
409 )
410 :
411  topoSetSource(mesh),
412  setName_(dict.lookupOrDefault<word>("set", "none")),
413  insidePoints_
414  (
415  dict.found("insidePoints")
416  ? dict.lookup("insidePoints")
417  : dict.lookup("insidePoint")
418  ),
419  nErode_(dict.lookupOrDefault("nErode", 0))
420 {}
421 
422 
423 // Construct from Istream
425 (
426  const polyMesh& mesh,
427  Istream& is
428 )
429 :
430  topoSetSource(mesh),
431  setName_(checkIs(is)),
432  insidePoints_(checkIs(is)),
433  nErode_(0)
434 {}
435 
436 
437 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
438 
440 {}
441 
442 
443 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
444 
446 (
447  const topoSetSource::setAction action,
448  topoSet& set
449 ) const
450 {
451  if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
452  {
453  Info<< " Adding all cells of connected region containing points "
454  << insidePoints_ << " ..." << endl;
455 
456  combine(set, true);
457  }
458  else if (action == topoSetSource::DELETE)
459  {
460  Info<< " Removing all cells of connected region containing points "
461  << insidePoints_ << " ..." << endl;
462 
463  combine(set, false);
464  }
465 }
466 
467 
468 // ************************************************************************* //
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:431
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
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
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:418
virtual void applyToSet(const topoSetSource::setAction action, topoSet &) const
Definition: regionToCell.C:446
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
AccessType combine(const List< T > &, AccessOp aop=accessOp< T >())
Combines sublists into one big list.
Definition: ListListOps.C:34
virtual ~regionToCell()
Destructor.
Definition: regionToCell.C:439
Macros for easy insertion into run-time selection tables.
UList< label > labelUList
Definition: UList.H:64
regionToCell(const polyMesh &mesh, const word &setName, const pointField &insidePoints, const label nErode)
Construct from components.
Definition: regionToCell.C:390
Base class of a source for a topoSet.
Definition: topoSetSource.H:63
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
A class for handling words, derived from string.
Definition: word.H:59
setAction
Enumeration defining the valid actions.
Definition: topoSetSource.H:82
List< label > labelList
A List of labels.
Definition: labelList.H:56
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
labelList f(nPoints)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
General set of labels of mesh quantity (points, cells, faces).
Definition: topoSet.H:61
static void swapBoundaryCellList(const polyMesh &mesh, const UList< T > &cellData, List< T > &neighbourCellData)
Swap to obtain neighbour cell values for all boundary faces.
label patchi
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.
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
messageStream Info
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
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