53 for (
auto idx : locations)
73 if (selectedElements.
empty() || subsetMap.
empty())
80 const bitSet selected(nElems, selectedElements);
86 if (selected[subsetMap[i]])
98 if (selected[subsetMap[i]])
100 subsettedElements[
n] = i;
105 return subsettedElements;
118 <<
"Mesh is not subsetted!" <<
nl 128 void Foam::fvMeshSubset::calcFaceFlipMap()
const 131 const labelList& subToBaseCell = cellMap();
133 faceFlipMapPtr_.reset(
new labelList(subToBaseFace.size()));
134 auto& faceFlipMap = *faceFlipMapPtr_;
138 const label subInt = subMesh().nInternalFaces();
140 const labelList& subOwn = subMesh().faceOwner();
141 const labelList& own = baseMesh_.faceOwner();
143 for (label subFaceI = 0; subFaceI < subInt; ++subFaceI)
145 faceFlipMap[subFaceI] = subToBaseFace[subFaceI]+1;
147 for (label subFaceI = subInt; subFaceI < subOwn.size(); ++subFaceI)
149 const label faceI = subToBaseFace[subFaceI];
150 if (subToBaseCell[subOwn[subFaceI]] == own[faceI])
152 faceFlipMap[subFaceI] = faceI+1;
156 faceFlipMap[subFaceI] = -faceI-1;
162 void Foam::fvMeshSubset::doCoupledPatches
171 const polyBoundaryMesh& oldPatches = baseMesh().
boundaryMesh();
173 label nUncoupled = 0;
177 PstreamBuffers pBufs;
180 if (!nCellsUsingFace.empty())
182 for (
const polyPatch&
pp : oldPatches)
184 const auto* procPatch = isA<processorPolyPatch>(
pp);
188 const label nbrProci = procPatch->neighbProcNo();
190 UOPstream toNbr(nbrProci, pBufs);
192 SubList<label>(nCellsUsingFace,
pp.size(),
pp.start());
197 pBufs.finishedSends();
200 for (
const polyPatch&
pp : oldPatches)
202 const auto* procPatch = isA<processorPolyPatch>(
pp);
206 const label nbrProci = procPatch->neighbProcNo();
208 if (!pBufs.recvDataCount(nbrProci))
216 UIPstream fromNbr(nbrProci, pBufs);
217 fromNbr >> nbrCellsUsingFace;
220 if (!nCellsUsingFace.empty() && !nbrCellsUsingFace.empty())
228 nCellsUsingFace[
pp.start()+i] == 1
229 && nbrCellsUsingFace[i] == 0
234 nCellsUsingFace[
pp.start()+i] = 3;
244 for (
const polyPatch&
pp : oldPatches)
246 const cyclicPolyPatch* cpp = isA<cyclicPolyPatch>(
pp);
248 if (cpp && !nCellsUsingFace.empty())
250 const auto& cycPatch = *cpp;
254 label thisFacei = cycPatch.start() + i;
255 label otherFacei = cycPatch.transformGlobalFace(thisFacei);
259 nCellsUsingFace[thisFacei] == 1
260 && nCellsUsingFace[otherFacei] == 0
263 nCellsUsingFace[thisFacei] = 3;
272 reduce(nUncoupled, sumOp<label>());
277 Info<<
"Uncoupled " << nUncoupled <<
" faces on coupled patches. " 278 <<
"(processorPolyPatch, cyclicPolyPatch)" <<
endl;
283 void Foam::fvMeshSubset::subsetZones()
291 auto& newSubMesh = subMeshPtr_();
297 List<pointZone*>
pZones(pointZones.size());
301 const pointZone& pz = pointZones[zonei];
303 pZones[zonei] =
new pointZone
308 newSubMesh.pointZones()
319 List<faceZone*> fZones(faceZones.size());
323 const faceZone& fz = faceZones[zonei];
332 zone[fz[j]] = (fz.flipMap()[j] ? 1 : -1);
349 const label meshFacei =
faceMap()[subFacei];
350 if (zone[meshFacei] != 0)
352 subAddressing[nSub] = subFacei;
353 const label subOwner = subMesh().faceOwner()[subFacei];
354 const label baseOwner = baseMesh().
faceOwner()[meshFacei];
356 const bool sameOwner = (cellMap()[subOwner] == baseOwner);
357 const bool flip = (zone[meshFacei] == 1);
358 subFlipStatus[nSub] = (sameOwner == flip);
364 fZones[zonei] =
new faceZone
370 newSubMesh.faceZones()
379 List<cellZone*> cZones(cellZones.size());
383 const cellZone& cz = cellZones[zonei];
385 cZones[zonei] =
new cellZone
390 newSubMesh.cellZones()
395 newSubMesh.addZones(
pZones, fZones, cZones);
404 subMeshPtr_(nullptr),
405 faceFlipMapPtr_(nullptr),
424 const bitSet& selectedCells,
431 reset(selectedCells,
patchID, syncPar);
445 reset(selectedCells,
patchID, syncPar);
459 reset(selectedCells,
patchID, syncPar);
474 reset(regioni, regions,
patchID, syncPar);
482 subMeshPtr_.reset(
nullptr);
483 faceFlipMapPtr_.reset(
nullptr);
507 subMeshPtr_.reset(std::move(subMeshPtr));
508 faceFlipMapPtr_.reset(
nullptr);
510 pointMap_ = std::move(pointMap);
512 cellMap_ = std::move(cellMap);
513 patchMap_ = std::move(patchMap);
535 baseMesh_.time().timeName(),
545 auto& newSubMesh = subMeshPtr_();
550 const polyBoundaryMesh& oldBoundary = baseMesh_.boundaryMesh();
551 const polyBoundaryMesh& newBoundary = newSubMesh.boundaryMesh();
555 patchMap_ =
identity(newPatches.size());
557 forAll(newPatches, patchi)
562 oldBoundary[patchi].clone
573 newSubMesh.addFvPatches(newPatches,
false);
584 const bitSet& selectedCells,
604 if (wantedPatchID == -1)
608 wantedPatchID = oldPatches.
findPatchID(exposedPatchName);
610 else if (wantedPatchID < 0 || wantedPatchID >= oldPatches.
size())
613 <<
"Non-existing patch index " << wantedPatchID <<
endl 614 <<
"Should be between 0 and " << oldPatches.
size()-1
624 label len = cellMap_.
size();
628 i >= 0 && (cellMap_[i] >= oldCells.
size());
634 cellMap_.resize(len);
650 label nFacesInSet = 0;
652 forAll(oldFaces, oldFacei)
654 bool faceUsed =
false;
656 if (selectedCells.
test(oldOwner[oldFacei]))
658 ++nCellsUsingFace[oldFacei];
665 && selectedCells.
test(oldNeighbour[oldFacei])
668 ++nCellsUsingFace[oldFacei];
677 faceMap_.resize(nFacesInSet);
680 doCoupledPatches(syncPar, nCellsUsingFace);
684 label oldInternalPatchID = 0;
687 label nextPatchID = oldPatches.
size();
693 label nbSize = oldPatches.
size();
695 if (wantedPatchID == -1)
701 forAll(oldPatches, patchi)
703 if (isA<processorPolyPatch>(oldPatches[patchi]))
705 nextPatchID = patchi;
708 ++oldInternalPatchID;
714 for (label oldPatchi = 0; oldPatchi < nextPatchID; oldPatchi++)
716 globalPatchMap[oldPatchi] = oldPatchi;
720 label oldPatchi = nextPatchID;
721 oldPatchi < oldPatches.
size();
725 globalPatchMap[oldPatchi] = oldPatchi+1;
730 oldInternalPatchID = wantedPatchID;
731 nextPatchID = wantedPatchID+1;
747 for (label oldFacei = 0; oldFacei < oldNInternalFaces; ++oldFacei)
749 if (nCellsUsingFace[oldFacei] == 2)
751 globalFaceMap[oldFacei] = facei;
752 faceMap_[facei++] = oldFacei;
755 markUsed(oldFaces[oldFacei], globalPointMap);
760 const label nInternalFaces = facei;
766 oldPatchi < oldPatches.
size()
767 && oldPatchi < nextPatchID;
771 const polyPatch& oldPatch = oldPatches[oldPatchi];
773 label oldFacei = oldPatch.
start();
777 if (nCellsUsingFace[oldFacei] == 1)
782 globalFaceMap[oldFacei] = facei;
783 faceMap_[facei++] = oldFacei;
786 markUsed(oldFaces[oldFacei], globalPointMap);
789 ++boundaryPatchSizes[globalPatchMap[oldPatchi]];
796 for (label oldFacei = 0; oldFacei < oldNInternalFaces; ++oldFacei)
798 if (nCellsUsingFace[oldFacei] == 1)
800 globalFaceMap[oldFacei] = facei;
801 faceMap_[facei++] = oldFacei;
804 markUsed(oldFaces[oldFacei], globalPointMap);
807 ++boundaryPatchSizes[oldInternalPatchID];
814 label oldFacei = oldNInternalFaces;
815 oldFacei < oldFaces.
size();
819 if (nCellsUsingFace[oldFacei] == 3)
821 globalFaceMap[oldFacei] = facei;
822 faceMap_[facei++] = oldFacei;
825 markUsed(oldFaces[oldFacei], globalPointMap);
828 ++boundaryPatchSizes[oldInternalPatchID];
835 label oldPatchi = nextPatchID;
836 oldPatchi < oldPatches.
size();
840 const polyPatch& oldPatch = oldPatches[oldPatchi];
842 label oldFacei = oldPatch.
start();
846 if (nCellsUsingFace[oldFacei] == 1)
851 globalFaceMap[oldFacei] = facei;
852 faceMap_[facei++] = oldFacei;
855 markUsed(oldFaces[oldFacei], globalPointMap);
858 ++boundaryPatchSizes[globalPatchMap[oldPatchi]];
864 if (facei != nFacesInSet)
872 label nPointsInSet = 0;
874 forAll(globalPointMap, pointi)
876 if (globalPointMap[pointi] != -1)
881 pointMap_.setSize(nPointsInSet);
885 forAll(globalPointMap, pointi)
887 if (globalPointMap[pointi] != -1)
889 pointMap_[nPointsInSet] = pointi;
890 globalPointMap[pointi] = nPointsInSet;
916 auto iter = newFaces.
begin();
917 const auto& renumbering = globalPointMap;
920 for (label facei = 0; facei < nInternalFaces; ++facei)
922 face& newItem = *iter;
925 const face& oldItem = oldFaces[faceMap_[facei]];
928 newItem.
resize(oldItem.size());
932 newItem[i] = renumbering[oldItem[i]];
937 for (label facei = nInternalFaces; facei < faceMap_.size(); ++facei)
939 const label oldFacei = faceMap_[facei];
941 face& newItem = *iter;
948 && selectedCells.
test(oldNeighbour[oldFacei])
949 && !selectedCells.
test(oldOwner[oldFacei])
952 ? oldFaces[oldFacei].reverseFace()
957 newItem.resize(oldItem.size());
961 newItem[i] = renumbering[oldItem[i]];
973 auto iter = newCells.
begin();
974 const auto& renumbering = globalFaceMap;
976 for (
const label oldCelli : cellMap_)
978 cell& newItem = *iter;
981 const labelList& oldItem = oldCells[oldCelli];
984 newItem.
resize(oldItem.size());
988 newItem[i] = renumbering[oldItem[i]];
1006 baseMesh_.time().timeName(),
1012 std::move(newPoints),
1013 std::move(newFaces),
1014 std::move(newCells),
1021 patchMap_.resize(nbSize);
1022 label nNewPatches = 0;
1023 label patchStart = nInternalFaces;
1031 labelList globalPatchSizes(boundaryPatchSizes);
1032 globalPatchSizes.setSize(nextPatchID);
1050 bool samePatches =
true;
1056 samePatches =
false;
1074 label oldPatchi = 0;
1075 oldPatchi < oldPatches.
size()
1076 && oldPatchi < nextPatchID;
1080 const label newSize = boundaryPatchSizes[globalPatchMap[oldPatchi]];
1082 if (oldInternalPatchID != oldPatchi)
1087 for (label patchFacei = 0; patchFacei < newSize; patchFacei++)
1089 const label facei = patchStart+patchFacei;
1090 const label oldFacei = faceMap_[facei];
1091 map[patchFacei] = oldPatches[oldPatchi].whichFace(oldFacei);
1097 oldPatches[oldPatchi].clone
1099 subMeshPtr_().boundaryMesh(),
1112 oldPatches[oldPatchi].clone
1114 subMeshPtr_().boundaryMesh(),
1122 patchStart += newSize;
1123 patchMap_[nNewPatches] = oldPatchi;
1129 if (wantedPatchID == -1)
1131 label oldInternalSize = boundaryPatchSizes[oldInternalPatchID];
1135 reduce(oldInternalSize, sumOp<label>());
1139 if (oldInternalSize > 0)
1147 boundaryPatchSizes[oldInternalPatchID],
1150 subMeshPtr_().boundaryMesh(),
1151 emptyPolyPatch::typeName
1160 patchStart += boundaryPatchSizes[oldInternalPatchID];
1161 patchMap_[nNewPatches] = -1;
1170 label oldPatchi = nextPatchID;
1171 oldPatchi < oldPatches.
size();
1175 const label newSize = boundaryPatchSizes[globalPatchMap[oldPatchi]];
1177 if (oldInternalPatchID != oldPatchi)
1182 for (label patchFacei = 0; patchFacei < newSize; patchFacei++)
1184 const label facei = patchStart+patchFacei;
1185 const label oldFacei = faceMap_[facei];
1186 map[patchFacei] = oldPatches[oldPatchi].whichFace(oldFacei);
1192 oldPatches[oldPatchi].clone
1194 subMeshPtr_().boundaryMesh(),
1207 oldPatches[oldPatchi].clone
1209 subMeshPtr_().boundaryMesh(),
1220 patchStart += newSize;
1221 patchMap_[nNewPatches] = oldPatchi;
1227 newBoundary.resize(nNewPatches);
1228 patchMap_.resize(nNewPatches);
1232 subMeshPtr_().addFvPatches(newBoundary, syncPar);
1273 const label regioni,
label findPatchID(const word &patchName, const bool allowNotFound=true) const
Find patch index given a name, return -1 if not found.
void size(const label n)
Older name for setAddressableSize.
ZoneMesh< faceZone, polyMesh > faceZoneMesh
A ZoneMesh with the type faceZone.
List< cell > cellList
List of cell.
void resize(const label len)
Adjust allocated size of list.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
virtual const labelList & faceNeighbour() const
Return face neighbour.
constexpr char nl
The newline '\n' character (0x0a)
label start() const noexcept
The start label of boundary faces in the polyMesh face list.
bool empty() const noexcept
True if List is empty (ie, size() is zero)
labelList sortedToc() const
The indices of the on bits as a sorted labelList.
Ostream & endl(Ostream &os)
Add newline and flush stream.
static bool & parRun() noexcept
Test if this a parallel run.
const cellList & cells() const
Ignore writing from objectRegistry::writeObject()
static int myProcNo(const label communicator=worldComm)
Rank of this process in the communicator (starting from masterNo()). Can be negative if the process i...
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
UList< label > labelUList
A UList of labels.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
virtual const pointField & points() const
Return raw points.
#define forAll(list, i)
Loop across all elements in list.
List< face > faceList
List of faces.
PtrList< polyPatch > polyPatchList
Store lists of polyPatch as a PtrList.
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator). It is 1 for serial run. ...
vectorField pointField
pointField is a vectorField.
void setSize(const label n)
Alias for resize()
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
bool checkHasSubMesh() const
FatalError if subset has not been performed.
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
A class for handling words, derived from Foam::string.
wordList names() const
Return a list of patch names.
label size() const noexcept
The number of entries in the list.
wordList patchNames(nPatches)
virtual const labelList & faceOwner() const
Return face owner.
void reset()
Reset subMesh and all maps. Same as clear()
label nInternalFaces() const noexcept
Number of internal faces.
ZoneMesh< pointZone, polyMesh > pointZoneMesh
A ZoneMesh with the type pointZone.
virtual const faceList & faces() const
Return raw faces.
errorManip< error > abort(error &err)
const T * set(const label i) const
Return const pointer to element (can be nullptr), or nullptr for out-of-range access (ie...
iterator begin() noexcept
Return an iterator to begin traversing the UList.
bool test(const label pos) const
Test for True value at specified position, never auto-vivify entries.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO...
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
Holds a reference to the original mesh (the baseMesh) and optionally to a subset of that mesh (the su...
static void allGatherList(UList< T > &values, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Gather data, but keep individual values separate. Uses MPI_Allgather or manual linear/tree communicat...
const pointZoneMesh & pointZones() const noexcept
Return point zone mesh.
static word exposedPatchName
Name for exposed internal faces (default: oldInternalFaces)
fvMeshSubset(const fvMeshSubset &)=delete
No copy construct.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Mesh data needed to do the Finite Volume discretisation.
static labelList subsetSubset(const label nElems, const labelUList &selectedElements, const labelUList &subsetMap)
const cellZoneMesh & cellZones() const noexcept
Return cell zone mesh.
messageStream Info
Information stream (stdout output on master, null elsewhere)
A class representing the concept of 0 (zero) that can be used to avoid manipulating objects known to ...
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
List< label > labelList
A List of labels.
static autoPtr< T > New(Args &&... args)
Construct autoPtr with forwarding arguments.
bitSet create(const label n, const labelHashSet &locations, const bool on=true)
Create a bitSet with length n with the specified on locations.
List< bool > boolList
A List of bools.
ZoneMesh< cellZone, polyMesh > cellZoneMesh
A ZoneMesh with the type cellZone.
void reduce(T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) using linear/tree communication schedule.
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
IOporosityModelList pZones(mesh)
void clear()
Reset subMesh and all maps.
static void listCombineReduce(List< T > &values, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Combines List elements. After completion all processors have the same data.
static constexpr const zero Zero
Global zero (0)