86 mesh.time().timeName(),
101 fld[celli] = elems[celli];
115 label
diff = neighbour[facei] - owner[facei];
129 const bool calculateIntersect,
135 scalar& sumSqrIntersect
143 label own = owner[facei];
144 label nei = neighbour[facei];
147 label
diff = nei-own;
148 cellBandwidth[nei] =
max(cellBandwidth[nei], diff);
151 bandwidth =
max(cellBandwidth);
155 forAll(cellBandwidth, celli)
157 profile += 1.0*cellBandwidth[celli];
160 sumSqrIntersect = 0.0;
161 if (calculateIntersect)
165 for (label colI = celli-cellBandwidth[celli]; colI <= celli; colI++)
167 nIntersect[colI] += 1.0;
192 forAll(cellOrder, newCelli)
194 label oldCelli = cellOrder[newCelli];
196 const cell& cFaces =
mesh.cells()[oldCelli];
203 label facei = cFaces[i];
205 if (
mesh.isInternalFace(facei))
208 label nbrCelli = reverseCellOrder[
mesh.faceNeighbour()[facei]];
209 if (nbrCelli == newCelli)
211 nbrCelli = reverseCellOrder[
mesh.faceOwner()[facei]];
214 if (newCelli < nbrCelli)
234 for (
const label index : order)
236 if (nbr[index] != -1)
238 oldToNewFace[cFaces[index]] = newFacei++;
244 for (label facei = newFacei; facei <
mesh.nFaces(); facei++)
246 oldToNewFace[facei] = facei;
251 forAll(oldToNewFace, facei)
253 if (oldToNewFace[facei] == -1)
256 <<
"Did not determine new position" <<
" for face " << facei
257 <<
abort(FatalError);
280 label prevRegion = -1;
282 forAll(cellOrder, newCelli)
284 label oldCelli = cellOrder[newCelli];
286 if (cellToRegion[oldCelli] != prevRegion)
288 prevRegion = cellToRegion[oldCelli];
291 const cell& cFaces =
mesh.cells()[oldCelli];
297 label facei = cFaces[i];
299 if (
mesh.isInternalFace(facei))
302 label nbrCelli = reverseCellOrder[
mesh.faceNeighbour()[facei]];
303 if (nbrCelli == newCelli)
305 nbrCelli = reverseCellOrder[
mesh.faceOwner()[facei]];
308 if (cellToRegion[oldCelli] != cellToRegion[cellOrder[nbrCelli]])
313 else if (newCelli < nbrCelli)
337 oldToNewFace[cFaces[nbr.indices()[i]]] = newFacei++;
343 label nRegions =
max(cellToRegion)+1;
348 for (label facei = 0; facei <
mesh.nInternalFaces(); facei++)
350 label ownRegion = cellToRegion[
mesh.faceOwner()[facei]];
351 label neiRegion = cellToRegion[
mesh.faceNeighbour()[facei]];
353 if (ownRegion != neiRegion)
356 min(ownRegion, neiRegion)*nRegions
357 +
max(ownRegion, neiRegion);
366 label
key = sortKey[i];
378 oldToNewFace[sortKey.indices()[i]] = newFacei++;
383 for (label facei = newFacei; facei <
mesh.nFaces(); facei++)
385 oldToNewFace[facei] = facei;
390 forAll(oldToNewFace, facei)
392 if (oldToNewFace[facei] == -1)
395 <<
"Did not determine new position"
396 <<
" for face " << facei
397 <<
abort(FatalError);
438 forAll(newNeighbour, facei)
440 label own = newOwner[facei];
441 label nei = newNeighbour[facei];
445 newFaces[facei].flip();
446 std::swap(newOwner[facei], newNeighbour[facei]);
447 flipFaceFlux.insert(facei);
459 patchSizes[patchi] =
patches[patchi].size();
460 patchStarts[patchi] =
patches[patchi].start();
461 oldPatchNMeshPoints[patchi] =
patches[patchi].nPoints();
488 label oldFacei = fZone[i];
489 newAddressing[i] = reverseFaceOrder[oldFacei];
490 if (flipFaceFlux.found(newAddressing[i]))
492 newFlipMap[i] = !fZone.
flipMap()[i];
496 newFlipMap[i] = fZone.
flipMap()[i];
565 Info<<
"Determining cell order:" <<
endl;
569 label nRegions =
max(cellToRegion)+1;
575 forAll(regionToCells, regioni)
577 Info<<
" region " << regioni <<
" starts at " << celli <<
endl;
580 const bool oldParRun = UPstream::parRun(
false);
585 const fvMesh& subMesh = subsetter.subMesh();
593 UPstream::parRun(oldParRun);
596 const labelList& cellMap = subsetter.cellMap();
600 cellOrder[celli++] = cellMap[subCellOrder[i]];
611int main(
int argc,
char *argv[])
615 "Renumber mesh cells to reduce the bandwidth"
622 argList::addOption(
"dict",
"file",
"Alternative renumberMeshDict");
624 argList::addBoolOption
627 "Calculate the rms of the front-width"
630 argList::noFunctionObjects();
639 Info<<
"renumberMesh built with zoltan support." <<
nl <<
endl;
640 (void)zoltanRenumber::typeName;
653 const bool readDict =
args.
found(
"dict");
654 const bool doFrontWidth =
args.
found(
"frontWidth");
655 const bool overwrite =
args.
found(
"overwrite");
663 const word oldInstance =
mesh.pointsInstance();
667 scalar sumSqrIntersect;
673 mesh.faceNeighbour(),
687 )/
mesh.globalData().nTotalCells()
691 <<
" size: " <<
mesh.globalData().nTotalCells() <<
nl
692 <<
"Before renumbering :" <<
nl
693 <<
" band : " << band <<
nl
694 <<
" profile : " << profile <<
nl;
698 Info<<
" rms frontwidth : " << rmsFrontwidth <<
nl;
703 bool sortCoupledFaceCells =
false;
704 bool writeMaps =
false;
705 bool orderPoints =
false;
722 renumberPtr = renumberMethod::New(renumberDict);
726 "sortCoupledFaceCells",
729 if (sortCoupledFaceCells)
731 Info<<
"Sorting cells on coupled boundaries to be last." <<
nl
738 Info<<
"Ordering cells into regions of size " << blockSize
739 <<
" (using decomposition);"
740 <<
" ordering faces into region-internal"
741 <<
" and region-external."
744 if (blockSize < 0 || blockSize >=
mesh.nCells())
747 <<
"Block size " << blockSize
748 <<
" should be positive integer"
749 <<
" and less than the number of cells in the mesh."
754 orderPoints = renumberDict.
getOrDefault(
"orderPoints",
false);
757 Info<<
"Ordering points into internal and boundary points."
762 renumberDict.
readEntry(
"writeMaps", writeMaps);
765 Info<<
"Writing renumber maps (new to old) to polyMesh." <<
nl
771 Info<<
"Using default renumberMethod." <<
nl <<
endl;
776 Info<<
"Selecting renumberMethod " << renumberPtr().type() <<
nl
786 "cellProcAddressing",
787 mesh.facesInstance(),
788 polyMesh::meshSubDir,
790 IOobject::READ_IF_PRESENT,
800 "faceProcAddressing",
801 mesh.facesInstance(),
802 polyMesh::meshSubDir,
804 IOobject::READ_IF_PRESENT,
813 "pointProcAddressing",
814 mesh.pointsInstance(),
815 polyMesh::meshSubDir,
817 IOobject::READ_IF_PRESENT,
826 "boundaryProcAddressing",
827 mesh.pointsInstance(),
828 polyMesh::meshSubDir,
830 IOobject::READ_IF_PRESENT,
920 label nBlocks =
mesh.nCells()/blockSize;
921 Info<<
"nBlocks = " << nBlocks <<
endl;
924 dictionary decomposeDict(renumberDictPtr().subDict(
"blockCoeffs"));
925 decomposeDict.set(
"numberOfSubdomains", nBlocks);
927 const bool oldParRun = UPstream::parRun(
false);
936 decomposePtr().decompose
943 UPstream::parRun(oldParRun);
954 Info<<
nl <<
"Written decomposition as volScalarField to "
955 <<
"cellDist for use in postprocessing."
959 cellOrder = regionRenumber(renumberPtr(),
mesh, cellToRegion);
962 faceOrder = getRegionFaceOrder
972 cellOrder = renumberPtr().renumber
978 if (sortCoupledFaceCells)
989 nBndCells += pbm[patchi].
size();
1007 if (reverseCellOrder[celli] != -1)
1009 bndCells[nBndCells] = celli;
1010 bndCellMap[nBndCells++] =
1011 reverseCellOrder[celli];
1012 reverseCellOrder[celli] = -1;
1018 bndCellMap.setSize(nBndCells);
1026 label sortedI =
mesh.nCells();
1029 label origCelli = bndCells[order[i]];
1030 newReverseCellOrder[origCelli] = --sortedI;
1033 Info<<
"Ordered all " << nBndCells
1034 <<
" cells with a coupled face"
1035 <<
" to the end of the cell list, starting at " << sortedI
1040 forAll(cellOrder, newCelli)
1042 label origCelli = cellOrder[newCelli];
1043 if (newReverseCellOrder[origCelli] == -1)
1045 newReverseCellOrder[origCelli] = sortedI++;
1050 cellOrder =
invert(
mesh.nCells(), newReverseCellOrder);
1055 faceOrder = getFaceOrder
1089 pointOrderMap().pointMap()
1094 pointOrderMap().reversePointMap(),
1095 const_cast<labelList&
>(map().reversePointMap())
1101 mesh.updateMesh(map());
1104 if (cellProcAddressing.headerOk())
1106 bool localOk = (cellProcAddressing.size() ==
mesh.nCells());
1110 Info<<
"Renumbering processor cell decomposition map "
1111 << cellProcAddressing.name() <<
endl;
1120 Info<<
"Not writing inconsistent processor cell decomposition"
1121 <<
" map " << cellProcAddressing.filePath() <<
endl;
1122 cellProcAddressing.writeOpt(IOobject::NO_WRITE);
1127 cellProcAddressing.writeOpt(IOobject::NO_WRITE);
1130 if (faceProcAddressing.headerOk())
1132 bool localOk = (faceProcAddressing.size() ==
mesh.nFaces());
1136 Info<<
"Renumbering processor face decomposition map "
1137 << faceProcAddressing.name() <<
endl;
1146 for (
const label facei : fff)
1148 label masterFacei = faceProcAddressing[facei];
1150 faceProcAddressing[facei] = -masterFacei;
1152 if (masterFacei == 0)
1155 <<
" masterFacei:" << masterFacei
1156 <<
exit(FatalError);
1162 Info<<
"Not writing inconsistent processor face decomposition"
1163 <<
" map " << faceProcAddressing.filePath() <<
endl;
1164 faceProcAddressing.writeOpt(IOobject::NO_WRITE);
1169 faceProcAddressing.writeOpt(IOobject::NO_WRITE);
1172 if (pointProcAddressing.headerOk())
1174 bool localOk = (pointProcAddressing.size() ==
mesh.nPoints());
1178 Info<<
"Renumbering processor point decomposition map "
1179 << pointProcAddressing.name() <<
endl;
1188 Info<<
"Not writing inconsistent processor point decomposition"
1189 <<
" map " << pointProcAddressing.filePath() <<
endl;
1190 pointProcAddressing.writeOpt(IOobject::NO_WRITE);
1195 pointProcAddressing.writeOpt(IOobject::NO_WRITE);
1198 if (boundaryProcAddressing.headerOk())
1202 boundaryProcAddressing.size()
1203 ==
mesh.boundaryMesh().size()
1211 Info<<
"Not writing inconsistent processor patch decomposition"
1212 <<
" map " << boundaryProcAddressing.filePath() <<
endl;
1213 boundaryProcAddressing.writeOpt(IOobject::NO_WRITE);
1218 boundaryProcAddressing.writeOpt(IOobject::NO_WRITE);
1225 if (map().hasMotionPoints())
1227 mesh.movePoints(map().preMotionPoints());
1234 scalar sumSqrIntersect;
1240 mesh.faceNeighbour(),
1253 )/
mesh.globalData().nTotalCells()
1256 Info<<
"After renumbering :" <<
nl
1257 <<
" band : " << band <<
nl
1258 <<
" profile : " << profile <<
nl;
1263 Info<<
" rms frontwidth : " << rmsFrontwidth <<
nl;
1282 mesh.nInternalPoints(),
1293 mesh.nInternalEdges(),
1298 mesh.nInternal0Edges(),
1303 mesh.nInternal1Edges(),
1308 <<
" total : " << nTotPoints <<
nl
1309 <<
" internal: " << nTotIntPoints <<
nl
1310 <<
" boundary: " << nTotPoints-nTotIntPoints <<
nl
1312 <<
" total : " << nTotEdges <<
nl
1313 <<
" internal: " << nTotIntEdges <<
nl
1314 <<
" internal using 0 boundary points: "
1315 << nTotInt0Edges <<
nl
1316 <<
" internal using 1 boundary points: "
1317 << nTotInt1Edges-nTotInt0Edges <<
nl
1318 <<
" internal using 2 boundary points: "
1319 << nTotIntEdges-nTotInt1Edges <<
nl
1320 <<
" boundary: " << nTotEdges-nTotIntEdges <<
nl
1327 mesh.setInstance(oldInstance);
1335 Info<<
"Writing mesh to " <<
mesh.facesInstance() <<
endl;
1338 processorMeshes::removeFiles(
mesh);
1346 mesh.facesInstance(),
1347 polyMesh::meshSubDir,
1349 IOobject::READ_IF_PRESENT,
1354 refData.updateMesh(map());
1358 topoSet::updateMesh(
mesh.facesInstance(), map(), cellSets);
1359 topoSet::updateMesh(
mesh.facesInstance(), map(), faceSets);
1360 topoSet::updateMesh(
mesh.facesInstance(), map(), pointSets);
1382 <<
"Written current cellID and origCellID as volScalarField"
1383 <<
" for use in postprocessing." <<
nl <<
endl;
1390 mesh.facesInstance(),
1391 polyMesh::meshSubDir,
1405 mesh.facesInstance(),
1406 polyMesh::meshSubDir,
1420 mesh.facesInstance(),
1421 polyMesh::meshSubDir,
Field reading functions for post-processing utilities.
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
reduce(hasMovingMesh, orOp< bool >())
Cuthill-McKee renumbering.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
List of IOobjects with searching and retrieving facilities.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
void setSize(const label n)
Alias for resize()
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
A list that is sorted upon construction or when explicitly requested with the sort() method.
A List with indirect addressing. Like IndirectList but does not store addressing.
void size(const label n)
Older name for setAddressableSize.
label size() const noexcept
The number of elements in the list.
void clearAddressing()
Clear addressing.
bool found(const word &optName) const
Return true if the named option is found.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
void reset(autoPtr< T > &&other) noexcept
Delete managed object and set to new given pointer.
A cell is defined as a list of faces with extra functionality.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, bool mandatory=true) const
Smooth ATC in cells next to a set of patches supplied by type.
A subset of mesh faces organised as a primitive patch.
virtual void resetAddressing(const labelUList &addr, const bool flipMapValue)
Reset addressing - use uniform flip map value.
const boolList & flipMap() const noexcept
Return face flip map.
Holds a reference to the original mesh (the baseMesh) and optionally to a subset of that mesh (the su...
Mesh data needed to do the Finite Volume discretisation.
Various for reading/decomposing/reconstructing/distributing refinement data.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
UPtrList< const labelUList > faceCells() const
Return a list of faceCells for each patch.
Mesh consisting of general polyhedral cells.
Direct mesh changes based on v1.3 polyTopoChange syntax.
Cell-face mesh analysis engine.
const vectorField & cellCentres() const
virtual bool write(const bool valid=true) const
Write using setting from DB.
Abstract base class for renumbering.
virtual labelList renumber(const pointField &) const
A class for managing temporary objects.
A class for handling words, derived from Foam::string.
const polyBoundaryMesh & patches
bool coupled(solutionDict.getOrDefault("coupledEnergyField", false))
Foam::PtrList< Foam::fvMesh > meshes(regionNames.size())
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const word dictName("faMeshDefinition")
void inplaceRenumber(const labelUList &oldToNew, IntListType &lists)
Inplace renumber the values (not the indices) of a list of lists.
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
wordList ReadFields(const typename GeoMesh::Mesh &mesh, const IOobjectList &objects, PtrList< GeometricField< Type, PatchField, GeoMesh > > &fields, const bool syncPar=true, const bool readOldTime=false)
Read Geometric fields of templated type.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar sqrt(const dimensionedScalar &ds)
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
void sort(UList< T > &list)
Sort the list.
errorManip< error > abort(error &err)
labelList invert(const label len, const labelUList &map)
Create an inverse one-to-one mapping.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
labelListList invertOneToMany(const label len, const labelUList &map)
Invert one-to-many map. Unmapped elements will be size 0.
errorManipArg< error, int > exit(error &err, const int errNo=1)
ListType reorder(const labelUList &oldToNew, const ListType &input, const bool prune=false)
Reorder the elements of a list.
constexpr char nl
The newline '\n' character (0x0a)
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
#define forAllReverse(list, i)
Reverse loop across all elements in list.
static const char *const typeName
The type name used in ensight case files.