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