mergeOrSplitBaffles.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) 2011-2019 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 Application
25  mergeOrSplitBaffles
26 
27 Description
28  Detects faces that share points (baffles). Either merge them or
29  duplicate the points.
30 
31  Notes:
32  - can only handle pairwise boundary faces. So three faces using
33  the same points is not handled (is illegal mesh anyway)
34 
35  - there is no option to only split/merge some baffles.
36 
37  - surfaces consisting of duplicate faces can be topologically split
38  if the points on the interior of the surface cannot walk to all the
39  cells that use them in one go.
40 
41  - Parallel operation (where duplicate face is perpendicular to a coupled
42  boundary) is supported but not really tested.
43  (Note that coupled faces themselves are not seen as duplicate faces)
44 
45 \*---------------------------------------------------------------------------*/
46 
47 #include "argList.H"
48 #include "Time.H"
49 #include "syncTools.H"
50 #include "faceSet.H"
51 #include "pointSet.H"
52 #include "meshTools.H"
53 #include "polyTopoChange.H"
54 #include "polyRemoveFace.H"
55 #include "polyModifyFace.H"
56 #include "indirectPrimitivePatch.H"
57 #include "processorPolyPatch.H"
58 #include "localPointRegion.H"
59 #include "duplicatePoints.H"
60 #include "ReadFields.H"
61 #include "volFields.H"
62 #include "surfaceFields.H"
63 #include "pointFields.H"
64 
65 using namespace Foam;
66 
67 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
68 
69 void insertDuplicateMerge
70 (
71  const polyMesh& mesh,
72  const labelList& duplicates,
73  polyTopoChange& meshMod
74 )
75 {
76  const faceList& faces = mesh.faces();
77  const labelList& faceOwner = mesh.faceOwner();
78  const faceZoneMesh& faceZones = mesh.faceZones();
79 
80  forAll(duplicates, bFacei)
81  {
82  label otherFacei = duplicates[bFacei];
83 
84  if (otherFacei != -1 && otherFacei > bFacei)
85  {
86  // Two duplicate faces. Merge.
87 
88  label face0 = mesh.nInternalFaces() + bFacei;
89  label face1 = mesh.nInternalFaces() + otherFacei;
90 
91  label own0 = faceOwner[face0];
92  label own1 = faceOwner[face1];
93 
94  if (own0 < own1)
95  {
96  // Use face0 as the new internal face.
97  label zoneID = faceZones.whichZone(face0);
98  bool zoneFlip = false;
99 
100  if (zoneID >= 0)
101  {
102  const faceZone& fZone = faceZones[zoneID];
103  zoneFlip = fZone.flipMap()[fZone.whichFace(face0)];
104  }
105 
106  meshMod.setAction(polyRemoveFace(face1));
107  meshMod.setAction
108  (
110  (
111  faces[face0], // modified face
112  face0, // label of face being modified
113  own0, // owner
114  own1, // neighbour
115  false, // face flip
116  -1, // patch for face
117  false, // remove from zone
118  zoneID, // zone for face
119  zoneFlip // face flip in zone
120  )
121  );
122  }
123  else
124  {
125  // Use face1 as the new internal face.
126  label zoneID = faceZones.whichZone(face1);
127  bool zoneFlip = false;
128 
129  if (zoneID >= 0)
130  {
131  const faceZone& fZone = faceZones[zoneID];
132  zoneFlip = fZone.flipMap()[fZone.whichFace(face1)];
133  }
134 
135  meshMod.setAction(polyRemoveFace(face0));
136  meshMod.setAction
137  (
139  (
140  faces[face1], // modified face
141  face1, // label of face being modified
142  own1, // owner
143  own0, // neighbour
144  false, // face flip
145  -1, // patch for face
146  false, // remove from zone
147  zoneID, // zone for face
148  zoneFlip // face flip in zone
149  )
150  );
151  }
152  }
153  }
154 }
155 
156 
157 labelList findBaffles(const polyMesh& mesh, const labelList& boundaryFaces)
158 {
159  // Get all duplicate face labels (in boundaryFaces indices!).
161  (
162  mesh,
163  boundaryFaces
164  );
165 
166 
167  // Check that none are on processor patches
168  const polyBoundaryMesh& patches = mesh.boundaryMesh();
169 
170  forAll(duplicates, bFacei)
171  {
172  if (duplicates[bFacei] != -1)
173  {
174  label facei = mesh.nInternalFaces() + bFacei;
175  label patchi = patches.whichPatch(facei);
176 
177  if (isA<processorPolyPatch>(patches[patchi]))
178  {
180  << "Duplicate face " << facei
181  << " is on a processorPolyPatch."
182  << "This is not allowed." << nl
183  << "Face:" << facei
184  << " is on patch:" << patches[patchi].name()
185  << abort(FatalError);
186  }
187  }
188  }
189 
190 
191  // Write to faceSet for ease of postprocessing.
192  {
193  faceSet duplicateSet
194  (
195  mesh,
196  "duplicateFaces",
197  (mesh.nFaces() - mesh.nInternalFaces())/256
198  );
199 
200  forAll(duplicates, bFacei)
201  {
202  label otherFacei = duplicates[bFacei];
203 
204  if (otherFacei != -1 && otherFacei > bFacei)
205  {
206  duplicateSet.insert(mesh.nInternalFaces() + bFacei);
207  duplicateSet.insert(mesh.nInternalFaces() + otherFacei);
208  }
209  }
210 
211  Pout<< "Writing " << duplicateSet.size()
212  << " duplicate faces to faceSet " << duplicateSet.objectPath()
213  << nl << endl;
214  duplicateSet.write();
215  }
216 
217  return duplicates;
218 }
219 
220 
221 
222 
223 int main(int argc, char *argv[])
224 {
226  (
227  "Detect faces that share points (baffles).\n"
228  "Merge them or duplicate the points."
229  );
230 
231  #include "addOverwriteOption.H"
232  #include "addRegionOption.H"
234  (
235  "detectOnly",
236  "find baffles only, but do not merge or split them"
237  );
239  (
240  "split",
241  "topologically split duplicate surfaces"
242  );
244  (
245  "fields",
246  "update fields"
247  );
248 
249  #include "setRootCase.H"
250  #include "createTime.H"
252  #include "createNamedMesh.H"
253 
254  const word oldInstance = mesh.pointsInstance();
255 
256  const bool split = args.optionFound("split");
257  const bool overwrite = args.optionFound("overwrite");
258  const bool detectOnly = args.optionFound("detectOnly");
259  const bool fields = args.optionFound("fields");
260 
261  // Collect all boundary faces
262  labelList boundaryFaces(mesh.nFaces() - mesh.nInternalFaces());
263 
264  forAll(boundaryFaces, i)
265  {
266  boundaryFaces[i] = i+mesh.nInternalFaces();
267  }
268 
269 
270  if (detectOnly)
271  {
272  findBaffles(mesh, boundaryFaces);
273  return 0;
274  }
275 
276 
277  // Read objects in time directory
279 
280  if (fields) Info<< "Reading geometric fields" << nl << endl;
281 
282  #include "readVolFields.H"
283  #include "readSurfaceFields.H"
284  #include "readPointFields.H"
285 
286  Info<< endl;
287 
288 
289  // Mesh change engine
290  polyTopoChange meshMod(mesh);
291 
292 
293  if (split)
294  {
295  Pout<< "Topologically splitting duplicate surfaces"
296  << ", i.e. duplicating points internal to duplicate surfaces."
297  << nl << endl;
298 
299  // Analyse which points need to be duplicated
301 
302  // Point duplication engine
303  duplicatePoints pointDuplicator(mesh);
304 
305  // Insert topo changes
306  pointDuplicator.setRefinement(regionSide, meshMod);
307  }
308  else
309  {
310  Pout<< "Merging duplicate faces."
311  << nl << endl;
312 
313  // Get all duplicate face labels (in boundaryFaces indices!).
314  labelList duplicates(findBaffles(mesh, boundaryFaces));
315 
316  // Merge into internal faces.
317  insertDuplicateMerge(mesh, duplicates, meshMod);
318  }
319 
320  if (!overwrite)
321  {
322  runTime++;
323  }
324 
325  // Change the mesh. No inflation.
326  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
327 
328  // Update fields
329  mesh.updateMesh(map);
330 
331  // Move mesh (since morphing does not do this)
332  if (map().hasMotionPoints())
333  {
334  mesh.movePoints(map().preMotionPoints());
335  }
336 
337  if (overwrite)
338  {
339  mesh.setInstance(oldInstance);
340  }
341  Pout<< "Writing mesh to time " << runTime.timeName() << endl;
342  mesh.write();
343 
344  // Dump duplicated points (if any)
345  if (split)
346  {
347  const labelList& pointMap = map().pointMap();
348 
349  labelList nDupPerPoint(map().nOldPoints(), 0);
350 
351  pointSet dupPoints(mesh, "duplicatedPoints", 100);
352 
353  forAll(pointMap, pointi)
354  {
355  label oldPointi = pointMap[pointi];
356 
357  nDupPerPoint[oldPointi]++;
358 
359  if (nDupPerPoint[oldPointi] > 1)
360  {
361  dupPoints.insert(map().reversePointMap()[oldPointi]);
362  dupPoints.insert(pointi);
363  }
364  }
365 
366  Pout<< "Writing " << dupPoints.size()
367  << " duplicated points to pointSet "
368  << dupPoints.objectPath() << nl << endl;
369 
370  dupPoints.write();
371  }
372 
373  Info<< "End\n" << endl;
374 
375  return 0;
376 }
377 
378 
379 // ************************************************************************* //
Foam::surfaceFields.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
Class containing data for face removal.
virtual tmp< scalarField > movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
Definition: polyMesh.C:1287
#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
const word & name() const
Return name.
Definition: IOobject.H:295
A list of face labels.
Definition: faceSet.H:48
List of IOobjects with searching and retrieving facilities.
Definition: IOobjectList.H:50
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:476
Class describing modification of a face.
void off()
Switch the function objects off.
A set of point labels.
Definition: pointSet.H:48
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
label nInternalFaces() const
label nFaces() const
const boolList & flipMap() const
Return face flip map.
Definition: faceZone.H:248
engineTime & runTime
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
bool optionFound(const word &opt) const
Return true if the named option is found.
Definition: argListI.H:108
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:306
Field reading functions for post-processing utilities.
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:626
Takes mesh with &#39;baffles&#39; (= boundary faces sharing points). Determines for selected points on bounda...
Determines the &#39;side&#39; for every face and connected to a singly-connected (through edges) region of fa...
Definition: regionSide.H:61
const fileName & pointsInstance() const
Return the current instance directory for points.
Definition: polyMesh.C:802
Duplicate points.
virtual void updateMesh(const mapPolyMesh &mpm)
Update the mesh corresponding to given map.
A class for handling words, derived from string.
Definition: word.H:59
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1169
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1156
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Foam::polyBoundaryMesh.
static const char nl
Definition: Ostream.H:265
objects
void setInstance(const fileName &)
Set the instance for mesh files.
Definition: polyMeshIO.C:32
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: ZoneMesh.C:221
label patchi
const functionObjectList & functionObjects() const
Return the list of function objects.
Definition: Time.H:410
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Direct mesh changes based on v1.3 polyTopoChange syntax.
static labelList findDuplicateFaces(const primitiveMesh &, const labelList &)
Helper routine to find baffles (two boundary faces using the.
messageStream Info
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
virtual Ostream & write(const token &)=0
Write next token to stream.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
static void addBoolOption(const word &opt, const string &usage="")
Add to a bool option to validOptions with usage information.
Definition: argList.C:117
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:64
static void addNote(const string &)
Add extra notes for the usage information.
Definition: argList.C:158
autoPtr< mapPolyMesh > changeMesh(polyMesh &mesh, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
virtual bool write(const bool write=true) const
Write using setting from DB.
Foam::argList args(argc, argv)
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Namespace for OpenFOAM.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.