28 #include "surfaceInterpolate.H" 50 const unsigned int val
78 if (protectedCell_.empty())
80 unrefineableCell.
clear();
84 const labelList& cellLevel = meshCutter_.cellLevel();
86 unrefineableCell = protectedCell_;
89 labelList neiLevel(nFaces()-nInternalFaces());
91 for (
label facei = nInternalFaces(); facei < nFaces(); facei++)
93 neiLevel[facei-nInternalFaces()] = cellLevel[faceOwner()[facei]];
103 forAll(faceNeighbour(), facei)
105 label own = faceOwner()[facei];
106 bool ownProtected = unrefineableCell.
get(own);
107 label nei = faceNeighbour()[facei];
108 bool neiProtected = unrefineableCell.
get(nei);
110 if (ownProtected && (cellLevel[nei] > cellLevel[own]))
112 seedFace[facei] =
true;
114 else if (neiProtected && (cellLevel[own] > cellLevel[nei]))
116 seedFace[facei] =
true;
119 for (
label facei = nInternalFaces(); facei < nFaces(); facei++)
121 label own = faceOwner()[facei];
122 bool ownProtected = unrefineableCell.
get(own);
126 && (neiLevel[facei-nInternalFaces()] > cellLevel[own])
129 seedFace[facei] =
true;
137 bool hasExtended =
false;
139 for (
label facei = 0; facei < nInternalFaces(); facei++)
143 label own = faceOwner()[facei];
144 if (unrefineableCell.
get(own) == 0)
146 unrefineableCell.
set(own, 1);
150 label nei = faceNeighbour()[facei];
151 if (unrefineableCell.
get(nei) == 0)
153 unrefineableCell.
set(nei, 1);
158 for (
label facei = nInternalFaces(); facei < nFaces(); facei++)
162 label own = faceOwner()[facei];
163 if (unrefineableCell.
get(own) == 0)
165 unrefineableCell.
set(own, 1);
194 ).optionalSubDict(typeName +
"Coeffs")
199 refineDict.
lookup(
"correctFluxes")
205 correctFluxes_.insert(fluxVelocities[i][0], fluxVelocities[i][1]);
223 meshCutter_.setRefinement(cellsToRefine, meshMod);
229 Info<<
"Refined from " 231 <<
" to " << globalData().nTotalCells() <<
" cells." <<
endl;
236 for (
label facei = 0; facei < nInternalFaces(); facei++)
238 label oldFacei = map().faceMap()[facei];
240 if (oldFacei >= nInternalFaces())
243 <<
"New internal face:" << facei
244 <<
" fc:" << faceCentres()[facei]
245 <<
" originates from boundary oldFace:" << oldFacei
278 const labelList& reverseFaceMap = map().reverseFaceMap();
287 label oldFacei = faceMap[facei];
291 label masterFacei = reverseFaceMap[oldFacei];
296 <<
"Problem: should not have removed faces" 300 else if (masterFacei != facei)
302 masterFaces.insert(masterFacei);
308 Pout<<
"Found " << masterFaces.size() <<
" split faces " <<
endl;
313 lookupClass<surfaceScalarField>()
317 if (!correctFluxes_.found(iter.key()))
320 <<
"Cannot find surfaceScalarField " << iter.key()
321 <<
" in user-provided flux mapping table " 322 << correctFluxes_ <<
endl 323 <<
" The flux mapping table is used to recreate the" 324 <<
" flux on newly created faces." <<
endl 325 <<
" Either add the entry if it is a flux or use (" 326 << iter.key() <<
" none) to suppress this warning." 331 const word& UName = correctFluxes_[iter.key()];
340 Pout<<
"Setting surfaceScalarField " << iter.key()
341 <<
" to NaN" <<
endl;
352 Pout<<
"Mapping flux " << iter.key()
353 <<
" using interpolated flux " << UName
362 lookupObject<volVectorField>(UName)
368 for (
label facei = 0; facei < nInternalFaces(); facei++)
370 label oldFacei = faceMap[facei];
375 phi[facei] = phiU[facei];
377 else if (reverseFaceMap[oldFacei] != facei)
380 phi[facei] = phiU[facei];
385 surfaceScalarField::Boundary& phiBf =
397 label oldFacei = faceMap[facei];
402 patchPhi[i] = patchPhiU[i];
404 else if (reverseFaceMap[oldFacei] != facei)
407 patchPhi[i] = patchPhiU[i];
417 label facei = iter.key();
419 if (isInternalFace(facei))
421 phi[facei] = phiU[facei];
433 patchPhi[i] = patchPhiU[i];
442 meshCutter_.updateMesh(map);
445 if (protectedCell_.size())
449 forAll(newProtectedCell, celli)
451 label oldCelli = map().cellMap()[celli];
452 newProtectedCell.
set(celli, protectedCell_.get(oldCelli));
454 protectedCell_.transfer(newProtectedCell);
458 meshCutter_.checkRefinementLevels(-1,
labelList(0));
473 meshCutter_.setUnrefinement(splitPoints, meshMod);
487 label pointi = splitPoints[i];
489 const labelList& pEdges = pointEdges()[pointi];
493 label otherPointi = edges()[pEdges[j]].otherVertex(pointi);
499 faceToSplitPoint.insert(pFaces[pFacei], otherPointi);
510 Info<<
"Unrefined from " 512 <<
" to " << globalData().nTotalCells() <<
" cells." 535 const labelList& reversePointMap = map().reversePointMap();
536 const labelList& reverseFaceMap = map().reverseFaceMap();
540 lookupClass<surfaceScalarField>()
544 if (!correctFluxes_.found(iter.key()))
547 <<
"Cannot find surfaceScalarField " << iter.key()
548 <<
" in user-provided flux mapping table " 549 << correctFluxes_ <<
endl 550 <<
" The flux mapping table is used to recreate the" 551 <<
" flux on newly created faces." <<
endl 552 <<
" Either add the entry if it is a flux or use (" 553 << iter.key() <<
" none) to suppress this warning." 558 const word& UName = correctFluxes_[iter.key()];
567 Info<<
"Mapping flux " << iter.key()
568 <<
" using interpolated flux " << UName
573 surfaceScalarField::Boundary& phiBf =
580 lookupObject<volVectorField>(UName)
588 label oldFacei = iter.key();
589 label oldPointi = iter();
591 if (reversePointMap[oldPointi] < 0)
594 label facei = reverseFaceMap[oldFacei];
598 if (isInternalFace(facei))
600 phi[facei] = phiU[facei];
610 patchPhi[i] = patchPhiU[i];
620 meshCutter_.updateMesh(map);
623 if (protectedCell_.size())
627 forAll(newProtectedCell, celli)
629 label oldCelli = map().cellMap()[celli];
632 newProtectedCell.
set(celli, protectedCell_.get(oldCelli));
635 protectedCell_.transfer(newProtectedCell);
639 meshCutter_.checkRefinementLevels(-1,
labelList(0));
650 forAll(pointCells(), pointi)
652 const labelList& pCells = pointCells()[pointi];
656 vFld[pCells[i]] =
max(vFld[pCells[i]], pFld[pointi]);
668 forAll(pointCells(), pointi)
670 const labelList& pCells = pointCells()[pointi];
674 pFld[pointi] =
max(pFld[pointi], vFld[pCells[i]]);
686 forAll(pointCells(), pointi)
688 const labelList& pCells = pointCells()[pointi];
693 sum += vFld[pCells[i]];
695 pFld[pointi] = sum/pCells.
size();
704 const scalar minLevel,
705 const scalar maxLevel
712 scalar err =
min(fld[i]-minLevel, maxLevel-fld[i]);
725 const scalar lowerRefineLevel,
726 const scalar upperRefineLevel,
749 if (cellError[celli] > 0)
751 candidateCell.
set(celli, 1);
759 const label maxCells,
760 const label maxRefinement,
765 label nTotToRefine = (maxCells - globalData().nTotalCells()) / 7;
767 const labelList& cellLevel = meshCutter_.cellLevel();
772 calculateProtectedCells(unrefineableCell);
775 label nLocalCandidates = count(candidateCell, 1);
781 if (nCandidates < nTotToRefine)
783 forAll(candidateCell, celli)
787 cellLevel[celli] < maxRefinement
788 && candidateCell.
get(celli)
790 unrefineableCell.
empty()
791 || !unrefineableCell.
get(celli)
802 for (
label level = 0; level < maxRefinement; level++)
804 forAll(candidateCell, celli)
808 cellLevel[celli] == level
809 && candidateCell.
get(celli)
811 unrefineableCell.
empty()
812 || !unrefineableCell.
get(celli)
830 meshCutter_.consistentRefinement
838 <<
" cells for refinement out of " << globalData().nTotalCells()
841 return consistentSet;
847 const scalar unrefineLevel,
853 const labelList splitPoints(meshCutter_.getSplitPoints());
859 label pointi = splitPoints[i];
861 if (pFld[pointi] < unrefineLevel)
864 const labelList& pCells = pointCells()[pointi];
866 bool hasMarked =
false;
870 if (markedCell.
get(pCells[pCelli]))
879 newSplitPoints.append(pointi);
885 newSplitPoints.shrink();
890 meshCutter_.consistentUnrefinement
897 <<
" split points out of a possible " 901 return consistentSet;
911 boolList markedFace(nFaces(),
false);
915 if (markedCell.
get(celli))
921 markedFace[cFaces[i]] =
true;
929 for (
label facei = 0; facei < nInternalFaces(); facei++)
931 if (markedFace[facei])
933 markedCell.
set(faceOwner()[facei], 1);
934 markedCell.
set(faceNeighbour()[facei], 1);
937 for (
label facei = nInternalFaces(); facei < nFaces(); facei++)
939 if (markedFace[facei])
941 markedCell.
set(faceOwner()[facei], 1);
953 const labelList& cellLevel = meshCutter_.cellLevel();
954 const labelList& pointLevel = meshCutter_.pointLevel();
958 forAll(pointLevel, pointi)
960 const labelList& pCells = pointCells(pointi);
964 label celli = pCells[pCelli];
966 if (pointLevel[pointi] <= cellLevel[celli])
969 if (nAnchorPoints[celli] == 8)
971 if (protectedCell.
set(celli,
true))
977 if (!protectedCell[celli])
979 nAnchorPoints[celli]++;
986 forAll(protectedCell, celli)
988 if (!protectedCell[celli] && nAnchorPoints[celli] != 8)
990 protectedCell.
set(celli,
true);
999 Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(
const IOobject& io)
1004 nRefinementIterations_(0),
1005 protectedCell_(nCells(), 0)
1024 label nProtected = 0;
1032 label celli = pCells[i];
1036 if (pointLevel[pointi] <= cellLevel[celli])
1040 if (nAnchors[celli] > 8)
1065 neiLevel[facei] = cellLevel[
faceOwner()[facei]];
1086 if (pointLevel[f[fp]] <= faceLevel)
1092 protectedFace[facei] =
true;
1103 if (protectedFace[facei])
1113 if (protectedFace[facei])
1125 if (cFaces.
size() < 6)
1159 cellSet protectedCells(*
this,
"protectedCells", nProtected);
1164 protectedCells.
insert(celli);
1169 <<
" cells that are protected from refinement." 1170 <<
" Writing these to cellSet " 1171 << protectedCells.
name()
1174 protectedCells.
write();
1210 bool hasChanged =
false;
1212 if (refineInterval == 0)
1218 else if (refineInterval < 0)
1221 <<
"Illegal refineInterval " << refineInterval <<
nl 1222 <<
"The refineInterval setting in the dynamicMeshDict should" 1223 <<
" be >= 1." <<
nl 1240 <<
"Illegal maximum number of cells " << maxCells <<
nl 1241 <<
"The maxCells setting in the dynamicMeshDict should" 1248 if (maxRefinement <= 0)
1251 <<
"Illegal maximum refinement level " << maxRefinement <<
nl 1252 <<
"The maxCells setting in the dynamicMeshDict should" 1257 const word fieldName(refineDict.
lookup(
"field"));
1259 const volScalarField& vFld = lookupObject<volScalarField>(fieldName);
1261 const scalar lowerRefineLevel =
1263 const scalar upperRefineLevel =
1270 const label nBufferLayers =
1304 if (nCellsToRefine > 0)
1312 const labelList& cellMap = map().cellMap();
1313 const labelList& reverseCellMap = map().reverseCellMap();
1319 label oldCelli = cellMap[celli];
1323 newRefineCell.set(celli, 1);
1325 else if (reverseCellMap[oldCelli] != celli)
1327 newRefineCell.set(celli, 1);
1331 newRefineCell.set(celli, refineCell.get(oldCelli));
1334 refineCell.
transfer(newRefineCell);
1339 for (
label i = 0; i < nBufferLayers; i++)
1363 pointsToUnrefine.
size(),
1367 if (nSplitPoints > 0)
1436 scalarCellLevel[celli] = cellLevel[celli];
1439 writeOk = writeOk && scalarCellLevel.
write();
const refinementHistory & history() const
#define forAll(list, i)
Loop across all elements in list.
virtual labelList selectRefineCells(const label maxCells, const label maxRefinement, const PackedBoolList &candidateCell) const
Subset candidate cells for refinement.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
const word & name() const
Return name.
scalarField maxPointField(const scalarField &) const
Get per cell max of connected point.
static label count(const PackedBoolList &, const unsigned int)
Count set/unset elements in packedlist.
bool moving() const
Is mesh moving.
const hexRef8 & meshCutter() const
Direct access to the refinement engine.
errorManipArg< error, int > exit(error &err, const int errNo=1)
PackedBoolList protectedCell_
Protected cells (usually since not hexes)
unsigned int get(const label) const
Get value at index I.
A face is a list of labels corresponding to mesh vertices.
const labelIOList & pointLevel() const
A list of keyword definitions, which are a keyword followed by any number of values (e...
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
label nInternalFaces() const
virtual const labelList & faceNeighbour() const
Return face neighbour.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
scalarField cellToPoint(const scalarField &vFld) const
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
void size(const label)
Override size to be inconsistent with allocated storage.
virtual bool writeObject(IOstream::streamFormat fmt, IOstream::versionNumber ver, IOstream::compressionType cmp, const bool valid) const
Write using given format, version and compression.
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface...
Ostream & endl(Ostream &os)
Add newline and flush stream.
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none.
void checkEightAnchorPoints(PackedBoolList &protectedCell, label &nProtected) const
Check all cells have 8 anchor points.
const cellList & cells() const
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
const Time & time() const
Return the top-level database.
void resize(const label)
Alias for setSize(const label)
bool insert(const Key &key)
Insert a new entry.
void calculateProtectedCells(PackedBoolList &unrefineableCell) const
Calculate cells that cannot be refined since would trigger.
label size() const
Return number of elements in table.
Macros for easy insertion into run-time selection tables.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Class to handle errors and exceptions in a simple, consistent stream-based manner.
virtual bool update()
Update the mesh for both mesh motion and topology change.
void set(const PackedList< 1 > &)
Set specified bits.
autoPtr< mapPolyMesh > unrefine(const labelList &)
Unrefine cells. Gets passed in centre points of cells to combine.
const dictionary & optionalSubDict(const word &) const
Find and return a sub-dictionary if found.
virtual bool writeObject(IOstream::streamFormat fmt, IOstream::versionNumber ver, IOstream::compressionType cmp, const bool valid) const
Write the underlying polyMesh and other data.
Refinement of (split) hexes using polyTopoChange.
A class for handling words, derived from string.
hexRef8 meshCutter_
Mesh cutting engine.
void extendMarkedCells(PackedBoolList &markedCell) const
Extend markedCell with cell-face-cell.
virtual const labelList & faceOwner() const
Return face owner.
streamFormat
Enumeration for the format of data in the stream.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
autoPtr< mapPolyMesh > refine(const labelList &)
Refine cells. Update mesh and fields.
const globalMeshData & globalData() const
Return parallel info.
List< label > labelList
A List of labels.
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if successful.
virtual const faceList & faces() const
Return raw faces.
static void fillNan(UList< scalar > &)
Fill block of data with NaN.
Container with cells to refine. Refinement given as single direction.
An STL-conforming hash table.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
errorManip< error > abort(error &err)
void clear()
Clear the list, i.e. set addressable size to zero.
label readLabel(Istream &is)
const labelListList & pointCells() const
compressionType
Enumeration for the format of data in the stream.
Internal::FieldType & primitiveFieldRef()
Return a reference to the internal field.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
A topoSetSource to select points based on usage in cells.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
label nRefinementIterations_
Number of refinement/unrefinement steps done so far.
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.
void setInstance(const fileName &)
Set the instance for mesh files.
virtual void selectRefineCandidates(const scalar lowerRefineLevel, const scalar upperRefineLevel, const scalarField &vFld, PackedBoolList &candidateCell) const
Select candidate cells for refinement.
const fvPatch & patch() const
Return patch.
void readDict()
Read the projection parameters from dictionary.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
#define WarningInFunction
Report a warning using Foam::Warning.
virtual labelList selectUnrefinePoints(const scalar unrefineLevel, const PackedBoolList &markedCell, const scalarField &pFld) const
Select points that can be unrefined.
bool write(const bool valid=true) const
Force writing refinement+history to polyMesh directory.
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
label start() const
Return start label of this patch in the polyMesh face list.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Abstract base class for geometry and/or topology changing fvMesh.
A cell is defined as a list of faces with extra functionality.
A collection of cell labels.
prefixOSstream Pout(cout, "Pout")
Direct mesh changes based on v1.3 polyTopoChange syntax.
const dimensionedScalar c
Speed of light in a vacuum.
bool topoChanging() const
Is mesh topology changing.
virtual ~dynamicRefineFvMesh()
Destructor.
bool empty() const
Return true if the list is empty (ie, size() is zero).
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
autoPtr< mapPolyMesh > changeMesh(polyMesh &mesh, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
All refinement history. Used in unrefinement.
virtual bool write(const bool valid=true) const
Write using setting from DB.
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]
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
An abstract base class with a fat-interface to all derived classes covering all possible ways in whic...
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
scalarField error(const scalarField &fld, const scalar minLevel, const scalar maxLevel) const
const labelIOList & cellLevel() const
Switch dumpLevel_
Dump cellLevel for postprocessing.
scalarField maxCellField(const volScalarField &) const
Get point max of connected cell.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.