addPatchCellLayer.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
27\*---------------------------------------------------------------------------*/
28
29#include "addPatchCellLayer.H"
30#include "polyMesh.H"
31#include "polyTopoChange.H"
32#include "meshTools.H"
33#include "mapPolyMesh.H"
34#include "syncTools.H"
35#include "polyAddPoint.H"
36#include "polyAddFace.H"
37#include "polyModifyFace.H"
38#include "polyAddCell.H"
39#include "globalIndex.H"
40#include "PatchTools.H"
41#include "dummyTransform.H"
42
43// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
44
45namespace Foam
46{
48
49 // Reduction class to get minimum value over face.
51 {
52 public:
53
54 void operator()(face& x, const face& y) const
55 {
56 if (x.size())
57 {
58 if (y.size())
59 {
60 label j = 0;
61 forAll(x, i)
62 {
63 x[i] = min(x[i], y[j]);
64
65 j = y.rcIndex(j);
66 }
67 }
68 }
69 else if (y.size())
70 {
71 x.setSize(y.size());
72 label j = 0;
73 forAll(x, i)
74 {
75 x[i] = y[j];
76 j = y.rcIndex(j);
77 }
78 }
79 }
80 };
81
82}
83
84
85// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
86
87Foam::label Foam::addPatchCellLayer::nbrFace
88(
89 const labelListList& edgeFaces,
90 const label edgei,
91 const label facei
92)
93{
94 const labelList& eFaces = edgeFaces[edgei];
95
96 if (eFaces.size() == 2)
97 {
98 return (eFaces[0] != facei ? eFaces[0] : eFaces[1]);
99 }
100 else
101 {
102 return -1;
103 }
104}
105
106
107void Foam::addPatchCellLayer::addVertex
108(
109 const label pointi,
110 face& f,
111 label& fp
112)
113{
114 if (fp == 0)
115 {
116 f[fp++] = pointi;
117 }
118 else
119 {
120 if (f[fp-1] != pointi && f[0] != pointi)
121 {
122 f[fp++] = pointi;
123 }
124 }
125}
126
127
128// Is edge to the same neighbour? (and needs extrusion and has not been
129// dealt with already)
130bool Foam::addPatchCellLayer::sameEdgeNeighbour
131(
132 const indirectPrimitivePatch& pp,
133 const labelListList& globalEdgeFaces,
134 const boolList& doneEdge,
135 const label thisGlobalFacei,
136 const label nbrGlobalFacei,
137 const label edgei
138) const
139{
140 const edge& e = pp.edges()[edgei];
141
142 return
143 !doneEdge[edgei] // not yet handled
144 && (
145 addedPoints_[e[0]].size() // is extruded
146 || addedPoints_[e[1]].size()
147 )
148 && (
149 nbrFace(globalEdgeFaces, edgei, thisGlobalFacei)
150 == nbrGlobalFacei // is to same neighbour
151 );
152}
153
154
155// Collect consecutive string of edges that connects the same two
156// (possibly coupled) faces. Returns -1 if no unvisited edge can be found.
157// Otherwise returns start and end index in face.
158Foam::labelPair Foam::addPatchCellLayer::getEdgeString
159(
160 const indirectPrimitivePatch& pp,
161 const labelListList& globalEdgeFaces,
162 const boolList& doneEdge,
163 const label patchFacei,
164 const label globalFacei
165) const
166{
167 const labelList& fEdges = pp.faceEdges()[patchFacei];
168
169 label startFp = -1;
170 label endFp = -1;
171
172 // Get edge that hasn't been done yet but needs extrusion
173 forAll(fEdges, fp)
174 {
175 label edgei = fEdges[fp];
176 const edge& e = pp.edges()[edgei];
177
178 if
179 (
180 !doneEdge[edgei]
181 && ( addedPoints_[e[0]].size() || addedPoints_[e[1]].size() )
182 )
183 {
184 startFp = fp;
185 break;
186 }
187 }
188
189 if (startFp != -1)
190 {
191 // We found an edge that needs extruding but hasn't been done yet.
192 // Now find the face on the other side
193 label nbrGlobalFacei = nbrFace
194 (
195 globalEdgeFaces,
196 fEdges[startFp],
197 globalFacei
198 );
199
200 if (nbrGlobalFacei == -1)
201 {
202 // Proper boundary edge. Only extrude single edge.
203 endFp = startFp;
204 }
205 else
206 {
207 // Search back for edge
208 // - which hasn't been handled yet
209 // - with same neighbour
210 // - that needs extrusion
211
212 const label initFp = startFp;
213 while (true)
214 {
215 label prevFp = fEdges.rcIndex(startFp);
216
217 if (prevFp == initFp)
218 {
219 const edge& e = pp.edges()[fEdges[initFp]];
220 const face& localF = pp.localFaces()[patchFacei];
221
223 << "On face:" << patchFacei
224 << " fc:" << pp.faceCentres()[patchFacei]
225 << " vertices:" << localF
226 << " points:"
227 << UIndirectList<point>(pp.points(), pp[patchFacei])
228 << " edges:" << fEdges
229 << " All edges of face seem to have same neighbour "
230 << nbrGlobalFacei
231 << " starting walking from edge " << e
232 << exit(FatalError);
233 }
234
235 if
236 (
237 !sameEdgeNeighbour
238 (
239 pp,
240 globalEdgeFaces,
241 doneEdge,
242 globalFacei,
243 nbrGlobalFacei,
244 fEdges[prevFp]
245 )
246 )
247 {
248 break;
249 }
250 startFp = prevFp;
251 }
252
253 // Search forward for end of string
254 endFp = startFp;
255 while (true)
256 {
257 label nextFp = fEdges.fcIndex(endFp);
258
259 if
260 (
261 !sameEdgeNeighbour
262 (
263 pp,
264 globalEdgeFaces,
265 doneEdge,
266 globalFacei,
267 nbrGlobalFacei,
268 fEdges[nextFp]
269 )
270 )
271 {
272 break;
273 }
274 endFp = nextFp;
275 }
276 }
277 }
278
279 return labelPair(startFp, endFp);
280}
281
282
283Foam::label Foam::addPatchCellLayer::addSideFace
284(
285 const indirectPrimitivePatch& pp,
286 const labelListList& addedCells, // per pp face the new extruded cell
287 const face& newFace,
288 const label newPatchID,
289 const label zonei,
290 const bool edgeFlip,
291 const label inflateFacei,
292
293 const label ownFacei, // pp face that provides owner
294 const label nbrFacei,
295 const label meshEdgei, // corresponding mesh edge
296 const label layeri, // layer
297 const label numEdgeFaces, // number of layers for edge
298 const labelList& meshFaces, // precalculated edgeFaces
299 polyTopoChange& meshMod
300) const
301{
302 // Adds a side face i.e. extrudes a patch edge.
303
304 label addedFacei = -1;
305
306
307 // Is patch edge external edge of indirectPrimitivePatch?
308 if (nbrFacei == -1)
309 {
310 // External edge so external face.
311
312 // Determine if different number of layer on owner and neighbour side
313 // (relevant only for coupled faces). See section for internal edge
314 // below.
315
316 label layerOwn;
317
318 if (addedCells[ownFacei].size() < numEdgeFaces)
319 {
320 label offset = numEdgeFaces - addedCells[ownFacei].size();
321 if (layeri <= offset)
322 {
323 layerOwn = 0;
324 }
325 else
326 {
327 layerOwn = layeri - offset;
328 }
329 }
330 else
331 {
332 layerOwn = layeri;
333 }
334
335
336 //Pout<< "Added boundary face:" << newFace
337 // << " atfc:" << newFace.centre(meshMod.points())
338 // << " n:" << newFace.unitNormal(meshMod.points())
339 // << " own:" << addedCells[ownFacei][layerOwn]
340 // << " patch:" << newPatchID
341 // << endl;
342
343 addedFacei = meshMod.setAction
344 (
345 polyAddFace
346 (
347 newFace, // face
348 addedCells[ownFacei][layerOwn], // owner
349 -1, // neighbour
350 -1, // master point
351 -1, // master edge
352 inflateFacei, // master face
353 false, // flux flip
354 newPatchID, // patch for face
355 zonei, // zone for face
356 edgeFlip // face zone flip
357 )
358 );
359 }
360 else
361 {
362 // When adding side faces we need to modify neighbour and owners
363 // in region where layer mesh is stopped. Determine which side
364 // has max number of faces and make sure layers match closest to
365 // original pp if there are different number of layers.
366
367 label layerNbr;
368 label layerOwn;
369
370 if (addedCells[ownFacei].size() > addedCells[nbrFacei].size())
371 {
372 label offset =
373 addedCells[ownFacei].size() - addedCells[nbrFacei].size();
374
375 layerOwn = layeri;
376
377 if (layeri <= offset)
378 {
379 layerNbr = 0;
380 }
381 else
382 {
383 layerNbr = layeri - offset;
384 }
385 }
386 else if (addedCells[nbrFacei].size() > addedCells[ownFacei].size())
387 {
388 label offset =
389 addedCells[nbrFacei].size() - addedCells[ownFacei].size();
390
391 layerNbr = layeri;
392
393 if (layeri <= offset)
394 {
395 layerOwn = 0;
396 }
397 else
398 {
399 layerOwn = layeri - offset;
400 }
401 }
402 else
403 {
404 // Same number of layers on both sides.
405 layerNbr = layeri;
406 layerOwn = layeri;
407 }
408
409
410 // Check mesh internal faces using edge to initialise
411 label inflateEdgei = -1;
412 if (addToMesh_)
413 {
414 forAll(meshFaces, i)
415 {
416 if (mesh_.isInternalFace(meshFaces[i]))
417 {
418 // meshEdge uses internal faces so ok to inflate from it
419 inflateEdgei = meshEdgei;
420 break;
421 }
422 }
423 }
424
425
426 addedFacei = meshMod.setAction
427 (
428 polyAddFace
429 (
430 newFace, // face
431 addedCells[ownFacei][layerOwn], // owner
432 addedCells[nbrFacei][layerNbr], // neighbour
433 -1, // master point
434 inflateEdgei, // master edge
435 -1, // master face
436 false, // flux flip
437 -1, // patch for face
438 zonei, // zone for face
439 edgeFlip // face zone flip
440 )
441 );
442
443 //Pout<< "Added internal face:" << newFace
444 // << " atfc:" << newFace.centre(meshMod.points())
445 // << " n:" << newFace.unitNormal(meshMod.points())
446 // << " own:" << addedCells[ownFacei][layerOwn]
447 // << " nei:" << addedCells[nbrFacei][layerNbr]
448 // << endl;
449 }
450
451 return addedFacei;
452}
453
454
455Foam::label Foam::addPatchCellLayer::findProcPatch
456(
457 const polyMesh& mesh,
458 const label nbrProcID
459)
460{
461 const polyBoundaryMesh& patches = mesh.boundaryMesh();
462
464 {
465 label patchi = mesh.globalData().processorPatches()[i];
466
467 if
468 (
469 refCast<const processorPolyPatch>(patches[patchi]).neighbProcNo()
470 == nbrProcID
471 )
472 {
473 return patchi;
474 }
475 }
476 return -1;
477}
478
479
480void Foam::addPatchCellLayer::setFaceProps
481(
482 const polyMesh& mesh,
483 const label facei,
484
485 label& patchi,
486 label& zonei,
487 bool& zoneFlip
488)
489{
490 patchi = mesh.boundaryMesh().whichPatch(facei);
491 zonei = mesh.faceZones().whichZone(facei);
492 if (zonei != -1)
493 {
494 label index = mesh.faceZones()[zonei].whichFace(facei);
495 zoneFlip = mesh.faceZones()[zonei].flipMap()[index];
496 }
497}
498
499
500void Foam::addPatchCellLayer::setFaceProps
501(
502 const polyMesh& mesh,
503 const indirectPrimitivePatch& pp,
504 const label ppEdgeI,
505 const label faceI,
506
507 label& patchI,
508 label& zoneI,
509 bool& zoneFlip,
510 label& inflateFaceI
511)
512{
513 setFaceProps
514 (
515 mesh,
516 faceI,
517
518 patchI,
519 zoneI,
520 zoneFlip
521 );
522
523 if (patchI != -1 || zoneI != -1)
524 {
525 inflateFaceI = faceI;
526 }
527
528 if (zoneI != -1)
529 {
530 // Correct flip for patch edge ordering
531 const edge& ppEdge = pp.edges()[ppEdgeI];
532 const edge patchEdge
533 (
534 pp.meshPoints()[ppEdge[0]],
535 pp.meshPoints()[ppEdge[1]]
536 );
537
538 const face& f = mesh.faces()[faceI];
539 bool found = false;
540 forAll(f, fp)
541 {
542 const edge e(f[fp], f.nextLabel(fp));
543 int stat = edge::compare(e, patchEdge);
544 if (stat == 1)
545 {
546 found = true;
547 break;
548 }
549 else if (stat == -1)
550 {
551 found = true;
552 zoneFlip = !zoneFlip;
553 break;
554 }
555 }
556
557 if (!found)
558 {
559 //FatalErrorInFunction
561 << "Problem: cannot find patch edge " << ppEdgeI
562 << " with mesh vertices " << patchEdge
563 << " at " << patchEdge.line(mesh.points())
564 << " in face " << faceI << " with mesh vertices "
565 << f
566 << " at " << pointField(mesh.points(), f)
567 << endl
568 << "Continuing with potentially incorrect faceZone orientation"
569 //<< exit(FatalError);
570 << endl;
571 }
572 }
573}
574
575
576void Foam::addPatchCellLayer::findZoneFace
577(
578 const bool useInternalFaces,
579 const bool useBoundaryFaces,
580
581 const polyMesh& mesh,
582 const indirectPrimitivePatch& pp,
583 const label ppEdgeI,
584 const labelUIndList& excludeFaces,
585 const labelList& meshFaces,
586
587 label& inflateFaceI,
588 label& patchI,
589 label& zoneI,
590 bool& zoneFlip
591)
592{
593 inflateFaceI = -1;
594 patchI = -1;
595 zoneI = -1;
596 zoneFlip = false;
597
598 forAll(meshFaces, k)
599 {
600 label faceI = meshFaces[k];
601
602 if
603 (
604 !excludeFaces.found(faceI)
605 && (
606 (mesh.isInternalFace(faceI) && useInternalFaces)
607 || (!mesh.isInternalFace(faceI) && useBoundaryFaces)
608 )
609 )
610 {
611 setFaceProps
612 (
613 mesh,
614 pp,
615 ppEdgeI,
616 faceI,
617
618 patchI,
619 zoneI,
620 zoneFlip,
621 inflateFaceI
622 );
623
624 if (zoneI != -1 || patchI != -1)
625 {
626 break;
627 }
628 }
629 }
630}
631
632
633// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
634
635// Construct from mesh
637(
638 const polyMesh& mesh,
639 const bool addToMesh
640)
641:
642 mesh_(mesh),
643 addToMesh_(addToMesh),
644 addedPoints_(0),
645 layerFaces_(0)
646{}
647
648
649// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
650
652(
653 const polyMesh& mesh,
654 const labelListList& layerFaces
655)
656{
657 labelListList layerCells(layerFaces.size());
658
659 forAll(layerFaces, patchFacei)
660 {
661 const labelList& faceLabels = layerFaces[patchFacei];
662
663 if (faceLabels.size())
664 {
665 labelList& added = layerCells[patchFacei];
666 added.setSize(faceLabels.size()-1);
667
668 for (label i = 0; i < faceLabels.size()-1; i++)
669 {
670 added[i] = mesh.faceNeighbour()[faceLabels[i]];
671 }
672 }
673 }
674 return layerCells;
675}
676
677
679{
680 return addedCells(mesh_, layerFaces_);
681}
682
683
685(
686 const polyMesh& mesh,
687 const globalIndex& globalFaces,
688 const indirectPrimitivePatch& pp
689)
690{
691 // Precalculate mesh edges for pp.edges.
692 const labelList meshEdges(pp.meshEdges(mesh.edges(), mesh.pointEdges()));
693
694 // From mesh edge to global face labels. Non-empty sublists only for
695 // pp edges.
696 labelListList globalEdgeFaces(mesh.nEdges());
697
698 const labelListList& edgeFaces = pp.edgeFaces();
699
700 forAll(edgeFaces, edgeI)
701 {
702 label meshEdgeI = meshEdges[edgeI];
703
704 const labelList& eFaces = edgeFaces[edgeI];
705
706 // Store face and processor as unique tag.
707 labelList& globalEFaces = globalEdgeFaces[meshEdgeI];
708 globalEFaces.setSize(eFaces.size());
709 forAll(eFaces, i)
710 {
711 globalEFaces[i] = globalFaces.toGlobal(pp.addressing()[eFaces[i]]);
712 }
713 }
714
715 // Synchronise across coupled edges.
717 (
718 mesh,
719 globalEdgeFaces,
721 labelList() // null value
722 );
723
724 // Extract pp part
725 return labelListList(UIndirectList<labelList>(globalEdgeFaces, meshEdges));
726}
727
728
729void Foam::addPatchCellLayer::markPatchEdges
730(
731 const polyMesh& mesh,
732 const indirectPrimitivePatch& pp,
733 const labelListList& edgeGlobalFaces,
734 const labelList& meshEdges,
735
736 bitSet& isPatchEdge,
737 bitSet& isPatchBoundaryEdge
738)
739{
740 // Mark (mesh) edges:
741 // - anywhere on extrusion
742 // - where the extrusion ends
743
744 isPatchEdge.setSize(mesh.nEdges());
745 isPatchEdge = false;
746 isPatchEdge.set(meshEdges);
747 // Synchronise across coupled edges
749 (
750 mesh,
751 isPatchEdge,
753 false // initial value
754 );
755
756 isPatchBoundaryEdge.setSize(mesh.nEdges());
757 isPatchBoundaryEdge = false;
758 forAll(edgeGlobalFaces, edgei)
759 {
760 // Test that edge has single global extruded face.
761 // Mark on processor that holds the face (since edgeGlobalFaces
762 // only gets filled from pp faces so if there is only one this
763 // is it)
764 if (edgeGlobalFaces[edgei].size() == 1)
765 {
766 isPatchBoundaryEdge.set(meshEdges[edgei]);
767 }
768 }
769 // Synchronise across coupled edges
771 (
772 mesh,
773 isPatchBoundaryEdge,
774 orEqOp<unsigned int>(),
775 false // initial value
776 );
777}
778
779
780void Foam::addPatchCellLayer::globalEdgeInfo
781(
782 const bool zoneFromAnyFace,
783
784 const polyMesh& mesh,
785 const globalIndex& globalFaces,
786 const labelListList& edgeGlobalFaces,
787 const indirectPrimitivePatch& pp,
788 const labelList& meshEdges,
789
790 labelList& patchEdgeToFace, // face (in globalFaces index)
791 labelList& patchEdgeToPatch, // patch on face (or -1 for internal faces)
792 labelList& patchEdgeToZone, // zone on face
793 bitSet& patchEdgeToFlip // flip orientation on face
794)
795{
796 // For every edge on the outside of the patch return a potential patch/
797 // faceZone to extrude into.
798
799 // Mark (mesh) edges on pp.
800 bitSet isExtrudeEdge;
801 bitSet isBoundaryEdge;
802 markPatchEdges
803 (
804 mesh,
805 pp,
806 edgeGlobalFaces,
807 meshEdges,
808
809 isExtrudeEdge,
810 isBoundaryEdge
811 );
812
813 // Build map of pp edges (in mesh point indexing). Note that this
814 // is now also on processors that do not have pp (but do have the edge)
815 EdgeMap<label> isBoundaryEdgeSet(pp.nEdges());
816 for (const label edgei : isBoundaryEdge)
817 {
818 isBoundaryEdgeSet.insert(mesh.edges()[edgei], edgei);
819 }
820 EdgeMap<label> isExtrudeEdgeSet(pp.nEdges());
821 for (const label edgei : isExtrudeEdge)
822 {
823 isExtrudeEdgeSet.insert(mesh.edges()[edgei], edgei);
824 }
825
826
827 const faceZoneMesh& fzs = mesh.faceZones();
828
829 // Extract zone info into mesh face indexing for ease of addressing
830 labelList faceToZone(mesh.nFaces(), -1);
831 bitSet faceToFlip(mesh.nFaces());
832 for (const faceZone& fz: fzs)
833 {
834 const labelList& addressing = fz;
835 UIndirectList<label>(faceToZone, addressing) = fz.index();
836
837 const boolList& fm = fz.flipMap();
838 forAll(addressing, i)
839 {
840 faceToFlip[addressing[i]] = fm[i];
841 }
842 }
843
844
845 // Storage (over all mesh edges)
846 // - face that data originates from (in globalFaces indexing)
847 labelList meshEdgeToFace(mesh.nEdges(), -1);
848 // - patch (for boundary faces)
849 labelList meshEdgeToPatch(mesh.nEdges(), -1);
850 // - faceZone
851 labelList meshEdgeToZone(mesh.nEdges(), -1);
852 // - faceZone orientation
853 bitSet meshEdgeToFlip(mesh.nEdges());
854
855 //if (useInternalFaces)
856 {
857 const bitSet isInternalOrCoupled
858 (
860 );
861
862 // Loop over edges of the face to find any faceZone.
863 // Edges kept as point pair so we don't construct mesh.faceEdges etc.
864
865 for (const label facei : isInternalOrCoupled)
866 {
867 const face& f = mesh.faces()[facei];
868
869 label prevPointi = f.last();
870 for (const label pointi : f)
871 {
872 const edge e(prevPointi, pointi);
873
874 // Check if edge is internal to extrusion. Take over faceZone
875 // etc from internal face.
876 const auto eFnd = isExtrudeEdgeSet.cfind(e);
877 if (eFnd)
878 {
879 const label edgei = eFnd();
880
881 if (faceToZone[facei] != -1)
882 {
883 // Found a zoned internal face. Use.
884 meshEdgeToFace[edgei] = globalFaces.toGlobal(facei);
885 meshEdgeToZone[edgei] = faceToZone[facei];
886 const edge& meshE = mesh.edges()[edgei];
887 const int d = edge::compare(e, meshE);
888 if (d == 1)
889 {
890 meshEdgeToFlip[edgei] = faceToFlip[facei];
891 }
892 else if (d == -1)
893 {
894 meshEdgeToFlip[edgei] = !faceToFlip[facei];
895 }
896 else
897 {
898 FatalErrorInFunction << "big problem"
899 << exit(FatalError);
900 }
901 }
902 }
903
904 prevPointi = pointi;
905 }
906 }
907 }
908
909
910 //if (useBoundaryFaces)
911 {
912 // Loop over all patches and find 'best' one (non-coupled,
913 // non-extrusion, non-constraint(?)). Note that logic is slightly
914 // different from internal faces loop above - first patch face
915 // is being used instead of first zoned face for internal faces
916
917 const polyBoundaryMesh& patches = mesh.boundaryMesh();
918
919 bitSet isPpFace(mesh.nFaces());
920 isPpFace.set(pp.addressing());
921 // Note: no need to sync ppFace since does not include processor patches
922
923 for (const polyPatch& pp : patches)
924 {
925 if (!pp.coupled())
926 {
927 // TBD. Check for constraint? This is usually a geometric
928 // constraint and not a topological one so should
929 // be handled in the extrusion vector calculation instead?
930
931 forAll(pp, i)
932 {
933 const label facei = pp.start()+i;
934
935 if (!isPpFace[facei])
936 {
937 const face& f = pp[i];
938
939 label prevPointi = f.last();
940 for (const label pointi : f)
941 {
942 const edge e(prevPointi, pointi);
943 const auto eFnd =
944 (
945 zoneFromAnyFace
946 ? isExtrudeEdgeSet.cfind(e)
947 : isBoundaryEdgeSet.cfind(e)
948 );
949 if (eFnd)
950 {
951 const label edgei = eFnd();
952 if (meshEdgeToFace[edgei] == -1)
953 {
954 // Found unassigned face. Use its
955 // information.
956 // Note that we use the lowest numbered
957 // patch face.
958
959 meshEdgeToFace[edgei] =
960 globalFaces.toGlobal(facei);
961 }
962
963 // Override any patch info. Note that
964 // meshEdgeToFace might be an internal face.
965 if (meshEdgeToPatch[edgei] == -1)
966 {
967 meshEdgeToPatch[edgei] = pp.index();
968 }
969
970 // Override any zone info
971 if (meshEdgeToZone[edgei] == -1)
972 {
973 meshEdgeToZone[edgei] =
974 faceToZone[facei];
975 const edge& meshE = mesh.edges()[edgei];
976 const int d = edge::compare(e, meshE);
977 if (d == 1)
978 {
979 meshEdgeToFlip[edgei] =
980 faceToFlip[facei];
981 }
982 else if (d == -1)
983 {
984 meshEdgeToFlip[edgei] =
985 !faceToFlip[facei];
986 }
987 else
988 {
990 << "big problem"
991 << exit(FatalError);
992 }
993 }
994 }
995
996 prevPointi = pointi;
997 }
998 }
999 }
1000 }
1001 }
1002 }
1003
1004
1005 // Synchronise across coupled edges. Max patch/face/faceZone wins
1007 (
1008 mesh,
1009 meshEdgeToFace,
1010 maxEqOp<label>(),
1011 label(-1)
1012 );
1014 (
1015 mesh,
1016 meshEdgeToPatch,
1017 maxEqOp<label>(),
1018 label(-1)
1019 );
1021 (
1022 mesh,
1023 meshEdgeToZone,
1024 maxEqOp<label>(),
1025 label(-1)
1026 );
1027// // Note: flipMap not yet done. Needs edge orientation. This is handled
1028// // later on.
1029// if (Pstream::parRun())
1030// {
1031// const globalMeshData& gd = mesh.globalData();
1032// const indirectPrimitivePatch& cpp = gd.coupledPatch();
1033//
1034// labelList patchEdges;
1035// labelList coupledEdges;
1036// bitSet sameEdgeOrientation;
1037// PatchTools::matchEdges
1038// (
1039// pp,
1040// cpp,
1041// patchEdges,
1042// coupledEdges,
1043// sameEdgeOrientation
1044// );
1045//
1046// // Convert data on pp edges to data on coupled patch
1047// labelList cppEdgeZoneID(cpp.nEdges(), -1);
1048// boolList cppEdgeFlip(cpp.nEdges(), false);
1049// forAll(coupledEdges, i)
1050// {
1051// label cppEdgei = coupledEdges[i];
1052// label ppEdgei = patchEdges[i];
1053//
1054// cppEdgeZoneID[cppEdgei] = edgeZoneID[ppEdgei];
1055// if (sameEdgeOrientation[i])
1056// {
1057// cppEdgeFlip[cppEdgei] = edgeFlip[ppEdgei];
1058// }
1059// else
1060// {
1061// cppEdgeFlip[cppEdgei] = !edgeFlip[ppEdgei];
1062// }
1063// }
1064//
1065// // Sync
1066// const globalIndexAndTransform& git = gd.globalTransforms();
1067// const mapDistribute& edgeMap = gd.globalEdgeSlavesMap();
1068//
1069// globalMeshData::syncData
1070// (
1071// cppEdgeZoneID,
1072// gd.globalEdgeSlaves(),
1073// gd.globalEdgeTransformedSlaves(),
1074// edgeMap,
1075// git,
1076// maxEqOp<label>(),
1077// dummyTransform()
1078// );
1079// globalMeshData::syncData
1080// (
1081// cppEdgeFlip,
1082// gd.globalEdgeSlaves(),
1083// gd.globalEdgeTransformedSlaves(),
1084// edgeMap,
1085// git,
1086// andEqOp<bool>(),
1087// dummyTransform()
1088// );
1089//
1090// // Convert data on coupled edges to pp edges
1091// forAll(coupledEdges, i)
1092// {
1093// label cppEdgei = coupledEdges[i];
1094// label ppEdgei = patchEdges[i];
1095//
1096// edgeZoneID[ppEdgei] = cppEdgeZoneID[cppEdgei];
1097// if (sameEdgeOrientation[i])
1098// {
1099// edgeFlip[ppEdgei] = cppEdgeFlip[cppEdgei];
1100// }
1101// else
1102// {
1103// edgeFlip[ppEdgei] = !cppEdgeFlip[cppEdgei];
1104// }
1105// }
1106// }
1107
1109 (
1110 mesh,
1111 meshEdgeToFlip,
1112 orEqOp<unsigned int>(),
1113 0
1114 );
1115
1116 // Extract pp info
1117 patchEdgeToFace = UIndirectList<label>(meshEdgeToFace, meshEdges);
1118 patchEdgeToPatch = UIndirectList<label>(meshEdgeToPatch, meshEdges);
1119 patchEdgeToZone = UIndirectList<label>(meshEdgeToZone, meshEdges);
1120 patchEdgeToFlip.setSize(meshEdges.size());
1121 patchEdgeToFlip = false;
1122 forAll(meshEdges, i)
1123 {
1124 patchEdgeToFlip[i] = meshEdgeToFlip[meshEdges[i]];
1125 }
1126}
1127
1128
1130(
1131 const bool zoneFromAnyFace,
1132
1133 const polyMesh& mesh,
1134 const globalIndex& globalFaces,
1135 const labelListList& globalEdgeFaces,
1136 const indirectPrimitivePatch& pp,
1137
1138 labelList& edgePatchID,
1139 label& nPatches,
1140 Map<label>& nbrProcToPatch,
1141 Map<label>& patchToNbrProc,
1142 labelList& edgeZoneID,
1143 boolList& edgeFlip,
1144 labelList& inflateFaceID
1145)
1146{
1148 const globalMeshData& gd = mesh.globalData();
1149
1150 // Precalculate mesh edges for pp.edges.
1151 const labelList meshEdges(pp.meshEdges(mesh.edges(), mesh.pointEdges()));
1152
1153 edgePatchID.setSize(pp.nEdges());
1154 edgePatchID = -1;
1155 nPatches = patches.size();
1156 edgeZoneID.setSize(pp.nEdges());
1157 edgeZoneID = -1;
1158 edgeFlip.setSize(pp.nEdges());
1159 edgeFlip = false;
1160 inflateFaceID.setSize(pp.nEdges(), -1);
1161
1162
1163 // Determine properties for faces extruded from edges
1164 // - edge inbetween two different processors:
1165 // - extrude as patch face on correct processor
1166 // - edge at end of patch (so edgeFaces.size() == 1):
1167 // - use mesh face that is a boundary face
1168 // - edge internal to patch (so edgeFaces.size() == 2):
1169
1170
1171 // Pass1:
1172 // For all edges inbetween two processors: see if matches to existing
1173 // processor patch or create interprocessor-patch if necessary.
1174 // Sets edgePatchID[edgeI] but none of the other quantities
1175
1176 forAll(globalEdgeFaces, edgei)
1177 {
1178 const labelList& eGlobalFaces = globalEdgeFaces[edgei];
1179 if
1180 (
1181 eGlobalFaces.size() == 2
1182 && pp.edgeFaces()[edgei].size() == 1
1183 )
1184 {
1185 // Locally but not globally a boundary edge. Hence a coupled
1186 // edge. Find the patch to use if on different processors.
1187
1188 label f0 = eGlobalFaces[0];
1189 label f1 = eGlobalFaces[1];
1190
1191 label otherProci = -1;
1192 if (globalFaces.isLocal(f0) && !globalFaces.isLocal(f1))
1193 {
1194 otherProci = globalFaces.whichProcID(f1);
1195 }
1196 else if (!globalFaces.isLocal(f0) && globalFaces.isLocal(f1))
1197 {
1198 otherProci = globalFaces.whichProcID(f0);
1199 }
1200
1201
1202 if (otherProci != -1)
1203 {
1204 // Use existing processorPolyPatch to otherProci?
1205
1206 label procPatchi =
1207 gd.topology().procPatchLookup(otherProci);
1208
1209 if (procPatchi < 0)
1210 {
1211 // No existing processorPolyPatch to otherProci.
1212 // See if already marked for addition
1213 procPatchi = nbrProcToPatch.lookup(otherProci, -1);
1214
1215 if (procPatchi < 0)
1216 {
1217 // Add new proc-patch, mark for addition.
1218
1219 procPatchi = nPatches;
1220 ++nPatches;
1221
1222 nbrProcToPatch.insert(otherProci, procPatchi);
1223 patchToNbrProc.insert(procPatchi, otherProci);
1224 }
1225 }
1226
1227 edgePatchID[edgei] = procPatchi;
1228 }
1229 }
1230 }
1231
1232
1233 // Pass2: determine face properties for all other edges
1234 // ----------------------------------------------------
1235
1236 // Info for edges of pp
1237 labelList edgeToFace;
1238 labelList edgeToPatch;
1239 labelList edgeToZone;
1240 bitSet edgeToFlip;
1241 globalEdgeInfo
1242 (
1243 zoneFromAnyFace, // internal edge info also from boundary faces
1244
1245 mesh,
1246 globalFaces,
1247 globalEdgeFaces,
1248 pp,
1249 meshEdges,
1250
1251 edgeToFace, // face (in globalFaces index)
1252 edgeToPatch, // patch on face (or -1 for internal faces)
1253 edgeToZone, // zone on face
1254 edgeToFlip // flip orientation on face
1255 );
1256
1257 const labelListList& edgeFaces = pp.edgeFaces();
1258
1259 DynamicList<label> dynMeshEdgeFaces;
1260
1261 forAll(edgeFaces, edgei)
1262 {
1263 if (edgePatchID[edgei] == -1)
1264 {
1265 if (edgeFaces[edgei].size() == 2)
1266 {
1267 // Internal edge. Look at any face (internal or boundary) to
1268 // determine extrusion properties. First one that has zone
1269 // info wins
1270 if (globalFaces.isLocal(edgeToFace[edgei]))
1271 {
1272 inflateFaceID[edgei] =
1273 globalFaces.toLocal(edgeToFace[edgei]);
1274 }
1275 edgeZoneID[edgei] = edgeToZone[edgei];
1276 edgeFlip[edgei] = edgeToFlip[edgei];
1277 }
1278 else
1279 {
1280 // Proper, uncoupled patch edge. Boundary to get info from
1281 // might be on a different processor!
1282
1283 if (globalFaces.isLocal(edgeToFace[edgei]))
1284 {
1285 inflateFaceID[edgei] =
1286 globalFaces.toLocal(edgeToFace[edgei]);
1287 }
1288 edgePatchID[edgei] = edgeToPatch[edgei];
1289 edgeZoneID[edgei] = edgeToZone[edgei];
1290 edgeFlip[edgei] = edgeToFlip[edgei];
1291 }
1292 }
1293 }
1294
1295
1296
1297 // Now hopefully every boundary edge has a edge patch. Check
1298 if (debug)
1299 {
1300 forAll(edgeFaces, edgei)
1301 {
1302 if (edgeFaces[edgei].size() == 1 && edgePatchID[edgei] == -1)
1303 {
1304 const edge& e = pp.edges()[edgei];
1306 << "Have no sidePatchID for edge " << edgei << " points "
1307 << pp.points()[pp.meshPoints()[e[0]]]
1308 << pp.points()[pp.meshPoints()[e[1]]]
1309 << endl;
1310 }
1311 }
1312 }
1313
1314
1315
1316 // Pass3: for any faces set in pass1 see if we can find a processor face
1317 // to inherit from (we only have a patch, not a patch face)
1318 forAll(edgeFaces, edgei)
1319 {
1320 if
1321 (
1322 edgeFaces[edgei].size() == 1
1323 && globalEdgeFaces[edgei].size() == 2
1324 && edgePatchID[edgei] != -1
1325 && inflateFaceID[edgei] == -1
1326 )
1327 {
1328 // 1. Do we have a local boundary face to inflate from
1329
1330 label myFaceI = pp.addressing()[edgeFaces[edgei][0]];
1331
1332 // Pick up any boundary face on this edge and use its properties
1333 label meshEdgei = meshEdges[edgei];
1334 const labelList& meshFaces = mesh.edgeFaces
1335 (
1336 meshEdgei,
1337 dynMeshEdgeFaces
1338 );
1339
1340 forAll(meshFaces, k)
1341 {
1342 label facei = meshFaces[k];
1343
1344 if (facei != myFaceI && !mesh.isInternalFace(facei))
1345 {
1346 if (patches.whichPatch(facei) == edgePatchID[edgei])
1347 {
1348 setFaceProps
1349 (
1350 mesh,
1351 pp,
1352 edgei,
1353 facei,
1354
1355 edgePatchID[edgei],
1356 edgeZoneID[edgei],
1357 edgeFlip[edgei],
1358 inflateFaceID[edgei]
1359 );
1360 break;
1361 }
1362 }
1363 }
1364 }
1365 }
1366
1367
1368
1369 // Sync all data:
1370 // - edgePatchID : might be local in case of processor patch. Do not
1371 // sync for now
1372 // - inflateFaceID: local. Do not sync
1373 // - edgeZoneID : global. sync.
1374 // - edgeFlip : global. sync.
1375
1376 if (Pstream::parRun())
1377 {
1378 const globalMeshData& gd = mesh.globalData();
1379 const indirectPrimitivePatch& cpp = gd.coupledPatch();
1380
1381 labelList patchEdges;
1382 labelList coupledEdges;
1383 bitSet sameEdgeOrientation;
1385 (
1386 pp,
1387 cpp,
1388 patchEdges,
1389 coupledEdges,
1390 sameEdgeOrientation
1391 );
1392
1393 // Convert data on pp edges to data on coupled patch
1394 labelList cppEdgeZoneID(cpp.nEdges(), -1);
1395 boolList cppEdgeFlip(cpp.nEdges(), false);
1396 forAll(coupledEdges, i)
1397 {
1398 label cppEdgei = coupledEdges[i];
1399 label ppEdgei = patchEdges[i];
1400
1401 cppEdgeZoneID[cppEdgei] = edgeZoneID[ppEdgei];
1402 if (sameEdgeOrientation[i])
1403 {
1404 cppEdgeFlip[cppEdgei] = edgeFlip[ppEdgei];
1405 }
1406 else
1407 {
1408 cppEdgeFlip[cppEdgei] = !edgeFlip[ppEdgei];
1409 }
1410 }
1411
1412 // Sync
1413 const globalIndexAndTransform& git = gd.globalTransforms();
1414 const mapDistribute& edgeMap = gd.globalEdgeSlavesMap();
1415
1417 (
1418 cppEdgeZoneID,
1419 gd.globalEdgeSlaves(),
1421 edgeMap,
1422 git,
1425 );
1427 (
1428 cppEdgeFlip,
1429 gd.globalEdgeSlaves(),
1431 edgeMap,
1432 git,
1433 andEqOp<bool>(),
1435 );
1436
1437 // Convert data on coupled edges to pp edges
1438 forAll(coupledEdges, i)
1439 {
1440 label cppEdgei = coupledEdges[i];
1441 label ppEdgei = patchEdges[i];
1442
1443 edgeZoneID[ppEdgei] = cppEdgeZoneID[cppEdgei];
1444 if (sameEdgeOrientation[i])
1445 {
1446 edgeFlip[ppEdgei] = cppEdgeFlip[cppEdgei];
1447 }
1448 else
1449 {
1450 edgeFlip[ppEdgei] = !cppEdgeFlip[cppEdgei];
1451 }
1452 }
1453 }
1454}
1455
1456
1458(
1459 const globalIndex& globalFaces,
1460 const labelListList& globalEdgeFaces,
1461 const scalarField& expansionRatio,
1462 const indirectPrimitivePatch& pp,
1463 const bitSet& ppFlip,
1464
1465 const labelList& edgePatchID,
1466 const labelList& edgeZoneID,
1467 const boolList& edgeFlip,
1468 const labelList& inflateFaceID,
1469
1470 const labelList& exposedPatchID,
1471 const labelList& nFaceLayers,
1472 const labelList& nPointLayers,
1473 const vectorField& firstLayerDisp,
1474 polyTopoChange& meshMod
1475)
1476{
1477 if (debug)
1478 {
1479 Pout<< "addPatchCellLayer::setRefinement : Adding up to "
1480 << gMax(nPointLayers)
1481 << " layers of cells to indirectPrimitivePatch with "
1482 << pp.nPoints() << " points" << endl;
1483 }
1484
1485 if
1486 (
1487 pp.nPoints() != firstLayerDisp.size()
1488 || pp.nPoints() != nPointLayers.size()
1489 || pp.size() != nFaceLayers.size()
1490 || pp.size() != ppFlip.size()
1491 )
1492 {
1494 << "Size of new points is not same as number of points used by"
1495 << " the face subset" << endl
1496 << " patch.nPoints:" << pp.nPoints()
1497 << " displacement:" << firstLayerDisp.size()
1498 << " nPointLayers:" << nPointLayers.size() << nl
1499 << " patch.nFaces:" << pp.size()
1500 << " flip map:" << ppFlip.size()
1501 << " nFaceLayers:" << nFaceLayers.size()
1502 << abort(FatalError);
1503 }
1504 if (!addToMesh_)
1505 {
1506 // flip map should be false
1507 if (ppFlip.count())
1508 {
1510 << "In generating stand-alone mesh the flip map should be empty"
1511 << ". Instead it is " << ppFlip.count()
1512 << abort(FatalError);
1513 }
1514 }
1515 else
1516 {
1517 // Maybe check for adding to neighbour of boundary faces? How about
1518 // coupled faces where the faceZone flipMap is negated
1519
1520 // For all boundary faces:
1521 // -1 : not extruded
1522 // 0 : extruded from owner outwards (flip = false)
1523 // 1 : extrude from neighbour outwards
1524 labelList stateAndFlip(mesh_.nBoundaryFaces(), 0);
1525 forAll(pp.addressing(), patchFacei)
1526 {
1527 if (nFaceLayers[patchFacei] > 0)
1528 {
1529 const label meshFacei = pp.addressing()[patchFacei];
1530 const label bFacei = meshFacei-mesh_.nInternalFaces();
1531 if (bFacei >= 0)
1532 {
1533 stateAndFlip[bFacei] = label(ppFlip[patchFacei]);
1534 }
1535 }
1536 }
1537 // Make sure uncoupled patches do not trigger the error below
1538 for (const auto& patch : mesh_.boundaryMesh())
1539 {
1540 if (!patch.coupled())
1541 {
1542 forAll(patch, i)
1543 {
1544 label& state = stateAndFlip[patch.offset()+i];
1545 state = (state == 0 ? 1 : 0);
1546 }
1547 }
1548 }
1549 syncTools::swapBoundaryFaceList(mesh_, stateAndFlip);
1550
1551 forAll(pp.addressing(), patchFacei)
1552 {
1553 if (nFaceLayers[patchFacei] > 0)
1554 {
1555 const label meshFacei = pp.addressing()[patchFacei];
1556 const label bFacei = meshFacei-mesh_.nInternalFaces();
1557 if (bFacei >= 0)
1558 {
1559 if (stateAndFlip[bFacei] == -1)
1560 {
1562 << "At extruded face:" << meshFacei
1563 << " at:" << mesh_.faceCentres()[meshFacei]
1564 << " locally have nLayers:"
1565 << nFaceLayers[patchFacei]
1566 << " but remotely have none" << exit(FatalError);
1567 }
1568 else if (stateAndFlip[bFacei] == label(ppFlip[patchFacei]))
1569 {
1571 << "At extruded face:" << meshFacei
1572 << " at:" << mesh_.faceCentres()[meshFacei]
1573 << " locally have flip:" << ppFlip[patchFacei]
1574 << " which is not the opposite of coupled version "
1575 << stateAndFlip[bFacei]
1576 << exit(FatalError);
1577 }
1578 }
1579 }
1580 }
1581 }
1582
1583
1584 forAll(nPointLayers, i)
1585 {
1586 if (nPointLayers[i] < 0)
1587 {
1589 << "Illegal number of layers " << nPointLayers[i]
1590 << " at patch point " << i << abort(FatalError);
1591 }
1592 }
1593 forAll(nFaceLayers, i)
1594 {
1595 if (nFaceLayers[i] < 0)
1596 {
1598 << "Illegal number of layers " << nFaceLayers[i]
1599 << " at patch face " << i << abort(FatalError);
1600 }
1601 }
1602
1603 forAll(globalEdgeFaces, edgei)
1604 {
1605 if (globalEdgeFaces[edgei].size() > 2)
1606 {
1607 const edge& e = pp.edges()[edgei];
1608
1609 if (nPointLayers[e[0]] > 0 || nPointLayers[e[1]] > 0)
1610 {
1612 << "Trying to extrude edge "
1613 << e.line(pp.localPoints())
1614 << " which is non-manifold (has "
1615 << globalEdgeFaces[edgei].size()
1616 << " local or coupled faces using it)"
1617 << " of which "
1618 << pp.edgeFaces()[edgei].size()
1619 << " local"
1620 << abort(FatalError);
1621 }
1622 }
1623 }
1624
1625
1626 const labelList& meshPoints = pp.meshPoints();
1627
1628
1629 // Determine which points are on which side of the extrusion
1630 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1631
1632 bitSet isBlockedFace(mesh_.nFaces());
1633 forAll(nFaceLayers, patchFacei)
1634 {
1635 if (nFaceLayers[patchFacei] > 0)
1636 {
1637 isBlockedFace.set(pp.addressing()[patchFacei]);
1638 }
1639 }
1640
1641 // Some storage for edge-face-addressing.
1643
1644 // Precalculate mesh edges for pp.edges.
1645 const labelList meshEdges(pp.meshEdges(mesh_.edges(), mesh_.pointEdges()));
1646
1647 if (debug)
1648 {
1649 // Check synchronisation
1650 // ~~~~~~~~~~~~~~~~~~~~~
1651
1652 {
1653 labelList n(mesh_.nPoints(), Zero);
1654 labelUIndList(n, meshPoints) = nPointLayers;
1655 syncTools::syncPointList(mesh_, n, maxEqOp<label>(), label(0));
1656
1657 // Non-synced
1658 forAll(meshPoints, i)
1659 {
1660 label meshPointi = meshPoints[i];
1661
1662 if (n[meshPointi] != nPointLayers[i])
1663 {
1665 << "At mesh point:" << meshPointi
1666 << " coordinate:" << mesh_.points()[meshPointi]
1667 << " specified nLayers:" << nPointLayers[i] << endl
1668 << "On coupled point a different nLayers:"
1669 << n[meshPointi] << " was specified."
1670 << abort(FatalError);
1671 }
1672 }
1673
1674
1675 // Check that nPointLayers equals the max layers of connected faces
1676 // (or 0). Anything else makes no sense.
1677 labelList nFromFace(mesh_.nPoints(), Zero);
1678 forAll(nFaceLayers, i)
1679 {
1680 const face& f = pp[i];
1681
1682 forAll(f, fp)
1683 {
1684 label pointi = f[fp];
1685
1686 nFromFace[pointi] = max(nFromFace[pointi], nFaceLayers[i]);
1687 }
1688 }
1690 (
1691 mesh_,
1692 nFromFace,
1694 label(0)
1695 );
1696
1697 forAll(nPointLayers, i)
1698 {
1699 label meshPointi = meshPoints[i];
1700
1701 if
1702 (
1703 nPointLayers[i] > 0
1704 && nPointLayers[i] != nFromFace[meshPointi]
1705 )
1706 {
1708 << "At mesh point:" << meshPointi
1709 << " coordinate:" << mesh_.points()[meshPointi]
1710 << " specified nLayers:" << nPointLayers[i] << endl
1711 << "but the max nLayers of surrounding faces is:"
1712 << nFromFace[meshPointi]
1713 << abort(FatalError);
1714 }
1715 }
1716 }
1717
1718 {
1719 pointField d(mesh_.nPoints(), vector::max);
1720 UIndirectList<point>(d, meshPoints) = firstLayerDisp;
1722 (
1723 mesh_,
1724 d,
1727 );
1728
1729 forAll(meshPoints, i)
1730 {
1731 label meshPointi = meshPoints[i];
1732
1733 if (mag(d[meshPointi] - firstLayerDisp[i]) > SMALL)
1734 {
1736 << "At mesh point:" << meshPointi
1737 << " coordinate:" << mesh_.points()[meshPointi]
1738 << " specified displacement:" << firstLayerDisp[i]
1739 << endl
1740 << "On coupled point a different displacement:"
1741 << d[meshPointi] << " was specified."
1742 << abort(FatalError);
1743 }
1744 }
1745 }
1746
1747 // Check that edges of pp (so ones that become boundary faces)
1748 // connect to only one boundary face. Guarantees uniqueness of
1749 // patch that they go into so if this is a coupled patch both
1750 // sides decide the same.
1751 // ~~~~~~~~~~~~~~~~~~~~~~
1752
1753 for (label edgei = pp.nInternalEdges(); edgei < pp.nEdges(); edgei++)
1754 {
1755 const edge& e = pp.edges()[edgei];
1756
1757 if (nPointLayers[e[0]] > 0 || nPointLayers[e[1]] > 0)
1758 {
1759 // Edge is to become a face
1760
1761 const labelList& eFaces = pp.edgeFaces()[edgei];
1762
1763 // First check: pp should be single connected.
1764 if (eFaces.size() != 1)
1765 {
1767 << "boundary-edge-to-be-extruded:"
1768 << pp.points()[meshPoints[e[0]]]
1769 << pp.points()[meshPoints[e[1]]]
1770 << " has more than two faces using it:" << eFaces
1771 << abort(FatalError);
1772 }
1773
1774 //label myFacei = pp.addressing()[eFaces[0]];
1775 //
1776 //label meshEdgei = meshEdges[edgei];
1777 //
1779 //const labelList& meshFaces = mesh_.edgeFaces(meshEdgei, ef);
1780 //
1782 //const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1783 //
1784 //label bFacei = -1;
1785 //
1786 //forAll(meshFaces, i)
1787 //{
1788 // label facei = meshFaces[i];
1789 //
1790 // if (facei != myFacei)
1791 // {
1792 // if (!mesh_.isInternalFace(facei))
1793 // {
1794 // if (bFacei == -1)
1795 // {
1796 // bFacei = facei;
1797 // }
1798 // else
1799 // {
1800 // //FatalErrorInFunction
1801 // WarningInFunction
1802 // << "boundary-edge-to-be-extruded:"
1803 // << pp.points()[meshPoints[e[0]]]
1804 // << pp.points()[meshPoints[e[1]]]
1805 // << " has more than one boundary faces"
1806 // << " using it:"
1807 // << bFacei << " fc:"
1808 // << mesh_.faceCentres()[bFacei]
1809 // << " patch:" << patches.whichPatch(bFacei)
1810 // << " and " << facei << " fc:"
1811 // << mesh_.faceCentres()[facei]
1812 // << " patch:" << patches.whichPatch(facei)
1813 // << endl;
1814 // //abort(FatalError);
1815 // }
1816 // }
1817 // }
1818 //}
1819 }
1820 }
1821 }
1822
1823
1824 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1825
1826 // Precalculated patchID for each patch face
1827 labelList patchID(pp.size());
1828
1829 forAll(pp, patchFacei)
1830 {
1831 label meshFacei = pp.addressing()[patchFacei];
1832
1833 patchID[patchFacei] = patches.whichPatch(meshFacei);
1834 }
1835
1836
1837 // From master point (in patch point label) to added points (in mesh point
1838 // label)
1839 addedPoints_.setSize(pp.nPoints());
1840
1841 // Mark points that do not get extruded by setting size of addedPoints_ to 0
1842 label nTruncated = 0;
1843
1844 forAll(nPointLayers, patchPointi)
1845 {
1846 if (nPointLayers[patchPointi] > 0)
1847 {
1848 addedPoints_[patchPointi].setSize(nPointLayers[patchPointi]);
1849 }
1850 else
1851 {
1852 nTruncated++;
1853 }
1854 }
1855
1856 if (debug)
1857 {
1858 Pout<< "Not adding points at " << nTruncated << " out of "
1859 << pp.nPoints() << " points" << endl;
1860 }
1861
1862
1863 // Store per face whether it uses the duplicated point or the original one
1864 // Also mark any affected cells. We could transport the duplicated point
1865 // itself but since it is a processor-local index only we only transport
1866 // a boolean.
1867
1868 // Per face, per index in face either labelMax or a valid index. Note:
1869 // most faces are not affected in which case the face will be zero size
1870 // and only have a nullptr and a size.
1871 faceList baseFaces(mesh_.nFaces());
1872 bitSet isAffectedCell(mesh_.nCells());
1873 {
1874 const faceList& localFaces = pp.localFaces();
1875 forAll(localFaces, patchFacei)
1876 {
1877 const face& f = localFaces[patchFacei];
1878 forAll(f, fp)
1879 {
1880 const label patchPointi = f[fp];
1881 if (nPointLayers[patchPointi] > 0)
1882 {
1883 const label meshFacei = pp.addressing()[patchFacei];
1884 face& baseF = baseFaces[meshFacei];
1885 // Initialise to labelMax if not yet sized
1886 baseF.setSize(f.size(), labelMax);
1887 baseF[fp] = pp.meshPoints()[patchPointi];
1888
1889 if (ppFlip[patchFacei])
1890 {
1891 // Neighbour stays. Affected points on the owner side.
1892 const label celli = mesh_.faceOwner()[meshFacei];
1893 isAffectedCell.set(celli);
1894 }
1895 else if (mesh_.isInternalFace(meshFacei))
1896 {
1897 // Owner unaffected. Unaffected points on neighbour side
1898 const label celli = mesh_.faceNeighbour()[meshFacei];
1899 isAffectedCell.set(celli);
1900 }
1901 }
1902 }
1903 }
1904 }
1905
1906 // Transport affected side across faces. Could do across edges: say we have
1907 // a loose cell edge-(but not face-)connected to face-to-be-extruded do
1908 // we want it to move with the extrusion or stay connected to the original?
1909 // For now just keep it connected to the original.
1910 {
1911 // Work space
1912 Map<label> minPointValue(128);
1913 faceList oldBoundaryFaces(mesh_.nBoundaryFaces());
1914
1915 while (true)
1916 {
1917 bitSet newIsAffectedCell(mesh_.nCells());
1918
1919 label nChanged = 0;
1920 for (const label celli : isAffectedCell)
1921 {
1922 const cell& cFaces = mesh_.cells()[celli];
1923
1924 // 1. Determine marked base points. Inside a single cell all
1925 // faces use the same 'instance' of a point.
1926 minPointValue.clear();
1927 for (const label facei : cFaces)
1928 {
1929 const face& baseF = baseFaces[facei];
1930 const face& f = mesh_.faces()[facei];
1931
1932 if (baseF.size())
1933 {
1934 forAll(f, fp)
1935 {
1936 if (baseF[fp] != labelMax)
1937 {
1938 // Could check here for inconsistent patchPoint
1939 // e.g. cell using both sides of a
1940 // face-to-be-extruded. Is not possible!
1941 minPointValue.insert(f[fp], baseF[fp]);
1942 }
1943 }
1944 }
1945 }
1946
1947 //Pout<< "For cell:" << celli
1948 // << " at:" << mesh_.cellCentres()[celli]
1949 // << " have minPointValue:" << minPointValue
1950 // << endl;
1951
1952 // 2. Transport marked points on all cell points
1953 for (const label facei : cFaces)
1954 {
1955 const face& f = mesh_.faces()[facei];
1956 face& baseF = baseFaces[facei];
1957
1958 const label oldNChanged = nChanged;
1959 forAll(f, fp)
1960 {
1961 const auto fnd = minPointValue.find(f[fp]);
1962 if (fnd.found())
1963 {
1964 baseF.setSize(f.size(), labelMax);
1965 if (baseF[fp] == labelMax)
1966 {
1967 baseF[fp] = fnd();
1968 nChanged++;
1969
1970 //Pout<< "For cell:" << celli
1971 // << " at:" << mesh_.cellCentres()[celli]
1972 // << " on face:" << facei
1973 // << " points:"
1974 // << UIndirectList<point>(mesh_.points(), f)
1975 // << " now have baseFace:" << baseF
1976 // << endl;
1977 }
1978 }
1979 }
1980
1981 if (!isBlockedFace(facei) && nChanged > oldNChanged)
1982 {
1983 // Mark neighbouring cells
1984 const label own = mesh_.faceOwner()[facei];
1985 if (!isAffectedCell[own])
1986 {
1987 newIsAffectedCell.set(own);
1988 }
1989 if (mesh_.isInternalFace(facei))
1990 {
1991 const label nei = mesh_.faceNeighbour()[facei];
1992 if (!isAffectedCell[nei])
1993 {
1994 newIsAffectedCell.set(nei);
1995 }
1996 }
1997 }
1998 }
1999 }
2000
2001 if (debug)
2002 {
2003 Pout<< "isAffectedCell:" << isAffectedCell.count() << endl;
2004 Pout<< "newIsAffectedCell:" << newIsAffectedCell.count()
2005 << endl;
2006 Pout<< "nChanged:" << nChanged << endl;
2007 }
2008
2009 if (returnReduce(nChanged, sumOp<label>()) == 0)
2010 {
2011 break;
2012 }
2013
2014
2015 // Transport minimum across coupled faces
2016 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2017
2019 (
2020 baseFaces,
2021 mesh_.nBoundaryFaces(),
2022 mesh_.nInternalFaces()
2023 );
2024 oldBoundaryFaces = l;
2026 (
2027 mesh_,
2028 l,
2029 minEqOpFace(),
2030 Foam::dummyTransform() // dummy transformation
2031 );
2032
2033 forAll(l, bFacei)
2034 {
2035 // Note: avoid special handling of comparing zero-sized faces
2036 // (see face::operator==). Review.
2037 const labelUList& baseVts = l[bFacei];
2038 const labelUList& oldVts = oldBoundaryFaces[bFacei];
2039 if (baseVts != oldVts)
2040 {
2041 const label facei = mesh_.nInternalFaces()+bFacei;
2042 const label own = mesh_.faceOwner()[facei];
2043 if (!isAffectedCell[own])
2044 {
2045 newIsAffectedCell.set(own);
2046 }
2047 }
2048 }
2049
2050 isAffectedCell = newIsAffectedCell;
2051 }
2052 }
2053
2054
2055 //
2056 // Create new points
2057 //
2058
2059 // If creating new mesh: copy existing patch points
2060 labelList copiedPatchPoints;
2061 if (!addToMesh_)
2062 {
2063 copiedPatchPoints.setSize(firstLayerDisp.size());
2064 forAll(firstLayerDisp, patchPointi)
2065 {
2066 if (addedPoints_[patchPointi].size())
2067 {
2068 label meshPointi = meshPoints[patchPointi];
2069 label zoneI = mesh_.pointZones().whichZone(meshPointi);
2070 copiedPatchPoints[patchPointi] = meshMod.setAction
2071 (
2073 (
2074 mesh_.points()[meshPointi], // point
2075 -1, // master point
2076 zoneI, // zone for point
2077 true // supports a cell
2078 )
2079 );
2080 }
2081 }
2082 }
2083
2084
2085 // Create points for additional layers
2086 forAll(firstLayerDisp, patchPointi)
2087 {
2088 if (addedPoints_[patchPointi].size())
2089 {
2090 const label meshPointi = meshPoints[patchPointi];
2091 const label zoneI = mesh_.pointZones().whichZone(meshPointi);
2092
2093 point pt = mesh_.points()[meshPointi];
2094
2095 vector disp = firstLayerDisp[patchPointi];
2096
2097 forAll(addedPoints_[patchPointi], i)
2098 {
2099 pt += disp;
2100
2101 const label addedVertI = meshMod.setAction
2102 (
2104 (
2105 pt, // point
2106 (addToMesh_ ? meshPointi : -1), // master point
2107 zoneI, // zone for point
2108 true // supports a cell
2109 )
2110 );
2111
2112
2113 //Pout<< "Adding point:" << addedVertI << " at:" << pt
2114 // << " from point:" << meshPointi
2115 // << " at:" << mesh_.points()[meshPointi]
2116 // << endl;
2117
2118 addedPoints_[patchPointi][i] = addedVertI;
2119
2120 disp *= expansionRatio[patchPointi];
2121 }
2122 }
2123 }
2124
2125
2126 //
2127 // Add cells to all boundaryFaces
2128 //
2129
2130 labelListList addedCells(pp.size());
2131
2132 forAll(pp, patchFacei)
2133 {
2134 if (nFaceLayers[patchFacei] > 0)
2135 {
2136 const label meshFacei = pp.addressing()[patchFacei];
2137
2138 label extrudeCelli = -2;
2139 label extrudeZonei;
2140 if (!addToMesh_)
2141 {
2142 extrudeCelli = -1;
2143 const label ownCelli = mesh_.faceOwner()[meshFacei];
2144 extrudeZonei = mesh_.cellZones().whichZone(ownCelli);
2145 }
2146 else if (!ppFlip[patchFacei])
2147 {
2148 // Normal: extrude from owner face
2149 extrudeCelli = mesh_.faceOwner()[meshFacei];
2150 extrudeZonei = mesh_.cellZones().whichZone(extrudeCelli);
2151 }
2152 else if (mesh_.isInternalFace(meshFacei))
2153 {
2154 // Extrude from neighbour face (if internal). Might be
2155 // that it is a coupled face and the other side is
2156 // extruded
2157 extrudeCelli = mesh_.faceNeighbour()[meshFacei];
2158 extrudeZonei = mesh_.cellZones().whichZone(extrudeCelli);
2159 }
2160
2161 if (extrudeCelli != -2)
2162 {
2163 addedCells[patchFacei].setSize(nFaceLayers[patchFacei]);
2164
2165 for (label i = 0; i < nFaceLayers[patchFacei]; i++)
2166 {
2167 // Note: add from cell (owner of patch face) or from face?
2168 // for now add from cell so we can map easily.
2169 addedCells[patchFacei][i] = meshMod.setAction
2170 (
2172 (
2173 -1, // master point
2174 -1, // master edge
2175 -1, // master face
2176 extrudeCelli, // master cell
2177 extrudeZonei // zone for cell
2178 )
2179 );
2180
2181 //Pout<< "Added cell:" << addedCells[patchFacei][i]
2182 // << " from master:" << extrudeCelli
2183 // << " at:" << mesh_.cellCentres()[extrudeCelli]
2184 // << endl;
2185 }
2186 }
2187 }
2188 }
2189
2190
2191
2192 // Create faces on top of the original patch faces.
2193 // These faces are created from original patch faces outwards so in order
2194 // of increasing cell number. So orientation should be same as original
2195 // patch face for them to have owner<neighbour.
2196
2197 layerFaces_.setSize(pp.size());
2198
2199 forAll(pp.localFaces(), patchFacei)
2200 {
2201 label meshFacei = pp.addressing()[patchFacei];
2202
2203 if (addedCells[patchFacei].size())
2204 {
2205 layerFaces_[patchFacei].setSize(addedCells[patchFacei].size() + 1);
2206
2207 // Get duplicated vertices on the patch face.
2208 const face& f = pp.localFaces()[patchFacei];
2209
2210 face newFace(f.size());
2211
2212 forAll(addedCells[patchFacei], i)
2213 {
2214 forAll(f, fp)
2215 {
2216 if (addedPoints_[f[fp]].empty())
2217 {
2218 // Keep original point
2219 newFace[fp] =
2220 (
2221 addToMesh_
2222 ? meshPoints[f[fp]]
2223 : copiedPatchPoints[f[fp]]
2224 );
2225 }
2226 else
2227 {
2228 // Get new outside point
2229 label offset =
2230 addedPoints_[f[fp]].size()
2231 - addedCells[patchFacei].size();
2232 newFace[fp] = addedPoints_[f[fp]][i+offset];
2233 }
2234 }
2235 //Pout<< " newFace:" << newFace << endl;
2236 //Pout<< " coords:"
2237 // << UIndirectList<point>(meshMod.points(), newFace)
2238 // << " normal:" << newFace.unitNormal(meshMod.points())
2239 // << endl;
2240
2241 // Get new neighbour
2242 label own = addedCells[patchFacei][i];
2243 label nei;
2244 label patchi;
2245 label zoneI = -1;
2246 bool flip = false;
2247 bool fluxFlip = false;
2248
2249 if (i == addedCells[patchFacei].size()-1)
2250 {
2251 // Top layer so is either patch face or connects to
2252 // the other cell
2253 patchi = patchID[patchFacei];
2254 if (patchi == -1)
2255 {
2256 // Internal face
2257 nei =
2258 (
2259 !ppFlip[patchFacei]
2260 ? mesh_.faceNeighbour()[meshFacei]
2261 : mesh_.faceOwner()[meshFacei]
2262 );
2263
2264 if (ppFlip[patchFacei])
2265 {
2266 newFace = newFace.reverseFace();
2267 }
2268 //Pout<< "** adding top (internal) face:"
2269 // << " at:" << mesh_.faceCentres()[meshFacei]
2270 // << " own:" << own << " nei:" << nei
2271 // << " patchi:" << patchi
2272 // << " newFace:" << newFace
2273 // << endl;
2274 }
2275 else
2276 {
2277 nei = -1;
2278 }
2279 zoneI = mesh_.faceZones().whichZone(meshFacei);
2280 if (zoneI != -1)
2281 {
2282 const faceZone& fz = mesh_.faceZones()[zoneI];
2283 flip = fz.flipMap()[fz.whichFace(meshFacei)];
2284 }
2285 }
2286 else
2287 {
2288 // Internal face between layer i and i+1
2289 nei = addedCells[patchFacei][i+1];
2290 patchi = -1;
2291 }
2292
2293 if (nei != -1 && nei < own)
2294 {
2295 // Wrongly oriented internal face
2296 newFace = newFace.reverseFace();
2297 std::swap(own, nei);
2298 flip = !flip;
2299 fluxFlip = true;
2300
2301 //Pout<< "Flipped newFace:"
2302 // << newFace.unitNormal(meshMod.points())
2303 // << " own:" << own
2304 // << " nei:" << nei
2305 // << endl;
2306 }
2307
2308
2309
2310 layerFaces_[patchFacei][i+1] = meshMod.setAction
2311 (
2313 (
2314 newFace, // face
2315 own, // owner
2316 nei, // neighbour
2317 -1, // master point
2318 -1, // master edge
2319 (addToMesh_ ? meshFacei : -1), // master face
2320 fluxFlip, // flux flip
2321 patchi, // patch for face
2322 zoneI, // zone for face
2323 flip // face zone flip
2324 )
2325 );
2326
2327 //Pout<< "added layer face:" << layerFaces_[patchFacei][i+1]
2328 // << " verts:" << newFace
2329 // << " newFc:" << newFace.centre(meshMod.points())
2330 // << " originalFc:" << mesh_.faceCentres()[meshFacei]
2331 // << nl
2332 // << " n:" << newFace.unitNormal(meshMod.points())
2333 // << " own:" << own << " nei:" << nei
2334 // << " patchi:" << patchi
2335 // << endl;
2336 }
2337 }
2338 }
2339
2340 //
2341 // Modify owner faces to have addedCells as neighbour
2342 //
2343
2344 if (addToMesh_)
2345 {
2346 forAll(pp, patchFacei)
2347 {
2348 if (addedCells[patchFacei].size())
2349 {
2350 label meshFacei = pp.addressing()[patchFacei];
2351
2352 layerFaces_[patchFacei][0] = meshFacei;
2353 const face& f = pp[patchFacei];
2354
2355 const label own =
2356 (
2357 !ppFlip[patchFacei]
2358 ? mesh_.faceOwner()[meshFacei]
2359 : mesh_.faceNeighbour()[meshFacei]
2360 );
2361 const label nei = addedCells[patchFacei][0];
2362
2363 meshMod.setAction
2364 (
2366 (
2367 (ppFlip[patchFacei] ? f.reverseFace() : f),// verts
2368 meshFacei, // label of face
2369 own, // owner
2370 nei, // neighbour
2371 ppFlip[patchFacei], // face flip
2372 -1, // patch for face
2373 true, //false, // remove from zone
2374 -1, //zoneI, // zone for face
2375 false // face flip in zone
2376 )
2377 );
2378
2379 //Pout<< "Modified bottom face " << meshFacei
2380 // << " at:" << mesh_.faceCentres()[meshFacei]
2381 // << " new own:" << own << " new nei:" << nei
2382 // << " verts:" << meshMod.faces()[meshFacei]
2383 // << " n:"
2384 // << meshMod.faces()[meshFacei].unitNormal(meshMod.points())
2385 // << endl;
2386 }
2387 }
2388 }
2389 else
2390 {
2391 // If creating new mesh: reverse original faces and put them
2392 // in the exposed patch ID.
2393 forAll(pp, patchFacei)
2394 {
2395 if (addedCells[patchFacei].size())
2396 {
2397 label meshFacei = pp.addressing()[patchFacei];
2398 label zoneI = mesh_.faceZones().whichZone(meshFacei);
2399 bool zoneFlip = false;
2400 if (zoneI != -1)
2401 {
2402 const faceZone& fz = mesh_.faceZones()[zoneI];
2403 zoneFlip = !fz.flipMap()[fz.whichFace(meshFacei)];
2404 }
2405
2406 // Reverse and renumber old patch face.
2407 face f(pp.localFaces()[patchFacei].reverseFace());
2408 forAll(f, fp)
2409 {
2410 f[fp] = copiedPatchPoints[f[fp]];
2411 }
2412
2413 layerFaces_[patchFacei][0] = meshMod.setAction
2414 (
2416 (
2417 f, // modified face
2418 addedCells[patchFacei][0], // owner
2419 -1, // neighbour
2420 -1, // masterPoint
2421 -1, // masterEdge
2422 -1, // masterFace
2423 true, // face flip
2424 exposedPatchID[patchFacei], // patch for face
2425 zoneI, // zone for face
2426 zoneFlip // face flip in zone
2427 )
2428 );
2429 }
2430 }
2431 }
2432
2433
2434
2435 //
2436 // Create 'side' faces, one per edge that is being extended.
2437 //
2438
2439 const labelListList& faceEdges = pp.faceEdges();
2440 const faceList& localFaces = pp.localFaces();
2441 const edgeList& edges = pp.edges();
2442
2443 // Get number of layers per edge. This is 0 if edge is not extruded;
2444 // max of connected faces otherwise.
2445 labelList edgeLayers(pp.nEdges());
2446
2447 {
2448 // Use list over mesh.nEdges() since syncTools does not yet support
2449 // partial list synchronisation.
2450 labelList meshEdgeLayers(mesh_.nEdges(), -1);
2451
2452 forAll(meshEdges, edgei)
2453 {
2454 const edge& e = edges[edgei];
2455
2456 label meshEdgei = meshEdges[edgei];
2457
2458 if ((nPointLayers[e[0]] == 0) && (nPointLayers[e[1]] == 0))
2459 {
2460 meshEdgeLayers[meshEdgei] = 0;
2461 }
2462 else
2463 {
2464 const labelList& eFaces = pp.edgeFaces()[edgei];
2465
2466 forAll(eFaces, i)
2467 {
2468 meshEdgeLayers[meshEdgei] = max
2469 (
2470 nFaceLayers[eFaces[i]],
2471 meshEdgeLayers[meshEdgei]
2472 );
2473 }
2474 }
2475 }
2476
2478 (
2479 mesh_,
2480 meshEdgeLayers,
2482 label(0) // initial value
2483 );
2484
2485 forAll(meshEdges, edgei)
2486 {
2487 edgeLayers[edgei] = meshEdgeLayers[meshEdges[edgei]];
2488 }
2489 }
2490
2491
2492 // Mark off which edges have been extruded
2493 boolList doneEdge(pp.nEdges(), false);
2494
2495
2496 // Create faces. Per face walk connected edges and find string of edges
2497 // between the same two faces and extrude string into a single face.
2498 forAll(pp, patchFacei)
2499 {
2500 const labelList& fEdges = faceEdges[patchFacei];
2501
2502 forAll(fEdges, fp)
2503 {
2504 // Get string of edges that needs to be extruded as a single face.
2505 // Returned as indices in fEdges.
2506 labelPair indexPair
2507 (
2508 getEdgeString
2509 (
2510 pp,
2511 globalEdgeFaces,
2512 doneEdge,
2513 patchFacei,
2514 globalFaces.toGlobal(pp.addressing()[patchFacei])
2515 )
2516 );
2517
2518 //Pout<< "Found unextruded edges in edges:" << fEdges
2519 // << " start:" << indexPair[0]
2520 // << " end:" << indexPair[1]
2521 // << endl;
2522
2523 const label startFp = indexPair[0];
2524 const label endFp = indexPair[1];
2525
2526 if (startFp != -1)
2527 {
2528 // Extrude edges from indexPair[0] up to indexPair[1]
2529 // (note indexPair = indices of edges. There is one more vertex
2530 // than edges)
2531 const face& f = localFaces[patchFacei];
2532
2533 labelList stringedVerts;
2534 if (endFp >= startFp)
2535 {
2536 stringedVerts.setSize(endFp-startFp+2);
2537 }
2538 else
2539 {
2540 stringedVerts.setSize(endFp+f.size()-startFp+2);
2541 }
2542
2543 label fp = startFp;
2544
2545 for (label i = 0; i < stringedVerts.size()-1; i++)
2546 {
2547 stringedVerts[i] = f[fp];
2548 doneEdge[fEdges[fp]] = true;
2549 fp = f.fcIndex(fp);
2550 }
2551 stringedVerts.last() = f[fp];
2552
2553
2554 // Now stringedVerts contains the vertices in order of face f.
2555 // This is consistent with the order if f becomes the owner cell
2556 // and nbrFacei the neighbour cell. Note that the cells get
2557 // added in order of pp so we can just use face ordering and
2558 // because we loop in incrementing order as well we will
2559 // always have nbrFacei > patchFacei.
2560
2561 label startEdgei = fEdges[startFp];
2562
2563 label meshEdgei = meshEdges[startEdgei];
2564
2565 label numEdgeSideFaces = edgeLayers[startEdgei];
2566
2567 for (label i = 0; i < numEdgeSideFaces; i++)
2568 {
2569 label vEnd = stringedVerts.last();
2570 label vStart = stringedVerts[0];
2571
2572 // calculate number of points making up a face
2573 label newFp = 2*stringedVerts.size();
2574
2575 if (i == 0)
2576 {
2577 // layer 0 gets all the truncation of neighbouring
2578 // faces with more layers.
2579 if (addedPoints_[vEnd].size())
2580 {
2581 newFp +=
2582 addedPoints_[vEnd].size() - numEdgeSideFaces;
2583 }
2584 if (addedPoints_[vStart].size())
2585 {
2586 newFp +=
2587 addedPoints_[vStart].size() - numEdgeSideFaces;
2588 }
2589 }
2590
2591 face newFace(newFp);
2592
2593 newFp = 0;
2594
2595 // For layer 0 get pp points, for all other layers get
2596 // points of layer-1.
2597 if (i == 0)
2598 {
2599 forAll(stringedVerts, stringedI)
2600 {
2601 label v = stringedVerts[stringedI];
2602 addVertex
2603 (
2604 (
2605 addToMesh_
2606 ? meshPoints[v]
2607 : copiedPatchPoints[v]
2608 ),
2609 newFace,
2610 newFp
2611 );
2612 }
2613 }
2614 else
2615 {
2616 forAll(stringedVerts, stringedI)
2617 {
2618 label v = stringedVerts[stringedI];
2619 if (addedPoints_[v].size())
2620 {
2621 label offset =
2622 addedPoints_[v].size() - numEdgeSideFaces;
2623 addVertex
2624 (
2625 addedPoints_[v][i+offset-1],
2626 newFace,
2627 newFp
2628 );
2629 }
2630 else
2631 {
2632 addVertex
2633 (
2634 (
2635 addToMesh_
2636 ? meshPoints[v]
2637 : copiedPatchPoints[v]
2638 ),
2639 newFace,
2640 newFp
2641 );
2642 }
2643 }
2644 }
2645
2646 // add points between stringed vertices (end)
2647 if (numEdgeSideFaces < addedPoints_[vEnd].size())
2648 {
2649 if (i == 0 && addedPoints_[vEnd].size())
2650 {
2651 label offset =
2652 addedPoints_[vEnd].size() - numEdgeSideFaces;
2653 for (label ioff = 0; ioff < offset; ioff++)
2654 {
2655 addVertex
2656 (
2657 addedPoints_[vEnd][ioff],
2658 newFace,
2659 newFp
2660 );
2661 }
2662 }
2663 }
2664
2665 forAllReverse(stringedVerts, stringedI)
2666 {
2667 label v = stringedVerts[stringedI];
2668 if (addedPoints_[v].size())
2669 {
2670 label offset =
2671 addedPoints_[v].size() - numEdgeSideFaces;
2672 addVertex
2673 (
2674 addedPoints_[v][i+offset],
2675 newFace,
2676 newFp
2677 );
2678 }
2679 else
2680 {
2681 addVertex
2682 (
2683 (
2684 addToMesh_
2685 ? meshPoints[v]
2686 : copiedPatchPoints[v]
2687 ),
2688 newFace,
2689 newFp
2690 );
2691 }
2692 }
2693
2694
2695 // add points between stringed vertices (start)
2696 if (numEdgeSideFaces < addedPoints_[vStart].size())
2697 {
2698 if (i == 0 && addedPoints_[vStart].size())
2699 {
2700 label offset =
2701 addedPoints_[vStart].size() - numEdgeSideFaces;
2702 for (label ioff = offset-1; ioff >= 0; ioff--)
2703 {
2704 addVertex
2705 (
2706 addedPoints_[vStart][ioff],
2707 newFace,
2708 newFp
2709 );
2710 }
2711 }
2712 }
2713
2714 if (newFp >= 3)
2715 {
2716 // Add face inbetween faces patchFacei and nbrFacei
2717 // (possibly -1 for external edges)
2718
2719 newFace.setSize(newFp);
2720
2721 // Walked edges as if owner face was extruded. Reverse
2722 // for neighbour face extrusion.
2723 if (ppFlip[patchFacei])
2724 {
2725 newFace = newFace.reverseFace();
2726 }
2727
2728 if (debug)
2729 {
2730 labelHashSet verts(2*newFace.size());
2731 forAll(newFace, fp)
2732 {
2733 if (!verts.insert(newFace[fp]))
2734 {
2736 << "Duplicate vertex in face"
2737 << " to be added." << nl
2738 << "newFace:" << newFace << nl
2739 << "points:"
2741 (
2742 meshMod.points(),
2743 newFace
2744 ) << nl
2745 << "Layer:" << i
2746 << " out of:" << numEdgeSideFaces << nl
2747 << "ExtrudeEdge:" << meshEdgei
2748 << " at:"
2749 << mesh_.edges()[meshEdgei].line
2750 (
2751 mesh_.points()
2752 ) << nl
2753 << "string:" << stringedVerts
2754 << "stringpoints:"
2756 (
2757 pp.localPoints(),
2758 stringedVerts
2759 ) << nl
2760 << "stringNLayers:"
2761 << labelUIndList
2762 (
2763 nPointLayers,
2764 stringedVerts
2765 ) << nl
2766 << abort(FatalError);
2767 }
2768 }
2769 }
2770
2771 label nbrFacei = nbrFace
2772 (
2773 pp.edgeFaces(),
2774 startEdgei,
2775 patchFacei
2776 );
2777
2778 const labelList& meshFaces = mesh_.edgeFaces
2779 (
2780 meshEdgei,
2781 ef
2782 );
2783
2784 // Because we walk in order of patch face and in order
2785 // of face edges so face orientation will be opposite
2786 // that of the patch edge
2787 bool zoneFlip = false;
2788 if (edgeZoneID[startEdgei] != -1)
2789 {
2790 zoneFlip = !edgeFlip[startEdgei];
2791 }
2792
2793 addSideFace
2794 (
2795 pp,
2796 addedCells,
2797
2798 newFace, // vertices of new face
2799 edgePatchID[startEdgei],// -1 or patch for face
2800 edgeZoneID[startEdgei],
2801 zoneFlip,
2802 inflateFaceID[startEdgei],
2803
2804 patchFacei,
2805 nbrFacei,
2806 meshEdgei, // (mesh) edge to inflate
2807 i, // layer
2808 numEdgeSideFaces, // num layers
2809 meshFaces, // edgeFaces
2810 meshMod
2811 );
2812 }
2813 }
2814 }
2815 }
2816 }
2817
2818
2819 // Adjust side faces if they're on the side where points were duplicated
2820 // (i.e. adding to internal faces). Like duplicatePoints::setRefinement.
2821 if (addToMesh_)
2822 {
2823 face newFace;
2824
2825 forAll(baseFaces, facei)
2826 {
2827 const face& f = mesh_.faces()[facei];
2828 const face& baseF = baseFaces[facei];
2829
2830 if (isBlockedFace(facei) || baseF.empty())
2831 {
2832 // Either part of patch or no duplicated points on face
2833 continue;
2834 }
2835
2836 // Start off from original face
2837 newFace = f;
2838 forAll(f, fp)
2839 {
2840 const label meshPointi = f[fp];
2841 if (baseF[fp] != labelMax)
2842 {
2843 // Duplicated point
2844 const label patchPointi = pp.meshPointMap()[meshPointi];
2845 const label addedPointi = addedPoints_[patchPointi].last();
2846
2847 //Pout<< " For point:" << meshPointi
2848 // << " at:" << mesh_.points()[meshPointi]
2849 // << " at:" << pp.localPoints()[patchPointi]
2850 // << " using addedpoint:" << addedPointi
2851 // << " at:" << meshMod.points()[addedPointi]
2852 // << endl;
2853
2854 newFace[fp] = addedPointi;
2855 }
2856 }
2857
2858 //Pout<< "for face:" << facei << nl
2859 // << " old:" << f
2860 // << " n:" << newFace.unitNormal(meshMod.points())
2861 // << " coords:" << UIndirectList<point>(mesh_.points(), f)
2862 // << nl
2863 // << " new:" << newFace
2864 // << " n:" << newFace.unitNormal(meshMod.points())
2865 // << " coords:"
2866 // << UIndirectList<point>(meshMod.points(), newFace)
2867 // << endl;
2868
2869 // Get current zone info
2870 label zoneID = mesh_.faceZones().whichZone(facei);
2871 bool zoneFlip = false;
2872 if (zoneID >= 0)
2873 {
2874 const faceZone& fZone = mesh_.faceZones()[zoneID];
2875 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
2876 }
2877
2878
2879 if (mesh_.isInternalFace(facei))
2880 {
2881 meshMod.modifyFace
2882 (
2883 newFace, // modified face
2884 facei, // label of modified face
2885 mesh_.faceOwner()[facei], // owner
2886 mesh_.faceNeighbour()[facei], // neighbour
2887 false, // face flip
2888 -1, // patch for face
2889 zoneID, // zone for face
2890 zoneFlip // face flip in zone
2891 );
2892 }
2893 else
2894 {
2895 meshMod.modifyFace
2896 (
2897 newFace, // modified face
2898 facei, // label of modified face
2899 mesh_.faceOwner()[facei], // owner
2900 -1, // neighbour
2901 false, // face flip
2902 patches.whichPatch(facei), // patch for face
2903 zoneID, // zone for face
2904 zoneFlip // face flip in zone
2905 );
2906 }
2907 }
2908 }
2909}
2910
2911
2913(
2914 const mapPolyMesh& morphMap,
2915 const labelList& faceMap, // new to old patch faces
2916 const labelList& pointMap // new to old patch points
2917)
2918{
2919 {
2920 labelListList newAddedPoints(pointMap.size());
2921
2922 forAll(newAddedPoints, newPointi)
2923 {
2924 label oldPointi = pointMap[newPointi];
2925
2926 const labelList& added = addedPoints_[oldPointi];
2927
2928 labelList& newAdded = newAddedPoints[newPointi];
2929 newAdded.setSize(added.size());
2930 label newI = 0;
2931
2932 forAll(added, i)
2933 {
2934 label newPointi = morphMap.reversePointMap()[added[i]];
2935
2936 if (newPointi >= 0)
2937 {
2938 newAdded[newI++] = newPointi;
2939 }
2940 }
2941 newAdded.setSize(newI);
2942 }
2943 addedPoints_.transfer(newAddedPoints);
2944 }
2945
2946 {
2947 labelListList newLayerFaces(faceMap.size());
2948
2949 forAll(newLayerFaces, newFacei)
2950 {
2951 label oldFacei = faceMap[newFacei];
2952
2953 const labelList& added = layerFaces_[oldFacei];
2954
2955 labelList& newAdded = newLayerFaces[newFacei];
2956 newAdded.setSize(added.size());
2957 label newI = 0;
2958
2959 forAll(added, i)
2960 {
2961 label newFacei = morphMap.reverseFaceMap()[added[i]];
2962
2963 if (newFacei >= 0)
2964 {
2965 newAdded[newI++] = newFacei;
2966 }
2967 }
2968 newAdded.setSize(newI);
2969 }
2970 layerFaces_.transfer(newLayerFaces);
2971 }
2972}
2973
2974
2975// ************************************************************************* //
scalar y
label k
bool found
label n
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:191
const T & lookup(const Key &key, const T &deflt) const
Return hashed entry if it exists, or return the given default.
Definition: HashTableI.H:224
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
Definition: HashTableI.H:180
iterator find(const Key &key)
Find and return an iterator set at the hashed entry.
Definition: HashTableI.H:114
void clear()
Clear all entries from table.
Definition: HashTable.C:678
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
void setSize(const label n, unsigned int val=0u)
Alias for resize()
Definition: PackedList.H:557
label size() const noexcept
Number of entries.
Definition: PackedListI.H:377
static void matchEdges(const PrimitivePatch< FaceList1, PointField1 > &p1, const PrimitivePatch< FaceList2, PointField2 > &p2, labelList &p1EdgeLabels, labelList &p2EdgeLabels, bitSet &sameOrientation)
Find corresponding edges on patches sharing the same points.
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.
label nInternalEdges() const
Number of internal edges.
labelList meshEdges(const edgeList &allEdges, const labelListList &cellEdges, const labelList &faceCells) const
const Map< label > & meshPointMap() const
Mesh point map.
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 > & points() const noexcept
Return reference to global points.
const labelListList & edgeFaces() const
Return edge-face addressing.
const labelListList & faceEdges() const
Return face-edge addressing.
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
A List obtained as a section of another List.
Definition: SubList.H:70
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
Definition: UListI.H:427
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 fcIndex(const label i) const noexcept
Definition: UListI.H:60
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: ZoneMesh.C:288
Adds layers of cells to outside of polyPatch. Can optionally create stand-alone extruded mesh (addToM...
void setRefinement(const globalIndex &globalFaces, const labelListList &globalEdgeFaces, const scalarField &expansionRatio, const indirectPrimitivePatch &pp, const bitSet &flip, const labelList &sidePatchID, const labelList &sideZoneID, const boolList &sideFlip, const labelList &inflateFaceID, const labelList &exposedPatchID, const labelList &nFaceLayers, const labelList &nPointLayers, const vectorField &firstLayerDisp, polyTopoChange &meshMod)
Play commands into polyTopoChange to create layers on top.
static void calcExtrudeInfo(const bool zoneFromAnyFace, const polyMesh &, const globalIndex &globalFaces, const labelListList &globalEdgeFaces, const indirectPrimitivePatch &pp, labelList &edgePatchID, label &nPatches, Map< label > &nbrProcToPatch, Map< label > &patchToNbrProc, labelList &edgeZoneID, boolList &edgeFlip, labelList &inflateFaceID)
Determine extrude information per patch edge:
labelListList addedCells() const
Added cells given current mesh & layerfaces.
static labelListList globalEdgeFaces(const polyMesh &, const globalIndex &globalFaces, const indirectPrimitivePatch &pp)
Per patch edge the pp faces (in global indices) using it.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
unsigned int count(const bool on=true) const
Count number of bits set.
Definition: bitSetI.H:500
void set(const bitSet &bitset)
Set specified bits from another bitset.
Definition: bitSetI.H:590
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
static int compare(const edge &a, const edge &b)
Compare edges.
Definition: edgeI.H:33
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:67
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:372
const boolList & flipMap() const noexcept
Return face flip map.
Definition: faceZone.H:272
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
face reverseFace() const
Return face with reverse direction.
Definition: face.C:636
Determination and storage of the possible independent transforms introduced by coupledPolyPatches,...
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
label toLocal(const label i) const
From global to local on current processor.
Definition: globalIndexI.H:325
bool isLocal(const label i) const
Is on local processor.
Definition: globalIndexI.H:244
Various mesh related information for a parallel run. Upon construction, constructs all info using par...
const labelListList & globalEdgeTransformedSlaves() const
const labelList & processorPatches() const noexcept
Return list of processor patch labels.
const mapDistribute & globalEdgeSlavesMap() const
static void syncData(List< Type > &elems, const labelListList &slaves, const labelListList &transformedSlaves, const mapDistribute &slavesMap, const globalIndexAndTransform &, const CombineOp &cop, const TransformOp &top)
Helper: synchronise data with transforms.
const globalIndexAndTransform & globalTransforms() const
Global transforms numbering.
const indirectPrimitivePatch & coupledPatch() const
Return patch of all coupled faces.
const labelListList & globalEdgeSlaves() const
const processorTopology & topology() const noexcept
The processor to processor topology.
Class containing processor-to-processor mapping information.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:162
const labelList & reversePointMap() const
Reverse point map.
Definition: mapPolyMesh.H:469
const labelList & reverseFaceMap() const
Reverse face map.
Definition: mapPolyMesh.H:501
void operator()(face &x, const face &y) const
void updateMesh()
Update for new mesh topology.
static const complex max
complex (VGREAT,VGREAT)
Definition: complex.H:293
Class containing data for cell addition.
Definition: polyAddCell.H:51
A face addition data class. A face can be inflated either from a point or from another face and can e...
Definition: polyAddFace.H:56
Class containing data for point addition.
Definition: polyAddPoint.H:52
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
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
Definition: polyMesh.H:498
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1108
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1310
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1127
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1083
Class describing modification of a face.
Direct mesh changes based on v1.3 polyTopoChange syntax.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
void modifyFace(const face &f, const label facei, const label own, const label nei, const bool flipFaceFlux, const label patchID, const label zoneID, const bool zoneFlip)
Modify vertices or cell of face.
const DynamicList< point > & points() const
Points. Shrunk after constructing mesh (or calling of compact())
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
const labelListList & pointEdges() const
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
label nFaces() const noexcept
Number of mesh faces.
const labelListList & edgeFaces() const
label nEdges() const
Number of mesh edges.
label procPatchLookup(const label proci) const
Which local boundary is attached to specified processor.
static bitSet getInternalOrCoupledFaces(const polyMesh &mesh)
Get per face whether it is internal or coupled.
Definition: syncTools.C:176
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &faceValues)
Swap coupled boundary face values. Uses eqOp.
Definition: syncTools.H:445
static void syncPointList(const polyMesh &mesh, List< T > &pointValues, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
static void syncEdgeList(const polyMesh &mesh, List< T > &edgeValues, const CombineOp &cop, const T &nullValue, const TransformOp &top, const FlipOp &fop)
Synchronize values on all mesh edges.
static void syncBoundaryFaceList(const polyMesh &mesh, UList< T > &faceValues, const CombineOp &cop, const TransformOp &top, const bool parRun=Pstream::parRun())
Synchronize values on boundary faces only.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
const polyBoundaryMesh & patches
dynamicFvMesh & mesh
Dummy transform to be used with syncTools.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const label nPatches
const labelIOList & zoneID
#define WarningInFunction
Report a warning using Foam::Warning.
Namespace for OpenFOAM.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
Pair< label > labelPair
A pair of labels.
Definition: Pair.H:57
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
List< label > labelList
A List of labels.
Definition: List.H:66
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
constexpr label labelMax
Definition: label.H:61
ZoneMesh< faceZone, polyMesh > faceZoneMesh
A ZoneMesh with the type faceZone.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
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
List< bool > boolList
A List of bools.
Definition: List.H:64
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
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.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Type gMax(const FieldField< Field, Type > &f)
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition: IndirectList.H:68
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
label newPointi
Definition: readKivaGrid.H:496
labelList f(nPoints)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: stdFoam.H:346
List helper to append y unique elements onto the end of x.
Definition: ListOps.H:570