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