extrudeToRegionMesh.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 extrudeToRegionMesh
29
30Group
31 grpMeshGenerationUtilities
32
33Description
34 Extrude faceZones (internal or boundary faces) or faceSets (boundary faces
35 only) into a separate mesh (as a different region).
36
37 - used to e.g. extrude baffles (extrude internal faces) or create
38 liquid film regions.
39 - if extruding internal faces:
40 - create baffles in original mesh with mappedWall patches
41 - if extruding boundary faces:
42 - convert boundary faces to mappedWall patches
43 - extrude edges of faceZone as a <zone>_sidePatch
44 - extrude edges inbetween different faceZones as a
45 (nonuniformTransform)cyclic <zoneA>_<zoneB>
46 - extrudes into master direction (i.e. away from the owner cell
47 if flipMap is false)
48
49\verbatim
50
51Internal face extrusion
52-----------------------
53
54 +-------------+
55 | |
56 | |
57 +---AAAAAAA---+
58 | |
59 | |
60 +-------------+
61
62 AAA=faceZone to extrude.
63
64
65For the case of no flipMap the extrusion starts at owner and extrudes
66into the space of the neighbour:
67
68 +CCCCCCC+
69 | | <= extruded mesh
70 +BBBBBBB+
71
72 +-------------+
73 | |
74 | (neighbour) |
75 |___CCCCCCC___| <= original mesh (with 'baffles' added)
76 | BBBBBBB |
77 |(owner side) |
78 | |
79 +-------------+
80
81 BBB=mapped between owner on original mesh and new extrusion.
82 (zero offset)
83 CCC=mapped between neighbour on original mesh and new extrusion
84 (offset due to the thickness of the extruded mesh)
85
86For the case of flipMap the extrusion is the other way around: from the
87neighbour side into the owner side.
88
89
90Boundary face extrusion
91-----------------------
92
93 +--AAAAAAA--+
94 | |
95 | |
96 +-----------+
97
98 AAA=faceZone to extrude. E.g. slave side is owner side (no flipmap)
99
100becomes
101
102 +CCCCCCC+
103 | | <= extruded mesh
104 +BBBBBBB+
105
106 +--BBBBBBB--+
107 | | <= original mesh
108 | |
109 +-----------+
110
111 BBB=mapped between original mesh and new extrusion
112 CCC=polypatch
113
114
115Notes:
116 - when extruding cyclics with only one cell inbetween it does not
117 detect this as a cyclic since the face is the same face. It will
118 only work if the coupled edge extrudes a different face so if there
119 are more than 1 cell inbetween.
120
121\endverbatim
122
123\*---------------------------------------------------------------------------*/
124
125#include "argList.H"
126#include "fvMesh.H"
127#include "polyTopoChange.H"
128#include "OFstream.H"
129#include "meshTools.H"
130#include "mappedWallPolyPatch.H"
131#include "createShellMesh.H"
132#include "syncTools.H"
133#include "cyclicPolyPatch.H"
134#include "wedgePolyPatch.H"
136#include "extrudeModel.H"
137#include "globalIndex.H"
138#include "faceSet.H"
139
140#include "volFields.H"
141#include "surfaceFields.H"
142#include "pointFields.H"
143//#include "ReadFields.H"
144#include "fvMeshTools.H"
145#include "OBJstream.H"
146#include "PatchTools.H"
147#include "processorMeshes.H"
148
149using namespace Foam;
150
151// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
152
153label findPatchID(const UList<polyPatch*>& newPatches, const word& name)
154{
155 forAll(newPatches, i)
156 {
157 if (newPatches[i]->name() == name)
158 {
159 return i;
160 }
161 }
162 return -1;
163}
164
165
166#ifdef FULLDEBUG
167void printPatches(Ostream& os, const UList<polyPatch*>& newPatches)
168{
169 for (const polyPatch* ppPtr : newPatches)
170 {
171 const polyPatch& pp = *ppPtr;
172 os << " name:" << pp.name() << " index:" << pp.index()
173 << " start:" << pp.start() << " size:" << pp.size()
174 << " type:" << pp.type() << nl;
175 }
176}
177#endif
178
179
180template<class PatchType>
181label addPatch
182(
184 const word& patchName,
185 DynamicList<polyPatch*>& newPatches
186)
187{
188 label patchi = findPatchID(newPatches, patchName);
189
190 if (patchi != -1)
191 {
192 if (isA<PatchType>(*newPatches[patchi]))
193 {
194 // Already there
195 return patchi;
196 }
197 else
198 {
200 << "Already have patch " << patchName
201 << " but of type " << newPatches[patchi]->type()
202 << exit(FatalError);
203 }
204 }
205
206
207 patchi = newPatches.size();
208
209 label startFacei = 0;
210 if (patchi > 0)
211 {
212 const polyPatch& pp = *newPatches.last();
213 startFacei = pp.start()+pp.size();
214 }
215
216 #ifdef FULLDEBUG
217 Pout<< "addPatch : starting newPatches:"
218 << " patch:" << patchi << " startFace:" << startFacei << nl;
219 printPatches(Pout, newPatches);
220 Pout<< "*** end of addPatch:" << endl;
221 #endif
222
223 newPatches.append
224 (
225 polyPatch::New
226 (
227 PatchType::typeName,
228 patchName,
229 0, // size
230 startFacei, // nFaces
231 patchi,
232 patches
233 ).ptr()
234 );
235
236 return patchi;
237}
238
239
240template<class PatchType>
241label addPatch
242(
244 const word& patchName,
245 const dictionary& dict,
246 DynamicList<polyPatch*>& newPatches
247)
248{
249 label patchi = findPatchID(newPatches, patchName);
250
251 if (patchi != -1)
252 {
253 if (isA<PatchType>(*newPatches[patchi]))
254 {
255 // Already there
256 return patchi;
257 }
258 else
259 {
261 << "Already have patch " << patchName
262 << " but of type " << newPatches[patchi]->type()
263 << exit(FatalError);
264 }
265 }
266
267
268 patchi = newPatches.size();
269
270 label startFacei = 0;
271 if (patchi > 0)
272 {
273 const polyPatch& pp = *newPatches.last();
274 startFacei = pp.start()+pp.size();
275 }
276
277
278 #ifdef FULLDEBUG
279 Pout<< "addPatch : starting newPatches:"
280 << " patch:" << patchi << " startFace:" << startFacei << nl;
281 printPatches(Pout, newPatches);
282 Pout<< "*** end of addPatch:" << endl;
283 #endif
284
285
286 dictionary patchDict(dict);
287 patchDict.set("type", PatchType::typeName);
288 patchDict.set("nFaces", 0);
289 patchDict.set("startFace", startFacei);
290
291 newPatches.append
292 (
293 polyPatch::New
294 (
295 patchName,
296 patchDict,
297 patchi,
298 patches
299 ).ptr()
300 );
301
302 return patchi;
303}
304
305
306// Remove zero-sized patches
307void deleteEmptyPatches(fvMesh& mesh)
308{
309 const polyBoundaryMesh& patches = mesh.boundaryMesh();
310
311 wordList masterNames;
312 if (Pstream::master())
313 {
314 masterNames = patches.names();
315 }
316 Pstream::broadcast(masterNames);
317
318
319 labelList oldToNew(patches.size(), -1);
320 label usedI = 0;
321 label notUsedI = patches.size();
322
323 // Add all the non-empty, non-processor patches
324 forAll(masterNames, masterI)
325 {
326 label patchi = patches.findPatchID(masterNames[masterI]);
327
328 if (patchi != -1)
329 {
330 if (isA<processorPolyPatch>(patches[patchi]))
331 {
332 // Similar named processor patch? Not 'possible'.
333 if (patches[patchi].empty())
334 {
335 Pout<< "Deleting processor patch " << patchi
336 << " name:" << patches[patchi].name()
337 << endl;
338 oldToNew[patchi] = --notUsedI;
339 }
340 else
341 {
342 oldToNew[patchi] = usedI++;
343 }
344 }
345 else
346 {
347 // Common patch.
348 if (returnReduce(patches[patchi].empty(), andOp<bool>()))
349 {
350 Pout<< "Deleting patch " << patchi
351 << " name:" << patches[patchi].name()
352 << endl;
353 oldToNew[patchi] = --notUsedI;
354 }
355 else
356 {
357 oldToNew[patchi] = usedI++;
358 }
359 }
360 }
361 }
362
363 // Add remaining patches at the end
364 forAll(patches, patchi)
365 {
366 if (oldToNew[patchi] == -1)
367 {
368 // Unique to this processor. Note: could check that these are
369 // only processor patches.
370 if (patches[patchi].empty())
371 {
372 Pout<< "Deleting processor patch " << patchi
373 << " name:" << patches[patchi].name()
374 << endl;
375 oldToNew[patchi] = --notUsedI;
376 }
377 else
378 {
379 oldToNew[patchi] = usedI++;
380 }
381 }
382 }
383
384 fvMeshTools::reorderPatches(mesh, oldToNew, usedI, true);
385}
386
387
388// Check zone either all internal or all external faces
389void checkZoneInside
390(
391 const polyMesh& mesh,
392 const wordList& zoneNames,
393 const labelList& zoneID,
394 const labelList& extrudeMeshFaces,
395 const boolList& isInternal
396)
397{
398 forAll(zoneNames, i)
399 {
400 if (isInternal[i])
401 {
402 Info<< "Zone " << zoneNames[i] << " has internal faces" << endl;
403 }
404 else
405 {
406 Info<< "Zone " << zoneNames[i] << " has boundary faces" << endl;
407 }
408 }
409
410 forAll(extrudeMeshFaces, i)
411 {
412 label facei = extrudeMeshFaces[i];
413 label zoneI = zoneID[i];
414 if (isInternal[zoneI] != mesh.isInternalFace(facei))
415 {
417 << "Zone " << zoneNames[zoneI]
418 << " is not consistently all internal or all boundary faces."
419 << " Face " << facei << " at " << mesh.faceCentres()[facei]
420 << " is the first occurrence."
421 << exit(FatalError);
422 }
423 }
424}
425
426
427// Calculate global pp faces per pp edge.
428labelListList globalEdgeFaces
429(
430 const polyMesh& mesh,
431 const globalIndex& globalFaces,
432 const primitiveFacePatch& pp,
433 const labelList& ppMeshEdges
434)
435{
436 // From mesh edge to global pp face labels.
437 labelListList globalEdgeFaces(ppMeshEdges.size());
438
439 const labelListList& edgeFaces = pp.edgeFaces();
440
441 forAll(edgeFaces, edgeI)
442 {
443 // Store pp face and processor as unique tag.
444 globalEdgeFaces[edgeI] = globalFaces.toGlobal(edgeFaces[edgeI]);
445 }
446
447 // Synchronise across coupled edges.
448 syncTools::syncEdgeList
449 (
450 mesh,
451 ppMeshEdges,
452 globalEdgeFaces,
454 labelList() // null value
455 );
456
457 return globalEdgeFaces;
458}
459
460
461// Find a patch face that is not extruded. Return -1 if not found.
462label findUncoveredPatchFace
463(
464 const fvMesh& mesh,
465 const labelUIndList& extrudeMeshFaces, // mesh faces that are extruded
466 const label meshEdgeI // mesh edge
467)
468{
469 // Make set of extruded faces.
470 labelHashSet extrudeFaceSet(extrudeMeshFaces.size());
471 extrudeFaceSet.insert(extrudeMeshFaces);
472
473 const polyBoundaryMesh& pbm = mesh.boundaryMesh();
474 const labelList& eFaces = mesh.edgeFaces()[meshEdgeI];
475 forAll(eFaces, i)
476 {
477 label facei = eFaces[i];
478 label patchi = pbm.whichPatch(facei);
479
480 if
481 (
482 patchi != -1
483 && !pbm[patchi].coupled()
484 && !extrudeFaceSet.found(facei)
485 )
486 {
487 return facei;
488 }
489 }
490
491 return -1;
492}
493
494
495// Same as findUncoveredPatchFace, except explicitly checks for cyclic faces
496label findUncoveredCyclicPatchFace
497(
498 const fvMesh& mesh,
499 const labelUIndList& extrudeMeshFaces, // mesh faces that are extruded
500 const label meshEdgeI // mesh edge
501)
502{
503 // Make set of extruded faces.
504 labelHashSet extrudeFaceSet(extrudeMeshFaces.size());
505 extrudeFaceSet.insert(extrudeMeshFaces);
506
507 const polyBoundaryMesh& pbm = mesh.boundaryMesh();
508 const labelList& eFaces = mesh.edgeFaces()[meshEdgeI];
509 forAll(eFaces, i)
510 {
511 label facei = eFaces[i];
512 label patchi = pbm.whichPatch(facei);
513
514 if
515 (
516 patchi != -1
517 && isA<cyclicPolyPatch>(pbm[patchi])
518 && !extrudeFaceSet.found(facei)
519 )
520 {
521 return facei;
522 }
523 }
524
525 return -1;
526}
527
528
529// Calculate per edge min and max zone
530void calcEdgeMinMaxZone
531(
532 const fvMesh& mesh,
533 const primitiveFacePatch& extrudePatch,
534 const labelList& extrudeMeshEdges,
535 const labelList& zoneID,
536 const mapDistribute& extrudeEdgeFacesMap,
537 const labelListList& extrudeEdgeGlobalFaces,
538
539 labelList& minZoneID,
540 labelList& maxZoneID
541)
542{
543 // Get zoneIDs in extrudeEdgeGlobalFaces order
544 labelList mappedZoneID(zoneID);
545 extrudeEdgeFacesMap.distribute(mappedZoneID);
546
547 // Get min and max zone per edge
548 minZoneID.setSize(extrudeEdgeGlobalFaces.size(), labelMax);
549 maxZoneID.setSize(extrudeEdgeGlobalFaces.size(), labelMin);
550
551 forAll(extrudeEdgeGlobalFaces, edgeI)
552 {
553 const labelList& eFaces = extrudeEdgeGlobalFaces[edgeI];
554 if (eFaces.size())
555 {
556 forAll(eFaces, i)
557 {
558 label zoneI = mappedZoneID[eFaces[i]];
559 minZoneID[edgeI] = min(minZoneID[edgeI], zoneI);
560 maxZoneID[edgeI] = max(maxZoneID[edgeI], zoneI);
561 }
562 }
563 }
564 syncTools::syncEdgeList
565 (
566 mesh,
567 extrudeMeshEdges,
568 minZoneID,
570 labelMax // null value
571 );
572 syncTools::syncEdgeList
573 (
574 mesh,
575 extrudeMeshEdges,
576 maxZoneID,
578 labelMin // null value
579 );
580}
581
582
583// Count the number of faces in patches that need to be created. Calculates:
584// zoneSidePatch[zoneI] : the number of side faces to be created
585// zoneZonePatch[zoneA,zoneB] : the number of faces inbetween zoneA and B
586// Since this only counts we're not taking the processor patches into
587// account.
588void countExtrudePatches
589(
590 const fvMesh& mesh,
591 const label nZones,
592 const primitiveFacePatch& extrudePatch,
593 const labelList& extrudeMeshFaces,
594 const labelList& extrudeMeshEdges,
595
596 const labelListList& extrudeEdgeGlobalFaces,
597 const labelList& minZoneID,
598 const labelList& maxZoneID,
599
600 labelList& zoneSidePatch,
601 labelList& zoneZonePatch
602)
603{
604 // Check on master edge for use of zones. Since we only want to know
605 // whether they are being used at all no need to accurately count on slave
606 // edge as well. Just add all together at the end of this routine so it
607 // gets detected at least.
608
609 forAll(extrudePatch.edgeFaces(), edgeI)
610 {
611 const labelList& eFaces = extrudePatch.edgeFaces()[edgeI];
612
613 if (eFaces.size() == 2)
614 {
615 // Internal edge - check if inbetween different zones.
616 if (minZoneID[edgeI] != maxZoneID[edgeI])
617 {
618 zoneZonePatch[minZoneID[edgeI]*nZones+maxZoneID[edgeI]]++;
619 }
620 }
621 else if
622 (
623 eFaces.size() == 1
624 && extrudeEdgeGlobalFaces[edgeI].size() == 2
625 )
626 {
627 // Coupled edge - check if inbetween different zones.
628 if (minZoneID[edgeI] != maxZoneID[edgeI])
629 {
630 const edge& e = extrudePatch.edges()[edgeI];
631 const pointField& pts = extrudePatch.localPoints();
633 << "Edge " << edgeI
634 << "at " << pts[e[0]] << pts[e[1]]
635 << " is a coupled edge and inbetween two different zones "
636 << minZoneID[edgeI] << " and " << maxZoneID[edgeI] << endl
637 << " This is currently not supported." << endl;
638
639 zoneZonePatch[minZoneID[edgeI]*nZones+maxZoneID[edgeI]]++;
640 }
641 }
642 else
643 {
644 // One or more than two edge-faces.
645 // Check whether we are on a mesh edge with external patches. If
646 // so choose any uncovered one. If none found put face in
647 // undetermined zone 'side' patch
648
649 label facei = findUncoveredPatchFace
650 (
651 mesh,
652 labelUIndList(extrudeMeshFaces, eFaces),
653 extrudeMeshEdges[edgeI]
654 );
655
656 if (facei == -1)
657 {
658 zoneSidePatch[minZoneID[edgeI]]++;
659 }
660 }
661 }
662 // Synchronise decision. Actual numbers are not important, just make
663 // sure that they're > 0 on all processors.
664 Pstream::listCombineAllGather(zoneSidePatch, plusEqOp<label>());
665 Pstream::listCombineAllGather(zoneZonePatch, plusEqOp<label>());
666}
667
668
669void addCouplingPatches
670(
671 const fvMesh& mesh,
672 const word& regionName,
673 const word& shellRegionName,
674 const wordList& zoneNames,
675 const wordList& zoneShadowNames,
676 const boolList& isInternal,
677 const labelList& zoneIDs,
678
679 DynamicList<polyPatch*>& newPatches,
680 labelList& interRegionTopPatch,
681 labelList& interRegionBottomPatch
682)
683{
684 Pout<< "Adding coupling patches:" << nl << nl
685 << "patchID\tpatch\ttype" << nl
686 << "-------\t-----\t----"
687 << endl;
688
689 interRegionTopPatch.setSize(zoneNames.size(), -1);
690 interRegionBottomPatch.setSize(zoneNames.size(), -1);
691
692 label nOldPatches = newPatches.size();
693 forAll(zoneNames, zoneI)
694 {
695 word interName
696 (
698 +"_to_"
699 +shellRegionName
700 +'_'
701 +zoneNames[zoneI]
702 );
703
704 if (isInternal[zoneI])
705 {
706 interRegionTopPatch[zoneI] = addPatch<mappedWallPolyPatch>
707 (
708 mesh.boundaryMesh(),
709 interName + "_top",
710 newPatches
711 );
712 Pout<< interRegionTopPatch[zoneI]
713 << '\t' << newPatches[interRegionTopPatch[zoneI]]->name()
714 << '\t' << newPatches[interRegionTopPatch[zoneI]]->type()
715 << nl;
716
717 interRegionBottomPatch[zoneI] = addPatch<mappedWallPolyPatch>
718 (
719 mesh.boundaryMesh(),
720 interName + "_bottom",
721 newPatches
722 );
723 Pout<< interRegionBottomPatch[zoneI]
724 << '\t' << newPatches[interRegionBottomPatch[zoneI]]->name()
725 << '\t' << newPatches[interRegionBottomPatch[zoneI]]->type()
726 << nl;
727 }
728 else if (zoneShadowNames.size() == 0)
729 {
730 interRegionTopPatch[zoneI] = addPatch<polyPatch>
731 (
732 mesh.boundaryMesh(),
733 zoneNames[zoneI] + "_top",
734 newPatches
735 );
736 Pout<< interRegionTopPatch[zoneI]
737 << '\t' << newPatches[interRegionTopPatch[zoneI]]->name()
738 << '\t' << newPatches[interRegionTopPatch[zoneI]]->type()
739 << nl;
740
741 interRegionBottomPatch[zoneI] = addPatch<mappedWallPolyPatch>
742 (
743 mesh.boundaryMesh(),
744 interName,
745 newPatches
746 );
747 Pout<< interRegionBottomPatch[zoneI]
748 << '\t' << newPatches[interRegionBottomPatch[zoneI]]->name()
749 << '\t' << newPatches[interRegionBottomPatch[zoneI]]->type()
750 << nl;
751 }
752 else //patch using shadow face zones.
753 {
754 interRegionTopPatch[zoneI] = addPatch<mappedWallPolyPatch>
755 (
756 mesh.boundaryMesh(),
757 zoneShadowNames[zoneI] + "_top",
758 newPatches
759 );
760 Pout<< interRegionTopPatch[zoneI]
761 << '\t' << newPatches[interRegionTopPatch[zoneI]]->name()
762 << '\t' << newPatches[interRegionTopPatch[zoneI]]->type()
763 << nl;
764
765 interRegionBottomPatch[zoneI] = addPatch<mappedWallPolyPatch>
766 (
767 mesh.boundaryMesh(),
768 interName,
769 newPatches
770 );
771 Pout<< interRegionBottomPatch[zoneI]
772 << '\t' << newPatches[interRegionBottomPatch[zoneI]]->name()
773 << '\t' << newPatches[interRegionBottomPatch[zoneI]]->type()
774 << nl;
775 }
776 }
777 Pout<< "Added " << newPatches.size()-nOldPatches
778 << " inter-region patches." << nl
779 << endl;
780}
781
782
783// Sets sidePatch[edgeI] to interprocessor or cyclic patch. Adds any
784// coupled patches if necessary.
785void addCoupledPatches
786(
787 const fvMesh& mesh,
788 const primitiveFacePatch& extrudePatch,
789 const labelList& extrudeMeshFaces,
790 const labelList& extrudeMeshEdges,
791 const mapDistribute& extrudeEdgeFacesMap,
792 const labelListList& extrudeEdgeGlobalFaces,
793
794 labelList& sidePatchID,
795 DynamicList<polyPatch*>& newPatches
796)
797{
798 // Calculate opposite processor for coupled edges (only if shared by
799 // two procs). Note: could have saved original globalEdgeFaces structure.
800
801 // Get procID in extrudeEdgeGlobalFaces order
802 labelList procID(extrudeEdgeGlobalFaces.size(), Pstream::myProcNo());
803 extrudeEdgeFacesMap.distribute(procID);
804
805 labelList minProcID(extrudeEdgeGlobalFaces.size(), labelMax);
806 labelList maxProcID(extrudeEdgeGlobalFaces.size(), labelMin);
807
808 forAll(extrudeEdgeGlobalFaces, edgeI)
809 {
810 const labelList& eFaces = extrudeEdgeGlobalFaces[edgeI];
811 if (eFaces.size())
812 {
813 forAll(eFaces, i)
814 {
815 label proci = procID[eFaces[i]];
816 minProcID[edgeI] = min(minProcID[edgeI], proci);
817 maxProcID[edgeI] = max(maxProcID[edgeI], proci);
818 }
819 }
820 }
821 syncTools::syncEdgeList
822 (
823 mesh,
824 extrudeMeshEdges,
825 minProcID,
827 labelMax // null value
828 );
829 syncTools::syncEdgeList
830 (
831 mesh,
832 extrudeMeshEdges,
833 maxProcID,
835 labelMin // null value
836 );
837
838 Pout<< "Adding processor or cyclic patches:" << nl << nl
839 << "patchID\tpatch" << nl
840 << "-------\t-----"
841 << endl;
842
843 label nOldPatches = newPatches.size();
844
845 sidePatchID.setSize(extrudePatch.edgeFaces().size(), -1);
846 forAll(extrudePatch.edgeFaces(), edgeI)
847 {
848 const labelList& eFaces = extrudePatch.edgeFaces()[edgeI];
849
850 if
851 (
852 eFaces.size() == 1
853 && extrudeEdgeGlobalFaces[edgeI].size() == 2
854 )
855 {
856 // coupled boundary edge. Find matching patch.
857 label nbrProci = minProcID[edgeI];
858 if (nbrProci == Pstream::myProcNo())
859 {
860 nbrProci = maxProcID[edgeI];
861 }
862
863
864 if (nbrProci == Pstream::myProcNo())
865 {
866 // Cyclic patch since both procs the same. This cyclic should
867 // already exist in newPatches so no adding necessary.
868
869 label facei = findUncoveredCyclicPatchFace
870 (
871 mesh,
872 labelUIndList(extrudeMeshFaces, eFaces),
873 extrudeMeshEdges[edgeI]
874 );
875
876 if (facei != -1)
877 {
878 const polyBoundaryMesh& patches = mesh.boundaryMesh();
879
880 label newPatchi = findPatchID
881 (
882 newPatches,
883 patches[patches.whichPatch(facei)].name()
884 );
885
886 sidePatchID[edgeI] = newPatchi;
887 }
888 else
889 {
891 << "Unable to determine coupled patch addressing"
892 << abort(FatalError);
893 }
894 }
895 else
896 {
897 // Processor patch
898 word name
899 (
900 processorPolyPatch::newName(Pstream::myProcNo(), nbrProci)
901 );
902
903 sidePatchID[edgeI] = findPatchID(newPatches, name);
904
905 if (sidePatchID[edgeI] == -1)
906 {
907 dictionary patchDict;
908 patchDict.add("myProcNo", Pstream::myProcNo());
909 patchDict.add("neighbProcNo", nbrProci);
910
911 sidePatchID[edgeI] = addPatch<processorPolyPatch>
912 (
913 mesh.boundaryMesh(),
914 name,
915 patchDict,
916 newPatches
917 );
918
919 Pout<< sidePatchID[edgeI] << '\t' << name
920 << nl;
921 }
922 }
923 }
924 }
925 Pout<< "Added " << newPatches.size()-nOldPatches
926 << " coupled patches." << nl
927 << endl;
928}
929
930
931void addZoneSidePatches
932(
933 const fvMesh& mesh,
934 const wordList& zoneNames,
935 const word& oneDPolyPatchType,
936
937 DynamicList<polyPatch*>& newPatches,
938 labelList& zoneSidePatch
939)
940{
941 Pout<< "Adding patches for sides on zones:" << nl << nl
942 << "patchID\tpatch" << nl
943 << "-------\t-----"
944 << endl;
945
946 label nOldPatches = newPatches.size();
947
948 forAll(zoneSidePatch, zoneI)
949 {
950 if (oneDPolyPatchType != word::null)
951 {
952 // Reuse single empty patch.
953 word patchName;
954 if (oneDPolyPatchType == "empty")
955 {
956 patchName = "oneDEmptyPatch";
957 zoneSidePatch[zoneI] = addPatch<emptyPolyPatch>
958 (
959 mesh.boundaryMesh(),
960 patchName,
961 newPatches
962 );
963 }
964 else if (oneDPolyPatchType == "wedge")
965 {
966 patchName = "oneDWedgePatch";
967 zoneSidePatch[zoneI] = addPatch<wedgePolyPatch>
968 (
969 mesh.boundaryMesh(),
970 patchName,
971 newPatches
972 );
973 }
974 else
975 {
977 << "Type " << oneDPolyPatchType << " does not exist "
978 << exit(FatalError);
979 }
980
981 Pout<< zoneSidePatch[zoneI] << '\t' << patchName << nl;
982 }
983 else if (zoneSidePatch[zoneI] > 0)
984 {
985 word patchName = zoneNames[zoneI] + "_" + "side";
986
987 zoneSidePatch[zoneI] = addPatch<polyPatch>
988 (
989 mesh.boundaryMesh(),
990 patchName,
991 newPatches
992 );
993
994 Pout<< zoneSidePatch[zoneI] << '\t' << patchName << nl;
995 }
996 }
997 Pout<< "Added " << newPatches.size()-nOldPatches << " zone-side patches."
998 << nl << endl;
999}
1000
1001
1002void addInterZonePatches
1003(
1004 const fvMesh& mesh,
1005 const wordList& zoneNames,
1006 const bool oneD,
1007
1008 labelList& zoneZonePatch_min,
1009 labelList& zoneZonePatch_max,
1010 DynamicList<polyPatch*>& newPatches
1011)
1012{
1013 Pout<< "Adding inter-zone patches:" << nl << nl
1014 << "patchID\tpatch" << nl
1015 << "-------\t-----"
1016 << endl;
1017
1018 dictionary transformDict;
1019 transformDict.add
1020 (
1021 "transform",
1022 cyclicPolyPatch::transformTypeNames[cyclicPolyPatch::NOORDERING]
1023 );
1024
1025 label nOldPatches = newPatches.size();
1026
1027 if (!oneD)
1028 {
1029 forAll(zoneZonePatch_min, minZone)
1030 {
1031 for (label maxZone = minZone; maxZone < zoneNames.size(); maxZone++)
1032 {
1033 label index = minZone*zoneNames.size()+maxZone;
1034
1035 if (zoneZonePatch_min[index] > 0)
1036 {
1037 word minToMax =
1038 zoneNames[minZone]
1039 + "_to_"
1040 + zoneNames[maxZone];
1041 word maxToMin =
1042 zoneNames[maxZone]
1043 + "_to_"
1044 + zoneNames[minZone];
1045
1046 {
1047 transformDict.set("neighbourPatch", maxToMin);
1048 zoneZonePatch_min[index] =
1049 addPatch<nonuniformTransformCyclicPolyPatch>
1050 (
1051 mesh.boundaryMesh(),
1052 minToMax,
1053 transformDict,
1054 newPatches
1055 );
1056 Pout<< zoneZonePatch_min[index] << '\t' << minToMax
1057 << nl;
1058 }
1059 {
1060 transformDict.set("neighbourPatch", minToMax);
1061 zoneZonePatch_max[index] =
1062 addPatch<nonuniformTransformCyclicPolyPatch>
1063 (
1064 mesh.boundaryMesh(),
1065 maxToMin,
1066 transformDict,
1067 newPatches
1068 );
1069 Pout<< zoneZonePatch_max[index] << '\t' << maxToMin
1070 << nl;
1071 }
1072
1073 }
1074 }
1075 }
1076 }
1077 Pout<< "Added " << newPatches.size()-nOldPatches << " inter-zone patches."
1078 << nl << endl;
1079}
1080
1081
1082tmp<pointField> calcOffset
1083(
1084 const primitiveFacePatch& extrudePatch,
1085 const createShellMesh& extruder,
1086 const polyPatch& pp
1087)
1088{
1090
1091 tmp<pointField> toffsets(new pointField(fc.size()));
1092 pointField& offsets = toffsets.ref();
1093
1094 forAll(fc, i)
1095 {
1096 label meshFacei = pp.start()+i;
1097 label patchFacei = mag(extruder.faceToFaceMap()[meshFacei])-1;
1098 point patchFc = extrudePatch[patchFacei].centre
1099 (
1100 extrudePatch.points()
1101 );
1102 offsets[i] = patchFc - fc[i];
1103 }
1104 return toffsets;
1105}
1106
1107
1108void setCouplingInfo
1109(
1110 fvMesh& mesh,
1111 const labelList& zoneToPatch,
1112 const word& sampleRegion,
1114 const List<pointField>& offsets
1115)
1116{
1117 const polyBoundaryMesh& patches = mesh.boundaryMesh();
1118
1119 List<polyPatch*> newPatches
1120 (
1121 patches.size(),
1122 static_cast<polyPatch*>(nullptr)
1123 );
1124
1125 forAll(zoneToPatch, zoneI)
1126 {
1127 const label patchi = zoneToPatch[zoneI];
1128
1129 if (patchi != -1)
1130 {
1131 const polyPatch& pp = patches[patchi];
1132
1133 if (isA<mappedWallPolyPatch>(pp))
1134 {
1135 const boundBox bb(pp.points(), pp.meshPoints(), true);
1136 const vector avgOffset = gAverage(offsets[zoneI]);
1137 const scalar mergeSqrDist =
1138 gMax(magSqr(offsets[zoneI]-avgOffset));
1139
1140 // Create with uniform offset initially
1141 auto mappedPtr = autoPtr<mappedWallPolyPatch>::New
1142 (
1143 pp.name(),
1144 pp.size(),
1145 pp.start(),
1146 patchi,
1147 sampleRegion, // sampleRegion
1148 mode, // sampleMode
1149 pp.name(), // samplePatch
1150
1151 avgOffset, // uniform offset
1152 patches
1153 );
1154
1155 Info<< "Adding on " << mesh.name() << " coupling patch "
1156 << pp.name() << " with ";
1157
1158 // Verify uniformity of offset
1159 // (same check as blockMesh geom merge)
1160 if (mergeSqrDist < magSqr(10*SMALL*bb.span()))
1161 {
1162 Info<< "uniform offset " << avgOffset << endl;
1163 }
1164 else
1165 {
1166 Info<< "non-uniform offset" << endl;
1167
1168 (*mappedPtr).setOffset(offsets[zoneI]);
1169 }
1170
1171 newPatches[patchi] = mappedPtr.release();
1172 }
1173 }
1174 }
1175
1176 forAll(newPatches, patchi)
1177 {
1178 if (!newPatches[patchi])
1179 {
1180 newPatches[patchi] = patches[patchi].clone(patches).ptr();
1181 //newPatches[patchi].index() = patchi;
1182 }
1183 }
1184
1185
1186 #ifdef FULLDEBUG
1187 Pout<< "*** setCouplingInfo addFvPAtches:" << nl;
1188 printPatches(Pout, newPatches);
1189 Pout<< "*** setCouplingInfo end of addFvPAtches:" << endl;
1190 #endif
1191
1192 mesh.removeFvBoundary();
1193 mesh.addFvPatches(newPatches, true);
1194}
1195
1196
1197// Extrude and write geometric properties
1198void extrudeGeometricProperties
1199(
1200 const polyMesh& mesh,
1201 const primitiveFacePatch& extrudePatch,
1202 const createShellMesh& extruder,
1203 const polyMesh& regionMesh,
1204 const extrudeModel& model
1205)
1206{
1207 const pointIOField patchFaceCentres
1208 (
1209 IOobject
1210 (
1211 "patchFaceCentres",
1212 mesh.pointsInstance(),
1213 mesh.meshSubDir,
1214 mesh,
1215 IOobject::MUST_READ
1216 )
1217 );
1218
1219 const pointIOField patchEdgeCentres
1220 (
1221 IOobject
1222 (
1223 "patchEdgeCentres",
1224 mesh.pointsInstance(),
1225 mesh.meshSubDir,
1226 mesh,
1227 IOobject::MUST_READ
1228 )
1229 );
1230
1231 //forAll(extrudePatch.edges(), edgeI)
1232 //{
1233 // const edge& e = extrudePatch.edges()[edgeI];
1234 // Pout<< "Edge:" << e.centre(extrudePatch.localPoints()) << nl
1235 // << "read:" << patchEdgeCentres[edgeI]
1236 // << endl;
1237 //}
1238
1239
1240 // Determine edge normals on original patch
1241 labelList patchEdges;
1242 labelList coupledEdges;
1243 bitSet sameEdgeOrientation;
1244 PatchTools::matchEdges
1245 (
1246 extrudePatch,
1247 mesh.globalData().coupledPatch(),
1248 patchEdges,
1249 coupledEdges,
1250 sameEdgeOrientation
1251 );
1252
1253 pointField patchEdgeNormals
1254 (
1255 PatchTools::edgeNormals
1256 (
1257 mesh,
1258 extrudePatch,
1259 patchEdges,
1260 coupledEdges
1261 )
1262 );
1263
1264
1265 pointIOField faceCentres
1266 (
1267 IOobject
1268 (
1269 "faceCentres",
1270 regionMesh.pointsInstance(),
1271 regionMesh.meshSubDir,
1272 regionMesh,
1273 IOobject::NO_READ,
1274 IOobject::NO_WRITE,
1275 false
1276 ),
1277 regionMesh.nFaces()
1278 );
1279
1280
1281 // Work out layers. Guaranteed in columns so no fancy parallel bits.
1282
1283
1284 forAll(extruder.faceToFaceMap(), facei)
1285 {
1286 if (extruder.faceToFaceMap()[facei] != 0)
1287 {
1288 // 'horizontal' face
1289 label patchFacei = mag(extruder.faceToFaceMap()[facei])-1;
1290
1291 label celli = regionMesh.faceOwner()[facei];
1292 if (regionMesh.isInternalFace(facei))
1293 {
1294 celli = max(celli, regionMesh.faceNeighbour()[facei]);
1295 }
1296
1297 // Calculate layer from cell numbering (see createShellMesh)
1298 label layerI = (celli % model.nLayers());
1299
1300 if
1301 (
1302 !regionMesh.isInternalFace(facei)
1303 && extruder.faceToFaceMap()[facei] > 0
1304 )
1305 {
1306 // Top face
1307 layerI++;
1308 }
1309
1310
1311 // Recalculate based on extrusion model
1312 faceCentres[facei] = model
1313 (
1314 patchFaceCentres[patchFacei],
1315 extrudePatch.faceNormals()[patchFacei],
1316 layerI
1317 );
1318 }
1319 else
1320 {
1321 // 'vertical face
1322 label patchEdgeI = extruder.faceToEdgeMap()[facei];
1323 label layerI =
1324 (
1325 regionMesh.faceOwner()[facei]
1326 % model.nLayers()
1327 );
1328
1329 // Extrude patch edge centre to this layer
1330 point pt0 = model
1331 (
1332 patchEdgeCentres[patchEdgeI],
1333 patchEdgeNormals[patchEdgeI],
1334 layerI
1335 );
1336 // Extrude patch edge centre to next layer
1337 point pt1 = model
1338 (
1339 patchEdgeCentres[patchEdgeI],
1340 patchEdgeNormals[patchEdgeI],
1341 layerI+1
1342 );
1343
1344 // Interpolate
1345 faceCentres[facei] = 0.5*(pt0+pt1);
1346 }
1347 }
1348
1349 pointIOField cellCentres
1350 (
1351 IOobject
1352 (
1353 "cellCentres",
1354 regionMesh.pointsInstance(),
1355 regionMesh.meshSubDir,
1356 regionMesh,
1357 IOobject::NO_READ,
1358 IOobject::NO_WRITE,
1359 false
1360 ),
1361 regionMesh.nCells()
1362 );
1363
1364 forAll(extruder.cellToFaceMap(), celli)
1365 {
1366 label patchFacei = extruder.cellToFaceMap()[celli];
1367
1368 // Calculate layer from cell numbering (see createShellMesh)
1369 label layerI = (celli % model.nLayers());
1370
1371 // Recalculate based on extrusion model
1372 point pt0 = model
1373 (
1374 patchFaceCentres[patchFacei],
1375 extrudePatch.faceNormals()[patchFacei],
1376 layerI
1377 );
1378 point pt1 = model
1379 (
1380 patchFaceCentres[patchFacei],
1381 extrudePatch.faceNormals()[patchFacei],
1382 layerI+1
1383 );
1384
1385 // Interpolate
1386 cellCentres[celli] = 0.5*(pt0+pt1);
1387 }
1388
1389
1390 // Bit of checking
1391 if (false)
1392 {
1393 OBJstream faceStr(regionMesh.time().path()/"faceCentres.obj");
1394 OBJstream cellStr(regionMesh.time().path()/"cellCentres.obj");
1395
1396 forAll(faceCentres, facei)
1397 {
1398 Pout<< "Model :" << faceCentres[facei] << endl
1399 << "regionMesh:" << regionMesh.faceCentres()[facei] << endl;
1400 faceStr.write
1401 (
1403 (
1404 faceCentres[facei],
1405 regionMesh.faceCentres()[facei]
1406 )
1407 );
1408 }
1409 forAll(cellCentres, celli)
1410 {
1411 Pout<< "Model :" << cellCentres[celli] << endl
1412 << "regionMesh:" << regionMesh.cellCentres()[celli] << endl;
1413 cellStr.write
1414 (
1416 (
1417 cellCentres[celli],
1418 regionMesh.cellCentres()[celli]
1419 )
1420 );
1421 }
1422 }
1423
1424
1425
1426 Info<< "Writing geometric properties for mesh " << regionMesh.name()
1427 << " to " << regionMesh.pointsInstance() << nl
1428 << endl;
1429
1430 bool ok = faceCentres.write() && cellCentres.write();
1431
1432 if (!ok)
1433 {
1435 << "Failed writing " << faceCentres.objectPath()
1436 << " and " << cellCentres.objectPath()
1437 << exit(FatalError);
1438 }
1439}
1440
1441
1442int main(int argc, char *argv[])
1443{
1444 argList::addNote
1445 (
1446 "Create region mesh by extruding a faceZone or faceSet"
1447 );
1448
1449 #include "addRegionOption.H"
1450 #include "addOverwriteOption.H"
1451
1452 argList::addOption
1453 (
1454 "dict", "file", "Alternative extrudeToRegionMeshDict"
1455 );
1456
1457 #include "setRootCase.H"
1458 #include "createTime.H"
1459 #include "createNamedMesh.H"
1460
1461 if (mesh.boundaryMesh().checkParallelSync(true))
1462 {
1463 List<wordList> allNames(Pstream::nProcs());
1464 allNames[Pstream::myProcNo()] = mesh.boundaryMesh().names();
1465 Pstream::allGatherList(allNames);
1466
1468 << "Patches are not synchronised on all processors."
1469 << " Per processor patches " << allNames
1470 << exit(FatalError);
1471 }
1472
1473
1474 const word oldInstance = mesh.pointsInstance();
1475 const bool overwrite = args.found("overwrite");
1476
1477 const word dictName("extrudeToRegionMeshDict");
1478
1480
1482
1483 // Point generator
1484 autoPtr<extrudeModel> model(extrudeModel::New(dict));
1485
1486 // Region
1487 const word shellRegionName(dict.get<word>("region"));
1488
1489 // Faces to extrude - either faceZones or faceSets (boundary faces only)
1490 wordList zoneNames;
1491 wordList zoneShadowNames;
1492
1493 const bool hasZones = dict.found("faceZones");
1494 if (hasZones)
1495 {
1496 dict.readEntry("faceZones", zoneNames);
1497 dict.readIfPresent("faceZonesShadow", zoneShadowNames);
1498
1499 // Check
1500 if (dict.found("faceSets"))
1501 {
1503 << "Please supply faces to extrude either through 'faceZones'"
1504 << " or 'faceSets' entry. Found both."
1505 << exit(FatalIOError);
1506 }
1507 }
1508 else
1509 {
1510 dict.readEntry("faceSets", zoneNames);
1511 dict.readIfPresent("faceSetsShadow", zoneShadowNames);
1512 }
1513
1514
1515 mappedPatchBase::sampleMode sampleMode =
1516 mappedPatchBase::sampleModeNames_.get("sampleMode", dict);
1517
1518 const bool oneD(dict.get<bool>("oneD"));
1519 bool oneDNonManifoldEdges(false);
1520 word oneDPatchType(emptyPolyPatch::typeName);
1521 if (oneD)
1522 {
1523 oneDNonManifoldEdges = dict.getOrDefault("nonManifold", false);
1524 oneDPatchType = dict.get<word>("oneDPolyPatchType");
1525 }
1526
1527 const bool adaptMesh(dict.get<bool>("adaptMesh"));
1528
1529 if (hasZones)
1530 {
1531 Info<< "Extruding zones " << zoneNames
1532 << " on mesh " << regionName
1533 << " into shell mesh " << shellRegionName
1534 << endl;
1535 }
1536 else
1537 {
1538 Info<< "Extruding faceSets " << zoneNames
1539 << " on mesh " << regionName
1540 << " into shell mesh " << shellRegionName
1541 << endl;
1542 }
1543
1544 if (shellRegionName == regionName)
1545 {
1547 << "Cannot extrude into same region as mesh." << endl
1548 << "Mesh region : " << regionName << endl
1549 << "Shell region : " << shellRegionName
1550 << exit(FatalIOError);
1551 }
1552
1553
1554 if (oneD)
1555 {
1556 if (oneDNonManifoldEdges)
1557 {
1558 Info<< "Extruding as 1D columns with sides in patch type "
1559 << oneDPatchType
1560 << " and connected points (except on non-manifold areas)."
1561 << endl;
1562 }
1563 else
1564 {
1565 Info<< "Extruding as 1D columns with sides in patch type "
1566 << oneDPatchType
1567 << " and duplicated points (overlapping volumes)."
1568 << endl;
1569 }
1570 }
1571
1572
1573
1574
1576 //IOobjectList objects(mesh, runTime.timeName());
1577 //
1579 //
1580 //PtrList<volScalarField> vsFlds;
1581 //ReadFields(mesh, objects, vsFlds);
1582 //
1583 //PtrList<volVectorField> vvFlds;
1584 //ReadFields(mesh, objects, vvFlds);
1585 //
1586 //PtrList<volSphericalTensorField> vstFlds;
1587 //ReadFields(mesh, objects, vstFlds);
1588 //
1589 //PtrList<volSymmTensorField> vsymtFlds;
1590 //ReadFields(mesh, objects, vsymtFlds);
1591 //
1592 //PtrList<volTensorField> vtFlds;
1593 //ReadFields(mesh, objects, vtFlds);
1594 //
1596 //
1597 //PtrList<surfaceScalarField> ssFlds;
1598 //ReadFields(mesh, objects, ssFlds);
1599 //
1600 //PtrList<surfaceVectorField> svFlds;
1601 //ReadFields(mesh, objects, svFlds);
1602 //
1603 //PtrList<surfaceSphericalTensorField> sstFlds;
1604 //ReadFields(mesh, objects, sstFlds);
1605 //
1606 //PtrList<surfaceSymmTensorField> ssymtFlds;
1607 //ReadFields(mesh, objects, ssymtFlds);
1608 //
1609 //PtrList<surfaceTensorField> stFlds;
1610 //ReadFields(mesh, objects, stFlds);
1611 //
1613 //
1614 //PtrList<pointScalarField> psFlds;
1615 //ReadFields(pointMesh::New(mesh), objects, psFlds);
1616 //
1617 //PtrList<pointVectorField> pvFlds;
1618 //ReadFields(pointMesh::New(mesh), objects, pvFlds);
1619
1620
1621
1622 // Create dummy fv* files
1623 fvMeshTools::createDummyFvMeshFiles(mesh, shellRegionName, true);
1624
1625
1626 word meshInstance;
1627 if (!overwrite)
1628 {
1629 ++runTime;
1630 meshInstance = runTime.timeName();
1631 }
1632 else
1633 {
1634 meshInstance = oldInstance;
1635 }
1636 Info<< "Writing meshes to " << meshInstance << nl << endl;
1637
1638
1639 const polyBoundaryMesh& patches = mesh.boundaryMesh();
1640
1641
1642 // Extract faces to extrude
1643 // ~~~~~~~~~~~~~~~~~~~~~~~~
1644 // Note: zoneID are regions of extrusion. They are not mesh.faceZones
1645 // indices.
1646
1647 // From extrude zone to mesh zone (or -1 if extruding faceSets)
1648 labelList meshZoneID;
1649 // Per extrude zone whether contains internal or external faces
1650 boolList isInternal(zoneNames.size(), false);
1651
1652 labelList extrudeMeshFaces;
1653 faceList zoneFaces;
1655 boolList zoneFlipMap;
1656 // Shadow
1657 labelList zoneShadowIDs; // from extrude shadow zone to mesh zone
1658 labelList extrudeMeshShadowFaces;
1659 boolList zoneShadowFlipMap;
1660 labelList zoneShadowID;
1661
1662 if (hasZones)
1663 {
1664 const faceZoneMesh& faceZones = mesh.faceZones();
1665
1666 meshZoneID.setSize(zoneNames.size());
1667 forAll(zoneNames, i)
1668 {
1669 meshZoneID[i] = faceZones.findZoneID(zoneNames[i]);
1670 if (meshZoneID[i] == -1)
1671 {
1673 << "Cannot find zone " << zoneNames[i] << endl
1674 << "Valid zones are " << faceZones.names()
1675 << exit(FatalIOError);
1676 }
1677 }
1678 // Collect per face information
1679 label nExtrudeFaces = 0;
1680 forAll(meshZoneID, i)
1681 {
1682 nExtrudeFaces += faceZones[meshZoneID[i]].size();
1683 }
1684 extrudeMeshFaces.setSize(nExtrudeFaces);
1685 zoneFaces.setSize(nExtrudeFaces);
1686 zoneID.setSize(nExtrudeFaces);
1687 zoneFlipMap.setSize(nExtrudeFaces);
1688 nExtrudeFaces = 0;
1689 forAll(meshZoneID, i)
1690 {
1691 const faceZone& fz = faceZones[meshZoneID[i]];
1692 const primitiveFacePatch& fzp = fz();
1693 forAll(fz, j)
1694 {
1695 extrudeMeshFaces[nExtrudeFaces] = fz[j];
1696 zoneFaces[nExtrudeFaces] = fzp[j];
1697 zoneID[nExtrudeFaces] = i;
1698 zoneFlipMap[nExtrudeFaces] = fz.flipMap()[j];
1699 nExtrudeFaces++;
1700
1701 if (mesh.isInternalFace(fz[j]))
1702 {
1703 isInternal[i] = true;
1704 }
1705 }
1706 }
1707
1708 // Shadow zone
1709 // ~~~~~~~~~~~
1710
1711 if (zoneShadowNames.size())
1712 {
1713 zoneShadowIDs.setSize(zoneShadowNames.size());
1714 forAll(zoneShadowNames, i)
1715 {
1716 zoneShadowIDs[i] = faceZones.findZoneID(zoneShadowNames[i]);
1717 if (zoneShadowIDs[i] == -1)
1718 {
1720 << "Cannot find zone " << zoneShadowNames[i] << endl
1721 << "Valid zones are " << faceZones.names()
1722 << exit(FatalIOError);
1723 }
1724 }
1725
1726 label nShadowFaces = 0;
1727 forAll(zoneShadowIDs, i)
1728 {
1729 nShadowFaces += faceZones[zoneShadowIDs[i]].size();
1730 }
1731
1732 extrudeMeshShadowFaces.setSize(nShadowFaces);
1733 zoneShadowFlipMap.setSize(nShadowFaces);
1734 zoneShadowID.setSize(nShadowFaces);
1735
1736 nShadowFaces = 0;
1737 forAll(zoneShadowIDs, i)
1738 {
1739 const faceZone& fz = faceZones[zoneShadowIDs[i]];
1740 forAll(fz, j)
1741 {
1742 extrudeMeshShadowFaces[nShadowFaces] = fz[j];
1743 zoneShadowFlipMap[nShadowFaces] = fz.flipMap()[j];
1744 zoneShadowID[nShadowFaces] = i;
1745 nShadowFaces++;
1746 }
1747 }
1748 }
1749 }
1750 else
1751 {
1752 meshZoneID.setSize(zoneNames.size(), -1);
1753 // Load faceSets
1754 PtrList<faceSet> zones(zoneNames.size());
1755 forAll(zoneNames, i)
1756 {
1757 Info<< "Loading faceSet " << zoneNames[i] << endl;
1758 zones.set(i, new faceSet(mesh, zoneNames[i]));
1759 }
1760
1761
1762 // Collect per face information
1763 label nExtrudeFaces = 0;
1764 forAll(zones, i)
1765 {
1766 nExtrudeFaces += zones[i].size();
1767 }
1768 extrudeMeshFaces.setSize(nExtrudeFaces);
1769 zoneFaces.setSize(nExtrudeFaces);
1770 zoneID.setSize(nExtrudeFaces);
1771 zoneFlipMap.setSize(nExtrudeFaces);
1772
1773 nExtrudeFaces = 0;
1774 forAll(zones, i)
1775 {
1776 const faceSet& fz = zones[i];
1777 for (const label facei : fz)
1778 {
1779 if (mesh.isInternalFace(facei))
1780 {
1782 << "faceSet " << fz.name()
1783 << "contains internal faces."
1784 << " This is not permitted."
1785 << exit(FatalIOError);
1786 }
1787 extrudeMeshFaces[nExtrudeFaces] = facei;
1788 zoneFaces[nExtrudeFaces] = mesh.faces()[facei];
1789 zoneID[nExtrudeFaces] = i;
1790 zoneFlipMap[nExtrudeFaces] = false;
1791 nExtrudeFaces++;
1792
1793 if (mesh.isInternalFace(facei))
1794 {
1795 isInternal[i] = true;
1796 }
1797 }
1798 }
1799
1800
1801 // Shadow zone
1802 // ~~~~~~~~~~~
1803
1804 PtrList<faceSet> shadowZones(zoneShadowNames.size());
1805 if (zoneShadowNames.size())
1806 {
1807 zoneShadowIDs.setSize(zoneShadowNames.size(), -1);
1808 forAll(zoneShadowNames, i)
1809 {
1810 shadowZones.set(i, new faceSet(mesh, zoneShadowNames[i]));
1811 }
1812
1813 label nShadowFaces = 0;
1814 for (const faceSet& fz : shadowZones)
1815 {
1816 nShadowFaces += fz.size();
1817 }
1818
1819 if (nExtrudeFaces != nShadowFaces)
1820 {
1822 << "Extruded faces " << nExtrudeFaces << endl
1823 << "is different from shadow faces. " << nShadowFaces
1824 << "This is not permitted " << endl
1825 << exit(FatalIOError);
1826 }
1827
1828 extrudeMeshShadowFaces.setSize(nShadowFaces);
1829 zoneShadowFlipMap.setSize(nShadowFaces);
1830 zoneShadowID.setSize(nShadowFaces);
1831
1832 nShadowFaces = 0;
1833 forAll(shadowZones, i)
1834 {
1835 const faceSet& fz = shadowZones[i];
1836 for (const label facei : fz)
1837 {
1838 if (mesh.isInternalFace(facei))
1839 {
1841 << "faceSet " << fz.name()
1842 << "contains internal faces."
1843 << " This is not permitted."
1844 << exit(FatalIOError);
1845 }
1846 extrudeMeshShadowFaces[nShadowFaces] = facei;
1847 zoneShadowFlipMap[nShadowFaces] = false;
1848 zoneShadowID[nShadowFaces] = i;
1849 nShadowFaces++;
1850 }
1851 }
1852 }
1853 }
1854 const primitiveFacePatch extrudePatch(std::move(zoneFaces), mesh.points());
1855
1856
1857 Pstream::listCombineAllGather(isInternal, orEqOp<bool>());
1858
1859 // Check zone either all internal or all external faces
1860 checkZoneInside(mesh, zoneNames, zoneID, extrudeMeshFaces, isInternal);
1861
1862
1863 const pointField& extrudePoints = extrudePatch.localPoints();
1864 const faceList& extrudeFaces = extrudePatch.localFaces();
1865 const labelListList& edgeFaces = extrudePatch.edgeFaces();
1866
1867
1868 Info<< "extrudePatch :"
1869 << " faces:" << extrudePatch.size()
1870 << " points:" << extrudePatch.nPoints()
1871 << " edges:" << extrudePatch.nEdges()
1872 << nl
1873 << endl;
1874
1875
1876 // Determine per-extrude-edge info
1877 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1878
1879 // Corresponding mesh edges
1880 const labelList extrudeMeshEdges
1881 (
1882 extrudePatch.meshEdges
1883 (
1884 mesh.edges(),
1885 mesh.pointEdges()
1886 )
1887 );
1888
1889 const globalIndex globalExtrudeFaces(extrudePatch.size());
1890
1891 // Global pp faces per pp edge.
1892 labelListList extrudeEdgeGlobalFaces
1893 (
1894 globalEdgeFaces
1895 (
1896 mesh,
1897 globalExtrudeFaces,
1898 extrudePatch,
1899 extrudeMeshEdges
1900 )
1901 );
1902 List<Map<label>> compactMap;
1903 const mapDistribute extrudeEdgeFacesMap
1904 (
1905 globalExtrudeFaces,
1906 extrudeEdgeGlobalFaces,
1907 compactMap
1908 );
1909
1910
1911 // Determine min and max zone per edge
1912 labelList edgeMinZoneID;
1913 labelList edgeMaxZoneID;
1914 calcEdgeMinMaxZone
1915 (
1916 mesh,
1917 extrudePatch,
1918 extrudeMeshEdges,
1919 zoneID,
1920 extrudeEdgeFacesMap,
1921 extrudeEdgeGlobalFaces,
1922
1923 edgeMinZoneID,
1924 edgeMaxZoneID
1925 );
1926
1927
1928
1929
1930 DynamicList<polyPatch*> regionPatches(patches.size());
1931 // Copy all non-local patches since these are used on boundary edges of
1932 // the extrusion
1933 forAll(patches, patchi)
1934 {
1935 if (!isA<processorPolyPatch>(patches[patchi]))
1936 {
1937 label newPatchi = regionPatches.size();
1938 regionPatches.append
1939 (
1940 patches[patchi].clone
1941 (
1942 patches,
1943 newPatchi,
1944 0, // size
1945 0 // start
1946 ).ptr()
1947 );
1948 }
1949 }
1950
1951
1952 // Add interface patches
1953 // ~~~~~~~~~~~~~~~~~~~~~
1954
1955 // From zone to interface patch (region side)
1956 labelList interRegionTopPatch;
1957 labelList interRegionBottomPatch;
1958
1959 addCouplingPatches
1960 (
1961 mesh,
1962 regionName,
1963 shellRegionName,
1964 zoneNames,
1965 zoneShadowNames,
1966 isInternal,
1967 meshZoneID,
1968
1969 regionPatches,
1970 interRegionTopPatch,
1971 interRegionBottomPatch
1972 );
1973
1974
1975 // From zone to interface patch (mesh side)
1976 labelList interMeshTopPatch;
1977 labelList interMeshBottomPatch;
1978
1979 if (adaptMesh)
1980 {
1981 // Add coupling patches to mesh
1982
1983 // 1. Clone existing global patches
1984 DynamicList<polyPatch*> newPatches(patches.size());
1985 forAll(patches, patchi)
1986 {
1987 if (!isA<processorPolyPatch>(patches[patchi]))
1988 {
1989 autoPtr<polyPatch> clonedPatch(patches[patchi].clone(patches));
1990 clonedPatch->index() = newPatches.size();
1991 newPatches.append(clonedPatch.ptr());
1992 }
1993 }
1994
1995 // 2. Add new patches
1996 addCouplingPatches
1997 (
1998 mesh,
1999 regionName,
2000 shellRegionName,
2001 zoneNames,
2002 zoneShadowNames,
2003 isInternal,
2004 meshZoneID,
2005
2006 newPatches,
2007 interMeshTopPatch,
2008 interMeshBottomPatch
2009 );
2010
2011 // 3. Clone processor patches
2012 forAll(patches, patchi)
2013 {
2014 if (isA<processorPolyPatch>(patches[patchi]))
2015 {
2016 autoPtr<polyPatch> clonedPatch(patches[patchi].clone(patches));
2017 clonedPatch->index() = newPatches.size();
2018 newPatches.append(clonedPatch.ptr());
2019 }
2020 }
2021
2022
2023 #ifdef FULLDEBUG
2024 Pout<< "*** adaptMesh : addFvPAtches:" << nl;
2025 printPatches(Pout, newPatches);
2026 Pout<< "*** end of adaptMesh : addFvPAtches:" << endl;
2027 #endif
2028
2029
2030 // Add to mesh
2031 mesh.clearOut();
2032 mesh.removeFvBoundary();
2033 mesh.addFvPatches(newPatches, true);
2034
2036 }
2037
2038
2039 // Patch per extruded face
2040 labelList extrudeTopPatchID(extrudePatch.size());
2041 labelList extrudeBottomPatchID(extrudePatch.size());
2042
2043 forAll(zoneID, facei)
2044 {
2045 extrudeTopPatchID[facei] = interRegionTopPatch[zoneID[facei]];
2046 extrudeBottomPatchID[facei] = interRegionBottomPatch[zoneID[facei]];
2047 }
2048
2049
2050
2051 // Count how many patches on special edges of extrudePatch are necessary
2052 // - zoneXXX_sides
2053 // - zoneXXX_zoneYYY
2054 labelList zoneSidePatch(zoneNames.size(), Zero);
2055 // Patch to use for minZone
2056 labelList zoneZonePatch_min(zoneNames.size()*zoneNames.size(), Zero);
2057 // Patch to use for maxZone
2058 labelList zoneZonePatch_max(zoneNames.size()*zoneNames.size(), Zero);
2059
2060 countExtrudePatches
2061 (
2062 mesh,
2063 zoneNames.size(),
2064
2065 extrudePatch, // patch
2066 extrudeMeshFaces, // mesh face per patch face
2067 extrudeMeshEdges, // mesh edge per patch edge
2068
2069 extrudeEdgeGlobalFaces, // global indexing per patch edge
2070 edgeMinZoneID, // minZone per patch edge
2071 edgeMaxZoneID, // maxZone per patch edge
2072
2073 zoneSidePatch, // per zone-side num edges that extrude into it
2074 zoneZonePatch_min // per zone-zone num edges that extrude into it
2075 );
2076
2077 // Now we'll have:
2078 // zoneSidePatch[zoneA] : number of faces needed on the side of zoneA
2079 // zoneZonePatch_min[zoneA,zoneB] : number of faces needed inbetween A,B
2080
2081
2082 // Add the zone-side patches.
2083 addZoneSidePatches
2084 (
2085 mesh,
2086 zoneNames,
2087 (oneD ? oneDPatchType : word::null),
2088
2089 regionPatches,
2090 zoneSidePatch
2091 );
2092
2093
2094 // Add the patches inbetween zones
2095 addInterZonePatches
2096 (
2097 mesh,
2098 zoneNames,
2099 oneD,
2100
2101 zoneZonePatch_min,
2102 zoneZonePatch_max,
2103 regionPatches
2104 );
2105
2106
2107 // Sets sidePatchID[edgeI] to interprocessor patch. Adds any
2108 // interprocessor or cyclic patches if necessary.
2109 labelList sidePatchID;
2110 addCoupledPatches
2111 (
2112 mesh,
2113 extrudePatch,
2114 extrudeMeshFaces,
2115 extrudeMeshEdges,
2116 extrudeEdgeFacesMap,
2117 extrudeEdgeGlobalFaces,
2118
2119 sidePatchID,
2120 regionPatches
2121 );
2122
2123
2124// // Add all the newPatches to the mesh and fields
2125// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2126// {
2127// forAll(newPatches, patchi)
2128// {
2129// Pout<< "Adding patch " << patchi
2130// << " name:" << newPatches[patchi]->name()
2131// << endl;
2132// }
2133// //label nOldPatches = mesh.boundary().size();
2134// mesh.clearOut();
2135// mesh.removeFvBoundary();
2136// mesh.addFvPatches(newPatches, true);
2137// //// Add calculated fvPatchFields for the added patches
2138// //for
2139// //(
2140// // label patchi = nOldPatches;
2141// // patchi < mesh.boundary().size();
2142// // patchi++
2143// //)
2144// //{
2145// // Pout<< "ADDing calculated to patch " << patchi
2146// // << endl;
2147// // addCalculatedPatchFields(mesh);
2148// //}
2149// //Pout<< "** Added " << mesh.boundary().size()-nOldPatches
2150// // << " patches." << endl;
2151// }
2152
2153
2154 // Set patches to use for edges to be extruded into boundary faces
2155 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2156 // In order of edgeFaces: per edge, per originating face the
2157 // patch to use for the side face (from the extruded edge).
2158 // If empty size create an internal face.
2159 labelListList extrudeEdgePatches(extrudePatch.nEdges());
2160
2161 // Is edge a non-manifold edge
2162 bitSet nonManifoldEdge(extrudePatch.nEdges());
2163
2164 // Note: logic has to be same as in countExtrudePatches.
2165 forAll(edgeFaces, edgeI)
2166 {
2167 const labelList& eFaces = edgeFaces[edgeI];
2168
2169 labelList& ePatches = extrudeEdgePatches[edgeI];
2170
2171 if (oneD)
2172 {
2173 ePatches.setSize(eFaces.size());
2174 forAll(eFaces, i)
2175 {
2176 ePatches[i] = zoneSidePatch[zoneID[eFaces[i]]];
2177 }
2178
2179 if (oneDNonManifoldEdges)
2180 {
2181 //- Set nonManifoldEdge[edgeI] for non-manifold edges only
2182 // The other option is to have non-manifold edges everywhere
2183 // and generate space overlapping columns of cells.
2184 if (eFaces.size() != 2)
2185 {
2186 nonManifoldEdge.set(edgeI);
2187 }
2188 }
2189 else
2190 {
2191 nonManifoldEdge.set(edgeI);
2192 }
2193 }
2194 else if (eFaces.size() == 2)
2195 {
2196 label zone0 = zoneID[eFaces[0]];
2197 label zone1 = zoneID[eFaces[1]];
2198
2199 if (zone0 != zone1) // || (cos(angle) > blabla))
2200 {
2201 label minZone = min(zone0,zone1);
2202 label maxZone = max(zone0,zone1);
2203 label index = minZone*zoneNames.size()+maxZone;
2204
2205 ePatches.setSize(eFaces.size());
2206
2207 if (zone0 == minZone)
2208 {
2209 ePatches[0] = zoneZonePatch_min[index];
2210 ePatches[1] = zoneZonePatch_max[index];
2211 }
2212 else
2213 {
2214 ePatches[0] = zoneZonePatch_max[index];
2215 ePatches[1] = zoneZonePatch_min[index];
2216 }
2217
2218 nonManifoldEdge.set(edgeI);
2219 }
2220 }
2221 else if (sidePatchID[edgeI] != -1)
2222 {
2223 // Coupled extrusion
2224 ePatches.setSize(eFaces.size());
2225 forAll(eFaces, i)
2226 {
2227 ePatches[i] = sidePatchID[edgeI];
2228 }
2229 }
2230 else
2231 {
2232 label facei = findUncoveredPatchFace
2233 (
2234 mesh,
2235 labelUIndList(extrudeMeshFaces, eFaces),
2236 extrudeMeshEdges[edgeI]
2237 );
2238
2239 if (facei != -1)
2240 {
2241 label newPatchi = findPatchID
2242 (
2243 regionPatches,
2244 patches[patches.whichPatch(facei)].name()
2245 );
2246 ePatches.setSize(eFaces.size(), newPatchi);
2247 }
2248 else
2249 {
2250 ePatches.setSize(eFaces.size());
2251 forAll(eFaces, i)
2252 {
2253 ePatches[i] = zoneSidePatch[zoneID[eFaces[i]]];
2254 }
2255 }
2256 nonManifoldEdge.set(edgeI);
2257 }
2258 }
2259
2260
2261
2262 // Assign point regions
2263 // ~~~~~~~~~~~~~~~~~~~~
2264
2265 // Per face, per point the region number.
2266 faceList pointGlobalRegions;
2267 faceList pointLocalRegions;
2268 labelList localToGlobalRegion;
2269
2270 createShellMesh::calcPointRegions
2271 (
2272 mesh.globalData(),
2273 extrudePatch,
2274 nonManifoldEdge,
2275 false, // keep cyclic separated regions apart
2276
2277 pointGlobalRegions,
2278 pointLocalRegions,
2279 localToGlobalRegion
2280 );
2281
2282 // Per local region an originating point
2283 labelList localRegionPoints(localToGlobalRegion.size());
2284 forAll(pointLocalRegions, facei)
2285 {
2286 const face& f = extrudePatch.localFaces()[facei];
2287 const face& pRegions = pointLocalRegions[facei];
2288 forAll(pRegions, fp)
2289 {
2290 localRegionPoints[pRegions[fp]] = f[fp];
2291 }
2292 }
2293
2294 // Calculate region normals by reducing local region normals
2295 pointField localRegionNormals(localToGlobalRegion.size());
2296 {
2297 pointField localSum(localToGlobalRegion.size(), Zero);
2298
2299 forAll(pointLocalRegions, facei)
2300 {
2301 const face& pRegions = pointLocalRegions[facei];
2302 forAll(pRegions, fp)
2303 {
2304 label localRegionI = pRegions[fp];
2305 localSum[localRegionI] += extrudePatch.faceNormals()[facei];
2306 }
2307 }
2308
2309 Map<point> globalSum(2*localToGlobalRegion.size());
2310
2311 forAll(localSum, localRegionI)
2312 {
2313 label globalRegionI = localToGlobalRegion[localRegionI];
2314 globalSum.insert(globalRegionI, localSum[localRegionI]);
2315 }
2316
2317 // Reduce
2318 Pstream::mapCombineAllGather(globalSum, plusEqOp<point>());
2319
2320 forAll(localToGlobalRegion, localRegionI)
2321 {
2322 label globalRegionI = localToGlobalRegion[localRegionI];
2323 localRegionNormals[localRegionI] = globalSum[globalRegionI];
2324 }
2325 localRegionNormals /= mag(localRegionNormals);
2326 }
2327
2328
2329 // For debugging: dump hedgehog plot of normals
2330 if (false)
2331 {
2332 OFstream str(runTime.path()/"localRegionNormals.obj");
2333 label vertI = 0;
2334
2335 scalar thickness = model().sumThickness(1);
2336
2337 forAll(pointLocalRegions, facei)
2338 {
2339 const face& f = extrudeFaces[facei];
2340
2341 forAll(f, fp)
2342 {
2343 label region = pointLocalRegions[facei][fp];
2344 const point& pt = extrudePoints[f[fp]];
2345
2346 meshTools::writeOBJ(str, pt);
2347 vertI++;
2348 meshTools::writeOBJ
2349 (
2350 str,
2351 pt+thickness*localRegionNormals[region]
2352 );
2353 vertI++;
2354 str << "l " << vertI-1 << ' ' << vertI << nl;
2355 }
2356 }
2357 }
2358
2359
2360 // Use model to create displacements of first layer
2361 vectorField firstDisp(localRegionNormals.size());
2362 forAll(firstDisp, regionI)
2363 {
2364 //const point& regionPt = regionCentres[regionI];
2365 const point& regionPt = extrudePatch.points()
2366 [
2367 extrudePatch.meshPoints()
2368 [
2369 localRegionPoints[regionI]
2370 ]
2371 ];
2372 const vector& n = localRegionNormals[regionI];
2373 firstDisp[regionI] = model()(regionPt, n, 1) - regionPt;
2374 }
2375
2376
2377 // Create a new mesh
2378 // ~~~~~~~~~~~~~~~~~
2379
2380 createShellMesh extruder
2381 (
2382 extrudePatch,
2383 pointLocalRegions,
2384 localRegionPoints
2385 );
2386
2387
2388 autoPtr<mapPolyMesh> shellMap;
2389 fvMesh regionMesh
2390 (
2391 IOobject
2392 (
2393 shellRegionName,
2394 meshInstance,
2395 runTime,
2396 IOobject::NO_READ,
2397 IOobject::AUTO_WRITE,
2398 false
2399 ),
2400 Zero,
2401 false
2402 );
2403
2404 // Add the new patches
2405 forAll(regionPatches, patchi)
2406 {
2407 polyPatch* ppPtr = regionPatches[patchi];
2408 regionPatches[patchi] = ppPtr->clone(regionMesh.boundaryMesh()).ptr();
2409 delete ppPtr;
2410 }
2411
2412 #ifdef FULLDEBUG
2413 Pout<< "*** regionPatches : regionPatches:" << nl;
2414 printPatches(Pout, regionPatches);
2415 Pout<< "*** end of regionPatches : regionPatches:" << endl;
2416 #endif
2417
2418
2419 regionMesh.clearOut();
2420 regionMesh.removeFvBoundary();
2421 regionMesh.addFvPatches(regionPatches, true);
2422
2423 {
2424 polyTopoChange meshMod(regionPatches.size());
2425
2426 extruder.setRefinement
2427 (
2428 firstDisp, // first displacement
2429 model().expansionRatio(),
2430 model().nLayers(), // nLayers
2431 extrudeTopPatchID,
2432 extrudeBottomPatchID,
2433 extrudeEdgePatches,
2434 meshMod
2435 );
2436
2437 // Enforce actual point positions according to extrudeModel (model)
2438 // (extruder.setRefinement only does fixed expansionRatio)
2439 // The regionPoints and nLayers are looped in the same way as in
2440 // createShellMesh
2441 DynamicList<point>& newPoints = const_cast<DynamicList<point>&>
2442 (
2443 meshMod.points()
2444 );
2445 label meshPointi = extrudePatch.localPoints().size();
2446 forAll(localRegionPoints, regionI)
2447 {
2448 label pointi = localRegionPoints[regionI];
2449 point pt = extrudePatch.localPoints()[pointi];
2450 const vector& n = localRegionNormals[regionI];
2451
2452 for (label layerI = 1; layerI <= model().nLayers(); layerI++)
2453 {
2454 newPoints[meshPointi++] = model()(pt, n, layerI);
2455 }
2456 }
2457
2458 shellMap = meshMod.changeMesh
2459 (
2460 regionMesh, // mesh to change
2461 false // inflate
2462 );
2463 }
2464
2465 // Necessary?
2466 regionMesh.setInstance(meshInstance);
2467
2468
2469 // Update numbering on extruder.
2470 extruder.updateMesh(shellMap());
2471
2472
2473 // Calculate offsets from shell mesh back to original mesh
2474 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2475
2476 List<pointField> topOffsets(zoneNames.size());
2477 List<pointField> bottomOffsets(zoneNames.size());
2478
2479 forAll(regionMesh.boundaryMesh(), patchi)
2480 {
2481 const polyPatch& pp = regionMesh.boundaryMesh()[patchi];
2482
2483 if (isA<mappedWallPolyPatch>(pp))
2484 {
2485 if (interRegionTopPatch.found(patchi))
2486 {
2487 label zoneI = interRegionTopPatch.find(patchi);
2488 topOffsets[zoneI] = calcOffset(extrudePatch, extruder, pp);
2489 }
2490 else if (interRegionBottomPatch.found(patchi))
2491 {
2492 label zoneI = interRegionBottomPatch.find(patchi);
2493 bottomOffsets[zoneI] = calcOffset(extrudePatch, extruder, pp);
2494 }
2495 }
2496 }
2497
2498
2499 // Change top and bottom boundary conditions on regionMesh
2500 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2501
2502 {
2503 // Correct top patches for offset
2504 setCouplingInfo
2505 (
2506 regionMesh,
2507 interRegionTopPatch,
2508 regionName, // name of main mesh
2509 sampleMode, // sampleMode
2510 topOffsets
2511 );
2512
2513 // Correct bottom patches for offset
2514 setCouplingInfo
2515 (
2516 regionMesh,
2517 interRegionBottomPatch,
2518 regionName,
2519 sampleMode, // sampleMode
2520 bottomOffsets
2521 );
2522
2523 // Remove any unused patches
2524 deleteEmptyPatches(regionMesh);
2525 }
2526
2527 // Change top and bottom boundary conditions on main mesh
2528 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2529
2530 if (adaptMesh)
2531 {
2532 // Correct top patches for offset
2533 setCouplingInfo
2534 (
2535 mesh,
2536 interMeshTopPatch,
2537 shellRegionName, // name of shell mesh
2538 sampleMode, // sampleMode
2539 -topOffsets
2540 );
2541
2542 // Correct bottom patches for offset
2543 setCouplingInfo
2544 (
2545 mesh,
2546 interMeshBottomPatch,
2547 shellRegionName,
2548 sampleMode,
2549 -bottomOffsets
2550 );
2551 }
2552
2553
2554
2555 // Write addressing from region mesh back to originating patch
2556 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2557
2558 labelIOList cellToPatchFaceAddressing
2559 (
2560 IOobject
2561 (
2562 "cellToPatchFaceAddressing",
2563 regionMesh.facesInstance(),
2564 regionMesh.meshSubDir,
2565 regionMesh,
2566 IOobject::NO_READ,
2567 IOobject::NO_WRITE,
2568 false
2569 ),
2570 extruder.cellToFaceMap()
2571 );
2572 cellToPatchFaceAddressing.note() = "cell to patch face addressing";
2573
2574 labelIOList faceToPatchFaceAddressing
2575 (
2576 IOobject
2577 (
2578 "faceToPatchFaceAddressing",
2579 regionMesh.facesInstance(),
2580 regionMesh.meshSubDir,
2581 regionMesh,
2582 IOobject::NO_READ,
2583 IOobject::NO_WRITE,
2584 false
2585 ),
2586 extruder.faceToFaceMap()
2587 );
2588 faceToPatchFaceAddressing.note() =
2589 "front/back face + turning index to patch face addressing";
2590
2591 labelIOList faceToPatchEdgeAddressing
2592 (
2593 IOobject
2594 (
2595 "faceToPatchEdgeAddressing",
2596 regionMesh.facesInstance(),
2597 regionMesh.meshSubDir,
2598 regionMesh,
2599 IOobject::NO_READ,
2600 IOobject::NO_WRITE,
2601 false
2602 ),
2603 extruder.faceToEdgeMap()
2604 );
2605 faceToPatchEdgeAddressing.note() =
2606 "side face to patch edge addressing";
2607
2608 labelIOList pointToPatchPointAddressing
2609 (
2610 IOobject
2611 (
2612 "pointToPatchPointAddressing",
2613 regionMesh.facesInstance(),
2614 regionMesh.meshSubDir,
2615 regionMesh,
2616 IOobject::NO_READ,
2617 IOobject::NO_WRITE,
2618 false
2619 ),
2620 extruder.pointToPointMap()
2621 );
2622 pointToPatchPointAddressing.note() =
2623 "point to patch point addressing";
2624
2625
2626 Info<< "Writing mesh " << regionMesh.name()
2627 << " to " << regionMesh.facesInstance() << nl
2628 << endl;
2629
2630 bool ok =
2631 regionMesh.write()
2632 && cellToPatchFaceAddressing.write()
2633 && faceToPatchFaceAddressing.write()
2634 && faceToPatchEdgeAddressing.write()
2635 && pointToPatchPointAddressing.write();
2636
2637 if (!ok)
2638 {
2640 << "Failed writing mesh " << regionMesh.name()
2641 << " at location " << regionMesh.facesInstance()
2642 << exit(FatalError);
2643 }
2644 topoSet::removeFiles(regionMesh);
2645 processorMeshes::removeFiles(regionMesh);
2646
2647
2648 // See if we need to extrude coordinates as well
2649 {
2650 autoPtr<pointIOField> patchFaceCentresPtr;
2651
2652 IOobject io
2653 (
2654 "patchFaceCentres",
2655 mesh.pointsInstance(),
2656 mesh.meshSubDir,
2657 mesh,
2658 IOobject::MUST_READ
2659 );
2660 if (io.typeHeaderOk<pointIOField>(true))
2661 {
2662 // Read patchFaceCentres and patchEdgeCentres
2663 Info<< "Reading patch face,edge centres : "
2664 << io.name() << " and patchEdgeCentres" << endl;
2665
2666 extrudeGeometricProperties
2667 (
2668 mesh,
2669 extrudePatch,
2670 extruder,
2671 regionMesh,
2672 model()
2673 );
2674 }
2675 }
2676
2677
2678
2679
2680 // Insert baffles into original mesh
2681 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2682
2683 autoPtr<mapPolyMesh> addBafflesMap;
2684
2685 if (adaptMesh)
2686 {
2687 polyTopoChange meshMod(mesh);
2688
2689 // Modify faces to be in bottom (= always coupled) patch
2690 forAll(extrudeMeshFaces, zoneFacei)
2691 {
2692 label meshFacei = extrudeMeshFaces[zoneFacei];
2693 label zoneI = zoneID[zoneFacei];
2694 bool flip = zoneFlipMap[zoneFacei];
2695 const face& f = mesh.faces()[meshFacei];
2696
2697 if (!flip)
2698 {
2699 meshMod.modifyFace
2700 (
2701 f, // modified face
2702 meshFacei, // label of face being modified
2703 mesh.faceOwner()[meshFacei],// owner
2704 -1, // neighbour
2705 false, // face flip
2706 interMeshBottomPatch[zoneI],// patch for face
2707 meshZoneID[zoneI], // zone for face
2708 flip // face flip in zone
2709 );
2710 }
2711 else if (mesh.isInternalFace(meshFacei))
2712 {
2713 meshMod.modifyFace
2714 (
2715 f.reverseFace(), // modified face
2716 meshFacei, // label of modified face
2717 mesh.faceNeighbour()[meshFacei],// owner
2718 -1, // neighbour
2719 true, // face flip
2720 interMeshBottomPatch[zoneI], // patch for face
2721 meshZoneID[zoneI], // zone for face
2722 !flip // face flip in zone
2723 );
2724 }
2725 }
2726
2727 if (zoneShadowNames.size() > 0) //if there is a top faceZone specified
2728 {
2729 forAll(extrudeMeshFaces, zoneFacei)
2730 {
2731 label meshFacei = extrudeMeshShadowFaces[zoneFacei];
2732 label zoneI = zoneShadowID[zoneFacei];
2733 bool flip = zoneShadowFlipMap[zoneFacei];
2734 const face& f = mesh.faces()[meshFacei];
2735
2736 if (!flip)
2737 {
2738 meshMod.modifyFace
2739 (
2740 f, // modified face
2741 meshFacei, // face being modified
2742 mesh.faceOwner()[meshFacei],// owner
2743 -1, // neighbour
2744 false, // face flip
2745 interMeshTopPatch[zoneI], // patch for face
2746 meshZoneID[zoneI], // zone for face
2747 flip // face flip in zone
2748 );
2749 }
2750 else if (mesh.isInternalFace(meshFacei))
2751 {
2752 meshMod.modifyFace
2753 (
2754 f.reverseFace(), // modified face
2755 meshFacei, // label modified face
2756 mesh.faceNeighbour()[meshFacei],// owner
2757 -1, // neighbour
2758 true, // face flip
2759 interMeshTopPatch[zoneI], // patch for face
2760 meshZoneID[zoneI], // zone for face
2761 !flip // face flip in zone
2762 );
2763 }
2764 }
2765 }
2766 else
2767 {
2768 // Add faces (using same points) to be in top patch
2769 forAll(extrudeMeshFaces, zoneFacei)
2770 {
2771 label meshFacei = extrudeMeshFaces[zoneFacei];
2772 label zoneI = zoneID[zoneFacei];
2773 bool flip = zoneFlipMap[zoneFacei];
2774 const face& f = mesh.faces()[meshFacei];
2775
2776 if (!flip)
2777 {
2778 if (mesh.isInternalFace(meshFacei))
2779 {
2780 meshMod.addFace
2781 (
2782 f.reverseFace(), // modified face
2783 mesh.faceNeighbour()[meshFacei],// owner
2784 -1, // neighbour
2785 -1, // master point
2786 -1, // master edge
2787 meshFacei, // master face
2788 true, // flip flux
2789 interMeshTopPatch[zoneI], // patch for face
2790 -1, // zone for face
2791 false //face flip in zone
2792 );
2793 }
2794 }
2795 else
2796 {
2797 meshMod.addFace
2798 (
2799 f, // face
2800 mesh.faceOwner()[meshFacei], // owner
2801 -1, // neighbour
2802 -1, // master point
2803 -1, // master edge
2804 meshFacei, // master face
2805 false, // flip flux
2806 interMeshTopPatch[zoneI], // patch for face
2807 -1, // zone for face
2808 false // zone flip
2809 );
2810 }
2811 }
2812 }
2813
2814 // Change the mesh. Change points directly (no inflation).
2815 addBafflesMap = meshMod.changeMesh(mesh, false);
2816
2817 // Update fields
2818 mesh.updateMesh(addBafflesMap());
2819
2820
2821//XXXXXX
2822// Update maps! e.g. faceToPatchFaceAddressing
2823//XXXXXX
2824
2825 // Move mesh (since morphing might not do this)
2826 if (addBafflesMap().hasMotionPoints())
2827 {
2828 mesh.movePoints(addBafflesMap().preMotionPoints());
2829 }
2830
2831 mesh.setInstance(meshInstance);
2832
2833 // Remove any unused patches
2834 deleteEmptyPatches(mesh);
2835
2836 Info<< "Writing mesh " << mesh.name()
2837 << " to " << mesh.facesInstance() << nl
2838 << endl;
2839
2840 if (!mesh.write())
2841 {
2843 << "Failed writing mesh " << mesh.name()
2844 << " at location " << mesh.facesInstance()
2845 << exit(FatalError);
2846 }
2847 topoSet::removeFiles(mesh);
2848 processorMeshes::removeFiles(mesh);
2849 }
2850
2851 Info << "End\n" << endl;
2852
2853 return 0;
2854}
2855
2856
2857// ************************************************************************* //
label n
Y[inertIndex] max(0.0)
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
A primitive field of type <T> with automated input and output.
Definition: IOField.H:58
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
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:65
label size() const noexcept
The number of elements in the list.
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
Definition: List.H:330
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
OFstream that keeps track of vertices.
Definition: OBJstream.H:61
Output to file stream, using an OSstream.
Definition: OFstream.H:57
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:62
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
A list of faces which address into the list of points.
label nEdges() const
Number of edges in patch.
label nPoints() const
Number of points supporting patch faces.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
labelList meshEdges(const edgeList &allEdges, const labelListList &cellEdges, const labelList &faceCells) const
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const Field< point_type > & localPoints() const
Return pointField of points in patch.
const Field< point_type > & faceNormals() const
Return face unit normals for patch.
const Field< point_type > & points() const noexcept
Return reference to global points.
const labelListList & edgeFaces() const
Return edge-face addressing.
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
SubField is a Field obtained as a section of another Field, without its own allocation....
Definition: SubField.H:62
fileName path() const
Return path.
Definition: Time.H:358
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: UList.H:94
bool found(const T &val, label pos=0) const
True if the value if found in the list.
Definition: UListI.H:265
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
Definition: UList.C:212
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 size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
label findZoneID(const word &zoneName) const
Find zone index by name, return -1 if not found.
Definition: ZoneMesh.C:525
wordList names() const
A list of the zone names.
Definition: ZoneMesh.C:304
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
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
Creates mesh by extruding a patch.
const labelList & faceToEdgeMap() const
From region side-face to patch edge. -1 for non-edge faces.
const labelList & pointToPointMap() const
From region point to patch point.
const labelList & faceToFaceMap() const
From region face to patch face. Contains turning index:
const labelList & cellToFaceMap() const
From region cell to patch face. Consecutively added so.
void updateMesh(const mapPolyMesh &)
Update any locally stored mesh information.
void setRefinement(const pointField &firstLayerThickness, const scalar expansionRatio, const label nLayers, const labelList &topPatchID, const labelList &bottomPatchID, const labelListList &extrudeEdgePatches, polyTopoChange &meshMod)
Play commands into polyTopoChange to create layer mesh.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
entry * add(entry *entryPtr, bool mergeEntry=false)
Add a new entry.
Definition: dictionary.C:640
entry * set(entry *entryPtr)
Assign a new entry, overwriting any existing entry.
Definition: dictionary.C:780
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
Top level extrusion model class.
Definition: extrudeModel.H:77
label nLayers() const
Return the number of layers.
Definition: extrudeModel.C:59
scalar sumThickness(const label layer) const
Helper: calculate cumulative relative thickness for layer.
Definition: extrudeModel.C:71
A list of face labels.
Definition: faceSet.H:54
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:67
const boolList & flipMap() const noexcept
Return face flip map.
Definition: faceZone.H:272
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:68
label toGlobal(const label i) const
From local to global index.
Definition: globalIndexI.H:260
A line primitive.
Definition: line.H:68
Class containing processor-to-processor mapping information.
void distribute(List< T > &fld, const bool dummyTransform=true, const int tag=UPstream::msgType()) const
Distribute data using default commsType.
sampleMode
Mesh items to sample.
const Time & time() const noexcept
Return time registry.
label index() const noexcept
The index of this patch in the boundaryMesh.
const word & name() const noexcept
The patch name.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
void setInstance(const fileName &instance, const IOobject::writeOption wOpt=IOobject::AUTO_WRITE)
Set the instance for mesh files.
Definition: polyMeshIO.C:36
const fileName & facesInstance() const
Return the current instance directory for faces.
Definition: polyMesh.C:866
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1121
const fileName & pointsInstance() const
Return the current instance directory for points.
Definition: polyMesh.C:860
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1127
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
Definition: polyMesh.H:324
void clearOut()
Clear all geometry and addressing unnecessary for CFD.
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
const vectorField::subField faceCentres() const
Return face centres.
Definition: polyPatch.C:321
Direct mesh changes based on v1.3 polyTopoChange syntax.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
const vectorField & faceCentres() const
const vectorField & cellCentres() const
label nCells() const noexcept
Number of mesh cells.
label nFaces() const noexcept
Number of mesh faces.
virtual bool write(const bool valid=true) const
Write using setting from DB.
A class for managing temporary objects.
Definition: tmp.H:65
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
const polyBoundaryMesh & patches
bool coupled(solutionDict.getOrDefault("coupledEnergyField", false))
dynamicFvMesh & mesh
engineTime & runTime
Foam::word regionName(Foam::polyMesh::defaultRegion)
Required Variables.
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:473
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
OBJstream os(runTime.globalPath()/outputName)
const word dictName("faMeshDefinition")
const labelIOList & zoneIDs
Definition: correctPhi.H:59
label nZones
const labelIOList & zoneID
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
#define WarningInFunction
Report a warning using Foam::Warning.
Namespace for OpenFOAM.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Type gAverage(const FieldField< Field, Type > &f)
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)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
labelList f(nPoints)
dictionary dict
IOobject dictIO
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
List helper to append y unique elements onto the end of x.
Definition: ListOps.H:570
Foam::surfaceFields.