mergeOrSplitBaffles.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-2016 OpenFOAM Foundation
9 Copyright (C) 2016-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 mergeOrSplitBaffles
29
30Group
31 grpMeshManipulationUtilities
32
33Description
34 Detects boundary faces that share points (baffles). Either merges them or
35 duplicate the points.
36
37Usage
38 \b mergeOrSplitBaffles [OPTION]
39
40 Options:
41 - \par -detectOnly
42 Detect baffles and write to faceSet duplicateFaces.
43
44 - \par -split
45 Detect baffles and duplicate the points (used so the two sides
46 can move independently)
47
48 - \par -dict <dictionary>
49 Specify a dictionary to read actions from.
50
51Note
52 - can only handle pairwise boundary faces. So three faces using
53 the same points is not handled (is illegal mesh anyway)
54
55 - surfaces consisting of duplicate faces can be topologically split
56 if the points on the interior of the surface cannot walk to all the
57 cells that use them in one go.
58
59 - Parallel operation (where duplicate face is perpendicular to a coupled
60 boundary) is supported but not really tested.
61 (Note that coupled faces themselves are not seen as duplicate faces)
62
63\*---------------------------------------------------------------------------*/
64
65#include "argList.H"
66#include "Time.H"
67#include "syncTools.H"
68#include "faceSet.H"
69#include "pointSet.H"
70#include "meshTools.H"
71#include "polyTopoChange.H"
72#include "polyRemoveFace.H"
73#include "polyModifyFace.H"
75#include "processorPolyPatch.H"
76#include "localPointRegion.H"
77#include "duplicatePoints.H"
78#include "ReadFields.H"
79#include "volFields.H"
80#include "surfaceFields.H"
81#include "processorMeshes.H"
82
83using namespace Foam;
84
85// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
86
87void insertDuplicateMerge
88(
89 const polyMesh& mesh,
90 const labelList& boundaryFaces,
91 const labelList& duplicates,
92 polyTopoChange& meshMod
93)
94{
95 const faceList& faces = mesh.faces();
96 const labelList& faceOwner = mesh.faceOwner();
97 const faceZoneMesh& faceZones = mesh.faceZones();
98
99 forAll(duplicates, bFacei)
100 {
101 label otherFacei = duplicates[bFacei];
102
103 if (otherFacei != -1 && otherFacei > bFacei)
104 {
105 // Two duplicate faces. Merge.
106
107 label face0 = boundaryFaces[bFacei];
108 label face1 = boundaryFaces[otherFacei];
109
110 label own0 = faceOwner[face0];
111 label own1 = faceOwner[face1];
112
113 if (own0 < own1)
114 {
115 // Use face0 as the new internal face.
116 label zoneID = faceZones.whichZone(face0);
117 bool zoneFlip = false;
118
119 if (zoneID >= 0)
120 {
121 const faceZone& fZone = faceZones[zoneID];
122 zoneFlip = fZone.flipMap()[fZone.whichFace(face0)];
123 }
124
125 meshMod.setAction(polyRemoveFace(face1));
126 meshMod.setAction
127 (
129 (
130 faces[face0], // modified face
131 face0, // label of face being modified
132 own0, // owner
133 own1, // neighbour
134 false, // face flip
135 -1, // patch for face
136 false, // remove from zone
137 zoneID, // zone for face
138 zoneFlip // face flip in zone
139 )
140 );
141 }
142 else
143 {
144 // Use face1 as the new internal face.
145 label zoneID = faceZones.whichZone(face1);
146 bool zoneFlip = false;
147
148 if (zoneID >= 0)
149 {
150 const faceZone& fZone = faceZones[zoneID];
151 zoneFlip = fZone.flipMap()[fZone.whichFace(face1)];
152 }
153
154 meshMod.setAction(polyRemoveFace(face0));
155 meshMod.setAction
156 (
158 (
159 faces[face1], // modified face
160 face1, // label of face being modified
161 own1, // owner
162 own0, // neighbour
163 false, // face flip
164 -1, // patch for face
165 false, // remove from zone
166 zoneID, // zone for face
167 zoneFlip // face flip in zone
168 )
169 );
170 }
171 }
172 }
173}
174
175
176label patchSize(const polyMesh& mesh, const labelList& patchIDs)
177{
178 const polyBoundaryMesh& patches = mesh.boundaryMesh();
179
180 label sz = 0;
181 forAll(patchIDs, i)
182 {
183 const polyPatch& pp = patches[patchIDs[i]];
184 sz += pp.size();
185 }
186 return sz;
187}
188
189
190labelList patchFaces(const polyMesh& mesh, const labelList& patchIDs)
191{
192 const polyBoundaryMesh& patches = mesh.boundaryMesh();
193
194 labelList faceIDs(patchSize(mesh, patchIDs));
195 label sz = 0;
196 forAll(patchIDs, i)
197 {
198 const polyPatch& pp = patches[patchIDs[i]];
199
200 forAll(pp, ppi)
201 {
202 faceIDs[sz++] = pp.start()+ppi;
203 }
204 }
205
206 if (faceIDs.size() != sz)
207 {
208 FatalErrorInFunction << exit(FatalError);
209 }
210
211 return faceIDs;
212}
213
214
215labelList findBaffles(const polyMesh& mesh, const labelList& boundaryFaces)
216{
217 // Get all duplicate face labels (in boundaryFaces indices!).
218 labelList duplicates = localPointRegion::findDuplicateFaces
219 (
220 mesh,
221 boundaryFaces
222 );
223
224
225 // Check that none are on processor patches
226 const polyBoundaryMesh& patches = mesh.boundaryMesh();
227
228 forAll(duplicates, bFacei)
229 {
230 if (duplicates[bFacei] != -1)
231 {
232 label facei = boundaryFaces[bFacei];
233 label patchi = patches.whichPatch(facei);
234
235 if (isA<processorPolyPatch>(patches[patchi]))
236 {
238 << "Duplicate face " << facei
239 << " is on a processorPolyPatch."
240 << "This is not allowed." << nl
241 << "Face:" << facei
242 << " is on patch:" << patches[patchi].name()
243 << abort(FatalError);
244 }
245 }
246 }
247
248
249 // Write to faceSet for ease of post-processing.
250 {
251 faceSet duplicateSet
252 (
253 mesh,
254 "duplicateFaces",
255 mesh.nBoundaryFaces()/256
256 );
257
258 forAll(duplicates, bFacei)
259 {
260 label otherFacei = duplicates[bFacei];
261
262 if (otherFacei != -1 && otherFacei > bFacei)
263 {
264 duplicateSet.insert(boundaryFaces[bFacei]);
265 duplicateSet.insert(boundaryFaces[otherFacei]);
266 }
267 }
268
269 Info<< "Writing " << returnReduce(duplicateSet.size(), sumOp<label>())
270 << " duplicate faces to faceSet " << duplicateSet.objectPath()
271 << nl << endl;
272 duplicateSet.write();
273 }
274
275 return duplicates;
276}
277
278
279int main(int argc, char *argv[])
280{
281 argList::addNote
282 (
283 "Detect faces that share points (baffles).\n"
284 "Merge them or duplicate the points."
285 );
286
287 #include "addOverwriteOption.H"
288 #include "addRegionOption.H"
289 argList::addOption
290 (
291 "dict",
292 "file",
293 "Specify a dictionary to read actions from"
294 );
295 argList::addBoolOption
296 (
297 "detectOnly",
298 "Find baffles only, but do not merge or split them"
299 );
300 argList::addBoolOption
301 (
302 "split",
303 "Topologically split duplicate surfaces"
304 );
305
306 argList::noFunctionObjects(); // Never use function objects
307
308 #include "setRootCase.H"
309 #include "createTime.H"
310 #include "createNamedMesh.H"
311
312 const word oldInstance = mesh.pointsInstance();
313 const polyBoundaryMesh& patches = mesh.boundaryMesh();
314
315 const bool readDict = args.found("dict");
316 const bool split = args.found("split");
317 const bool overwrite = args.found("overwrite");
318 const bool detectOnly = args.found("detectOnly");
319
320 if (readDict && (split || detectOnly))
321 {
323 << "Use of dictionary for settings not compatible with"
324 << " using command line arguments for \"split\""
325 << " or \"detectOnly\"" << exit(FatalError);
326 }
327
328
329 labelList detectPatchIDs;
330 labelList splitPatchIDs;
331 labelList mergePatchIDs;
332
333 if (readDict)
334 {
335 const word dictName;
337
338 Info<< "Reading " << dictIO.name() << nl << endl;
340
341 if (dict.found("detect"))
342 {
343 detectPatchIDs = patches.patchSet
344 (
345 dict.subDict("detect").get<wordRes>("patches")
346 ).sortedToc();
347
348 Info<< "Detecting baffles on " << detectPatchIDs.size()
349 << " patches with "
350 << returnReduce(patchSize(mesh, detectPatchIDs), sumOp<label>())
351 << " faces" << endl;
352 }
353 if (dict.found("merge"))
354 {
355 mergePatchIDs = patches.patchSet
356 (
357 dict.subDict("merge").get<wordRes>("patches")
358 ).sortedToc();
359
360 Info<< "Detecting baffles on " << mergePatchIDs.size()
361 << " patches with "
362 << returnReduce(patchSize(mesh, mergePatchIDs), sumOp<label>())
363 << " faces" << endl;
364 }
365 if (dict.found("split"))
366 {
367 splitPatchIDs = patches.patchSet
368 (
369 dict.subDict("split").get<wordRes>("patches")
370 ).sortedToc();
371
372 Info<< "Detecting baffles on " << splitPatchIDs.size()
373 << " patches with "
374 << returnReduce(patchSize(mesh, splitPatchIDs), sumOp<label>())
375 << " faces" << endl;
376 }
377 }
378 else
379 {
380 if (detectOnly)
381 {
382 detectPatchIDs = identity(patches.size());
383 }
384 else if (split)
385 {
386 splitPatchIDs = identity(patches.size());
387 }
388 else
389 {
390 mergePatchIDs = identity(patches.size());
391 }
392 }
393
394
395 if (detectPatchIDs.size())
396 {
397 findBaffles(mesh, patchFaces(mesh, detectPatchIDs));
398
399 if (detectOnly)
400 {
401 return 0;
402 }
403 }
404
405
406
407 // Read objects in time directory
408 IOobjectList objects(mesh, runTime.timeName());
409
410 // Read vol fields.
411
413 ReadFields(mesh, objects, vsFlds);
414
416 ReadFields(mesh, objects, vvFlds);
417
419 ReadFields(mesh, objects, vstFlds);
420
422 ReadFields(mesh, objects, vsymtFlds);
423
425 ReadFields(mesh, objects, vtFlds);
426
427 // Read surface fields.
428
430 ReadFields(mesh, objects, ssFlds);
431
433 ReadFields(mesh, objects, svFlds);
434
436 ReadFields(mesh, objects, sstFlds);
437
439 ReadFields(mesh, objects, ssymtFlds);
440
442 ReadFields(mesh, objects, stFlds);
443
444
445
446 if (mergePatchIDs.size())
447 {
448 Info<< "Merging duplicate faces" << nl << endl;
449
450 // Mesh change engine
451 polyTopoChange meshMod(mesh);
452
453 const labelList boundaryFaces(patchFaces(mesh, mergePatchIDs));
454
455 // Get all duplicate face pairs (in boundaryFaces indices!).
456 labelList duplicates(findBaffles(mesh, boundaryFaces));
457
458 // Merge into internal faces.
459 insertDuplicateMerge(mesh, boundaryFaces, duplicates, meshMod);
460
461 if (!overwrite)
462 {
463 ++runTime;
464 }
465
466 // Change the mesh. No inflation.
467 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
468
469 // Update fields
470 mesh.updateMesh(map());
471
472 // Move mesh (since morphing does not do this)
473 if (map().hasMotionPoints())
474 {
475 mesh.movePoints(map().preMotionPoints());
476 }
477
478 if (overwrite)
479 {
480 mesh.setInstance(oldInstance);
481 }
482 Info<< "Writing mesh to time " << runTime.timeName() << endl;
483 mesh.write();
484 }
485
486
487 if (splitPatchIDs.size())
488 {
489 Info<< "Topologically splitting duplicate surfaces"
490 << ", i.e. duplicating points internal to duplicate surfaces"
491 << nl << endl;
492
493 // Determine points on split patches
494 DynamicList<label> candidates;
495 {
496 label sz = 0;
497 forAll(splitPatchIDs, i)
498 {
499 sz += patches[splitPatchIDs[i]].nPoints();
500 }
501 candidates.setCapacity(sz);
502
503 bitSet isCandidate(mesh.nPoints());
504 forAll(splitPatchIDs, i)
505 {
506 const labelList& mp = patches[splitPatchIDs[i]].meshPoints();
507 forAll(mp, mpi)
508 {
509 label pointi = mp[mpi];
510 if (isCandidate.set(pointi))
511 {
512 candidates.append(pointi);
513 }
514 }
515 }
516 }
517
518
519 // Analyse which points need to be duplicated
520 localPointRegion regionSide(mesh, candidates);
521
522 // Point duplication engine
523 duplicatePoints pointDuplicator(mesh);
524
525 // Mesh change engine
526 polyTopoChange meshMod(mesh);
527
528 // Insert topo changes
529 pointDuplicator.setRefinement(regionSide, meshMod);
530
531 if (!overwrite)
532 {
533 ++runTime;
534 }
535
536 // Change the mesh. No inflation.
537 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
538
539 // Update fields
540 mesh.updateMesh(map());
541
542 // Move mesh (since morphing does not do this)
543 if (map().hasMotionPoints())
544 {
545 mesh.movePoints(map().preMotionPoints());
546 }
547
548 if (overwrite)
549 {
550 mesh.setInstance(oldInstance);
551 }
552 Info<< "Writing mesh to time " << runTime.timeName() << endl;
553 mesh.write();
554
555 topoSet::removeFiles(mesh);
556 processorMeshes::removeFiles(mesh);
557
558 // Dump duplicated points (if any)
559 const labelList& pointMap = map().pointMap();
560
561 labelList nDupPerPoint(map().nOldPoints(), Zero);
562
563 pointSet dupPoints(mesh, "duplicatedPoints", 100);
564
565 forAll(pointMap, pointi)
566 {
567 label oldPointi = pointMap[pointi];
568
569 nDupPerPoint[oldPointi]++;
570
571 if (nDupPerPoint[oldPointi] > 1)
572 {
573 dupPoints.insert(map().reversePointMap()[oldPointi]);
574 dupPoints.insert(pointi);
575 }
576 }
577
578 Info<< "Writing " << returnReduce(dupPoints.size(), sumOp<label>())
579 << " duplicated points to pointSet "
580 << dupPoints.objectPath() << nl << endl;
581
582 dupPoints.write();
583 }
584
585 Info<< "End\n" << endl;
586
587 return 0;
588}
589
590
591// ************************************************************************* //
static bool split(const std::string &line, std::string &key, std::string &val)
Definition: cpuInfo.C:39
Field reading functions for post-processing utilities.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
void setCapacity(const label len)
Alter the size of the underlying storage.
Definition: DynamicListI.H:303
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
List of IOobjects with searching and retrieving facilities.
Definition: IOobjectList.H:59
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: ZoneMesh.C:288
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
Duplicate points.
A list of face labels.
Definition: faceSet.H:54
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:67
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:372
const boolList & flipMap() const noexcept
Return face flip map.
Definition: faceZone.H:272
Takes mesh with 'baffles' (= boundary faces sharing points). Determines for selected points on bounda...
A set of point labels.
Definition: pointSet.H:54
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
Class describing modification of a face.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:364
Class containing data for face removal.
Direct mesh changes based on v1.3 polyTopoChange syntax.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
autoPtr< mapPolyMesh > changeMesh(polyMesh &mesh, const labelUList &patchMap, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
Determines the 'side' for every face and connected to a singly-connected (through edges) region of fa...
Definition: regionSide.H:64
A List of wordRe with additional matching capabilities.
Definition: wordRes.H:54
A class for handling words, derived from Foam::string.
Definition: word.H:68
const polyBoundaryMesh & patches
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const word dictName("faMeshDefinition")
const labelIOList & zoneID
const dimensionedScalar mp
Proton mass.
Namespace for OpenFOAM.
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.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
errorManip< error > abort(error &err)
Definition: errorManip.H:144
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
dictionary dict
IOobject dictIO
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Foam::surfaceFields.