Merge branch 'feature-overlapping_zones' into 'develop'

Feature overlapping zones

See merge request Development/openfoam!670
This commit is contained in:
Andrew Heather
2024-03-13 14:44:14 +00:00
8 changed files with 1164 additions and 187 deletions

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2016-2023 OpenCFD Ltd.
Copyright (C) 2016-2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -83,17 +83,77 @@ void Foam::ZoneMesh<ZoneType, MeshType>::calcZoneMap() const
map.reserve(this->totalSize());
const PtrList<ZoneType>& zones = *this;
label zonei = 0;
for (const ZoneType& zn : zones)
{
for (const label id : static_cast<const labelList&>(zn))
{
map.insert(id, zonei);
//map.insert(id, zonei);
const auto fnd = map.cfind(id);
if (!fnd)
{
map.insert(id, zonei);
}
else if (fnd() != zonei)
{
// Multiple zones for same id
if (!additionalMapPtr_)
{
// First occurrence
label maxIndex = -1;
for (const ZoneType& zn : zones)
{
for
(
const label id
: static_cast<const labelList&>(zn)
)
{
maxIndex = max(maxIndex, id);
}
}
additionalMapPtr_.reset(new labelListList(maxIndex+1));
}
auto& additionalMap = *additionalMapPtr_;
additionalMap[id].push_uniq(zonei);
}
}
++zonei;
}
// Sort such that map contains lowest zoneID
if (additionalMapPtr_)
{
auto& additionalMap = *additionalMapPtr_;
forAll(additionalMap, id)
{
labelList& zones = additionalMap[id];
if (zones.size())
{
stableSort(zones);
const label zonei = map[id];
const label index = findLower(zones, zonei);
if (index == -1)
{
// Already first
}
else
{
map.set(id, zones[0]);
for (label i = 0; i < zones.size() && i <= index; i++)
{
zones[i] = zones[i+1];
}
zones[index+1] = zonei;
}
}
}
}
}
}
@ -358,6 +418,35 @@ Foam::label Foam::ZoneMesh<ZoneType, MeshType>::whichZone
}
template<class ZoneType, class MeshType>
Foam::label Foam::ZoneMesh<ZoneType, MeshType>::whichZones
(
const label objectIndex,
DynamicList<label>& zones
) const
{
zones.clear();
const auto fnd = zoneMap().find(objectIndex);
if (fnd)
{
// Add main element
zones.push_back(fnd());
if (additionalMapPtr_)
{
const auto& additionalMap = *additionalMapPtr_;
if (objectIndex < additionalMap.size())
{
for (const label zonei : additionalMap[objectIndex])
{
zones.push_uniq(zonei);
}
}
}
}
return zones.size();
}
template<class ZoneType, class MeshType>
Foam::wordList Foam::ZoneMesh<ZoneType, MeshType>::types() const
{
@ -814,6 +903,7 @@ template<class ZoneType, class MeshType>
void Foam::ZoneMesh<ZoneType, MeshType>::clearLocalAddressing()
{
zoneMapPtr_.reset(nullptr);
additionalMapPtr_.reset(nullptr);
groupIDsPtr_.reset(nullptr);
}

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2016-2023 OpenCFD Ltd.
Copyright (C) 2016-2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -76,6 +76,10 @@ class ZoneMesh
//- Demand-driven: map of zone labels for given element
mutable std::unique_ptr<Map<label>> zoneMapPtr_;
//- Demand-driven: map of additional zone labels for given element.
// Only populated if overlapping zones
mutable std::unique_ptr<labelListList> additionalMapPtr_;
//- Demand-driven: list of zone ids per group
mutable std::unique_ptr<HashTable<labelList>> groupIDsPtr_;
@ -190,6 +194,15 @@ public:
// If object does not belong to any zones, return -1
label whichZone(const label objectIndex) const;
//- Given a global object index, return (in argument) the zones it is
// in. Returns number of zones (0 if object does not belong to any
// zones)
label whichZones
(
const label objectIndex,
DynamicList<label>& zones
) const;
//- Return a list of zone types
wordList types() const;

View File

@ -2698,37 +2698,14 @@ void Foam::polyMeshAdder::add
const polyMesh& mesh = meshes[meshi];
const cellZoneMesh& cellZones = mesh.cellZones();
// Precalc offset zones
labelList newZoneID(mesh.nCells(), -1);
forAll(cellZones, zonei)
{
const labelList& cellLabels = cellZones[zonei];
for (const label celli : cellLabels)
{
if (newZoneID[celli] != -1)
{
WarningInFunction
<< "Cell:" << celli
<< " centre:" << mesh.cellCentres()[celli]
<< " is in two zones:"
<< cellZones[newZoneID[celli]].name()
<< " and " << cellZones[zonei].name() << endl
<< " This is not supported."
<< " Continuing with first zone only." << endl;
}
else
{
newZoneID[celli] = cellZoneMap[meshi][zonei];
}
}
}
// Add cells in mesh order
cellProcAddressing[meshi].setSize(mesh.nCells());
DynamicList<label> zones;
for (label celli = 0; celli < mesh.nCells(); celli++)
{
// Get current zones and renumber
cellZones.whichZones(celli, zones);
inplaceRenumber(cellZoneMap[meshi], zones);
// Add cell from cell
cellProcAddressing[meshi][celli] = meshMod.addCell
(
@ -2736,7 +2713,7 @@ void Foam::polyMeshAdder::add
-1,
-1,
celli,
newZoneID[celli]
zones
);
}
}
@ -2755,30 +2732,21 @@ void Foam::polyMeshAdder::add
const pointField& points = mesh.points();
const pointZoneMesh& pointZones = mesh.pointZones();
// Precalc offset zones
labelList newZoneID(points.size(), -1);
forAll(pointZones, zonei)
{
const labelList& pointLabels = pointZones[zonei];
for (const label pointi : pointLabels)
{
newZoneID[pointi] = pointZoneMap[meshi][zonei];
}
}
// Add points in mesh order
const labelList& myUniquePoints = uniquePoints[meshi];
DynamicList<label> zones;
for (const label pointi : myUniquePoints)
{
// Get current zones and renumber
pointZones.whichZones(pointi, zones);
inplaceRenumber(pointZoneMap[meshi], zones);
// Rewrite the addressing (should already be compact) just in
// case the polyTopoChange does some special ordering
pointProcAddressing[meshi][pointi] = meshMod.addPoint
(
points[pointi],
pointi,
newZoneID[pointi],
zones,
true
);
}
@ -2813,21 +2781,8 @@ void Foam::polyMeshAdder::add
const labelList& faceNeighbour = mesh.faceNeighbour();
const faceZoneMesh& faceZones = mesh.faceZones();
// Precalc offset zones
labelList newZoneID(faces.size(), -1);
boolList zoneFlip(faces.size(), false);
forAll(faceZones, zonei)
{
const labelList& faceLabels = faceZones[zonei];
const boolList& flipMap = faceZones[zonei].flipMap();
forAll(faceLabels, facei)
{
newZoneID[faceLabels[facei]] = faceZoneMap[meshi][zonei];
zoneFlip[faceLabels[facei]] = flipMap[facei];
}
}
DynamicList<label> zones;
DynamicList<bool> flips;
// Add faces in mesh order
@ -2846,14 +2801,28 @@ void Foam::polyMeshAdder::add
{
newFace[fp] = pointProcAddressing[meshi][f[fp]];
}
// Get current zones and renumber
faceZones.whichZones(facei, zones);
flips.clear();
forAll(zones, i)
{
const label zonei = zones[i];
const label index = faceZones[zonei].localID(facei);
flips.append(faceZones[zonei].flipMap()[index]);
}
inplaceRenumber(faceZoneMap[meshi], zones);
bool flipFaceFlux = false;
bool flip = zoneFlip[facei];
if (newNei < newOwn)
{
std::swap(newOwn, newNei);
newFace = newFace.reverseFace();
flipFaceFlux = !flipFaceFlux;
flip = !flip;
for (bool& flip : flips)
{
flip = !flip;
}
}
const label newFacei = meshMod.addFace
@ -2866,8 +2835,8 @@ void Foam::polyMeshAdder::add
facei, // masterFaceID
flipFaceFlux, // flipFaceFlux
-1, // patchID
newZoneID[facei], // zoneID
flip // zoneFlip
zones, // zoneIDs
flips // zoneFlips
);
faceProcAddressing[meshi][facei] = newFacei;
@ -2892,6 +2861,17 @@ void Foam::polyMeshAdder::add
const label newOwn =
cellProcAddressing[meshi][faceOwner[facei]];
faceZones.whichZones(facei, zones);
flips.clear();
forAll(zones, i)
{
const label zonei = zones[i];
const label index = faceZones[zonei].localID(facei);
flips.append(faceZones[zonei].flipMap()[index]);
}
inplaceRenumber(faceZoneMap[meshi], zones);
// Neighbour mesh face
const auto& nbrMesh = meshes[nbrMeshi];
const label nbrFacei = nbrMesh.nInternalFaces()+nbrBFacei;
@ -2912,7 +2892,6 @@ void Foam::polyMeshAdder::add
newFace[fp] = pointProcAddressing[meshi][f[fp]];
}
const bool flipFaceFlux = false;
const bool flip = zoneFlip[facei];
const label newFacei = meshMod.addFace
(
newFace,
@ -2923,8 +2902,8 @@ void Foam::polyMeshAdder::add
facei, // masterFaceID
flipFaceFlux, // flipFaceFlux
-1, // patchID
newZoneID[facei], // zoneID
flip // zoneFlip
zones, // zoneIDs
flips // zoneFlips
);
faceProcAddressing[meshi][facei] = newFacei;
@ -2960,6 +2939,18 @@ void Foam::polyMeshAdder::add
newFace[fp] = pointProcAddressing[meshi][f[fp]];
}
faceZones.whichZones(facei, zones);
flips.clear();
forAll(zones, i)
{
const label zonei = zones[i];
const label index =
faceZones[zonei].localID(facei);
flips.append(faceZones[zonei].flipMap()[index]);
}
inplaceRenumber(faceZoneMap[meshi], zones);
const label newFacei = meshMod.addFace
(
newFace,
@ -2970,8 +2961,8 @@ void Foam::polyMeshAdder::add
facei, // masterFaceID
false, // flipFaceFlux
newPatchi, // patchID
newZoneID[facei], // zoneID
zoneFlip[facei] // zoneFlip
zones, // zoneID
flips // zoneFlip
);
faceProcAddressing[meshi][facei] = newFacei;

