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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2017 OpenFOAM Foundation
9  Copyright (C) 2018-2020 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "dynamicRefineFvMesh.H"
31 #include "surfaceInterpolate.H"
32 #include "volFields.H"
33 #include "polyTopoChange.H"
34 #include "surfaceFields.H"
35 #include "syncTools.H"
36 #include "pointFields.H"
37 #include "sigFpe.H"
38 #include "cellSet.H"
39 #include "HashOps.H"
40 
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42 
43 namespace Foam
44 {
45  defineTypeNameAndDebug(dynamicRefineFvMesh, 0);
46  addToRunTimeSelectionTable(dynamicFvMesh, dynamicRefineFvMesh, IOobject);
47  addToRunTimeSelectionTable(dynamicFvMesh, dynamicRefineFvMesh, doInit);
48 }
49 
50 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
51 
53 (
54  bitSet& unrefineableCell
55 ) const
56 {
57  if (protectedCell_.empty())
58  {
59  unrefineableCell.clear();
60  return;
61  }
62 
63  const labelList& cellLevel = meshCutter_.cellLevel();
64 
65  unrefineableCell = protectedCell_;
66 
67  // Get neighbouring cell level
68  labelList neiLevel(nFaces()-nInternalFaces());
69 
70  for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
71  {
72  neiLevel[facei-nInternalFaces()] = cellLevel[faceOwner()[facei]];
73  }
74  syncTools::swapBoundaryFaceList(*this, neiLevel);
75 
76 
77  bitSet seedFace;
78 
79  while (true)
80  {
81  // Pick up faces on border of protected cells
82  seedFace.reset();
83  seedFace.resize(nFaces());
84 
85  for (label facei = 0; facei < nInternalFaces(); ++facei)
86  {
87  const label own = faceOwner()[facei];
88  const label nei = faceNeighbour()[facei];
89 
90  if
91  (
92  // Protected owner
93  (
94  unrefineableCell.test(own)
95  && (cellLevel[nei] > cellLevel[own])
96  )
97  ||
98  // Protected neighbour
99  (
100  unrefineableCell.test(nei)
101  && (cellLevel[own] > cellLevel[nei])
102  )
103  )
104  {
105  seedFace.set(facei);
106  }
107  }
108  for (label facei = nInternalFaces(); facei < nFaces(); facei++)
109  {
110  const label own = faceOwner()[facei];
111 
112  if
113  (
114  // Protected owner
115  (
116  unrefineableCell.test(own)
117  && (neiLevel[facei-nInternalFaces()] > cellLevel[own])
118  )
119  )
120  {
121  seedFace.set(facei);
122  }
123  }
124 
125  syncTools::syncFaceList(*this, seedFace, orEqOp<unsigned int>());
126 
127 
128  // Extend unrefineableCell
129  bool hasExtended = false;
130 
131  for (label facei = 0; facei < nInternalFaces(); ++facei)
132  {
133  if (seedFace.test(facei))
134  {
135  if (unrefineableCell.set(faceOwner()[facei]))
136  {
137  hasExtended = true;
138  }
139  if (unrefineableCell.set(faceNeighbour()[facei]))
140  {
141  hasExtended = true;
142  }
143  }
144  }
145  for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
146  {
147  if (seedFace.test(facei))
148  {
149  const label own = faceOwner()[facei];
150 
151  if (unrefineableCell.set(own))
152  {
153  hasExtended = true;
154  }
155  }
156  }
157 
158  if (!returnReduce(hasExtended, orOp<bool>()))
159  {
160  break;
161  }
162  }
163 }
164 
165 
167 {
168  dictionary refineDict
169  (
171  (
172  IOobject
173  (
174  "dynamicMeshDict",
175  time().constant(),
176  *this,
179  false
180  )
181  ).optionalSubDict(typeName + "Coeffs")
182  );
183 
184  auto fluxVelocities = refineDict.get<List<Pair<word>>>("correctFluxes");
185 
186  // Rework into hashtable.
187  correctFluxes_.resize(fluxVelocities.size());
188  for (const auto& pr : fluxVelocities)
189  {
190  correctFluxes_.insert(pr.first(), pr.second());
191  }
192 
193  refineDict.readEntry("dumpLevel", dumpLevel_);
194 }
195 
196 
198 {
200 
201  // Correct the flux for modified/added faces. All the faces which only
202  // have been renumbered will already have been handled by the mapping.
203  {
204  const labelList& faceMap = mpm.faceMap();
205  const labelList& reverseFaceMap = mpm.reverseFaceMap();
206 
207  // Storage for any master faces. These will be the original faces
208  // on the coarse cell that get split into four (or rather the
209  // master face gets modified and three faces get added from the master)
210  // Estimate number of faces created
211 
212  bitSet masterFaces(nFaces());
213 
214  forAll(faceMap, facei)
215  {
216  const label oldFacei = faceMap[facei];
217 
218  if (oldFacei >= 0)
219  {
220  const label masterFacei = reverseFaceMap[oldFacei];
221 
222  if (masterFacei < 0)
223  {
225  << "Problem: should not have removed faces"
226  << " when refining."
227  << nl << "face:" << facei << endl
228  << abort(FatalError);
229  }
230  else if (masterFacei != facei)
231  {
232  masterFaces.set(masterFacei);
233  }
234  }
235  }
236 
237  if (debug)
238  {
239  Pout<< "Found " << masterFaces.count() << " split faces " << endl;
240  }
241 
243  (
244  lookupClass<surfaceScalarField>()
245  );
246 
247  // Remove surfaceInterpolation to allow re-calculation on demand
248  // This could be done in fvMesh::updateMesh but some dynamicFvMesh
249  // might need the old interpolation fields (weights, etc).
251 
252  forAllIters(fluxes, iter)
253  {
254  if (!correctFluxes_.found(iter.key()))
255  {
257  << "Cannot find surfaceScalarField " << iter.key()
258  << " in user-provided flux mapping table "
259  << correctFluxes_ << endl
260  << " The flux mapping table is used to recreate the"
261  << " flux on newly created faces." << endl
262  << " Either add the entry if it is a flux or use ("
263  << iter.key() << " none) to suppress this warning."
264  << endl;
265  continue;
266  }
267 
268  const word& UName = correctFluxes_[iter.key()];
269 
270  if (UName == "none")
271  {
272  continue;
273  }
274 
275  surfaceScalarField& phi = *iter();
276 
277  if (UName == "NaN")
278  {
279  Pout<< "Setting surfaceScalarField " << iter.key()
280  << " to NaN" << endl;
281 
283 
284  continue;
285  }
286 
287  if (debug)
288  {
289  Pout<< "Mapping flux " << iter.key()
290  << " using interpolated flux " << UName
291  << endl;
292  }
293 
294  const surfaceScalarField phiU
295  (
297  (
298  lookupObject<volVectorField>(UName)
299  )
300  & Sf()
301  );
302 
303  // Recalculate new internal faces.
304  for (label facei = 0; facei < nInternalFaces(); ++facei)
305  {
306  const label oldFacei = faceMap[facei];
307 
308  if (oldFacei == -1)
309  {
310  // Inflated/appended
311  phi[facei] = phiU[facei];
312  }
313  else if (reverseFaceMap[oldFacei] != facei)
314  {
315  // face-from-masterface
316  phi[facei] = phiU[facei];
317  }
318  }
319 
320  // Recalculate new boundary faces.
321  surfaceScalarField::Boundary& phiBf = phi.boundaryFieldRef();
322 
323  forAll(phiBf, patchi)
324  {
325  fvsPatchScalarField& patchPhi = phiBf[patchi];
326  const fvsPatchScalarField& patchPhiU =
327  phiU.boundaryField()[patchi];
328 
329  label facei = patchPhi.patch().start();
330 
331  forAll(patchPhi, i)
332  {
333  const label oldFacei = faceMap[facei];
334 
335  if (oldFacei == -1)
336  {
337  // Inflated/appended
338  patchPhi[i] = patchPhiU[i];
339  }
340  else if (reverseFaceMap[oldFacei] != facei)
341  {
342  // face-from-masterface
343  patchPhi[i] = patchPhiU[i];
344  }
345 
346  ++facei;
347  }
348  }
349 
350  // Update master faces
351  for (const label facei : masterFaces)
352  {
353  if (isInternalFace(facei))
354  {
355  phi[facei] = phiU[facei];
356  }
357  else
358  {
359  const label patchi = boundaryMesh().whichPatch(facei);
360  const label i = facei - boundaryMesh()[patchi].start();
361 
362  const fvsPatchScalarField& patchPhiU =
363  phiU.boundaryField()[patchi];
364 
365  fvsPatchScalarField& patchPhi = phiBf[patchi];
366 
367  patchPhi[i] = patchPhiU[i];
368  }
369  }
370  }
371  }
372 
373  // Correct the flux for injected faces - these are the faces which have
374  // no correspondence to the old mesh (i.e. added without a masterFace, edge
375  // or point). An example is the internal faces from hexRef8.
376  {
377  const labelList& faceMap = mpm.faceMap();
378 
379  mapNewInternalFaces<scalar>(this->Sf(), this->magSf(), faceMap);
380  mapNewInternalFaces<vector>(this->Sf(), this->magSf(), faceMap);
381 
382  // No oriented fields of more complex type
383  mapNewInternalFaces<sphericalTensor>(faceMap);
384  mapNewInternalFaces<symmTensor>(faceMap);
385  mapNewInternalFaces<tensor>(faceMap);
386  }
387 }
388 
389 
390 // Refines cells, maps fields and recalculates (an approximate) flux
393 (
394  const labelList& cellsToRefine
395 )
396 {
397  // Mesh changing engine.
398  polyTopoChange meshMod(*this);
399 
400  // Play refinement commands into mesh changer.
401  meshCutter_.setRefinement(cellsToRefine, meshMod);
402 
403  // Create mesh (with inflation), return map from old to new mesh.
404  //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
405  autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
406 
407  Info<< "Refined from "
408  << returnReduce(map().nOldCells(), sumOp<label>())
409  << " to " << globalData().nTotalCells() << " cells." << endl;
410 
411  if (debug)
412  {
413  // Check map.
414  for (label facei = 0; facei < nInternalFaces(); ++facei)
415  {
416  const label oldFacei = map().faceMap()[facei];
417 
418  if (oldFacei >= nInternalFaces())
419  {
421  << "New internal face:" << facei
422  << " fc:" << faceCentres()[facei]
423  << " originates from boundary oldFace:" << oldFacei
424  << abort(FatalError);
425  }
426  }
427  }
428 
429  // // Remove the stored tet base points
430  // tetBasePtIsPtr_.clear();
431  // // Remove the cell tree
432  // cellTreePtr_.clear();
433 
434  // Update fields
435  updateMesh(*map);
436 
437 
438  // Move mesh
439  /*
440  pointField newPoints;
441  if (map().hasMotionPoints())
442  {
443  newPoints = map().preMotionPoints();
444  }
445  else
446  {
447  newPoints = points();
448  }
449  movePoints(newPoints);
450  */
451 
452 
453 
454  // Update numbering of cells/vertices.
455  meshCutter_.updateMesh(*map);
456 
457  // Update numbering of protectedCell_
458  if (protectedCell_.size())
459  {
460  bitSet newProtectedCell(nCells());
461 
462  forAll(newProtectedCell, celli)
463  {
464  const label oldCelli = map().cellMap()[celli];
465  if (protectedCell_.test(oldCelli))
466  {
467  newProtectedCell.set(celli);
468  }
469  }
470  protectedCell_.transfer(newProtectedCell);
471  }
472 
473  // Debug: Check refinement levels (across faces only)
474  meshCutter_.checkRefinementLevels(-1, labelList());
475 
476  return map;
477 }
478 
479 
482 (
483  const labelList& splitPoints
484 )
485 {
486  polyTopoChange meshMod(*this);
487 
488  // Play refinement commands into mesh changer.
489  meshCutter_.setUnrefinement(splitPoints, meshMod);
490 
491 
492  // Save information on faces that will be combined
493  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
494 
495  // Find the faceMidPoints on cells to be combined.
496  // for each face resulting of split of face into four store the
497  // midpoint
498  Map<label> faceToSplitPoint(3*splitPoints.size());
499 
500  {
501  for (const label pointi : splitPoints)
502  {
503  const labelList& pEdges = pointEdges()[pointi];
504 
505  for (const label edgei : pEdges)
506  {
507  const label otherPointi = edges()[edgei].otherVertex(pointi);
508 
509  const labelList& pFaces = pointFaces()[otherPointi];
510 
511  for (const label facei : pFaces)
512  {
513  faceToSplitPoint.insert(facei, otherPointi);
514  }
515  }
516  }
517  }
518 
519 
520  // Change mesh and generate map.
521  //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
522  autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
523 
524  Info<< "Unrefined from "
525  << returnReduce(map().nOldCells(), sumOp<label>())
526  << " to " << globalData().nTotalCells() << " cells."
527  << endl;
528 
529  // Update fields
530  updateMesh(*map);
531 
532 
533  // Move mesh
534  /*
535  pointField newPoints;
536  if (map().hasMotionPoints())
537  {
538  newPoints = map().preMotionPoints();
539  }
540  else
541  {
542  newPoints = points();
543  }
544  movePoints(newPoints);
545  */
546 
547  // Correct the flux for modified faces.
548  {
549  const labelList& reversePointMap = map().reversePointMap();
550  const labelList& reverseFaceMap = map().reverseFaceMap();
551 
553  (
554  lookupClass<surfaceScalarField>()
555  );
556  forAllIters(fluxes, iter)
557  {
558  if (!correctFluxes_.found(iter.key()))
559  {
561  << "Cannot find surfaceScalarField " << iter.key()
562  << " in user-provided flux mapping table "
563  << correctFluxes_ << endl
564  << " The flux mapping table is used to recreate the"
565  << " flux on newly created faces." << endl
566  << " Either add the entry if it is a flux or use ("
567  << iter.key() << " none) to suppress this warning."
568  << endl;
569  continue;
570  }
571 
572  const word& UName = correctFluxes_[iter.key()];
573 
574  if (UName == "none")
575  {
576  continue;
577  }
578 
579  DebugInfo
580  << "Mapping flux " << iter.key()
581  << " using interpolated flux " << UName
582  << endl;
583 
584 
585  surfaceScalarField& phi = *iter();
586  surfaceScalarField::Boundary& phiBf =
588 
589  const surfaceScalarField phiU
590  (
592  (
593  lookupObject<volVectorField>(UName)
594  )
595  & Sf()
596  );
597 
598 
599  forAllConstIters(faceToSplitPoint, iter)
600  {
601  const label oldFacei = iter.key();
602  const label oldPointi = iter.val();
603 
604  if (reversePointMap[oldPointi] < 0)
605  {
606  // midpoint was removed. See if face still exists.
607  const label facei = reverseFaceMap[oldFacei];
608 
609  if (facei >= 0)
610  {
611  if (isInternalFace(facei))
612  {
613  phi[facei] = phiU[facei];
614  }
615  else
616  {
617  label patchi = boundaryMesh().whichPatch(facei);
618  label i = facei - boundaryMesh()[patchi].start();
619 
620  const fvsPatchScalarField& patchPhiU =
621  phiU.boundaryField()[patchi];
622  fvsPatchScalarField& patchPhi = phiBf[patchi];
623  patchPhi[i] = patchPhiU[i];
624  }
625  }
626  }
627  }
628  }
629  }
630 
631 
632  // Update numbering of cells/vertices.
633  meshCutter_.updateMesh(*map);
634 
635  // Update numbering of protectedCell_
636  if (protectedCell_.size())
637  {
638  bitSet newProtectedCell(nCells());
639 
640  forAll(newProtectedCell, celli)
641  {
642  const label oldCelli = map().cellMap()[celli];
643  if (protectedCell_.test(oldCelli))
644  {
645  newProtectedCell.set(celli);
646  }
647  }
648  protectedCell_.transfer(newProtectedCell);
649  }
650 
651  // Debug: Check refinement levels (across faces only)
652  meshCutter_.checkRefinementLevels(-1, labelList());
653 
654  return map;
655 }
656 
657 
660 {
661  scalarField vFld(nCells(), -GREAT);
662 
663  forAll(pointCells(), pointi)
664  {
665  const labelList& pCells = pointCells()[pointi];
666 
667  for (const label celli : pCells)
668  {
669  vFld[celli] = max(vFld[celli], pFld[pointi]);
670  }
671  }
672  return vFld;
673 }
674 
675 
678 {
679  scalarField pFld(nPoints(), -GREAT);
680 
681  forAll(pointCells(), pointi)
682  {
683  const labelList& pCells = pointCells()[pointi];
684 
685  for (const label celli : pCells)
686  {
687  pFld[pointi] = max(pFld[pointi], vFld[celli]);
688  }
689  }
690  return pFld;
691 }
692 
693 
696 {
697  scalarField pFld(nPoints());
698 
699  forAll(pointCells(), pointi)
700  {
701  const labelList& pCells = pointCells()[pointi];
702 
703  scalar sum = 0.0;
704  for (const label celli : pCells)
705  {
706  sum += vFld[celli];
707  }
708  pFld[pointi] = sum/pCells.size();
709  }
710  return pFld;
711 }
712 
713 
715 (
716  const scalarField& fld,
717  const scalar minLevel,
718  const scalar maxLevel
719 ) const
720 {
721  scalarField c(fld.size(), scalar(-1));
722 
723  forAll(fld, i)
724  {
725  scalar err = min(fld[i]-minLevel, maxLevel-fld[i]);
726 
727  if (err >= 0)
728  {
729  c[i] = err;
730  }
731  }
732  return c;
733 }
734 
735 
737 (
738  const scalar lowerRefineLevel,
739  const scalar upperRefineLevel,
740  const scalarField& vFld,
741  bitSet& candidateCell
742 ) const
743 {
744  // Get error per cell. Is -1 (not to be refined) to >0 (to be refined,
745  // higher more desirable to be refined).
746  scalarField cellError
747  (
748  maxPointField
749  (
750  error
751  (
752  cellToPoint(vFld),
753  lowerRefineLevel,
754  upperRefineLevel
755  )
756  )
757  );
758 
759  // Mark cells that are candidates for refinement.
760  forAll(cellError, celli)
761  {
762  if (cellError[celli] > 0)
763  {
764  candidateCell.set(celli);
765  }
766  }
767 }
768 
769 
771 (
772  const label maxCells,
773  const label maxRefinement,
774  const bitSet& candidateCell
775 ) const
776 {
777  // Every refined cell causes 7 extra cells
778  label nTotToRefine = (maxCells - globalData().nTotalCells()) / 7;
779 
780  const labelList& cellLevel = meshCutter_.cellLevel();
781 
782  // Mark cells that cannot be refined since they would trigger refinement
783  // of protected cells (since 2:1 cascade)
784  bitSet unrefineableCell;
785  calculateProtectedCells(unrefineableCell);
786 
787  // Count current selection
788  label nLocalCandidates = candidateCell.count();
789  label nCandidates = returnReduce(nLocalCandidates, sumOp<label>());
790 
791  // Collect all cells
792  DynamicList<label> candidates(nLocalCandidates);
793 
794  if (nCandidates < nTotToRefine)
795  {
796  for (const label celli : candidateCell)
797  {
798  if
799  (
800  (!unrefineableCell.test(celli))
801  && cellLevel[celli] < maxRefinement
802  )
803  {
804  candidates.append(celli);
805  }
806  }
807  }
808  else
809  {
810  // Sort by error? For now just truncate.
811  for (label level = 0; level < maxRefinement; ++level)
812  {
813  for (const label celli : candidateCell)
814  {
815  if
816  (
817  (!unrefineableCell.test(celli))
818  && cellLevel[celli] == level
819  )
820  {
821  candidates.append(celli);
822  }
823  }
824 
825  if (returnReduce(candidates.size(), sumOp<label>()) > nTotToRefine)
826  {
827  break;
828  }
829  }
830  }
831 
832  // Guarantee 2:1 refinement after refinement
833  labelList consistentSet
834  (
835  meshCutter_.consistentRefinement
836  (
837  candidates.shrink(),
838  true // Add to set to guarantee 2:1
839  )
840  );
841 
842  Info<< "Selected " << returnReduce(consistentSet.size(), sumOp<label>())
843  << " cells for refinement out of " << globalData().nTotalCells()
844  << "." << endl;
845 
846  return consistentSet;
847 }
848 
849 
851 (
852  const scalar unrefineLevel,
853  const bitSet& markedCell,
854  const scalarField& pFld
855 ) const
856 {
857  // All points that can be unrefined
858  const labelList splitPoints(meshCutter_.getSplitPoints());
859 
860 
861  const labelListList& pointCells = this->pointCells();
862 
863  // If we have any protected cells make sure they also are not being
864  // unrefined
865 
866  bitSet protectedPoint(nPoints());
867 
868  if (protectedCell_.size())
869  {
870  // Get all points on a protected cell
871  forAll(pointCells, pointi)
872  {
873  for (const label celli : pointCells[pointi])
874  {
875  if (protectedCell_.test(celli))
876  {
877  protectedPoint.set(pointi);
878  break;
879  }
880  }
881  }
882 
884  (
885  *this,
886  protectedPoint,
888  0u
889  );
890 
891  DebugInfo<< "From "
892  << returnReduce(protectedCell_.count(), sumOp<label>())
893  << " protected cells found "
894  << returnReduce(protectedPoint.count(), sumOp<label>())
895  << " protected points." << endl;
896  }
897 
898 
899  DynamicList<label> newSplitPoints(splitPoints.size());
900 
901  for (const label pointi : splitPoints)
902  {
903  if (!protectedPoint[pointi] && pFld[pointi] < unrefineLevel)
904  {
905  // Check that all cells are not marked
906  bool hasMarked = false;
907 
908  for (const label celli : pointCells[pointi])
909  {
910  if (markedCell.test(celli))
911  {
912  hasMarked = true;
913  break;
914  }
915  }
916 
917  if (!hasMarked)
918  {
919  newSplitPoints.append(pointi);
920  }
921  }
922  }
923 
924 
925  newSplitPoints.shrink();
926 
927  // Guarantee 2:1 refinement after unrefinement
928  labelList consistentSet
929  (
930  meshCutter_.consistentUnrefinement
931  (
932  newSplitPoints,
933  false
934  )
935  );
936  Info<< "Selected " << returnReduce(consistentSet.size(), sumOp<label>())
937  << " split points out of a possible "
938  << returnReduce(splitPoints.size(), sumOp<label>())
939  << "." << endl;
940 
941  return consistentSet;
942 }
943 
944 
946 (
947  bitSet& markedCell
948 ) const
949 {
950  // Mark faces using any marked cell
951  bitSet markedFace(nFaces());
952 
953  for (const label celli : markedCell)
954  {
955  markedFace.set(cells()[celli]); // set multiple faces
956  }
957 
958  syncTools::syncFaceList(*this, markedFace, orEqOp<unsigned int>());
959 
960  // Update cells using any markedFace
961  for (label facei = 0; facei < nInternalFaces(); ++facei)
962  {
963  if (markedFace.test(facei))
964  {
965  markedCell.set(faceOwner()[facei]);
966  markedCell.set(faceNeighbour()[facei]);
967  }
968  }
969  for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
970  {
971  if (markedFace.test(facei))
972  {
973  markedCell.set(faceOwner()[facei]);
974  }
975  }
976 }
977 
978 
980 (
981  bitSet& protectedCell
982 ) const
983 {
984  const labelList& cellLevel = meshCutter_.cellLevel();
985  const labelList& pointLevel = meshCutter_.pointLevel();
986 
987  labelList nAnchorPoints(nCells(), Zero);
988 
989  forAll(pointLevel, pointi)
990  {
991  const labelList& pCells = pointCells(pointi);
992 
993  for (const label celli : pCells)
994  {
995  if (pointLevel[pointi] <= cellLevel[celli])
996  {
997  // Check if cell has already 8 anchor points -> protect cell
998  if (nAnchorPoints[celli] == 8)
999  {
1000  protectedCell.set(celli);
1001  }
1002 
1003  if (!protectedCell.test(celli))
1004  {
1005  ++nAnchorPoints[celli];
1006  }
1007  }
1008  }
1009  }
1010 
1011 
1012  forAll(protectedCell, celli)
1013  {
1014  if (nAnchorPoints[celli] != 8)
1015  {
1016  protectedCell.set(celli);
1017  }
1018  }
1019 }
1020 
1021 
1022 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1023 
1024 Foam::dynamicRefineFvMesh::dynamicRefineFvMesh
1026  const IOobject& io,
1027  const bool doInit
1028 )
1029 :
1030  dynamicFvMesh(io, doInit),
1031  meshCutter_(*this)
1032 {
1033  if (doInit)
1034  {
1035  init(false); // do not initialise lower levels
1036  }
1037 }
1038 
1039 
1040 bool Foam::dynamicRefineFvMesh::init(const bool doInit)
1041 {
1042  if (doInit)
1043  {
1044  dynamicFvMesh::init(doInit);
1045  }
1046 
1047  protectedCell_.setSize(nCells());
1048  nRefinementIterations_ = 0;
1049  dumpLevel_ = false;
1050 
1051  // Read static part of dictionary
1052  readDict();
1053 
1054 
1055  const labelList& cellLevel = meshCutter_.cellLevel();
1056  const labelList& pointLevel = meshCutter_.pointLevel();
1057 
1058  // Set cells that should not be refined.
1059  // This is currently any cell which does not have 8 anchor points or
1060  // uses any face which does not have 4 anchor points.
1061  // Note: do not use cellPoint addressing
1062 
1063  // Count number of points <= cellLevel
1064  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1065 
1066  labelList nAnchors(nCells(), Zero);
1067 
1068  forAll(pointCells(), pointi)
1069  {
1070  const labelList& pCells = pointCells()[pointi];
1071 
1072  for (const label celli : pCells)
1073  {
1074  if (!protectedCell_.test(celli))
1075  {
1076  if (pointLevel[pointi] <= cellLevel[celli])
1077  {
1078  ++nAnchors[celli];
1079 
1080  if (nAnchors[celli] > 8)
1081  {
1082  protectedCell_.set(celli);
1083  }
1084  }
1085  }
1086  }
1087  }
1088 
1089 
1090  // Count number of points <= faceLevel
1091  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1092  // Bit tricky since proc face might be one more refined than the owner since
1093  // the coupled one is refined.
1094 
1095  {
1096  labelList neiLevel(nFaces());
1097 
1098  for (label facei = 0; facei < nInternalFaces(); ++facei)
1099  {
1100  neiLevel[facei] = cellLevel[faceNeighbour()[facei]];
1101  }
1102  for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
1103  {
1104  neiLevel[facei] = cellLevel[faceOwner()[facei]];
1105  }
1106  syncTools::swapFaceList(*this, neiLevel);
1107 
1108 
1109  bitSet protectedFace(nFaces());
1110 
1111  forAll(faceOwner(), facei)
1112  {
1113  const label faceLevel = max
1114  (
1115  cellLevel[faceOwner()[facei]],
1116  neiLevel[facei]
1117  );
1118 
1119  const face& f = faces()[facei];
1120 
1121  label nAnchors = 0;
1122 
1123  for (const label pointi : f)
1124  {
1125  if (pointLevel[pointi] <= faceLevel)
1126  {
1127  ++nAnchors;
1128 
1129  if (nAnchors > 4)
1130  {
1131  protectedFace.set(facei);
1132  break;
1133  }
1134  }
1135  }
1136  }
1137 
1138  syncTools::syncFaceList(*this, protectedFace, orEqOp<unsigned int>());
1139 
1140  for (label facei = 0; facei < nInternalFaces(); ++facei)
1141  {
1142  if (protectedFace.test(facei))
1143  {
1144  protectedCell_.set(faceOwner()[facei]);
1145  protectedCell_.set(faceNeighbour()[facei]);
1146  }
1147  }
1148  for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
1149  {
1150  if (protectedFace.test(facei))
1151  {
1152  protectedCell_.set(faceOwner()[facei]);
1153  }
1154  }
1155 
1156  // Also protect any cells that are less than hex
1157  forAll(cells(), celli)
1158  {
1159  const cell& cFaces = cells()[celli];
1160 
1161  if (cFaces.size() < 6)
1162  {
1163  protectedCell_.set(celli);
1164  }
1165  else
1166  {
1167  for (const label cfacei : cFaces)
1168  {
1169  if (faces()[cfacei].size() < 4)
1170  {
1171  protectedCell_.set(celli);
1172  break;
1173  }
1174  }
1175  }
1176  }
1177 
1178  // Check cells for 8 corner points
1179  checkEightAnchorPoints(protectedCell_);
1180  }
1181 
1182  if (!returnReduce(protectedCell_.any(), orOp<bool>()))
1183  {
1184  protectedCell_.clear();
1185  }
1186  else
1187  {
1188  cellSet protectedCells
1189  (
1190  *this,
1191  "protectedCells",
1192  HashSetOps::used(protectedCell_)
1193  );
1194 
1195  Info<< "Detected "
1196  << returnReduce(protectedCells.size(), sumOp<label>())
1197  << " cells that are protected from refinement."
1198  << " Writing these to cellSet "
1199  << protectedCells.name()
1200  << "." << endl;
1201 
1202  protectedCells.write();
1203  }
1204 
1205  return true;
1206 }
1207 
1208 
1209 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1210 
1212 {
1213  // Re-read dictionary. Chosen since usually -small so trivial amount
1214  // of time compared to actual refinement. Also very useful to be able
1215  // to modify on-the-fly.
1216  dictionary refineDict
1217  (
1218  IOdictionary
1219  (
1220  IOobject
1221  (
1222  "dynamicMeshDict",
1223  time().constant(),
1224  *this,
1227  false
1228  )
1229  ).optionalSubDict(typeName + "Coeffs")
1230  );
1231 
1232  const label refineInterval = refineDict.get<label>("refineInterval");
1233 
1234  bool hasChanged = false;
1235 
1236  if (refineInterval == 0)
1237  {
1238  topoChanging(hasChanged);
1239 
1240  return false;
1241  }
1242  else if (refineInterval < 0)
1243  {
1245  << "Illegal refineInterval " << refineInterval << nl
1246  << "The refineInterval setting in the dynamicMeshDict should"
1247  << " be >= 1." << nl
1248  << exit(FatalError);
1249  }
1250 
1251 
1252  // Note: cannot refine at time 0 since no V0 present since mesh not
1253  // moved yet.
1254 
1255  if (time().timeIndex() > 0 && time().timeIndex() % refineInterval == 0)
1256  {
1257  const label maxCells = refineDict.get<label>("maxCells");
1258 
1259  if (maxCells <= 0)
1260  {
1262  << "Illegal maximum number of cells " << maxCells << nl
1263  << "The maxCells setting in the dynamicMeshDict should"
1264  << " be > 0." << nl
1265  << exit(FatalError);
1266  }
1267 
1268  const label maxRefinement = refineDict.get<label>("maxRefinement");
1269 
1270  if (maxRefinement <= 0)
1271  {
1273  << "Illegal maximum refinement level " << maxRefinement << nl
1274  << "The maxCells setting in the dynamicMeshDict should"
1275  << " be > 0." << nl
1276  << exit(FatalError);
1277  }
1278 
1279  const word fieldName(refineDict.get<word>("field"));
1280 
1281  const volScalarField& vFld = lookupObject<volScalarField>(fieldName);
1282 
1283  const scalar lowerRefineLevel =
1284  refineDict.get<scalar>("lowerRefineLevel");
1285  const scalar upperRefineLevel =
1286  refineDict.get<scalar>("upperRefineLevel");
1287  const scalar unrefineLevel = refineDict.getOrDefault<scalar>
1288  (
1289  "unrefineLevel",
1290  GREAT
1291  );
1292  const label nBufferLayers = refineDict.get<label>("nBufferLayers");
1293 
1294  // Cells marked for refinement or otherwise protected from unrefinement.
1295  bitSet refineCell(nCells());
1296 
1297  // Determine candidates for refinement (looking at field only)
1298  selectRefineCandidates
1299  (
1300  lowerRefineLevel,
1301  upperRefineLevel,
1302  vFld,
1303  refineCell
1304  );
1305 
1306  if (globalData().nTotalCells() < maxCells)
1307  {
1308  // Select subset of candidates. Take into account max allowable
1309  // cells, refinement level, protected cells.
1310  labelList cellsToRefine
1311  (
1312  selectRefineCells
1313  (
1314  maxCells,
1315  maxRefinement,
1316  refineCell
1317  )
1318  );
1319 
1320  const label nCellsToRefine = returnReduce
1321  (
1322  cellsToRefine.size(), sumOp<label>()
1323  );
1324 
1325  if (nCellsToRefine > 0)
1326  {
1327  // Refine/update mesh and map fields
1328  autoPtr<mapPolyMesh> map = refine(cellsToRefine);
1329 
1330  // Update refineCell. Note that some of the marked ones have
1331  // not been refined due to constraints.
1332  {
1333  const labelList& cellMap = map().cellMap();
1334  const labelList& reverseCellMap = map().reverseCellMap();
1335 
1336  bitSet newRefineCell(cellMap.size());
1337 
1338  forAll(cellMap, celli)
1339  {
1340  const label oldCelli = cellMap[celli];
1341 
1342  if
1343  (
1344  (oldCelli < 0)
1345  || (reverseCellMap[oldCelli] != celli)
1346  || (refineCell.test(oldCelli))
1347  )
1348  {
1349  newRefineCell.set(celli);
1350  }
1351  }
1352  refineCell.transfer(newRefineCell);
1353  }
1354 
1355  // Extend with a buffer layer to prevent neighbouring points
1356  // being unrefined.
1357  for (label i = 0; i < nBufferLayers; ++i)
1358  {
1359  extendMarkedCells(refineCell);
1360  }
1361 
1362  hasChanged = true;
1363  }
1364  }
1365 
1366 
1367  {
1368  // Select unrefineable points that are not marked in refineCell
1369  labelList pointsToUnrefine
1370  (
1371  selectUnrefinePoints
1372  (
1373  unrefineLevel,
1374  refineCell,
1375  maxCellField(vFld)
1376  )
1377  );
1378 
1379  const label nSplitPoints = returnReduce
1380  (
1381  pointsToUnrefine.size(),
1382  sumOp<label>()
1383  );
1384 
1385  if (nSplitPoints > 0)
1386  {
1387  // Refine/update mesh
1388  unrefine(pointsToUnrefine);
1389 
1390  hasChanged = true;
1391  }
1392  }
1393 
1394 
1395  if ((nRefinementIterations_ % 10) == 0)
1396  {
1397  // Compact refinement history occasionally (how often?).
1398  // Unrefinement causes holes in the refinementHistory.
1399  const_cast<refinementHistory&>(meshCutter().history()).compact();
1400  }
1401  nRefinementIterations_++;
1402  }
1403 
1404  topoChanging(hasChanged);
1405  if (hasChanged)
1406  {
1407  // Reset moving flag (if any). If not using inflation we'll not move,
1408  // if are using inflation any follow on movePoints will set it.
1409  moving(false);
1410  }
1411 
1412  return hasChanged;
1413 }
1414 
1415 
1418  IOstreamOption streamOpt,
1419  const bool valid
1420 ) const
1421 {
1422  // Force refinement data to go to the current time directory.
1423  const_cast<hexRef8&>(meshCutter_).setInstance(time().timeName());
1424 
1425  bool writeOk =
1426  (
1427  dynamicFvMesh::writeObject(streamOpt, valid)
1428  && meshCutter_.write(valid)
1429  );
1430 
1431  if (dumpLevel_)
1432  {
1433  volScalarField scalarCellLevel
1434  (
1435  IOobject
1436  (
1437  "cellLevel",
1438  time().timeName(),
1439  *this,
1442  false
1443  ),
1444  *this,
1446  );
1447 
1448  const labelList& cellLevel = meshCutter_.cellLevel();
1449 
1450  forAll(cellLevel, celli)
1451  {
1452  scalarCellLevel[celli] = cellLevel[celli];
1453  }
1454 
1455  writeOk = writeOk && scalarCellLevel.write();
1456  }
1457 
1458  return writeOk;
1459 }
1460 
1461 
1462 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
Foam::IOobject::NO_WRITE
Definition: IOobject.H:195
Foam::surfaceInterpolation::clearOut
void clearOut()
Clear all geometry and addressing.
Definition: surfaceInterpolation.C:49
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:67
volFields.H
Foam::IOdictionary
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:54
Foam::PackedList::reset
void reset()
Clear all bits but do not adjust the addressable size.
Definition: PackedListI.H:505
Foam::PackedList::resize
void resize(const label numElem, const unsigned int val=0u)
Reset addressable list size, does not shrink the allocated size.
Definition: PackedListI.H:409
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:169
Foam::meshCutter
Cuts (splits) cells.
Definition: meshCutter.H:138
Foam::IOobject::AUTO_WRITE
Definition: IOobject.H:194
Foam::dynamicRefineFvMesh::error
scalarField error(const scalarField &fld, const scalar minLevel, const scalar maxLevel) const
Definition: dynamicRefineFvMesh.C:715
Foam::faceMap
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
Definition: blockMeshMergeTopological.C:94
Foam::dynamicRefineFvMesh::extendMarkedCells
void extendMarkedCells(bitSet &markedCell) const
Extend markedCell with cell-face-cell.
Definition: dynamicRefineFvMesh.C:946
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
Foam::returnReduce
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Definition: PstreamReduceOps.H:94
Foam::bitSet
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:63
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::DynamicList< label >
Foam::fvPatch::start
virtual label start() const
Return start label of this patch in the polyMesh face list.
Definition: fvPatch.H:173
Foam::cellToPoint
A topoSetPointSource to select all the points from given cellSet(s).
Definition: cellToPoint.H:173
UName
const word UName(propsDict.getOrDefault< word >("U", "U"))
Foam::dynamicRefineFvMesh::init
virtual bool init(const bool doInit)
Initialise all non-demand-driven data.
Definition: dynamicRefineFvMesh.C:1040
polyTopoChange.H
Foam::fvMesh::mapFields
virtual void mapFields(const mapPolyMesh &mpm)
Map all fields in time using given map.
Definition: fvMesh.C:714
Foam::fvsPatchField
An abstract base class with a fat-interface to all derived classes covering all possible ways in whic...
Definition: fvsPatchField.H:68
Foam::polyTopoChange
Direct mesh changes based on v1.3 polyTopoChange syntax.
Definition: polyTopoChange.H:99
Foam::dynamicFvMesh::init
virtual bool init(const bool doInit)
Initialise all non-demand-driven data.
Definition: dynamicFvMesh.C:90
Foam::dynamicFvMesh
Abstract base class for geometry and/or topology changing fvMesh.
Definition: dynamicFvMesh.H:78
Foam::dynamicRefineFvMesh::mapFields
virtual void mapFields(const mapPolyMesh &mpm)
Map all fields in time using given map.
Definition: dynamicRefineFvMesh.C:197
Foam::dynamicRefineFvMesh::update
virtual bool update()
Update the mesh for both mesh motion and topology change.
Definition: dynamicRefineFvMesh.C:1211
Foam::dynamicRefineFvMesh::selectRefineCells
virtual labelList selectRefineCells(const label maxCells, const label maxRefinement, const bitSet &candidateCell) const
Subset candidate cells for refinement.
Definition: dynamicRefineFvMesh.C:771
Foam::bitSet::set
void set(const bitSet &bitset)
Set specified bits from another bitset.
Definition: bitSetI.H:574
Foam::Map< label >
Foam::orEqOp
Definition: ops.H:86
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
surfaceFields.H
Foam::surfaceFields.
Foam::dynamicRefineFvMesh::correctFluxes_
HashTable< word > correctFluxes_
Fluxes to map.
Definition: dynamicRefineFvMesh.H:100
Foam::dictionary::get
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:107
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::bitSet::test
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
Definition: bitSetI.H:520
Foam::DynamicList::shrink
DynamicList< T, SizeMin > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:434
syncTools.H
Foam::bitSet::count
unsigned int count(const bool on=true) const
Count number of bits set.
Definition: bitSetI.H:499
Foam::dynamicRefineFvMesh::readDict
void readDict()
Read the projection parameters from dictionary.
Definition: dynamicRefineFvMesh.C:166
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::syncTools::syncPointList
static void syncPointList(const polyMesh &mesh, List< T > &pointValues, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
Definition: syncToolsTemplates.C:721
Foam::sumOp
Definition: ops.H:213
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::dynamicRefineFvMesh::writeObject
virtual bool writeObject(IOstreamOption streamOpt, const bool valid) const
Write using stream options.
Definition: dynamicRefineFvMesh.C:1417
nPoints
label nPoints
Definition: gmvOutputHeader.H:2
Foam::refinementHistory
All refinement history. Used in unrefinement.
Definition: refinementHistory.H:102
Foam::syncTools::syncFaceList
static void syncFaceList(const polyMesh &mesh, UList< T > &faceValues, const CombineOp &cop)
Synchronize values on all mesh faces.
Definition: syncTools.H:396
Foam::dynamicRefineFvMesh::checkEightAnchorPoints
void checkEightAnchorPoints(bitSet &protectedCell) const
Check all cells have 8 anchor points.
Definition: dynamicRefineFvMesh.C:980
Foam::Field< scalar >
Foam::regIOobject::write
virtual bool write(const bool valid=true) const
Write using setting from DB.
Definition: regIOobjectWrite.C:132
Foam::sigFpe::fillNan
static void fillNan(UList< scalar > &list)
Fill data block with NaN values.
Definition: sigFpe.C:225
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::DynamicList::append
DynamicList< T, SizeMin > & append(const T &val)
Append an element to the end of this list.
Definition: DynamicListI.H:511
Foam::dynamicRefineFvMesh::selectRefineCandidates
virtual void selectRefineCandidates(const scalar lowerRefineLevel, const scalar upperRefineLevel, const scalarField &vFld, bitSet &candidateCell) const
Select candidate cells for refinement.
Definition: dynamicRefineFvMesh.C:737
init
mesh init(true)
Foam::IOstreamOption
The IOstreamOption is a simple container for options an IOstream can normally have.
Definition: IOstreamOption.H:63
Foam::dictionary::readEntry
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, bool mandatory=true) const
Definition: dictionaryTemplates.C:302
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition: dimensionedScalarFwd.H:42
phi
surfaceScalarField & phi
Definition: setRegionFluidFields.H:8
Foam::dynamicRefineFvMesh::maxCellField
scalarField maxCellField(const volScalarField &) const
Get point max of connected cell.
Definition: dynamicRefineFvMesh.C:677
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
fld
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< ' ';}gmvFile<< nl;for(const word &name :lagrangianScalarNames){ IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Definition: gmvOutputLagrangian.H:23
Foam::polyTopoChange::changeMesh
autoPtr< mapPolyMesh > changeMesh(polyMesh &mesh, const labelUList &patchMap, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
Definition: polyTopoChange.C:2980
forAllIters
#define forAllIters(container, iter)
Iterate across all elements in the container object.
Definition: stdFoam.H:223
Foam::dynamicRefineFvMesh::refine
virtual autoPtr< mapPolyMesh > refine(const labelList &)
Refine cells. Update mesh and fields.
Definition: dynamicRefineFvMesh.C:393
timeName
word timeName
Definition: getTimeIndex.H:3
Foam::FatalError
error FatalError
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
addToRunTimeSelectionTable.H
Macros for easy insertion into run-time selection tables.
dynamicRefineFvMesh.H
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Foam::fvc::interpolate
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
Foam::cellSet
A collection of cell labels.
Definition: cellSet.H:51
Foam::dynamicRefineFvMesh::selectUnrefinePoints
virtual labelList selectUnrefinePoints(const scalar unrefineLevel, const bitSet &markedCell, const scalarField &pFld) const
Select points that can be unrefined.
Definition: dynamicRefineFvMesh.C:851
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Foam::GeometricField::primitiveFieldRef
Internal::FieldType & primitiveFieldRef(const bool updateAccessTime=true)
Return a reference to the internal field.
Definition: GeometricField.C:766
Foam::HashTable
A HashTable similar to std::unordered_map.
Definition: HashTable.H:105
Foam::mapPolyMesh::reverseFaceMap
const labelList & reverseFaceMap() const
Reverse face map.
Definition: mapPolyMesh.H:501
Foam::IOobject::name
const word & name() const noexcept
Return name.
Definition: IOobjectI.H:65
Foam::autoPtr
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: HashPtrTable.H:53
Foam::dynamicRefineFvMesh::calculateProtectedCells
void calculateProtectedCells(bitSet &unrefineableCell) const
Calculate cells that cannot be refined since would trigger.
Definition: dynamicRefineFvMesh.C:53
Foam::refineCell
Container with cells to refine. Refinement given as single direction.
Definition: refineCell.H:56
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
DebugInfo
#define DebugInfo
Report an information message using Foam::Info.
Definition: messageStream.H:382
Foam::nl
constexpr char nl
Definition: Ostream.H:404
forAllConstIters
forAllConstIters(mixture.phases(), phase)
Definition: pEqn.H:28
Foam::fvsPatchField::patch
const fvPatch & patch() const
Return patch.
Definition: fvsPatchField.H:281
Foam::hexRef8
Refinement of (split) hexes using polyTopoChange.
Definition: hexRef8.H:67
f
labelList f(nPoints)
Foam::GeometricField::boundaryFieldRef
Boundary & boundaryFieldRef(const bool updateAccessTime=true)
Return a reference to the boundary field.
Definition: GeometricField.C:783
Foam::syncTools::swapFaceList
static void swapFaceList(const polyMesh &mesh, UList< T > &faceValues)
Swap coupled face values. Uses eqOp.
Definition: syncTools.H:478
Foam::List< label >
Foam::HashSetOps::used
labelHashSet used(const bitSet &select)
Convert a bitset to a labelHashSet of the indices used.
Definition: HashOps.C:35
Foam::dynamicRefineFvMesh::dumpLevel_
bool dumpLevel_
Dump cellLevel for post-processing.
Definition: dynamicRefineFvMesh.H:109
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
Foam::IOobject::MUST_READ_IF_MODIFIED
Definition: IOobject.H:186
Foam::mapPolyMesh::faceMap
const labelList & faceMap() const
Old face map.
Definition: mapPolyMesh.H:410
Foam::dictionary::optionalSubDict
const dictionary & optionalSubDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary, otherwise return this dictionary.
Definition: dictionary.C:577
Foam::List::set
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
Definition: List.H:341
Foam::mapPolyMesh
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:161
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Foam::fvMesh::time
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:280
Foam::boundaryMesh
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface....
Definition: boundaryMesh.H:62
Foam::dynamicRefineFvMesh::unrefine
virtual autoPtr< mapPolyMesh > unrefine(const labelList &)
Unrefine cells. Gets passed in centre points of cells to combine.
Definition: dynamicRefineFvMesh.C:482
Foam::dynamicRefineFvMesh::cellToPoint
scalarField cellToPoint(const scalarField &vFld) const
Definition: dynamicRefineFvMesh.C:695
Foam::face
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:72
timeIndex
label timeIndex
Definition: getTimeIndex.H:30
Foam::pointCells
Smooth ATC in cells having a point to a set of patches supplied by type.
Definition: pointCells.H:56
cells
const cellShapeList & cells
Definition: gmvOutputHeader.H:3
cellSet.H
Foam::dictionary::getOrDefault
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:148
Foam::GeometricField< scalar, fvsPatchField, surfaceMesh >
Foam::IOobject::NO_READ
Definition: IOobject.H:188
Foam::orOp
Definition: ops.H:234
constant
constant condensation/saturation model.
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
Foam::dynamicRefineFvMesh::maxPointField
scalarField maxPointField(const scalarField &) const
Get per cell max of connected point.
Definition: dynamicRefineFvMesh.C:659
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
Foam::cell
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:54
Foam::error
Class to handle errors and exceptions in a simple, consistent stream-based manner.
Definition: error.H:73
pointFields.H
pFaces
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::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].reset(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
Foam::PackedList::clear
void clear()
Clear the list, i.e. set addressable size to zero.
Definition: PackedListI.H:512
HashOps.H
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Definition: GeometricFieldI.H:62
Foam::dimless
const dimensionSet dimless
Dimensionless.
Definition: dimensionSets.C:189
Foam::fvMesh::writeObject
virtual bool writeObject(IOstreamOption streamOpt, const bool valid) const
Write the underlying polyMesh and other data.
Definition: fvMesh.C:1018