FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 <OpenFOAM/polyMesh.H>
29 #include "polyTopoChange.H"
30 #include <meshTools/meshTools.H>
35 #include <OpenFOAM/syncTools.H>
36 #include <meshTools/faceSet.H>
37 #include <meshTools/cellSet.H>
38 #include <meshTools/pointSet.H>
39 #include <OpenFOAM/OFstream.H>
40 #include <OpenFOAM/Time.H>
41 #include <OpenFOAM/FaceCellWave.H>
43 #include "refinementData.H"
44 #include "refinementDistanceData.H"
45 
46 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
47 
48 namespace Foam
49 {
50  defineTypeNameAndDebug(hexRef8, 0);
51 
52  //- Reduction class. If x and y are not equal assign value.
53  template<int value>
54  class ifEqEqOp
55  {
56  public:
57  void operator()(label& x, const label& y) const
58  {
59  x = (x==y) ? x : value;
60  }
61  };
62 }
63 
64 
65 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
66 
67 void Foam::hexRef8::reorder
68 (
69  const labelList& map,
70  const label len,
71  const label null,
72  labelList& elems
73 )
74 {
75  labelList newElems(len, null);
76 
77  forAll(elems, i)
78  {
79  label newI = map[i];
80 
81  if (newI >= len)
82  {
83  FatalErrorIn("hexRef8::reorder(..)") << abort(FatalError);
84  }
85 
86  if (newI >= 0)
87  {
88  newElems[newI] = elems[i];
89  }
90  }
91 
92  elems.transfer(newElems);
93 }
94 
95 
96 void Foam::hexRef8::getFaceInfo
97 (
98  const label faceI,
99  label& patchID,
100  label& zoneID,
101  label& zoneFlip
102 ) const
103 {
104  patchID = -1;
105 
106  if (!mesh_.isInternalFace(faceI))
107  {
108  patchID = mesh_.boundaryMesh().whichPatch(faceI);
109  }
110 
111  zoneID = mesh_.faceZones().whichZone(faceI);
112 
113  zoneFlip = false;
114 
115  if (zoneID >= 0)
116  {
117  const faceZone& fZone = mesh_.faceZones()[zoneID];
118 
119  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
120  }
121 }
122 
123 
124 // Adds a face on top of existing faceI.
125 Foam::label Foam::hexRef8::addFace
126 (
127  polyTopoChange& meshMod,
128  const label faceI,
129  const face& newFace,
130  const label own,
131  const label nei
132 ) const
133 {
134  label patchID, zoneID, zoneFlip;
135 
136  getFaceInfo(faceI, patchID, zoneID, zoneFlip);
137 
138  label newFaceI = -1;
139 
140  if ((nei == -1) || (own < nei))
141  {
142  // Ordering ok.
143  newFaceI = meshMod.setAction
144  (
145  polyAddFace
146  (
147  newFace, // face
148  own, // owner
149  nei, // neighbour
150  -1, // master point
151  -1, // master edge
152  faceI, // master face for addition
153  false, // flux flip
154  patchID, // patch for face
155  zoneID, // zone for face
156  zoneFlip // face zone flip
157  )
158  );
159  }
160  else
161  {
162  // Reverse owner/neighbour
163  newFaceI = meshMod.setAction
164  (
165  polyAddFace
166  (
167  newFace.reverseFace(), // face
168  nei, // owner
169  own, // neighbour
170  -1, // master point
171  -1, // master edge
172  faceI, // master face for addition
173  false, // flux flip
174  patchID, // patch for face
175  zoneID, // zone for face
176  zoneFlip // face zone flip
177  )
178  );
179  }
180  return newFaceI;
181 }
182 
183 
184 // Adds an internal face from an edge. Assumes orientation correct.
185 // Problem is that the face is between four new vertices. So what do we provide
186 // as master? The only existing mesh item we have is the edge we have split.
187 // Have to be careful in only using it if it has internal faces since otherwise
188 // polyMeshMorph will complain (because it cannot generate a sensible mapping
189 // for the face)
190 Foam::label Foam::hexRef8::addInternalFace
191 (
192  polyTopoChange& meshMod,
193  const label meshFaceI,
194  const label meshPointI,
195  const face& newFace,
196  const label own,
197  const label nei
198 ) const
199 {
200  if (mesh_.isInternalFace(meshFaceI))
201  {
202  return meshMod.setAction
203  (
204  polyAddFace
205  (
206  newFace, // face
207  own, // owner
208  nei, // neighbour
209  -1, // master point
210  -1, // master edge
211  meshFaceI, // master face for addition
212  false, // flux flip
213  -1, // patch for face
214  -1, // zone for face
215  false // face zone flip
216  )
217  );
218  }
219  else
220  {
221  // Two choices:
222  // - append (i.e. create out of nothing - will not be mapped)
223  // problem: field does not get mapped.
224  // - inflate from point.
225  // problem: does interpolative mapping which constructs full
226  // volPointInterpolation!
227 
228  // For now create out of nothing
229 
230  return meshMod.setAction
231  (
232  polyAddFace
233  (
234  newFace, // face
235  own, // owner
236  nei, // neighbour
237  -1, // master point
238  -1, // master edge
239  -1, // master face for addition
240  false, // flux flip
241  -1, // patch for face
242  -1, // zone for face
243  false // face zone flip
244  )
245  );
246 
247 
250  //label masterPointI = -1;
251  //
252  //const labelList& pFaces = mesh_.pointFaces()[meshPointI];
253  //
254  //forAll(pFaces, i)
255  //{
256  // if (mesh_.isInternalFace(pFaces[i]))
257  // {
258  // // meshPoint uses internal faces so ok to inflate from it
259  // masterPointI = meshPointI;
260  //
261  // break;
262  // }
263  //}
264  //
265  //return meshMod.setAction
266  //(
267  // polyAddFace
268  // (
269  // newFace, // face
270  // own, // owner
271  // nei, // neighbour
272  // masterPointI, // master point
273  // -1, // master edge
274  // -1, // master face for addition
275  // false, // flux flip
276  // -1, // patch for face
277  // -1, // zone for face
278  // false // face zone flip
279  // )
280  //);
281  }
282 }
283 
284 
285 // Modifies existing faceI for either new owner/neighbour or new face points.
286 void Foam::hexRef8::modFace
287 (
288  polyTopoChange& meshMod,
289  const label faceI,
290  const face& newFace,
291  const label own,
292  const label nei
293 ) const
294 {
295  label patchID, zoneID, zoneFlip;
296 
297  getFaceInfo(faceI, patchID, zoneID, zoneFlip);
298 
299  if
300  (
301  (own != mesh_.faceOwner()[faceI])
302  || (
303  mesh_.isInternalFace(faceI)
304  && (nei != mesh_.faceNeighbour()[faceI])
305  )
306  || (newFace != mesh_.faces()[faceI])
307  )
308  {
309  if ((nei == -1) || (own < nei))
310  {
311  meshMod.setAction
312  (
313  polyModifyFace
314  (
315  newFace, // modified face
316  faceI, // label of face being modified
317  own, // owner
318  nei, // neighbour
319  false, // face flip
320  patchID, // patch for face
321  false, // remove from zone
322  zoneID, // zone for face
323  zoneFlip // face flip in zone
324  )
325  );
326  }
327  else
328  {
329  meshMod.setAction
330  (
331  polyModifyFace
332  (
333  newFace.reverseFace(), // modified face
334  faceI, // label of face being modified
335  nei, // owner
336  own, // neighbour
337  false, // face flip
338  patchID, // patch for face
339  false, // remove from zone
340  zoneID, // zone for face
341  zoneFlip // face flip in zone
342  )
343  );
344  }
345  }
346 }
347 
348 
349 // Bit complex way to determine the unrefined edge length.
350 Foam::scalar Foam::hexRef8::getLevel0EdgeLength() const
351 {
352  if (cellLevel_.size() != mesh_.nCells())
353  {
355  (
356  "hexRef8::getLevel0EdgeLength() const"
357  ) << "Number of cells in mesh:" << mesh_.nCells()
358  << " does not equal size of cellLevel:" << cellLevel_.size()
359  << endl
360  << "This might be because of a restart with inconsistent cellLevel."
361  << abort(FatalError);
362  }
363 
364  // Determine minimum edge length per refinement level
365  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
366 
367  const scalar GREAT2 = sqr(GREAT);
368 
369  label nLevels = gMax(cellLevel_)+1;
370 
371  scalarField typEdgeLenSqr(nLevels, GREAT2);
372 
373 
374  // 1. Look only at edges surrounded by cellLevel cells only.
375  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
376 
377  {
378  // Per edge the cellLevel of connected cells. -1 if not set,
379  // labelMax if different levels, otherwise levels of connected cells.
380  labelList edgeLevel(mesh_.nEdges(), -1);
381 
382  forAll(cellLevel_, cellI)
383  {
384  const label cLevel = cellLevel_[cellI];
385 
386  const labelList& cEdges = mesh_.cellEdges(cellI);
387 
388  forAll(cEdges, i)
389  {
390  label edgeI = cEdges[i];
391 
392  if (edgeLevel[edgeI] == -1)
393  {
394  edgeLevel[edgeI] = cLevel;
395  }
396  else if (edgeLevel[edgeI] == labelMax)
397  {
398  // Already marked as on different cellLevels
399  }
400  else if (edgeLevel[edgeI] != cLevel)
401  {
402  edgeLevel[edgeI] = labelMax;
403  }
404  }
405  }
406 
407  // Make sure that edges with different levels on different processors
408  // are also marked. Do the same test (edgeLevel != cLevel) on coupled
409  // edges.
411  (
412  mesh_,
413  edgeLevel,
414  ifEqEqOp<labelMax>(),
415  labelMin,
416  false // no separation
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.
438  Pstream::listCombineGather(typEdgeLenSqr, minEqOp<scalar>());
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 
473  Pstream::listCombineGather(maxEdgeLenSqr, maxEqOp<scalar>());
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.
798 Foam::label Foam::hexRef8::getAnchorLevel(const label faceI) const
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 nessecary 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 = newFace.reverseFace();
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, 1);
1587  }
1588  else
1589  {
1590  refineCell.set(own, 0);
1591  }
1592  nChanged++;
1593  }
1594  else if (neiLevel > (ownLevel+1))
1595  {
1596  if (maxSet)
1597  {
1598  refineCell.set(own, 1);
1599  }
1600  else
1601  {
1602  refineCell.set(nei, 0);
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, false);
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.set(own, 0);
1634  nChanged++;
1635  }
1636  }
1637  else if (neiLevel[i] > (ownLevel+1))
1638  {
1639  if (maxSet)
1640  {
1641  refineCell.set(own, 1);
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], 1);
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, false);
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  history_.instance() = inst;
1750 }
1751 
1752 
1753 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1754 
1755 // Construct from mesh, read refinement data
1756 Foam::hexRef8::hexRef8(const polyMesh& mesh)
1757 :
1758  mesh_(mesh),
1759  cellLevel_
1760  (
1761  IOobject
1762  (
1763  "cellLevel",
1764  mesh_.facesInstance(),
1765  polyMesh::meshSubDir,
1766  mesh_,
1767  IOobject::READ_IF_PRESENT,
1768  IOobject::AUTO_WRITE
1769  ),
1770  labelList(mesh_.nCells(), 0)
1771  ),
1772  pointLevel_
1773  (
1774  IOobject
1775  (
1776  "pointLevel",
1777  mesh_.facesInstance(),
1778  polyMesh::meshSubDir,
1779  mesh_,
1780  IOobject::READ_IF_PRESENT,
1781  IOobject::AUTO_WRITE
1782  ),
1783  labelList(mesh_.nPoints(), 0)
1784  ),
1785  level0Edge_(getLevel0EdgeLength()),
1786  history_
1787  (
1788  IOobject
1789  (
1790  "refinementHistory",
1791  mesh_.facesInstance(),
1792  polyMesh::meshSubDir,
1793  mesh_,
1794  IOobject::READ_IF_PRESENT,
1795  IOobject::AUTO_WRITE
1796  ),
1797  mesh_.nCells() // All cells visible if could not be read
1798  ),
1799  faceRemover_(mesh_, GREAT), // merge boundary faces wherever possible
1800  savedPointLevel_(0),
1801  savedCellLevel_(0)
1802 {
1803  if (history_.active() && history_.visibleCells().size() != mesh_.nCells())
1804  {
1805  FatalErrorIn
1806  (
1807  "hexRef8::hexRef8(const polyMesh&)"
1808  ) << "History enabled but number of visible cells in it "
1809  << history_.visibleCells().size()
1810  << " is not equal to the number of cells in the mesh "
1811  << mesh_.nCells()
1812  << abort(FatalError);
1813  }
1814 
1815  if
1816  (
1817  cellLevel_.size() != mesh_.nCells()
1818  || pointLevel_.size() != mesh_.nPoints()
1819  )
1820  {
1821  FatalErrorIn
1822  (
1823  "hexRef8::hexRef8(const polyMesh&)"
1824  ) << "Restarted from inconsistent cellLevel or pointLevel files."
1825  << endl
1826  << "Number of cells in mesh:" << mesh_.nCells()
1827  << " does not equal size of cellLevel:" << cellLevel_.size() << endl
1828  << "Number of points in mesh:" << mesh_.nPoints()
1829  << " does not equal size of pointLevel:" << pointLevel_.size()
1830  << abort(FatalError);
1831  }
1832 
1833 
1834  // Check refinement levels for consistency
1836 
1837 
1838  // Check initial mesh for consistency
1839 
1840  //if (debug)
1841  {
1842  checkMesh();
1843  }
1844 }
1845 
1846 
1847 // Construct from components
1848 Foam::hexRef8::hexRef8
1850  const polyMesh& mesh,
1851  const labelList& cellLevel,
1852  const labelList& pointLevel,
1853  const refinementHistory& history
1854 )
1855 :
1856  mesh_(mesh),
1857  cellLevel_
1858  (
1859  IOobject
1860  (
1861  "cellLevel",
1862  mesh_.facesInstance(),
1864  mesh_,
1867  ),
1868  cellLevel
1869  ),
1870  pointLevel_
1871  (
1872  IOobject
1873  (
1874  "pointLevel",
1875  mesh_.facesInstance(),
1877  mesh_,
1880  ),
1881  pointLevel
1882  ),
1883  level0Edge_(getLevel0EdgeLength()),
1884  history_
1885  (
1886  IOobject
1887  (
1888  "refinementHistory",
1889  mesh_.facesInstance(),
1891  mesh_,
1894  ),
1895  history
1896  ),
1897  faceRemover_(mesh_, GREAT), // merge boundary faces wherever possible
1898  savedPointLevel_(0),
1899  savedCellLevel_(0)
1900 {
1901  if (history_.active() && history_.visibleCells().size() != mesh_.nCells())
1902  {
1903  FatalErrorIn
1904  (
1905  "hexRef8::hexRef8(const polyMesh&, const labelList&"
1906  ", const labelList&, const refinementHistory&)"
1907  ) << "History enabled but number of visible cells in it "
1908  << history_.visibleCells().size()
1909  << " is not equal to the number of cells in the mesh "
1910  << mesh_.nCells() << abort(FatalError);
1911  }
1912 
1913  if
1914  (
1915  cellLevel_.size() != mesh_.nCells()
1916  || pointLevel_.size() != mesh_.nPoints()
1917  )
1918  {
1919  FatalErrorIn
1920  (
1921  "hexRef8::hexRef8(const polyMesh&, const labelList&"
1922  ", const labelList&, const refinementHistory&)"
1923  ) << "Incorrect cellLevel or pointLevel size." << endl
1924  << "Number of cells in mesh:" << mesh_.nCells()
1925  << " does not equal size of cellLevel:" << cellLevel_.size() << endl
1926  << "Number of points in mesh:" << mesh_.nPoints()
1927  << " does not equal size of pointLevel:" << pointLevel_.size()
1928  << abort(FatalError);
1929  }
1930 
1931  // Check refinement levels for consistency
1932  checkRefinementLevels(-1, labelList(0));
1933 
1934 
1935  // Check initial mesh for consistency
1936 
1937  //if (debug)
1938  {
1939  checkMesh();
1940  }
1941 }
1942 
1943 
1944 // Construct from components
1945 Foam::hexRef8::hexRef8
1947  const polyMesh& mesh,
1948  const labelList& cellLevel,
1949  const labelList& pointLevel
1950 )
1951 :
1952  mesh_(mesh),
1953  cellLevel_
1954  (
1955  IOobject
1956  (
1957  "cellLevel",
1958  mesh_.facesInstance(),
1960  mesh_,
1963  ),
1964  cellLevel
1965  ),
1966  pointLevel_
1967  (
1968  IOobject
1969  (
1970  "pointLevel",
1971  mesh_.facesInstance(),
1973  mesh_,
1976  ),
1977  pointLevel
1978  ),
1979  level0Edge_(getLevel0EdgeLength()),
1980  history_
1981  (
1982  IOobject
1983  (
1984  "refinementHistory",
1985  mesh_.facesInstance(),
1987  mesh_,
1990  ),
1992  labelList(0)
1993  ),
1994  faceRemover_(mesh_, GREAT), // merge boundary faces wherever possible
1995  savedPointLevel_(0),
1996  savedCellLevel_(0)
1997 {
1998  if
1999  (
2000  cellLevel_.size() != mesh_.nCells()
2001  || pointLevel_.size() != mesh_.nPoints()
2002  )
2003  {
2004  FatalErrorIn
2005  (
2006  "hexRef8::hexRef8(const polyMesh&, const labelList&"
2007  ", const labelList&)"
2008  ) << "Incorrect cellLevel or pointLevel size." << endl
2009  << "Number of cells in mesh:" << mesh_.nCells()
2010  << " does not equal size of cellLevel:" << cellLevel_.size() << endl
2011  << "Number of points in mesh:" << mesh_.nPoints()
2012  << " does not equal size of pointLevel:" << pointLevel_.size()
2013  << abort(FatalError);
2014  }
2015 
2016  // Check refinement levels for consistency
2017  checkRefinementLevels(-1, labelList(0));
2018 
2019  // Check initial mesh for consistency
2020 
2021  //if (debug)
2022  {
2023  checkMesh();
2024  }
2025 }
2026 
2027 
2028 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2029 
2032  const labelList& cellsToRefine,
2033  const bool maxSet
2034 ) const
2035 {
2036  // Loop, modifying cellsToRefine, until no more changes to due to 2:1
2037  // conflicts.
2038  // maxSet = false : unselect cells to refine
2039  // maxSet = true : select cells to refine
2040 
2041  // Go to straight boolList.
2042  PackedBoolList refineCell(mesh_.nCells());
2043  forAll(cellsToRefine, i)
2044  {
2045  refineCell.set(cellsToRefine[i], 1);
2046  }
2047 
2048  while (true)
2049  {
2050  label nChanged = faceConsistentRefinement(maxSet, refineCell);
2051 
2052  reduce(nChanged, sumOp<label>());
2053 
2054  if (debug)
2055  {
2056  Pout<< "hexRef8::consistentRefinement : Changed " << nChanged
2057  << " refinement levels due to 2:1 conflicts."
2058  << endl;
2059  }
2060 
2061  if (nChanged == 0)
2062  {
2063  break;
2064  }
2065  }
2066 
2067 
2068  // Convert back to labelList.
2069  label nRefined = 0;
2070 
2071  forAll(refineCell, cellI)
2072  {
2073  if (refineCell.get(cellI) == 1)
2074  {
2075  nRefined++;
2076  }
2077  }
2078 
2079  labelList newCellsToRefine(nRefined);
2080  nRefined = 0;
2081 
2082  forAll(refineCell, cellI)
2083  {
2084  if (refineCell.get(cellI) == 1)
2085  {
2086  newCellsToRefine[nRefined++] = cellI;
2087  }
2088  }
2089 
2090  if (debug)
2091  {
2092  checkWantedRefinementLevels(newCellsToRefine);
2093  }
2094 
2095  return newCellsToRefine;
2096 }
2097 
2098 
2099 // Given a list of cells to refine determine additional cells to refine
2100 // such that the overall refinement:
2101 // - satisfies maxFaceDiff (e.g. 2:1) across neighbouring faces
2102 // - satisfies maxPointDiff (e.g. 4:1) across selected point connected
2103 // cells. This is used to ensure that e.g. cells on the surface are not
2104 // point connected to cells which are 8 times smaller.
2107  const label maxFaceDiff,
2108  const labelList& cellsToRefine,
2109  const labelList& facesToCheck,
2110  const label maxPointDiff,
2111  const labelList& pointsToCheck
2112 ) const
2113 {
2114  const labelList& faceOwner = mesh_.faceOwner();
2115  const labelList& faceNeighbour = mesh_.faceNeighbour();
2116 
2117 
2118  if (maxFaceDiff <= 0)
2119  {
2120  FatalErrorIn
2121  (
2122  "hexRef8::consistentSlowRefinement"
2123  "(const label, const labelList&, const labelList&"
2124  ", const label, const labelList&)"
2125  ) << "Illegal maxFaceDiff " << maxFaceDiff << nl
2126  << "Value should be >= 1" << exit(FatalError);
2127  }
2128 
2129 
2130  // Bit tricky. Say we want a distance of three cells between two
2131  // consecutive refinement levels. This is done by using FaceCellWave to
2132  // transport out the new refinement level. It gets decremented by one
2133  // every cell it crosses so if we initialize it to maxFaceDiff
2134  // we will get a field everywhere that tells us whether an unselected cell
2135  // needs refining as well.
2136 
2137 
2138  // Initial information about (distance to) cellLevel on all cells
2139  List<refinementData> allCellInfo(mesh_.nCells());
2140 
2141  // Initial information about (distance to) cellLevel on all faces
2142  List<refinementData> allFaceInfo(mesh_.nFaces());
2143 
2144  forAll(allCellInfo, cellI)
2145  {
2146  // maxFaceDiff since refinementData counts both
2147  // faces and cells.
2148  allCellInfo[cellI] = refinementData
2149  (
2150  maxFaceDiff*(cellLevel_[cellI]+1),// when cell is to be refined
2151  maxFaceDiff*cellLevel_[cellI] // current level
2152  );
2153  }
2154 
2155  // Cells to be refined will have cellLevel+1
2156  forAll(cellsToRefine, i)
2157  {
2158  label cellI = cellsToRefine[i];
2159 
2160  allCellInfo[cellI].count() = allCellInfo[cellI].refinementCount();
2161  }
2162 
2163 
2164  // Labels of seed faces
2165  DynamicList<label> seedFaces(mesh_.nFaces()/100);
2166  // refinementLevel data on seed faces
2167  DynamicList<refinementData> seedFacesInfo(mesh_.nFaces()/100);
2168 
2169 
2170  // Additional buffer layer thickness by changing initial count. Usually
2171  // this happens on boundary faces. Bit tricky. Use allFaceInfo to mark
2172  // off thus marked faces so they're skipped in the next loop.
2173  forAll(facesToCheck, i)
2174  {
2175  label faceI = facesToCheck[i];
2176 
2177  if (allFaceInfo[faceI].valid())
2178  {
2179  // Can only occur if face has already gone through loop below.
2180  FatalErrorIn
2181  (
2182  "hexRef8::consistentSlowRefinement"
2183  "(const label, const labelList&, const labelList&"
2184  ", const label, const labelList&)"
2185  ) << "Argument facesToCheck seems to have duplicate entries!"
2186  << endl
2187  << "face:" << faceI << " occurs at positions "
2188  << findIndices(facesToCheck, faceI)
2189  << abort(FatalError);
2190  }
2191 
2192 
2193  const refinementData& ownData = allCellInfo[faceOwner[faceI]];
2194 
2195  if (mesh_.isInternalFace(faceI))
2196  {
2197  // Seed face if neighbouring cell (after possible refinement)
2198  // will be refined one more than the current owner or neighbour.
2199 
2200  const refinementData& neiData = allCellInfo[faceNeighbour[faceI]];
2201 
2202  label faceCount;
2203  label faceRefineCount;
2204  if (neiData.count() > ownData.count())
2205  {
2206  faceCount = neiData.count() + maxFaceDiff;
2207  faceRefineCount = faceCount + maxFaceDiff;
2208  }
2209  else
2210  {
2211  faceCount = ownData.count() + maxFaceDiff;
2212  faceRefineCount = faceCount + maxFaceDiff;
2213  }
2214 
2215  seedFaces.append(faceI);
2216  seedFacesInfo.append
2217  (
2219  (
2220  faceRefineCount,
2221  faceCount
2222  )
2223  );
2224  allFaceInfo[faceI] = seedFacesInfo[seedFacesInfo.size()-1];
2225  }
2226  else
2227  {
2228  label faceCount = ownData.count() + maxFaceDiff;
2229  label faceRefineCount = faceCount + maxFaceDiff;
2230 
2231  seedFaces.append(faceI);
2232  seedFacesInfo.append
2233  (
2235  (
2236  faceRefineCount,
2237  faceCount
2238  )
2239  );
2240  allFaceInfo[faceI] = seedFacesInfo[seedFacesInfo.size()-1];
2241  }
2242  }
2243 
2244 
2245  // Just seed with all faces inbetween different refinement levels for now
2246  // (alternatively only seed faces on cellsToRefine but that gives problems
2247  // if no cells to refine)
2248  forAll(faceNeighbour, faceI)
2249  {
2250  // Check if face already handled in loop above
2251  if (!allFaceInfo[faceI].valid())
2252  {
2253  label own = faceOwner[faceI];
2254  label nei = faceNeighbour[faceI];
2255 
2256  // Seed face with transported data from highest cell.
2257 
2258  if (allCellInfo[own].count() > allCellInfo[nei].count())
2259  {
2260  allFaceInfo[faceI].updateFace
2261  (
2262  mesh_,
2263  faceI,
2264  own,
2265  allCellInfo[own],
2267  );
2268  seedFaces.append(faceI);
2269  seedFacesInfo.append(allFaceInfo[faceI]);
2270  }
2271  else if (allCellInfo[own].count() < allCellInfo[nei].count())
2272  {
2273  allFaceInfo[faceI].updateFace
2274  (
2275  mesh_,
2276  faceI,
2277  nei,
2278  allCellInfo[nei],
2280  );
2281  seedFaces.append(faceI);
2282  seedFacesInfo.append(allFaceInfo[faceI]);
2283  }
2284  }
2285  }
2286 
2287  // Seed all boundary faces with owner value. This is to make sure that
2288  // they are visited (probably only important for coupled faces since
2289  // these need to be visited from both sides)
2290  for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
2291  {
2292  // Check if face already handled in loop above
2293  if (!allFaceInfo[faceI].valid())
2294  {
2295  label own = faceOwner[faceI];
2296 
2297  // Seed face with transported data from owner.
2298  refinementData faceData;
2299  faceData.updateFace
2300  (
2301  mesh_,
2302  faceI,
2303  own,
2304  allCellInfo[own],
2306  );
2307  seedFaces.append(faceI);
2308  seedFacesInfo.append(faceData);
2309  }
2310  }
2311 
2312 
2313  // face-cell-face transport engine
2315  (
2316  mesh_,
2317  allFaceInfo,
2318  allCellInfo
2319  );
2320 
2321  while(true)
2322  {
2323  if (debug)
2324  {
2325  Pout<< "hexRef8::consistentSlowRefinement : Seeded "
2326  << seedFaces.size() << " faces between cells with different"
2327  << " refinement level." << endl;
2328  }
2329 
2330  // Set seed faces
2331  levelCalc.setFaceInfo(seedFaces.shrink(), seedFacesInfo.shrink());
2332  seedFaces.clear();
2333  seedFacesInfo.clear();
2334 
2335  // Iterate until no change. Now 2:1 face difference should be satisfied
2336  levelCalc.iterate(mesh_.globalData().nTotalFaces());
2337 
2338 
2339  // Now check point-connected cells (face-connected cells already ok):
2340  // - get per point max of connected cells
2341  // - sync across coupled points
2342  // - check cells against above point max
2343 
2344  if (maxPointDiff == -1)
2345  {
2346  // No need to do any point checking.
2347  break;
2348  }
2349 
2350  // Determine per point the max cell level. (done as count, not
2351  // as cell level purely for ease)
2352  labelList maxPointCount(mesh_.nPoints(), 0);
2353 
2354  forAll(maxPointCount, pointI)
2355  {
2356  label& pLevel = maxPointCount[pointI];
2357 
2358  const labelList& pCells = mesh_.pointCells(pointI);
2359 
2360  forAll(pCells, i)
2361  {
2362  pLevel = max(pLevel, allCellInfo[pCells[i]].count());
2363  }
2364  }
2365 
2366  // Sync maxPointCount to neighbour
2368  (
2369  mesh_,
2370  maxPointCount,
2371  maxEqOp<label>(),
2372  labelMin, // null value
2373  false // no separation
2374  );
2375 
2376  // Update allFaceInfo from maxPointCount for all points to check
2377  // (usually on boundary faces)
2378 
2379  // Per face the new refinement data
2380  Map<refinementData> changedFacesInfo(pointsToCheck.size());
2381 
2382  forAll(pointsToCheck, i)
2383  {
2384  label pointI = pointsToCheck[i];
2385 
2386  // Loop over all cells using the point and check whether their
2387  // refinement level is much less than the maximum.
2388 
2389  const labelList& pCells = mesh_.pointCells(pointI);
2390 
2391  forAll(pCells, pCellI)
2392  {
2393  label cellI = pCells[pCellI];
2394 
2395  refinementData& cellInfo = allCellInfo[cellI];
2396 
2397  if
2398  (
2399  !cellInfo.isRefined()
2400  && (
2401  maxPointCount[pointI]
2402  > cellInfo.count() + maxFaceDiff*maxPointDiff
2403  )
2404  )
2405  {
2406  // Mark cell for refinement
2407  cellInfo.count() = cellInfo.refinementCount();
2408 
2409  // Insert faces of cell as seed faces.
2410  const cell& cFaces = mesh_.cells()[cellI];
2411 
2412  forAll(cFaces, cFaceI)
2413  {
2414  label faceI = cFaces[cFaceI];
2415 
2416  refinementData faceData;
2417  faceData.updateFace
2418  (
2419  mesh_,
2420  faceI,
2421  cellI,
2422  cellInfo,
2424  );
2425 
2426  if (faceData.count() > allFaceInfo[faceI].count())
2427  {
2428  changedFacesInfo.insert(faceI, faceData);
2429  }
2430  }
2431  }
2432  }
2433  }
2434 
2435  label nChanged = changedFacesInfo.size();
2436  reduce(nChanged, sumOp<label>());
2437 
2438  if (nChanged == 0)
2439  {
2440  break;
2441  }
2442 
2443 
2444  // Transfer into seedFaces, seedFacesInfo
2445  seedFaces.setCapacity(changedFacesInfo.size());
2446  seedFacesInfo.setCapacity(changedFacesInfo.size());
2447 
2448  forAllConstIter(Map<refinementData>, changedFacesInfo, iter)
2449  {
2450  seedFaces.append(iter.key());
2451  seedFacesInfo.append(iter());
2452  }
2453  }
2454 
2455 
2456  if (debug)
2457  {
2458  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
2459  {
2460  label own = mesh_.faceOwner()[faceI];
2461  label ownLevel =
2462  cellLevel_[own]
2463  + (allCellInfo[own].isRefined() ? 1 : 0);
2464 
2465  label nei = mesh_.faceNeighbour()[faceI];
2466  label neiLevel =
2467  cellLevel_[nei]
2468  + (allCellInfo[nei].isRefined() ? 1 : 0);
2469 
2470  if (mag(ownLevel-neiLevel) > 1)
2471  {
2472  dumpCell(own);
2473  dumpCell(nei);
2474  FatalErrorIn
2475  (
2476  "hexRef8::consistentSlowRefinement"
2477  "(const label, const labelList&, const labelList&"
2478  ", const label, const labelList&)"
2479  ) << "cell:" << own
2480  << " current level:" << cellLevel_[own]
2481  << " current refData:" << allCellInfo[own]
2482  << " level after refinement:" << ownLevel
2483  << nl
2484  << "neighbour cell:" << nei
2485  << " current level:" << cellLevel_[nei]
2486  << " current refData:" << allCellInfo[nei]
2487  << " level after refinement:" << neiLevel
2488  << nl
2489  << "which does not satisfy 2:1 constraints anymore." << nl
2490  << "face:" << faceI << " faceRefData:" << allFaceInfo[faceI]
2491  << abort(FatalError);
2492  }
2493  }
2494 
2495 
2496  // Coupled faces. Swap owner level to get neighbouring cell level.
2497  // (only boundary faces of neiLevel used)
2498 
2499  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2500  labelList neiCount(mesh_.nFaces()-mesh_.nInternalFaces());
2501  labelList neiRefCount(mesh_.nFaces()-mesh_.nInternalFaces());
2502 
2503  forAll(neiLevel, i)
2504  {
2505  label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
2506  neiLevel[i] = cellLevel_[own];
2507  neiCount[i] = allCellInfo[own].count();
2508  neiRefCount[i] = allCellInfo[own].refinementCount();
2509  }
2510 
2511  // Swap to neighbour
2512  syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
2513  syncTools::swapBoundaryFaceList(mesh_, neiCount, false);
2514  syncTools::swapBoundaryFaceList(mesh_, neiRefCount, false);
2515 
2516  // Now we have neighbour value see which cells need refinement
2517  forAll(neiLevel, i)
2518  {
2519  label faceI = i+mesh_.nInternalFaces();
2520 
2521  label own = mesh_.faceOwner()[faceI];
2522  label ownLevel =
2523  cellLevel_[own]
2524  + (allCellInfo[own].isRefined() ? 1 : 0);
2525 
2526  label nbrLevel =
2527  neiLevel[i]
2528  + ((neiCount[i] >= neiRefCount[i]) ? 1 : 0);
2529 
2530  if (mag(ownLevel - nbrLevel) > 1)
2531  {
2532  dumpCell(own);
2533  label patchI = mesh_.boundaryMesh().whichPatch(faceI);
2534 
2535  FatalErrorIn
2536  (
2537  "hexRef8::consistentSlowRefinement"
2538  "(const label, const labelList&, const labelList&"
2539  ", const label, const labelList&)"
2540  ) << "Celllevel does not satisfy 2:1 constraint."
2541  << " On coupled face "
2542  << faceI
2543  << " refData:" << allFaceInfo[faceI]
2544  << " on patch " << patchI << " "
2545  << mesh_.boundaryMesh()[patchI].name() << nl
2546  << "owner cell " << own
2547  << " current level:" << cellLevel_[own]
2548  << " current count:" << allCellInfo[own].count()
2549  << " current refCount:"
2550  << allCellInfo[own].refinementCount()
2551  << " level after refinement:" << ownLevel
2552  << nl
2553  << "(coupled) neighbour cell"
2554  << " has current level:" << neiLevel[i]
2555  << " current count:" << neiCount[i]
2556  << " current refCount:" << neiRefCount[i]
2557  << " level after refinement:" << nbrLevel
2558  << abort(FatalError);
2559  }
2560  }
2561  }
2562 
2563  // Convert back to labelList of cells to refine.
2564 
2565  label nRefined = 0;
2566 
2567  forAll(allCellInfo, cellI)
2568  {
2569  if (allCellInfo[cellI].isRefined())
2570  {
2571  nRefined++;
2572  }
2573  }
2574 
2575  // Updated list of cells to refine
2576  labelList newCellsToRefine(nRefined);
2577  nRefined = 0;
2578 
2579  forAll(allCellInfo, cellI)
2580  {
2581  if (allCellInfo[cellI].isRefined())
2582  {
2583  newCellsToRefine[nRefined++] = cellI;
2584  }
2585  }
2586 
2587  if (debug)
2588  {
2589  Pout<< "hexRef8::consistentSlowRefinement : From "
2590  << cellsToRefine.size() << " to " << newCellsToRefine.size()
2591  << " cells to refine." << endl;
2592  }
2593 
2594  return newCellsToRefine;
2595 }
2596 
2597 
2600  const label maxFaceDiff,
2601  const labelList& cellsToRefine,
2602  const labelList& facesToCheck
2603 ) const
2604 {
2605  const labelList& faceOwner = mesh_.faceOwner();
2606  const labelList& faceNeighbour = mesh_.faceNeighbour();
2607 
2608  if (maxFaceDiff <= 0)
2609  {
2610  FatalErrorIn
2611  (
2612  "hexRef8::consistentSlowRefinement2"
2613  "(const label, const labelList&, const labelList&)"
2614  ) << "Illegal maxFaceDiff " << maxFaceDiff << nl
2615  << "Value should be >= 1" << exit(FatalError);
2616  }
2617 
2618  const scalar level0Size = 2*maxFaceDiff*level0Edge_;
2619 
2620 
2621  // Bit tricky. Say we want a distance of three cells between two
2622  // consecutive refinement levels. This is done by using FaceCellWave to
2623  // transport out the 'refinement shell'. Anything inside the refinement
2624  // shell (given by a distance) gets marked for refinement.
2625 
2626  // Initial information about (distance to) cellLevel on all cells
2627  List<refinementDistanceData> allCellInfo(mesh_.nCells());
2628 
2629  // Initial information about (distance to) cellLevel on all faces
2630  List<refinementDistanceData> allFaceInfo(mesh_.nFaces());
2631 
2632 
2633  // Mark cells with wanted refinement level
2634  forAll(cellsToRefine, i)
2635  {
2636  label cellI = cellsToRefine[i];
2637 
2638  allCellInfo[cellI] = refinementDistanceData
2639  (
2640  level0Size,
2641  mesh_.cellCentres()[cellI],
2642  cellLevel_[cellI]+1 // wanted refinement
2643  );
2644  }
2645  // Mark all others with existing refinement level
2646  forAll(allCellInfo, cellI)
2647  {
2648  if (!allCellInfo[cellI].valid())
2649  {
2650  allCellInfo[cellI] = refinementDistanceData
2651  (
2652  level0Size,
2653  mesh_.cellCentres()[cellI],
2654  cellLevel_[cellI] // wanted refinement
2655  );
2656  }
2657  }
2658 
2659 
2660  // Labels of seed faces
2661  DynamicList<label> seedFaces(mesh_.nFaces()/100);
2662  // refinementLevel data on seed faces
2663  DynamicList<refinementDistanceData> seedFacesInfo(mesh_.nFaces()/100);
2664 
2665 
2666  const pointField& cc = mesh_.cellCentres();
2667 
2668  forAll(facesToCheck, i)
2669  {
2670  label faceI = facesToCheck[i];
2671 
2672  if (allFaceInfo[faceI].valid())
2673  {
2674  // Can only occur if face has already gone through loop below.
2675  FatalErrorIn
2676  (
2677  "hexRef8::consistentSlowRefinement2"
2678  "(const label, const labelList&, const labelList&)"
2679  ) << "Argument facesToCheck seems to have duplicate entries!"
2680  << endl
2681  << "face:" << faceI << " occurs at positions "
2682  << findIndices(facesToCheck, faceI)
2683  << abort(FatalError);
2684  }
2685 
2686  label own = faceOwner[faceI];
2687 
2688  label ownLevel =
2689  (
2690  allCellInfo[own].valid()
2691  ? allCellInfo[own].originLevel()
2692  : cellLevel_[own]
2693  );
2694 
2695  if (!mesh_.isInternalFace(faceI))
2696  {
2697  // Do as if boundary face would have neighbour with one higher
2698  // refinement level.
2699  const point& fc = mesh_.faceCentres()[faceI];
2700 
2701  refinementDistanceData neiData
2702  (
2703  level0Size,
2704  2*fc - cc[own], // est'd cell centre
2705  ownLevel+1
2706  );
2707 
2708  allFaceInfo[faceI].updateFace
2709  (
2710  mesh_,
2711  faceI,
2712  own, // not used (should be nei)
2713  neiData,
2715  );
2716  }
2717  else
2718  {
2719  label nei = faceNeighbour[faceI];
2720 
2721  label neiLevel =
2722  (
2723  allCellInfo[nei].valid()
2724  ? allCellInfo[nei].originLevel()
2725  : cellLevel_[nei]
2726  );
2727 
2728  if (ownLevel == neiLevel)
2729  {
2730  // Fake as if nei>own or own>nei (whichever one 'wins')
2731  allFaceInfo[faceI].updateFace
2732  (
2733  mesh_,
2734  faceI,
2735  nei,
2736  refinementDistanceData(level0Size, cc[nei], neiLevel+1),
2738  );
2739  allFaceInfo[faceI].updateFace
2740  (
2741  mesh_,
2742  faceI,
2743  own,
2744  refinementDistanceData(level0Size, cc[own], ownLevel+1),
2746  );
2747  }
2748  else
2749  {
2750  // Difference in level anyway.
2751  allFaceInfo[faceI].updateFace
2752  (
2753  mesh_,
2754  faceI,
2755  nei,
2756  refinementDistanceData(level0Size, cc[nei], neiLevel),
2758  );
2759  allFaceInfo[faceI].updateFace
2760  (
2761  mesh_,
2762  faceI,
2763  own,
2764  refinementDistanceData(level0Size, cc[own], ownLevel),
2766  );
2767  }
2768  }
2769  seedFaces.append(faceI);
2770  seedFacesInfo.append(allFaceInfo[faceI]);
2771  }
2772 
2773 
2774  // Create some initial seeds to start walking from. This is only if there
2775  // are no facesToCheck.
2776  // Just seed with all faces inbetween different refinement levels for now
2777  forAll(faceNeighbour, faceI)
2778  {
2779  // Check if face already handled in loop above
2780  if (!allFaceInfo[faceI].valid())
2781  {
2782  label own = faceOwner[faceI];
2783 
2784  label ownLevel =
2785  (
2786  allCellInfo[own].valid()
2787  ? allCellInfo[own].originLevel()
2788  : cellLevel_[own]
2789  );
2790 
2791  label nei = faceNeighbour[faceI];
2792 
2793  label neiLevel =
2794  (
2795  allCellInfo[nei].valid()
2796  ? allCellInfo[nei].originLevel()
2797  : cellLevel_[nei]
2798  );
2799 
2800  if (ownLevel > neiLevel)
2801  {
2802  // Set face to owner data. (since face not yet would be copy)
2803  seedFaces.append(faceI);
2804  allFaceInfo[faceI].updateFace
2805  (
2806  mesh_,
2807  faceI,
2808  own,
2809  refinementDistanceData(level0Size, cc[own], ownLevel),
2811  );
2812  seedFacesInfo.append(allFaceInfo[faceI]);
2813  }
2814  else if (neiLevel > ownLevel)
2815  {
2816  seedFaces.append(faceI);
2817  allFaceInfo[faceI].updateFace
2818  (
2819  mesh_,
2820  faceI,
2821  nei,
2822  refinementDistanceData(level0Size, cc[nei], neiLevel),
2824  );
2825  seedFacesInfo.append(allFaceInfo[faceI]);
2826  }
2827  }
2828  }
2829 
2830  seedFaces.shrink();
2831  seedFacesInfo.shrink();
2832 
2833  // face-cell-face transport engine
2835  (
2836  mesh_,
2837  seedFaces,
2838  seedFacesInfo,
2839  allFaceInfo,
2840  allCellInfo,
2841  mesh_.globalData().nTotalCells()
2842  );
2843 
2844 
2845  //if (debug)
2846  //{
2847  // // Dump wanted level
2848  // volScalarField wantedLevel
2849  // (
2850  // IOobject
2851  // (
2852  // "wantedLevel",
2853  // fMesh.time().timeName(),
2854  // fMesh,
2855  // IOobject::NO_READ,
2856  // IOobject::AUTO_WRITE,
2857  // false
2858  // ),
2859  // fMesh,
2860  // dimensionedScalar("zero", dimless, 0)
2861  // );
2862  //
2863  // forAll(wantedLevel, cellI)
2864  // {
2865  // wantedLevel[cellI] = allCellInfo[cellI].wantedLevel(cc[cellI]);
2866  // }
2867  //
2868  // Pout<< "Writing " << wantedLevel.objectPath() << endl;
2869  // wantedLevel.write();
2870  //}
2871 
2872 
2873  // Convert back to labelList of cells to refine.
2874  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2875 
2876  // 1. Force original refinement cells to be picked up by setting the
2877  // originLevel of input cells to be a very large level (but within range
2878  // of 1<< shift inside refinementDistanceData::wantedLevel)
2879  forAll(cellsToRefine, i)
2880  {
2881  label cellI = cellsToRefine[i];
2882 
2883  allCellInfo[cellI].originLevel() = sizeof(label)*8-2;
2884  allCellInfo[cellI].origin() = cc[cellI];
2885  }
2886 
2887  // 2. Extend to 2:1. I don't understand yet why this is not done
2888  // 2. Extend to 2:1. For non-cube cells the scalar distance does not work
2889  // so make sure it at least provides 2:1.
2890  PackedBoolList refineCell(mesh_.nCells());
2891  forAll(allCellInfo, cellI)
2892  {
2893  label wanted = allCellInfo[cellI].wantedLevel(cc[cellI]);
2894 
2895  if (wanted > cellLevel_[cellI]+1)
2896  {
2897  refineCell.set(cellI, 1);
2898  }
2899  }
2900  faceConsistentRefinement(true, refineCell);
2901 
2902  while (true)
2903  {
2904  label nChanged = faceConsistentRefinement(true, refineCell);
2905 
2906  reduce(nChanged, sumOp<label>());
2907 
2908  if (debug)
2909  {
2910  Pout<< "hexRef8::consistentSlowRefinement2 : Changed " << nChanged
2911  << " refinement levels due to 2:1 conflicts."
2912  << endl;
2913  }
2914 
2915  if (nChanged == 0)
2916  {
2917  break;
2918  }
2919  }
2920 
2921  // 3. Convert back to labelList.
2922  label nRefined = 0;
2923 
2924  forAll(refineCell, cellI)
2925  {
2926  if (refineCell.get(cellI) == 1)
2927  {
2928  nRefined++;
2929  }
2930  }
2931 
2932  labelList newCellsToRefine(nRefined);
2933  nRefined = 0;
2934 
2935  forAll(refineCell, cellI)
2936  {
2937  if (refineCell.get(cellI) == 1)
2938  {
2939  newCellsToRefine[nRefined++] = cellI;
2940  }
2941  }
2942 
2943  if (debug)
2944  {
2945  Pout<< "hexRef8::consistentSlowRefinement2 : From "
2946  << cellsToRefine.size() << " to " << newCellsToRefine.size()
2947  << " cells to refine." << endl;
2948 
2949  // Check that newCellsToRefine obeys at least 2:1.
2950 
2951  {
2952  cellSet cellsIn(mesh_, "cellsToRefineIn", cellsToRefine);
2953  Pout<< "hexRef8::consistentSlowRefinement2 : writing "
2954  << cellsIn.size() << " to cellSet "
2955  << cellsIn.objectPath() << endl;
2956  cellsIn.write();
2957  }
2958  {
2959  cellSet cellsOut(mesh_, "cellsToRefineOut", newCellsToRefine);
2960  Pout<< "hexRef8::consistentSlowRefinement2 : writing "
2961  << cellsOut.size() << " to cellSet "
2962  << cellsOut.objectPath() << endl;
2963  cellsOut.write();
2964  }
2965 
2966  // Extend to 2:1
2967  PackedBoolList refineCell(mesh_.nCells());
2968  forAll(newCellsToRefine, i)
2969  {
2970  refineCell.set(newCellsToRefine[i], 1);
2971  }
2972  const PackedBoolList savedRefineCell(refineCell);
2973 
2974  label nChanged = faceConsistentRefinement(true, refineCell);
2975 
2976  {
2977  cellSet cellsOut2
2978  (
2979  mesh_, "cellsToRefineOut2", newCellsToRefine.size()
2980  );
2981  forAll(refineCell, cellI)
2982  {
2983  if (refineCell.get(cellI) == 1)
2984  {
2985  cellsOut2.insert(cellI);
2986  }
2987  }
2988  Pout<< "hexRef8::consistentSlowRefinement2 : writing "
2989  << cellsOut2.size() << " to cellSet "
2990  << cellsOut2.objectPath() << endl;
2991  cellsOut2.write();
2992  }
2993 
2994  if (nChanged > 0)
2995  {
2996  forAll(refineCell, cellI)
2997  {
2998  if
2999  (
3000  refineCell.get(cellI) == 1
3001  && savedRefineCell.get(cellI) == 0
3002  )
3003  {
3004  dumpCell(cellI);
3005  FatalErrorIn
3006  (
3007  "hexRef8::consistentSlowRefinement2"
3008  "(const label, const labelList&, const labelList&)"
3009  ) << "Cell:" << cellI << " cc:"
3010  << mesh_.cellCentres()[cellI]
3011  << " was not marked for refinement but does not obey"
3012  << " 2:1 constraints."
3013  << abort(FatalError);
3014  }
3015  }
3016  }
3017  }
3018 
3019  return newCellsToRefine;
3020 }
3021 
3022 
3023 // Top level driver to insert topo changes to do all refinement.
3026  const labelList& cellLabels,
3027  polyTopoChange& meshMod
3028 )
3029 {
3030  if (debug)
3031  {
3032  Pout<< "hexRef8::setRefinement :"
3033  << " Checking initial mesh just to make sure" << endl;
3034 
3035  checkMesh();
3036  // Cannot call checkRefinementlevels since hanging points might
3037  // get triggered by the mesher after subsetting.
3038  //checkRefinementLevels(-1, labelList(0));
3039  }
3040 
3041  // Clear any saved point/cell data.
3042  savedPointLevel_.clear();
3043  savedCellLevel_.clear();
3044 
3045 
3046  // New point/cell level. Copy of pointLevel for existing points.
3047  DynamicList<label> newCellLevel(cellLevel_.size());
3048  forAll(cellLevel_, cellI)
3049  {
3050  newCellLevel.append(cellLevel_[cellI]);
3051  }
3052  DynamicList<label> newPointLevel(pointLevel_.size());
3053  forAll(pointLevel_, pointI)
3054  {
3055  newPointLevel.append(pointLevel_[pointI]);
3056  }
3057 
3058 
3059  if (debug)
3060  {
3061  Pout<< "hexRef8::setRefinement :"
3062  << " Allocating " << cellLabels.size() << " cell midpoints."
3063  << endl;
3064  }
3065 
3066 
3067  // Mid point per refined cell.
3068  // -1 : not refined
3069  // >=0: label of mid point.
3070  labelList cellMidPoint(mesh_.nCells(), -1);
3071 
3072  forAll(cellLabels, i)
3073  {
3074  label cellI = cellLabels[i];
3075 
3076  label anchorPointI = mesh_.faces()[mesh_.cells()[cellI][0]][0];
3077 
3078  cellMidPoint[cellI] = meshMod.setAction
3079  (
3080  polyAddPoint
3081  (
3082  mesh_.cellCentres()[cellI], // point
3083  anchorPointI, // master point
3084  -1, // zone for point
3085  true // supports a cell
3086  )
3087  );
3088 
3089  newPointLevel(cellMidPoint[cellI]) = cellLevel_[cellI]+1;
3090  }
3091 
3092 
3093  if (debug)
3094  {
3095  cellSet splitCells(mesh_, "splitCells", cellLabels.size());
3096 
3097  forAll(cellMidPoint, cellI)
3098  {
3099  if (cellMidPoint[cellI] >= 0)
3100  {
3101  splitCells.insert(cellI);
3102  }
3103  }
3104 
3105  Pout<< "hexRef8::setRefinement : Dumping " << splitCells.size()
3106  << " cells to split to cellSet " << splitCells.objectPath()
3107  << endl;
3108 
3109  splitCells.write();
3110  }
3111 
3112 
3113 
3114  // Split edges
3115  // ~~~~~~~~~~~
3116 
3117  if (debug)
3118  {
3119  Pout<< "hexRef8::setRefinement :"
3120  << " Allocating edge midpoints."
3121  << endl;
3122  }
3123 
3124  // Unrefined edges are ones between cellLevel or lower points.
3125  // If any cell using this edge gets split then the edge needs to be split.
3126 
3127  // -1 : no need to split edge
3128  // >=0 : label of introduced mid point
3129  labelList edgeMidPoint(mesh_.nEdges(), -1);
3130 
3131  // Note: Loop over cells to be refined or edges?
3132  forAll(cellMidPoint, cellI)
3133  {
3134  if (cellMidPoint[cellI] >= 0)
3135  {
3136  const labelList& cEdges = mesh_.cellEdges(cellI);
3137 
3138  forAll(cEdges, i)
3139  {
3140  label edgeI = cEdges[i];
3141 
3142  const edge& e = mesh_.edges()[edgeI];
3143 
3144  if
3145  (
3146  pointLevel_[e[0]] <= cellLevel_[cellI]
3147  && pointLevel_[e[1]] <= cellLevel_[cellI]
3148  )
3149  {
3150  edgeMidPoint[edgeI] = 12345; // mark need for splitting
3151  }
3152  }
3153  }
3154  }
3155 
3156  // Synchronize edgeMidPoint across coupled patches. Take max so that
3157  // any split takes precedence.
3159  (
3160  mesh_,
3161  edgeMidPoint,
3162  maxEqOp<label>(),
3163  labelMin,
3164  false // no separation
3165  );
3166 
3167 
3168  // Introduce edge points
3169  // ~~~~~~~~~~~~~~~~~~~~~
3170 
3171  {
3172  // Phase 1: calculate midpoints and sync.
3173  // This needs doing for if people do not write binary and we slowly
3174  // get differences.
3175 
3176  pointField edgeMids(mesh_.nEdges(), point(-GREAT, -GREAT, -GREAT));
3177 
3178  forAll(edgeMidPoint, edgeI)
3179  {
3180  if (edgeMidPoint[edgeI] >= 0)
3181  {
3182  // Edge marked to be split.
3183  edgeMids[edgeI] = mesh_.edges()[edgeI].centre(mesh_.points());
3184  }
3185  }
3187  (
3188  mesh_,
3189  edgeMids,
3190  maxEqOp<vector>(),
3191  point(-GREAT, -GREAT, -GREAT),
3192  true // apply separation
3193  );
3194 
3195 
3196  // Phase 2: introduce points at the synced locations.
3197  forAll(edgeMidPoint, edgeI)
3198  {
3199  if (edgeMidPoint[edgeI] >= 0)
3200  {
3201  // Edge marked to be split. Replace edgeMidPoint with actual
3202  // point label.
3203 
3204  const edge& e = mesh_.edges()[edgeI];
3205 
3206  edgeMidPoint[edgeI] = meshMod.setAction
3207  (
3208  polyAddPoint
3209  (
3210  edgeMids[edgeI], // point
3211  e[0], // master point
3212  -1, // zone for point
3213  true // supports a cell
3214  )
3215  );
3216 
3217  newPointLevel(edgeMidPoint[edgeI]) =
3218  max
3219  (
3220  pointLevel_[e[0]],
3221  pointLevel_[e[1]]
3222  )
3223  + 1;
3224  }
3225  }
3226  }
3227 
3228  if (debug)
3229  {
3230  OFstream str(mesh_.time().path()/"edgeMidPoint.obj");
3231 
3232  forAll(edgeMidPoint, edgeI)
3233  {
3234  if (edgeMidPoint[edgeI] >= 0)
3235  {
3236  const edge& e = mesh_.edges()[edgeI];
3237 
3238  meshTools::writeOBJ(str, e.centre(mesh_.points()));
3239  }
3240  }
3241 
3242  Pout<< "hexRef8::setRefinement :"
3243  << " Dumping edge centres to split to file " << str.name() << endl;
3244  }
3245 
3246 
3247  // Calculate face level
3248  // ~~~~~~~~~~~~~~~~~~~~
3249  // (after splitting)
3250 
3251  if (debug)
3252  {
3253  Pout<< "hexRef8::setRefinement :"
3254  << " Allocating face midpoints."
3255  << endl;
3256  }
3257 
3258  // Face anchor level. There are guaranteed 4 points with level
3259  // <= anchorLevel. These are the corner points.
3260  labelList faceAnchorLevel(mesh_.nFaces());
3261 
3262  for (label faceI = 0; faceI < mesh_.nFaces(); faceI++)
3263  {
3264  faceAnchorLevel[faceI] = getAnchorLevel(faceI);
3265  }
3266 
3267  // -1 : no need to split face
3268  // >=0 : label of introduced mid point
3269  labelList faceMidPoint(mesh_.nFaces(), -1);
3270 
3271 
3272  // Internal faces: look at cells on both sides. Uniquely determined since
3273  // face itself guaranteed to be same level as most refined neighbour.
3274  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
3275  {
3276  if (faceAnchorLevel[faceI] >= 0)
3277  {
3278  label own = mesh_.faceOwner()[faceI];
3279  label ownLevel = cellLevel_[own];
3280  label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3281 
3282  label nei = mesh_.faceNeighbour()[faceI];
3283  label neiLevel = cellLevel_[nei];
3284  label newNeiLevel = neiLevel + (cellMidPoint[nei] >= 0 ? 1 : 0);
3285 
3286  if
3287  (
3288  newOwnLevel > faceAnchorLevel[faceI]
3289  || newNeiLevel > faceAnchorLevel[faceI]
3290  )
3291  {
3292  faceMidPoint[faceI] = 12345; // mark to be split
3293  }
3294  }
3295  }
3296 
3297  // Coupled patches handled like internal faces except now all information
3298  // from neighbour comes from across processor.
3299  // Boundary faces are more complicated since the boundary face can
3300  // be more refined than its owner (or neighbour for coupled patches)
3301  // (does not happen if refining/unrefining only, but does e.g. when
3302  // refinining and subsetting)
3303 
3304  {
3305  labelList newNeiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
3306 
3307  forAll(newNeiLevel, i)
3308  {
3309  label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
3310  label ownLevel = cellLevel_[own];
3311  label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3312 
3313  newNeiLevel[i] = newOwnLevel;
3314  }
3315 
3316  // Swap.
3317  syncTools::swapBoundaryFaceList(mesh_, newNeiLevel, false);
3318 
3319  // So now we have information on the neighbour.
3320 
3321  forAll(newNeiLevel, i)
3322  {
3323  label faceI = i+mesh_.nInternalFaces();
3324 
3325  if (faceAnchorLevel[faceI] >= 0)
3326  {
3327  label own = mesh_.faceOwner()[faceI];
3328  label ownLevel = cellLevel_[own];
3329  label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3330 
3331  if
3332  (
3333  newOwnLevel > faceAnchorLevel[faceI]
3334  || newNeiLevel[i] > faceAnchorLevel[faceI]
3335  )
3336  {
3337  faceMidPoint[faceI] = 12345; // mark to be split
3338  }
3339  }
3340  }
3341  }
3342 
3343 
3344  // Synchronize faceMidPoint across coupled patches. (logical or)
3346  (
3347  mesh_,
3348  faceMidPoint,
3349  maxEqOp<label>(),
3350  false
3351  );
3352 
3353 
3354 
3355  // Introduce face points
3356  // ~~~~~~~~~~~~~~~~~~~~~
3357 
3358  {
3359  // Phase 1: determine mid points and sync. See comment for edgeMids
3360  // above
3361  pointField bFaceMids
3362  (
3363  mesh_.nFaces()-mesh_.nInternalFaces(),
3364  point(-GREAT, -GREAT, -GREAT)
3365  );
3366 
3367  forAll(bFaceMids, i)
3368  {
3369  label faceI = i+mesh_.nInternalFaces();
3370 
3371  if (faceMidPoint[faceI] >= 0)
3372  {
3373  bFaceMids[i] = mesh_.faceCentres()[faceI];
3374  }
3375  }
3377  (
3378  mesh_,
3379  bFaceMids,
3380  maxEqOp<vector>(),
3381  true // apply separation
3382  );
3383 
3384  forAll(faceMidPoint, faceI)
3385  {
3386  if (faceMidPoint[faceI] >= 0)
3387  {
3388  // Face marked to be split. Replace faceMidPoint with actual
3389  // point label.
3390 
3391  const face& f = mesh_.faces()[faceI];
3392 
3393  faceMidPoint[faceI] = meshMod.setAction
3394  (
3395  polyAddPoint
3396  (
3397  (
3398  faceI < mesh_.nInternalFaces()
3399  ? mesh_.faceCentres()[faceI]
3400  : bFaceMids[faceI-mesh_.nInternalFaces()]
3401  ), // point
3402  f[0], // master point
3403  -1, // zone for point
3404  true // supports a cell
3405  )
3406  );
3407 
3408  // Determine the level of the corner points and midpoint will
3409  // be one higher.
3410  newPointLevel(faceMidPoint[faceI]) = faceAnchorLevel[faceI]+1;
3411  }
3412  }
3413  }
3414 
3415  if (debug)
3416  {
3417  faceSet splitFaces(mesh_, "splitFaces", cellLabels.size());
3418 
3419  forAll(faceMidPoint, faceI)
3420  {
3421  if (faceMidPoint[faceI] >= 0)
3422  {
3423  splitFaces.insert(faceI);
3424  }
3425  }
3426 
3427  Pout<< "hexRef8::setRefinement : Dumping " << splitFaces.size()
3428  << " faces to split to faceSet " << splitFaces.objectPath() << endl;
3429 
3430  splitFaces.write();
3431  }
3432 
3433 
3434  // Information complete
3435  // ~~~~~~~~~~~~~~~~~~~~
3436  // At this point we have all the information we need. We should no
3437  // longer reference the cellLabels to refine. All the information is:
3438  // - cellMidPoint >= 0 : cell needs to be split
3439  // - faceMidPoint >= 0 : face needs to be split
3440  // - edgeMidPoint >= 0 : edge needs to be split
3441 
3442 
3443 
3444  // Get the corner/anchor points
3445  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3446 
3447  if (debug)
3448  {
3449  Pout<< "hexRef8::setRefinement :"
3450  << " Finding cell anchorPoints (8 per cell)"
3451  << endl;
3452  }
3453 
3454  // There will always be 8 points on the hex that have were introduced
3455  // with the hex and will have the same or lower refinement level.
3456 
3457  // Per cell the 8 corner points.
3458  labelListList cellAnchorPoints(mesh_.nCells());
3459 
3460  {
3461  labelList nAnchorPoints(mesh_.nCells(), 0);
3462 
3463  forAll(cellMidPoint, cellI)
3464  {
3465  if (cellMidPoint[cellI] >= 0)
3466  {
3467  cellAnchorPoints[cellI].setSize(8);
3468  }
3469  }
3470 
3471  forAll(pointLevel_, pointI)
3472  {
3473  const labelList& pCells = mesh_.pointCells(pointI);
3474 
3475  forAll(pCells, pCellI)
3476  {
3477  label cellI = pCells[pCellI];
3478 
3479  if
3480  (
3481  cellMidPoint[cellI] >= 0
3482  && pointLevel_[pointI] <= cellLevel_[cellI]
3483  )
3484  {
3485  if (nAnchorPoints[cellI] == 8)
3486  {
3487  dumpCell(cellI);
3488  FatalErrorIn
3489  (
3490  "hexRef8::setRefinement(const labelList&"
3491  ", polyTopoChange&)"
3492  ) << "cell " << cellI
3493  << " of level " << cellLevel_[cellI]
3494  << " uses more than 8 points of equal or"
3495  << " lower level" << nl
3496  << "Points so far:" << cellAnchorPoints[cellI]
3497  << abort(FatalError);
3498  }
3499 
3500  cellAnchorPoints[cellI][nAnchorPoints[cellI]++] = pointI;
3501  }
3502  }
3503  }
3504 
3505  forAll(cellMidPoint, cellI)
3506  {
3507  if (cellMidPoint[cellI] >= 0)
3508  {
3509  if (nAnchorPoints[cellI] != 8)
3510  {
3511  dumpCell(cellI);
3512 
3513  const labelList& cPoints = mesh_.cellPoints(cellI);
3514 
3515  FatalErrorIn
3516  (
3517  "hexRef8::setRefinement(const labelList&"
3518  ", polyTopoChange&)"
3519  ) << "cell " << cellI
3520  << " of level " << cellLevel_[cellI]
3521  << " does not seem to have 8 points of equal or"
3522  << " lower level" << endl
3523  << "cellPoints:" << cPoints << endl
3524  << "pointLevels:"
3525  << UIndirectList<label>(pointLevel_, cPoints)() << endl
3526  << abort(FatalError);
3527  }
3528  }
3529  }
3530  }
3531 
3532 
3533  // Add the cells
3534  // ~~~~~~~~~~~~~
3535 
3536  if (debug)
3537  {
3538  Pout<< "hexRef8::setRefinement :"
3539  << " Adding cells (1 per anchorPoint)"
3540  << endl;
3541  }
3542 
3543  // Per cell the 7 added cells (+ original cell)
3544  labelListList cellAddedCells(mesh_.nCells());
3545 
3546  forAll(cellAnchorPoints, cellI)
3547  {
3548  const labelList& cAnchors = cellAnchorPoints[cellI];
3549 
3550  if (cAnchors.size() == 8)
3551  {
3552  labelList& cAdded = cellAddedCells[cellI];
3553  cAdded.setSize(8);
3554 
3555  // Original cell at 0
3556  cAdded[0] = cellI;
3557 
3558  // Update cell level
3559  newCellLevel[cellI] = cellLevel_[cellI]+1;
3560 
3561 
3562  for (label i = 1; i < 8; i++)
3563  {
3564  cAdded[i] = meshMod.setAction
3565  (
3566  polyAddCell
3567  (
3568  -1, // master point
3569  -1, // master edge
3570  -1, // master face
3571  cellI, // master cell
3572  mesh_.cellZones().whichZone(cellI) // zone for cell
3573  )
3574  );
3575 
3576  newCellLevel(cAdded[i]) = cellLevel_[cellI]+1;
3577  }
3578  }
3579  }
3580 
3581 
3582  // Faces
3583  // ~~~~~
3584  // 1. existing faces that get split (into four always)
3585  // 2. existing faces that do not get split but only edges get split
3586  // 3. existing faces that do not get split but get new owner/neighbour
3587  // 4. new internal faces inside split cells.
3588 
3589  if (debug)
3590  {
3591  Pout<< "hexRef8::setRefinement :"
3592  << " Marking faces to be handled"
3593  << endl;
3594  }
3595 
3596  // Get all affected faces.
3597  PackedBoolList affectedFace(mesh_.nFaces());
3598 
3599  {
3600  forAll(cellMidPoint, cellI)
3601  {
3602  if (cellMidPoint[cellI] >= 0)
3603  {
3604  const cell& cFaces = mesh_.cells()[cellI];
3605 
3606  forAll(cFaces, i)
3607  {
3608  affectedFace.set(cFaces[i], 1);
3609  }
3610  }
3611  }
3612 
3613  forAll(faceMidPoint, faceI)
3614  {
3615  if (faceMidPoint[faceI] >= 0)
3616  {
3617  affectedFace.set(faceI, 1);
3618  }
3619  }
3620 
3621  forAll(edgeMidPoint, edgeI)
3622  {
3623  if (edgeMidPoint[edgeI] >= 0)
3624  {
3625  const labelList& eFaces = mesh_.edgeFaces(edgeI);
3626 
3627  forAll(eFaces, i)
3628  {
3629  affectedFace.set(eFaces[i], 1);
3630  }
3631  }
3632  }
3633  }
3634 
3635 
3636  // 1. Faces that get split
3637  // ~~~~~~~~~~~~~~~~~~~~~~~
3638 
3639  if (debug)
3640  {
3641  Pout<< "hexRef8::setRefinement : Splitting faces" << endl;
3642  }
3643 
3644  forAll(faceMidPoint, faceI)
3645  {
3646  if (faceMidPoint[faceI] >= 0 && affectedFace.get(faceI) == 1)
3647  {
3648  // Face needs to be split and hasn't yet been done in some way
3649  // (affectedFace - is impossible since this is first change but
3650  // just for completeness)
3651 
3652  const face& f = mesh_.faces()[faceI];
3653 
3654  // Has original faceI been used (three faces added, original gets
3655  // modified)
3656  bool modifiedFace = false;
3657 
3658  label anchorLevel = faceAnchorLevel[faceI];
3659 
3660  face newFace(4);
3661 
3662  forAll(f, fp)
3663  {
3664  label pointI = f[fp];
3665 
3666  if (pointLevel_[pointI] <= anchorLevel)
3667  {
3668  // point is anchor. Start collecting face.
3669 
3670  DynamicList<label> faceVerts(4);
3671 
3672  faceVerts.append(pointI);
3673 
3674  // Walk forward to mid point.
3675  // - if next is +2 midpoint is +1
3676  // - if next is +1 it is midpoint
3677  // - if next is +0 there has to be edgeMidPoint
3678 
3679  walkFaceToMid
3680  (
3681  edgeMidPoint,
3682  anchorLevel,
3683  faceI,
3684  fp,
3685  faceVerts
3686  );
3687 
3688  faceVerts.append(faceMidPoint[faceI]);
3689 
3690  walkFaceFromMid
3691  (
3692  edgeMidPoint,
3693  anchorLevel,
3694  faceI,
3695  fp,
3696  faceVerts
3697  );
3698 
3699  // Convert dynamiclist to face.
3700  newFace.transfer(faceVerts);
3701 
3702  //Pout<< "Split face:" << faceI << " verts:" << f
3703  // << " into quad:" << newFace << endl;
3704 
3705  // Get new owner/neighbour
3706  label own, nei;
3707  getFaceNeighbours
3708  (
3709  cellAnchorPoints,
3710  cellAddedCells,
3711  faceI,
3712  pointI, // Anchor point
3713 
3714  own,
3715  nei
3716  );
3717 
3718 
3719  if (debug)
3720  {
3721  if (mesh_.isInternalFace(faceI))
3722  {
3723  label oldOwn = mesh_.faceOwner()[faceI];
3724  label oldNei = mesh_.faceNeighbour()[faceI];
3725 
3726  checkInternalOrientation
3727  (
3728  meshMod,
3729  oldOwn,
3730  faceI,
3731  mesh_.cellCentres()[oldOwn],
3732  mesh_.cellCentres()[oldNei],
3733  newFace
3734  );
3735  }
3736  else
3737  {
3738  label oldOwn = mesh_.faceOwner()[faceI];
3739 
3740  checkBoundaryOrientation
3741  (
3742  meshMod,
3743  oldOwn,
3744  faceI,
3745  mesh_.cellCentres()[oldOwn],
3746  mesh_.faceCentres()[faceI],
3747  newFace
3748  );
3749  }
3750  }
3751 
3752 
3753  if (!modifiedFace)
3754  {
3755  modifiedFace = true;
3756 
3757  modFace(meshMod, faceI, newFace, own, nei);
3758  }
3759  else
3760  {
3761  addFace(meshMod, faceI, newFace, own, nei);
3762  }
3763  }
3764  }
3765 
3766  // Mark face as having been handled
3767  affectedFace.set(faceI, 0);
3768  }
3769  }
3770 
3771 
3772  // 2. faces that do not get split but use edges that get split
3773  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3774 
3775  if (debug)
3776  {
3777  Pout<< "hexRef8::setRefinement :"
3778  << " Adding edge splits to unsplit faces"
3779  << endl;
3780  }
3781 
3782  DynamicList<label> eFacesStorage;
3783  DynamicList<label> fEdgesStorage;
3784 
3785  forAll(edgeMidPoint, edgeI)
3786  {
3787  if (edgeMidPoint[edgeI] >= 0)
3788  {
3789  // Split edge. Check that face not already handled above.
3790 
3791  const labelList& eFaces = mesh_.edgeFaces(edgeI, eFacesStorage);
3792 
3793  forAll(eFaces, i)
3794  {
3795  label faceI = eFaces[i];
3796 
3797  if (faceMidPoint[faceI] < 0 && affectedFace.get(faceI) == 1)
3798  {
3799  // Unsplit face. Add edge splits to face.
3800 
3801  const face& f = mesh_.faces()[faceI];
3802  const labelList& fEdges = mesh_.faceEdges
3803  (
3804  faceI,
3805  fEdgesStorage
3806  );
3807 
3808  DynamicList<label> newFaceVerts(f.size());
3809 
3810  forAll(f, fp)
3811  {
3812  newFaceVerts.append(f[fp]);
3813 
3814  label edgeI = fEdges[fp];
3815 
3816  if (edgeMidPoint[edgeI] >= 0)
3817  {
3818  newFaceVerts.append(edgeMidPoint[edgeI]);
3819  }
3820  }
3821 
3822  face newFace;
3823  newFace.transfer(newFaceVerts);
3824 
3825  // The point with the lowest level should be an anchor
3826  // point of the neighbouring cells.
3827  label anchorFp = findMinLevel(f);
3828 
3829  label own, nei;
3830  getFaceNeighbours
3831  (
3832  cellAnchorPoints,
3833  cellAddedCells,
3834  faceI,
3835  f[anchorFp], // Anchor point
3836 
3837  own,
3838  nei
3839  );
3840 
3841 
3842  if (debug)
3843  {
3844  if (mesh_.isInternalFace(faceI))
3845  {
3846  label oldOwn = mesh_.faceOwner()[faceI];
3847  label oldNei = mesh_.faceNeighbour()[faceI];
3848 
3849  checkInternalOrientation
3850  (
3851  meshMod,
3852  oldOwn,
3853  faceI,
3854  mesh_.cellCentres()[oldOwn],
3855  mesh_.cellCentres()[oldNei],
3856  newFace
3857  );
3858  }
3859  else
3860  {
3861  label oldOwn = mesh_.faceOwner()[faceI];
3862 
3863  checkBoundaryOrientation
3864  (
3865  meshMod,
3866  oldOwn,
3867  faceI,
3868  mesh_.cellCentres()[oldOwn],
3869  mesh_.faceCentres()[faceI],
3870  newFace
3871  );
3872  }
3873  }
3874 
3875  modFace(meshMod, faceI, newFace, own, nei);
3876 
3877  // Mark face as having been handled
3878  affectedFace.set(faceI, 0);
3879  }
3880  }
3881  }
3882  }
3883 
3884 
3885  // 3. faces that do not get split but whose owner/neighbour change
3886  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3887 
3888  if (debug)
3889  {
3890  Pout<< "hexRef8::setRefinement :"
3891  << " Changing owner/neighbour for otherwise unaffected faces"
3892  << endl;
3893  }
3894 
3895  forAll(affectedFace, faceI)
3896  {
3897  if (affectedFace.get(faceI) == 1)
3898  {
3899  const face& f = mesh_.faces()[faceI];
3900 
3901  // The point with the lowest level should be an anchor
3902  // point of the neighbouring cells.
3903  label anchorFp = findMinLevel(f);
3904 
3905  label own, nei;
3906  getFaceNeighbours
3907  (
3908  cellAnchorPoints,
3909  cellAddedCells,
3910  faceI,
3911  f[anchorFp], // Anchor point
3912 
3913  own,
3914  nei
3915  );
3916 
3917  modFace(meshMod, faceI, f, own, nei);
3918 
3919  // Mark face as having been handled
3920  affectedFace.set(faceI, 0);
3921  }
3922  }
3923 
3924 
3925  // 4. new internal faces inside split cells.
3926  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3927 
3928 
3929  // This is the hard one. We have to find the splitting points between
3930  // the anchor points. But the edges between the anchor points might have
3931  // been split (into two,three or four edges).
3932 
3933  if (debug)
3934  {
3935  Pout<< "hexRef8::setRefinement :"
3936  << " Create new internal faces for split cells"
3937  << endl;
3938  }
3939 
3940  forAll(cellMidPoint, cellI)
3941  {
3942  if (cellMidPoint[cellI] >= 0)
3943  {
3944  createInternalFaces
3945  (
3946  cellAnchorPoints,
3947  cellAddedCells,
3948  cellMidPoint,
3949  faceMidPoint,
3950  faceAnchorLevel,
3951  edgeMidPoint,
3952  cellI,
3953  meshMod
3954  );
3955  }
3956  }
3957 
3958  // Extend pointLevels and cellLevels for the new cells. Could also be done
3959  // in updateMesh but saves passing cellAddedCells out of this routine.
3960 
3961  // Check
3962  if (debug)
3963  {
3964  label minPointI = labelMax;
3965  label maxPointI = labelMin;
3966 
3967  forAll(cellMidPoint, cellI)
3968  {
3969  if (cellMidPoint[cellI] >= 0)
3970  {
3971  minPointI = min(minPointI, cellMidPoint[cellI]);
3972  maxPointI = max(maxPointI, cellMidPoint[cellI]);
3973  }
3974  }
3975  forAll(faceMidPoint, faceI)
3976  {
3977  if (faceMidPoint[faceI] >= 0)
3978  {
3979  minPointI = min(minPointI, faceMidPoint[faceI]);
3980  maxPointI = max(maxPointI, faceMidPoint[faceI]);
3981  }
3982  }
3983  forAll(edgeMidPoint, edgeI)
3984  {
3985  if (edgeMidPoint[edgeI] >= 0)
3986  {
3987  minPointI = min(minPointI, edgeMidPoint[edgeI]);
3988  maxPointI = max(maxPointI, edgeMidPoint[edgeI]);
3989  }
3990  }
3991 
3992  if (minPointI != labelMax && minPointI != mesh_.nPoints())
3993  {
3994  FatalErrorIn("hexRef8::setRefinement(..)")
3995  << "Added point labels not consecutive to existing mesh points."
3996  << nl
3997  << "mesh_.nPoints():" << mesh_.nPoints()
3998  << " minPointI:" << minPointI
3999  << " maxPointI:" << maxPointI
4000  << abort(FatalError);
4001  }
4002  }
4003 
4004  pointLevel_.transfer(newPointLevel);
4005  cellLevel_.transfer(newCellLevel);
4006 
4007  // Mark files as changed
4008  setInstance(mesh_.facesInstance());
4009 
4010 
4011  // Update the live split cells tree.
4012  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4013 
4014  // New unrefinement structure
4015  if (history_.active())
4016  {
4017  if (debug)
4018  {
4019  Pout<< "hexRef8::setRefinement :"
4020  << " Updating refinement history to " << cellLevel_.size()
4021  << " cells" << endl;
4022  }
4023 
4024  // Extend refinement history for new cells
4025  history_.resize(cellLevel_.size());
4026 
4027  forAll(cellAddedCells, cellI)
4028  {
4029  const labelList& addedCells = cellAddedCells[cellI];
4030 
4031  if (addedCells.size())
4032  {
4033  // Cell was split.
4034  history_.storeSplit(cellI, addedCells);
4035  }
4036  }
4037  }
4038 
4039  // Compact cellAddedCells.
4040 
4041  labelListList refinedCells(cellLabels.size());
4042 
4043  forAll(cellLabels, i)
4044  {
4045  label cellI = cellLabels[i];
4046 
4047  refinedCells[i].transfer(cellAddedCells[cellI]);
4048  }
4049 
4050  return refinedCells;
4051 }
4052 
4053 
4056  const labelList& pointsToStore,
4057  const labelList& facesToStore,
4058  const labelList& cellsToStore
4059 )
4060 {
4061  savedPointLevel_.resize(2*pointsToStore.size());
4062  forAll(pointsToStore, i)
4063  {
4064  label pointI = pointsToStore[i];
4065  savedPointLevel_.insert(pointI, pointLevel_[pointI]);
4066  }
4067 
4068  savedCellLevel_.resize(2*cellsToStore.size());
4069  forAll(cellsToStore, i)
4070  {
4071  label cellI = cellsToStore[i];
4072  savedCellLevel_.insert(cellI, cellLevel_[cellI]);
4073  }
4074 }
4075 
4076 
4077 // Gets called after the mesh change. setRefinement will already have made
4078 // sure the pointLevel_ and cellLevel_ are the size of the new mesh so we
4079 // only need to account for reordering.
4081 {
4082  Map<label> dummyMap(0);
4083 
4084  updateMesh(map, dummyMap, dummyMap, dummyMap);
4085 }
4086 
4087 
4088 // Gets called after the mesh change. setRefinement will already have made
4089 // sure the pointLevel_ and cellLevel_ are the size of the new mesh so we
4090 // only need to account for reordering.
4093  const mapPolyMesh& map,
4094  const Map<label>& pointsToRestore,
4095  const Map<label>& facesToRestore,
4096  const Map<label>& cellsToRestore
4097 )
4098 {
4099  // Update celllevel
4100  if (debug)
4101  {
4102  Pout<< "hexRef8::updateMesh :"
4103  << " Updating various lists"
4104  << endl;
4105  }
4106 
4107  {
4108  const labelList& reverseCellMap = map.reverseCellMap();
4109 
4110  if (debug)
4111  {
4112  Pout<< "hexRef8::updateMesh :"
4113  << " reverseCellMap:" << map.reverseCellMap().size()
4114  << " cellMap:" << map.cellMap().size()
4115  << " nCells:" << mesh_.nCells()
4116  << " nOldCells:" << map.nOldCells()
4117  << " cellLevel_:" << cellLevel_.size()
4118  << " reversePointMap:" << map.reversePointMap().size()
4119  << " pointMap:" << map.pointMap().size()
4120  << " nPoints:" << mesh_.nPoints()
4121  << " nOldPoints:" << map.nOldPoints()
4122  << " pointLevel_:" << pointLevel_.size()
4123  << endl;
4124  }
4125 
4126  if (reverseCellMap.size() == cellLevel_.size())
4127  {
4128  // Assume it is after hexRef8 that this routine is called.
4129  // Just account for reordering. We cannot use cellMap since
4130  // then cells created from cells would get cellLevel_ of
4131  // cell they were created from.
4132  reorder(reverseCellMap, mesh_.nCells(), -1, cellLevel_);
4133  }
4134  else
4135  {
4136  // Map data
4137  const labelList& cellMap = map.cellMap();
4138 
4139  labelList newCellLevel(cellMap.size());
4140  forAll(cellMap, newCellI)
4141  {
4142  label oldCellI = cellMap[newCellI];
4143 
4144  if (oldCellI == -1)
4145  {
4146  newCellLevel[newCellI] = -1;
4147  }
4148  else
4149  {
4150  newCellLevel[newCellI] = cellLevel_[oldCellI];
4151  }
4152  }
4153  cellLevel_.transfer(newCellLevel);
4154  }
4155 
4156  // See if any cells to restore. This will be for some new cells
4157  // the corresponding old cell.
4158  forAllConstIter(Map<label>, cellsToRestore, iter)
4159  {
4160  label newCellI = iter.key();
4161  label storedCellI = iter();
4162 
4163  Map<label>::iterator fnd = savedCellLevel_.find(storedCellI);
4164 
4165  if (fnd == savedCellLevel_.end())
4166  {
4167  FatalErrorIn("hexRef8::updateMesh(const mapPolyMesh&)")
4168  << "Problem : trying to restore old value for new cell "
4169  << newCellI << " but cannot find old cell " << storedCellI
4170  << " in map of stored values " << savedCellLevel_
4171  << abort(FatalError);
4172  }
4173  cellLevel_[newCellI] = fnd();
4174  }
4175 
4176  //if (findIndex(cellLevel_, -1) != -1)
4177  //{
4178  // WarningIn("hexRef8::updateMesh(const mapPolyMesh&)")
4179  // << "Problem : "
4180  // << "cellLevel_ contains illegal value -1 after mapping
4181  // << " at cell " << findIndex(cellLevel_, -1) << endl
4182  // << "This means that another program has inflated cells"
4183  // << " (created cells out-of-nothing) and hence we don't know"
4184  // << " their cell level. Continuing with illegal value."
4185  // << abort(FatalError);
4186  //}
4187  }
4188 
4189 
4190  // Update pointlevel
4191  {
4192  const labelList& reversePointMap = map.reversePointMap();
4193 
4194  if (reversePointMap.size() == pointLevel_.size())
4195  {
4196  // Assume it is after hexRef8 that this routine is called.
4197  reorder(reversePointMap, mesh_.nPoints(), -1, pointLevel_);
4198  }
4199  else
4200  {
4201  // Map data
4202  const labelList& pointMap = map.pointMap();
4203 
4204  labelList newPointLevel(pointMap.size());
4205 
4206  forAll(pointMap, newPointI)
4207  {
4208  label oldPointI = pointMap[newPointI];
4209 
4210  if (oldPointI == -1)
4211  {
4212  //FatalErrorIn("hexRef8::updateMesh(const mapPolyMesh&)")
4213  // << "Problem : point " << newPointI
4214  // << " at " << mesh_.points()[newPointI]
4215  // << " does not originate from another point"
4216  // << " (i.e. is inflated)." << nl
4217  // << "Hence we cannot determine the new pointLevel"
4218  // << " for it." << abort(FatalError);
4219  newPointLevel[newPointI] = -1;
4220  }
4221  else
4222  {
4223  newPointLevel[newPointI] = pointLevel_[oldPointI];
4224  }
4225  }
4226  pointLevel_.transfer(newPointLevel);
4227  }
4228 
4229  // See if any points to restore. This will be for some new points
4230  // the corresponding old point (the one from the call to storeData)
4231  forAllConstIter(Map<label>, pointsToRestore, iter)
4232  {
4233  label newPointI = iter.key();
4234  label storedPointI = iter();
4235 
4236  Map<label>::iterator fnd = savedPointLevel_.find(storedPointI);
4237 
4238  if (fnd == savedPointLevel_.end())
4239  {
4240  FatalErrorIn("hexRef8::updateMesh(const mapPolyMesh&)")
4241  << "Problem : trying to restore old value for new point "
4242  << newPointI << " but cannot find old point "
4243  << storedPointI
4244  << " in map of stored values " << savedPointLevel_
4245  << abort(FatalError);
4246  }
4247  pointLevel_[newPointI] = fnd();
4248  }
4249 
4250  //if (findIndex(pointLevel_, -1) != -1)
4251  //{
4252  // WarningIn("hexRef8::updateMesh(const mapPolyMesh&)")
4253  // << "Problem : "
4254  // << "pointLevel_ contains illegal value -1 after mapping"
4255  // << " at point" << findIndex(pointLevel_, -1) << endl
4256  // << "This means that another program has inflated points"
4257  // << " (created points out-of-nothing) and hence we don't know"
4258  // << " their point level. Continuing with illegal value."
4259  // //<< abort(FatalError);
4260  //}
4261  }
4262 
4263  // Update refinement tree
4264  if (history_.active())
4265  {
4266  history_.updateMesh(map);
4267  }
4268 
4269  // Mark files as changed
4270  setInstance(mesh_.facesInstance());
4271 
4272  // Update face removal engine
4273  faceRemover_.updateMesh(map);
4274 }
4275 
4276 
4277 // Gets called after mesh subsetting. Maps are from new back to old.
4280  const labelList& pointMap,
4281  const labelList& faceMap,
4282  const labelList& cellMap
4283 )
4284 {
4285  // Update celllevel
4286  if (debug)
4287  {
4288  Pout<< "hexRef8::subset :"
4289  << " Updating various lists"
4290  << endl;
4291  }
4292 
4293  if (history_.active())
4294  {
4295  WarningIn
4296  (
4297  "hexRef8::subset(const labelList&, const labelList&"
4298  ", const labelList&)"
4299  ) << "Subsetting will not work in combination with unrefinement."
4300  << nl
4301  << "Proceed at your own risk." << endl;
4302  }
4303 
4304 
4305  // Update celllevel
4306  {
4307  labelList newCellLevel(cellMap.size());
4308 
4309  forAll(cellMap, newCellI)
4310  {
4311  newCellLevel[newCellI] = cellLevel_[cellMap[newCellI]];
4312  }
4313 
4314  cellLevel_.transfer(newCellLevel);
4315 
4316  if (findIndex(cellLevel_, -1) != -1)
4317  {
4318  FatalErrorIn("hexRef8::subset(..)")
4319  << "Problem : "
4320  << "cellLevel_ contains illegal value -1 after mapping:"
4321  << cellLevel_
4322  << abort(FatalError);
4323  }
4324  }
4325 
4326  // Update pointlevel
4327  {
4328  labelList newPointLevel(pointMap.size());
4329 
4330  forAll(pointMap, newPointI)
4331  {
4332  newPointLevel[newPointI] = pointLevel_[pointMap[newPointI]];
4333  }
4334 
4335  pointLevel_.transfer(newPointLevel);
4336 
4337  if (findIndex(pointLevel_, -1) != -1)
4338  {
4339  FatalErrorIn("hexRef8::subset(..)")
4340  << "Problem : "
4341  << "pointLevel_ contains illegal value -1 after mapping:"
4342  << pointLevel_
4343  << abort(FatalError);
4344  }
4345  }
4346 
4347  // Update refinement tree
4348  if (history_.active())
4349  {
4350  history_.subset(pointMap, faceMap, cellMap);
4351  }
4352 
4353  // Mark files as changed
4354  setInstance(mesh_.facesInstance());
4355 
4356  // Nothing needs doing to faceRemover.
4357  //faceRemover_.subset(pointMap, faceMap, cellMap);
4358 }
4359 
4360 
4361 // Gets called after the mesh distribution
4363 {
4364  if (debug)
4365  {
4366  Pout<< "hexRef8::distribute :"
4367  << " Updating various lists"
4368  << endl;
4369  }
4370 
4371  // Update celllevel
4372  map.distributeCellData(cellLevel_);
4373  // Update pointlevel
4374  map.distributePointData(pointLevel_);
4375 
4376  // Update refinement tree
4377  if (history_.active())
4378  {
4379  history_.distribute(map);
4380  }
4381 
4382  // Update face removal engine
4383  faceRemover_.distribute(map);
4384 }
4385 
4386 
4388 {
4389  const scalar smallDim = 1E-6 * mesh_.globalData().bb().mag();
4390 
4391  if (debug)
4392  {
4393  Pout<< "hexRef8::checkMesh : Using matching tolerance smallDim:"
4394  << smallDim << endl;
4395  }
4396 
4397  // Check owner on coupled faces.
4398  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4399 
4400  // There should be only one coupled face between two cells. Why? Since
4401  // otherwise mesh redistribution might cause multiple faces between two
4402  // cells
4403  {
4404  labelList nei(mesh_.nFaces()-mesh_.nInternalFaces());
4405  forAll(nei, i)
4406  {
4407  nei[i] = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4408  }
4409 
4410  // Replace data on coupled patches with their neighbour ones.
4411  syncTools::swapBoundaryFaceList(mesh_, nei, false);
4412 
4413  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
4414 
4415  forAll(patches, patchI)
4416  {
4417  const polyPatch& pp = patches[patchI];
4418 
4419  if (pp.coupled())
4420  {
4421  // Check how many faces between owner and neighbour. Should
4422  // be only one.
4424  cellToFace(2*pp.size());
4425 
4426  label faceI = pp.start();
4427 
4428  forAll(pp, i)
4429  {
4430  label own = mesh_.faceOwner()[faceI];
4431  label bFaceI = faceI-mesh_.nInternalFaces();
4432 
4433  if (!cellToFace.insert(labelPair(own, nei[bFaceI]), faceI))
4434  {
4435  dumpCell(own);
4436  FatalErrorIn("hexRef8::checkMesh()")
4437  << "Faces do not seem to be correct across coupled"
4438  << " boundaries" << endl
4439  << "Coupled face " << faceI
4440  << " between owner " << own
4441  << " on patch " << pp.name()
4442  << " and coupled neighbour " << nei[bFaceI]
4443  << " has two faces connected to it:"
4444  << faceI << " and "
4445  << cellToFace[labelPair(own, nei[bFaceI])]
4446  << abort(FatalError);
4447  }
4448 
4449  faceI++;
4450  }
4451  }
4452  }
4453  }
4454 
4455  // Check face areas.
4456  // ~~~~~~~~~~~~~~~~~
4457 
4458  {
4459  scalarField neiFaceAreas(mesh_.nFaces()-mesh_.nInternalFaces());
4460  forAll(neiFaceAreas, i)
4461  {
4462  neiFaceAreas[i] = mag(mesh_.faceAreas()[i+mesh_.nInternalFaces()]);
4463  }
4464 
4465  // Replace data on coupled patches with their neighbour ones.
4466  syncTools::swapBoundaryFaceList(mesh_, neiFaceAreas, false);
4467 
4468  forAll(neiFaceAreas, i)
4469  {
4470  label faceI = i+mesh_.nInternalFaces();
4471 
4472  const scalar magArea = mag(mesh_.faceAreas()[faceI]);
4473 
4474  if (mag(magArea - neiFaceAreas[i]) > smallDim)
4475  {
4476  const face& f = mesh_.faces()[faceI];
4477  label patchI = mesh_.boundaryMesh().whichPatch(faceI);
4478 
4479  dumpCell(mesh_.faceOwner()[faceI]);
4480 
4481  FatalErrorIn("hexRef8::checkMesh()")
4482  << "Faces do not seem to be correct across coupled"
4483  << " boundaries" << endl
4484  << "Coupled face " << faceI
4485  << " on patch " << patchI
4486  << " " << mesh_.boundaryMesh()[patchI].name()
4487  << " coords:" << UIndirectList<point>(mesh_.points(), f)()
4488  << " has face area:" << magArea
4489  << " (coupled) neighbour face area differs:"
4490  << neiFaceAreas[i]
4491  << " to within tolerance " << smallDim
4492  << abort(FatalError);
4493  }
4494  }
4495  }
4496 
4497 
4498  // Check number of points on faces.
4499  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4500  {
4501  labelList nVerts(mesh_.nFaces()-mesh_.nInternalFaces());
4502 
4503  forAll(nVerts, i)
4504  {
4505  nVerts[i] = mesh_.faces()[i+mesh_.nInternalFaces()].size();
4506  }
4507 
4508  // Replace data on coupled patches with their neighbour ones.
4509  syncTools::swapBoundaryFaceList(mesh_, nVerts, false);
4510 
4511  forAll(nVerts, i)
4512  {
4513  label faceI = i+mesh_.nInternalFaces();
4514 
4515  const face& f = mesh_.faces()[faceI];
4516 
4517  if (f.size() != nVerts[i])
4518  {
4519  dumpCell(mesh_.faceOwner()[faceI]);
4520 
4521  label patchI = mesh_.boundaryMesh().whichPatch(faceI);
4522 
4523  FatalErrorIn("hexRef8::checkMesh()")
4524  << "Faces do not seem to be correct across coupled"
4525  << " boundaries" << endl
4526  << "Coupled face " << faceI
4527  << " on patch " << patchI
4528  << " " << mesh_.boundaryMesh()[patchI].name()
4529  << " coords:" << UIndirectList<point>(mesh_.points(), f)()
4530  << " has size:" << f.size()
4531  << " (coupled) neighbour face has size:"
4532  << nVerts[i]
4533  << abort(FatalError);
4534  }
4535  }
4536  }
4537 
4538 
4539  // Check points of face
4540  // ~~~~~~~~~~~~~~~~~~~~
4541  {
4542  // Anchor points.
4543  pointField anchorPoints(mesh_.nFaces()-mesh_.nInternalFaces());
4544 
4545  forAll(anchorPoints, i)
4546  {
4547  label faceI = i+mesh_.nInternalFaces();
4548  const point& fc = mesh_.faceCentres()[faceI];
4549  const face& f = mesh_.faces()[faceI];
4550  const vector anchorVec(mesh_.points()[f[0]] - fc);
4551 
4552  anchorPoints[i] = anchorVec;
4553  }
4554 
4555  // Replace data on coupled patches with their neighbour ones. Apply
4556  // rotation transformation (but not separation since is relative vector
4557  // to point on same face.
4558  syncTools::swapBoundaryFaceList(mesh_, anchorPoints, false);
4559 
4560  forAll(anchorPoints, i)
4561  {
4562  label faceI = i+mesh_.nInternalFaces();
4563  const point& fc = mesh_.faceCentres()[faceI];
4564  const face& f = mesh_.faces()[faceI];
4565  const vector anchorVec(mesh_.points()[f[0]] - fc);
4566 
4567  if (mag(anchorVec - anchorPoints[i]) > smallDim)
4568  {
4569  dumpCell(mesh_.faceOwner()[faceI]);
4570 
4571  label patchI = mesh_.boundaryMesh().whichPatch(faceI);
4572 
4573  FatalErrorIn("hexRef8::checkMesh()")
4574  << "Faces do not seem to be correct across coupled"
4575  << " boundaries" << endl
4576  << "Coupled face " << faceI
4577  << " on patch " << patchI
4578  << " " << mesh_.boundaryMesh()[patchI].name()
4579  << " coords:" << UIndirectList<point>(mesh_.points(), f)()
4580  << " has anchor vector:" << anchorVec
4581  << " (coupled) neighbour face anchor vector differs:"
4582  << anchorPoints[i]
4583  << " to within tolerance " << smallDim
4584  << abort(FatalError);
4585  }
4586  }
4587  }
4588 
4589  if (debug)
4590  {
4591  Pout<< "hexRef8::checkMesh : Returning" << endl;
4592  }
4593 }
4594 
4595 
4598  const label maxPointDiff,
4599  const labelList& pointsToCheck
4600 ) const
4601 {
4602  if (debug)
4603  {
4604  Pout<< "hexRef8::checkRefinementLevels :"
4605  << " Checking 2:1 refinement level" << endl;
4606  }
4607 
4608  if
4609  (
4610  cellLevel_.size() != mesh_.nCells()
4611  || pointLevel_.size() != mesh_.nPoints()
4612  )
4613  {
4614  FatalErrorIn("hexRef8::checkRefinementLevels(const label)")
4615  << "cellLevel size should be number of cells"
4616  << " and pointLevel size should be number of points."<< nl
4617  << "cellLevel:" << cellLevel_.size()
4618  << " mesh.nCells():" << mesh_.nCells() << nl
4619  << "pointLevel:" << pointLevel_.size()
4620  << " mesh.nPoints():" << mesh_.nPoints()
4621  << abort(FatalError);
4622  }
4623 
4624 
4625  // Check 2:1 consistency.
4626  // ~~~~~~~~~~~~~~~~~~~~~~
4627 
4628  {
4629  // Internal faces.
4630  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
4631  {
4632  label own = mesh_.faceOwner()[faceI];
4633  label nei = mesh_.faceNeighbour()[faceI];
4634 
4635  if (mag(cellLevel_[own] - cellLevel_[nei]) > 1)
4636  {
4637  dumpCell(own);
4638  dumpCell(nei);
4639 
4640  FatalErrorIn
4641  (
4642  "hexRef8::checkRefinementLevels(const label)"
4643  ) << "Celllevel does not satisfy 2:1 constraint." << nl
4644  << "On face " << faceI << " owner cell " << own
4645  << " has refinement " << cellLevel_[own]
4646  << " neighbour cell " << nei << " has refinement "
4647  << cellLevel_[nei]
4648  << abort(FatalError);
4649  }
4650  }
4651 
4652  // Coupled faces. Get neighbouring value
4653  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
4654 
4655  forAll(neiLevel, i)
4656  {
4657  label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4658 
4659  neiLevel[i] = cellLevel_[own];
4660  }
4661 
4662  // No separation
4663  syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
4664 
4665  forAll(neiLevel, i)
4666  {
4667  label faceI = i+mesh_.nInternalFaces();
4668 
4669  label own = mesh_.faceOwner()[faceI];
4670 
4671  if (mag(cellLevel_[own] - neiLevel[i]) > 1)
4672  {
4673  dumpCell(own);
4674 
4675  label patchI = mesh_.boundaryMesh().whichPatch(faceI);
4676 
4677  FatalErrorIn
4678  (
4679  "hexRef8::checkRefinementLevels(const label)"
4680  ) << "Celllevel does not satisfy 2:1 constraint."
4681  << " On coupled face " << faceI
4682  << " on patch " << patchI << " "
4683  << mesh_.boundaryMesh()[patchI].name()
4684  << " owner cell " << own << " has refinement "
4685  << cellLevel_[own]
4686  << " (coupled) neighbour cell has refinement "
4687  << neiLevel[i]
4688  << abort(FatalError);
4689  }
4690  }
4691  }
4692 
4693 
4694  // Check pointLevel is synchronized
4695  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4696  {
4697  labelList syncPointLevel(pointLevel_);
4698 
4699  // Get min level
4701  (
4702  mesh_,
4703  syncPointLevel,
4704  minEqOp<label>(),
4705  labelMax,
4706  false // no separation
4707  );
4708 
4709 
4710  forAll(syncPointLevel, pointI)
4711  {
4712  if (pointLevel_[pointI] != syncPointLevel[pointI])
4713  {
4714  FatalErrorIn
4715  (
4716  "hexRef8::checkRefinementLevels(const label)"
4717  ) << "PointLevel is not consistent across coupled patches."
4718  << endl
4719  << "point:" << pointI << " coord:" << mesh_.points()[pointI]
4720  << " has level " << pointLevel_[pointI]
4721  << " whereas the coupled point has level "
4722  << syncPointLevel[pointI]
4723  << abort(FatalError);
4724  }
4725  }
4726  }
4727 
4728 
4729  // Check 2:1 across points (instead of faces)
4730  if (maxPointDiff != -1)
4731  {
4732  // Determine per point the max cell level.
4733  labelList maxPointLevel(mesh_.nPoints(), 0);
4734 
4735  forAll(maxPointLevel, pointI)
4736  {
4737  const labelList& pCells = mesh_.pointCells(pointI);
4738 
4739  label& pLevel = maxPointLevel[pointI];
4740 
4741  forAll(pCells, i)
4742  {
4743  pLevel = max(pLevel, cellLevel_[pCells[i]]);
4744  }
4745  }
4746 
4747  // Sync maxPointLevel to neighbour
4749  (
4750  mesh_,
4751  maxPointLevel,
4752  maxEqOp<label>(),
4753  labelMin, // null value
4754  false // no separation
4755  );
4756 
4757  // Check 2:1 across boundary points
4758  forAll(pointsToCheck, i)
4759  {
4760  label pointI = pointsToCheck[i];
4761 
4762  const labelList& pCells = mesh_.pointCells(pointI);
4763 
4764  forAll(pCells, i)
4765  {
4766  label cellI = pCells[i];
4767 
4768  if
4769  (
4770  mag(cellLevel_[cellI]-maxPointLevel[pointI])
4771  > maxPointDiff
4772  )
4773  {
4774  dumpCell(cellI);
4775 
4776  FatalErrorIn
4777  (
4778  "hexRef8::checkRefinementLevels(const label)"
4779  ) << "Too big a difference between"
4780  << " point-connected cells." << nl
4781  << "cell:" << cellI
4782  << " cellLevel:" << cellLevel_[cellI]
4783  << " uses point:" << pointI
4784  << " coord:" << mesh_.points()[pointI]
4785  << " which is also used by a cell with level:"
4786  << maxPointLevel[pointI]
4787  << abort(FatalError);
4788  }
4789  }
4790  }
4791  }
4792 
4793 
4794  //- Gives problems after first splitting off inside mesher.
4796  //{
4797  // // Any patches with points having only two edges.
4798  //
4799  // boolList isHangingPoint(mesh_.nPoints(), false);
4800  //
4801  // const polyBoundaryMesh& patches = mesh_.boundaryMesh();
4802  //
4803  // forAll(patches, patchI)
4804  // {
4805  // const polyPatch& pp = patches[patchI];
4806  //
4807  // const labelList& meshPoints = pp.meshPoints();
4808  //
4809  // forAll(meshPoints, i)
4810  // {
4811  // label pointI = meshPoints[i];
4812  //
4813  // const labelList& pEdges = mesh_.pointEdges()[pointI];
4814  //
4815  // if (pEdges.size() == 2)
4816  // {
4817  // isHangingPoint[pointI] = true;
4818  // }
4819  // }
4820  // }
4821  //
4822  // syncTools::syncPointList
4823  // (
4824  // mesh_,
4825  // isHangingPoint,
4826  // andEqOp<bool>(), // only if all decide it is hanging point
4827  // true, // null
4828  // false // no separation
4829  // );
4830  //
4831  // //OFstream str(mesh_.time().path()/"hangingPoints.obj");
4832  //
4833  // label nHanging = 0;
4834  //
4835  // forAll(isHangingPoint, pointI)
4836  // {
4837  // if (isHangingPoint[pointI])
4838  // {
4839  // nHanging++;
4840  //
4841  // Pout<< "Hanging boundary point " << pointI
4842  // << " at " << mesh_.points()[pointI]
4843  // << endl;
4844  // //meshTools::writeOBJ(str, mesh_.points()[pointI]);
4845  // }
4846  // }
4847  //
4848  // if (returnReduce(nHanging, sumOp<label>()) > 0)
4849  // {
4850  // FatalErrorIn
4851  // (
4852  // "hexRef8::checkRefinementLevels(const label)"
4853  // ) << "Detected a point used by two edges only (hanging point)"
4854  // << nl << "This is not allowed"
4855  // << abort(FatalError);
4856  // }
4857  //}
4858 }
4859 
4860 
4861 
4862 //
4863 // Unrefinement
4864 // ~~~~~~~~~~~~
4865 //
4866 
4867 
4869 {
4870  if (debug)
4871  {
4872  checkRefinementLevels(-1, labelList(0));
4873  }
4874 
4875  if (debug)
4876  {
4877  Pout<< "hexRef8::getSplitPoints :"
4878  << " Calculating unrefineable points" << endl;
4879  }
4880 
4881 
4882  if (!history_.active())
4883  {
4884  FatalErrorIn("hexRef8::getSplitPoints()")
4885  << "Only call if constructed with history capability"
4886  << abort(FatalError);
4887  }
4888 
4889  // Master cell
4890  // -1 undetermined
4891  // -2 certainly not split point
4892  // >= label of master cell
4893  labelList splitMaster(mesh_.nPoints(), -1);
4894  labelList splitMasterLevel(mesh_.nPoints(), 0);
4895 
4896  // Unmark all with not 8 cells
4897  //const labelListList& pointCells = mesh_.pointCells();
4898 
4899  for (label pointI = 0; pointI < mesh_.nPoints(); pointI++)
4900  {
4901  const labelList& pCells = mesh_.pointCells(pointI);
4902 
4903  if (pCells.size() != 8)
4904  {
4905  splitMaster[pointI] = -2;
4906  }
4907  }
4908 
4909  // Unmark all with different master cells
4910  const labelList& visibleCells = history_.visibleCells();
4911 
4912  forAll(visibleCells, cellI)
4913  {
4914  const labelList& cPoints = mesh_.cellPoints(cellI);
4915 
4916  if (visibleCells[cellI] != -1 && history_.parentIndex(cellI) >= 0)
4917  {
4918  label parentIndex = history_.parentIndex(cellI);
4919 
4920  // Check same master.
4921  forAll(cPoints, i)
4922  {
4923  label pointI = cPoints[i];
4924 
4925  label masterCellI = splitMaster[pointI];
4926 
4927  if (masterCellI == -1)
4928  {
4929  // First time visit of point. Store parent cell and
4930  // level of the parent cell (with respect to cellI). This
4931  // is additional guarantee that we're referring to the
4932  // same master at the same refinement level.
4933 
4934  splitMaster[pointI] = parentIndex;
4935  splitMasterLevel[pointI] = cellLevel_[cellI] - 1;
4936  }
4937  else if (masterCellI == -2)
4938  {
4939  // Already decided that point is not splitPoint
4940  }
4941  else if
4942  (
4943  (masterCellI != parentIndex)
4944  || (splitMasterLevel[pointI] != cellLevel_[cellI] - 1)
4945  )
4946  {
4947  // Different masters so point is on two refinement
4948  // patterns
4949  splitMaster[pointI] = -2;
4950  }
4951  }
4952  }
4953  else
4954  {
4955  // Either not visible or is unrefined cell
4956  forAll(cPoints, i)
4957  {
4958  label pointI = cPoints[i];
4959 
4960  splitMaster[pointI] = -2;
4961  }
4962  }
4963  }
4964 
4965  // Unmark boundary faces
4966  for
4967  (
4968  label faceI = mesh_.nInternalFaces();
4969  faceI < mesh_.nFaces();
4970  faceI++
4971  )
4972  {
4973  const face& f = mesh_.faces()[faceI];
4974 
4975  forAll(f, fp)
4976  {
4977  splitMaster[f[fp]] = -2;
4978  }
4979  }
4980 
4981 
4982  // Collect into labelList
4983 
4984  label nSplitPoints = 0;
4985 
4986  forAll(splitMaster, pointI)
4987  {
4988  if (splitMaster[pointI] >= 0)
4989  {
4990  nSplitPoints++;
4991  }
4992  }
4993 
4994  labelList splitPoints(nSplitPoints);
4995  nSplitPoints = 0;
4996 
4997  forAll(splitMaster, pointI)
4998  {
4999  if (splitMaster[pointI] >= 0)
5000  {
5001  splitPoints[nSplitPoints++] = pointI;
5002  }
5003  }
5004 
5005  return splitPoints;
5006 }
5007 
5008 
5009 //void Foam::hexRef8::markIndex
5010 //(
5011 // const label maxLevel,
5012 // const label level,
5013 // const label index,
5014 // const label markValue,
5015 // labelList& indexValues
5016 //) const
5017 //{
5018 // if (level < maxLevel && indexValues[index] == -1)
5019 // {
5020 // // Mark
5021 // indexValues[index] = markValue;
5022 //
5023 // // Mark parent
5024 // const splitCell8& split = history_.splitCells()[index];
5025 //
5026 // if (split.parent_ >= 0)
5027 // {
5028 // markIndex
5029 // (
5030 // maxLevel, level+1, split.parent_, markValue, indexValues);
5031 // )
5032 // }
5033 // }
5034 //}
5035 //
5036 //
5041 //void Foam::hexRef8::markCellClusters
5042 //(
5043 // const label maxLevel,
5044 // labelList& cluster
5045 //) const
5046 //{
5047 // cluster.setSize(mesh_.nCells());
5048 // cluster = -1;
5049 //
5050 // const DynamicList<splitCell8>& splitCells = history_.splitCells();
5051 //
5052 // // Mark all splitCells down to level maxLevel with a cell originating from
5053 // // it.
5054 //
5055 // labelList indexLevel(splitCells.size(), -1);
5056 //
5057 // forAll(visibleCells, cellI)
5058 // {
5059 // label index = visibleCells[cellI];
5060 //
5061 // if (index >= 0)
5062 // {
5063 // markIndex(maxLevel, 0, index, cellI, indexLevel);
5064 // }
5065 // }
5066 //
5067 // // Mark cells with splitCell
5068 //}
5069 
5070 
5073  const labelList& pointsToUnrefine,
5074  const bool maxSet
5075 ) const
5076 {
5077  if (debug)
5078  {
5079  Pout<< "hexRef8::consistentUnrefinement :"
5080  << " Determining 2:1 consistent unrefinement" << endl;
5081  }
5082 
5083  if (maxSet)
5084  {
5085  FatalErrorIn
5086  (
5087  "hexRef8::consistentUnrefinement(const labelList&, const bool"
5088  ) << "maxSet not implemented yet."
5089  << abort(FatalError);
5090  }
5091 
5092  // Loop, modifying pointsToUnrefine, until no more changes to due to 2:1
5093  // conflicts.
5094  // maxSet = false : unselect points to refine
5095  // maxSet = true: select points to refine
5096 
5097  // Maintain boolList for pointsToUnrefine and cellsToUnrefine
5098  PackedBoolList unrefinePoint(mesh_.nPoints());
5099 
5100  forAll(pointsToUnrefine, i)
5101  {
5102  label pointI = pointsToUnrefine[i];
5103 
5104  unrefinePoint.set(pointI, 1);
5105  }
5106 
5107 
5108  while (true)
5109  {
5110  // Construct cells to unrefine
5111  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
5112 
5113  PackedBoolList unrefineCell(mesh_.nCells());
5114 
5115  forAll(unrefinePoint, pointI)
5116  {
5117  if (unrefinePoint.get(pointI) == 1)
5118  {
5119  const labelList& pCells = mesh_.pointCells(pointI);
5120 
5121  forAll(pCells, j)
5122  {
5123  unrefineCell.set(pCells[j], 1);
5124  }
5125  }
5126  }
5127 
5128 
5129  label nChanged = 0;
5130 
5131 
5132  // Check 2:1 consistency taking refinement into account
5133  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5134 
5135  // Internal faces.
5136  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
5137  {
5138  label own = mesh_.faceOwner()[faceI];
5139  label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5140 
5141  label nei = mesh_.faceNeighbour()[faceI];
5142  label neiLevel = cellLevel_[nei] - unrefineCell.get(nei);
5143 
5144  if (ownLevel < (neiLevel-1))
5145  {
5146  // Since was 2:1 this can only occur if own is marked for
5147  // unrefinement.
5148 
5149  if (maxSet)
5150  {
5151  unrefineCell.set(nei, 1);
5152  }
5153  else
5154  {
5155  if (unrefineCell.get(own) == 0)
5156  {
5157  FatalErrorIn("hexRef8::consistentUnrefinement(..)")
5158  << "problem" << abort(FatalError);
5159  }
5160 
5161  unrefineCell.set(own, 0);
5162  }
5163  nChanged++;
5164  }
5165  else if (neiLevel < (ownLevel-1))
5166  {
5167  if (maxSet)
5168  {
5169  unrefineCell.set(own, 1);
5170  }
5171  else
5172  {
5173  if (unrefineCell.get(nei) == 0)
5174  {
5175  FatalErrorIn("hexRef8::consistentUnrefinement(..)")
5176  << "problem" << abort(FatalError);
5177  }
5178 
5179  unrefineCell.set(nei, 0);
5180  }
5181  nChanged++;
5182  }
5183  }
5184 
5185 
5186  // Coupled faces. Swap owner level to get neighbouring cell level.
5187  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
5188 
5189  forAll(neiLevel, i)
5190  {
5191  label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
5192 
5193  neiLevel[i] = cellLevel_[own] - unrefineCell.get(own);
5194  }
5195 
5196  // Swap to neighbour
5197  syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
5198 
5199  forAll(neiLevel, i)
5200  {
5201  label faceI = i+mesh_.nInternalFaces();
5202  label own = mesh_.faceOwner()[faceI];
5203  label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5204 
5205  if (ownLevel < (neiLevel[i]-1))
5206  {
5207  if (!maxSet)
5208  {
5209  if (unrefineCell.get(own) == 0)
5210  {
5211  FatalErrorIn("hexRef8::consistentUnrefinement(..)")
5212  << "problem" << abort(FatalError);
5213  }
5214 
5215  unrefineCell.set(own, 0);
5216  nChanged++;
5217  }
5218  }
5219  else if (neiLevel[i] < (ownLevel-1))
5220  {
5221  if (maxSet)
5222  {
5223  if (unrefineCell.get(own) == 1)
5224  {
5225  FatalErrorIn("hexRef8::consistentUnrefinement(..)")
5226  << "problem" << abort(FatalError);
5227  }
5228 
5229  unrefineCell.set(own, 1);
5230  nChanged++;
5231  }
5232  }
5233  }
5234 
5235  reduce(nChanged, sumOp<label>());
5236 
5237  if (debug)
5238  {
5239  Pout<< "hexRef8::consistentUnrefinement :"
5240  << " Changed " << nChanged
5241  << " refinement levels due to 2:1 conflicts."
5242  << endl;
5243  }
5244 
5245  if (nChanged == 0)
5246  {
5247  break;
5248  }
5249 
5250 
5251  // Convert cellsToUnrefine back into points to unrefine
5252  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5253 
5254  // Knock out any point whose cell neighbour cannot be unrefined.
5255  forAll(unrefinePoint, pointI)
5256  {
5257  if (unrefinePoint.get(pointI) == 1)
5258  {
5259  const labelList& pCells = mesh_.pointCells(pointI);
5260 
5261  forAll(pCells, j)
5262  {
5263  if (unrefineCell.get(pCells[j]) == 0)
5264  {
5265  unrefinePoint.set(pointI, 0);
5266  break;
5267  }
5268  }
5269  }
5270  }
5271  }
5272 
5273 
5274  // Convert back to labelList.
5275  label nSet = 0;
5276 
5277  forAll(unrefinePoint, pointI)
5278  {
5279  if (unrefinePoint.get(pointI) == 1)
5280  {
5281  nSet++;
5282  }
5283  }
5284 
5285  labelList newPointsToUnrefine(nSet);
5286  nSet = 0;
5287 
5288  forAll(unrefinePoint, pointI)
5289  {
5290  if (unrefinePoint.get(pointI) == 1)
5291  {
5292  newPointsToUnrefine[nSet++] = pointI;
5293  }
5294  }
5295 
5296  return newPointsToUnrefine;
5297 }
5298 
5299 
5302  const labelList& splitPointLabels,
5303  polyTopoChange& meshMod
5304 )
5305 {
5306  if (!history_.active())
5307  {
5308  FatalErrorIn
5309  (
5310  "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
5311  ) << "Only call if constructed with history capability"
5312  << abort(FatalError);
5313  }
5314 
5315  if (debug)
5316  {
5317  Pout<< "hexRef8::setUnrefinement :"
5318  << " Checking initial mesh just to make sure" << endl;
5319 
5320  checkMesh();
5321 
5322  forAll(cellLevel_, cellI)
5323  {
5324  if (cellLevel_[cellI] < 0)
5325  {
5326  FatalErrorIn
5327  (
5328  "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
5329  ) << "Illegal cell level " << cellLevel_[cellI]
5330  << " for cell " << cellI
5331  << abort(FatalError);
5332  }
5333  }
5334 
5335 
5336  // Write to sets.
5337  pointSet pSet(mesh_, "splitPoints", splitPointLabels);
5338  pSet.write();
5339 
5340  cellSet cSet(mesh_, "splitPointCells", splitPointLabels.size());
5341 
5342  forAll(splitPointLabels, i)
5343  {
5344  const labelList& pCells = mesh_.pointCells(splitPointLabels[i]);
5345 
5346  forAll(pCells, j)
5347  {
5348  cSet.insert(pCells[j]);
5349  }
5350  }
5351  cSet.write();
5352 
5353  Pout<< "hexRef8::setRefinement : Dumping " << pSet.size()
5354  << " points and "
5355  << cSet.size() << " cells for unrefinement to" << nl
5356  << " pointSet " << pSet.objectPath() << nl
5357  << " cellSet " << cSet.objectPath()
5358  << endl;
5359  }
5360 
5361 
5362  labelList cellRegion;
5363  labelList cellRegionMaster;
5364  labelList facesToRemove;
5365 
5366  {
5367  labelHashSet splitFaces(12*splitPointLabels.size());
5368 
5369  forAll(splitPointLabels, i)
5370  {
5371  const labelList& pFaces = mesh_.pointFaces()[splitPointLabels[i]];
5372 
5373  forAll(pFaces, j)
5374  {
5375  splitFaces.insert(pFaces[j]);
5376  }
5377  }
5378 
5379  // Check with faceRemover what faces will get removed. Note that this
5380  // can be more (but never less) than splitFaces provided.
5381  faceRemover_.compatibleRemoves
5382  (
5383  splitFaces.toc(), // pierced faces
5384  cellRegion, // per cell -1 or region it is merged into
5385  cellRegionMaster, // per region the master cell
5386  facesToRemove // new faces to be removed.
5387  );
5388 
5389  if (facesToRemove.size() != splitFaces.size())
5390  {
5391  FatalErrorIn
5392  (
5393  "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
5394  ) << "Ininitial set of split points to unrefine does not"
5395  << " seem to be consistent or not mid points of refined cells"
5396  << abort(FatalError);
5397  }
5398  }
5399 
5400  // Redo the region master so it is consistent with our master.
5401  // This will guarantee that the new cell (for which faceRemover uses
5402  // the region master) is already compatible with our refinement structure.
5403 
5404  forAll(splitPointLabels, i)
5405  {
5406  label pointI = splitPointLabels[i];
5407 
5408  // Get original cell label
5409 
5410  const labelList& pCells = mesh_.pointCells(pointI);
5411 
5412  // Check
5413  if (pCells.size() != 8)
5414  {
5415  FatalErrorIn
5416  (
5417  "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
5418  ) << "splitPoint " << pointI
5419  << " should have 8 cells using it. It has " << pCells
5420  << abort(FatalError);
5421  }
5422 
5423 
5424  // Check that the lowest numbered pCells is the master of the region
5425  // (should be guaranteed by directRemoveFaces)
5426  //if (debug)
5427  {
5428  label masterCellI = min(pCells);
5429 
5430  forAll(pCells, j)
5431  {
5432  label cellI = pCells[j];
5433 
5434  label region = cellRegion[cellI];
5435 
5436  if (region == -1)
5437  {
5438  FatalErrorIn("hexRef8::setUnrefinement(..)")
5439  << "Ininitial set of split points to unrefine does not"
5440  << " seem to be consistent or not mid points"
5441  << " of refined cells" << nl
5442  << "cell:" << cellI << " on splitPoint " << pointI
5443  << " has no region to be merged into"
5444  << abort(FatalError);
5445  }
5446 
5447  if (masterCellI != cellRegionMaster[region])
5448  {
5449  FatalErrorIn("hexRef8::setUnrefinement(..)")
5450  << "cell:" << cellI << " on splitPoint:" << pointI
5451  << " in region " << region
5452  << " has master:" << cellRegionMaster[region]
5453  << " which is not the lowest numbered cell"
5454  << " among the pointCells:" << pCells
5455  << abort(FatalError);
5456  }
5457  }
5458  }
5459  }
5460 
5461  // Insert all commands to combine cells. Never fails so don't have to
5462  // test for success.
5463  faceRemover_.setRefinement
5464  (
5465  facesToRemove,
5466  cellRegion,
5467  cellRegionMaster,
5468  meshMod
5469  );
5470 
5471  // Remove the 8 cells that originated from merging around the split point
5472  // and adapt cell levels (not that pointLevels stay the same since points
5473  // either get removed or stay at the same position.
5474  forAll(splitPointLabels, i)
5475  {
5476  label pointI = splitPointLabels[i];
5477 
5478  const labelList& pCells = mesh_.pointCells(pointI);
5479 
5480  label masterCellI = min(pCells);
5481 
5482  forAll(pCells, j)
5483  {
5484  cellLevel_[pCells[j]]--;
5485  }
5486 
5487  history_.combineCells(masterCellI, pCells);
5488  }
5489 
5490  // Mark files as changed
5491  setInstance(mesh_.facesInstance());
5492 
5493  // history_.updateMesh will take care of truncating.
5494 }
5495 
5496 
5497 // Write refinement to polyMesh directory.
5499 {
5500  bool writeOk = cellLevel_.write() && pointLevel_.write();
5501 
5502  if (history_.active())
5503  {
5504  writeOk = writeOk && history_.write();
5505  }
5506 
5507  return writeOk;
5508 }
5509 
5510 
5511 // ************************ vim: set sw=4 sts=4 et: ************************ //