planeToFaceZone.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2020 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 "planeToFaceZone.H"
27 #include "polyMesh.H"
28 #include "faceZoneSet.H"
30 #include "PatchTools.H"
31 #include "syncTools.H"
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38  template<>
40  {
41  "all",
42  "closest"
43  };
44 }
45 
48 
49 
50 namespace Foam
51 {
55 }
56 
57 
58 Foam::topoSetSource::addToUsageTable Foam::planeToFaceZone::usage_
59 (
60  planeToFaceZone::typeName,
61  "\n Usage: planeToFaceZone (px py pz) (nx ny nz) include\n\n"
62  " Select faces for which the adjacent cell centres lie on opposite "
63  " of a plane\n\n"
64 );
65 
66 
67 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
68 
69 void Foam::planeToFaceZone::combine(faceZoneSet& fzSet, const bool add) const
70 {
71  // Mark all cells with centres above the plane
72  boolList cellIsAbovePlane(mesh_.nCells());
73  forAll(mesh_.cells(), celli)
74  {
75  cellIsAbovePlane[celli] =
76  ((mesh_.cellCentres()[celli] - point_) & normal_) > 0;
77  }
78 
79  // Mark all faces that sit between cells above and below the plane
80  boolList faceIsOnPlane(mesh_.nFaces());
81  forAll(mesh_.faceNeighbour(), facei)
82  {
83  faceIsOnPlane[facei] =
84  cellIsAbovePlane[mesh_.faceOwner()[facei]]
85  != cellIsAbovePlane[mesh_.faceNeighbour()[facei]];
86  }
87  forAll(mesh_.boundaryMesh(), patchi)
88  {
89  const polyPatch& patch = mesh_.boundaryMesh()[patchi];
90  forAll(patch, patchFacei)
91  {
92  const label facei = patch.start() + patchFacei;
93  faceIsOnPlane[facei] =
94  patch.coupled() && cellIsAbovePlane[mesh_.faceOwner()[facei]];
95  }
96  }
97  syncTools::syncFaceList(mesh_, faceIsOnPlane, notEqOp<bool>());
98 
99  // Convert marked faces to a list of indices
100  labelList newSetFaces(findIndices(faceIsOnPlane, true));
101 
102  // If constructing a single contiguous set, remove all faces except those
103  // connected to the contiguous region closest to the specified point
104  if (include_ == include::closest)
105  {
106  // Step 1: Get locally contiguous regions for the new face set and the
107  // total number of regions across all processors.
108  labelList newSetFaceRegions(newSetFaces.size(), -1);
109  label nRegions = -1;
110  {
111  // Create a patch of the set faces
112  const uindirectPrimitivePatch newSetPatch
113  (
114  UIndirectList<face>(mesh_.faces(), newSetFaces),
115  mesh_.points()
116  );
117 
118  // Get the region ID-s and store the total number of regions on
119  // each processor
120  labelList procNRegions(Pstream::nProcs(), -1);
121  procNRegions[Pstream::myProcNo()] =
123  (
124  newSetPatch,
125  boolList(newSetPatch.nEdges(), false),
126  newSetFaceRegions
127  );
128  Pstream::gatherList(procNRegions);
129  Pstream::scatterList(procNRegions);
130 
131  // Cumulative sum the number of regions on each processor to get an
132  // offset which makes the local region ID-s globally unique
133  labelList procRegionOffset(Pstream::nProcs(), 0);
134  for (label proci = 1; proci < Pstream::nProcs(); ++ proci)
135  {
136  procRegionOffset[proci] +=
137  procRegionOffset[proci - 1]
138  + procNRegions[proci - 1];
139  }
140 
141  // Apply the offset
142  forAll(newSetFaces, newSetFacei)
143  {
144  newSetFaceRegions[newSetFacei] +=
145  procRegionOffset[Pstream::myProcNo()];
146  }
147 
148  // Store the total number of regions across all processors
149  nRegions = procRegionOffset.last() + procNRegions.last();
150  }
151 
152  // Step 2: Create a region map which combines regions which are
153  // connected across coupled interfaces
154  labelList regionMap(identity(nRegions));
155  {
156  // Put region labels on connected boundary edges and synchronise to
157  // create a list of all regions connected to a given edge
158  labelListList meshEdgeRegions(mesh_.nEdges(), labelList());
159  forAll(newSetFaces, newSetFacei)
160  {
161  const label facei = newSetFaces[newSetFacei];
162  const label regioni = newSetFaceRegions[newSetFacei];
163  forAll(mesh_.faceEdges()[facei], faceEdgei)
164  {
165  const label edgei = mesh_.faceEdges()[facei][faceEdgei];
166  meshEdgeRegions[edgei] = labelList(1, regioni);
167  }
168  }
170  (
171  mesh_,
172  meshEdgeRegions,
174  labelList()
175  );
176 
177  // Combine edge regions to create a list of what regions a given
178  // region is connected to
179  List<labelHashSet> regionRegions(nRegions);
180  forAll(newSetFaces, newSetFacei)
181  {
182  const label facei = newSetFaces[newSetFacei];
183  const label regioni = newSetFaceRegions[newSetFacei];
184  forAll(mesh_.faceEdges()[facei], faceEdgei)
185  {
186  const label edgei = mesh_.faceEdges()[facei][faceEdgei];
187  forAll(meshEdgeRegions[edgei], edgeRegioni)
188  {
189  if (meshEdgeRegions[edgei][edgeRegioni] != regioni)
190  {
191  regionRegions[regioni].insert
192  (
193  meshEdgeRegions[edgei][edgeRegioni]
194  );
195  }
196  }
197  }
198  }
200  Pstream::listCombineScatter(regionRegions);
201 
202  // Collapse the region connections into a map between each region
203  // and the lowest numbered region that it connects to
204  forAll(regionRegions, regioni)
205  {
206  forAllConstIter(labelHashSet, regionRegions[regioni], iter)
207  {
208  regionMap[iter.key()] =
209  min(regionMap[iter.key()], regionMap[regioni]);
210  }
211  }
212  }
213 
214  // Step 3: Combine connected regions
215  labelList regionNFaces;
216  {
217  // Remove duplicates from the region map
218  label regioni0 = 0;
219  forAll(regionMap, regioni)
220  {
221  if (regionMap[regioni] > regioni0)
222  {
223  ++ regioni0;
224  regionMap[regioni] = regioni0;
225  }
226  }
227 
228  // Recompute the number of regions
229  nRegions = regioni0 + 1;
230 
231  // Renumber the face region ID-s
232  newSetFaceRegions =
233  IndirectList<label>(regionMap, newSetFaceRegions);
234 
235  // Report the final number and size of the regions
236  regionNFaces = labelList(nRegions, 0);
237  forAll(newSetFaces, newSetFacei)
238  {
239  regionNFaces[newSetFaceRegions[newSetFacei]] ++;
240  }
242  Pstream::listCombineScatter(regionNFaces);
243  Info<< " Found " << nRegions << " contiguous regions with "
244  << regionNFaces << " faces" << endl;
245  }
246 
247  // Step 4: Choose the closest region to output
248  label selectedRegioni = -1;
249  {
250  // Compute the region centres
251  scalarField regionMagAreas(nRegions, 0);
252  pointField regionCentres(nRegions, Zero);
253  forAll(newSetFaces, newSetFacei)
254  {
255  const label facei = newSetFaces[newSetFacei];
256  const label regioni = newSetFaceRegions[newSetFacei];
257 
258  const vector& a = mesh_.faceAreas()[facei];
259  const point& c = mesh_.faceCentres()[facei];
260 
261  regionMagAreas[regioni] += mag(a);
262  regionCentres[regioni] += mag(a)*c;
263  }
264  Pstream::listCombineGather(regionMagAreas, plusEqOp<scalar>());
266  Pstream::listCombineScatter(regionMagAreas);
267  Pstream::listCombineScatter(regionCentres);
268  regionCentres /= regionMagAreas;
269 
270  // Find the region centroid closest to the reference point
271  selectedRegioni =
273  (
274  findMin(mag(regionCentres - point_)()),
275  minOp<label>()
276  );
277 
278  // Report the selection
279  Info<< " Selecting region " << selectedRegioni << " with "
280  << regionNFaces[selectedRegioni]
281  << " faces as the closest to point " << point_ << endl;
282  }
283 
284  // Step 5: Remove any faces from the set list not in the selected region
285  {
286  // Remove faces from the list by shuffling up and resizing
287  label newSetFacei0 = 0;
288  forAll(newSetFaces, newSetFacei)
289  {
290  newSetFaces[newSetFacei0] = newSetFaces[newSetFacei];
291 
292  if (newSetFaceRegions[newSetFacei] == selectedRegioni)
293  {
294  ++ newSetFacei0;
295  }
296  }
297  newSetFaces.resize(newSetFacei0);
298  }
299  }
300 
301  // Modify the face zone set
302  DynamicList<label> newAddressing;
303  DynamicList<bool> newFlipMap;
304  if (add)
305  {
306  // Start from copy
307  newAddressing = DynamicList<label>(fzSet.addressing());
308  newFlipMap = DynamicList<bool>(fzSet.flipMap());
309 
310  // Add anything from the new set that is not already in the zone set
311  forAll(newSetFaces, newSetFacei)
312  {
313  const label facei = newSetFaces[newSetFacei];
314 
315  if (!fzSet.found(facei))
316  {
317  newAddressing.append(facei);
318  newFlipMap.append(cellIsAbovePlane[mesh_.faceOwner()[facei]]);
319  }
320  }
321  }
322  else
323  {
324  // Start from empty
325  newAddressing = DynamicList<label>(fzSet.addressing().size());
326  newFlipMap = DynamicList<bool>(fzSet.flipMap().size());
327 
328  // Add everything from the zone set that is not also in the new set
329  labelHashSet newSet(newSetFaces);
330  forAll(fzSet.addressing(), i)
331  {
332  const label facei = fzSet.addressing()[i];
333 
334  if (!newSet.found(facei))
335  {
336  newAddressing.append(facei);
337  newFlipMap.append(cellIsAbovePlane[mesh_.faceOwner()[facei]]);
338  }
339  }
340  }
341  fzSet.addressing().transfer(newAddressing);
342  fzSet.flipMap().transfer(newFlipMap);
343  fzSet.updateSet();
344 }
345 
346 
347 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
348 
350 (
351  const polyMesh& mesh,
352  const dictionary& dict
353 )
354 :
355  topoSetSource(mesh),
356  point_(dict.lookup<vector>("point")),
357  normal_(dict.lookup<vector>("normal")),
358  include_
359  (
360  includeNames_
361  [
362  dict.lookupOrDefault<word>
363  (
364  "include",
365  includeNames_[include::all]
366  )
367  ]
368  )
369 {}
370 
371 
373 (
374  const polyMesh& mesh,
375  Istream& is
376 )
377 :
378  topoSetSource(mesh),
379  point_(checkIs(is)),
380  normal_(checkIs(is)),
381  include_(includeNames_[word(checkIs(is))])
382 {}
383 
384 
385 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
386 
388 {}
389 
390 
391 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
392 
394 (
395  const topoSetSource::setAction action,
396  topoSet& set
397 ) const
398 {
399  if (!isA<faceZoneSet>(set))
400  {
402  << "Operation only allowed on a faceZoneSet." << endl;
403  }
404  else
405  {
406  faceZoneSet& fzSet = refCast<faceZoneSet>(set);
407 
408  if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
409  {
410  Info<< " Adding faces which form a plane at " << point_
411  << " with normal " << normal_ << endl;
412 
413  combine(fzSet, true);
414  }
415  else if (action == topoSetSource::DELETE)
416  {
417  Info<< " Removing faces which form a plane at " << point_
418  << " with normal " << normal_ << endl;
419 
420  combine(fzSet, false);
421  }
422  }
423 }
424 
425 
426 // ************************************************************************* //
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
const boolList & flipMap() const
Definition: faceZoneSet.H:118
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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 list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
const polyBoundaryMesh & boundaryMesh() const
Return boundaryMesh reference.
Definition: polyPatch.C:291
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
label findMin(const ListType &, const label start=0)
Find index of min element (and less than given element).
void resize(const label)
Alias for setSize(const label)
Definition: ListI.H:138
labelList findIndices(const ListType &, typename ListType::const_reference, const label start=0)
Find all occurrences of given element. Linear search.
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:51
AccessType combine(const List< T > &, AccessOp aop=accessOp< T >())
Combines sublists into one big list.
Definition: ListListOps.C:34
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
const dimensionedScalar c
Speed of light in a vacuum.
Macros for easy insertion into run-time selection tables.
static label markZones(const PrimitivePatch< FaceList, PointField > &, const BoolListType &borderEdge, labelList &faceZone)
Size and fills faceZone with zone of face.
Base class of a source for a topoSet.
Definition: topoSetSource.H:63
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
static const NamedEnum< include, 2 > includeNames_
Included region names.
A list of faces which address into the list of points.
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:319
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:113
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronise values on all mesh edges.
A class for handling words, derived from string.
Definition: word.H:59
setAction
Enumeration defining the valid actions.
Definition: topoSetSource.H:82
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
virtual ~planeToFaceZone()
Destructor.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
void add(FieldField< Field1, typename typeOfSum< Type1, Type2 >::type > &f, const FieldField< Field1, Type1 > &f1, const FieldField< Field2, Type2 > &f2)
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void updateSet()
Sort addressing and make faceSet part consistent with addressing.
Definition: faceZoneSet.C:51
Like faceSet but -reads data from faceZone -updates faceZone when writing.
Definition: faceZoneSet.H:49
General set of labels of mesh quantity (points, cells, faces).
Definition: topoSet.H:61
static void syncFaceList(const polyMesh &mesh, UList< T > &l, const CombineOp &cop)
Synchronise values on all mesh faces.
Definition: syncTools.H:381
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
label patchi
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
#define WarningInFunction
Report a warning using Foam::Warning.
Class with constructor to add usage string to table.
A topoSetSource to select faces based on the adjacent cell centres spanning a given plane...
A List with indirect addressing.
Definition: fvMatrix.H:106
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:309
messageStream Info
planeToFaceZone(const polyMesh &mesh, const dictionary &dict)
Construct from dictionary.
dimensioned< scalar > mag(const dimensioned< Type > &)
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
virtual void applyToSet(const topoSetSource::setAction action, topoSet &) const
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
T & last()
Return the last element of the list.
Definition: UListI.H:128
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
const labelList & addressing() const
Definition: faceZoneSet.H:107
A List with indirect addressing.
Definition: IndirectList.H:101
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:844