File diff suppressed because it is too large Load Diff

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2018-2022 OpenCFD Ltd.
Copyright (C) 2018-2022,2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -55,6 +55,14 @@ Description
- coupled patches: the reorderCoupledFaces routine (borrowed from
the couplePatches utility) reorders coupled patch faces and
uses the cyclicPolyPatch,processorPolyPatch functionality.
- zones are assumed to be non-overlapping by default. If desired to be
overlapping either set the multiZone to true when calling
modifyCell|Face|Point or use the variants of addCell|Face|Point
and modifyCell|Face|Point that take a list of zones.
- if overlapping zones:
- 'main' zone is the lowest numbered zone. -1 means no zones.
- 'additional' zones are stored in incremental ordering (and cannot
contain -1)
SourceFiles
polyTopoChange.C
@ -127,6 +135,9 @@ class polyTopoChange
//- Zone of point
Map<label> pointZone_;
//- Additional zones of points
DynamicList<labelList> pointAdditionalZones_;
//- Retired points
labelHashSet retiredPoints_;
@ -170,6 +181,9 @@ class polyTopoChange
//- Orientation of face in zone
bitSet faceZoneFlip_;
//- Additional zones of face (stored as signed face)
DynamicList<labelList> faceAdditionalZones_;
//- Active faces
label nActiveFaces_;
@ -193,9 +207,12 @@ class polyTopoChange
//- Cells added from face
Map<label> cellFromFace_;
//- Zone of cell
//- First zone of cell (or -1)
DynamicList<label> cellZone_;
//- Additional zones of cell
DynamicList<labelList> cellAdditionalZones_;
// Private Member Functions
@ -207,13 +224,6 @@ class polyTopoChange
DynamicList<Type>& lst
);
template<class Type>
static void reorder
(
const labelUList& oldToNew,
List<DynamicList<Type>>& lst
);
template<class Type>
static void renumberKey
(
@ -533,6 +543,18 @@ public:
const bool inCell
);
//- Add point. Return new point label.
// Notes:
// - masterPointID can be < 0 (appended points)
// - inCell = false: add retired point (to end of point list)
label addPoint
(
const point& pt,
const label masterPointID,
const labelUList& zoneIDs,
const bool inCell
);
//- Modify coordinate.
// Notes:
// - zoneID = +ve (add to zoneID), -ve (remove from zones)
@ -542,12 +564,32 @@ public:
const label pointi,
const point& pt,
const label zoneID,
const bool inCell,
const bool multiZone = false
);
//- Modify coordinate.
// Notes:
// - zoneIDs = set pointZones
// - inCell = false: add retired point (to end of point list)
void modifyPoint
(
const label pointi,
const point& pt,
const labelUList& zoneIDs,
const bool inCell
);
//- Remove/merge point.
void removePoint(const label pointi, const label mergePointi);
//- Get current cellZone(s). Return number of zones.
label pointZones
(
const label pointi,
DynamicList<label>& zones
) const;
//- Add face to cells. Return new face label.
// own,nei<0, zoneID>=0 : add inactive face (to end of face list)
label addFace
@ -564,6 +606,22 @@ public:
const bool zoneFlip
);
//- Add face to cells. Return new face label.
// own,nei<0, zoneID>=0 : add inactive face (to end of face list)
label addFace
(
const face& f,
const label own,
const label nei,
const label masterPointID,
const label masterEdgeID,
const label masterFaceID,
const bool flipFaceFlux,
const label patchID,
const labelUList& zoneIDs,
const UList<bool>& zoneFlips
);
//- Modify vertices or cell of face.
void modifyFace
(
@ -574,12 +632,34 @@ public:
const bool flipFaceFlux,
const label patchID,
const label zoneID,
const bool zoneFlip
const bool zoneFlip,
const bool multiZone = false
);
//- Modify vertices or cell of face.
void modifyFace
(
const face& f,
const label facei,
const label own,
const label nei,
const bool flipFaceFlux,
const label patchID,
const labelUList& zoneIDs,
const UList<bool>& zoneFlips
);
//- Remove/merge face.
void removeFace(const label facei, const label mergeFacei);
//- Get current faceZone(s). Return number of zones.
label faceZones
(
const label facei,
DynamicList<label>& zones,
DynamicList<bool>& flips
) const;
//- Add cell. Return new cell label.
label addCell
(
@ -590,12 +670,34 @@ public:
const label zoneID
);
//- Modify zone of cell
void modifyCell(const label celli, const label zoneID);
//- Add zoned cell (zones cannot be -1). Return new cell label.
label addCell
(
const label masterPointID,
const label masterEdgeID,
const label masterFaceID,
const label masterCellID,
const labelUList& zoneIDs
);
//- Modify zone of cell. Optionally add to zone.
void modifyCell
(
const label celli,
const label zoneID,
const bool multiZone = false
);
//- Set zones of cell
void modifyCell(const label celli, const labelUList& zoneIDs);
//- Remove/merge cell.
void removeCell(const label celli, const label mergeCelli);
//- Get current cellZone(s). Return number of zones.
label cellZones(const label celli, DynamicList<label>& zones) const;
//- Explicitly set the number of patches if construct-without-mesh
//- used.
inline void setNumPatches(const label nPatches);

