extrudeMesh.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) 2015-2022 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 extrudeMesh
29
30Group
31 grpMeshGenerationUtilities
32
33Description
34 Extrude mesh from existing patch (by default outwards facing normals;
35 optional flips faces) or from patch read from file.
36
37 Note: Merges close points so be careful.
38
39 Type of extrusion prescribed by run-time selectable model.
40
41\*---------------------------------------------------------------------------*/
42
43#include "argList.H"
44#include "Time.H"
45#include "polyTopoChange.H"
46#include "polyTopoChanger.H"
47#include "edgeCollapser.H"
48#include "perfectInterface.H"
49#include "addPatchCellLayer.H"
50#include "fvMesh.H"
51#include "MeshedSurfaces.H"
52#include "globalIndex.H"
53#include "cellSet.H"
54#include "fvMeshTools.H"
55
56#include "extrudedMesh.H"
57#include "extrudeModel.H"
58
59#include "wedge.H"
60#include "wedgePolyPatch.H"
61#include "planeExtrusion.H"
62#include "emptyPolyPatch.H"
63#include "processorPolyPatch.H"
64#include "processorMeshes.H"
65#include "hexRef8Data.H"
66
67using namespace Foam;
68
69// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
70
71enum ExtrudeMode
72{
73 MESH,
74 PATCH,
75 SURFACE
76};
77
78static const Enum<ExtrudeMode> ExtrudeModeNames
79{
80 { ExtrudeMode::MESH, "mesh" },
81 { ExtrudeMode::PATCH, "patch" },
82 { ExtrudeMode::SURFACE, "surface" },
83};
84
85
86label findPatchID(const polyBoundaryMesh& patches, const word& name)
87{
88 const label patchID = patches.findPatchID(name);
89
90 if (patchID == -1)
91 {
93 << "Cannot find patch " << name
94 << " in the source mesh.\n"
95 << "Valid patch names are " << patches.names()
96 << exit(FatalError);
97 }
98 return patchID;
99}
100
101
102labelList patchFaces(const polyBoundaryMesh& patches, const wordList& names)
103{
104 label n = 0;
105
106 forAll(names, i)
107 {
108 const polyPatch& pp = patches[findPatchID(patches, names[i])];
109
110 n += pp.size();
111 }
112 labelList faceLabels(n);
113 n = 0;
114 forAll(names, i)
115 {
116 const polyPatch& pp = patches[findPatchID(patches, names[i])];
117
118 forAll(pp, j)
119 {
120 faceLabels[n++] = pp.start()+j;
121 }
122 }
123
124 return faceLabels;
125}
126
127
128void zoneFaces
129(
130 const faceZoneMesh& fzs,
131 const wordList& names,
132 labelList& faceLabels,
133 bitSet& faceFlip
134)
135{
136 label n = 0;
137
138 forAll(names, i)
139 {
140 const auto& pp = fzs[fzs.findZoneID(names[i])];
141 n += pp.size();
142 }
143 faceLabels.setSize(n);
144 faceFlip.setSize(n);
145 n = 0;
146 forAll(names, i)
147 {
148 const auto& pp = fzs[fzs.findZoneID(names[i])];
149 const boolList& ppFlip = pp.flipMap();
150 forAll(pp, i)
151 {
152 faceLabels[n] = pp[i];
153 faceFlip[n] = ppFlip[i];
154 n++;
155 }
156 }
157}
158
159
160void updateFaceLabels(const mapPolyMesh& map, labelList& faceLabels)
161{
162 const labelList& reverseMap = map.reverseFaceMap();
163
164 labelList newFaceLabels(faceLabels.size());
165 label newI = 0;
166
167 forAll(faceLabels, i)
168 {
169 label oldFacei = faceLabels[i];
170
171 if (reverseMap[oldFacei] >= 0)
172 {
173 newFaceLabels[newI++] = reverseMap[oldFacei];
174 }
175 }
176 newFaceLabels.setSize(newI);
177 faceLabels.transfer(newFaceLabels);
178}
179
180
181void updateCellSet(const mapPolyMesh& map, labelHashSet& cellLabels)
182{
183 const labelList& reverseMap = map.reverseCellMap();
184
185 labelHashSet newCellLabels(2*cellLabels.size());
186
187 forAll(cellLabels, i)
188 {
189 label oldCelli = cellLabels[i];
190
191 if (reverseMap[oldCelli] >= 0)
192 {
193 newCellLabels.insert(reverseMap[oldCelli]);
194 }
195 }
196 cellLabels.transfer(newCellLabels);
197}
198
199
200template<class PatchType>
201void changeFrontBackPatches
202(
203 polyMesh& mesh,
204 const word& frontPatchName,
205 const word& backPatchName
206)
207{
208 const polyBoundaryMesh& patches = mesh.boundaryMesh();
209
210 label frontPatchi = findPatchID(patches, frontPatchName);
211 label backPatchi = findPatchID(patches, backPatchName);
212
213 DynamicList<polyPatch*> newPatches(patches.size());
214
215 forAll(patches, patchi)
216 {
217 const polyPatch& pp(patches[patchi]);
218
219 if (patchi == frontPatchi || patchi == backPatchi)
220 {
221 newPatches.append
222 (
223 new PatchType
224 (
225 pp.name(),
226 pp.size(),
227 pp.start(),
228 pp.index(),
229 mesh.boundaryMesh(),
230 PatchType::typeName
231 )
232 );
233 }
234 else
235 {
236 newPatches.append(pp.clone(mesh.boundaryMesh()).ptr());
237 }
238 }
239
240 // Edit patches
241 mesh.removeBoundary();
242 mesh.addPatches(newPatches, true);
243}
244
245
246int main(int argc, char *argv[])
247{
248 argList::addNote
249 (
250 "Extrude mesh from existing patch."
251 );
252
253 #include "addRegionOption.H"
254 argList::addOption("dict", "file", "Alternative extrudeMeshDict");
255 #include "setRootCase.H"
256 #include "createTimeExtruded.H"
257
258 // Get optional regionName
259 word regionName(polyMesh::defaultRegion);
260 if (args.readIfPresent("region", regionName))
261 {
262 Info<< "Create mesh " << regionName << " for time = "
263 << runTimeExtruded.timeName() << nl << endl;
264 }
265 else
266 {
267 Info<< "Create mesh for time = "
268 << runTimeExtruded.timeName() << nl << endl;
269 }
270
271
272 const IOdictionary dict
273 (
274 IOobject::selectIO
275 (
277 (
278 "extrudeMeshDict",
279 runTimeExtruded.system(),
280 runTimeExtruded,
281 IOobject::MUST_READ_IF_MODIFIED,
282 IOobject::NO_WRITE
283 ),
284 args.getOrDefault<fileName>("dict", "")
285 )
286 );
287
288 // Point generator
289 autoPtr<extrudeModel> model(extrudeModel::New(dict));
290
291 // Whether to flip normals
292 const bool flipNormals(dict.get<bool>("flipNormals"));
293
294 // What to extrude
295 const ExtrudeMode mode = ExtrudeModeNames.get("constructFrom", dict);
296
297 // Any merging of small edges
298 const scalar mergeTol(dict.getOrDefault<scalar>("mergeTol", 1e-4));
299
300 Info<< "Extruding from " << ExtrudeModeNames[mode]
301 << " using model " << model().type() << endl;
302 if (flipNormals)
303 {
304 Info<< "Flipping normals before extruding" << endl;
305 }
306 if (mergeTol > 0)
307 {
308 Info<< "Collapsing edges < " << mergeTol << " of bounding box" << endl;
309 }
310 else
311 {
312 Info<< "Not collapsing any edges after extrusion" << endl;
313 }
314 Info<< endl;
315
316
317 // Generated mesh (one of either)
318 autoPtr<fvMesh> meshFromMesh;
319 autoPtr<polyMesh> meshFromSurface;
320
321 // Faces on front and back for stitching (in case of mergeFaces)
322 word frontPatchName;
323 labelList frontPatchFaces;
324 word backPatchName;
325 labelList backPatchFaces;
326
327 // Optional added cells (get written to cellSet)
328 labelHashSet addedCellsSet;
329
330 // Optional refinement data
331 autoPtr<hexRef8Data> refDataPtr;
332
333 if (mode == PATCH || mode == MESH)
334 {
335 if (flipNormals && mode == MESH)
336 {
338 << "Flipping normals not supported for extrusions from mesh."
339 << exit(FatalError);
340 }
341
342 fileName sourceCasePath(dict.get<fileName>("sourceCase").expand());
343 fileName sourceRootDir = sourceCasePath.path();
344 fileName sourceCaseDir = sourceCasePath.name();
345 if (Pstream::parRun())
346 {
347 sourceCaseDir =
348 sourceCaseDir/("processor" + Foam::name(Pstream::myProcNo()));
349 }
350 wordList sourcePatches;
351 wordList sourceFaceZones;
352 if
353 (
354 dict.readIfPresent
355 (
356 "sourcePatches",
357 sourcePatches,
358 keyType::LITERAL
359 )
360 )
361 {
362 if (sourcePatches.size() == 1)
363 {
364 frontPatchName = sourcePatches[0];
365 }
366
367 Info<< "Extruding patches " << sourcePatches
368 << " on mesh " << sourceCasePath << nl
369 << endl;
370 }
371 else
372 {
373 dict.readEntry("sourceFaceZones", sourceFaceZones);
374 Info<< "Extruding faceZones " << sourceFaceZones
375 << " on mesh " << sourceCasePath << nl
376 << endl;
377 }
378
379
381 (
382 Time::controlDictName,
383 sourceRootDir,
384 sourceCaseDir
385 );
386
387 #include "createNamedMesh.H"
388
389 const polyBoundaryMesh& patches = mesh.boundaryMesh();
390
391
392 // Extrusion engine. Either adding to existing mesh or
393 // creating separate mesh.
394 addPatchCellLayer layerExtrude(mesh, (mode == MESH));
395
396 labelList meshFaces;
397 bitSet faceFlips;
398 if (sourceFaceZones.size())
399 {
400 zoneFaces(mesh.faceZones(), sourceFaceZones, meshFaces, faceFlips);
401 }
402 else
403 {
404 meshFaces = patchFaces(patches, sourcePatches);
405 faceFlips.setSize(meshFaces.size());
406 }
407
408 if (mode == PATCH && flipNormals)
409 {
410 // Cheat. Flip patch faces in mesh. This invalidates the
411 // mesh (open cells) but does produce the correct extrusion.
412 polyTopoChange meshMod(mesh);
413 forAll(meshFaces, i)
414 {
415 label meshFacei = meshFaces[i];
416
417 label patchi = patches.whichPatch(meshFacei);
418 label own = mesh.faceOwner()[meshFacei];
419 label nei = -1;
420 if (patchi == -1)
421 {
422 nei = mesh.faceNeighbour()[meshFacei];
423 }
424
425 label zoneI = mesh.faceZones().whichZone(meshFacei);
426 bool zoneFlip = false;
427 if (zoneI != -1)
428 {
429 label index = mesh.faceZones()[zoneI].whichFace(meshFacei);
430 zoneFlip = mesh.faceZones()[zoneI].flipMap()[index];
431 }
432
433 meshMod.modifyFace
434 (
435 mesh.faces()[meshFacei].reverseFace(), // modified face
436 meshFacei, // label of face
437 own, // owner
438 nei, // neighbour
439 true, // face flip
440 patchi, // patch for face
441 zoneI, // zone for face
442 zoneFlip // face flip in zone
443 );
444 }
445
446 // Change the mesh. No inflation.
447 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
448
449 // Update fields
450 mesh.updateMesh(map());
451
452 // Move mesh (since morphing does not do this)
453 if (map().hasMotionPoints())
454 {
455 mesh.movePoints(map().preMotionPoints());
456 }
457 }
458
459
460 indirectPrimitivePatch extrudePatch
461 (
463 (
464 mesh.faces(),
465 meshFaces
466 ),
467 mesh.points()
468 );
469
470 // Determine extrudePatch normal
471 pointField extrudePatchPointNormals
472 (
473 PatchTools::pointNormals(mesh, extrudePatch, faceFlips)
474 );
475
476
477 // Precalculate mesh edges for pp.edges.
478 const labelList meshEdges
479 (
480 extrudePatch.meshEdges
481 (
482 mesh.edges(),
483 mesh.pointEdges()
484 )
485 );
486
487 // Global face indices engine
488 const globalIndex globalFaces(mesh.nFaces());
489
490 // Determine extrudePatch.edgeFaces in global numbering (so across
491 // coupled patches)
492 labelListList edgeGlobalFaces
493 (
494 addPatchCellLayer::globalEdgeFaces
495 (
496 mesh,
497 globalFaces,
498 extrudePatch
499 )
500 );
501
502
503 // Determine what patches boundary edges need to get extruded into.
504 // This might actually cause edge-connected processors to become
505 // face-connected so might need to introduce new processor boundaries.
506 // Calculates:
507 // - per pp.edge the patch to extrude into
508 // - any additional processor boundaries (patchToNbrProc = map from
509 // new patchID to neighbour processor)
510 // - number of new patches (nPatches)
511
512 labelList edgePatchID;
513 labelList edgeZoneID;
514 boolList edgeFlip;
515 labelList inflateFaceID;
516 label nPatches;
517 Map<label> nbrProcToPatch;
518 Map<label> patchToNbrProc;
519 addPatchCellLayer::calcExtrudeInfo
520 (
521 true, // for internal edges get zone info from any face
522
523 mesh,
524 globalFaces,
525 edgeGlobalFaces,
526 extrudePatch,
527
528 edgePatchID,
529 nPatches,
530 nbrProcToPatch,
531 patchToNbrProc,
532
533 edgeZoneID,
534 edgeFlip,
535 inflateFaceID
536 );
537
538
539 // Add any patches.
540
541 label nAdded = nPatches - mesh.boundaryMesh().size();
542 reduce(nAdded, sumOp<label>());
543
544 Info<< "Adding overall " << nAdded << " processor patches." << endl;
545
546 if (nAdded > 0)
547 {
549 forAll(mesh.boundaryMesh(), patchi)
550 {
551 newPatches.append
552 (
553 mesh.boundaryMesh()[patchi].clone
554 (
555 mesh.boundaryMesh()
556 ).ptr()
557 );
558 }
559 for
560 (
561 label patchi = mesh.boundaryMesh().size();
562 patchi < nPatches;
563 patchi++
564 )
565 {
566 label nbrProci = patchToNbrProc[patchi];
567
568 Pout<< "Adding patch " << patchi
569 << " between " << Pstream::myProcNo()
570 << " and " << nbrProci
571 << endl;
572
573 newPatches.append
574 (
576 (
577 0, // size
578 mesh.nFaces(), // start
579 patchi, // index
580 mesh.boundaryMesh(),// polyBoundaryMesh
581 Pstream::myProcNo(),// myProcNo
582 nbrProci // neighbProcNo
583 )
584 );
585 }
586
587 // Add patches. Do no parallel updates.
588 mesh.removeFvBoundary();
589 mesh.addFvPatches(newPatches, true);
590 }
591
592
593
594 // Only used for addPatchCellLayer into new mesh
595 labelList exposedPatchID;
596 if (mode == PATCH)
597 {
598 dict.readEntry("exposedPatchName", backPatchName);
599 exposedPatchID.setSize
600 (
601 extrudePatch.size(),
602 findPatchID(patches, backPatchName)
603 );
604 }
605
606 // Determine points and extrusion
607 pointField layer0Points(extrudePatch.nPoints());
608 pointField displacement(extrudePatch.nPoints());
609 forAll(displacement, pointi)
610 {
611 const vector& patchNormal = extrudePatchPointNormals[pointi];
612
613 // layer0 point
614 layer0Points[pointi] = model()
615 (
616 extrudePatch.localPoints()[pointi],
617 patchNormal,
618 0
619 );
620 // layerN point
621 point extrudePt = model()
622 (
623 extrudePatch.localPoints()[pointi],
624 patchNormal,
625 model().nLayers()
626 );
627 displacement[pointi] = extrudePt - layer0Points[pointi];
628 }
629
630
631 // Check if wedge (has layer0 different from original patch points)
632 // If so move the mesh to starting position.
633 if (gMax(mag(layer0Points-extrudePatch.localPoints())) > SMALL)
634 {
635 Info<< "Moving mesh to layer0 points since differ from original"
636 << " points - this can happen for wedge extrusions." << nl
637 << endl;
638
639 pointField newPoints(mesh.points());
640 forAll(extrudePatch.meshPoints(), i)
641 {
642 newPoints[extrudePatch.meshPoints()[i]] = layer0Points[i];
643 }
644 mesh.movePoints(newPoints);
645 }
646
647
648 // Layers per face
649 labelList nFaceLayers(extrudePatch.size(), model().nLayers());
650
651 // Layers per point
652 labelList nPointLayers(extrudePatch.nPoints(), model().nLayers());
653
654 // Displacement for first layer
655 vectorField firstLayerDisp(displacement*model().sumThickness(1));
656
657 // Expansion ratio not used.
658 scalarField ratio(extrudePatch.nPoints(), 1.0);
659
660
661 // Load any refinement data
662 if (mode == MESH)
663 {
664 // Tricky: register hexRef8 onto the database
665 // since the mesh does not yet exist. It should not be registered
666 // onto the input mesh.
667 refDataPtr.reset
668 (
669 new hexRef8Data
670 (
672 (
673 "dummy",
674 mesh.facesInstance(),
675 polyMesh::meshSubDir,
676 runTimeExtruded, //mesh,
677 IOobject::READ_IF_PRESENT,
678 IOobject::NO_WRITE,
679 false
680 )
681 )
682 );
683 }
684
685
686 // Topo change container. Either copy an existing mesh or start
687 // with empty storage (number of patches only needed for checking)
689 (
690 (
691 mode == MESH
692 ? new polyTopoChange(mesh)
693 : new polyTopoChange(patches.size())
694 )
695 );
696
697 layerExtrude.setRefinement
698 (
699 globalFaces,
700 edgeGlobalFaces,
701
702 ratio, // expansion ratio
703 extrudePatch, // patch faces to extrude
704 faceFlips, // side to extrude (for internal/coupled faces)
705
706 edgePatchID, // if boundary edge: patch for extruded face
707 edgeZoneID, // optional zone for extruded face
708 edgeFlip,
709 inflateFaceID, // mesh face that zone/patch info is from
710
711 exposedPatchID, // if new mesh: patches for exposed faces
712 nFaceLayers,
713 nPointLayers,
714 firstLayerDisp,
715 meshMod()
716 );
717
718 // Reset points according to extrusion model
719 forAll(layerExtrude.addedPoints(), pointi)
720 {
721 const labelList& pPoints = layerExtrude.addedPoints()[pointi];
722 forAll(pPoints, pPointi)
723 {
724 label meshPointi = pPoints[pPointi];
725
726 point modelPt
727 (
728 model()
729 (
730 extrudePatch.localPoints()[pointi],
731 extrudePatchPointNormals[pointi],
732 pPointi+1 // layer
733 )
734 );
735
736 const_cast<DynamicList<point>&>
737 (
738 meshMod().points()
739 )[meshPointi] = modelPt;
740 }
741 }
742
743 // Store faces on front and exposed patch (if mode=patch there are
744 // only added faces so cannot used map to old faces)
745 const labelListList& layerFaces = layerExtrude.layerFaces();
746 backPatchFaces.setSize(layerFaces.size());
747 frontPatchFaces.setSize(layerFaces.size());
748 forAll(backPatchFaces, patchFacei)
749 {
750 backPatchFaces[patchFacei] = layerFaces[patchFacei].first();
751 frontPatchFaces[patchFacei] = layerFaces[patchFacei].last();
752 }
753
754
755 // Create dummy fvSchemes, fvSolution
756 fvMeshTools::createDummyFvMeshFiles
757 (
758 mesh,
759 polyMesh::regionName(regionName),
760 true
761 );
762
763 // Create actual mesh from polyTopoChange container
764 autoPtr<mapPolyMesh> map = meshMod().makeMesh
765 (
766 meshFromMesh,
768 (
770 runTimeExtruded.constant(),
771 runTimeExtruded,
772 IOobject::READ_IF_PRESENT, // Read fv* if present
773 IOobject::AUTO_WRITE,
774 false
775 ),
776 mesh
777 );
778
779 layerExtrude.updateMesh
780 (
781 map(),
782 identity(extrudePatch.size()),
783 identity(extrudePatch.nPoints())
784 );
785
786 // Calculate face labels for front and back.
787 frontPatchFaces = renumber
788 (
789 map().reverseFaceMap(),
790 frontPatchFaces
791 );
792 backPatchFaces = renumber
793 (
794 map().reverseFaceMap(),
795 backPatchFaces
796 );
797
798 // Update
799 if (refDataPtr)
800 {
801 refDataPtr->updateMesh(map());
802 }
803
804 // Store added cells
805 if (mode == MESH)
806 {
807 const labelListList addedCells
808 (
809 layerExtrude.addedCells
810 (
811 *meshFromMesh,
812 layerExtrude.layerFaces()
813 )
814 );
815 forAll(addedCells, facei)
816 {
817 const labelList& aCells = addedCells[facei];
818 addedCellsSet.insert(aCells);
819 }
820 }
821 }
822 else
823 {
824 // Read from surface
825 fileName surfName(dict.get<fileName>("surface").expand());
826
827 Info<< "Extruding surfaceMesh read from file " << surfName << nl
828 << endl;
829
830 MeshedSurface<face> fMesh(surfName);
831
832 if (flipNormals)
833 {
834 Info<< "Flipping faces." << nl << endl;
835 faceList& faces = const_cast<faceList&>(fMesh.surfFaces());
836 forAll(faces, i)
837 {
838 faces[i] = fMesh[i].reverseFace();
839 }
840 }
841
842 Info<< "Extruding surface with :" << nl
843 << " points : " << fMesh.points().size() << nl
844 << " faces : " << fMesh.size() << nl
845 << " normals[0] : " << fMesh.faceNormals()[0]
846 << nl
847 << endl;
848
849 meshFromSurface.reset
850 (
851 new extrudedMesh
852 (
854 (
855 extrudedMesh::defaultRegion,
856 runTimeExtruded.constant(),
857 runTimeExtruded
858 ),
859 fMesh,
860 model()
861 )
862 );
863
864
865 // Get the faces on front and back
866 frontPatchName = "originalPatch";
867 frontPatchFaces = patchFaces
868 (
869 meshFromSurface().boundaryMesh(),
870 wordList(1, frontPatchName)
871 );
872 backPatchName = "otherSide";
873 backPatchFaces = patchFaces
874 (
875 meshFromSurface().boundaryMesh(),
876 wordList(1, backPatchName)
877 );
878 }
879
880
881 polyMesh& mesh =
882 (
883 meshFromMesh
884 ? *meshFromMesh
885 : *meshFromSurface
886 );
887
888
889 const boundBox& bb = mesh.bounds();
890 const vector span = bb.span();
891 const scalar mergeDim = mergeTol * bb.minDim();
892
893 Info<< "Mesh bounding box : " << bb << nl
894 << " with span : " << span << nl
895 << "Merge distance : " << mergeDim << nl
896 << endl;
897
898
899 // Collapse edges
900 // ~~~~~~~~~~~~~~
901
902 if (mergeDim > 0)
903 {
904 Info<< "Collapsing edges < " << mergeDim << " ..." << nl << endl;
905
906 // Edge collapsing engine
907 edgeCollapser collapser(mesh);
908
909 const edgeList& edges = mesh.edges();
910 const pointField& points = mesh.points();
911
912 bitSet collapseEdge(mesh.nEdges());
913 Map<point> collapsePointToLocation(mesh.nPoints());
914
915 forAll(edges, edgeI)
916 {
917 const edge& e = edges[edgeI];
918
919 scalar d = e.mag(points);
920
921 if (d < mergeDim)
922 {
923 Info<< "Merging edge " << e << " since length " << d
924 << " << " << mergeDim << nl;
925
926 collapseEdge.set(edgeI);
927 collapsePointToLocation.set(e[1], points[e[0]]);
928 }
929 }
930
931 List<pointEdgeCollapse> allPointInfo;
932 const globalIndex globalPoints(mesh.nPoints());
933 labelList pointPriority(mesh.nPoints(), Zero);
934
935 collapser.consistentCollapse
936 (
938 pointPriority,
939 collapsePointToLocation,
941 allPointInfo
942 );
943
944 // Topo change container
945 polyTopoChange meshMod(mesh);
946
947 // Put all modifications into meshMod
948 bool anyChange = collapser.setRefinement(allPointInfo, meshMod);
949 reduce(anyChange, orOp<bool>());
950
951 if (anyChange)
952 {
953 // Construct new mesh from polyTopoChange.
954 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
955
956 // Update fields
957 mesh.updateMesh(map());
958
959 // Update stored data
960 updateFaceLabels(map(), frontPatchFaces);
961 updateFaceLabels(map(), backPatchFaces);
962 updateCellSet(map(), addedCellsSet);
963
964 if (refDataPtr)
965 {
966 refDataPtr->updateMesh(map());
967 }
968
969 // Move mesh (if inflation used)
970 if (map().hasMotionPoints())
971 {
972 mesh.movePoints(map().preMotionPoints());
973 }
974 }
975 }
976
977
978 // Change the front and back patch types as required
979 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
980
981 word frontBackType(word::null);
982
983 if (isType<extrudeModels::wedge>(model()))
984 {
985 changeFrontBackPatches<wedgePolyPatch>
986 (
987 mesh,
988 frontPatchName,
989 backPatchName
990 );
991 }
992 else if (isType<extrudeModels::plane>(model()))
993 {
994 changeFrontBackPatches<emptyPolyPatch>
995 (
996 mesh,
997 frontPatchName,
998 backPatchName
999 );
1000 }
1001
1002
1003 // Merging front and back patch faces
1004 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1005
1006 if (dict.get<bool>("mergeFaces"))
1007 {
1008 if (mode == MESH)
1009 {
1011 << "Cannot stitch front and back of extrusion since"
1012 << " in 'mesh' mode (extrusion appended to mesh)."
1013 << exit(FatalError);
1014 }
1015
1016 Info<< "Assuming full 360 degree axisymmetric case;"
1017 << " stitching faces on patches "
1018 << frontPatchName << " and "
1019 << backPatchName << " together ..." << nl << endl;
1020
1021 if (frontPatchFaces.size() != backPatchFaces.size())
1022 {
1024 << "Differing number of faces on front ("
1025 << frontPatchFaces.size() << ") and back ("
1026 << backPatchFaces.size() << ")"
1027 << exit(FatalError);
1028 }
1029
1030
1031 polyTopoChanger stitcher(mesh);
1032 stitcher.setSize(1);
1033
1034 const word cutZoneName("originalCutFaceZone");
1035
1037 (
1038 1,
1039 new faceZone
1040 (
1041 cutZoneName,
1042 frontPatchFaces,
1043 false, // none are flipped
1044 0,
1045 mesh.faceZones()
1046 )
1047 );
1048
1049 mesh.addZones(List<pointZone*>(0), fz, List<cellZone*>(0));
1050
1051 // Add the perfect interface mesh modifier
1052 perfectInterface perfectStitcher
1053 (
1054 "couple",
1055 0,
1056 stitcher,
1057 cutZoneName,
1058 word::null, // dummy patch name
1059 word::null // dummy patch name
1060 );
1061
1062 // Topo change container
1063 polyTopoChange meshMod(mesh);
1064
1065 perfectStitcher.setRefinement
1066 (
1068 (
1070 (
1071 mesh.faces(),
1072 frontPatchFaces
1073 ),
1074 mesh.points()
1075 ),
1077 (
1079 (
1080 mesh.faces(),
1081 backPatchFaces
1082 ),
1083 mesh.points()
1084 ),
1085 meshMod
1086 );
1087
1088 // Construct new mesh from polyTopoChange.
1089 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
1090
1091 // Update fields
1092 mesh.updateMesh(map());
1093
1094 // Update local data
1095 updateCellSet(map(), addedCellsSet);
1096
1097 if (refDataPtr)
1098 {
1099 refDataPtr->updateMesh(map());
1100 }
1101
1102 // Move mesh (if inflation used)
1103 if (map().hasMotionPoints())
1104 {
1105 mesh.movePoints(map().preMotionPoints());
1106 }
1107 }
1108
1109 mesh.setInstance(runTimeExtruded.constant());
1110 Info<< "Writing mesh to " << mesh.objectPath() << nl << endl;
1111
1112 if (!mesh.write())
1113 {
1115 << exit(FatalError);
1116 }
1117 // Remove any left-over files
1118 processorMeshes::removeFiles(mesh);
1119
1120 // Need writing cellSet
1121 label nAdded = returnReduce(addedCellsSet.size(), sumOp<label>());
1122 if (nAdded > 0)
1123 {
1124 cellSet addedCells(mesh, "addedCells", addedCellsSet);
1125 Info<< "Writing added cells to cellSet " << addedCells.name()
1126 << nl << endl;
1127 if (!addedCells.write())
1128 {
1130 << addedCells.name()
1131 << exit(FatalError);
1132 }
1133 }
1134
1135 if (refDataPtr)
1136 {
1137 refDataPtr->write();
1138 }
1139
1140
1141 Info<< "End\n" << endl;
1142
1143 return 0;
1144}
1145
1146
1147// ************************************************************************* //
label n
reduce(hasMovingMesh, orOp< bool >())
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
Definition: Enum.H:61
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:191
void transfer(HashTable< T, Key, Hash > &rhs)
Transfer contents into this table.
Definition: HashTable.C:719
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
A List with indirect addressing.
Definition: IndirectList.H:119
void transfer(List< T > &list)
Definition: List.C:447
void setSize(const label n)
Alias for resize()
Definition: List.H:218
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
A surface geometry mesh with zone information, not to be confused with the similarly named surfaceMes...
Definition: MeshedSurface.H:99
Bookkeeping for patch definitions.
Definition: PDRpatchDef.H:54
void setSize(const label n, unsigned int val=0u)
Alias for resize()
Definition: PackedList.H:557
A list of faces which address into the list of points.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:80
T & first()
Return the first element of the list.
Definition: UListI.H:202
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
T & last()
Return the last element of the list.
Definition: UListI.H:216
label findZoneID(const word &zoneName) const
Find zone index by name, return -1 if not found.
Definition: ZoneMesh.C:525
Adds layers of cells to outside of polyPatch. Can optionally create stand-alone extruded mesh (addToM...
bool readIfPresent(const word &optName, T &val) const
Read a value from the named option if present.
Definition: argListI.H:323
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Definition: argListI.H:307
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
void reset(autoPtr< T > &&other) noexcept
Delete managed object and set to new given pointer.
Definition: autoPtrI.H:117
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
scalar minDim() const
Smallest length/height/width dimension.
Definition: boundBoxI.H:145
vector span() const
The bounding box span (from minimum to maximum)
Definition: boundBoxI.H:127
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface....
Definition: boundaryMesh.H:63
A collection of cell labels.
Definition: cellSet.H:54
Does polyTopoChanges to remove edges. Can remove faces due to edge collapse but can not remove cells ...
Definition: edgeCollapser.H:69
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:67
A class for handling file names.
Definition: fileName.H:76
static std::string path(const std::string &str)
Return directory path name (part before last /)
Definition: fileNameI.H:176
static std::string name(const std::string &str)
Return basename (part beyond last /), including its extension.
Definition: fileNameI.H:199
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:68
Calculates points shared by more than two processor patches or cyclic patches.
Definition: globalPoints.H:103
Various for reading/decomposing/reconstructing/distributing refinement data.
Definition: hexRef8Data.H:61
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:162
const labelList & reverseCellMap() const
Reverse cell map.
Definition: mapPolyMesh.H:532
const labelList & reverseFaceMap() const
Reverse face map.
Definition: mapPolyMesh.H:501
label index() const noexcept
The index of this patch in the boundaryMesh.
const word & name() const noexcept
The patch name.
Hack of attachDetach to couple patches when they perfectly align. Does not decouple....
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
virtual autoPtr< polyPatch > clone(const labelList &faceCells) const
Construct and return a clone, setting faceCells.
Definition: polyPatch.H:238
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:364
Direct mesh changes based on v1.3 polyTopoChange syntax.
List of mesh modifiers defining the mesh dynamics.
Neighbour processor patch.
string & expand(const bool allowEmpty=false)
Definition: string.C:173
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
Definition: word.H:68
label collapseEdge(triSurface &surf, const scalar minLen)
Keep collapsing all edges < minLen.
const polyBoundaryMesh & patches
dynamicFvMesh & mesh
engineTime & runTime
Foam::word regionName(Foam::polyMesh::defaultRegion)
Required Variables.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const label nPatches
const pointField & points
Namespace for OpenFOAM.
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)
IntListType renumber(const labelUList &oldToNew, const IntListType &input)
Renumber the values (not the indices) of a list.
mode_t mode(const fileName &name, const bool followLink=true)
Return the file mode, normally following symbolic links.
Definition: MSwindows.C:572
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Type gMax(const FieldField< Field, Type > &f)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
dictionary dict
Foam::argList args(argc, argv)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333