FreeFOAM The Cross-Platform CFD Toolkit
dynamicRefineFvMesh.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "dynamicRefineFvMesh.H"
28 #include <finiteVolume/volFields.H>
31 #include <finiteVolume/fvCFD.H>
32 #include <OpenFOAM/syncTools.H>
33 #include <OpenFOAM/pointFields.H>
34 
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39 
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41 
43 
45 
46 
47 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48 
50 (
51  const PackedBoolList& l,
52  const unsigned int val
53 )
54 {
55  label n = 0;
56  forAll(l, i)
57  {
58  if (l.get(i) == val)
59  {
60  n++;
61  }
62  }
63  return n;
64 }
65 
66 
68 (
69  PackedBoolList& unrefineableCell
70 ) const
71 {
72  if (protectedCell_.empty())
73  {
74  unrefineableCell.clear();
75  return;
76  }
77 
78  const labelList& cellLevel = meshCutter_.cellLevel();
79 
80  unrefineableCell = protectedCell_;
81 
82  // Get neighbouring cell level
83  labelList neiLevel(nFaces()-nInternalFaces());
84 
85  for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
86  {
87  neiLevel[faceI-nInternalFaces()] = cellLevel[faceOwner()[faceI]];
88  }
89  syncTools::swapBoundaryFaceList(*this, neiLevel, false);
90 
91 
92  while (true)
93  {
94  // Pick up faces on border of protected cells
95  boolList seedFace(nFaces(), false);
96 
97  forAll(faceNeighbour(), faceI)
98  {
99  label own = faceOwner()[faceI];
100  bool ownProtected = (unrefineableCell.get(own) == 1);
101  label nei = faceNeighbour()[faceI];
102  bool neiProtected = (unrefineableCell.get(nei) == 1);
103 
104  if (ownProtected && (cellLevel[nei] > cellLevel[own]))
105  {
106  seedFace[faceI] = true;
107  }
108  else if (neiProtected && (cellLevel[own] > cellLevel[nei]))
109  {
110  seedFace[faceI] = true;
111  }
112  }
113  for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
114  {
115  label own = faceOwner()[faceI];
116  bool ownProtected = (unrefineableCell.get(own) == 1);
117  if
118  (
119  ownProtected
120  && (neiLevel[faceI-nInternalFaces()] > cellLevel[own])
121  )
122  {
123  seedFace[faceI] = true;
124  }
125  }
126 
127  syncTools::syncFaceList(*this, seedFace, orEqOp<bool>(), false);
128 
129 
130  // Extend unrefineableCell
131  bool hasExtended = false;
132 
133  for (label faceI = 0; faceI < nInternalFaces(); faceI++)
134  {
135  if (seedFace[faceI])
136  {
137  label own = faceOwner()[faceI];
138  if (unrefineableCell.get(own) == 0)
139  {
140  unrefineableCell.set(own, 1);
141  hasExtended = true;
142  }
143 
144  label nei = faceNeighbour()[faceI];
145  if (unrefineableCell.get(nei) == 0)
146  {
147  unrefineableCell.set(nei, 1);
148  hasExtended = true;
149  }
150  }
151  }
152  for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
153  {
154  if (seedFace[faceI])
155  {
156  label own = faceOwner()[faceI];
157  if (unrefineableCell.get(own) == 0)
158  {
159  unrefineableCell.set(own, 1);
160  hasExtended = true;
161  }
162  }
163  }
164 
165  if (!returnReduce(hasExtended, orOp<bool>()))
166  {
167  break;
168  }
169  }
170 }
171 
172 
174 {
175  dictionary refineDict
176  (
178  (
179  IOobject
180  (
181  "dynamicMeshDict",
182  time().constant(),
183  *this,
186  false
187  )
188  ).subDict(typeName + "Coeffs")
189  );
190 
191  correctFluxes_ = List<Pair<word> >(refineDict.lookup("correctFluxes"));
192 
193  dumpLevel_ = Switch(refineDict.lookup("dumpLevel"));
194 }
195 
196 
197 // Refines cells, maps fields and recalculates (an approximate) flux
199 (
200  const labelList& cellsToRefine
201 )
202 {
203  // Mesh changing engine.
204  polyTopoChange meshMod(*this);
205 
206  // Play refinement commands into mesh changer.
207  meshCutter_.setRefinement(cellsToRefine, meshMod);
208 
209  // Create mesh (with inflation), return map from old to new mesh.
210  //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
211  autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
212 
213  Info<< "Refined from "
214  << returnReduce(map().nOldCells(), sumOp<label>())
215  << " to " << globalData().nTotalCells() << " cells." << endl;
216 
217  if (debug)
218  {
219  // Check map.
220  for (label faceI = 0; faceI < nInternalFaces(); faceI++)
221  {
222  label oldFaceI = map().faceMap()[faceI];
223 
224  if (oldFaceI >= nInternalFaces())
225  {
226  FatalErrorIn("dynamicRefineFvMesh::refine(const labelList&)")
227  << "New internal face:" << faceI
228  << " fc:" << faceCentres()[faceI]
229  << " originates from boundary oldFace:" << oldFaceI
230  << abort(FatalError);
231  }
232  }
233  }
234 
235 
236  // Update fields
237  updateMesh(map);
238 
239  // Move mesh
240  /*
241  pointField newPoints;
242  if (map().hasMotionPoints())
243  {
244  newPoints = map().preMotionPoints();
245  }
246  else
247  {
248  newPoints = points();
249  }
250  movePoints(newPoints);
251  */
252 
253  // Correct the flux for modified/added faces. All the faces which only
254  // have been renumbered will already have been handled by the mapping.
255  {
256  const labelList& faceMap = map().faceMap();
257  const labelList& reverseFaceMap = map().reverseFaceMap();
258 
259  // Storage for any master faces. These will be the original faces
260  // on the coarse cell that get split into four (or rather the
261  // master face gets modified and three faces get added from the master)
262  labelHashSet masterFaces(4*cellsToRefine.size());
263 
264  forAll(faceMap, faceI)
265  {
266  label oldFaceI = faceMap[faceI];
267 
268  if (oldFaceI >= 0)
269  {
270  label masterFaceI = reverseFaceMap[oldFaceI];
271 
272  if (masterFaceI < 0)
273  {
275  (
276  "dynamicRefineFvMesh::refine(const labelList&)"
277  ) << "Problem: should not have removed faces"
278  << " when refining."
279  << nl << "face:" << faceI << abort(FatalError);
280  }
281  else if (masterFaceI != faceI)
282  {
283  masterFaces.insert(masterFaceI);
284  }
285  }
286  }
287  if (debug)
288  {
289  Info<< "Found " << returnReduce(masterFaces.size(), sumOp<label>())
290  << " split faces " << endl;
291  }
292 
293  forAll(correctFluxes_, i)
294  {
295  if (debug)
296  {
297  Info<< "Mapping flux " << correctFluxes_[i][0]
298  << " using interpolated flux " << correctFluxes_[i][1]
299  << endl;
300  }
302  (
303  lookupObject<surfaceScalarField>(correctFluxes_[i][0])
304  );
307  (
308  lookupObject<volVectorField>(correctFluxes_[i][1])
309  )
310  & Sf();
311 
312  // Recalculate new internal faces.
313  for (label faceI = 0; faceI < nInternalFaces(); faceI++)
314  {
315  label oldFaceI = faceMap[faceI];
316 
317  if (oldFaceI == -1)
318  {
319  // Inflated/appended
320  phi[faceI] = phiU[faceI];
321  }
322  else if (reverseFaceMap[oldFaceI] != faceI)
323  {
324  // face-from-masterface
325  phi[faceI] = phiU[faceI];
326  }
327  }
328 
329  // Recalculate new boundary faces.
330  forAll(phi.boundaryField(), patchI)
331  {
332  fvsPatchScalarField& patchPhi = phi.boundaryField()[patchI];
333  const fvsPatchScalarField& patchPhiU =
334  phiU.boundaryField()[patchI];
335 
336  label faceI = patchPhi.patch().patch().start();
337 
338  forAll(patchPhi, i)
339  {
340  label oldFaceI = faceMap[faceI];
341 
342  if (oldFaceI == -1)
343  {
344  // Inflated/appended
345  patchPhi[i] = patchPhiU[i];
346  }
347  else if (reverseFaceMap[oldFaceI] != faceI)
348  {
349  // face-from-masterface
350  patchPhi[i] = patchPhiU[i];
351  }
352 
353  faceI++;
354  }
355  }
356 
357  // Update master faces
358  forAllConstIter(labelHashSet, masterFaces, iter)
359  {
360  label faceI = iter.key();
361 
362  if (isInternalFace(faceI))
363  {
364  phi[faceI] = phiU[faceI];
365  }
366  else
367  {
368  label patchI = boundaryMesh().whichPatch(faceI);
369  label i = faceI - boundaryMesh()[patchI].start();
370 
371  const fvsPatchScalarField& patchPhiU =
372  phiU.boundaryField()[patchI];
373 
374  fvsPatchScalarField& patchPhi =
375  phi.boundaryField()[patchI];
376 
377  patchPhi[i] = patchPhiU[i];
378  }
379  }
380  }
381  }
382 
383 
384 
385  // Update numbering of cells/vertices.
386  meshCutter_.updateMesh(map);
387 
388  // Update numbering of protectedCell_
389  if (protectedCell_.size())
390  {
391  PackedBoolList newProtectedCell(nCells());
392 
393  forAll(newProtectedCell, cellI)
394  {
395  label oldCellI = map().cellMap()[cellI];
396  newProtectedCell.set(cellI, protectedCell_.get(oldCellI));
397  }
398  protectedCell_.transfer(newProtectedCell);
399  }
400 
401  // Debug: Check refinement levels (across faces only)
402  meshCutter_.checkRefinementLevels(-1, labelList(0));
403 
404  return map;
405 }
406 
407 
408 // Combines previously split cells, maps fields and recalculates
409 // (an approximate) flux
411 (
412  const labelList& splitPoints
413 )
414 {
415  polyTopoChange meshMod(*this);
416 
417  // Play refinement commands into mesh changer.
418  meshCutter_.setUnrefinement(splitPoints, meshMod);
419 
420 
421  // Save information on faces that will be combined
422  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
423 
424  // Find the faceMidPoints on cells to be combined.
425  // for each face resulting of split of face into four store the
426  // midpoint
427  Map<label> faceToSplitPoint(3*splitPoints.size());
428 
429  {
430  forAll(splitPoints, i)
431  {
432  label pointI = splitPoints[i];
433 
434  const labelList& pEdges = pointEdges()[pointI];
435 
436  forAll(pEdges, j)
437  {
438  label otherPointI = edges()[pEdges[j]].otherVertex(pointI);
439 
440  const labelList& pFaces = pointFaces()[otherPointI];
441 
442  forAll(pFaces, pFaceI)
443  {
444  faceToSplitPoint.insert(pFaces[pFaceI], otherPointI);
445  }
446  }
447  }
448  }
449 
450 
451  // Change mesh and generate map.
452  //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
453  autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
454 
455  Info<< "Unrefined from "
456  << returnReduce(map().nOldCells(), sumOp<label>())
457  << " to " << globalData().nTotalCells() << " cells."
458  << endl;
459 
460  // Update fields
461  updateMesh(map);
462 
463 
464  // Move mesh
465  /*
466  pointField newPoints;
467  if (map().hasMotionPoints())
468  {
469  newPoints = map().preMotionPoints();
470  }
471  else
472  {
473  newPoints = points();
474  }
475  movePoints(newPoints);
476  */
477 
478  // Correct the flux for modified faces.
479  {
480  const labelList& reversePointMap = map().reversePointMap();
481  const labelList& reverseFaceMap = map().reverseFaceMap();
482 
483  forAll(correctFluxes_, i)
484  {
485  if (debug)
486  {
487  Info<< "Mapping flux " << correctFluxes_[i][0]
488  << " using interpolated flux " << correctFluxes_[i][1]
489  << endl;
490  }
492  (
493  lookupObject<surfaceScalarField>(correctFluxes_[i][0])
494  );
497  (
498  lookupObject<volVectorField>(correctFluxes_[i][1])
499  )
500  & Sf();
501 
502  forAllConstIter(Map<label>, faceToSplitPoint, iter)
503  {
504  label oldFaceI = iter.key();
505  label oldPointI = iter();
506 
507  if (reversePointMap[oldPointI] < 0)
508  {
509  // midpoint was removed. See if face still exists.
510  label faceI = reverseFaceMap[oldFaceI];
511 
512  if (faceI >= 0)
513  {
514  if (isInternalFace(faceI))
515  {
516  phi[faceI] = phiU[faceI];
517  }
518  else
519  {
520  label patchI = boundaryMesh().whichPatch(faceI);
521  label i = faceI - boundaryMesh()[patchI].start();
522 
523  const fvsPatchScalarField& patchPhiU =
524  phiU.boundaryField()[patchI];
525 
526  fvsPatchScalarField& patchPhi =
527  phi.boundaryField()[patchI];
528 
529  patchPhi[i] = patchPhiU[i];
530  }
531  }
532  }
533  }
534  }
535  }
536 
537 
538  // Update numbering of cells/vertices.
539  meshCutter_.updateMesh(map);
540 
541  // Update numbering of protectedCell_
542  if (protectedCell_.size())
543  {
544  PackedBoolList newProtectedCell(nCells());
545 
546  forAll(newProtectedCell, cellI)
547  {
548  label oldCellI = map().cellMap()[cellI];
549  if (oldCellI >= 0)
550  {
551  newProtectedCell.set(cellI, protectedCell_.get(oldCellI));
552  }
553  }
554  protectedCell_.transfer(newProtectedCell);
555  }
556 
557  // Debug: Check refinement levels (across faces only)
558  meshCutter_.checkRefinementLevels(-1, labelList(0));
559 
560  return map;
561 }
562 
563 
564 // Get max of connected point
566 {
567  scalarField vFld(nCells(), -GREAT);
568 
569  forAll(pointCells(), pointI)
570  {
571  const labelList& pCells = pointCells()[pointI];
572 
573  forAll(pCells, i)
574  {
575  vFld[pCells[i]] = max(vFld[pCells[i]], pFld[pointI]);
576  }
577  }
578  return vFld;
579 }
580 
581 
582 // Get min of connected cell
584 {
585  scalarField pFld(nPoints(), GREAT);
586 
587  forAll(pointCells(), pointI)
588  {
589  const labelList& pCells = pointCells()[pointI];
590 
591  forAll(pCells, i)
592  {
593  pFld[pointI] = min(pFld[pointI], vFld[pCells[i]]);
594  }
595  }
596  return pFld;
597 }
598 
599 
600 // Simple (non-parallel) interpolation by averaging.
602 {
603  scalarField pFld(nPoints());
604 
605  forAll(pointCells(), pointI)
606  {
607  const labelList& pCells = pointCells()[pointI];
608 
609  scalar sum = 0.0;
610  forAll(pCells, i)
611  {
612  sum += vFld[pCells[i]];
613  }
614  pFld[pointI] = sum/pCells.size();
615  }
616  return pFld;
617 }
618 
619 
620 // Calculate error. Is < 0 or distance to minLevel, maxLevel
622 (
623  const scalarField& fld,
624  const scalar minLevel,
625  const scalar maxLevel
626 ) const
627 {
628  scalarField c(fld.size(), -1);
629 
630  forAll(fld, i)
631  {
632  scalar err = min(fld[i]-minLevel, maxLevel-fld[i]);
633 
634  if (err >= 0)
635  {
636  c[i] = err;
637  }
638  }
639  return c;
640 }
641 
642 
644 (
645  const scalar lowerRefineLevel,
646  const scalar upperRefineLevel,
647  const scalarField& vFld,
648  PackedBoolList& candidateCell
649 ) const
650 {
651  // Get error per cell. Is -1 (not to be refined) to >0 (to be refined,
652  // higher more desirable to be refined).
653  scalarField cellError
654  (
655  maxPointField
656  (
657  error
658  (
659  cellToPoint(vFld),
660  lowerRefineLevel,
661  upperRefineLevel
662  )
663  )
664  );
665 
666  // Mark cells that are candidates for refinement.
667  forAll(cellError, cellI)
668  {
669  if (cellError[cellI] > 0)
670  {
671  candidateCell.set(cellI, 1);
672  }
673  }
674 }
675 
676 
678 (
679  const label maxCells,
680  const label maxRefinement,
681  const PackedBoolList& candidateCell
682 ) const
683 {
684  // Every refined cell causes 7 extra cells
685  label nTotToRefine = (maxCells - globalData().nTotalCells()) / 7;
686 
687  const labelList& cellLevel = meshCutter_.cellLevel();
688 
689  // Mark cells that cannot be refined since they would trigger refinement
690  // of protected cells (since 2:1 cascade)
691  PackedBoolList unrefineableCell;
692  calculateProtectedCells(unrefineableCell);
693 
694  // Count current selection
695  label nCandidates = returnReduce(count(candidateCell, 1), sumOp<label>());
696 
697  // Collect all cells
698  DynamicList<label> candidates(nCells());
699 
700  if (nCandidates < nTotToRefine)
701  {
702  forAll(candidateCell, cellI)
703  {
704  if
705  (
706  cellLevel[cellI] < maxRefinement
707  && candidateCell.get(cellI) == 1
708  && (
709  unrefineableCell.empty()
710  || unrefineableCell.get(cellI) == 0
711  )
712  )
713  {
714  candidates.append(cellI);
715  }
716  }
717  }
718  else
719  {
720  // Sort by error? For now just truncate.
721  for (label level = 0; level < maxRefinement; level++)
722  {
723  forAll(candidateCell, cellI)
724  {
725  if
726  (
727  cellLevel[cellI] == level
728  && candidateCell.get(cellI) == 1
729  && (
730  unrefineableCell.empty()
731  || unrefineableCell.get(cellI) == 0
732  )
733  )
734  {
735  candidates.append(cellI);
736  }
737  }
738 
739  if (returnReduce(candidates.size(), sumOp<label>()) > nTotToRefine)
740  {
741  break;
742  }
743  }
744  }
745 
746  // Guarantee 2:1 refinement after refinement
747  labelList consistentSet
748  (
749  meshCutter_.consistentRefinement
750  (
751  candidates.shrink(),
752  true // Add to set to guarantee 2:1
753  )
754  );
755 
756  Info<< "Selected " << returnReduce(consistentSet.size(), sumOp<label>())
757  << " cells for refinement out of " << globalData().nTotalCells()
758  << "." << endl;
759 
760  return consistentSet;
761 }
762 
763 
765 (
766  const scalar unrefineLevel,
767  const PackedBoolList& markedCell,
768  const scalarField& pFld
769 ) const
770 {
771  // All points that can be unrefined
772  const labelList splitPoints(meshCutter_.getSplitPoints());
773 
774  DynamicList<label> newSplitPoints(splitPoints.size());
775 
776  forAll(splitPoints, i)
777  {
778  label pointI = splitPoints[i];
779 
780  if (pFld[pointI] < unrefineLevel)
781  {
782  // Check that all cells are not marked
783  const labelList& pCells = pointCells()[pointI];
784 
785  bool hasMarked = false;
786 
787  forAll(pCells, pCellI)
788  {
789  if (markedCell.get(pCells[pCellI]) == 1)
790  {
791  hasMarked = true;
792  break;
793  }
794  }
795 
796  if (!hasMarked)
797  {
798  newSplitPoints.append(pointI);
799  }
800  }
801  }
802 
803 
804  newSplitPoints.shrink();
805 
806  // Guarantee 2:1 refinement after unrefinement
807  labelList consistentSet
808  (
809  meshCutter_.consistentUnrefinement
810  (
811  newSplitPoints,
812  false
813  )
814  );
815  Info<< "Selected " << returnReduce(consistentSet.size(), sumOp<label>())
816  << " split points out of a possible "
817  << returnReduce(splitPoints.size(), sumOp<label>())
818  << "." << endl;
819 
820  return consistentSet;
821 }
822 
823 
825 {
826  // Mark faces using any marked cell
827  boolList markedFace(nFaces(), false);
828 
829  forAll(markedCell, cellI)
830  {
831  if (markedCell.get(cellI) == 1)
832  {
833  const cell& cFaces = cells()[cellI];
834 
835  forAll(cFaces, i)
836  {
837  markedFace[cFaces[i]] = true;
838  }
839  }
840  }
841 
842  syncTools::syncFaceList(*this, markedFace, orEqOp<bool>(), false);
843 
844  // Update cells using any markedFace
845  for (label faceI = 0; faceI < nInternalFaces(); faceI++)
846  {
847  if (markedFace[faceI])
848  {
849  markedCell.set(faceOwner()[faceI], 1);
850  markedCell.set(faceNeighbour()[faceI], 1);
851  }
852  }
853  for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
854  {
855  if (markedFace[faceI])
856  {
857  markedCell.set(faceOwner()[faceI], 1);
858  }
859  }
860 }
861 
862 
863 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
864 
865 dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
866 :
867  dynamicFvMesh(io),
868  meshCutter_(*this),
869  dumpLevel_(false),
870  nRefinementIterations_(0),
871  protectedCell_(nCells(), 0)
872 {
873  // Read static part of dictionary
874  readDict();
875 
876 
877  const labelList& cellLevel = meshCutter_.cellLevel();
878  const labelList& pointLevel = meshCutter_.pointLevel();
879 
880  // Set cells that should not be refined.
881  // This is currently any cell which does not have 8 anchor points or
882  // uses any face which does not have 4 anchor points.
883  // Note: do not use cellPoint addressing
884 
885  // Count number of points <= cellLevel
886  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
887 
888  label nProtected = 0;
889 
890  forAll(cellLevel, cellI)
891  {
892  const labelList& cPoints = cellPoints(cellI);
893 
894  label nAnchors = 0;
895  forAll(cPoints, cPointI)
896  {
897  label pointI = cPoints[cPointI];
898  if (pointLevel[pointI] <= cellLevel[cellI])
899  {
900  nAnchors++;
901  }
902  }
903  if (nAnchors != 8)
904  {
905  protectedCell_.set(cellI, 1);
906  nProtected++;
907  }
908  }
909 
910 
911  // Count number of points <= faceLevel
912  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
913  // Bit tricky since proc face might be one more refined than the owner since
914  // the coupled one is refined.
915 
916  {
917  labelList neiLevel(nFaces());
918 
919  for (label faceI = 0; faceI < nInternalFaces(); faceI++)
920  {
921  neiLevel[faceI] = cellLevel[faceNeighbour()[faceI]];
922  }
923  for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
924  {
925  neiLevel[faceI] = cellLevel[faceOwner()[faceI]];
926  }
927  syncTools::swapFaceList(*this, neiLevel, false);
928 
929 
930  boolList protectedFace(nFaces(), false);
931 
932  forAll(faceOwner(), faceI)
933  {
934  label faceLevel = max
935  (
936  cellLevel[faceOwner()[faceI]],
937  neiLevel[faceI]
938  );
939 
940  const face& f = faces()[faceI];
941 
942  label nAnchors = 0;
943 
944  forAll(f, fp)
945  {
946  if (pointLevel[f[fp]] <= faceLevel)
947  {
948  nAnchors++;
949  }
950  }
951 
952  if (nAnchors != 4)
953  {
954  protectedFace[faceI] = true;
955  }
956  }
957 
959  (
960  *this,
961  protectedFace,
962  orEqOp<bool>(),
963  false
964  );
965 
966  for (label faceI = 0; faceI < nInternalFaces(); faceI++)
967  {
968  if (protectedFace[faceI])
969  {
970  if (protectedCell_.set(faceOwner()[faceI], 1))
971  {
972  nProtected++;
973  }
974  if (protectedCell_.set(faceNeighbour()[faceI], 1))
975  {
976  nProtected++;
977  }
978  }
979  }
980  for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
981  {
982  if (protectedFace[faceI])
983  {
984  if (protectedCell_.set(faceOwner()[faceI], 1))
985  {
986  nProtected++;
987  }
988  }
989  }
990  }
991 
992  reduce(nProtected, sumOp<label>());
993 
994  //Info<< "Protecting " << nProtected << " out of "
995  // << returnReduce(nCells(), sumOp<label>())
996  // << endl;
997 
998 
999 
1000  if (nProtected == 0)
1001  {
1003  }
1004 }
1005 
1006 
1007 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
1008 
1010 {}
1011 
1012 
1013 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1014 
1016 {
1017  // Re-read dictionary. Choosen since usually -small so trivial amount
1018  // of time compared to actual refinement. Also very useful to be able
1019  // to modify on-the-fly.
1020  dictionary refineDict
1021  (
1022  IOdictionary
1023  (
1024  IOobject
1025  (
1026  "dynamicMeshDict",
1027  time().constant(),
1028  *this,
1031  false
1032  )
1033  ).subDict(typeName + "Coeffs")
1034  );
1035 
1036  label refineInterval = readLabel(refineDict.lookup("refineInterval"));
1037 
1038  bool hasChanged = false;
1039 
1040  if (refineInterval == 0)
1041  {
1042  changing(hasChanged);
1043 
1044  return false;
1045  }
1046  else if (refineInterval < 0)
1047  {
1048  FatalErrorIn("dynamicRefineFvMesh::update()")
1049  << "Illegal refineInterval " << refineInterval << nl
1050  << "The refineInterval setting in the dynamicMeshDict should"
1051  << " be >= 1." << nl
1052  << exit(FatalError);
1053  }
1054 
1055 
1056 
1057 
1058  // Note: cannot refine at time 0 since no V0 present since mesh not
1059  // moved yet.
1060 
1061  if (time().timeIndex() > 0 && time().timeIndex() % refineInterval == 0)
1062  {
1063  label maxCells = readLabel(refineDict.lookup("maxCells"));
1064 
1065  if (maxCells <= 0)
1066  {
1067  FatalErrorIn("dynamicRefineFvMesh::update()")
1068  << "Illegal maximum number of cells " << maxCells << nl
1069  << "The maxCells setting in the dynamicMeshDict should"
1070  << " be > 0." << nl
1071  << exit(FatalError);
1072  }
1073 
1074  label maxRefinement = readLabel(refineDict.lookup("maxRefinement"));
1075 
1076  if (maxRefinement <= 0)
1077  {
1078  FatalErrorIn("dynamicRefineFvMesh::update()")
1079  << "Illegal maximum refinement level " << maxRefinement << nl
1080  << "The maxCells setting in the dynamicMeshDict should"
1081  << " be > 0." << nl
1082  << exit(FatalError);
1083  }
1084 
1085  word field(refineDict.lookup("field"));
1086 
1087  const volScalarField& vFld = lookupObject<volScalarField>(field);
1088 
1089  const scalar lowerRefineLevel =
1090  readScalar(refineDict.lookup("lowerRefineLevel"));
1091  const scalar upperRefineLevel =
1092  readScalar(refineDict.lookup("upperRefineLevel"));
1093  const scalar unrefineLevel =
1094  readScalar(refineDict.lookup("unrefineLevel"));
1095  const label nBufferLayers =
1096  readLabel(refineDict.lookup("nBufferLayers"));
1097 
1098  // Cells marked for refinement or otherwise protected from unrefinement.
1100 
1101  if (globalData().nTotalCells() < maxCells)
1102  {
1103  // Determine candidates for refinement (looking at field only)
1105  (
1106  lowerRefineLevel,
1107  upperRefineLevel,
1108  vFld,
1109  refineCell
1110  );
1111 
1112  // Select subset of candidates. Take into account max allowable
1113  // cells, refinement level, protected cells.
1114  labelList cellsToRefine
1115  (
1117  (
1118  maxCells,
1119  maxRefinement,
1120  refineCell
1121  )
1122  );
1123 
1124  label nCellsToRefine = returnReduce
1125  (
1126  cellsToRefine.size(), sumOp<label>()
1127  );
1128 
1129  if (nCellsToRefine > 0)
1130  {
1131  // Refine/update mesh and map fields
1132  autoPtr<mapPolyMesh> map = refine(cellsToRefine);
1133 
1134  // Update refineCell. Note that some of the marked ones have
1135  // not been refined due to constraints.
1136  {
1137  const labelList& cellMap = map().cellMap();
1138  const labelList& reverseCellMap = map().reverseCellMap();
1139 
1140  PackedBoolList newRefineCell(cellMap.size());
1141 
1142  forAll(cellMap, cellI)
1143  {
1144  label oldCellI = cellMap[cellI];
1145 
1146  if (oldCellI < 0)
1147  {
1148  newRefineCell.set(cellI, 1);
1149  }
1150  else if (reverseCellMap[oldCellI] != cellI)
1151  {
1152  newRefineCell.set(cellI, 1);
1153  }
1154  else
1155  {
1156  newRefineCell.set(cellI, refineCell.get(oldCellI));
1157  }
1158  }
1159  refineCell.transfer(newRefineCell);
1160  }
1161 
1162  // Extend with a buffer layer to prevent neighbouring points
1163  // being unrefined.
1164  for (label i = 0; i < nBufferLayers; i++)
1165  {
1166  extendMarkedCells(refineCell);
1167  }
1168 
1169  hasChanged = true;
1170  }
1171  }
1172 
1173 
1174  {
1175  // Select unrefineable points that are not marked in refineCell
1176  labelList pointsToUnrefine
1177  (
1179  (
1180  unrefineLevel,
1181  refineCell,
1182  minCellField(vFld)
1183  )
1184  );
1185 
1186  label nSplitPoints = returnReduce
1187  (
1188  pointsToUnrefine.size(),
1189  sumOp<label>()
1190  );
1191 
1192  if (nSplitPoints > 0)
1193  {
1194  // Refine/update mesh
1195  unrefine(pointsToUnrefine);
1196 
1197  hasChanged = true;
1198  }
1199  }
1200 
1201 
1202  if ((nRefinementIterations_ % 10) == 0)
1203  {
1204  // Compact refinement history occassionally (how often?).
1205  // Unrefinement causes holes in the refinementHistory.
1206  const_cast<refinementHistory&>(meshCutter().history()).compact();
1207  }
1209  }
1210 
1211  changing(hasChanged);
1212 
1213  return hasChanged;
1214 }
1215 
1216 
1222 ) const
1223 {
1224  // Force refinement data to go to the current time directory.
1225  const_cast<hexRef8&>(meshCutter_).setInstance(time().timeName());
1226 
1227  bool writeOk =
1228  dynamicFvMesh::writeObjects(fmt, ver, cmp)
1229  && meshCutter_.write();
1230 
1231  if (dumpLevel_)
1232  {
1233  volScalarField scalarCellLevel
1234  (
1235  IOobject
1236  (
1237  "cellLevel",
1238  time().timeName(),
1239  *this,
1242  false
1243  ),
1244  *this,
1245  dimensionedScalar("level", dimless, 0)
1246  );
1247 
1248  const labelList& cellLevel = meshCutter_.cellLevel();
1249 
1250  forAll(cellLevel, cellI)
1251  {
1252  scalarCellLevel[cellI] = cellLevel[cellI];
1253  }
1254 
1255  writeOk = writeOk && scalarCellLevel.write();
1256  }
1257 
1258  return writeOk;
1259 }
1260 
1261 
1262 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1263 
1264 } // End namespace Foam
1265 
1266 // ************************ vim: set sw=4 sts=4 et: ************************ //