FreeFOAM The Cross-Platform CFD Toolkit
stitchMesh.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) 1991-2010 OpenCFD Ltd.
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
42  Is just a front-end to perfectInterface/slidingInterface.
43 
44  Comparable to running a meshModifier of the form
45  (if masterPatch is called "M" and slavePatch "S"):
46 
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 
59 
60 Usage
61 
62  - stitchMesh [OPTIONS] <masterPatch> <slavePatch>
63 
64  @param <masterPatch> \n
65  @todo Detailed description of argument.
66 
67  @param <slavePatch> \n
68  @todo Detailed description of argument.
69 
70  @param -perfect \n
71  The patches, vertices and face centres are perfectly aligned.
72 
73  @param -partial \n
74  The surface overlap only partialy.
75 
76  @param -overwrite \n
77  Overwrite existing data.
78 
79  @param -case <dir>\n
80  Case directory.
81 
82  @param -help \n
83  Display help message.
84 
85  @param -doc \n
86  Display Doxygen API documentation page for this application.
87 
88  @param -srcDoc \n
89  Display Doxygen source documentation page for this application.
90 
91 \*---------------------------------------------------------------------------*/
92 
93 #include <finiteVolume/fvCFD.H>
95 #include <OpenFOAM/mapPolyMesh.H>
96 #include <OpenFOAM/ListOps.H>
99 #include <OpenFOAM/IOobjectList.H>
100 #include <OpenFOAM/ReadFields.H>
101 
102 
103 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
104 
105 label addPointZone(const polyMesh& mesh, const word& name)
106 {
107  label zoneID = mesh.pointZones().findZoneID(name);
108 
109  if (zoneID != -1)
110  {
111  Info<< "Reusing existing pointZone "
112  << mesh.pointZones()[zoneID].name()
113  << " at index " << zoneID << endl;
114  }
115  else
116  {
117  pointZoneMesh& pointZones = const_cast<polyMesh&>(mesh).pointZones();
118  zoneID = pointZones.size();
119  Info<< "Adding pointZone " << name << " at index " << zoneID << endl;
120 
121  pointZones.setSize(zoneID+1);
122  pointZones.set
123  (
124  zoneID,
125  new pointZone
126  (
127  name,
128  labelList(0),
129  zoneID,
130  pointZones
131  )
132  );
133  }
134  return zoneID;
135 }
136 
137 
138 label addFaceZone(const polyMesh& mesh, const word& name)
139 {
140  label zoneID = mesh.faceZones().findZoneID(name);
141 
142  if (zoneID != -1)
143  {
144  Info<< "Reusing existing faceZone " << mesh.faceZones()[zoneID].name()
145  << " at index " << zoneID << endl;
146  }
147  else
148  {
149  faceZoneMesh& faceZones = const_cast<polyMesh&>(mesh).faceZones();
150  zoneID = faceZones.size();
151  Info<< "Adding faceZone " << name << " at index " << zoneID << endl;
152 
153  faceZones.setSize(zoneID+1);
154  faceZones.set
155  (
156  zoneID,
157  new faceZone
158  (
159  name,
160  labelList(0),
161  boolList(),
162  zoneID,
163  faceZones
164  )
165  );
166  }
167  return zoneID;
168 }
169 
170 
171 label addCellZone(const polyMesh& mesh, const word& name)
172 {
173  label zoneID = mesh.cellZones().findZoneID(name);
174 
175  if (zoneID != -1)
176  {
177  Info<< "Reusing existing cellZone " << mesh.cellZones()[zoneID].name()
178  << " at index " << zoneID << endl;
179  }
180  else
181  {
182  cellZoneMesh& cellZones = const_cast<polyMesh&>(mesh).cellZones();
183  zoneID = cellZones.size();
184  Info<< "Adding cellZone " << name << " at index " << zoneID << endl;
185 
186  cellZones.setSize(zoneID+1);
187  cellZones.set
188  (
189  zoneID,
190  new cellZone
191  (
192  name,
193  labelList(0),
194  zoneID,
195  cellZones
196  )
197  );
198  }
199  return zoneID;
200 }
201 
202 
203 // Checks whether patch present
204 void checkPatch(const polyBoundaryMesh& bMesh, const word& name)
205 {
206  label patchI = bMesh.findPatchID(name);
207 
208  if (patchI == -1)
209  {
210  FatalErrorIn("checkPatch(const polyBoundaryMesh&, const word&)")
211  << "Cannot find patch " << name << endl
212  << "It should be present and of non-zero size" << endl
213  << "Valid patches are " << bMesh.names()
214  << exit(FatalError);
215  }
216 
217  if (bMesh[patchI].empty())
218  {
219  FatalErrorIn("checkPatch(const polyBoundaryMesh&, const word&)")
220  << "Patch " << name << " is present but zero size"
221  << exit(FatalError);
222  }
223 }
224 
225 
226 // Main program:
227 
228 int main(int argc, char *argv[])
229 {
231 # include <OpenFOAM/addRegionOption.H>
232  Foam::argList::validArgs.append("masterPatch");
233  Foam::argList::validArgs.append("slavePatch");
234 
235  Foam::argList::validOptions.insert("partial", "");
236  Foam::argList::validOptions.insert("perfect", "");
237 
238  Foam::argList::validOptions.insert("overwrite", "");
239 
240  Foam::argList::validOptions.insert("toleranceDict", "file with tolerances");
241 
242 # include <OpenFOAM/setRootCase.H>
243 # include <OpenFOAM/createTime.H>
244  runTime.functionObjects().off();
245  #include <OpenFOAM/createNamedMesh.H>
246  const word oldInstance = mesh.pointsInstance();
247 
248 
249  word masterPatchName(args.additionalArgs()[0]);
250  word slavePatchName(args.additionalArgs()[1]);
251 
252  bool partialCover = args.optionFound("partial");
253  bool perfectCover = args.optionFound("perfect");
254  bool overwrite = args.optionFound("overwrite");
255 
256  if (partialCover && perfectCover)
257  {
259  << "Cannot both supply partial and perfect." << endl
260  << "Use perfect match option if the patches perfectly align"
261  << " (both vertex positions and face centres)" << endl
262  << exit(FatalError);
263  }
264 
265 
266  const word mergePatchName(masterPatchName + slavePatchName);
267  const word cutZoneName(mergePatchName + "CutFaceZone");
268 
269  slidingInterface::typeOfMatch tom = slidingInterface::INTEGRAL;
270 
271  if (partialCover)
272  {
273  Info<< "Coupling partially overlapping patches "
274  << masterPatchName << " and " << slavePatchName << nl
275  << "Resulting internal faces will be in faceZone " << cutZoneName
276  << nl
277  << "Any uncovered faces will remain in their patch"
278  << endl;
279 
280  tom = slidingInterface::PARTIAL;
281  }
282  else if (perfectCover)
283  {
284  Info<< "Coupling perfectly aligned patches "
285  << masterPatchName << " and " << slavePatchName << nl
286  << "Resulting (internal) faces will be in faceZone " << cutZoneName
287  << nl << nl
288  << "Note: both patches need to align perfectly." << nl
289  << "Both the vertex"
290  << " positions and the face centres need to align to within" << nl
291  << "a tolerance given by the minimum edge length on the patch"
292  << endl;
293  }
294  else
295  {
296  Info<< "Coupling patches " << masterPatchName << " and "
297  << slavePatchName << nl
298  << "Resulting (internal) faces will be in faceZone " << cutZoneName
299  << nl << nl
300  << "Note: the overall area covered by both patches should be"
301  << " identical (\"integral\" interface)." << endl
302  << "If this is not the case use the -partial option" << nl << endl;
303  }
304 
305  // set up the tolerances for the sliding mesh
306  dictionary slidingTolerances;
307  if (args.options().found("toleranceDict"))
308  {
309  IOdictionary toleranceFile(
310  IOobject(
311  args.options()["toleranceDict"],
312  runTime.constant(),
313  mesh,
314  IOobject::MUST_READ,
315  IOobject::NO_WRITE
316  )
317  );
318  slidingTolerances += toleranceFile;
319  }
320 
321  // Check for non-empty master and slave patches
322  checkPatch(mesh.boundaryMesh(), masterPatchName);
323  checkPatch(mesh.boundaryMesh(), slavePatchName);
324 
325  // Create and add face zones and mesh modifiers
326 
327  // Master patch
328  const polyPatch& masterPatch =
329  mesh.boundaryMesh()
330  [
331  mesh.boundaryMesh().findPatchID(masterPatchName)
332  ];
333 
334  // Make list of masterPatch faces
335  labelList isf(masterPatch.size());
336 
337  forAll (isf, i)
338  {
339  isf[i] = masterPatch.start() + i;
340  }
341 
342  polyTopoChanger stitcher(mesh);
343  stitcher.setSize(1);
344 
345  mesh.pointZones().clearAddressing();
346  mesh.faceZones().clearAddressing();
347  mesh.cellZones().clearAddressing();
348 
349  if (perfectCover)
350  {
351  // Add empty zone for resulting internal faces
352  label cutZoneID = addFaceZone(mesh, cutZoneName);
353 
354  mesh.faceZones()[cutZoneID].resetAddressing
355  (
356  isf,
357  boolList(masterPatch.size(), false)
358  );
359 
360  // Add the perfect interface mesh modifier
361  stitcher.set
362  (
363  0,
364  new perfectInterface
365  (
366  "couple",
367  0,
368  stitcher,
369  cutZoneName,
370  masterPatchName,
371  slavePatchName
372  )
373  );
374  }
375  else
376  {
377  label pointZoneID = addPointZone(mesh, mergePatchName + "CutPointZone");
378  mesh.pointZones()[pointZoneID] = labelList(0);
379 
380  label masterZoneID = addFaceZone(mesh, mergePatchName + "MasterZone");
381 
382  mesh.faceZones()[masterZoneID].resetAddressing
383  (
384  isf,
385  boolList(masterPatch.size(), false)
386  );
387 
388  // Slave patch
389  const polyPatch& slavePatch =
390  mesh.boundaryMesh()
391  [
392  mesh.boundaryMesh().findPatchID(slavePatchName)
393  ];
394 
395  labelList osf(slavePatch.size());
396 
397  forAll (osf, i)
398  {
399  osf[i] = slavePatch.start() + i;
400  }
401 
402  label slaveZoneID = addFaceZone(mesh, mergePatchName + "SlaveZone");
403  mesh.faceZones()[slaveZoneID].resetAddressing
404  (
405  osf,
406  boolList(slavePatch.size(), false)
407  );
408 
409  // Add empty zone for cut faces
410  label cutZoneID = addFaceZone(mesh, cutZoneName);
411  mesh.faceZones()[cutZoneID].resetAddressing
412  (
413  labelList(0),
414  boolList(0, false)
415  );
416 
417 
418  // Add the sliding interface mesh modifier
419  stitcher.set
420  (
421  0,
422  new slidingInterface
423  (
424  "couple",
425  0,
426  stitcher,
427  mergePatchName + "MasterZone",
428  mergePatchName + "SlaveZone",
429  mergePatchName + "CutPointZone",
430  cutZoneName,
431  masterPatchName,
432  slavePatchName,
433  tom, // integral or partial
434  true // couple/decouple mode
435  )
436  );
437  static_cast<slidingInterface&>(stitcher[0]).setTolerances
438  (
439  slidingTolerances,
440  true
441  );
442  }
443 
444 
445  // Search for list of objects for this time
446  IOobjectList objects(mesh, runTime.timeName());
447 
448  // Read all current fvFields so they will get mapped
449  Info<< "Reading all current volfields" << endl;
450  PtrList<volScalarField> volScalarFields;
451  ReadFields(mesh, objects, volScalarFields);
452 
453  PtrList<volVectorField> volVectorFields;
454  ReadFields(mesh, objects, volVectorFields);
455 
456  PtrList<volSphericalTensorField> volSphericalTensorFields;
457  ReadFields(mesh, objects, volSphericalTensorFields);
458 
459  PtrList<volSymmTensorField> volSymmTensorFields;
460  ReadFields(mesh, objects, volSymmTensorFields);
461 
462  PtrList<volTensorField> volTensorFields;
463  ReadFields(mesh, objects, volTensorFields);
464 
465  //- uncomment if you want to interpolate surface fields (usually bad idea)
466  //Info<< "Reading all current surfaceFields" << endl;
467  //PtrList<surfaceScalarField> surfaceScalarFields;
468  //ReadFields(mesh, objects, surfaceScalarFields);
469  //
470  //PtrList<surfaceVectorField> surfaceVectorFields;
471  //ReadFields(mesh, objects, surfaceVectorFields);
472  //
473  //PtrList<surfaceTensorField> surfaceTensorFields;
474  //ReadFields(mesh, objects, surfaceTensorFields);
475 
476  if (!overwrite)
477  {
478  runTime++;
479  }
480 
481  // Execute all polyMeshModifiers
482  autoPtr<mapPolyMesh> morphMap = stitcher.changeMesh(true);
483 
484  mesh.movePoints(morphMap->preMotionPoints());
485 
486  // Write mesh
487  if (overwrite)
488  {
489  mesh.setInstance(oldInstance);
490  stitcher.instance() = oldInstance;
491  }
492  Info << nl << "Writing polyMesh to time " << runTime.timeName() << endl;
493 
494  IOstream::defaultPrecision(10);
495 
496  // Bypass runTime write (since only writes at outputTime)
497  if
498  (
499  !runTime.objectRegistry::writeObject
500  (
501  runTime.writeFormat(),
502  IOstream::currentVersion,
503  runTime.writeCompression()
504  )
505  )
506  {
508  << "Failed writing polyMesh."
509  << exit(FatalError);
510  }
511 
512  mesh.faceZones().write();
513  mesh.pointZones().write();
514  mesh.cellZones().write();
515 
516  // Write fields
517  runTime.write();
518 
519  Info<< nl << "end" << endl;
520 
521  return 0;
522 }
523 
524 
525 // ************************ vim: set sw=4 sts=4 et: ************************ //