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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2017-2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 stitchMesh
29
30Group
31 grpMeshManipulationUtilities
32
33Description
34 'Stitches' a mesh.
35
36 Takes a mesh and two patches and merges the faces on the two patches
37 (if geometrically possible) so the faces become internal.
38
39 Can do
40 - 'perfect' match: faces and points on patches align exactly. Order might
41 be different though.
42 - 'integral' match: where the surfaces on both patches exactly
43 match but the individual faces not
44 - 'partial' match: where the non-overlapping part of the surface remains
45 in the respective patch.
46
47 Note : Is just a front-end to perfectInterface/slidingInterface.
48
49 Comparable to running a meshModifier of the form
50 (if masterPatch is called "M" and slavePatch "S"):
51
52 \verbatim
53 couple
54 {
55 type slidingInterface;
56 masterFaceZoneName MSMasterZone
57 slaveFaceZoneName MSSlaveZone
58 cutPointZoneName MSCutPointZone
59 cutFaceZoneName MSCutFaceZone
60 masterPatchName M;
61 slavePatchName S;
62 typeOfMatch partial or integral
63 }
64 \endverbatim
65
66
67\*---------------------------------------------------------------------------*/
68
69#include "fvCFD.H"
70#include "polyTopoChanger.H"
71#include "mapPolyMesh.H"
72#include "slidingInterface.H"
73#include "perfectInterface.H"
74#include "IOobjectList.H"
75#include "ReadFields.H"
76#include <numeric>
77
78// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
79
80// Checks whether patch present and non-zero
81bool checkPatch(const polyBoundaryMesh& bMesh, const word& name)
82{
83 const label patchi = bMesh.findPatchID(name);
84
85 if (patchi == -1)
86 {
87 Info<< "No patch " << name << " in mesh" << nl
88 << "Known patches: " << bMesh.names() << endl;
89
90 return false;
91 }
92
93 if (bMesh[patchi].empty())
94 {
95 Info<< "Patch " << name << " has zero size" << nl;
96
97 return false;
98 }
99
100 return true;
101}
102
103
104int main(int argc, char *argv[])
105{
106 argList::addNote
107 (
108 "Merge the faces on specified patches (if geometrically possible)"
109 " so that the\n"
110 "faces become internal.\n"
111 "This utility can be called without arguments (uses stitchMeshDict)"
112 " or with\n"
113 "two arguments (master/slave patch names)."
114 );
115
116 argList::noParallel();
117 argList::noFunctionObjects(); // Never use function objects
118
119 #include "addOverwriteOption.H"
120 #include "addRegionOption.H"
121
122 argList::addOption("dict", "file", "Alternative stitchMeshDict");
123
124 argList::addBoolOption
125 (
126 "integral",
127 "Couple integral master/slave patches (2 argument mode: default)"
128 );
129 argList::addBoolOption
130 (
131 "partial",
132 "Couple partially overlapping master/slave patches (2 argument mode)"
133 );
134 argList::addBoolOption
135 (
136 "perfect",
137 "Couple perfectly aligned master/slave patches (2 argument mode)"
138 );
139 argList::addBoolOption
140 (
141 "intermediate",
142 "Write intermediate stages, not just the final result"
143 );
144 argList::addOption
145 (
146 "toleranceDict",
147 "file",
148 "Dictionary file with tolerances"
149 );
150
151 // The arguments are optional (non-mandatory) when using dictionary mode
152 argList::noMandatoryArgs();
153 argList::addArgument
154 (
155 "master",
156 "The master patch name (non-dictionary mode)"
157 );
158 argList::addArgument
159 (
160 "slave",
161 "The slave patch name (non-dictionary mode)"
162 );
163
164 #include "setRootCase.H"
165
166 // We now handle checking args and general sanity etc.
167 const bool useCommandArgs = (args.size() > 1);
168
169 if (useCommandArgs)
170 {
171 if (args.found("dict"))
172 {
174 << "Cannot specify both dictionary and command-line arguments"
175 << nl
176 << endl;
177 }
178
179 // If we have arguments - we require all arguments!
180 if (!args.check(true, false))
181 {
183 }
184 }
185 else
186 {
187 // Carp about inapplicable options
188
189 if (args.found("integral"))
190 {
192 << "Only specify -integral with command-line arguments"
193 << endl;
194 }
195
196 if (args.found("partial"))
197 {
199 << "Only specify -partial with command-line arguments"
200 << endl;
201 }
202
203 if (args.found("perfect"))
204 {
206 << "Only specify -perfect with command-line arguments"
207 << endl;
208 }
209 }
210
211 #include "createTime.H"
212 #include "createNamedMesh.H"
213
214 const word oldInstance = mesh.pointsInstance();
215
216 const bool intermediate = args.found("intermediate");
217 const bool overwrite = args.found("overwrite");
218
219 const word dictName("stitchMeshDict");
220
221 // A validated input dictionary
222 dictionary validatedDict;
223
224 if (useCommandArgs)
225 {
226 // Command argument driven:
227 const int integralCover = args.found("integral");
228 const int partialCover = args.found("partial");
229 const int perfectCover = args.found("perfect");
230
231 if ((integralCover + partialCover + perfectCover) > 1)
232 {
234 << "Can only specify one of -integral | -partial | -perfect."
235 << nl
236 << "Use perfect match option if the patches perfectly align"
237 << " (both vertex positions and face centres)" << endl
238 << exit(FatalError);
239 }
240
241 // Patch names
242 const word masterPatchName(args[1]);
243 const word slavePatchName(args[2]);
244
245 // Patch names
246 Info<< " " << masterPatchName
247 << " / " << slavePatchName << nl;
248
249 // Bail out if either patch has problems
250 if
251 (
252 !checkPatch(mesh.boundaryMesh(), masterPatchName)
253 || !checkPatch(mesh.boundaryMesh(), slavePatchName)
254 )
255 {
257 << "Cannot continue"
258 << exit(FatalError);
259
260 return 1;
261 }
262
263 // Input was validated
264 dictionary dict;
265
266 if (perfectCover)
267 {
268 dict.add("match", word("perfect"));
269 }
270 else if (partialCover)
271 {
272 dict.add
273 (
274 "match",
275 slidingInterface::typeOfMatchNames[slidingInterface::PARTIAL]
276 );
277 }
278 else
279 {
280 dict.add
281 (
282 "match",
283 slidingInterface::typeOfMatchNames[slidingInterface::INTEGRAL]
284 );
285 }
286
287 // Patch names
288 dict.add("master", masterPatchName);
289 dict.add("slave", slavePatchName);
290
291 validatedDict.add("stitchMesh", dict);
292 }
293 else
294 {
295 // dictionary-driven:
296
298
299 Info<< "Reading " << dictIO.name() << flush;
300
301 IOdictionary stitchDict(dictIO);
302
303 Info<< " with " << stitchDict.size() << " entries" << nl;
304
305 // Suppress duplicate names
306 wordHashSet requestedPatches;
307
308 for (const entry& dEntry : stitchDict)
309 {
310 if (!dEntry.isDict())
311 {
312 Info<< "Ignoring non-dictionary entry: "
313 << dEntry.keyword() << nl;
314 continue;
315 }
316
317 const keyType& key = dEntry.keyword();
318 const dictionary& dict = dEntry.dict();
319
320 // Match type
321 word matchName;
322 if (dict.readIfPresent("match", matchName))
323 {
324 if
325 (
326 matchName != "perfect"
327 && !slidingInterface::typeOfMatchNames.found(matchName)
328 )
329 {
330 Info<< "Error: unknown match type - " << matchName
331 << " should be one of "
332 << slidingInterface::typeOfMatchNames.toc() << nl;
333 continue;
334 }
335 }
336
337 // Patch names
338 const word masterPatchName(dict.get<word>("master"));
339 const word slavePatchName(dict.get<word>("slave"));
340
341 // Patch names
342 Info<< " " << masterPatchName
343 << " / " << slavePatchName << nl;
344
345 if (!requestedPatches.insert(masterPatchName))
346 {
347 Info<< "Error: patch specified multiple times - "
348 << masterPatchName << nl;
349 continue;
350 }
351
352 if (!requestedPatches.insert(slavePatchName))
353 {
354 Info<< "Error: patch specified multiple times - "
355 << slavePatchName << nl;
356
357 requestedPatches.erase(masterPatchName);
358 continue;
359 }
360
361 // Bail out if either patch has problems
362 if
363 (
364 !checkPatch(mesh.boundaryMesh(), masterPatchName)
365 || !checkPatch(mesh.boundaryMesh(), slavePatchName)
366 )
367 {
368 requestedPatches.erase(masterPatchName);
369 requestedPatches.erase(slavePatchName);
370 continue;
371 }
372
373 // Input was validated
374
375 validatedDict.add(key, dict);
376 }
377 }
378
379 const label nActions = validatedDict.size();
380
381 Info<< nl << nActions << " validated actions" << endl;
382
383 if (!nActions)
384 {
385 Info<<"\nStopping" << nl << endl;
386 return 1;
387 }
388
389
390 // ------------------------------------------
391 // This is where the real work begins
392
393 // set up the tolerances for the sliding mesh
394 dictionary slidingTolerances;
395 if (args.found("toleranceDict"))
396 {
397 IOdictionary toleranceFile
398 (
399 IOobject
400 (
401 args["toleranceDict"],
402 runTime.constant(),
403 mesh,
404 IOobject::MUST_READ_IF_MODIFIED,
405 IOobject::NO_WRITE
406 )
407 );
408 slidingTolerances += toleranceFile;
409 }
410
411
412 // Search for list of objects for this time
413 IOobjectList objects(mesh, runTime.timeName());
414
415 // Read all current fvFields so they will get mapped
416 Info<< "Reading all current volfields" << endl;
417 PtrList<volScalarField> volScalarFields;
418 ReadFields(mesh, objects, volScalarFields);
419
420 PtrList<volVectorField> volVectorFields;
421 ReadFields(mesh, objects, volVectorFields);
422
423 PtrList<volSphericalTensorField> volSphericalTensorFields;
424 ReadFields(mesh, objects, volSphericalTensorFields);
425
426 PtrList<volSymmTensorField> volSymmTensorFields;
427 ReadFields(mesh, objects, volSymmTensorFields);
428
429 PtrList<volTensorField> volTensorFields;
430 ReadFields(mesh, objects, volTensorFields);
431
432 //- Uncomment if you want to interpolate surface fields (usually a bad idea)
433 //Info<< "Reading all current surfaceFields" << endl;
434 //PtrList<surfaceScalarField> surfaceScalarFields;
435 //ReadFields(mesh, objects, surfaceScalarFields);
436 //
437 //PtrList<surfaceVectorField> surfaceVectorFields;
438 //ReadFields(mesh, objects, surfaceVectorFields);
439 //
440 //PtrList<surfaceTensorField> surfaceTensorFields;
441 //ReadFields(mesh, objects, surfaceTensorFields);
442
443 // Increase precision for output mesh points
444 IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
445
446 polyTopoChanger stitcher(mesh, IOobject::NO_READ);
447
448 // Step through the topology changes
449 label actioni = 0;
450 for (const entry& dEntry : validatedDict)
451 {
452 const dictionary& dict = dEntry.dict();
453
454 // Match type
455 bool perfect = false;
456 slidingInterface::typeOfMatch matchType = slidingInterface::PARTIAL;
457
458 word matchName;
459 if (dict.readIfPresent("match", matchName))
460 {
461 if (matchName == "perfect")
462 {
463 perfect = true;
464 }
465 else
466 {
467 matchType = slidingInterface::typeOfMatchNames[matchName];
468 }
469 }
470
471 // Patch names
472 const word masterPatchName(dict.get<word>("master"));
473 const word slavePatchName(dict.get<word>("slave"));
474
475 // Zone names
476 const word mergePatchName(masterPatchName + slavePatchName);
477 const word cutZoneName(mergePatchName + "CutFaceZone");
478
479 Info<< nl << "========================================" << nl;
480
481 // Information messages
482 if (perfect)
483 {
484 Info<< "Coupling PERFECTLY aligned patches "
485 << masterPatchName << " / " << slavePatchName << nl << nl
486 << "Resulting (internal) faces in faceZone "
487 << cutZoneName << nl << nl
488 << "The patch vertices and face centres must align within a"
489 << " tolerance relative to the minimum edge length on the patch"
490 << nl << endl;
491 }
492 else if (matchType == slidingInterface::INTEGRAL)
493 {
494 Info<< "Coupling INTEGRALLY matching of patches "
495 << masterPatchName << " / " << slavePatchName << nl << nl
496 << "Resulting (internal) faces in faceZone "
497 << cutZoneName << nl << nl
498 << "The overall area covered by both patches should be"
499 << " identical!" << endl
500 << "If this is not the case use partial"
501 << nl << endl;
502 }
503 else
504 {
505 Info<< "Coupling PARTIALLY overlapping patches "
506 << masterPatchName << " / " << slavePatchName << nl
507 << "Resulting internal faces in faceZone "
508 << cutZoneName << nl
509 << "Uncovered faces remain in their patch"
510 << nl << endl;
511 }
512
513
514 // Master/slave patches
515 const polyPatch& masterPatch = mesh.boundaryMesh()[masterPatchName];
516 const polyPatch& slavePatch = mesh.boundaryMesh()[slavePatchName];
517
518 mesh.pointZones().clearAddressing();
519 mesh.faceZones().clearAddressing();
520 mesh.cellZones().clearAddressing();
521
522 // Lists of master and slave faces:
523 labelList faceIds;
524
525 // Markup master face ids
526 faceIds.setSize(masterPatch.size());
527 std::iota(faceIds.begin(), faceIds.end(), masterPatch.start());
528
529 stitcher.clear();
530 stitcher.setSize(1);
531
532 if (perfect)
533 {
534 // Add new (empty) zone for resulting internal faces
535 mesh.faceZones()
536 (
537 cutZoneName,
538 true // verbose
539 ).resetAddressing(std::move(faceIds), false);
540
541
542 // Add the perfect interface mesh modifier
543 stitcher.set
544 (
545 0,
546 new perfectInterface
547 (
548 "couple" + Foam::name(actioni),
549 0,
550 stitcher,
551 cutZoneName,
552 masterPatchName,
553 slavePatchName
554 )
555 );
556 }
557 else
558 {
559 mesh.pointZones()
560 (
561 mergePatchName + "CutPointZone",
562 true // verbose
563 ) = labelList();
564
565 mesh.faceZones()
566 (
567 mergePatchName + "MasterZone",
568 true // verbose
569 ).resetAddressing(std::move(faceIds), false);
570
571 // Markup slave face ids
572 faceIds.setSize(slavePatch.size());
573 std::iota(faceIds.begin(), faceIds.end(), slavePatch.start());
574
575 mesh.faceZones()
576 (
577 mergePatchName + "SlaveZone",
578 true // verbose
579 ).resetAddressing(std::move(faceIds), false);
580
581 // Add empty zone for cut faces
582 mesh.faceZones()
583 (
584 cutZoneName,
585 true // verbose
586 ).resetAddressing(labelList(), false);
587
588
589 // Add the sliding interface mesh modifier
590 stitcher.set
591 (
592 0,
593 new slidingInterface
594 (
595 "couple" + Foam::name(actioni),
596 0,
597 stitcher,
598 mergePatchName + "MasterZone",
599 mergePatchName + "SlaveZone",
600 mergePatchName + "CutPointZone",
601 cutZoneName,
602 masterPatchName,
603 slavePatchName,
604 matchType, // integral or partial
605 true // couple/decouple mode
606 )
607 );
608
609 static_cast<slidingInterface&>(stitcher[0]).setTolerances
610 (
611 slidingTolerances,
612 true
613 );
614 }
615
616 ++actioni;
617
618 // Advance time for intermediate results or only on final
619 if (!overwrite && (intermediate || actioni == nActions))
620 {
621 ++runTime;
622 }
623
624 // Execute all polyMeshModifiers
625 autoPtr<mapPolyMesh> morphMap = stitcher.changeMesh(true);
626
627 mesh.movePoints(morphMap->preMotionPoints());
628
629 // Write mesh
630 if (overwrite)
631 {
632 mesh.setInstance(oldInstance);
633 stitcher.instance() = oldInstance;
634 }
635
636 if (intermediate || actioni == nActions)
637 {
638 Info<< nl << "Writing polyMesh to time "
639 << runTime.timeName() << endl;
640
641 // Bypass runTime write (since only writes at writeTime)
642 if
643 (
644 !runTime.objectRegistry::writeObject
645 (
646 IOstreamOption
647 (
648 runTime.writeFormat(),
649 runTime.writeCompression()
650 ),
651 true
652 )
653 )
654 {
656 << "Failed writing polyMesh."
657 << exit(FatalError);
658 }
659
660 mesh.faceZones().write();
661 mesh.pointZones().write();
662 mesh.cellZones().write();
663
664 // Write fields
665 runTime.write();
666 }
667 }
668
669 Info<< "\nEnd\n" << endl;
670
671 return 0;
672}
673
674
675// ************************************************************************* //
Field reading functions for post-processing utilities.
bool found
Y[inertIndex] max(0.0)
void setSize(const label n)
Alias for resize()
Definition: List.H:218
bool check(bool checkArgs=argList::argsMandatory(), bool checkOpts=true) const
Definition: argList.C:1913
label size() const noexcept
The number of arguments.
Definition: argListI.H:146
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
void exit(const int errNo=1)
Exit : can be called for any error to exit program.
Definition: error.C:331
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const word dictName("faMeshDefinition")
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
Definition: foamGltfBase.H:108
wordList ReadFields(const typename GeoMesh::Mesh &mesh, const IOobjectList &objects, PtrList< GeometricField< Type, PatchField, GeoMesh > > &fields, const bool syncPar=true, const bool readOldTime=false)
Read Geometric fields of templated type.
List< label > labelList
A List of labels.
Definition: List.H:66
void checkPatch(const bool allGeometry, const word &name, const PatchType &pp, pointSet &points)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
PrimitivePatch< List< face >, const pointField > bMesh
Holder of faceList and points. (v.s. e.g. primitivePatch which references points)
Definition: bMesh.H:48
HashSet< word, Hash< word > > wordHashSet
A HashSet of words, uses string hasher.
Definition: HashSet.H:82
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Ostream & flush(Ostream &os)
Flush stream.
Definition: Ostream.H:364
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
dictionary dict
IOobject dictIO
Foam::argList args(argc, argv)