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