hexRef8.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2016-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 "hexRef8.H"
30
31#include "polyMesh.H"
32#include "polyTopoChange.H"
33#include "meshTools.H"
34#include "polyAddFace.H"
35#include "polyAddPoint.H"
36#include "polyAddCell.H"
37#include "polyModifyFace.H"
38#include "syncTools.H"
39#include "faceSet.H"
40#include "cellSet.H"
41#include "pointSet.H"
42#include "labelPairHashes.H"
43#include "OFstream.H"
44#include "Time.H"
45#include "FaceCellWave.H"
47#include "refinementData.H"
49#include "degenerateMatcher.H"
50
51// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
52
53namespace Foam
54{
56
57 //- Reduction class. If x and y are not equal assign value.
58 template<label value>
59 struct ifEqEqOp
60 {
61 void operator()(label& x, const label y) const
62 {
63 x = (x == y) ? x : value;
64 }
65 };
66}
67
68
69// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
70
71void Foam::hexRef8::reorder
72(
73 const labelList& map,
74 const label len,
75 const label null,
76 labelList& elems
77)
78{
79 labelList newElems(len, null);
80
81 forAll(elems, i)
82 {
83 label newI = map[i];
84
85 if (newI >= len)
86 {
88 }
89
90 if (newI >= 0)
91 {
92 newElems[newI] = elems[i];
93 }
94 }
95
96 elems.transfer(newElems);
97}
98
99
100void Foam::hexRef8::getFaceInfo
101(
102 const label facei,
103 label& patchID,
104 label& zoneID,
105 label& zoneFlip
106) const
107{
108 patchID = -1;
109
110 if (!mesh_.isInternalFace(facei))
111 {
112 patchID = mesh_.boundaryMesh().whichPatch(facei);
113 }
114
115 zoneID = mesh_.faceZones().whichZone(facei);
116
117 zoneFlip = false;
118
119 if (zoneID >= 0)
120 {
121 const faceZone& fZone = mesh_.faceZones()[zoneID];
122
123 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
124 }
125}
126
127
128// Adds a face on top of existing facei.
129Foam::label Foam::hexRef8::addFace
130(
131 polyTopoChange& meshMod,
132 const label facei,
133 const face& newFace,
134 const label own,
135 const label nei
136) const
137{
138 label patchID, zoneID, zoneFlip;
139
140 getFaceInfo(facei, patchID, zoneID, zoneFlip);
141
142 label newFacei = -1;
143
144 if ((nei == -1) || (own < nei))
145 {
146 // Ordering ok.
147 newFacei = meshMod.setAction
148 (
149 polyAddFace
150 (
151 newFace, // face
152 own, // owner
153 nei, // neighbour
154 -1, // master point
155 -1, // master edge
156 facei, // master face for addition
157 false, // flux flip
158 patchID, // patch for face
159 zoneID, // zone for face
160 zoneFlip // face zone flip
161 )
162 );
163 }
164 else
165 {
166 // Reverse owner/neighbour
167 newFacei = meshMod.setAction
168 (
169 polyAddFace
170 (
171 newFace.reverseFace(), // face
172 nei, // owner
173 own, // neighbour
174 -1, // master point
175 -1, // master edge
176 facei, // master face for addition
177 false, // flux flip
178 patchID, // patch for face
179 zoneID, // zone for face
180 zoneFlip // face zone flip
181 )
182 );
183 }
184 return newFacei;
185}
186
187
188// Adds an internal face from an edge. Assumes orientation correct.
189// Problem is that the face is between four new vertices. So what do we provide
190// as master? The only existing mesh item we have is the edge we have split.
191// Have to be careful in only using it if it has internal faces since otherwise
192// polyMeshMorph will complain (because it cannot generate a sensible mapping
193// for the face)
194Foam::label Foam::hexRef8::addInternalFace
195(
196 polyTopoChange& meshMod,
197 const label meshFacei,
198 const label meshPointi,
199 const face& newFace,
200 const label own,
201 const label nei
202) const
203{
204 if (mesh_.isInternalFace(meshFacei))
205 {
206 return meshMod.setAction
207 (
208 polyAddFace
209 (
210 newFace, // face
211 own, // owner
212 nei, // neighbour
213 -1, // master point
214 -1, // master edge
215 -1, // master face for addition
216 false, // flux flip
217 -1, // patch for face
218 -1, // zone for face
219 false // face zone flip
220 )
221 );
222 }
223 else
224 {
225 // Two choices:
226 // - append (i.e. create out of nothing - will not be mapped)
227 // problem: field does not get mapped.
228 // - inflate from point.
229 // problem: does interpolative mapping which constructs full
230 // volPointInterpolation!
231
232 // For now create out of nothing
233
234 return meshMod.setAction
235 (
236 polyAddFace
237 (
238 newFace, // face
239 own, // owner
240 nei, // neighbour
241 -1, // master point
242 -1, // master edge
243 -1, // master face for addition
244 false, // flux flip
245 -1, // patch for face
246 -1, // zone for face
247 false // face zone flip
248 )
249 );
250
251
254 //label masterPointi = -1;
255 //
256 //const labelList& pFaces = mesh_.pointFaces()[meshPointi];
257 //
258 //forAll(pFaces, i)
259 //{
260 // if (mesh_.isInternalFace(pFaces[i]))
261 // {
262 // // meshPoint uses internal faces so ok to inflate from it
263 // masterPointi = meshPointi;
264 //
265 // break;
266 // }
267 //}
268 //
269 //return meshMod.setAction
270 //(
271 // polyAddFace
272 // (
273 // newFace, // face
274 // own, // owner
275 // nei, // neighbour
276 // masterPointi, // master point
277 // -1, // master edge
278 // -1, // master face for addition
279 // false, // flux flip
280 // -1, // patch for face
281 // -1, // zone for face
282 // false // face zone flip
283 // )
284 //);
285 }
286}
287
288
289// Modifies existing facei for either new owner/neighbour or new face points.
290void Foam::hexRef8::modFace
291(
292 polyTopoChange& meshMod,
293 const label facei,
294 const face& newFace,
295 const label own,
296 const label nei
297) const
298{
299 label patchID, zoneID, zoneFlip;
300
301 getFaceInfo(facei, patchID, zoneID, zoneFlip);
302
303 if
304 (
305 (own != mesh_.faceOwner()[facei])
306 || (
307 mesh_.isInternalFace(facei)
308 && (nei != mesh_.faceNeighbour()[facei])
309 )
310 || (newFace != mesh_.faces()[facei])
311 )
312 {
313 if ((nei == -1) || (own < nei))
314 {
315 meshMod.setAction
316 (
317 polyModifyFace
318 (
319 newFace, // modified face
320 facei, // label of face being modified
321 own, // owner
322 nei, // neighbour
323 false, // face flip
324 patchID, // patch for face
325 false, // remove from zone
326 zoneID, // zone for face
327 zoneFlip // face flip in zone
328 )
329 );
330 }
331 else
332 {
333 meshMod.setAction
334 (
335 polyModifyFace
336 (
337 newFace.reverseFace(), // modified face
338 facei, // label of face being modified
339 nei, // owner
340 own, // neighbour
341 false, // face flip
342 patchID, // patch for face
343 false, // remove from zone
344 zoneID, // zone for face
345 zoneFlip // face flip in zone
346 )
347 );
348 }
349 }
350}
351
352
353// Bit complex way to determine the unrefined edge length.
354Foam::scalar Foam::hexRef8::getLevel0EdgeLength() const
355{
356 if (cellLevel_.size() != mesh_.nCells())
357 {
359 << "Number of cells in mesh:" << mesh_.nCells()
360 << " does not equal size of cellLevel:" << cellLevel_.size()
361 << endl
362 << "This might be because of a restart with inconsistent cellLevel."
363 << abort(FatalError);
364 }
365
366 // Determine minimum edge length per refinement level
367 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
368
369 const scalar GREAT2 = sqr(GREAT);
370
371 label nLevels = gMax(cellLevel_)+1;
372
373 scalarField typEdgeLenSqr(nLevels, GREAT2);
374
375
376 // 1. Look only at edges surrounded by cellLevel cells only.
377 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
378
379 {
380 // Per edge the cellLevel of connected cells. -1 if not set,
381 // labelMax if different levels, otherwise levels of connected cells.
382 labelList edgeLevel(mesh_.nEdges(), -1);
383
384 forAll(cellLevel_, celli)
385 {
386 const label cLevel = cellLevel_[celli];
387
388 const labelList& cEdges = mesh_.cellEdges(celli);
389
390 forAll(cEdges, i)
391 {
392 label edgeI = cEdges[i];
393
394 if (edgeLevel[edgeI] == -1)
395 {
396 edgeLevel[edgeI] = cLevel;
397 }
398 else if (edgeLevel[edgeI] == labelMax)
399 {
400 // Already marked as on different cellLevels
401 }
402 else if (edgeLevel[edgeI] != cLevel)
403 {
404 edgeLevel[edgeI] = labelMax;
405 }
406 }
407 }
408
409 // Make sure that edges with different levels on different processors
410 // are also marked. Do the same test (edgeLevel != cLevel) on coupled
411 // edges.
413 (
414 mesh_,
415 edgeLevel,
416 ifEqEqOp<labelMax>(),
418 );
419
420 // Now use the edgeLevel with a valid value to determine the
421 // length per level.
422 forAll(edgeLevel, edgeI)
423 {
424 const label eLevel = edgeLevel[edgeI];
425
426 if (eLevel >= 0 && eLevel < labelMax)
427 {
428 const edge& e = mesh_.edges()[edgeI];
429
430 scalar edgeLenSqr = magSqr(e.vec(mesh_.points()));
431
432 typEdgeLenSqr[eLevel] = min(typEdgeLenSqr[eLevel], edgeLenSqr);
433 }
434 }
435 }
436
437 // Get the minimum per level over all processors. Note minimum so if
438 // cells are not cubic we use the smallest edge side.
439 Pstream::listCombineAllGather(typEdgeLenSqr, minEqOp<scalar>());
440
441 if (debug)
442 {
443 Pout<< "hexRef8::getLevel0EdgeLength() :"
444 << " After phase1: Edgelengths (squared) per refinementlevel:"
445 << typEdgeLenSqr << endl;
446 }
447
448
449 // 2. For any levels where we haven't determined a valid length yet
450 // use any surrounding cell level. Here we use the max so we don't
451 // pick up levels between celllevel and higher celllevel (will have
452 // edges sized according to highest celllevel)
453 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
454
455 scalarField maxEdgeLenSqr(nLevels, -GREAT2);
456
457 forAll(cellLevel_, celli)
458 {
459 const label cLevel = cellLevel_[celli];
460
461 const labelList& cEdges = mesh_.cellEdges(celli);
462
463 forAll(cEdges, i)
464 {
465 const edge& e = mesh_.edges()[cEdges[i]];
466
467 scalar edgeLenSqr = magSqr(e.vec(mesh_.points()));
468
469 maxEdgeLenSqr[cLevel] = max(maxEdgeLenSqr[cLevel], edgeLenSqr);
470 }
471 }
472
473 Pstream::listCombineAllGather(maxEdgeLenSqr, maxEqOp<scalar>());
474
475 if (debug)
476 {
477 Pout<< "hexRef8::getLevel0EdgeLength() :"
478 << " Poor Edgelengths (squared) per refinementlevel:"
479 << maxEdgeLenSqr << endl;
480 }
481
482
483 // 3. Combine the two sets of lengths
484 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
485
486 forAll(typEdgeLenSqr, levelI)
487 {
488 if (typEdgeLenSqr[levelI] == GREAT2 && maxEdgeLenSqr[levelI] >= 0)
489 {
490 typEdgeLenSqr[levelI] = maxEdgeLenSqr[levelI];
491 }
492 }
493
494 if (debug)
495 {
496 Pout<< "hexRef8::getLevel0EdgeLength() :"
497 << " Final Edgelengths (squared) per refinementlevel:"
498 << typEdgeLenSqr << endl;
499 }
500
501 // Find lowest level present
502 scalar level0Size = -1;
503
504 forAll(typEdgeLenSqr, levelI)
505 {
506 scalar lenSqr = typEdgeLenSqr[levelI];
507
508 if (lenSqr < GREAT2)
509 {
510 level0Size = Foam::sqrt(lenSqr)*(1<<levelI);
511
512 if (debug)
513 {
514 Pout<< "hexRef8::getLevel0EdgeLength() :"
515 << " For level:" << levelI
516 << " have edgeLen:" << Foam::sqrt(lenSqr)
517 << " with equivalent level0 len:" << level0Size
518 << endl;
519 }
520 break;
521 }
522 }
523
524 if (level0Size == -1)
525 {
527 << "Problem : typEdgeLenSqr:" << typEdgeLenSqr << abort(FatalError);
528 }
529
530 return level0Size;
531}
532
533
534// Check whether pointi is an anchor on celli.
535// If it is not check whether any other point on the face is an anchor cell.
536Foam::label Foam::hexRef8::getAnchorCell
537(
538 const labelListList& cellAnchorPoints,
539 const labelListList& cellAddedCells,
540 const label celli,
541 const label facei,
542 const label pointi
543) const
544{
545 if (cellAnchorPoints[celli].size())
546 {
547 label index = cellAnchorPoints[celli].find(pointi);
548
549 if (index != -1)
550 {
551 return cellAddedCells[celli][index];
552 }
553
554
555 // pointi is not an anchor cell.
556 // Maybe we are already a refined face so check all the face
557 // vertices.
558 const face& f = mesh_.faces()[facei];
559
560 forAll(f, fp)
561 {
562 label index = cellAnchorPoints[celli].find(f[fp]);
563
564 if (index != -1)
565 {
566 return cellAddedCells[celli][index];
567 }
568 }
569
570 // Problem.
571 dumpCell(celli);
572 Perr<< "cell:" << celli << " anchorPoints:" << cellAnchorPoints[celli]
573 << endl;
574
576 << "Could not find point " << pointi
577 << " in the anchorPoints for cell " << celli << endl
578 << "Does your original mesh obey the 2:1 constraint and"
579 << " did you use consistentRefinement to make your cells to refine"
580 << " obey this constraint as well?"
581 << abort(FatalError);
582
583 return -1;
584 }
585 else
586 {
587 return celli;
588 }
589}
590
591
592// Get new owner and neighbour
593void Foam::hexRef8::getFaceNeighbours
594(
595 const labelListList& cellAnchorPoints,
596 const labelListList& cellAddedCells,
597 const label facei,
598 const label pointi,
599
600 label& own,
601 label& nei
602) const
603{
604 // Is owner split?
605 own = getAnchorCell
606 (
607 cellAnchorPoints,
608 cellAddedCells,
609 mesh_.faceOwner()[facei],
610 facei,
611 pointi
612 );
613
614 if (mesh_.isInternalFace(facei))
615 {
616 nei = getAnchorCell
617 (
618 cellAnchorPoints,
619 cellAddedCells,
620 mesh_.faceNeighbour()[facei],
621 facei,
622 pointi
623 );
624 }
625 else
626 {
627 nei = -1;
628 }
629}
630
631
632// Get point with the lowest pointLevel
633Foam::label Foam::hexRef8::findMinLevel(const labelList& f) const
634{
635 label minLevel = labelMax;
636 label minFp = -1;
637
638 forAll(f, fp)
639 {
640 label level = pointLevel_[f[fp]];
641
642 if (level < minLevel)
643 {
644 minLevel = level;
645 minFp = fp;
646 }
647 }
648
649 return minFp;
650}
651
652
653// Get point with the highest pointLevel
654Foam::label Foam::hexRef8::findMaxLevel(const labelList& f) const
655{
656 label maxLevel = labelMin;
657 label maxFp = -1;
658
659 forAll(f, fp)
660 {
661 label level = pointLevel_[f[fp]];
662
663 if (level > maxLevel)
664 {
665 maxLevel = level;
666 maxFp = fp;
667 }
668 }
669
670 return maxFp;
671}
672
673
674Foam::label Foam::hexRef8::countAnchors
675(
676 const labelList& f,
677 const label anchorLevel
678) const
679{
680 label nAnchors = 0;
681
682 forAll(f, fp)
683 {
684 if (pointLevel_[f[fp]] <= anchorLevel)
685 {
686 nAnchors++;
687 }
688 }
689 return nAnchors;
690}
691
692
693void Foam::hexRef8::dumpCell(const label celli) const
694{
695 OFstream str(mesh_.time().path()/"cell_" + Foam::name(celli) + ".obj");
696 Pout<< "hexRef8 : Dumping cell as obj to " << str.name() << endl;
697
698 const cell& cFaces = mesh_.cells()[celli];
699
700 Map<label> pointToObjVert;
701 label objVertI = 0;
702
703 forAll(cFaces, i)
704 {
705 const face& f = mesh_.faces()[cFaces[i]];
706
707 forAll(f, fp)
708 {
709 if (pointToObjVert.insert(f[fp], objVertI))
710 {
711 meshTools::writeOBJ(str, mesh_.points()[f[fp]]);
712 objVertI++;
713 }
714 }
715 }
716
717 forAll(cFaces, i)
718 {
719 const face& f = mesh_.faces()[cFaces[i]];
720
721 forAll(f, fp)
722 {
723 label pointi = f[fp];
724 label nexPointi = f[f.fcIndex(fp)];
725
726 str << "l " << pointToObjVert[pointi]+1
727 << ' ' << pointToObjVert[nexPointi]+1 << nl;
728 }
729 }
730}
731
732
733// Find point with certain pointLevel. Skip any higher levels.
734Foam::label Foam::hexRef8::findLevel
735(
736 const label facei,
737 const face& f,
738 const label startFp,
739 const bool searchForward,
740 const label wantedLevel
741) const
742{
743 label fp = startFp;
744
745 forAll(f, i)
746 {
747 label pointi = f[fp];
748
749 if (pointLevel_[pointi] < wantedLevel)
750 {
751 dumpCell(mesh_.faceOwner()[facei]);
752 if (mesh_.isInternalFace(facei))
753 {
754 dumpCell(mesh_.faceNeighbour()[facei]);
755 }
756
758 << "face:" << f
759 << " level:" << labelUIndList(pointLevel_, f)
760 << " startFp:" << startFp
761 << " wantedLevel:" << wantedLevel
762 << abort(FatalError);
763 }
764 else if (pointLevel_[pointi] == wantedLevel)
765 {
766 return fp;
767 }
768
769 if (searchForward)
770 {
771 fp = f.fcIndex(fp);
772 }
773 else
774 {
775 fp = f.rcIndex(fp);
776 }
777 }
778
779 dumpCell(mesh_.faceOwner()[facei]);
780 if (mesh_.isInternalFace(facei))
781 {
782 dumpCell(mesh_.faceNeighbour()[facei]);
783 }
784
786 << "face:" << f
787 << " level:" << labelUIndList(pointLevel_, f)
788 << " startFp:" << startFp
789 << " wantedLevel:" << wantedLevel
790 << abort(FatalError);
791
792 return -1;
793}
794
795
796// Gets cell level such that the face has four points <= level.
797Foam::label Foam::hexRef8::faceLevel(const label facei) const
798{
799 const face& f = mesh_.faces()[facei];
800
801 if (f.size() <= 4)
802 {
803 return pointLevel_[f[findMaxLevel(f)]];
804 }
805 else
806 {
807 label ownLevel = cellLevel_[mesh_.faceOwner()[facei]];
808
809 if (countAnchors(f, ownLevel) == 4)
810 {
811 return ownLevel;
812 }
813 else if (countAnchors(f, ownLevel+1) == 4)
814 {
815 return ownLevel+1;
816 }
817 else
818 {
819 return -1;
820 }
821 }
822}
823
824
825void Foam::hexRef8::checkInternalOrientation
826(
827 polyTopoChange& meshMod,
828 const label celli,
829 const label facei,
830 const point& ownPt,
831 const point& neiPt,
832 const face& newFace
833)
834{
835 face compactFace(identity(newFace.size()));
836 pointField compactPoints(meshMod.points(), newFace);
837
838 const vector areaNorm(compactFace.areaNormal(compactPoints));
839
840 const vector dir(neiPt - ownPt);
841
842 if ((dir & areaNorm) < 0)
843 {
845 << "cell:" << celli << " old face:" << facei
846 << " newFace:" << newFace << endl
847 << " coords:" << compactPoints
848 << " ownPt:" << ownPt
849 << " neiPt:" << neiPt
850 << abort(FatalError);
851 }
852
853 const vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
854
855 const scalar s = (fcToOwn & areaNorm) / (dir & areaNorm);
856
857 if (s < 0.1 || s > 0.9)
858 {
860 << "cell:" << celli << " old face:" << facei
861 << " newFace:" << newFace << endl
862 << " coords:" << compactPoints
863 << " ownPt:" << ownPt
864 << " neiPt:" << neiPt
865 << " s:" << s
866 << abort(FatalError);
867 }
868}
869
870
871void Foam::hexRef8::checkBoundaryOrientation
872(
873 polyTopoChange& meshMod,
874 const label celli,
875 const label facei,
876 const point& ownPt,
877 const point& boundaryPt,
878 const face& newFace
879)
880{
881 face compactFace(identity(newFace.size()));
882 pointField compactPoints(meshMod.points(), newFace);
883
884 const vector areaNorm(compactFace.areaNormal(compactPoints));
885
886 const vector dir(boundaryPt - ownPt);
887
888 if ((dir & areaNorm) < 0)
889 {
891 << "cell:" << celli << " old face:" << facei
892 << " newFace:" << newFace
893 << " coords:" << compactPoints
894 << " ownPt:" << ownPt
895 << " boundaryPt:" << boundaryPt
896 << abort(FatalError);
897 }
898
899 const vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
900
901 const scalar s = (fcToOwn & dir) / magSqr(dir);
902
903 if (s < 0.7 || s > 1.3)
904 {
906 << "cell:" << celli << " old face:" << facei
907 << " newFace:" << newFace
908 << " coords:" << compactPoints
909 << " ownPt:" << ownPt
910 << " boundaryPt:" << boundaryPt
911 << " s:" << s
912 << endl;
913 //<< abort(FatalError);
914 }
915}
916
917
918// If p0 and p1 are existing vertices check if edge is split and insert
919// splitPoint.
920void Foam::hexRef8::insertEdgeSplit
921(
922 const labelList& edgeMidPoint,
923 const label p0,
924 const label p1,
925 DynamicList<label>& verts
926) const
927{
928 if (p0 < mesh_.nPoints() && p1 < mesh_.nPoints())
929 {
930 label edgeI = meshTools::findEdge(mesh_, p0, p1);
931
932 if (edgeI != -1 && edgeMidPoint[edgeI] != -1)
933 {
934 verts.append(edgeMidPoint[edgeI]);
935 }
936 }
937}
938
939
940// Internal faces are one per edge between anchor points. So one per midPoint
941// between the anchor points. Here we store the information on the midPoint
942// and if we have enough information:
943// - two anchors
944// - two face mid points
945// we add the face. Note that this routine can get called anywhere from
946// two times (two unrefined faces) to four times (two refined faces) so
947// the first call that adds the information creates the face.
948Foam::label Foam::hexRef8::storeMidPointInfo
949(
950 const labelListList& cellAnchorPoints,
951 const labelListList& cellAddedCells,
952 const labelList& cellMidPoint,
953 const labelList& edgeMidPoint,
954 const label celli,
955 const label facei,
956 const bool faceOrder,
957 const label edgeMidPointi,
958 const label anchorPointi,
959 const label faceMidPointi,
960
961 Map<edge>& midPointToAnchors,
962 Map<edge>& midPointToFaceMids,
963 polyTopoChange& meshMod
964) const
965{
966 // See if need to store anchors.
967
968 bool changed = false;
969 bool haveTwoAnchors = false;
970
971 auto edgeMidFnd = midPointToAnchors.find(edgeMidPointi);
972
973 if (!edgeMidFnd.found())
974 {
975 midPointToAnchors.insert(edgeMidPointi, edge(anchorPointi, -1));
976 }
977 else
978 {
979 edge& e = edgeMidFnd.val();
980
981 if (anchorPointi != e[0])
982 {
983 if (e[1] == -1)
984 {
985 e[1] = anchorPointi;
986 changed = true;
987 }
988 }
989
990 if (e[0] != -1 && e[1] != -1)
991 {
992 haveTwoAnchors = true;
993 }
994 }
995
996 bool haveTwoFaceMids = false;
997
998 auto faceMidFnd = midPointToFaceMids.find(edgeMidPointi);
999
1000 if (!faceMidFnd.found())
1001 {
1002 midPointToFaceMids.insert(edgeMidPointi, edge(faceMidPointi, -1));
1003 }
1004 else
1005 {
1006 edge& e = faceMidFnd.val();
1007
1008 if (faceMidPointi != e[0])
1009 {
1010 if (e[1] == -1)
1011 {
1012 e[1] = faceMidPointi;
1013 changed = true;
1014 }
1015 }
1016
1017 if (e[0] != -1 && e[1] != -1)
1018 {
1019 haveTwoFaceMids = true;
1020 }
1021 }
1022
1023 // Check if this call of storeMidPointInfo is the one that completed all
1024 // the necessary information.
1025
1026 if (changed && haveTwoAnchors && haveTwoFaceMids)
1027 {
1028 const edge& anchors = midPointToAnchors[edgeMidPointi];
1029 const edge& faceMids = midPointToFaceMids[edgeMidPointi];
1030
1031 label otherFaceMidPointi = faceMids.otherVertex(faceMidPointi);
1032
1033 // Create face consistent with anchorI being the owner.
1034 // Note that the edges between the edge mid point and the face mids
1035 // might be marked for splitting. Note that these edge splits cannot
1036 // be between cellMid and face mids.
1037
1038 DynamicList<label> newFaceVerts(4);
1039 if (faceOrder == (mesh_.faceOwner()[facei] == celli))
1040 {
1041 newFaceVerts.append(faceMidPointi);
1042
1043 // Check & insert edge split if any
1044 insertEdgeSplit
1045 (
1046 edgeMidPoint,
1047 faceMidPointi, // edge between faceMid and
1048 edgeMidPointi, // edgeMid
1049 newFaceVerts
1050 );
1051
1052 newFaceVerts.append(edgeMidPointi);
1053
1054 insertEdgeSplit
1055 (
1056 edgeMidPoint,
1057 edgeMidPointi,
1058 otherFaceMidPointi,
1059 newFaceVerts
1060 );
1061
1062 newFaceVerts.append(otherFaceMidPointi);
1063 newFaceVerts.append(cellMidPoint[celli]);
1064 }
1065 else
1066 {
1067 newFaceVerts.append(otherFaceMidPointi);
1068
1069 insertEdgeSplit
1070 (
1071 edgeMidPoint,
1072 otherFaceMidPointi,
1073 edgeMidPointi,
1074 newFaceVerts
1075 );
1076
1077 newFaceVerts.append(edgeMidPointi);
1078
1079 insertEdgeSplit
1080 (
1081 edgeMidPoint,
1082 edgeMidPointi,
1083 faceMidPointi,
1084 newFaceVerts
1085 );
1086
1087 newFaceVerts.append(faceMidPointi);
1088 newFaceVerts.append(cellMidPoint[celli]);
1089 }
1090
1091 face newFace;
1092 newFace.transfer(newFaceVerts);
1093
1094 label anchorCell0 = getAnchorCell
1095 (
1096 cellAnchorPoints,
1097 cellAddedCells,
1098 celli,
1099 facei,
1100 anchorPointi
1101 );
1102 label anchorCell1 = getAnchorCell
1103 (
1104 cellAnchorPoints,
1105 cellAddedCells,
1106 celli,
1107 facei,
1108 anchors.otherVertex(anchorPointi)
1109 );
1110
1111
1112 label own, nei;
1113 point ownPt, neiPt;
1114
1115 if (anchorCell0 < anchorCell1)
1116 {
1117 own = anchorCell0;
1118 nei = anchorCell1;
1119
1120 ownPt = mesh_.points()[anchorPointi];
1121 neiPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1122
1123 }
1124 else
1125 {
1126 own = anchorCell1;
1127 nei = anchorCell0;
1128 newFace.flip();
1129
1130 ownPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1131 neiPt = mesh_.points()[anchorPointi];
1132 }
1133
1134 if (debug)
1135 {
1136 point ownPt, neiPt;
1137
1138 if (anchorCell0 < anchorCell1)
1139 {
1140 ownPt = mesh_.points()[anchorPointi];
1141 neiPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1142 }
1143 else
1144 {
1145 ownPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1146 neiPt = mesh_.points()[anchorPointi];
1147 }
1148
1149 checkInternalOrientation
1150 (
1151 meshMod,
1152 celli,
1153 facei,
1154 ownPt,
1155 neiPt,
1156 newFace
1157 );
1158 }
1159
1160 return addInternalFace
1161 (
1162 meshMod,
1163 facei,
1164 anchorPointi,
1165 newFace,
1166 own,
1167 nei
1168 );
1169 }
1170 else
1171 {
1172 return -1;
1173 }
1174}
1175
1176
1177// Creates all the 12 internal faces for celli.
1178void Foam::hexRef8::createInternalFaces
1179(
1180 const labelListList& cellAnchorPoints,
1181 const labelListList& cellAddedCells,
1182 const labelList& cellMidPoint,
1183 const labelList& faceMidPoint,
1184 const labelList& faceAnchorLevel,
1185 const labelList& edgeMidPoint,
1186 const label celli,
1187
1188 polyTopoChange& meshMod
1189) const
1190{
1191 // Find in every face the cellLevel+1 points (from edge subdivision)
1192 // and the anchor points.
1193
1194 const cell& cFaces = mesh_.cells()[celli];
1195 const label cLevel = cellLevel_[celli];
1196
1197 // From edge mid to anchor points
1198 Map<edge> midPointToAnchors(24);
1199 // From edge mid to face mids
1200 Map<edge> midPointToFaceMids(24);
1201
1202 // Storage for on-the-fly addressing
1203 DynamicList<label> storage;
1204
1205
1206 // Running count of number of internal faces added so far.
1207 label nFacesAdded = 0;
1208
1209 forAll(cFaces, i)
1210 {
1211 label facei = cFaces[i];
1212
1213 const face& f = mesh_.faces()[facei];
1214 const labelList& fEdges = mesh_.faceEdges(facei, storage);
1215
1216 // We are on the celli side of face f. The face will have 1 or 4
1217 // cLevel points and lots of higher numbered ones.
1218
1219 label faceMidPointi = -1;
1220
1221 label nAnchors = countAnchors(f, cLevel);
1222
1223 if (nAnchors == 1)
1224 {
1225 // Only one anchor point. So the other side of the face has already
1226 // been split using cLevel+1 and cLevel+2 points.
1227
1228 // Find the one anchor.
1229 label anchorFp = -1;
1230
1231 forAll(f, fp)
1232 {
1233 if (pointLevel_[f[fp]] <= cLevel)
1234 {
1235 anchorFp = fp;
1236 break;
1237 }
1238 }
1239
1240 // Now the face mid point is the second cLevel+1 point
1241 label edgeMid = findLevel
1242 (
1243 facei,
1244 f,
1245 f.fcIndex(anchorFp),
1246 true,
1247 cLevel+1
1248 );
1249 label faceMid = findLevel
1250 (
1251 facei,
1252 f,
1253 f.fcIndex(edgeMid),
1254 true,
1255 cLevel+1
1256 );
1257
1258 faceMidPointi = f[faceMid];
1259 }
1260 else if (nAnchors == 4)
1261 {
1262 // There is no face middle yet but the face will be marked for
1263 // splitting.
1264
1265 faceMidPointi = faceMidPoint[facei];
1266 }
1267 else
1268 {
1269 dumpCell(mesh_.faceOwner()[facei]);
1270 if (mesh_.isInternalFace(facei))
1271 {
1272 dumpCell(mesh_.faceNeighbour()[facei]);
1273 }
1274
1276 << "nAnchors:" << nAnchors
1277 << " facei:" << facei
1278 << abort(FatalError);
1279 }
1280
1281
1282
1283 // Now loop over all the anchors (might be just one) and store
1284 // the edge mids connected to it. storeMidPointInfo will collect
1285 // all the info and combine it all.
1286
1287 forAll(f, fp0)
1288 {
1289 label point0 = f[fp0];
1290
1291 if (pointLevel_[point0] <= cLevel)
1292 {
1293 // Anchor.
1294
1295 // Walk forward
1296 // ~~~~~~~~~~~~
1297 // to cLevel+1 or edgeMidPoint of this level.
1298
1299
1300 label edgeMidPointi = -1;
1301
1302 label fp1 = f.fcIndex(fp0);
1303
1304 if (pointLevel_[f[fp1]] <= cLevel)
1305 {
1306 // Anchor. Edge will be split.
1307 label edgeI = fEdges[fp0];
1308
1309 edgeMidPointi = edgeMidPoint[edgeI];
1310
1311 if (edgeMidPointi == -1)
1312 {
1313 dumpCell(celli);
1314
1315 const labelList& cPoints = mesh_.cellPoints(celli);
1316
1318 << "cell:" << celli << " cLevel:" << cLevel
1319 << " cell points:" << cPoints
1320 << " pointLevel:"
1321 << labelUIndList(pointLevel_, cPoints)
1322 << " face:" << facei
1323 << " f:" << f
1324 << " pointLevel:"
1325 << labelUIndList(pointLevel_, f)
1326 << " faceAnchorLevel:" << faceAnchorLevel[facei]
1327 << " faceMidPoint:" << faceMidPoint[facei]
1328 << " faceMidPointi:" << faceMidPointi
1329 << " fp:" << fp0
1330 << abort(FatalError);
1331 }
1332 }
1333 else
1334 {
1335 // Search forward in face to clevel+1
1336 label edgeMid = findLevel(facei, f, fp1, true, cLevel+1);
1337
1338 edgeMidPointi = f[edgeMid];
1339 }
1340
1341 label newFacei = storeMidPointInfo
1342 (
1343 cellAnchorPoints,
1344 cellAddedCells,
1345 cellMidPoint,
1346 edgeMidPoint,
1347
1348 celli,
1349 facei,
1350 true, // mid point after anchor
1351 edgeMidPointi, // edgemid
1352 point0, // anchor
1353 faceMidPointi,
1354
1355 midPointToAnchors,
1356 midPointToFaceMids,
1357 meshMod
1358 );
1359
1360 if (newFacei != -1)
1361 {
1362 nFacesAdded++;
1363
1364 if (nFacesAdded == 12)
1365 {
1366 break;
1367 }
1368 }
1369
1370
1371
1372 // Walk backward
1373 // ~~~~~~~~~~~~~
1374
1375 label fpMin1 = f.rcIndex(fp0);
1376
1377 if (pointLevel_[f[fpMin1]] <= cLevel)
1378 {
1379 // Anchor. Edge will be split.
1380 label edgeI = fEdges[fpMin1];
1381
1382 edgeMidPointi = edgeMidPoint[edgeI];
1383
1384 if (edgeMidPointi == -1)
1385 {
1386 dumpCell(celli);
1387
1388 const labelList& cPoints = mesh_.cellPoints(celli);
1389
1391 << "cell:" << celli << " cLevel:" << cLevel
1392 << " cell points:" << cPoints
1393 << " pointLevel:"
1394 << labelUIndList(pointLevel_, cPoints)
1395 << " face:" << facei
1396 << " f:" << f
1397 << " pointLevel:"
1398 << labelUIndList(pointLevel_, f)
1399 << " faceAnchorLevel:" << faceAnchorLevel[facei]
1400 << " faceMidPoint:" << faceMidPoint[facei]
1401 << " faceMidPointi:" << faceMidPointi
1402 << " fp:" << fp0
1403 << abort(FatalError);
1404 }
1405 }
1406 else
1407 {
1408 // Search back to clevel+1
1409 label edgeMid = findLevel
1410 (
1411 facei,
1412 f,
1413 fpMin1,
1414 false,
1415 cLevel+1
1416 );
1417
1418 edgeMidPointi = f[edgeMid];
1419 }
1420
1421 newFacei = storeMidPointInfo
1422 (
1423 cellAnchorPoints,
1424 cellAddedCells,
1425 cellMidPoint,
1426 edgeMidPoint,
1427
1428 celli,
1429 facei,
1430 false, // mid point before anchor
1431 edgeMidPointi, // edgemid
1432 point0, // anchor
1433 faceMidPointi,
1434
1435 midPointToAnchors,
1436 midPointToFaceMids,
1437 meshMod
1438 );
1439
1440 if (newFacei != -1)
1441 {
1442 nFacesAdded++;
1443
1444 if (nFacesAdded == 12)
1445 {
1446 break;
1447 }
1448 }
1449 } // done anchor
1450 } // done face
1451
1452 if (nFacesAdded == 12)
1453 {
1454 break;
1455 }
1456 }
1457}
1458
1459
1460void Foam::hexRef8::walkFaceToMid
1461(
1462 const labelList& edgeMidPoint,
1463 const label cLevel,
1464 const label facei,
1465 const label startFp,
1466 DynamicList<label>& faceVerts
1467) const
1468{
1469 const face& f = mesh_.faces()[facei];
1470 const labelList& fEdges = mesh_.faceEdges(facei);
1471
1472 label fp = startFp;
1473
1474 // Starting from fp store all (1 or 2) vertices until where the face
1475 // gets split
1476
1477 while (true)
1478 {
1479 if (edgeMidPoint[fEdges[fp]] >= 0)
1480 {
1481 faceVerts.append(edgeMidPoint[fEdges[fp]]);
1482 }
1483
1484 fp = f.fcIndex(fp);
1485
1486 if (pointLevel_[f[fp]] <= cLevel)
1487 {
1488 // Next anchor. Have already append split point on edge in code
1489 // above.
1490 return;
1491 }
1492 else if (pointLevel_[f[fp]] == cLevel+1)
1493 {
1494 // Mid level
1495 faceVerts.append(f[fp]);
1496
1497 return;
1498 }
1499 else if (pointLevel_[f[fp]] == cLevel+2)
1500 {
1501 // Store and continue to cLevel+1.
1502 faceVerts.append(f[fp]);
1503 }
1504 }
1505}
1506
1507
1508// Same as walkFaceToMid but now walk back.
1509void Foam::hexRef8::walkFaceFromMid
1510(
1511 const labelList& edgeMidPoint,
1512 const label cLevel,
1513 const label facei,
1514 const label startFp,
1515 DynamicList<label>& faceVerts
1516) const
1517{
1518 const face& f = mesh_.faces()[facei];
1519 const labelList& fEdges = mesh_.faceEdges(facei);
1520
1521 label fp = f.rcIndex(startFp);
1522
1523 while (true)
1524 {
1525 if (pointLevel_[f[fp]] <= cLevel)
1526 {
1527 // anchor.
1528 break;
1529 }
1530 else if (pointLevel_[f[fp]] == cLevel+1)
1531 {
1532 // Mid level
1533 faceVerts.append(f[fp]);
1534 break;
1535 }
1536 else if (pointLevel_[f[fp]] == cLevel+2)
1537 {
1538 // Continue to cLevel+1.
1539 }
1540 fp = f.rcIndex(fp);
1541 }
1542
1543 // Store
1544 while (true)
1545 {
1546 if (edgeMidPoint[fEdges[fp]] >= 0)
1547 {
1548 faceVerts.append(edgeMidPoint[fEdges[fp]]);
1549 }
1550
1551 fp = f.fcIndex(fp);
1552
1553 if (fp == startFp)
1554 {
1555 break;
1556 }
1557 faceVerts.append(f[fp]);
1558 }
1559}
1560
1561
1562// Updates refineCell (cells marked for refinement) so across all faces
1563// there will be 2:1 consistency after refinement.
1564Foam::label Foam::hexRef8::faceConsistentRefinement
1565(
1566 const bool maxSet,
1567 const labelUList& cellLevel,
1568 bitSet& refineCell
1569) const
1570{
1571 label nChanged = 0;
1572
1573 // Internal faces.
1574 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
1575 {
1576 label own = mesh_.faceOwner()[facei];
1577 label ownLevel = cellLevel[own] + refineCell.get(own);
1578
1579 label nei = mesh_.faceNeighbour()[facei];
1580 label neiLevel = cellLevel[nei] + refineCell.get(nei);
1581
1582 if (ownLevel > (neiLevel+1))
1583 {
1584 if (maxSet)
1585 {
1586 refineCell.set(nei);
1587 }
1588 else
1589 {
1590 refineCell.unset(own);
1591 }
1592 nChanged++;
1593 }
1594 else if (neiLevel > (ownLevel+1))
1595 {
1596 if (maxSet)
1597 {
1598 refineCell.set(own);
1599 }
1600 else
1601 {
1602 refineCell.unset(nei);
1603 }
1604 nChanged++;
1605 }
1606 }
1607
1608
1609 // Coupled faces. Swap owner level to get neighbouring cell level.
1610 // (only boundary faces of neiLevel used)
1611 labelList neiLevel(mesh_.nBoundaryFaces());
1612
1613 forAll(neiLevel, i)
1614 {
1615 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1616
1617 neiLevel[i] = cellLevel[own] + refineCell.get(own);
1618 }
1619
1620 // Swap to neighbour
1621 syncTools::swapBoundaryFaceList(mesh_, neiLevel);
1622
1623 // Now we have neighbour value see which cells need refinement
1624 forAll(neiLevel, i)
1625 {
1626 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1627 label ownLevel = cellLevel[own] + refineCell.get(own);
1628
1629 if (ownLevel > (neiLevel[i]+1))
1630 {
1631 if (!maxSet)
1632 {
1633 refineCell.unset(own);
1634 nChanged++;
1635 }
1636 }
1637 else if (neiLevel[i] > (ownLevel+1))
1638 {
1639 if (maxSet)
1640 {
1641 refineCell.set(own);
1642 nChanged++;
1643 }
1644 }
1645 }
1646
1647 return nChanged;
1648}
1649
1650
1651// Debug: check if wanted refinement is compatible with 2:1
1652void Foam::hexRef8::checkWantedRefinementLevels
1653(
1654 const labelUList& cellLevel,
1655 const labelList& cellsToRefine
1656) const
1657{
1658 bitSet refineCell(mesh_.nCells(), cellsToRefine);
1659
1660 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
1661 {
1662 label own = mesh_.faceOwner()[facei];
1663 label ownLevel = cellLevel[own] + refineCell.get(own);
1664
1665 label nei = mesh_.faceNeighbour()[facei];
1666 label neiLevel = cellLevel[nei] + refineCell.get(nei);
1667
1668 if (mag(ownLevel-neiLevel) > 1)
1669 {
1670 dumpCell(own);
1671 dumpCell(nei);
1673 << "cell:" << own
1674 << " current level:" << cellLevel[own]
1675 << " level after refinement:" << ownLevel
1676 << nl
1677 << "neighbour cell:" << nei
1678 << " current level:" << cellLevel[nei]
1679 << " level after refinement:" << neiLevel
1680 << nl
1681 << "which does not satisfy 2:1 constraints anymore."
1682 << abort(FatalError);
1683 }
1684 }
1685
1686 // Coupled faces. Swap owner level to get neighbouring cell level.
1687 // (only boundary faces of neiLevel used)
1688 labelList neiLevel(mesh_.nBoundaryFaces());
1689
1690 forAll(neiLevel, i)
1691 {
1692 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1693
1694 neiLevel[i] = cellLevel[own] + refineCell.get(own);
1695 }
1696
1697 // Swap to neighbour
1698 syncTools::swapBoundaryFaceList(mesh_, neiLevel);
1699
1700 // Now we have neighbour value see which cells need refinement
1701 forAll(neiLevel, i)
1702 {
1703 label facei = i + mesh_.nInternalFaces();
1704
1705 label own = mesh_.faceOwner()[facei];
1706 label ownLevel = cellLevel[own] + refineCell.get(own);
1707
1708 if (mag(ownLevel - neiLevel[i]) > 1)
1709 {
1710 label patchi = mesh_.boundaryMesh().whichPatch(facei);
1711
1712 dumpCell(own);
1714 << "Celllevel does not satisfy 2:1 constraint."
1715 << " On coupled face "
1716 << facei
1717 << " on patch " << patchi << " "
1718 << mesh_.boundaryMesh()[patchi].name()
1719 << " owner cell " << own
1720 << " current level:" << cellLevel[own]
1721 << " level after refinement:" << ownLevel
1722 << nl
1723 << " (coupled) neighbour cell will get refinement "
1724 << neiLevel[i]
1725 << abort(FatalError);
1726 }
1727 }
1728}
1729
1730
1731// Set instance for mesh files
1733{
1734 if (debug)
1735 {
1736 Pout<< "hexRef8::setInstance(const fileName& inst) : "
1737 << "Resetting file instance to " << inst << endl;
1738 }
1739
1740 cellLevel_.instance() = inst;
1741 pointLevel_.instance() = inst;
1742 level0Edge_.instance() = inst;
1743 history_.instance() = inst;
1744}
1745
1746
1747void Foam::hexRef8::collectLevelPoints
1748(
1749 const labelList& f,
1750 const label level,
1752) const
1753{
1754 forAll(f, fp)
1755 {
1756 if (pointLevel_[f[fp]] <= level)
1757 {
1758 points.append(f[fp]);
1759 }
1760 }
1761}
1762
1763
1764void Foam::hexRef8::collectLevelPoints
1765(
1766 const labelList& meshPoints,
1767 const labelList& f,
1768 const label level,
1769 DynamicList<label>& points
1770) const
1771{
1772 forAll(f, fp)
1773 {
1774 label pointi = meshPoints[f[fp]];
1775 if (pointLevel_[pointi] <= level)
1776 {
1777 points.append(pointi);
1778 }
1779 }
1780}
1781
1782
1783// Return true if we've found 6 quads. faces guaranteed to be outwards pointing.
1784bool Foam::hexRef8::matchHexShape
1785(
1786 const label celli,
1787 const label cellLevel,
1788 DynamicList<face>& quads
1789) const
1790{
1791 const cell& cFaces = mesh_.cells()[celli];
1792
1793 // Work arrays
1794 DynamicList<label> verts(4);
1795 quads.clear();
1796
1797
1798 // 1. pick up any faces with four cellLevel points
1799
1800 forAll(cFaces, i)
1801 {
1802 label facei = cFaces[i];
1803 const face& f = mesh_.faces()[facei];
1804
1805 verts.clear();
1806 collectLevelPoints(f, cellLevel, verts);
1807 if (verts.size() == 4)
1808 {
1809 if (mesh_.faceOwner()[facei] != celli)
1810 {
1811 reverse(verts);
1812 }
1813 quads.append(face(0));
1814 labelList& quadVerts = quads.last();
1815 quadVerts.transfer(verts);
1816 }
1817 }
1818
1819
1820 if (quads.size() < 6)
1821 {
1822 Map<labelList> pointFaces(2*cFaces.size());
1823
1824 forAll(cFaces, i)
1825 {
1826 label facei = cFaces[i];
1827 const face& f = mesh_.faces()[facei];
1828
1829 // Pick up any faces with only one level point.
1830 // See if there are four of these where the common point
1831 // is a level+1 point. This common point is then the mid of
1832 // a split face.
1833
1834 verts.clear();
1835 collectLevelPoints(f, cellLevel, verts);
1836 if (verts.size() == 1)
1837 {
1838 // Add to pointFaces for any level+1 point (this might be
1839 // a midpoint of a split face)
1840 forAll(f, fp)
1841 {
1842 label pointi = f[fp];
1843 if (pointLevel_[pointi] == cellLevel+1)
1844 {
1845 auto iter = pointFaces.find(pointi);
1846
1847 if (iter.found())
1848 {
1849 labelList& pFaces = iter.val();
1850 pFaces.appendUniq(facei);
1851 }
1852 else
1853 {
1854 pointFaces.insert
1855 (
1856 pointi,
1857 labelList(one{}, facei)
1858 );
1859 }
1860 }
1861 }
1862 }
1863 }
1864
1865 // 2. Check if we've collected any midPoints.
1866 forAllConstIters(pointFaces, iter)
1867 {
1868 const labelList& pFaces = iter.val();
1869
1870 if (pFaces.size() == 4)
1871 {
1872 // Collect and orient.
1873 faceList fourFaces(pFaces.size());
1874 forAll(pFaces, pFacei)
1875 {
1876 label facei = pFaces[pFacei];
1877 const face& f = mesh_.faces()[facei];
1878 if (mesh_.faceOwner()[facei] == celli)
1879 {
1880 fourFaces[pFacei] = f;
1881 }
1882 else
1883 {
1884 fourFaces[pFacei] = f.reverseFace();
1885 }
1886 }
1887
1888 primitivePatch bigFace
1889 (
1890 SubList<face>(fourFaces),
1891 mesh_.points()
1892 );
1893 const labelListList& edgeLoops = bigFace.edgeLoops();
1894
1895 if (edgeLoops.size() == 1)
1896 {
1897 // Collect the 4 cellLevel points
1898 verts.clear();
1899 collectLevelPoints
1900 (
1901 bigFace.meshPoints(),
1902 bigFace.edgeLoops()[0],
1903 cellLevel,
1904 verts
1905 );
1906
1907 if (verts.size() == 4)
1908 {
1909 quads.append(face(0));
1910 labelList& quadVerts = quads.last();
1911 quadVerts.transfer(verts);
1912 }
1913 }
1914 }
1915 }
1916 }
1917
1918 return (quads.size() == 6);
1919}
1920
1921
1922// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1923
1924// Construct from mesh, read refinement data
1925Foam::hexRef8::hexRef8(const polyMesh& mesh, const bool readHistory)
1926:
1927 mesh_(mesh),
1928 cellLevel_
1929 (
1930 IOobject
1931 (
1932 "cellLevel",
1933 mesh_.facesInstance(),
1934 polyMesh::meshSubDir,
1935 mesh_,
1936 IOobject::READ_IF_PRESENT,
1937 IOobject::NO_WRITE
1938 ),
1939 labelList(mesh_.nCells(), Zero)
1940 ),
1941 pointLevel_
1942 (
1943 IOobject
1944 (
1945 "pointLevel",
1946 mesh_.facesInstance(),
1947 polyMesh::meshSubDir,
1948 mesh_,
1949 IOobject::READ_IF_PRESENT,
1950 IOobject::NO_WRITE
1951 ),
1952 labelList(mesh_.nPoints(), Zero)
1953 ),
1954 level0Edge_
1955 (
1956 IOobject
1957 (
1958 "level0Edge",
1959 mesh_.facesInstance(),
1960 polyMesh::meshSubDir,
1961 mesh_,
1962 IOobject::READ_IF_PRESENT,
1963 IOobject::NO_WRITE
1964 ),
1965 // Needs name:
1966 dimensionedScalar("level0Edge", dimLength, getLevel0EdgeLength())
1967 ),
1968 history_
1969 (
1970 IOobject
1971 (
1972 "refinementHistory",
1973 mesh_.facesInstance(),
1974 polyMesh::meshSubDir,
1975 mesh_,
1976 IOobject::NO_READ,
1977 IOobject::NO_WRITE
1978 ),
1979 // All cells visible if not read or readHistory = false
1980 (readHistory ? mesh_.nCells() : 0)
1981 ),
1982 faceRemover_(mesh_, GREAT), // merge boundary faces wherever possible
1983 savedPointLevel_(0),
1984 savedCellLevel_(0)
1985{
1986 if (readHistory)
1987 {
1989 if (history_.typeHeaderOk<refinementHistory>(true))
1990 {
1991 history_.read();
1992 }
1993 }
1994
1995 if (history_.active() && history_.visibleCells().size() != mesh_.nCells())
1996 {
1998 << "History enabled but number of visible cells "
1999 << history_.visibleCells().size() << " in "
2000 << history_.objectPath()
2001 << " is not equal to the number of cells in the mesh "
2002 << mesh_.nCells()
2003 << abort(FatalError);
2004 }
2005
2006 if
2007 (
2008 cellLevel_.size() != mesh_.nCells()
2009 || pointLevel_.size() != mesh_.nPoints()
2010 )
2011 {
2013 << "Restarted from inconsistent cellLevel or pointLevel files."
2014 << endl
2015 << "cellLevel file " << cellLevel_.objectPath() << endl
2016 << "pointLevel file " << pointLevel_.objectPath() << endl
2017 << "Number of cells in mesh:" << mesh_.nCells()
2018 << " does not equal size of cellLevel:" << cellLevel_.size() << endl
2019 << "Number of points in mesh:" << mesh_.nPoints()
2020 << " does not equal size of pointLevel:" << pointLevel_.size()
2021 << abort(FatalError);
2022 }
2023
2024
2025 // Check refinement levels for consistency
2027
2028
2029 // Check initial mesh for consistency
2030
2031 //if (debug)
2032 {
2033 checkMesh();
2034 }
2035}
2036
2037
2039(
2040 const polyMesh& mesh,
2041 const labelList& cellLevel,
2042 const labelList& pointLevel,
2043 const refinementHistory& history,
2044 const scalar level0Edge
2045)
2046:
2047 mesh_(mesh),
2048 cellLevel_
2049 (
2050 IOobject
2051 (
2052 "cellLevel",
2053 mesh_.facesInstance(),
2054 polyMesh::meshSubDir,
2055 mesh_,
2056 IOobject::NO_READ,
2057 IOobject::NO_WRITE
2058 ),
2059 cellLevel
2060 ),
2061 pointLevel_
2062 (
2063 IOobject
2064 (
2065 "pointLevel",
2066 mesh_.facesInstance(),
2067 polyMesh::meshSubDir,
2068 mesh_,
2069 IOobject::NO_READ,
2070 IOobject::NO_WRITE
2071 ),
2072 pointLevel
2073 ),
2074 level0Edge_
2075 (
2076 IOobject
2077 (
2078 "level0Edge",
2079 mesh_.facesInstance(),
2080 polyMesh::meshSubDir,
2081 mesh_,
2082 IOobject::NO_READ,
2083 IOobject::NO_WRITE
2084 ),
2085 // Needs name:
2087 (
2088 "level0Edge",
2089 dimLength,
2090 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2091 )
2092 ),
2093 history_
2094 (
2095 IOobject
2096 (
2097 "refinementHistory",
2098 mesh_.facesInstance(),
2099 polyMesh::meshSubDir,
2100 mesh_,
2101 IOobject::NO_READ,
2102 IOobject::NO_WRITE
2103 ),
2104 history
2105 ),
2106 faceRemover_(mesh_, GREAT), // merge boundary faces wherever possible
2107 savedPointLevel_(0),
2108 savedCellLevel_(0)
2109{
2110 if (history_.active() && history_.visibleCells().size() != mesh_.nCells())
2111 {
2113 << "History enabled but number of visible cells in it "
2114 << history_.visibleCells().size()
2115 << " is not equal to the number of cells in the mesh "
2116 << mesh_.nCells() << abort(FatalError);
2117 }
2118
2119 if
2120 (
2121 cellLevel_.size() != mesh_.nCells()
2122 || pointLevel_.size() != mesh_.nPoints()
2123 )
2124 {
2126 << "Incorrect cellLevel or pointLevel size." << endl
2127 << "Number of cells in mesh:" << mesh_.nCells()
2128 << " does not equal size of cellLevel:" << cellLevel_.size() << endl
2129 << "Number of points in mesh:" << mesh_.nPoints()
2130 << " does not equal size of pointLevel:" << pointLevel_.size()
2131 << abort(FatalError);
2132 }
2133
2134 // Check refinement levels for consistency
2136
2137
2138 // Check initial mesh for consistency
2139
2140 //if (debug)
2141 {
2142 checkMesh();
2143 }
2144}
2145
2146
2148(
2149 const polyMesh& mesh,
2150 const labelList& cellLevel,
2151 const labelList& pointLevel,
2152 const scalar level0Edge
2153)
2154:
2155 mesh_(mesh),
2156 cellLevel_
2157 (
2158 IOobject
2159 (
2160 "cellLevel",
2161 mesh_.facesInstance(),
2162 polyMesh::meshSubDir,
2163 mesh_,
2164 IOobject::NO_READ,
2165 IOobject::NO_WRITE
2166 ),
2167 cellLevel
2168 ),
2169 pointLevel_
2170 (
2171 IOobject
2172 (
2173 "pointLevel",
2174 mesh_.facesInstance(),
2175 polyMesh::meshSubDir,
2176 mesh_,
2177 IOobject::NO_READ,
2178 IOobject::NO_WRITE
2179 ),
2180 pointLevel
2181 ),
2182 level0Edge_
2183 (
2184 IOobject
2185 (
2186 "level0Edge",
2187 mesh_.facesInstance(),
2188 polyMesh::meshSubDir,
2189 mesh_,
2190 IOobject::NO_READ,
2191 IOobject::NO_WRITE
2192 ),
2193 // Needs name:
2195 (
2196 "level0Edge",
2197 dimLength,
2198 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2199 )
2200 ),
2201 history_
2202 (
2203 IOobject
2204 (
2205 "refinementHistory",
2206 mesh_.facesInstance(),
2207 polyMesh::meshSubDir,
2208 mesh_,
2209 IOobject::NO_READ,
2210 IOobject::NO_WRITE
2211 ),
2212 List<refinementHistory::splitCell8>(0),
2213 labelList(0),
2214 false
2215 ),
2216 faceRemover_(mesh_, GREAT), // merge boundary faces wherever possible
2217 savedPointLevel_(0),
2218 savedCellLevel_(0)
2219{
2220 if
2221 (
2222 cellLevel_.size() != mesh_.nCells()
2223 || pointLevel_.size() != mesh_.nPoints()
2224 )
2225 {
2227 << "Incorrect cellLevel or pointLevel size." << endl
2228 << "Number of cells in mesh:" << mesh_.nCells()
2229 << " does not equal size of cellLevel:" << cellLevel_.size() << endl
2230 << "Number of points in mesh:" << mesh_.nPoints()
2231 << " does not equal size of pointLevel:" << pointLevel_.size()
2232 << abort(FatalError);
2233 }
2234
2235 // Check refinement levels for consistency
2237
2238 // Check initial mesh for consistency
2239
2240 //if (debug)
2241 {
2242 checkMesh();
2243 }
2244}
2245
2246
2247// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2248
2250(
2251 const labelUList& cellLevel,
2252 const labelList& cellsToRefine,
2253 const bool maxSet
2254) const
2255{
2256 // Loop, modifying cellsToRefine, until no more changes to due to 2:1
2257 // conflicts.
2258 // maxSet = false : unselect cells to refine
2259 // maxSet = true : select cells to refine
2260
2261 bitSet refineCell(mesh_.nCells(), cellsToRefine);
2262
2263 while (true)
2264 {
2265 label nChanged = faceConsistentRefinement
2266 (
2267 maxSet,
2268 cellLevel,
2270 );
2271
2272 reduce(nChanged, sumOp<label>());
2273
2274 if (debug)
2275 {
2276 Pout<< "hexRef8::consistentRefinement : Changed " << nChanged
2277 << " refinement levels due to 2:1 conflicts."
2278 << endl;
2279 }
2280
2281 if (nChanged == 0)
2282 {
2283 break;
2284 }
2285 }
2286
2287 // Convert back to labelList.
2288 labelList newCellsToRefine(refineCell.toc());
2289
2290 if (debug)
2291 {
2292 checkWantedRefinementLevels(cellLevel, newCellsToRefine);
2293 }
2294
2295 return newCellsToRefine;
2296}
2297
2298
2299// Given a list of cells to refine determine additional cells to refine
2300// such that the overall refinement:
2301// - satisfies maxFaceDiff (e.g. 2:1) across neighbouring faces
2302// - satisfies maxPointDiff (e.g. 4:1) across selected point connected
2303// cells. This is used to ensure that e.g. cells on the surface are not
2304// point connected to cells which are 8 times smaller.
2306(
2307 const label maxFaceDiff,
2308 const labelList& cellsToRefine,
2309 const labelList& facesToCheck,
2310 const label maxPointDiff,
2311 const labelList& pointsToCheck
2312) const
2313{
2314 const labelList& faceOwner = mesh_.faceOwner();
2315 const labelList& faceNeighbour = mesh_.faceNeighbour();
2316
2317
2318 if (maxFaceDiff <= 0)
2319 {
2321 << "Illegal maxFaceDiff " << maxFaceDiff << nl
2322 << "Value should be >= 1" << exit(FatalError);
2323 }
2324
2325
2326 // Bit tricky. Say we want a distance of three cells between two
2327 // consecutive refinement levels. This is done by using FaceCellWave to
2328 // transport out the new refinement level. It gets decremented by one
2329 // every cell it crosses so if we initialize it to maxFaceDiff
2330 // we will get a field everywhere that tells us whether an unselected cell
2331 // needs refining as well.
2332
2333
2334 // Initial information about (distance to) cellLevel on all cells
2335 List<refinementData> allCellInfo(mesh_.nCells());
2336
2337 // Initial information about (distance to) cellLevel on all faces
2338 List<refinementData> allFaceInfo(mesh_.nFaces());
2339
2340 forAll(allCellInfo, celli)
2341 {
2342 // maxFaceDiff since refinementData counts both
2343 // faces and cells.
2344 allCellInfo[celli] = refinementData
2345 (
2346 maxFaceDiff*(cellLevel_[celli]+1),// when cell is to be refined
2347 maxFaceDiff*cellLevel_[celli] // current level
2348 );
2349 }
2350
2351 // Cells to be refined will have cellLevel+1
2352 forAll(cellsToRefine, i)
2353 {
2354 label celli = cellsToRefine[i];
2355
2356 allCellInfo[celli].count() = allCellInfo[celli].refinementCount();
2357 }
2358
2359
2360 // Labels of seed faces
2361 DynamicList<label> seedFaces(mesh_.nFaces()/100);
2362 // refinementLevel data on seed faces
2363 DynamicList<refinementData> seedFacesInfo(mesh_.nFaces()/100);
2364
2365 // Dummy additional info for FaceCellWave
2366 int dummyTrackData = 0;
2367
2368
2369 // Additional buffer layer thickness by changing initial count. Usually
2370 // this happens on boundary faces. Bit tricky. Use allFaceInfo to mark
2371 // off thus marked faces so they're skipped in the next loop.
2372 forAll(facesToCheck, i)
2373 {
2374 label facei = facesToCheck[i];
2375
2376 if (allFaceInfo[facei].valid(dummyTrackData))
2377 {
2378 // Can only occur if face has already gone through loop below.
2380 << "Argument facesToCheck seems to have duplicate entries!"
2381 << endl
2382 << "face:" << facei << " occurs at positions "
2383 << findIndices(facesToCheck, facei)
2384 << abort(FatalError);
2385 }
2386
2387
2388 const refinementData& ownData = allCellInfo[faceOwner[facei]];
2389
2390 if (mesh_.isInternalFace(facei))
2391 {
2392 // Seed face if neighbouring cell (after possible refinement)
2393 // will be refined one more than the current owner or neighbour.
2394
2395 const refinementData& neiData = allCellInfo[faceNeighbour[facei]];
2396
2397 label faceCount;
2398 label faceRefineCount;
2399 if (neiData.count() > ownData.count())
2400 {
2401 faceCount = neiData.count() + maxFaceDiff;
2402 faceRefineCount = faceCount + maxFaceDiff;
2403 }
2404 else
2405 {
2406 faceCount = ownData.count() + maxFaceDiff;
2407 faceRefineCount = faceCount + maxFaceDiff;
2408 }
2409
2410 seedFaces.append(facei);
2411 seedFacesInfo.append
2412 (
2414 (
2415 faceRefineCount,
2416 faceCount
2417 )
2418 );
2419 allFaceInfo[facei] = seedFacesInfo.last();
2420 }
2421 else
2422 {
2423 label faceCount = ownData.count() + maxFaceDiff;
2424 label faceRefineCount = faceCount + maxFaceDiff;
2425
2426 seedFaces.append(facei);
2427 seedFacesInfo.append
2428 (
2430 (
2431 faceRefineCount,
2432 faceCount
2433 )
2434 );
2435 allFaceInfo[facei] = seedFacesInfo.last();
2436 }
2437 }
2438
2439
2440 // Just seed with all faces inbetween different refinement levels for now
2441 // (alternatively only seed faces on cellsToRefine but that gives problems
2442 // if no cells to refine)
2443 forAll(faceNeighbour, facei)
2444 {
2445 // Check if face already handled in loop above
2446 if (!allFaceInfo[facei].valid(dummyTrackData))
2447 {
2448 label own = faceOwner[facei];
2449 label nei = faceNeighbour[facei];
2450
2451 // Seed face with transported data from highest cell.
2452
2453 if (allCellInfo[own].count() > allCellInfo[nei].count())
2454 {
2455 allFaceInfo[facei].updateFace
2456 (
2457 mesh_,
2458 facei,
2459 own,
2460 allCellInfo[own],
2462 dummyTrackData
2463 );
2464 seedFaces.append(facei);
2465 seedFacesInfo.append(allFaceInfo[facei]);
2466 }
2467 else if (allCellInfo[own].count() < allCellInfo[nei].count())
2468 {
2469 allFaceInfo[facei].updateFace
2470 (
2471 mesh_,
2472 facei,
2473 nei,
2474 allCellInfo[nei],
2476 dummyTrackData
2477 );
2478 seedFaces.append(facei);
2479 seedFacesInfo.append(allFaceInfo[facei]);
2480 }
2481 }
2482 }
2483
2484 // Seed all boundary faces with owner value. This is to make sure that
2485 // they are visited (probably only important for coupled faces since
2486 // these need to be visited from both sides)
2487 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
2488 {
2489 // Check if face already handled in loop above
2490 if (!allFaceInfo[facei].valid(dummyTrackData))
2491 {
2492 label own = faceOwner[facei];
2493
2494 // Seed face with transported data from owner.
2495 refinementData faceData;
2496 faceData.updateFace
2497 (
2498 mesh_,
2499 facei,
2500 own,
2501 allCellInfo[own],
2503 dummyTrackData
2504 );
2505 seedFaces.append(facei);
2506 seedFacesInfo.append(faceData);
2507 }
2508 }
2509
2510
2511 // face-cell-face transport engine
2513 (
2514 mesh_,
2515 allFaceInfo,
2516 allCellInfo,
2517 dummyTrackData
2518 );
2519
2520 while (true)
2521 {
2522 if (debug)
2523 {
2524 Pout<< "hexRef8::consistentSlowRefinement : Seeded "
2525 << seedFaces.size() << " faces between cells with different"
2526 << " refinement level." << endl;
2527 }
2528
2529 // Set seed faces
2530 levelCalc.setFaceInfo(seedFaces.shrink(), seedFacesInfo.shrink());
2531 seedFaces.clear();
2532 seedFacesInfo.clear();
2533
2534 // Iterate until no change. Now 2:1 face difference should be satisfied
2535 levelCalc.iterate(mesh_.globalData().nTotalFaces()+1);
2536
2537
2538 // Now check point-connected cells (face-connected cells already ok):
2539 // - get per point max of connected cells
2540 // - sync across coupled points
2541 // - check cells against above point max
2542
2543 if (maxPointDiff == -1)
2544 {
2545 // No need to do any point checking.
2546 break;
2547 }
2548
2549 // Determine per point the max cell level. (done as count, not
2550 // as cell level purely for ease)
2551 labelList maxPointCount(mesh_.nPoints(), Zero);
2552
2553 forAll(maxPointCount, pointi)
2554 {
2555 label& pLevel = maxPointCount[pointi];
2556
2557 const labelList& pCells = mesh_.pointCells(pointi);
2558
2559 forAll(pCells, i)
2560 {
2561 pLevel = max(pLevel, allCellInfo[pCells[i]].count());
2562 }
2563 }
2564
2565 // Sync maxPointCount to neighbour
2567 (
2568 mesh_,
2569 maxPointCount,
2571 labelMin // null value
2572 );
2573
2574 // Update allFaceInfo from maxPointCount for all points to check
2575 // (usually on boundary faces)
2576
2577 // Per face the new refinement data
2578 Map<refinementData> changedFacesInfo(pointsToCheck.size());
2579
2580 forAll(pointsToCheck, i)
2581 {
2582 label pointi = pointsToCheck[i];
2583
2584 // Loop over all cells using the point and check whether their
2585 // refinement level is much less than the maximum.
2586
2587 const labelList& pCells = mesh_.pointCells(pointi);
2588
2589 forAll(pCells, pCelli)
2590 {
2591 label celli = pCells[pCelli];
2592
2593 refinementData& cellInfo = allCellInfo[celli];
2594
2595 if
2596 (
2597 !cellInfo.isRefined()
2598 && (
2599 maxPointCount[pointi]
2600 > cellInfo.count() + maxFaceDiff*maxPointDiff
2601 )
2602 )
2603 {
2604 // Mark cell for refinement
2605 cellInfo.count() = cellInfo.refinementCount();
2606
2607 // Insert faces of cell as seed faces.
2608 const cell& cFaces = mesh_.cells()[celli];
2609
2610 forAll(cFaces, cFacei)
2611 {
2612 label facei = cFaces[cFacei];
2613
2614 refinementData faceData;
2615 faceData.updateFace
2616 (
2617 mesh_,
2618 facei,
2619 celli,
2620 cellInfo,
2622 dummyTrackData
2623 );
2624
2625 if (faceData.count() > allFaceInfo[facei].count())
2626 {
2627 changedFacesInfo.insert(facei, faceData);
2628 }
2629 }
2630 }
2631 }
2632 }
2633
2634 label nChanged = changedFacesInfo.size();
2635 reduce(nChanged, sumOp<label>());
2636
2637 if (nChanged == 0)
2638 {
2639 break;
2640 }
2641
2642
2643 // Transfer into seedFaces, seedFacesInfo
2644 seedFaces.setCapacity(changedFacesInfo.size());
2645 seedFacesInfo.setCapacity(changedFacesInfo.size());
2646
2647 forAllConstIters(changedFacesInfo, iter)
2648 {
2649 seedFaces.append(iter.key());
2650 seedFacesInfo.append(iter.val());
2651 }
2652 }
2653
2654
2655 if (debug)
2656 {
2657 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
2658 {
2659 label own = mesh_.faceOwner()[facei];
2660 label ownLevel =
2661 cellLevel_[own]
2662 + (allCellInfo[own].isRefined() ? 1 : 0);
2663
2664 label nei = mesh_.faceNeighbour()[facei];
2665 label neiLevel =
2666 cellLevel_[nei]
2667 + (allCellInfo[nei].isRefined() ? 1 : 0);
2668
2669 if (mag(ownLevel-neiLevel) > 1)
2670 {
2671 dumpCell(own);
2672 dumpCell(nei);
2674 << "cell:" << own
2675 << " current level:" << cellLevel_[own]
2676 << " current refData:" << allCellInfo[own]
2677 << " level after refinement:" << ownLevel
2678 << nl
2679 << "neighbour cell:" << nei
2680 << " current level:" << cellLevel_[nei]
2681 << " current refData:" << allCellInfo[nei]
2682 << " level after refinement:" << neiLevel
2683 << nl
2684 << "which does not satisfy 2:1 constraints anymore." << nl
2685 << "face:" << facei << " faceRefData:" << allFaceInfo[facei]
2686 << abort(FatalError);
2687 }
2688 }
2689
2690
2691 // Coupled faces. Swap owner level to get neighbouring cell level.
2692 // (only boundary faces of neiLevel used)
2693
2694 labelList neiLevel(mesh_.nBoundaryFaces());
2695 labelList neiCount(mesh_.nBoundaryFaces());
2696 labelList neiRefCount(mesh_.nBoundaryFaces());
2697
2698 forAll(neiLevel, i)
2699 {
2700 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
2701 neiLevel[i] = cellLevel_[own];
2702 neiCount[i] = allCellInfo[own].count();
2703 neiRefCount[i] = allCellInfo[own].refinementCount();
2704 }
2705
2706 // Swap to neighbour
2707 syncTools::swapBoundaryFaceList(mesh_, neiLevel);
2708 syncTools::swapBoundaryFaceList(mesh_, neiCount);
2709 syncTools::swapBoundaryFaceList(mesh_, neiRefCount);
2710
2711 // Now we have neighbour value see which cells need refinement
2712 forAll(neiLevel, i)
2713 {
2714 label facei = i+mesh_.nInternalFaces();
2715
2716 label own = mesh_.faceOwner()[facei];
2717 label ownLevel =
2718 cellLevel_[own]
2719 + (allCellInfo[own].isRefined() ? 1 : 0);
2720
2721 label nbrLevel =
2722 neiLevel[i]
2723 + ((neiCount[i] >= neiRefCount[i]) ? 1 : 0);
2724
2725 if (mag(ownLevel - nbrLevel) > 1)
2726 {
2727 dumpCell(own);
2728 label patchi = mesh_.boundaryMesh().whichPatch(facei);
2729
2731 << "Celllevel does not satisfy 2:1 constraint."
2732 << " On coupled face "
2733 << facei
2734 << " refData:" << allFaceInfo[facei]
2735 << " on patch " << patchi << " "
2736 << mesh_.boundaryMesh()[patchi].name() << nl
2737 << "owner cell " << own
2738 << " current level:" << cellLevel_[own]
2739 << " current count:" << allCellInfo[own].count()
2740 << " current refCount:"
2741 << allCellInfo[own].refinementCount()
2742 << " level after refinement:" << ownLevel
2743 << nl
2744 << "(coupled) neighbour cell"
2745 << " has current level:" << neiLevel[i]
2746 << " current count:" << neiCount[i]
2747 << " current refCount:" << neiRefCount[i]
2748 << " level after refinement:" << nbrLevel
2749 << abort(FatalError);
2750 }
2751 }
2752 }
2753
2754 // Convert back to labelList of cells to refine.
2755
2756 label nRefined = 0;
2757
2758 forAll(allCellInfo, celli)
2759 {
2760 if (allCellInfo[celli].isRefined())
2761 {
2762 nRefined++;
2763 }
2764 }
2765
2766 // Updated list of cells to refine
2767 labelList newCellsToRefine(nRefined);
2768 nRefined = 0;
2769
2770 forAll(allCellInfo, celli)
2771 {
2772 if (allCellInfo[celli].isRefined())
2773 {
2774 newCellsToRefine[nRefined++] = celli;
2775 }
2776 }
2777
2778 if (debug)
2779 {
2780 Pout<< "hexRef8::consistentSlowRefinement : From "
2781 << cellsToRefine.size() << " to " << newCellsToRefine.size()
2782 << " cells to refine." << endl;
2783 }
2784
2785 return newCellsToRefine;
2786}
2787
2788
2790(
2791 const label maxFaceDiff,
2792 const labelList& cellsToRefine,
2793 const labelList& facesToCheck
2794) const
2795{
2796 const labelList& faceOwner = mesh_.faceOwner();
2797 const labelList& faceNeighbour = mesh_.faceNeighbour();
2798
2799 if (maxFaceDiff <= 0)
2800 {
2802 << "Illegal maxFaceDiff " << maxFaceDiff << nl
2803 << "Value should be >= 1" << exit(FatalError);
2804 }
2805
2806 const scalar level0Size = 2*maxFaceDiff*level0EdgeLength();
2807
2808
2809 // Bit tricky. Say we want a distance of three cells between two
2810 // consecutive refinement levels. This is done by using FaceCellWave to
2811 // transport out the 'refinement shell'. Anything inside the refinement
2812 // shell (given by a distance) gets marked for refinement.
2813
2814 // Initial information about (distance to) cellLevel on all cells
2815 List<refinementDistanceData> allCellInfo(mesh_.nCells());
2816
2817 // Initial information about (distance to) cellLevel on all faces
2818 List<refinementDistanceData> allFaceInfo(mesh_.nFaces());
2819
2820 // Dummy additional info for FaceCellWave
2821 int dummyTrackData = 0;
2822
2823
2824 // Mark cells with wanted refinement level
2825 forAll(cellsToRefine, i)
2826 {
2827 label celli = cellsToRefine[i];
2828
2829 allCellInfo[celli] = refinementDistanceData
2830 (
2831 level0Size,
2832 mesh_.cellCentres()[celli],
2833 cellLevel_[celli]+1 // wanted refinement
2834 );
2835 }
2836 // Mark all others with existing refinement level
2837 forAll(allCellInfo, celli)
2838 {
2839 if (!allCellInfo[celli].valid(dummyTrackData))
2840 {
2841 allCellInfo[celli] = refinementDistanceData
2842 (
2843 level0Size,
2844 mesh_.cellCentres()[celli],
2845 cellLevel_[celli] // wanted refinement
2846 );
2847 }
2848 }
2849
2850
2851 // Labels of seed faces
2852 DynamicList<label> seedFaces(mesh_.nFaces()/100);
2853 // refinementLevel data on seed faces
2854 DynamicList<refinementDistanceData> seedFacesInfo(mesh_.nFaces()/100);
2855
2856 const pointField& cc = mesh_.cellCentres();
2857
2858 forAll(facesToCheck, i)
2859 {
2860 label facei = facesToCheck[i];
2861
2862 if (allFaceInfo[facei].valid(dummyTrackData))
2863 {
2864 // Can only occur if face has already gone through loop below.
2866 << "Argument facesToCheck seems to have duplicate entries!"
2867 << endl
2868 << "face:" << facei << " occurs at positions "
2869 << findIndices(facesToCheck, facei)
2870 << abort(FatalError);
2871 }
2872
2873 label own = faceOwner[facei];
2874
2875 label ownLevel =
2876 (
2877 allCellInfo[own].valid(dummyTrackData)
2878 ? allCellInfo[own].originLevel()
2879 : cellLevel_[own]
2880 );
2881
2882 if (!mesh_.isInternalFace(facei))
2883 {
2884 // Do as if boundary face would have neighbour with one higher
2885 // refinement level.
2886 const point& fc = mesh_.faceCentres()[facei];
2887
2889 (
2890 level0Size,
2891 2*fc - cc[own], // est'd cell centre
2892 ownLevel+1
2893 );
2894
2895 allFaceInfo[facei].updateFace
2896 (
2897 mesh_,
2898 facei,
2899 own, // not used (should be nei)
2900 neiData,
2902 dummyTrackData
2903 );
2904 }
2905 else
2906 {
2907 label nei = faceNeighbour[facei];
2908
2909 label neiLevel =
2910 (
2911 allCellInfo[nei].valid(dummyTrackData)
2912 ? allCellInfo[nei].originLevel()
2913 : cellLevel_[nei]
2914 );
2915
2916 if (ownLevel == neiLevel)
2917 {
2918 // Fake as if nei>own or own>nei (whichever one 'wins')
2919 allFaceInfo[facei].updateFace
2920 (
2921 mesh_,
2922 facei,
2923 nei,
2924 refinementDistanceData(level0Size, cc[nei], neiLevel+1),
2926 dummyTrackData
2927 );
2928 allFaceInfo[facei].updateFace
2929 (
2930 mesh_,
2931 facei,
2932 own,
2933 refinementDistanceData(level0Size, cc[own], ownLevel+1),
2935 dummyTrackData
2936 );
2937 }
2938 else
2939 {
2940 // Difference in level anyway.
2941 allFaceInfo[facei].updateFace
2942 (
2943 mesh_,
2944 facei,
2945 nei,
2946 refinementDistanceData(level0Size, cc[nei], neiLevel),
2948 dummyTrackData
2949 );
2950 allFaceInfo[facei].updateFace
2951 (
2952 mesh_,
2953 facei,
2954 own,
2955 refinementDistanceData(level0Size, cc[own], ownLevel),
2957 dummyTrackData
2958 );
2959 }
2960 }
2961 seedFaces.append(facei);
2962 seedFacesInfo.append(allFaceInfo[facei]);
2963 }
2964
2965
2966 // Create some initial seeds to start walking from. This is only if there
2967 // are no facesToCheck.
2968 // Just seed with all faces inbetween different refinement levels for now
2969 forAll(faceNeighbour, facei)
2970 {
2971 // Check if face already handled in loop above
2972 if (!allFaceInfo[facei].valid(dummyTrackData))
2973 {
2974 label own = faceOwner[facei];
2975
2976 label ownLevel =
2977 (
2978 allCellInfo[own].valid(dummyTrackData)
2979 ? allCellInfo[own].originLevel()
2980 : cellLevel_[own]
2981 );
2982
2983 label nei = faceNeighbour[facei];
2984
2985 label neiLevel =
2986 (
2987 allCellInfo[nei].valid(dummyTrackData)
2988 ? allCellInfo[nei].originLevel()
2989 : cellLevel_[nei]
2990 );
2991
2992 if (ownLevel > neiLevel)
2993 {
2994 // Set face to owner data. (since face not yet would be copy)
2995 seedFaces.append(facei);
2996 allFaceInfo[facei].updateFace
2997 (
2998 mesh_,
2999 facei,
3000 own,
3001 refinementDistanceData(level0Size, cc[own], ownLevel),
3003 dummyTrackData
3004 );
3005 seedFacesInfo.append(allFaceInfo[facei]);
3006 }
3007 else if (neiLevel > ownLevel)
3008 {
3009 seedFaces.append(facei);
3010 allFaceInfo[facei].updateFace
3011 (
3012 mesh_,
3013 facei,
3014 nei,
3015 refinementDistanceData(level0Size, cc[nei], neiLevel),
3017 dummyTrackData
3018 );
3019 seedFacesInfo.append(allFaceInfo[facei]);
3020 }
3021 }
3022 }
3023
3024 seedFaces.shrink();
3025 seedFacesInfo.shrink();
3026
3027 // face-cell-face transport engine
3029 (
3030 mesh_,
3031 seedFaces,
3032 seedFacesInfo,
3033 allFaceInfo,
3034 allCellInfo,
3035 mesh_.globalData().nTotalCells()+1,
3036 dummyTrackData
3037 );
3038
3039
3040 //if (debug)
3041 //{
3042 // // Dump wanted level
3043 // volScalarField wantedLevel
3044 // (
3045 // IOobject
3046 // (
3047 // "wantedLevel",
3048 // fMesh.time().timeName(),
3049 // fMesh,
3050 // IOobject::NO_READ,
3051 // IOobject::NO_WRITE,
3052 // false
3053 // ),
3054 // fMesh,
3055 // dimensionedScalar(dimless, Zero)
3056 // );
3057 //
3058 // forAll(wantedLevel, celli)
3059 // {
3060 // wantedLevel[celli] = allCellInfo[celli].wantedLevel(cc[celli]);
3061 // }
3062 //
3063 // Pout<< "Writing " << wantedLevel.objectPath() << endl;
3064 // wantedLevel.write();
3065 //}
3066
3067
3068 // Convert back to labelList of cells to refine.
3069 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3070
3071 // 1. Force original refinement cells to be picked up by setting the
3072 // originLevel of input cells to be a very large level (but within range
3073 // of 1<< shift inside refinementDistanceData::wantedLevel)
3074 forAll(cellsToRefine, i)
3075 {
3076 label celli = cellsToRefine[i];
3077
3078 allCellInfo[celli].originLevel() = sizeof(label)*8-2;
3079 allCellInfo[celli].origin() = cc[celli];
3080 }
3081
3082 // 2. Extend to 2:1. I don't understand yet why this is not done
3083 // 2. Extend to 2:1. For non-cube cells the scalar distance does not work
3084 // so make sure it at least provides 2:1.
3085 bitSet refineCell(mesh_.nCells());
3086 forAll(allCellInfo, celli)
3087 {
3088 label wanted = allCellInfo[celli].wantedLevel(cc[celli]);
3089
3090 if (wanted > cellLevel_[celli]+1)
3091 {
3092 refineCell.set(celli);
3093 }
3094 }
3095 faceConsistentRefinement(true, cellLevel_, refineCell);
3096
3097 while (true)
3098 {
3099 label nChanged = faceConsistentRefinement(true, cellLevel_, refineCell);
3100
3101 reduce(nChanged, sumOp<label>());
3102
3103 if (debug)
3104 {
3105 Pout<< "hexRef8::consistentSlowRefinement2 : Changed " << nChanged
3106 << " refinement levels due to 2:1 conflicts."
3107 << endl;
3108 }
3109
3110 if (nChanged == 0)
3111 {
3112 break;
3113 }
3114 }
3115
3116 // 3. Convert back to labelList.
3117 labelList newCellsToRefine(refineCell.toc());
3118
3119 if (debug)
3120 {
3121 Pout<< "hexRef8::consistentSlowRefinement2 : From "
3122 << cellsToRefine.size() << " to " << newCellsToRefine.size()
3123 << " cells to refine." << endl;
3124
3125 // Check that newCellsToRefine obeys at least 2:1.
3126
3127 {
3128 cellSet cellsIn(mesh_, "cellsToRefineIn", cellsToRefine);
3129 Pout<< "hexRef8::consistentSlowRefinement2 : writing "
3130 << cellsIn.size() << " to cellSet "
3131 << cellsIn.objectPath() << endl;
3132 cellsIn.write();
3133 }
3134 {
3135 cellSet cellsOut(mesh_, "cellsToRefineOut", newCellsToRefine);
3136 Pout<< "hexRef8::consistentSlowRefinement2 : writing "
3137 << cellsOut.size() << " to cellSet "
3138 << cellsOut.objectPath() << endl;
3139 cellsOut.write();
3140 }
3141
3142 // Extend to 2:1
3143 bitSet refineCell(mesh_.nCells(), newCellsToRefine);
3144
3145 const bitSet savedRefineCell(refineCell);
3146 label nChanged = faceConsistentRefinement(true, cellLevel_, refineCell);
3147
3148 {
3149 cellSet cellsOut2
3150 (
3151 mesh_, "cellsToRefineOut2", newCellsToRefine.size()
3152 );
3153 forAll(refineCell, celli)
3154 {
3155 if (refineCell.test(celli))
3156 {
3157 cellsOut2.insert(celli);
3158 }
3159 }
3160 Pout<< "hexRef8::consistentSlowRefinement2 : writing "
3161 << cellsOut2.size() << " to cellSet "
3162 << cellsOut2.objectPath() << endl;
3163 cellsOut2.write();
3164 }
3165
3166 if (nChanged > 0)
3167 {
3168 forAll(refineCell, celli)
3169 {
3170 if (refineCell.test(celli) && !savedRefineCell.test(celli))
3171 {
3172 dumpCell(celli);
3174 << "Cell:" << celli << " cc:"
3175 << mesh_.cellCentres()[celli]
3176 << " was not marked for refinement but does not obey"
3177 << " 2:1 constraints."
3178 << abort(FatalError);
3179 }
3180 }
3181 }
3182 }
3183
3184 return newCellsToRefine;
3185}
3186
3187
3188// Top level driver to insert topo changes to do all refinement.
3190(
3191 const labelList& cellLabels,
3192 polyTopoChange& meshMod
3193)
3194{
3195 if (debug)
3196 {
3197 Pout<< "hexRef8::setRefinement :"
3198 << " Checking initial mesh just to make sure" << endl;
3199
3200 checkMesh();
3201 // Cannot call checkRefinementlevels since hanging points might
3202 // get triggered by the mesher after subsetting.
3203 //checkRefinementLevels(-1, labelList(0));
3204 }
3205
3206 // Clear any saved point/cell data.
3207 savedPointLevel_.clear();
3208 savedCellLevel_.clear();
3209
3210
3211 // New point/cell level. Copy of pointLevel for existing points.
3212 DynamicList<label> newCellLevel(cellLevel_.size());
3213 forAll(cellLevel_, celli)
3214 {
3215 newCellLevel.append(cellLevel_[celli]);
3216 }
3217 DynamicList<label> newPointLevel(pointLevel_.size());
3218 forAll(pointLevel_, pointi)
3219 {
3220 newPointLevel.append(pointLevel_[pointi]);
3221 }
3222
3223
3224 if (debug)
3225 {
3226 Pout<< "hexRef8::setRefinement :"
3227 << " Allocating " << cellLabels.size() << " cell midpoints."
3228 << endl;
3229 }
3230
3231
3232 // Mid point per refined cell.
3233 // -1 : not refined
3234 // >=0: label of mid point.
3235 labelList cellMidPoint(mesh_.nCells(), -1);
3236
3237 forAll(cellLabels, i)
3238 {
3239 label celli = cellLabels[i];
3240
3241 label anchorPointi = mesh_.faces()[mesh_.cells()[celli][0]][0];
3242
3243 cellMidPoint[celli] = meshMod.setAction
3244 (
3246 (
3247 mesh_.cellCentres()[celli], // point
3248 anchorPointi, // master point
3249 -1, // zone for point
3250 true // supports a cell
3251 )
3252 );
3253
3254 newPointLevel(cellMidPoint[celli]) = cellLevel_[celli]+1;
3255 }
3256
3257
3258 if (debug)
3259 {
3260 cellSet splitCells(mesh_, "splitCells", cellLabels.size());
3261
3262 forAll(cellMidPoint, celli)
3263 {
3264 if (cellMidPoint[celli] >= 0)
3265 {
3266 splitCells.insert(celli);
3267 }
3268 }
3269
3270 Pout<< "hexRef8::setRefinement : Dumping " << splitCells.size()
3271 << " cells to split to cellSet " << splitCells.objectPath()
3272 << endl;
3273
3274 splitCells.write();
3275 }
3276
3277
3278
3279 // Split edges
3280 // ~~~~~~~~~~~
3281
3282 if (debug)
3283 {
3284 Pout<< "hexRef8::setRefinement :"
3285 << " Allocating edge midpoints."
3286 << endl;
3287 }
3288
3289 // Unrefined edges are ones between cellLevel or lower points.
3290 // If any cell using this edge gets split then the edge needs to be split.
3291
3292 // -1 : no need to split edge
3293 // >=0 : label of introduced mid point
3294 labelList edgeMidPoint(mesh_.nEdges(), -1);
3295
3296 // Note: Loop over cells to be refined or edges?
3297 forAll(cellMidPoint, celli)
3298 {
3299 if (cellMidPoint[celli] >= 0)
3300 {
3301 const labelList& cEdges = mesh_.cellEdges(celli);
3302
3303 forAll(cEdges, i)
3304 {
3305 label edgeI = cEdges[i];
3306
3307 const edge& e = mesh_.edges()[edgeI];
3308
3309 if
3310 (
3311 pointLevel_[e[0]] <= cellLevel_[celli]
3312 && pointLevel_[e[1]] <= cellLevel_[celli]
3313 )
3314 {
3315 edgeMidPoint[edgeI] = 12345; // mark need for splitting
3316 }
3317 }
3318 }
3319 }
3320
3321 // Synchronize edgeMidPoint across coupled patches. Take max so that
3322 // any split takes precedence.
3324 (
3325 mesh_,
3326 edgeMidPoint,
3328 labelMin
3329 );
3330
3331
3332 // Introduce edge points
3333 // ~~~~~~~~~~~~~~~~~~~~~
3334
3335 {
3336 // Phase 1: calculate midpoints and sync.
3337 // This needs doing for if people do not write binary and we slowly
3338 // get differences.
3339
3340 pointField edgeMids(mesh_.nEdges(), point(-GREAT, -GREAT, -GREAT));
3341
3342 forAll(edgeMidPoint, edgeI)
3343 {
3344 if (edgeMidPoint[edgeI] >= 0)
3345 {
3346 // Edge marked to be split.
3347 edgeMids[edgeI] = mesh_.edges()[edgeI].centre(mesh_.points());
3348 }
3349 }
3351 (
3352 mesh_,
3353 edgeMids,
3355 point(-GREAT, -GREAT, -GREAT)
3356 );
3357
3358
3359 // Phase 2: introduce points at the synced locations.
3360 forAll(edgeMidPoint, edgeI)
3361 {
3362 if (edgeMidPoint[edgeI] >= 0)
3363 {
3364 // Edge marked to be split. Replace edgeMidPoint with actual
3365 // point label.
3366
3367 const edge& e = mesh_.edges()[edgeI];
3368
3369 edgeMidPoint[edgeI] = meshMod.setAction
3370 (
3372 (
3373 edgeMids[edgeI], // point
3374 e[0], // master point
3375 -1, // zone for point
3376 true // supports a cell
3377 )
3378 );
3379
3380 newPointLevel(edgeMidPoint[edgeI]) =
3381 max
3382 (
3383 pointLevel_[e[0]],
3384 pointLevel_[e[1]]
3385 )
3386 + 1;
3387 }
3388 }
3389 }
3390
3391 if (debug)
3392 {
3393 OFstream str(mesh_.time().path()/"edgeMidPoint.obj");
3394
3395 forAll(edgeMidPoint, edgeI)
3396 {
3397 if (edgeMidPoint[edgeI] >= 0)
3398 {
3399 const edge& e = mesh_.edges()[edgeI];
3400
3401 meshTools::writeOBJ(str, e.centre(mesh_.points()));
3402 }
3403 }
3404
3405 Pout<< "hexRef8::setRefinement :"
3406 << " Dumping edge centres to split to file " << str.name() << endl;
3407 }
3408
3409
3410 // Calculate face level
3411 // ~~~~~~~~~~~~~~~~~~~~
3412 // (after splitting)
3413
3414 if (debug)
3415 {
3416 Pout<< "hexRef8::setRefinement :"
3417 << " Allocating face midpoints."
3418 << endl;
3419 }
3420
3421 // Face anchor level. There are guaranteed 4 points with level
3422 // <= anchorLevel. These are the corner points.
3423 labelList faceAnchorLevel(mesh_.nFaces());
3424
3425 for (label facei = 0; facei < mesh_.nFaces(); facei++)
3426 {
3427 faceAnchorLevel[facei] = faceLevel(facei);
3428 }
3429
3430 // -1 : no need to split face
3431 // >=0 : label of introduced mid point
3432 labelList faceMidPoint(mesh_.nFaces(), -1);
3433
3434
3435 // Internal faces: look at cells on both sides. Uniquely determined since
3436 // face itself guaranteed to be same level as most refined neighbour.
3437 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
3438 {
3439 if (faceAnchorLevel[facei] >= 0)
3440 {
3441 label own = mesh_.faceOwner()[facei];
3442 label ownLevel = cellLevel_[own];
3443 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3444
3445 label nei = mesh_.faceNeighbour()[facei];
3446 label neiLevel = cellLevel_[nei];
3447 label newNeiLevel = neiLevel + (cellMidPoint[nei] >= 0 ? 1 : 0);
3448
3449 if
3450 (
3451 newOwnLevel > faceAnchorLevel[facei]
3452 || newNeiLevel > faceAnchorLevel[facei]
3453 )
3454 {
3455 faceMidPoint[facei] = 12345; // mark to be split
3456 }
3457 }
3458 }
3459
3460 // Coupled patches handled like internal faces except now all information
3461 // from neighbour comes from across processor.
3462 // Boundary faces are more complicated since the boundary face can
3463 // be more refined than its owner (or neighbour for coupled patches)
3464 // (does not happen if refining/unrefining only, but does e.g. when
3465 // refinining and subsetting)
3466
3467 {
3468 labelList newNeiLevel(mesh_.nBoundaryFaces());
3469
3470 forAll(newNeiLevel, i)
3471 {
3472 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
3473 label ownLevel = cellLevel_[own];
3474 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3475
3476 newNeiLevel[i] = newOwnLevel;
3477 }
3478
3479 // Swap.
3480 syncTools::swapBoundaryFaceList(mesh_, newNeiLevel);
3481
3482 // So now we have information on the neighbour.
3483
3484 forAll(newNeiLevel, i)
3485 {
3486 label facei = i+mesh_.nInternalFaces();
3487
3488 if (faceAnchorLevel[facei] >= 0)
3489 {
3490 label own = mesh_.faceOwner()[facei];
3491 label ownLevel = cellLevel_[own];
3492 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3493
3494 if
3495 (
3496 newOwnLevel > faceAnchorLevel[facei]
3497 || newNeiLevel[i] > faceAnchorLevel[facei]
3498 )
3499 {
3500 faceMidPoint[facei] = 12345; // mark to be split
3501 }
3502 }
3503 }
3504 }
3505
3506
3507 // Synchronize faceMidPoint across coupled patches. (logical or)
3509 (
3510 mesh_,
3511 faceMidPoint,
3513 );
3514
3515
3516
3517 // Introduce face points
3518 // ~~~~~~~~~~~~~~~~~~~~~
3519
3520 {
3521 // Phase 1: determine mid points and sync. See comment for edgeMids
3522 // above
3523 pointField bFaceMids
3524 (
3525 mesh_.nBoundaryFaces(),
3526 point(-GREAT, -GREAT, -GREAT)
3527 );
3528
3529 forAll(bFaceMids, i)
3530 {
3531 label facei = i+mesh_.nInternalFaces();
3532
3533 if (faceMidPoint[facei] >= 0)
3534 {
3535 bFaceMids[i] = mesh_.faceCentres()[facei];
3536 }
3537 }
3539 (
3540 mesh_,
3541 bFaceMids,
3543 );
3544
3545 forAll(faceMidPoint, facei)
3546 {
3547 if (faceMidPoint[facei] >= 0)
3548 {
3549 // Face marked to be split. Replace faceMidPoint with actual
3550 // point label.
3551
3552 const face& f = mesh_.faces()[facei];
3553
3554 faceMidPoint[facei] = meshMod.setAction
3555 (
3557 (
3558 (
3559 facei < mesh_.nInternalFaces()
3560 ? mesh_.faceCentres()[facei]
3561 : bFaceMids[facei-mesh_.nInternalFaces()]
3562 ), // point
3563 f[0], // master point
3564 -1, // zone for point
3565 true // supports a cell
3566 )
3567 );
3568
3569 // Determine the level of the corner points and midpoint will
3570 // be one higher.
3571 newPointLevel(faceMidPoint[facei]) = faceAnchorLevel[facei]+1;
3572 }
3573 }
3574 }
3575
3576 if (debug)
3577 {
3578 faceSet splitFaces(mesh_, "splitFaces", cellLabels.size());
3579
3580 forAll(faceMidPoint, facei)
3581 {
3582 if (faceMidPoint[facei] >= 0)
3583 {
3584 splitFaces.insert(facei);
3585 }
3586 }
3587
3588 Pout<< "hexRef8::setRefinement : Dumping " << splitFaces.size()
3589 << " faces to split to faceSet " << splitFaces.objectPath() << endl;
3590
3591 splitFaces.write();
3592 }
3593
3594
3595 // Information complete
3596 // ~~~~~~~~~~~~~~~~~~~~
3597 // At this point we have all the information we need. We should no
3598 // longer reference the cellLabels to refine. All the information is:
3599 // - cellMidPoint >= 0 : cell needs to be split
3600 // - faceMidPoint >= 0 : face needs to be split
3601 // - edgeMidPoint >= 0 : edge needs to be split
3602
3603
3604
3605 // Get the corner/anchor points
3606 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3607
3608 if (debug)
3609 {
3610 Pout<< "hexRef8::setRefinement :"
3611 << " Finding cell anchorPoints (8 per cell)"
3612 << endl;
3613 }
3614
3615 // There will always be 8 points on the hex that have were introduced
3616 // with the hex and will have the same or lower refinement level.
3617
3618 // Per cell the 8 corner points.
3619 labelListList cellAnchorPoints(mesh_.nCells());
3620
3621 {
3622 labelList nAnchorPoints(mesh_.nCells(), Zero);
3623
3624 forAll(cellMidPoint, celli)
3625 {
3626 if (cellMidPoint[celli] >= 0)
3627 {
3628 cellAnchorPoints[celli].setSize(8);
3629 }
3630 }
3631
3632 forAll(pointLevel_, pointi)
3633 {
3634 const labelList& pCells = mesh_.pointCells(pointi);
3635
3636 forAll(pCells, pCelli)
3637 {
3638 label celli = pCells[pCelli];
3639
3640 if
3641 (
3642 cellMidPoint[celli] >= 0
3643 && pointLevel_[pointi] <= cellLevel_[celli]
3644 )
3645 {
3646 if (nAnchorPoints[celli] == 8)
3647 {
3648 dumpCell(celli);
3650 << "cell " << celli
3651 << " of level " << cellLevel_[celli]
3652 << " uses more than 8 points of equal or"
3653 << " lower level" << nl
3654 << "Points so far:" << cellAnchorPoints[celli]
3655 << abort(FatalError);
3656 }
3657
3658 cellAnchorPoints[celli][nAnchorPoints[celli]++] = pointi;
3659 }
3660 }
3661 }
3662
3663 forAll(cellMidPoint, celli)
3664 {
3665 if (cellMidPoint[celli] >= 0)
3666 {
3667 if (nAnchorPoints[celli] != 8)
3668 {
3669 dumpCell(celli);
3670
3671 const labelList& cPoints = mesh_.cellPoints(celli);
3672
3674 << "cell " << celli
3675 << " of level " << cellLevel_[celli]
3676 << " does not seem to have 8 points of equal or"
3677 << " lower level" << endl
3678 << "cellPoints:" << cPoints << endl
3679 << "pointLevels:"
3680 << labelUIndList(pointLevel_, cPoints) << endl
3681 << abort(FatalError);
3682 }
3683 }
3684 }
3685 }
3686
3687
3688 // Add the cells
3689 // ~~~~~~~~~~~~~
3690
3691 if (debug)
3692 {
3693 Pout<< "hexRef8::setRefinement :"
3694 << " Adding cells (1 per anchorPoint)"
3695 << endl;
3696 }
3697
3698 // Per cell the 7 added cells (+ original cell)
3699 labelListList cellAddedCells(mesh_.nCells());
3700
3701 forAll(cellAnchorPoints, celli)
3702 {
3703 const labelList& cAnchors = cellAnchorPoints[celli];
3704
3705 if (cAnchors.size() == 8)
3706 {
3707 labelList& cAdded = cellAddedCells[celli];
3708 cAdded.setSize(8);
3709
3710 // Original cell at 0
3711 cAdded[0] = celli;
3712
3713 // Update cell level
3714 newCellLevel[celli] = cellLevel_[celli]+1;
3715
3716
3717 for (label i = 1; i < 8; i++)
3718 {
3719 cAdded[i] = meshMod.setAction
3720 (
3722 (
3723 -1, // master point
3724 -1, // master edge
3725 -1, // master face
3726 celli, // master cell
3727 mesh_.cellZones().whichZone(celli) // zone for cell
3728 )
3729 );
3730
3731 newCellLevel(cAdded[i]) = cellLevel_[celli]+1;
3732 }
3733 }
3734 }
3735
3736
3737 // Faces
3738 // ~~~~~
3739 // 1. existing faces that get split (into four always)
3740 // 2. existing faces that do not get split but only edges get split
3741 // 3. existing faces that do not get split but get new owner/neighbour
3742 // 4. new internal faces inside split cells.
3743
3744 if (debug)
3745 {
3746 Pout<< "hexRef8::setRefinement :"
3747 << " Marking faces to be handled"
3748 << endl;
3749 }
3750
3751 // Get all affected faces.
3752 bitSet affectedFace(mesh_.nFaces());
3753
3754 {
3755 forAll(cellMidPoint, celli)
3756 {
3757 if (cellMidPoint[celli] >= 0)
3758 {
3759 const cell& cFaces = mesh_.cells()[celli];
3760
3761 affectedFace.set(cFaces);
3762 }
3763 }
3764
3765 forAll(faceMidPoint, facei)
3766 {
3767 if (faceMidPoint[facei] >= 0)
3768 {
3769 affectedFace.set(facei);
3770 }
3771 }
3772
3773 forAll(edgeMidPoint, edgeI)
3774 {
3775 if (edgeMidPoint[edgeI] >= 0)
3776 {
3777 const labelList& eFaces = mesh_.edgeFaces(edgeI);
3778
3779 affectedFace.set(eFaces);
3780 }
3781 }
3782 }
3783
3784
3785 // 1. Faces that get split
3786 // ~~~~~~~~~~~~~~~~~~~~~~~
3787
3788 if (debug)
3789 {
3790 Pout<< "hexRef8::setRefinement : Splitting faces" << endl;
3791 }
3792
3793 forAll(faceMidPoint, facei)
3794 {
3795 if (faceMidPoint[facei] >= 0 && affectedFace.test(facei))
3796 {
3797 // Face needs to be split and hasn't yet been done in some way
3798 // (affectedFace - is impossible since this is first change but
3799 // just for completeness)
3800
3801 const face& f = mesh_.faces()[facei];
3802
3803 // Has original facei been used (three faces added, original gets
3804 // modified)
3805 bool modifiedFace = false;
3806
3807 label anchorLevel = faceAnchorLevel[facei];
3808
3809 face newFace(4);
3810
3811 forAll(f, fp)
3812 {
3813 label pointi = f[fp];
3814
3815 if (pointLevel_[pointi] <= anchorLevel)
3816 {
3817 // point is anchor. Start collecting face.
3818
3819 DynamicList<label> faceVerts(4);
3820
3821 faceVerts.append(pointi);
3822
3823 // Walk forward to mid point.
3824 // - if next is +2 midpoint is +1
3825 // - if next is +1 it is midpoint
3826 // - if next is +0 there has to be edgeMidPoint
3827
3828 walkFaceToMid
3829 (
3830 edgeMidPoint,
3831 anchorLevel,
3832 facei,
3833 fp,
3834 faceVerts
3835 );
3836
3837 faceVerts.append(faceMidPoint[facei]);
3838
3839 walkFaceFromMid
3840 (
3841 edgeMidPoint,
3842 anchorLevel,
3843 facei,
3844 fp,
3845 faceVerts
3846 );
3847
3848 // Convert dynamiclist to face.
3849 newFace.transfer(faceVerts);
3850
3851 //Pout<< "Split face:" << facei << " verts:" << f
3852 // << " into quad:" << newFace << endl;
3853
3854 // Get new owner/neighbour
3855 label own, nei;
3856 getFaceNeighbours
3857 (
3858 cellAnchorPoints,
3859 cellAddedCells,
3860 facei,
3861 pointi, // Anchor point
3862
3863 own,
3864 nei
3865 );
3866
3867
3868 if (debug)
3869 {
3870 if (mesh_.isInternalFace(facei))
3871 {
3872 label oldOwn = mesh_.faceOwner()[facei];
3873 label oldNei = mesh_.faceNeighbour()[facei];
3874
3875 checkInternalOrientation
3876 (
3877 meshMod,
3878 oldOwn,
3879 facei,
3880 mesh_.cellCentres()[oldOwn],
3881 mesh_.cellCentres()[oldNei],
3882 newFace
3883 );
3884 }
3885 else
3886 {
3887 label oldOwn = mesh_.faceOwner()[facei];
3888
3889 checkBoundaryOrientation
3890 (
3891 meshMod,
3892 oldOwn,
3893 facei,
3894 mesh_.cellCentres()[oldOwn],
3895 mesh_.faceCentres()[facei],
3896 newFace
3897 );
3898 }
3899 }
3900
3901
3902 if (!modifiedFace)
3903 {
3904 modifiedFace = true;
3905
3906 modFace(meshMod, facei, newFace, own, nei);
3907 }
3908 else
3909 {
3910 addFace(meshMod, facei, newFace, own, nei);
3911 }
3912 }
3913 }
3914
3915 // Mark face as having been handled
3916 affectedFace.unset(facei);
3917 }
3918 }
3919
3920
3921 // 2. faces that do not get split but use edges that get split
3922 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3923
3924 if (debug)
3925 {
3926 Pout<< "hexRef8::setRefinement :"
3927 << " Adding edge splits to unsplit faces"
3928 << endl;
3929 }
3930
3931 DynamicList<label> eFacesStorage;
3932 DynamicList<label> fEdgesStorage;
3933
3934 forAll(edgeMidPoint, edgeI)
3935 {
3936 if (edgeMidPoint[edgeI] >= 0)
3937 {
3938 // Split edge. Check that face not already handled above.
3939
3940 const labelList& eFaces = mesh_.edgeFaces(edgeI, eFacesStorage);
3941
3942 forAll(eFaces, i)
3943 {
3944 label facei = eFaces[i];
3945
3946 if (faceMidPoint[facei] < 0 && affectedFace.test(facei))
3947 {
3948 // Unsplit face. Add edge splits to face.
3949
3950 const face& f = mesh_.faces()[facei];
3951 const labelList& fEdges = mesh_.faceEdges
3952 (
3953 facei,
3954 fEdgesStorage
3955 );
3956
3957 DynamicList<label> newFaceVerts(f.size());
3958
3959 forAll(f, fp)
3960 {
3961 newFaceVerts.append(f[fp]);
3962
3963 label edgeI = fEdges[fp];
3964
3965 if (edgeMidPoint[edgeI] >= 0)
3966 {
3967 newFaceVerts.append(edgeMidPoint[edgeI]);
3968 }
3969 }
3970
3971 face newFace;
3972 newFace.transfer(newFaceVerts);
3973
3974 // The point with the lowest level should be an anchor
3975 // point of the neighbouring cells.
3976 label anchorFp = findMinLevel(f);
3977
3978 label own, nei;
3979 getFaceNeighbours
3980 (
3981 cellAnchorPoints,
3982 cellAddedCells,
3983 facei,
3984 f[anchorFp], // Anchor point
3985
3986 own,
3987 nei
3988 );
3989
3990
3991 if (debug)
3992 {
3993 if (mesh_.isInternalFace(facei))
3994 {
3995 label oldOwn = mesh_.faceOwner()[facei];
3996 label oldNei = mesh_.faceNeighbour()[facei];
3997
3998 checkInternalOrientation
3999 (
4000 meshMod,
4001 oldOwn,
4002 facei,
4003 mesh_.cellCentres()[oldOwn],
4004 mesh_.cellCentres()[oldNei],
4005 newFace
4006 );
4007 }
4008 else
4009 {
4010 label oldOwn = mesh_.faceOwner()[facei];
4011
4012 checkBoundaryOrientation
4013 (
4014 meshMod,
4015 oldOwn,
4016 facei,
4017 mesh_.cellCentres()[oldOwn],
4018 mesh_.faceCentres()[facei],
4019 newFace
4020 );
4021 }
4022 }
4023
4024 modFace(meshMod, facei, newFace, own, nei);
4025
4026 // Mark face as having been handled
4027 affectedFace.unset(facei);
4028 }
4029 }
4030 }
4031 }
4032
4033
4034 // 3. faces that do not get split but whose owner/neighbour change
4035 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4036
4037 if (debug)
4038 {
4039 Pout<< "hexRef8::setRefinement :"
4040 << " Changing owner/neighbour for otherwise unaffected faces"
4041 << endl;
4042 }
4043
4044 forAll(affectedFace, facei)
4045 {
4046 if (affectedFace.test(facei))
4047 {
4048 const face& f = mesh_.faces()[facei];
4049
4050 // The point with the lowest level should be an anchor
4051 // point of the neighbouring cells.
4052 label anchorFp = findMinLevel(f);
4053
4054 label own, nei;
4055 getFaceNeighbours
4056 (
4057 cellAnchorPoints,
4058 cellAddedCells,
4059 facei,
4060 f[anchorFp], // Anchor point
4061
4062 own,
4063 nei
4064 );
4065
4066 modFace(meshMod, facei, f, own, nei);
4067
4068 // Mark face as having been handled
4069 affectedFace.unset(facei);
4070 }
4071 }
4072
4073
4074 // 4. new internal faces inside split cells.
4075 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4076
4077
4078 // This is the hard one. We have to find the splitting points between
4079 // the anchor points. But the edges between the anchor points might have
4080 // been split (into two,three or four edges).
4081
4082 if (debug)
4083 {
4084 Pout<< "hexRef8::setRefinement :"
4085 << " Create new internal faces for split cells"
4086 << endl;
4087 }
4088
4089 forAll(cellMidPoint, celli)
4090 {
4091 if (cellMidPoint[celli] >= 0)
4092 {
4093 createInternalFaces
4094 (
4095 cellAnchorPoints,
4096 cellAddedCells,
4097 cellMidPoint,
4098 faceMidPoint,
4099 faceAnchorLevel,
4100 edgeMidPoint,
4101 celli,
4102 meshMod
4103 );
4104 }
4105 }
4106
4107 // Extend pointLevels and cellLevels for the new cells. Could also be done
4108 // in updateMesh but saves passing cellAddedCells out of this routine.
4109
4110 // Check
4111 if (debug)
4112 {
4113 label minPointi = labelMax;
4114 label maxPointi = labelMin;
4115
4116 forAll(cellMidPoint, celli)
4117 {
4118 if (cellMidPoint[celli] >= 0)
4119 {
4120 minPointi = min(minPointi, cellMidPoint[celli]);
4121 maxPointi = max(maxPointi, cellMidPoint[celli]);
4122 }
4123 }
4124 forAll(faceMidPoint, facei)
4125 {
4126 if (faceMidPoint[facei] >= 0)
4127 {
4128 minPointi = min(minPointi, faceMidPoint[facei]);
4129 maxPointi = max(maxPointi, faceMidPoint[facei]);
4130 }
4131 }
4132 forAll(edgeMidPoint, edgeI)
4133 {
4134 if (edgeMidPoint[edgeI] >= 0)
4135 {
4136 minPointi = min(minPointi, edgeMidPoint[edgeI]);
4137 maxPointi = max(maxPointi, edgeMidPoint[edgeI]);
4138 }
4139 }
4140
4141 if (minPointi != labelMax && minPointi != mesh_.nPoints())
4142 {
4144 << "Added point labels not consecutive to existing mesh points."
4145 << nl
4146 << "mesh_.nPoints():" << mesh_.nPoints()
4147 << " minPointi:" << minPointi
4148 << " maxPointi:" << maxPointi
4149 << abort(FatalError);
4150 }
4151 }
4152
4153 pointLevel_.transfer(newPointLevel);
4154 cellLevel_.transfer(newCellLevel);
4155
4156 // Mark files as changed
4157 setInstance(mesh_.facesInstance());
4158
4159
4160 // Update the live split cells tree.
4161 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4162
4163 // New unrefinement structure
4164 if (history_.active())
4165 {
4166 if (debug)
4167 {
4168 Pout<< "hexRef8::setRefinement :"
4169 << " Updating refinement history to " << cellLevel_.size()
4170 << " cells" << endl;
4171 }
4172
4173 // Extend refinement history for new cells
4174 history_.resize(cellLevel_.size());
4175
4176 forAll(cellAddedCells, celli)
4177 {
4178 const labelList& addedCells = cellAddedCells[celli];
4179
4180 if (addedCells.size())
4181 {
4182 // Cell was split.
4183 history_.storeSplit(celli, addedCells);
4184 }
4185 }
4186 }
4187
4188 // Compact cellAddedCells.
4189
4190 labelListList refinedCells(cellLabels.size());
4191
4192 forAll(cellLabels, i)
4193 {
4194 label celli = cellLabels[i];
4195
4196 refinedCells[i].transfer(cellAddedCells[celli]);
4197 }
4198
4199 return refinedCells;
4200}
4201
4202
4204(
4205 const labelList& pointsToStore,
4206 const labelList& facesToStore,
4207 const labelList& cellsToStore
4208)
4209{
4210 savedPointLevel_.clear();
4211 savedPointLevel_.resize(2*pointsToStore.size());
4212 forAll(pointsToStore, i)
4213 {
4214 label pointi = pointsToStore[i];
4215 savedPointLevel_.insert(pointi, pointLevel_[pointi]);
4216 }
4217
4218 savedCellLevel_.clear();
4219 savedCellLevel_.resize(2*cellsToStore.size());
4220 forAll(cellsToStore, i)
4221 {
4222 label celli = cellsToStore[i];
4223 savedCellLevel_.insert(celli, cellLevel_[celli]);
4224 }
4225}
4226
4227
4228// Gets called after the mesh change. setRefinement will already have made
4229// sure the pointLevel_ and cellLevel_ are the size of the new mesh so we
4230// only need to account for reordering.
4232{
4233 Map<label> dummyMap(0);
4234
4235 updateMesh(map, dummyMap, dummyMap, dummyMap);
4236}
4237
4238
4239// Gets called after the mesh change. setRefinement will already have made
4240// sure the pointLevel_ and cellLevel_ are the size of the new mesh so we
4241// only need to account for reordering.
4243(
4244 const mapPolyMesh& map,
4245 const Map<label>& pointsToRestore,
4246 const Map<label>& facesToRestore,
4247 const Map<label>& cellsToRestore
4248)
4249{
4250 // Update celllevel
4251 if (debug)
4252 {
4253 Pout<< "hexRef8::updateMesh :"
4254 << " Updating various lists"
4255 << endl;
4256 }
4257
4258 {
4259 const labelList& reverseCellMap = map.reverseCellMap();
4260
4261 if (debug)
4262 {
4263 Pout<< "hexRef8::updateMesh :"
4264 << " reverseCellMap:" << map.reverseCellMap().size()
4265 << " cellMap:" << map.cellMap().size()
4266 << " nCells:" << mesh_.nCells()
4267 << " nOldCells:" << map.nOldCells()
4268 << " cellLevel_:" << cellLevel_.size()
4269 << " reversePointMap:" << map.reversePointMap().size()
4270 << " pointMap:" << map.pointMap().size()
4271 << " nPoints:" << mesh_.nPoints()
4272 << " nOldPoints:" << map.nOldPoints()
4273 << " pointLevel_:" << pointLevel_.size()
4274 << endl;
4275 }
4276
4277 if (reverseCellMap.size() == cellLevel_.size())
4278 {
4279 // Assume it is after hexRef8 that this routine is called.
4280 // Just account for reordering. We cannot use cellMap since
4281 // then cells created from cells would get cellLevel_ of
4282 // cell they were created from.
4283 reorder(reverseCellMap, mesh_.nCells(), -1, cellLevel_);
4284 }
4285 else
4286 {
4287 // Map data
4288 const labelList& cellMap = map.cellMap();
4289
4290 labelList newCellLevel(cellMap.size());
4291 forAll(cellMap, newCelli)
4292 {
4293 label oldCelli = cellMap[newCelli];
4294
4295 if (oldCelli == -1)
4296 {
4297 newCellLevel[newCelli] = -1;
4298 }
4299 else
4300 {
4301 newCellLevel[newCelli] = cellLevel_[oldCelli];
4302 }
4303 }
4304 cellLevel_.transfer(newCellLevel);
4305 }
4306
4307 // See if any cells to restore. This will be for some new cells
4308 // the corresponding old cell.
4309 forAllConstIters(cellsToRestore, iter)
4310 {
4311 const label newCelli = iter.key();
4312 const label storedCelli = iter.val();
4313
4314 const auto fnd = savedCellLevel_.cfind(storedCelli);
4315
4316 if (!fnd.found())
4317 {
4319 << "Problem : trying to restore old value for new cell "
4320 << newCelli << " but cannot find old cell " << storedCelli
4321 << " in map of stored values " << savedCellLevel_
4322 << abort(FatalError);
4323 }
4324 cellLevel_[newCelli] = fnd.val();
4325 }
4326
4327 //if (cellLevel_.found(-1))
4328 //{
4329 // WarningInFunction
4330 // << "Problem : "
4331 // << "cellLevel_ contains illegal value -1 after mapping
4332 // << " at cell " << cellLevel_.find(-1) << endl
4333 // << "This means that another program has inflated cells"
4334 // << " (created cells out-of-nothing) and hence we don't know"
4335 // << " their cell level. Continuing with illegal value."
4336 // << abort(FatalError);
4337 //}
4338 }
4339
4340
4341 // Update pointlevel
4342 {
4343 const labelList& reversePointMap = map.reversePointMap();
4344
4345 if (reversePointMap.size() == pointLevel_.size())
4346 {
4347 // Assume it is after hexRef8 that this routine is called.
4348 reorder(reversePointMap, mesh_.nPoints(), -1, pointLevel_);
4349 }
4350 else
4351 {
4352 // Map data
4353 const labelList& pointMap = map.pointMap();
4354
4355 labelList newPointLevel(pointMap.size());
4356
4357 forAll(pointMap, newPointi)
4358 {
4359 const label oldPointi = pointMap[newPointi];
4360
4361 if (oldPointi == -1)
4362 {
4363 //FatalErrorInFunction
4364 // << "Problem : point " << newPointi
4365 // << " at " << mesh_.points()[newPointi]
4366 // << " does not originate from another point"
4367 // << " (i.e. is inflated)." << nl
4368 // << "Hence we cannot determine the new pointLevel"
4369 // << " for it." << abort(FatalError);
4370 newPointLevel[newPointi] = -1;
4371 }
4372 else
4373 {
4374 newPointLevel[newPointi] = pointLevel_[oldPointi];
4375 }
4376 }
4377 pointLevel_.transfer(newPointLevel);
4378 }
4379
4380 // See if any points to restore. This will be for some new points
4381 // the corresponding old point (the one from the call to storeData)
4382 forAllConstIters(pointsToRestore, iter)
4383 {
4384 const label newPointi = iter.key();
4385 const label storedPointi = iter.val();
4386
4387 auto fnd = savedPointLevel_.find(storedPointi);
4388
4389 if (!fnd.found())
4390 {
4392 << "Problem : trying to restore old value for new point "
4393 << newPointi << " but cannot find old point "
4394 << storedPointi
4395 << " in map of stored values " << savedPointLevel_
4396 << abort(FatalError);
4397 }
4398 pointLevel_[newPointi] = fnd.val();
4399 }
4400
4401 //if (pointLevel_.found(-1))
4402 //{
4403 // WarningInFunction
4404 // << "Problem : "
4405 // << "pointLevel_ contains illegal value -1 after mapping"
4406 // << " at point" << pointLevel_.find(-1) << endl
4407 // << "This means that another program has inflated points"
4408 // << " (created points out-of-nothing) and hence we don't know"
4409 // << " their point level. Continuing with illegal value."
4410 // //<< abort(FatalError);
4411 //}
4412 }
4413
4414 // Update refinement tree
4415 if (history_.active())
4416 {
4417 history_.updateMesh(map);
4418 }
4419
4420 // Mark files as changed
4421 setInstance(mesh_.facesInstance());
4422
4423 // Update face removal engine
4424 faceRemover_.updateMesh(map);
4425
4426 // Clear cell shapes
4427 cellShapesPtr_.clear();
4428}
4429
4430
4431// Gets called after mesh subsetting. Maps are from new back to old.
4433(
4434 const labelList& pointMap,
4435 const labelList& faceMap,
4436 const labelList& cellMap
4437)
4438{
4439 // Update celllevel
4440 if (debug)
4441 {
4442 Pout<< "hexRef8::subset :"
4443 << " Updating various lists"
4444 << endl;
4445 }
4446
4447 if (history_.active())
4448 {
4450 << "Subsetting will not work in combination with unrefinement."
4451 << nl
4452 << "Proceed at your own risk." << endl;
4453 }
4454
4455
4456 // Update celllevel
4457 {
4458 labelList newCellLevel(cellMap.size());
4459
4460 forAll(cellMap, newCelli)
4461 {
4462 newCellLevel[newCelli] = cellLevel_[cellMap[newCelli]];
4463 }
4464
4465 cellLevel_.transfer(newCellLevel);
4466
4467 if (cellLevel_.found(-1))
4468 {
4470 << "Problem : "
4471 << "cellLevel_ contains illegal value -1 after mapping:"
4472 << cellLevel_
4473 << abort(FatalError);
4474 }
4475 }
4476
4477 // Update pointlevel
4478 {
4479 labelList newPointLevel(pointMap.size());
4480
4481 forAll(pointMap, newPointi)
4482 {
4483 newPointLevel[newPointi] = pointLevel_[pointMap[newPointi]];
4484 }
4485
4486 pointLevel_.transfer(newPointLevel);
4487
4488 if (pointLevel_.found(-1))
4489 {
4491 << "Problem : "
4492 << "pointLevel_ contains illegal value -1 after mapping:"
4493 << pointLevel_
4494 << abort(FatalError);
4495 }
4496 }
4497
4498 // Update refinement tree
4499 if (history_.active())
4500 {
4501 history_.subset(pointMap, faceMap, cellMap);
4502 }
4503
4504 // Mark files as changed
4505 setInstance(mesh_.facesInstance());
4506
4507 // Nothing needs doing to faceRemover.
4508 //faceRemover_.subset(pointMap, faceMap, cellMap);
4509
4510 // Clear cell shapes
4511 cellShapesPtr_.clear();
4512}
4513
4514
4515// Gets called after the mesh distribution
4517{
4518 if (debug)
4519 {
4520 Pout<< "hexRef8::distribute :"
4521 << " Updating various lists"
4522 << endl;
4523 }
4524
4525 // Update celllevel
4526 map.distributeCellData(cellLevel_);
4527 // Update pointlevel
4528 map.distributePointData(pointLevel_);
4529
4530 // Update refinement tree
4531 if (history_.active())
4532 {
4533 history_.distribute(map);
4534 }
4535
4536 // Update face removal engine
4537 faceRemover_.distribute(map);
4538
4539 // Clear cell shapes
4540 cellShapesPtr_.clear();
4541}
4542
4543
4545{
4546 const scalar smallDim = 1e-6 * mesh_.bounds().mag();
4547
4548 if (debug)
4549 {
4550 Pout<< "hexRef8::checkMesh : Using matching tolerance smallDim:"
4551 << smallDim << endl;
4552 }
4553
4554 // Check owner on coupled faces.
4555 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4556
4557 // There should be only one coupled face between two cells. Why? Since
4558 // otherwise mesh redistribution might cause multiple faces between two
4559 // cells
4560 {
4561 labelList nei(mesh_.nBoundaryFaces());
4562 forAll(nei, i)
4563 {
4564 nei[i] = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4565 }
4566
4567 // Replace data on coupled patches with their neighbour ones.
4569
4570 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
4571
4572 forAll(patches, patchi)
4573 {
4574 const polyPatch& pp = patches[patchi];
4575
4576 if (pp.coupled())
4577 {
4578 // Check how many faces between owner and neighbour.
4579 // Should be only one.
4580 labelPairLookup cellToFace(2*pp.size());
4581
4582 label facei = pp.start();
4583
4584 forAll(pp, i)
4585 {
4586 label own = mesh_.faceOwner()[facei];
4587 label bFacei = facei-mesh_.nInternalFaces();
4588
4589 if (!cellToFace.insert(labelPair(own, nei[bFacei]), facei))
4590 {
4591 dumpCell(own);
4593 << "Faces do not seem to be correct across coupled"
4594 << " boundaries" << endl
4595 << "Coupled face " << facei
4596 << " between owner " << own
4597 << " on patch " << pp.name()
4598 << " and coupled neighbour " << nei[bFacei]
4599 << " has two faces connected to it:"
4600 << facei << " and "
4601 << cellToFace[labelPair(own, nei[bFacei])]
4602 << abort(FatalError);
4603 }
4604
4605 facei++;
4606 }
4607 }
4608 }
4609 }
4610
4611 // Check face areas.
4612 // ~~~~~~~~~~~~~~~~~
4613
4614 {
4615 scalarField neiFaceAreas(mesh_.nBoundaryFaces());
4616 forAll(neiFaceAreas, i)
4617 {
4618 neiFaceAreas[i] = mag(mesh_.faceAreas()[i+mesh_.nInternalFaces()]);
4619 }
4620
4621 // Replace data on coupled patches with their neighbour ones.
4622 syncTools::swapBoundaryFaceList(mesh_, neiFaceAreas);
4623
4624 forAll(neiFaceAreas, i)
4625 {
4626 label facei = i+mesh_.nInternalFaces();
4627
4628 const scalar magArea = mag(mesh_.faceAreas()[facei]);
4629
4630 if (mag(magArea - neiFaceAreas[i]) > smallDim)
4631 {
4632 const face& f = mesh_.faces()[facei];
4633 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4634
4635 dumpCell(mesh_.faceOwner()[facei]);
4636
4638 << "Faces do not seem to be correct across coupled"
4639 << " boundaries" << endl
4640 << "Coupled face " << facei
4641 << " on patch " << patchi
4642 << " " << mesh_.boundaryMesh()[patchi].name()
4643 << " coords:" << UIndirectList<point>(mesh_.points(), f)
4644 << " has face area:" << magArea
4645 << " (coupled) neighbour face area differs:"
4646 << neiFaceAreas[i]
4647 << " to within tolerance " << smallDim
4648 << abort(FatalError);
4649 }
4650 }
4651 }
4652
4653
4654 // Check number of points on faces.
4655 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4656 {
4657 labelList nVerts(mesh_.nBoundaryFaces());
4658
4659 forAll(nVerts, i)
4660 {
4661 nVerts[i] = mesh_.faces()[i+mesh_.nInternalFaces()].size();
4662 }
4663
4664 // Replace data on coupled patches with their neighbour ones.
4665 syncTools::swapBoundaryFaceList(mesh_, nVerts);
4666
4667 forAll(nVerts, i)
4668 {
4669 label facei = i+mesh_.nInternalFaces();
4670
4671 const face& f = mesh_.faces()[facei];
4672
4673 if (f.size() != nVerts[i])
4674 {
4675 dumpCell(mesh_.faceOwner()[facei]);
4676
4677 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4678
4680 << "Faces do not seem to be correct across coupled"
4681 << " boundaries" << endl
4682 << "Coupled face " << facei
4683 << " on patch " << patchi
4684 << " " << mesh_.boundaryMesh()[patchi].name()
4685 << " coords:" << UIndirectList<point>(mesh_.points(), f)
4686 << " has size:" << f.size()
4687 << " (coupled) neighbour face has size:"
4688 << nVerts[i]
4689 << abort(FatalError);
4690 }
4691 }
4692 }
4693
4694
4695 // Check points of face
4696 // ~~~~~~~~~~~~~~~~~~~~
4697 {
4698 // Anchor points.
4699 pointField anchorPoints(mesh_.nBoundaryFaces());
4700
4701 forAll(anchorPoints, i)
4702 {
4703 label facei = i+mesh_.nInternalFaces();
4704 const point& fc = mesh_.faceCentres()[facei];
4705 const face& f = mesh_.faces()[facei];
4706 const vector anchorVec(mesh_.points()[f[0]] - fc);
4707
4708 anchorPoints[i] = anchorVec;
4709 }
4710
4711 // Replace data on coupled patches with their neighbour ones. Apply
4712 // rotation transformation (but not separation since is relative vector
4713 // to point on same face.
4714 syncTools::swapBoundaryFaceList(mesh_, anchorPoints);
4715
4716 forAll(anchorPoints, i)
4717 {
4718 label facei = i+mesh_.nInternalFaces();
4719 const point& fc = mesh_.faceCentres()[facei];
4720 const face& f = mesh_.faces()[facei];
4721 const vector anchorVec(mesh_.points()[f[0]] - fc);
4722
4723 if (mag(anchorVec - anchorPoints[i]) > smallDim)
4724 {
4725 dumpCell(mesh_.faceOwner()[facei]);
4726
4727 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4728
4730 << "Faces do not seem to be correct across coupled"
4731 << " boundaries" << endl
4732 << "Coupled face " << facei
4733 << " on patch " << patchi
4734 << " " << mesh_.boundaryMesh()[patchi].name()
4735 << " coords:" << UIndirectList<point>(mesh_.points(), f)
4736 << " has anchor vector:" << anchorVec
4737 << " (coupled) neighbour face anchor vector differs:"
4738 << anchorPoints[i]
4739 << " to within tolerance " << smallDim
4740 << abort(FatalError);
4741 }
4742 }
4743 }
4744
4745 if (debug)
4746 {
4747 Pout<< "hexRef8::checkMesh : Returning" << endl;
4748 }
4749}
4750
4751
4753(
4754 const label maxPointDiff,
4755 const labelList& pointsToCheck
4756) const
4757{
4758 if (debug)
4759 {
4760 Pout<< "hexRef8::checkRefinementLevels :"
4761 << " Checking 2:1 refinement level" << endl;
4762 }
4763
4764 if
4765 (
4766 cellLevel_.size() != mesh_.nCells()
4767 || pointLevel_.size() != mesh_.nPoints()
4768 )
4769 {
4771 << "cellLevel size should be number of cells"
4772 << " and pointLevel size should be number of points."<< nl
4773 << "cellLevel:" << cellLevel_.size()
4774 << " mesh.nCells():" << mesh_.nCells() << nl
4775 << "pointLevel:" << pointLevel_.size()
4776 << " mesh.nPoints():" << mesh_.nPoints()
4777 << abort(FatalError);
4778 }
4779
4780
4781 // Check 2:1 consistency.
4782 // ~~~~~~~~~~~~~~~~~~~~~~
4783
4784 {
4785 // Internal faces.
4786 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
4787 {
4788 label own = mesh_.faceOwner()[facei];
4789 label nei = mesh_.faceNeighbour()[facei];
4790
4791 if (mag(cellLevel_[own] - cellLevel_[nei]) > 1)
4792 {
4793 dumpCell(own);
4794 dumpCell(nei);
4795
4797 << "Celllevel does not satisfy 2:1 constraint." << nl
4798 << "On face " << facei << " owner cell " << own
4799 << " has refinement " << cellLevel_[own]
4800 << " neighbour cell " << nei << " has refinement "
4801 << cellLevel_[nei]
4802 << abort(FatalError);
4803 }
4804 }
4805
4806 // Coupled faces. Get neighbouring value
4807 labelList neiLevel(mesh_.nBoundaryFaces());
4808
4809 forAll(neiLevel, i)
4810 {
4811 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4812
4813 neiLevel[i] = cellLevel_[own];
4814 }
4815
4816 // No separation
4817 syncTools::swapBoundaryFaceList(mesh_, neiLevel);
4818
4819 forAll(neiLevel, i)
4820 {
4821 label facei = i+mesh_.nInternalFaces();
4822
4823 label own = mesh_.faceOwner()[facei];
4824
4825 if (mag(cellLevel_[own] - neiLevel[i]) > 1)
4826 {
4827 dumpCell(own);
4828
4829 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4830
4832 << "Celllevel does not satisfy 2:1 constraint."
4833 << " On coupled face " << facei
4834 << " on patch " << patchi << " "
4835 << mesh_.boundaryMesh()[patchi].name()
4836 << " owner cell " << own << " has refinement "
4837 << cellLevel_[own]
4838 << " (coupled) neighbour cell has refinement "
4839 << neiLevel[i]
4840 << abort(FatalError);
4841 }
4842 }
4843 }
4844
4845
4846 // Check pointLevel is synchronized
4847 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4848 {
4849 labelList syncPointLevel(pointLevel_);
4850
4851 // Get min level
4853 (
4854 mesh_,
4855 syncPointLevel,
4857 labelMax
4858 );
4859
4860
4861 forAll(syncPointLevel, pointi)
4862 {
4863 if (pointLevel_[pointi] != syncPointLevel[pointi])
4864 {
4866 << "PointLevel is not consistent across coupled patches."
4867 << endl
4868 << "point:" << pointi << " coord:" << mesh_.points()[pointi]
4869 << " has level " << pointLevel_[pointi]
4870 << " whereas the coupled point has level "
4871 << syncPointLevel[pointi]
4872 << abort(FatalError);
4873 }
4874 }
4875 }
4876
4877
4878 // Check 2:1 across points (instead of faces)
4879 if (maxPointDiff != -1)
4880 {
4881 // Determine per point the max cell level.
4882 labelList maxPointLevel(mesh_.nPoints(), Zero);
4883
4884 forAll(maxPointLevel, pointi)
4885 {
4886 const labelList& pCells = mesh_.pointCells(pointi);
4887
4888 label& pLevel = maxPointLevel[pointi];
4889
4890 forAll(pCells, i)
4891 {
4892 pLevel = max(pLevel, cellLevel_[pCells[i]]);
4893 }
4894 }
4895
4896 // Sync maxPointLevel to neighbour
4898 (
4899 mesh_,
4900 maxPointLevel,
4902 labelMin // null value
4903 );
4904
4905 // Check 2:1 across boundary points
4906 forAll(pointsToCheck, i)
4907 {
4908 label pointi = pointsToCheck[i];
4909
4910 const labelList& pCells = mesh_.pointCells(pointi);
4911
4912 forAll(pCells, i)
4913 {
4914 label celli = pCells[i];
4915
4916 if
4917 (
4918 mag(cellLevel_[celli]-maxPointLevel[pointi])
4919 > maxPointDiff
4920 )
4921 {
4922 dumpCell(celli);
4923
4925 << "Too big a difference between"
4926 << " point-connected cells." << nl
4927 << "cell:" << celli
4928 << " cellLevel:" << cellLevel_[celli]
4929 << " uses point:" << pointi
4930 << " coord:" << mesh_.points()[pointi]
4931 << " which is also used by a cell with level:"
4932 << maxPointLevel[pointi]
4933 << abort(FatalError);
4934 }
4935 }
4936 }
4937 }
4938
4939
4940 //- Gives problems after first splitting off inside mesher.
4942 //{
4943 // // Any patches with points having only two edges.
4944 //
4945 // boolList isHangingPoint(mesh_.nPoints(), false);
4946 //
4947 // const polyBoundaryMesh& patches = mesh_.boundaryMesh();
4948 //
4949 // forAll(patches, patchi)
4950 // {
4951 // const polyPatch& pp = patches[patchi];
4952 //
4953 // const labelList& meshPoints = pp.meshPoints();
4954 //
4955 // forAll(meshPoints, i)
4956 // {
4957 // label pointi = meshPoints[i];
4958 //
4959 // const labelList& pEdges = mesh_.pointEdges()[pointi];
4960 //
4961 // if (pEdges.size() == 2)
4962 // {
4963 // isHangingPoint[pointi] = true;
4964 // }
4965 // }
4966 // }
4967 //
4968 // syncTools::syncPointList
4969 // (
4970 // mesh_,
4971 // isHangingPoint,
4972 // andEqOp<bool>(), // only if all decide it is hanging point
4973 // true, // null
4974 // false // no separation
4975 // );
4976 //
4977 // //OFstream str(mesh_.time().path()/"hangingPoints.obj");
4978 //
4979 // label nHanging = 0;
4980 //
4981 // forAll(isHangingPoint, pointi)
4982 // {
4983 // if (isHangingPoint[pointi])
4984 // {
4985 // nHanging++;
4986 //
4987 // Pout<< "Hanging boundary point " << pointi
4988 // << " at " << mesh_.points()[pointi]
4989 // << endl;
4990 // //meshTools::writeOBJ(str, mesh_.points()[pointi]);
4991 // }
4992 // }
4993 //
4994 // if (returnReduce(nHanging, sumOp<label>()) > 0)
4995 // {
4996 // FatalErrorInFunction
4997 // << "Detected a point used by two edges only (hanging point)"
4998 // << nl << "This is not allowed"
4999 // << abort(FatalError);
5000 // }
5001 //}
5002}
5003
5004
5006{
5007 if (!cellShapesPtr_)
5008 {
5009 if (debug)
5010 {
5011 Pout<< "hexRef8::cellShapes() : calculating splitHex cellShapes."
5012 << " cellLevel:" << cellLevel_.size()
5013 << " pointLevel:" << pointLevel_.size()
5014 << endl;
5015 }
5016
5017 const cellShapeList& meshShapes = mesh_.cellShapes();
5018 cellShapesPtr_.reset(new cellShapeList(meshShapes));
5019
5020 label nSplitHex = 0;
5021 label nUnrecognised = 0;
5022
5023 forAll(cellLevel_, celli)
5024 {
5025 if (meshShapes[celli].model().index() == 0)
5026 {
5027 label level = cellLevel_[celli];
5028
5029 // Return true if we've found 6 quads
5030 DynamicList<face> quads;
5031 bool haveQuads = matchHexShape
5032 (
5033 celli,
5034 level,
5035 quads
5036 );
5037
5038 if (haveQuads)
5039 {
5040 faceList faces(std::move(quads));
5041 cellShapesPtr_()[celli] = degenerateMatcher::match(faces);
5042 nSplitHex++;
5043 }
5044 else
5045 {
5046 nUnrecognised++;
5047 }
5048 }
5049 }
5050 if (debug)
5051 {
5052 Pout<< "hexRef8::cellShapes() :"
5053 << " nCells:" << mesh_.nCells() << " of which" << nl
5054 << " primitive:" << (mesh_.nCells()-nSplitHex-nUnrecognised)
5055 << nl
5056 << " split-hex:" << nSplitHex << nl
5057 << " poly :" << nUnrecognised << nl
5058 << endl;
5059 }
5060 }
5061 return *cellShapesPtr_;
5062}
5063
5064
5066{
5067 if (debug)
5068 {
5069 checkRefinementLevels(-1, labelList(0));
5070 }
5071
5072 if (debug)
5073 {
5074 Pout<< "hexRef8::getSplitPoints :"
5075 << " Calculating unrefineable points" << endl;
5076 }
5077
5078
5079 if (!history_.active())
5080 {
5082 << "Only call if constructed with history capability"
5083 << abort(FatalError);
5084 }
5085
5086 // Master cell
5087 // -1 undetermined
5088 // -2 certainly not split point
5089 // >= label of master cell
5090 labelList splitMaster(mesh_.nPoints(), -1);
5091 labelList splitMasterLevel(mesh_.nPoints(), Zero);
5092
5093 // Unmark all with not 8 cells
5094 //const labelListList& pointCells = mesh_.pointCells();
5095
5096 for (label pointi = 0; pointi < mesh_.nPoints(); pointi++)
5097 {
5098 const labelList& pCells = mesh_.pointCells(pointi);
5099
5100 if (pCells.size() != 8)
5101 {
5102 splitMaster[pointi] = -2;
5103 }
5104 }
5105
5106 // Unmark all with different master cells
5107 const labelList& visibleCells = history_.visibleCells();
5108
5109 forAll(visibleCells, celli)
5110 {
5111 const labelList& cPoints = mesh_.cellPoints(celli);
5112
5113 if (visibleCells[celli] != -1 && history_.parentIndex(celli) >= 0)
5114 {
5115 label parentIndex = history_.parentIndex(celli);
5116
5117 // Check same master.
5118 forAll(cPoints, i)
5119 {
5120 label pointi = cPoints[i];
5121
5122 label masterCelli = splitMaster[pointi];
5123
5124 if (masterCelli == -1)
5125 {
5126 // First time visit of point. Store parent cell and
5127 // level of the parent cell (with respect to celli). This
5128 // is additional guarantee that we're referring to the
5129 // same master at the same refinement level.
5130
5131 splitMaster[pointi] = parentIndex;
5132 splitMasterLevel[pointi] = cellLevel_[celli] - 1;
5133 }
5134 else if (masterCelli == -2)
5135 {
5136 // Already decided that point is not splitPoint
5137 }
5138 else if
5139 (
5140 (masterCelli != parentIndex)
5141 || (splitMasterLevel[pointi] != cellLevel_[celli] - 1)
5142 )
5143 {
5144 // Different masters so point is on two refinement
5145 // patterns
5146 splitMaster[pointi] = -2;
5147 }
5148 }
5149 }
5150 else
5151 {
5152 // Either not visible or is unrefined cell
5153 forAll(cPoints, i)
5154 {
5155 label pointi = cPoints[i];
5156
5157 splitMaster[pointi] = -2;
5158 }
5159 }
5160 }
5161
5162 // Unmark boundary faces
5163 for
5164 (
5165 label facei = mesh_.nInternalFaces();
5166 facei < mesh_.nFaces();
5167 facei++
5168 )
5169 {
5170 const face& f = mesh_.faces()[facei];
5171
5172 forAll(f, fp)
5173 {
5174 splitMaster[f[fp]] = -2;
5175 }
5176 }
5177
5178
5179 // Collect into labelList
5180
5181 label nSplitPoints = 0;
5182
5183 forAll(splitMaster, pointi)
5184 {
5185 if (splitMaster[pointi] >= 0)
5186 {
5187 nSplitPoints++;
5188 }
5189 }
5190
5191 labelList splitPoints(nSplitPoints);
5192 nSplitPoints = 0;
5193
5194 forAll(splitMaster, pointi)
5195 {
5196 if (splitMaster[pointi] >= 0)
5197 {
5198 splitPoints[nSplitPoints++] = pointi;
5199 }
5200 }
5201
5202 return splitPoints;
5203}
5204
5205
5206//void Foam::hexRef8::markIndex
5207//(
5208// const label maxLevel,
5209// const label level,
5210// const label index,
5211// const label markValue,
5212// labelList& indexValues
5213//) const
5214//{
5215// if (level < maxLevel && indexValues[index] == -1)
5216// {
5217// // Mark
5218// indexValues[index] = markValue;
5219//
5220// // Mark parent
5221// const splitCell8& split = history_.splitCells()[index];
5222//
5223// if (split.parent_ >= 0)
5224// {
5225// markIndex
5226// (
5227// maxLevel, level+1, split.parent_, markValue, indexValues);
5228// )
5229// }
5230// }
5231//}
5232//
5233//
5238//void Foam::hexRef8::markCellClusters
5239//(
5240// const label maxLevel,
5241// labelList& cluster
5242//) const
5243//{
5244// cluster.setSize(mesh_.nCells());
5245// cluster = -1;
5246//
5247// const DynamicList<splitCell8>& splitCells = history_.splitCells();
5248//
5249// // Mark all splitCells down to level maxLevel with a cell originating from
5250// // it.
5251//
5252// labelList indexLevel(splitCells.size(), -1);
5253//
5254// forAll(visibleCells, celli)
5255// {
5256// label index = visibleCells[celli];
5257//
5258// if (index >= 0)
5259// {
5260// markIndex(maxLevel, 0, index, celli, indexLevel);
5261// }
5262// }
5263//
5264// // Mark cells with splitCell
5265//}
5266
5267
5269(
5270 const labelList& pointsToUnrefine,
5271 const bool maxSet
5272) const
5273{
5274 if (debug)
5275 {
5276 Pout<< "hexRef8::consistentUnrefinement :"
5277 << " Determining 2:1 consistent unrefinement" << endl;
5278 }
5279
5280 if (maxSet)
5281 {
5283 << "maxSet not implemented yet."
5284 << abort(FatalError);
5285 }
5286
5287 // Loop, modifying pointsToUnrefine, until no more changes to due to 2:1
5288 // conflicts.
5289 // maxSet = false : unselect points to refine
5290 // maxSet = true: select points to refine
5291
5292 // Maintain bitset for pointsToUnrefine and cellsToUnrefine
5293 bitSet unrefinePoint(mesh_.nPoints(), pointsToUnrefine);
5294
5295 while (true)
5296 {
5297 // Construct cells to unrefine
5298 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
5299
5300 bitSet unrefineCell(mesh_.nCells());
5301
5302 forAll(unrefinePoint, pointi)
5303 {
5304 if (unrefinePoint.test(pointi))
5305 {
5306 const labelList& pCells = mesh_.pointCells(pointi);
5307
5308 unrefineCell.set(pCells);
5309 }
5310 }
5311
5312
5313 label nChanged = 0;
5314
5315
5316 // Check 2:1 consistency taking refinement into account
5317 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5318
5319 // Internal faces.
5320 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
5321 {
5322 label own = mesh_.faceOwner()[facei];
5323 label nei = mesh_.faceNeighbour()[facei];
5324
5325 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5326 label neiLevel = cellLevel_[nei] - unrefineCell.get(nei);
5327
5328 if (ownLevel < (neiLevel-1))
5329 {
5330 // Since was 2:1 this can only occur if own is marked for
5331 // unrefinement.
5332
5333 if (maxSet)
5334 {
5335 unrefineCell.set(nei);
5336 }
5337 else
5338 {
5339 // could also combine with unset:
5340 // if (!unrefineCell.unset(own))
5341 // {
5342 // FatalErrorInFunction
5343 // << "problem cell already unset"
5344 // << abort(FatalError);
5345 // }
5346 if (!unrefineCell.test(own))
5347 {
5349 << "problem" << abort(FatalError);
5350 }
5351
5352 unrefineCell.unset(own);
5353 }
5354 nChanged++;
5355 }
5356 else if (neiLevel < (ownLevel-1))
5357 {
5358 if (maxSet)
5359 {
5360 unrefineCell.set(own);
5361 }
5362 else
5363 {
5364 if (!unrefineCell.test(nei))
5365 {
5367 << "problem" << abort(FatalError);
5368 }
5369
5370 unrefineCell.unset(nei);
5371 }
5372 nChanged++;
5373 }
5374 }
5375
5376
5377 // Coupled faces. Swap owner level to get neighbouring cell level.
5378 labelList neiLevel(mesh_.nBoundaryFaces());
5379
5380 forAll(neiLevel, i)
5381 {
5382 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
5383
5384 neiLevel[i] = cellLevel_[own] - unrefineCell.get(own);
5385 }
5386
5387 // Swap to neighbour
5388 syncTools::swapBoundaryFaceList(mesh_, neiLevel);
5389
5390 forAll(neiLevel, i)
5391 {
5392 label facei = i+mesh_.nInternalFaces();
5393 label own = mesh_.faceOwner()[facei];
5394 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5395
5396 if (ownLevel < (neiLevel[i]-1))
5397 {
5398 if (!maxSet)
5399 {
5400 if (!unrefineCell.test(own))
5401 {
5403 << "problem" << abort(FatalError);
5404 }
5405
5406 unrefineCell.unset(own);
5407 nChanged++;
5408 }
5409 }
5410 else if (neiLevel[i] < (ownLevel-1))
5411 {
5412 if (maxSet)
5413 {
5414 if (unrefineCell.test(own))
5415 {
5417 << "problem" << abort(FatalError);
5418 }
5419
5420 unrefineCell.set(own);
5421 nChanged++;
5422 }
5423 }
5424 }
5425
5426 reduce(nChanged, sumOp<label>());
5427
5428 if (debug)
5429 {
5430 Pout<< "hexRef8::consistentUnrefinement :"
5431 << " Changed " << nChanged
5432 << " refinement levels due to 2:1 conflicts."
5433 << endl;
5434 }
5435
5436 if (nChanged == 0)
5437 {
5438 break;
5439 }
5440
5441
5442 // Convert cellsToUnrefine back into points to unrefine
5443 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5444
5445 // Knock out any point whose cell neighbour cannot be unrefined.
5446 forAll(unrefinePoint, pointi)
5447 {
5448 if (unrefinePoint.test(pointi))
5449 {
5450 const labelList& pCells = mesh_.pointCells(pointi);
5451
5452 forAll(pCells, j)
5453 {
5454 if (!unrefineCell.test(pCells[j]))
5455 {
5456 unrefinePoint.unset(pointi);
5457 break;
5458 }
5459 }
5460 }
5461 }
5462 }
5463
5464
5465 // Convert back to labelList.
5466 label nSet = 0;
5467
5468 forAll(unrefinePoint, pointi)
5469 {
5470 if (unrefinePoint.test(pointi))
5471 {
5472 nSet++;
5473 }
5474 }
5475
5476 labelList newPointsToUnrefine(nSet);
5477 nSet = 0;
5478
5479 forAll(unrefinePoint, pointi)
5480 {
5481 if (unrefinePoint.test(pointi))
5482 {
5483 newPointsToUnrefine[nSet++] = pointi;
5484 }
5485 }
5486
5487 return newPointsToUnrefine;
5488}
5489
5490
5492(
5493 const labelList& splitPointLabels,
5494 polyTopoChange& meshMod
5495)
5496{
5497 if (!history_.active())
5498 {
5500 << "Only call if constructed with history capability"
5501 << abort(FatalError);
5502 }
5503
5504 if (debug)
5505 {
5506 Pout<< "hexRef8::setUnrefinement :"
5507 << " Checking initial mesh just to make sure" << endl;
5508
5509 checkMesh();
5510
5511 forAll(cellLevel_, celli)
5512 {
5513 if (cellLevel_[celli] < 0)
5514 {
5516 << "Illegal cell level " << cellLevel_[celli]
5517 << " for cell " << celli
5518 << abort(FatalError);
5519 }
5520 }
5521
5522
5523 // Write to sets.
5524 pointSet pSet(mesh_, "splitPoints", splitPointLabels);
5525 pSet.write();
5526
5527 cellSet cSet(mesh_, "splitPointCells", splitPointLabels.size());
5528
5529 forAll(splitPointLabels, i)
5530 {
5531 const labelList& pCells = mesh_.pointCells(splitPointLabels[i]);
5532
5533 cSet.insert(pCells);
5534 }
5535 cSet.write();
5536
5537 Pout<< "hexRef8::setRefinement : Dumping " << pSet.size()
5538 << " points and "
5539 << cSet.size() << " cells for unrefinement to" << nl
5540 << " pointSet " << pSet.objectPath() << nl
5541 << " cellSet " << cSet.objectPath()
5542 << endl;
5543 }
5544
5545
5546 labelList cellRegion;
5547 labelList cellRegionMaster;
5548 labelList facesToRemove;
5549
5550 {
5551 labelHashSet splitFaces(12*splitPointLabels.size());
5552
5553 forAll(splitPointLabels, i)
5554 {
5555 const labelList& pFaces = mesh_.pointFaces()[splitPointLabels[i]];
5556
5557 splitFaces.insert(pFaces);
5558 }
5559
5560 // Check with faceRemover what faces will get removed. Note that this
5561 // can be more (but never less) than splitFaces provided.
5562 faceRemover_.compatibleRemoves
5563 (
5564 splitFaces.toc(), // pierced faces
5565 cellRegion, // per cell -1 or region it is merged into
5566 cellRegionMaster, // per region the master cell
5567 facesToRemove // new faces to be removed.
5568 );
5569
5570 if (facesToRemove.size() != splitFaces.size())
5571 {
5572 // Dump current mesh
5573 {
5574 const_cast<polyMesh&>(mesh_).setInstance
5575 (
5576 mesh_.time().timeName()
5577 );
5578 mesh_.write();
5579 pointSet pSet(mesh_, "splitPoints", splitPointLabels);
5580 pSet.write();
5581 faceSet fSet(mesh_, "splitFaces", splitFaces);
5582 fSet.write();
5583 faceSet removeSet(mesh_, "facesToRemove", facesToRemove);
5584 removeSet.write();
5585 }
5586
5588 << "Ininitial set of split points to unrefine does not"
5589 << " seem to be consistent or not mid points of refined cells"
5590 << " splitPoints:" << splitPointLabels.size()
5591 << " splitFaces:" << splitFaces.size()
5592 << " facesToRemove:" << facesToRemove.size()
5593 << abort(FatalError);
5594 }
5595 }
5596
5597 // Redo the region master so it is consistent with our master.
5598 // This will guarantee that the new cell (for which faceRemover uses
5599 // the region master) is already compatible with our refinement structure.
5600
5601 forAll(splitPointLabels, i)
5602 {
5603 label pointi = splitPointLabels[i];
5604
5605 // Get original cell label
5606
5607 const labelList& pCells = mesh_.pointCells(pointi);
5608
5609 // Check
5610 if (pCells.size() != 8)
5611 {
5613 << "splitPoint " << pointi
5614 << " should have 8 cells using it. It has " << pCells
5615 << abort(FatalError);
5616 }
5617
5618
5619 // Check that the lowest numbered pCells is the master of the region
5620 // (should be guaranteed by directRemoveFaces)
5621 //if (debug)
5622 {
5623 label masterCelli = min(pCells);
5624
5625 forAll(pCells, j)
5626 {
5627 label celli = pCells[j];
5628
5629 label region = cellRegion[celli];
5630
5631 if (region == -1)
5632 {
5634 << "Ininitial set of split points to unrefine does not"
5635 << " seem to be consistent or not mid points"
5636 << " of refined cells" << nl
5637 << "cell:" << celli << " on splitPoint " << pointi
5638 << " has no region to be merged into"
5639 << abort(FatalError);
5640 }
5641
5642 if (masterCelli != cellRegionMaster[region])
5643 {
5645 << "cell:" << celli << " on splitPoint:" << pointi
5646 << " in region " << region
5647 << " has master:" << cellRegionMaster[region]
5648 << " which is not the lowest numbered cell"
5649 << " among the pointCells:" << pCells
5650 << abort(FatalError);
5651 }
5652 }
5653 }
5654 }
5655
5656 // Insert all commands to combine cells. Never fails so don't have to
5657 // test for success.
5658 faceRemover_.setRefinement
5659 (
5660 facesToRemove,
5661 cellRegion,
5662 cellRegionMaster,
5663 meshMod
5664 );
5665
5666 // Remove the 8 cells that originated from merging around the split point
5667 // and adapt cell levels (not that pointLevels stay the same since points
5668 // either get removed or stay at the same position.
5669 forAll(splitPointLabels, i)
5670 {
5671 label pointi = splitPointLabels[i];
5672
5673 const labelList& pCells = mesh_.pointCells(pointi);
5674
5675 label masterCelli = min(pCells);
5676
5677 forAll(pCells, j)
5678 {
5679 cellLevel_[pCells[j]]--;
5680 }
5681
5682 history_.combineCells(masterCelli, pCells);
5683 }
5684
5685 // Mark files as changed
5686 setInstance(mesh_.facesInstance());
5687
5688 // history_.updateMesh will take care of truncating.
5689}
5690
5691
5692// Write refinement to polyMesh directory.
5693bool Foam::hexRef8::write(const bool valid) const
5694{
5695 bool writeOk =
5696 cellLevel_.write(valid)
5697 && pointLevel_.write(valid)
5698 && level0Edge_.write(valid);
5699
5700 if (history_.active())
5701 {
5702 writeOk = writeOk && history_.write(valid);
5703 }
5704 else
5705 {
5707 }
5708
5709 return writeOk;
5710}
5711
5712
5714{
5715 IOobject io
5716 (
5717 "dummy",
5720 mesh
5721 );
5722 fileName setsDir(io.path());
5723
5724 if (topoSet::debug) DebugVar(setsDir);
5725
5726 if (exists(setsDir/"cellLevel"))
5727 {
5728 rm(setsDir/"cellLevel");
5729 }
5730 if (exists(setsDir/"pointLevel"))
5731 {
5732 rm(setsDir/"pointLevel");
5733 }
5734 if (exists(setsDir/"level0Edge"))
5735 {
5736 rm(setsDir/"level0Edge");
5737 }
5739}
5740
5741
5742// ************************************************************************* //
scalar y
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void clear() noexcept
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:391
DynamicList< T, SizeMin > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:434
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
void setCapacity(const label len)
Alter the size of the underlying storage.
Definition: DynamicListI.H:303
Wave propagation of information through grid. Every iteration information goes through one layer of c...
Definition: FaceCellWave.H:81
static scalar propagationTol()
Access to tolerance.
Definition: FaceCellWave.H:271
virtual label iterate(const label maxIter)
Iterate until no changes or maxIter reached.
void setFaceInfo(const label facei, const Type &faceInfo)
Set single initial changed face.
Definition: FaceCellWave.C:318
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:191
List< Key > toc() const
The table of contents (the keys) in unsorted order.
Definition: HashTable.C:122
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
Definition: HashTableI.H:180
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
readOption readOpt() const noexcept
The read option.
Definition: IOobjectI.H:164
bool typeHeaderOk(const bool checkType=true, const bool search=true, const bool verbose=true)
Read header (uses typeFilePath to find file) and check its info.
fileName path() const
The complete path.
Definition: IOobject.C:524
fileName objectPath() const
The complete path + object name.
Definition: IOobjectI.H:214
virtual const fileName & name() const
Return the name of the stream.
Definition: IOstream.C:40
void transfer(List< T > &list)
Definition: List.C:447
void setSize(const label n)
Alias for resize()
Definition: List.H:218
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:116
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
Output to file stream, using an OSstream.
Definition: OFstream.H:57
virtual const fileName & name() const
Read/write access to the name of the stream.
Definition: OSstream.H:107
virtual const fileName & name() const
Get the name of the stream.
Definition: OSstream.H:107
unsigned int get(const label i) const
Get value at index i or 0 for out-of-range.
Definition: PackedListI.H:630
static void listCombineAllGather(const List< commsStruct > &comms, List< T > &values, const CombineOp &cop, const int tag, const label comm)
After completion all processors have the same data.
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
Definition: UList.C:212
label rcIndex(const label i) const noexcept
Definition: UListI.H:67
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
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
void set(const bitSet &bitset)
Set specified bits from another bitset.
Definition: bitSetI.H:590
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
Definition: bitSetI.H:521
bitSet & unset(const bitSet &other)
Definition: bitSetI.H:628
Holds information regarding type of cell. Used in inside/outside determination in cellClassification.
Definition: cellInfo.H:65
A collection of cell labels.
Definition: cellSet.H:54
A topoSetFaceSource to select all the faces from given cellSet(s).
Definition: cellToFace.H:177
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
void removeFiles() const
Remove all files from mesh instance()
Definition: faMesh.C:718
A list of face labels.
Definition: faceSet.H:54
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
A class for handling file names.
Definition: fileName.H:76
virtual bool write()
Write the output fields.
Refinement of (split) hexes using polyTopoChange.
Definition: hexRef8.H:68
void checkRefinementLevels(const label maxPointDiff, const labelList &pointsToCheck) const
Debug: Check 2:1 consistency across faces.
Definition: hexRef8.C:4753
labelList consistentSlowRefinement(const label maxFaceDiff, const labelList &cellsToRefine, const labelList &facesToCheck, const label maxPointDiff, const labelList &pointsToCheck) const
Like consistentRefinement but slower:
Definition: hexRef8.C:2306
void checkMesh() const
Debug: Check coupled mesh for correctness.
Definition: hexRef8.C:4544
labelListList setRefinement(const labelList &cells, polyTopoChange &)
Insert refinement. All selected cells will be split into 8.
Definition: hexRef8.C:3190
labelList consistentRefinement(const labelUList &cellLevel, const labelList &cellsToRefine, const bool maxSet) const
Given valid mesh and current cell level and proposed.
Definition: hexRef8.C:2250
void distribute(const mapDistributePolyMesh &)
Update local numbering for mesh redistribution.
Definition: hexRef8.C:4516
const cellShapeList & cellShapes() const
Utility: get hexes as cell shapes.
Definition: hexRef8.C:5005
void setUnrefinement(const labelList &splitPointLabels, polyTopoChange &)
Remove some refinement. Needs to be supplied output of.
Definition: hexRef8.C:5492
labelList consistentUnrefinement(const labelList &pointsToUnrefine, const bool maxSet) const
Given proposed.
Definition: hexRef8.C:5269
labelList getSplitPoints() const
Return the points at the centre of top-level split cells.
Definition: hexRef8.C:5065
label faceLevel(const label facei) const
Gets level such that the face has four points <= level.
Definition: hexRef8.C:797
void storeData(const labelList &pointsToStore, const labelList &facesToStore, const labelList &cellsToStore)
Signal points/face/cells for which to store data.
Definition: hexRef8.C:4204
labelList consistentSlowRefinement2(const label maxFaceDiff, const labelList &cellsToRefine, const labelList &facesToCheck) const
Like consistentSlowRefinement but uses different meshWave.
Definition: hexRef8.C:2790
void subset(const labelList &pointMap, const labelList &faceMap, const labelList &cellMap)
Update local numbering for subsetted mesh.
Definition: hexRef8.C:4433
void setInstance(const fileName &inst)
Definition: hexRef8.C:1732
Class containing mesh-to-mesh mapping information after a mesh distribution where we send parts of me...
void distributeCellData(List< T > &values) const
Distribute list of cell data.
void distributePointData(List< T > &values) const
Distribute list of point data.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:162
label nOldCells() const
Number of old cells.
Definition: mapPolyMesh.H:387
const labelList & reverseCellMap() const
Reverse cell map.
Definition: mapPolyMesh.H:532
const labelList & cellMap() const
Old cell map.
Definition: mapPolyMesh.H:435
const labelList & pointMap() const
Old point map.
Definition: mapPolyMesh.H:396
label nOldPoints() const
Number of old points.
Definition: mapPolyMesh.H:369
const labelList & reversePointMap() const
Reverse point map.
Definition: mapPolyMesh.H:469
void updateMesh()
Update for new mesh topology.
const word & name() const noexcept
The patch name.
A set of point labels.
Definition: pointSet.H:54
Class containing data for cell addition.
Definition: polyAddCell.H:51
Class containing data for point addition.
Definition: polyAddPoint.H:52
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
const fileName & facesInstance() const
Return the current instance directory for faces.
Definition: polyMesh.C:866
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
Definition: polyMesh.H:324
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:380
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:364
Direct mesh changes based on v1.3 polyTopoChange syntax.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
const DynamicList< point > & points() const
Points. Shrunk after constructing mesh (or calling of compact())
label nPoints() const noexcept
Number of mesh points.
label nCells() const noexcept
Number of mesh cells.
Container with cells to refine. Refinement given as single direction.
Definition: refineCell.H:57
Transfers refinement levels such that slow transition between levels is maintained....
label count() const
bool updateFace(const polyMesh &, const label thisFacei, const label neighbourCelli, const refinementData &neighbourInfo, const scalar tol, TrackingData &td)
Influence of neighbouring cell.
Transfers refinement levels such that slow transition between levels is maintained....
All refinement history. Used in unrefinement.
const labelList & visibleCells() const
bool active() const
Is there unrefinement history?
virtual bool read()
Read object. If global number of visible cells > 0 becomes active.
virtual bool write(const bool valid=true) const
Write using setting from DB.
static void syncEdgePositions(const polyMesh &mesh, List< point > &positions, const CombineOp &cop, const point &nullValue)
Synchronize locations on all mesh edges.
Definition: syncTools.H:289
static void syncFaceList(const polyMesh &mesh, UList< T > &faceValues, const CombineOp &cop)
Synchronize values on all mesh faces.
Definition: syncTools.H:396
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 syncBoundaryFacePositions(const polyMesh &mesh, UList< point > &positions, const CombineOp &cop)
Synchronize locations on boundary faces only.
Definition: syncTools.H:378
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.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
const volScalarField & p0
Definition: EEqn.H:36
const polyBoundaryMesh & patches
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const pointField & points
label nPoints
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const labelIOList & zoneID
A HashTable to objects of type <T> with a labelPair key. The hashing is based on labelPair (FixedList...
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
#define WarningInFunction
Report a warning using Foam::Warning.
#define DebugVar(var)
Report a variable name and value.
label findEdge(const edgeList &edges, const labelList &candidates, const label v0, const label v1)
Return edge among candidates that uses the two vertices.
Definition: meshTools.C:359
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of a point.
Definition: meshTools.C:203
Namespace for OpenFOAM.
bool rm(const fileName &file)
Remove a file (or its gz equivalent), returning true if successful.
Definition: MSwindows.C:1012
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
List< cellShape > cellShapeList
List of cellShapes and PtrList of List of cellShape.
Definition: cellShapeList.H:45
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
bool exists(const fileName &name, const bool checkGzip=true, const bool followLink=true)
Does the name exist (as DIRECTORY or FILE) in the file system?
Definition: MSwindows.C:633
List< label > labelList
A List of labels.
Definition: List.H:66
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
constexpr label labelMin
Definition: label.H:60
prefixOSstream Perr
OSstream wrapped stderr (std::cerr) with parallel prefix.
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
Definition: UListI.H:449
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
vector point
Point is a vector.
Definition: point.H:43
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
constexpr label labelMax
Definition: label.H:61
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
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
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
UList< label > labelUList
A UList of labels.
Definition: UList.H:85
Type gMax(const FieldField< Field, Type > &f)
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition: IndirectList.H:68
ListType reorder(const labelUList &oldToNew, const ListType &input, const bool prune=false)
Reorder the elements of a list.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
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 forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278
Reduction class. If x and y are not equal assign value.
Definition: hexRef8.C:60
void operator()(label &x, const label y) const
Definition: hexRef8.C:61