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-2023 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 {
54 }
55 
56 
57 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
58 
59 void Foam::planeToFaceZone::combine(faceZoneSet& fzSet, const bool add) const
60 {
61  // Mark all cells with centres above the plane
62  boolList cellIsAbovePlane(mesh_.nCells());
63  forAll(mesh_.cells(), celli)
64  {
65  cellIsAbovePlane[celli] =
66  ((mesh_.cellCentres()[celli] - point_) & normal_) > 0;
67  }
68 
69  // Mark all faces that sit between cells above and below the plane
70  boolList faceIsOnPlane(mesh_.nFaces());
71  forAll(mesh_.faceNeighbour(), facei)
72  {
73  faceIsOnPlane[facei] =
74  cellIsAbovePlane[mesh_.faceOwner()[facei]]
75  != cellIsAbovePlane[mesh_.faceNeighbour()[facei]];
76  }
78  {
79  const polyPatch& patch = mesh_.boundaryMesh()[patchi];
80  forAll(patch, patchFacei)
81  {
82  const label facei = patch.start() + patchFacei;
83  faceIsOnPlane[facei] =
84  patch.coupled() && cellIsAbovePlane[mesh_.faceOwner()[facei]];
85  }
86  }
87  syncTools::syncFaceList(mesh_, faceIsOnPlane, notEqOp<bool>());
88 
89  // Convert marked faces to a list of indices
90  labelList newSetFaces(findIndices(faceIsOnPlane, true));
91 
92  // If constructing a single contiguous set, remove all faces except those
93  // connected to the contiguous region closest to the specified point
94  if (include_ == include::closest)
95  {
96  // Step 1: Get locally contiguous regions for the new face set and the
97  // total number of regions across all processors.
98  labelList newSetFaceRegions(newSetFaces.size(), -1);
99  label nRegions = -1;
100  {
101  // Create a patch of the set faces
102  const uindirectPrimitivePatch newSetPatch
103  (
104  UIndirectList<face>(mesh_.faces(), newSetFaces),
105  mesh_.points()
106  );
107 
108  // Get the region ID-s and store the total number of regions on
109  // each processor
110  labelList procNRegions(Pstream::nProcs(), -1);
111  procNRegions[Pstream::myProcNo()] =
113  (
114  newSetPatch,
115  boolList(newSetPatch.nEdges(), false),
116  newSetFaceRegions
117  );
118  Pstream::gatherList(procNRegions);
119  Pstream::scatterList(procNRegions);
120 
121  // Cumulative sum the number of regions on each processor to get an
122  // offset which makes the local region ID-s globally unique
123  labelList procRegionOffset(Pstream::nProcs(), 0);
124  for (label proci = 1; proci < Pstream::nProcs(); ++ proci)
125  {
126  procRegionOffset[proci] +=
127  procRegionOffset[proci - 1]
128  + procNRegions[proci - 1];
129  }
130 
131  // Apply the offset
132  forAll(newSetFaces, newSetFacei)
133  {
134  newSetFaceRegions[newSetFacei] +=
135  procRegionOffset[Pstream::myProcNo()];
136  }
137 
138  // Store the total number of regions across all processors
139  nRegions = procRegionOffset.last() + procNRegions.last();
140  }
141 
142  // Step 2: Create a region map which combines regions which are
143  // connected across coupled interfaces
144  labelList regionMap(identityMap(nRegions));
145  {
146  // Put region labels on connected boundary edges and synchronise to
147  // create a list of all regions connected to a given edge
148  labelListList meshEdgeRegions(mesh_.nEdges(), labelList());
149  forAll(newSetFaces, newSetFacei)
150  {
151  const label facei = newSetFaces[newSetFacei];
152  const label regioni = newSetFaceRegions[newSetFacei];
153  forAll(mesh_.faceEdges()[facei], faceEdgei)
154  {
155  const label edgei = mesh_.faceEdges()[facei][faceEdgei];
156  meshEdgeRegions[edgei] = labelList(1, regioni);
157  }
158  }
160  (
161  mesh_,
162  meshEdgeRegions,
163  globalMeshData::ListPlusEqOp<labelList>(),
164  labelList()
165  );
166 
167  // Combine edge regions to create a list of what regions a given
168  // region is connected to
169  List<labelHashSet> regionRegions(nRegions);
170  forAll(newSetFaces, newSetFacei)
171  {
172  const label facei = newSetFaces[newSetFacei];
173  const label regioni = newSetFaceRegions[newSetFacei];
174  forAll(mesh_.faceEdges()[facei], faceEdgei)
175  {
176  const label edgei = mesh_.faceEdges()[facei][faceEdgei];
177  forAll(meshEdgeRegions[edgei], edgeRegioni)
178  {
179  if (meshEdgeRegions[edgei][edgeRegioni] != regioni)
180  {
181  regionRegions[regioni].insert
182  (
183  meshEdgeRegions[edgei][edgeRegioni]
184  );
185  }
186  }
187  }
188  }
189  Pstream::listCombineGather(regionRegions, plusEqOp<labelHashSet>());
190  Pstream::listCombineScatter(regionRegions);
191 
192  // Collapse the region connections into a map between each region
193  // and the lowest numbered region that it connects to
194  forAll(regionRegions, regioni)
195  {
196  forAllConstIter(labelHashSet, regionRegions[regioni], iter)
197  {
198  regionMap[iter.key()] =
199  min(regionMap[iter.key()], regionMap[regioni]);
200  }
201  }
202  }
203 
204  // Step 3: Combine connected regions
205  labelList regionNFaces;
206  {
207  // Remove duplicates from the region map
208  label regioni0 = 0;
209  forAll(regionMap, regioni)
210  {
211  if (regionMap[regioni] > regioni0)
212  {
213  ++ regioni0;
214  regionMap[regioni] = regioni0;
215  }
216  }
217 
218  // Recompute the number of regions
219  nRegions = regioni0 + 1;
220 
221  // Renumber the face region ID-s
222  newSetFaceRegions =
223  IndirectList<label>(regionMap, newSetFaceRegions);
224 
225  // Report the final number and size of the regions
226  regionNFaces = labelList(nRegions, 0);
227  forAll(newSetFaces, newSetFacei)
228  {
229  regionNFaces[newSetFaceRegions[newSetFacei]] ++;
230  }
231  Pstream::listCombineGather(regionNFaces, plusEqOp<label>());
232  Pstream::listCombineScatter(regionNFaces);
233  Info<< " Found " << nRegions << " contiguous regions with "
234  << regionNFaces << " faces" << endl;
235  }
236 
237  // Step 4: Choose the closest region to output
238  label selectedRegioni = -1;
239  {
240  // Compute the region centres
241  scalarField regionMagAreas(nRegions, 0);
242  pointField regionCentres(nRegions, Zero);
243  forAll(newSetFaces, newSetFacei)
244  {
245  const label facei = newSetFaces[newSetFacei];
246  const label regioni = newSetFaceRegions[newSetFacei];
247 
248  const vector& a = mesh_.faceAreas()[facei];
249  const point& c = mesh_.faceCentres()[facei];
250 
251  regionMagAreas[regioni] += mag(a);
252  regionCentres[regioni] += mag(a)*c;
253  }
254  Pstream::listCombineGather(regionMagAreas, plusEqOp<scalar>());
255  Pstream::listCombineGather(regionCentres, plusEqOp<point>());
256  Pstream::listCombineScatter(regionMagAreas);
257  Pstream::listCombineScatter(regionCentres);
258  regionCentres /= regionMagAreas;
259 
260  // Find the region centroid closest to the reference point
261  selectedRegioni =
263  (
264  findMin(mag(regionCentres - point_)()),
265  minOp<label>()
266  );
267 
268  // Report the selection
269  Info<< " Selecting region " << selectedRegioni << " with "
270  << regionNFaces[selectedRegioni]
271  << " faces as the closest to point " << point_ << endl;
272  }
273 
274  // Step 5: Remove any faces from the set list not in the selected region
275  {
276  // Remove faces from the list by shuffling up and resizing
277  label newSetFacei0 = 0;
278  forAll(newSetFaces, newSetFacei)
279  {
280  newSetFaces[newSetFacei0] = newSetFaces[newSetFacei];
281 
282  if (newSetFaceRegions[newSetFacei] == selectedRegioni)
283  {
284  ++ newSetFacei0;
285  }
286  }
287  newSetFaces.resize(newSetFacei0);
288  }
289  }
290 
291  // Modify the face zone set
292  DynamicList<label> newAddressing;
293  DynamicList<bool> newFlipMap;
294  if (add)
295  {
296  // Start from copy
297  newAddressing = DynamicList<label>(fzSet.addressing());
298  newFlipMap = DynamicList<bool>(fzSet.flipMap());
299 
300  // Add anything from the new set that is not already in the zone set
301  forAll(newSetFaces, newSetFacei)
302  {
303  const label facei = newSetFaces[newSetFacei];
304 
305  if (!fzSet.found(facei))
306  {
307  newAddressing.append(facei);
308  newFlipMap.append(cellIsAbovePlane[mesh_.faceOwner()[facei]]);
309  }
310  }
311  }
312  else
313  {
314  // Start from empty
315  newAddressing = DynamicList<label>(fzSet.addressing().size());
316  newFlipMap = DynamicList<bool>(fzSet.flipMap().size());
317 
318  // Add everything from the zone set that is not also in the new set
319  labelHashSet newSet(newSetFaces);
320  forAll(fzSet.addressing(), i)
321  {
322  const label facei = fzSet.addressing()[i];
323 
324  if (!newSet.found(facei))
325  {
326  newAddressing.append(facei);
327  newFlipMap.append(cellIsAbovePlane[mesh_.faceOwner()[facei]]);
328  }
329  }
330  }
331  fzSet.addressing().transfer(newAddressing);
332  fzSet.flipMap().transfer(newFlipMap);
333  fzSet.updateSet();
334 }
335 
336 
337 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
338 
340 (
341  const polyMesh& mesh,
342  const dictionary& dict
343 )
344 :
345  topoSetSource(mesh),
346  point_(dict.lookup<vector>("point")),
347  normal_(dict.lookup<vector>("normal")),
348  include_
349  (
350  includeNames_
351  [
352  dict.lookupOrDefault<word>
353  (
354  "include",
355  includeNames_[include::all]
356  )
357  ]
358  )
359 {}
360 
361 
362 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
363 
365 {}
366 
367 
368 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
369 
371 (
372  const topoSetSource::setAction action,
373  topoSet& set
374 ) const
375 {
376  if (!isA<faceZoneSet>(set))
377  {
379  << "Operation only allowed on a faceZoneSet." << endl;
380  }
381  else
382  {
383  faceZoneSet& fzSet = refCast<faceZoneSet>(set);
384 
385  if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
386  {
387  Info<< " Adding faces which form a plane at " << point_
388  << " with normal " << normal_ << endl;
389 
390  combine(fzSet, true);
391  }
392  else if (action == topoSetSource::DELETE)
393  {
394  Info<< " Removing faces which form a plane at " << point_
395  << " with normal " << normal_ << endl;
396 
397  combine(fzSet, false);
398  }
399  }
400 }
401 
402 
403 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
Macros for easy insertion into run-time selection tables.
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:54
static label markZones(const PrimitivePatch< FaceList, PointField > &, const BoolListType &borderEdge, labelList &faceZone)
Size and fills faceZone with zone of face.
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:160
Like faceSet but -reads data from faceZone -updates faceZone when writing.
Definition: faceZoneSet.H:52
A topoSetSource to select faces based on the adjacent cell centres spanning a given plane....
virtual void applyToSet(const topoSetSource::setAction action, topoSet &) const
static const NamedEnum< include, 2 > includeNames_
Included region names.
include
Enumeration for what to include.
virtual ~planeToFaceZone()
Destructor.
planeToFaceZone(const polyMesh &mesh, const dictionary &dict)
Construct from dictionary.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1373
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1386
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:405
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1392
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1360
label nEdges() const
label nCells() const
const vectorField & faceCentres() const
const vectorField & cellCentres() const
label nFaces() const
const labelListList & faceEdges() const
const vectorField & faceAreas() const
const cellList & cells() const
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronise values on all mesh edges.
static void syncFaceList(const polyMesh &mesh, UList< T > &l, const CombineOp &cop)
Synchronise values on all mesh faces.
Definition: syncTools.H:387
Base class of a source for a topoSet.
Definition: topoSetSource.H:64
setAction
Enumeration defining the valid actions.
Definition: topoSetSource.H:83
const polyMesh & mesh_
Definition: topoSetSource.H:99
General set of labels of mesh quantity (points, cells, faces).
Definition: topoSet.H:65
A class for handling words, derived from string.
Definition: word.H:62
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionedScalar c
Speed of light in a vacuum.
Namespace for OpenFOAM.
static const zero Zero
Definition: zero.H:97
List< label > labelList
A List of labels.
Definition: labelList.H:56
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
labelList findIndices(const ListType &, typename ListType::const_reference, const label start=0)
Find all occurrences of given element. Linear search.
messageStream Info
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
vector point
Point is a vector.
Definition: point.H:41
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
PrimitivePatch< UIndirectList< face >, const pointField & > uindirectPrimitivePatch
Foam::uindirectPrimitivePatch.
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
void add(FieldField< Field1, typename typeOfSum< Type1, Type2 >::type > &f, const FieldField< Field1, Type1 > &f1, const FieldField< Field2, Type2 > &f2)
dimensioned< scalar > mag(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
label findMin(const ListType &, const label start=0)
Find index of min element (and less than given element).
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:211
labelList identityMap(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
treeBoundBox combine(const treeBoundBox &a, const treeBoundBox &b)
Definition: patchToPatch.C:78
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
dictionary dict