View File

@ -43,35 +43,13 @@ void Foam::polyTopoChange::reorder
// Create copy
DynamicList<Type> oldLst(lst);
forAll(oldToNew, i)
forAll(oldLst, i)
{
const label newIdx = oldToNew[i];
if (newIdx >= 0)
{
lst[newIdx] = oldLst[i];
}
}
}
template<class Type>
void Foam::polyTopoChange::reorder
(
const labelUList& oldToNew,
List<DynamicList<Type>>& lst
)
{
// Create copy
List<DynamicList<Type>> oldLst(lst);
forAll(oldToNew, i)
{
const label newIdx = oldToNew[i];
if (newIdx >= 0)
{
lst[newIdx].transfer(oldLst[i]);
lst[newIdx] = std::move(oldLst[i]);
}
}
}

View File

@ -42,6 +42,29 @@ filter1
}
heater1
{
// Add some additional water inside a cellZone that is nested inside
// the filter. This is just to demonstrate overlapping cellZones
type scalarSemiImplicitSource;
active yes;
timeStart 0.0;
duration 0.4;
volumeMode specific;
selectionMode cellZone;
cellZone heater;
sources
{
rho (1e-2 0); // kg/s/m^3
H2O (1e-2 0); // kg/s/m^3
}
}
massSource1
{
type scalarSemiImplicitSource;

View File

@ -32,7 +32,6 @@ actions
set filterCells;
}
{
name leftFluidCells;
type cellSet;
@ -114,6 +113,41 @@ actions
source setToFaceZone;
faceSet cycRightFaces; // name of faceSet
}
// Additional zones inside filter
// heater inside filter
{
name heaterCells;
type cellSet;
action new;
source boxToCell;
box (1.7 -10 -10) (1.8 10 10);
}
// Note: could be done directly using boxToCell
{
name heater;
type cellZoneSet;
action new;
source setToCellZone;
set heaterCells;
}
{
name heaterOutsideFaces;
type faceSet;
action new;
source cellToFace;
set heaterCells;
option outside;
}
{
name heaterOutside;
type faceZoneSet;
action new;
source setToFaceZone;
faceSet heaterOutsideFaces; // name of faceSet
}
);