stitchMesh.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  stitchMesh
26 
27 Description
28  'Stitches' a mesh.
29 
30  Takes a mesh and two patches and merges the faces on the two patches
31  (if geometrically possible) so the faces become internal.
32 
33  Can do
34  - 'perfect' match: faces and points on patches align exactly. Order might
35  be different though.
36  - 'integral' match: where the surfaces on both patches exactly
37  match but the individual faces not
38  - 'partial' match: where the non-overlapping part of the surface remains
39  in the respective patch.
40 
41  Note : Is just a front-end to perfectInterface/slidingInterface.
42 
43  Comparable to running a meshModifier of the form
44  (if masterPatch is called "M" and slavePatch "S"):
45 
46  \verbatim
47  couple
48  {
49  type slidingInterface;
50  masterFaceZoneName MSMasterZone
51  slaveFaceZoneName MSSlaveZone
52  cutPointZoneName MSCutPointZone
53  cutFaceZoneName MSCutFaceZone
54  masterPatchName M;
55  slavePatchName S;
56  typeOfMatch partial or integral
57  }
58  \endverbatim
59 
60 
61 \*---------------------------------------------------------------------------*/
62 
63 #include "argList.H"
64 #include "polyTopoChanger.H"
65 #include "mapPolyMesh.H"
66 #include "slidingInterface.H"
67 #include "perfectInterface.H"
68 #include "ReadFields.H"
69 #include "volFields.H"
70 #include "surfaceFields.H"
71 #include "pointFields.H"
72 
73 using namespace Foam;
74 
75 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
76 
77 label addPointZone(const polyMesh& mesh, const word& name)
78 {
79  label zoneID = mesh.pointZones().findZoneID(name);
80 
81  if (zoneID != -1)
82  {
83  Info<< "Reusing existing pointZone "
84  << mesh.pointZones()[zoneID].name()
85  << " at index " << zoneID << endl;
86  }
87  else
88  {
89  pointZoneMesh& pointZones = const_cast<polyMesh&>(mesh).pointZones();
90  zoneID = pointZones.size();
91  Info<< "Adding pointZone " << name << " at index " << zoneID << endl;
92 
93  pointZones.setSize(zoneID+1);
94  pointZones.set
95  (
96  zoneID,
97  new pointZone
98  (
99  name,
100  labelList(0),
101  zoneID,
102  pointZones
103  )
104  );
105  }
106  return zoneID;
107 }
108 
109 
110 label addFaceZone(const polyMesh& mesh, const word& name)
111 {
112  label zoneID = mesh.faceZones().findZoneID(name);
113 
114  if (zoneID != -1)
115  {
116  Info<< "Reusing existing faceZone " << mesh.faceZones()[zoneID].name()
117  << " at index " << zoneID << endl;
118  }
119  else
120  {
121  faceZoneMesh& faceZones = const_cast<polyMesh&>(mesh).faceZones();
122  zoneID = faceZones.size();
123  Info<< "Adding faceZone " << name << " at index " << zoneID << endl;
124 
125  faceZones.setSize(zoneID+1);
126  faceZones.set
127  (
128  zoneID,
129  new faceZone
130  (
131  name,
132  labelList(0),
133  boolList(),
134  zoneID,
135  faceZones
136  )
137  );
138  }
139  return zoneID;
140 }
141 
142 
143 label addCellZone(const polyMesh& mesh, const word& name)
144 {
145  label zoneID = mesh.cellZones().findZoneID(name);
146 
147  if (zoneID != -1)
148  {
149  Info<< "Reusing existing cellZone " << mesh.cellZones()[zoneID].name()
150  << " at index " << zoneID << endl;
151  }
152  else
153  {
154  cellZoneMesh& cellZones = const_cast<polyMesh&>(mesh).cellZones();
155  zoneID = cellZones.size();
156  Info<< "Adding cellZone " << name << " at index " << zoneID << endl;
157 
158  cellZones.setSize(zoneID+1);
159  cellZones.set
160  (
161  zoneID,
162  new cellZone
163  (
164  name,
165  labelList(0),
166  zoneID,
167  cellZones
168  )
169  );
170  }
171  return zoneID;
172 }
173 
174 
175 // Checks whether patch present
176 void checkPatch(const polyBoundaryMesh& bMesh, const word& name)
177 {
178  const label patchi = bMesh.findPatchID(name);
179 
180  if (patchi == -1)
181  {
183  << "Cannot find patch " << name << endl
184  << "It should be present and of non-zero size" << endl
185  << "Valid patches are " << bMesh.names()
186  << exit(FatalError);
187  }
188 
189  if (bMesh[patchi].empty())
190  {
192  << "Patch " << name << " is present but zero size"
193  << exit(FatalError);
194  }
195 }
196 
197 
198 
199 int main(int argc, char *argv[])
200 {
202  (
203  "Merge the faces on the specified patches (if geometrically possible)\n"
204  "so the faces become internal.\n"
205  "Integral matching is used when the options -partial and -perfect are "
206  "omitted.\n"
207  );
208 
210  #include "addOverwriteOption.H"
211  #include "addRegionOption.H"
213  (
214  "noFields",
215  "do not update fields"
216  );
217 
218  argList::validArgs.append("masterPatch");
219  argList::validArgs.append("slavePatch");
220 
222  (
223  "partial",
224  "couple partially overlapping patches (optional)"
225  );
227  (
228  "perfect",
229  "couple perfectly aligned patches (optional)"
230  );
232  (
233  "toleranceDict",
234  "file",
235  "dictionary file with tolerances"
236  );
237 
238  #include "setRootCase.H"
239  #include "createTime.H"
241  #include "createNamedMesh.H"
242 
243  const word oldInstance = mesh.pointsInstance();
244 
245  const word masterPatchName = args[1];
246  const word slavePatchName = args[2];
247 
248  const bool partialCover = args.optionFound("partial");
249  const bool perfectCover = args.optionFound("perfect");
250  const bool overwrite = args.optionFound("overwrite");
251  const bool fields = !args.optionFound("noFields");
252 
253  if (partialCover && perfectCover)
254  {
256  << "Cannot supply both partial and perfect." << endl
257  << "Use perfect match option if the patches perfectly align"
258  << " (both vertex positions and face centres)" << endl
259  << exit(FatalError);
260  }
261 
262 
263  const word mergePatchName(masterPatchName + slavePatchName);
264  const word cutZoneName(mergePatchName + "CutFaceZone");
265 
267 
268  if (partialCover)
269  {
270  Info<< "Coupling partially overlapping patches "
271  << masterPatchName << " and " << slavePatchName << nl
272  << "Resulting internal faces will be in faceZone " << cutZoneName
273  << nl
274  << "Any uncovered faces will remain in their patch"
275  << endl;
276 
278  }
279  else if (perfectCover)
280  {
281  Info<< "Coupling perfectly aligned patches "
282  << masterPatchName << " and " << slavePatchName << nl
283  << "Resulting (internal) faces will be in faceZone " << cutZoneName
284  << nl << nl
285  << "Note: both patches need to align perfectly." << nl
286  << "Both the vertex"
287  << " positions and the face centres need to align to within" << nl
288  << "a tolerance given by the minimum edge length on the patch"
289  << endl;
290  }
291  else
292  {
293  Info<< "Coupling patches " << masterPatchName << " and "
294  << slavePatchName << nl
295  << "Resulting (internal) faces will be in faceZone " << cutZoneName
296  << nl << nl
297  << "Note: the overall area covered by both patches should be"
298  << " identical (\"integral\" interface)." << endl
299  << "If this is not the case use the -partial option" << nl << endl;
300  }
301 
302  // set up the tolerances for the sliding mesh
303  dictionary slidingTolerances;
304  if (args.options().found("toleranceDict"))
305  {
306  IOdictionary toleranceFile
307  (
308  IOobject
309  (
310  args.options()["toleranceDict"],
311  runTime.constant(),
312  mesh,
315  )
316  );
317  slidingTolerances += toleranceFile;
318  }
319 
320  // Check for non-empty master and slave patches
321  checkPatch(mesh.boundaryMesh(), masterPatchName);
322  checkPatch(mesh.boundaryMesh(), slavePatchName);
323 
324  // Create and add face zones and mesh modifiers
325 
326  // Master patch
327  const polyPatch& masterPatch = mesh.boundaryMesh()[masterPatchName];
328 
329  // Make list of masterPatch faces
330  labelList isf(masterPatch.size());
331 
332  forAll(isf, i)
333  {
334  isf[i] = masterPatch.start() + i;
335  }
336 
337  polyTopoChanger stitcher(mesh);
338  stitcher.setSize(1);
339 
340  mesh.pointZones().clearAddressing();
341  mesh.faceZones().clearAddressing();
342  mesh.cellZones().clearAddressing();
343 
344  if (perfectCover)
345  {
346  // Add empty zone for resulting internal faces
347  label cutZoneID = addFaceZone(mesh, cutZoneName);
348 
349  mesh.faceZones()[cutZoneID].resetAddressing
350  (
351  isf,
352  boolList(masterPatch.size(), false)
353  );
354 
355  // Add the perfect interface mesh modifier
356  stitcher.set
357  (
358  0,
359  new perfectInterface
360  (
361  "couple",
362  0,
363  stitcher,
364  cutZoneName,
365  masterPatchName,
366  slavePatchName
367  )
368  );
369  }
370  else
371  {
372  label pointZoneID = addPointZone(mesh, mergePatchName + "CutPointZone");
373  mesh.pointZones()[pointZoneID] = labelList(0);
374 
375  label masterZoneID = addFaceZone(mesh, mergePatchName + "MasterZone");
376 
377  mesh.faceZones()[masterZoneID].resetAddressing
378  (
379  isf,
380  boolList(masterPatch.size(), false)
381  );
382 
383  // Slave patch
384  const polyPatch& slavePatch = mesh.boundaryMesh()[slavePatchName];
385 
386  labelList osf(slavePatch.size());
387 
388  forAll(osf, i)
389  {
390  osf[i] = slavePatch.start() + i;
391  }
392 
393  label slaveZoneID = addFaceZone(mesh, mergePatchName + "SlaveZone");
394  mesh.faceZones()[slaveZoneID].resetAddressing
395  (
396  osf,
397  boolList(slavePatch.size(), false)
398  );
399 
400  // Add empty zone for cut faces
401  label cutZoneID = addFaceZone(mesh, cutZoneName);
402  mesh.faceZones()[cutZoneID].resetAddressing
403  (
404  labelList(0),
405  boolList(0, false)
406  );
407 
408 
409  // Add the sliding interface mesh modifier
410  stitcher.set
411  (
412  0,
413  new slidingInterface
414  (
415  "couple",
416  0,
417  stitcher,
418  mergePatchName + "MasterZone",
419  mergePatchName + "SlaveZone",
420  mergePatchName + "CutPointZone",
421  cutZoneName,
422  masterPatchName,
423  slavePatchName,
424  tom, // integral or partial
425  true // couple/decouple mode
426  )
427  );
428  static_cast<slidingInterface&>(stitcher[0]).setTolerances
429  (
430  slidingTolerances,
431  true
432  );
433  }
434 
435 
436  // Search for list of objects for this time
438 
439  if (fields) Info<< "Reading geometric fields" << nl << endl;
440 
441  #include "readVolFields.H"
442  // #include "readSurfaceFields.H"
443  #include "readPointFields.H"
444 
445  if (!overwrite)
446  {
447  runTime++;
448  }
449 
450  // Execute all polyMeshModifiers
451  autoPtr<mapPolyMesh> morphMap = stitcher.changeMesh(true);
452 
453  mesh.movePoints(morphMap->preMotionPoints());
454 
455  // Write mesh
456  if (overwrite)
457  {
458  mesh.setInstance(oldInstance);
459  stitcher.instance() = oldInstance;
460  stitcher.writeOpt() = IOobject::NO_WRITE;
461  }
462  Info<< nl << "Writing polyMesh to time " << runTime.timeName() << endl;
463 
465 
466  // Bypass runTime write (since only writes at writeTime)
467  if
468  (
469  !runTime.objectRegistry::writeObject
470  (
474  true
475  )
476  )
477  {
479  << "Failed writing polyMesh."
480  << exit(FatalError);
481  }
482 
483  mesh.faceZones().write();
484  mesh.pointZones().write();
485  mesh.cellZones().write();
486 
487  // Write fields
488  runTime.write();
489 
490  Info<< "End\n" << endl;
491 
492  return 0;
493 }
494 
495 
496 // ************************************************************************* //
Foam::surfaceFields.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
DynamicID< pointZoneMesh > pointZoneID
Foam::pointZoneID.
Definition: ZoneIDs.H:47
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
void clearAddressing()
Clear addressing.
Definition: ZoneMesh.C:387
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
bool set(const label) const
Is element set.
Definition: PtrListI.H:65
List of IOobjects with searching and retrieving facilities.
Definition: IOobjectList.H:50
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:476
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
void off()
Switch the function objects off.
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
static unsigned int defaultPrecision()
Return the default precision.
Definition: IOstream.H:461
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 findPatchID(const word &patchName) const
Find patch index given a name.
static void noParallel()
Remove the parallel options.
Definition: argList.C:174
static SLList< string > validArgs
A list of valid (mandatory) arguments.
Definition: argList.H:153
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
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:52
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
const fileName & pointsInstance() const
Return the current instance directory for points.
Definition: polyMesh.C:802
const Foam::HashTable< string > & options() const
Return options.
Definition: argListI.H:90
dynamicFvMesh & mesh
label findZoneID(const word &zoneName) const
Find zone index given a name.
Definition: ZoneMesh.C:341
List of mesh modifiers defining the mesh dynamics.
const pointZoneMesh & pointZones() const
Return point zone mesh.
Definition: polyMesh.H:470
A class for handling words, derived from string.
Definition: word.H:59
const cellZoneMesh & cellZones() const
Return cell zone mesh.
Definition: polyMesh.H:482
wordList names() const
Return a list of patch names.
const word & constant() const
Return constant name.
Definition: TimePaths.H:124
static void addOption(const word &opt, const string &param="", const string &usage="")
Add to an option to validOptions with usage information.
Definition: argList.C:127
IOstream::compressionType writeCompression() const
Default write compression.
Definition: Time.H:300
List< label > labelList
A List of labels.
Definition: labelList.H:56
void setSize(const label)
Reset size of PtrList. If extending the PtrList, new entries are.
Definition: PtrList.C:131
Foam::polyBoundaryMesh.
static const char nl
Definition: Ostream.H:265
Sliding interface mesh modifier. Given two face zones, couple the master and slave side using a cutti...
IOstream::streamFormat writeFormat() const
Default write format.
Definition: Time.H:288
objects
A subset of mesh cells.
Definition: cellZone.H:61
typeOfMatch
Type of match.
void setInstance(const fileName &)
Set the instance for mesh files.
Definition: polyMeshIO.C:32
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
A subset of mesh points. The labels of points in the zone can be obtained from the addressing() list...
Definition: pointZone.H:62
const functionObjectList & functionObjects() const
Return the list of function objects.
Definition: Time.H:410
static const versionNumber currentVersion
Current version number.
Definition: IOstream.H:206
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:303
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
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
Hack of attachDetach to couple patches when they perfectly align. Does not decouple. Used by stitchMesh app. Does geometric matching.
virtual bool write(const bool write=true) const
Write using setting from DB.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
Foam::argList args(argc, argv)
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
const pointField & preMotionPoints() const
Pre-motion point positions.
Definition: mapPolyMesh.H:602
Namespace for OpenFOAM.