Skip to content

Commit

Permalink
GRIDEDIT-778 Changed made after review
Browse files Browse the repository at this point in the history
  • Loading branch information
@BillSenior
BillSenior committed May 23, 2024
1 parent b995bd7 commit 6ab5ae0
Show file tree
Hide file tree
Showing 5 changed files with 78 additions and 89 deletions.
35 changes: 19 additions & 16 deletions libs/MeshKernel/include/MeshKernel/CasulliDeRefinement.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,17 +57,25 @@ namespace meshkernel
std::vector<Point> ElementsToDelete(const Mesh2D& mesh, const Polygons& polygon) const;

private:
/// @brief Element mask values used in the de-refinement procedure.
/// @brief Maximum number of iterations allowed when initialising the element mask
static const UInt maxIterationCount = 1000;

/// @brief Maximum reserve size for arrays used in the de-refinement
static const UInt maximumSize = 1000;

/// @brief Label for the element in the de-refined grid.
///
/// Enumeration comments are from the original Fortran code.
enum class ElementMask
/// The prefix, e.g. WasNode, indicates the original mesh entity around/from which the de-refined element was constructed.
/// First and after suffix, indicates the order in which the elements are to be processed.
/// Enumeration values and comments are from the original Fortran code.
enum class ElementType
{
A = 1, //< front, 'A' cell (used to be node, delete it): 1
B = 2, //< front, 'B' cell (used to be link, keep it): 2
C = 3, //< 'C' cell (used to be cell, keep it): 3
NotA = -1, //< not in front, 'A' cell: -1
NotB = -2, //< not in front, 'B' cell: -2
Unassigned = 0 //< not assigned a value 0
WasNodeFirst = 1, //< front, 'A' cell (used to be node, delete it): 1
WasEdgeFirst = 2, //< front, 'B' cell (used to be link, keep it): 2
WasCell = 3, //< 'C' cell (used to be cell, keep it): 3
WasNodeAfter = -1, //< not in front, 'A' cell: -1
WasEdgeAfter = -2, //< not in front, 'B' cell: -2
Unassigned = 0 //< not assigned a value 0
};

/// @brief Indicate if the element can be a seed element or not.
Expand Down Expand Up @@ -104,7 +112,7 @@ namespace meshkernel
std::vector<std::array<int, 2>>& kne) const;

/// @brief Initialise the element mask.
std::vector<ElementMask> InitialiseElementMask(const Mesh2D& mesh,
std::vector<ElementType> InitialiseElementType(const Mesh2D& mesh,
const std::vector<int>& nodeTypes) const;

/// \brief Determine if the element can be deleted from the mesh or not.
Expand Down Expand Up @@ -132,7 +140,7 @@ namespace meshkernel
/// @brief Add element id to the list of id's
///
/// only added is it is not already on the list and the element is a quadrilateral
void AddElementToList(const Mesh& mesh, const std::vector<UInt>& frontList, std::vector<UInt>& frontListCopy, const UInt kNew) const;
void AddElementToList(const Mesh& mesh, const std::vector<UInt>& frontList, std::vector<UInt>& frontListCopy, const UInt elementId) const;

/// @brief Redirect nodes of connected cells, deactivate polygons of degree smaller than three
void RedirectNodesOfConnectedElements(Mesh2D& mesh, const UInt elementId, const UInt nodeId, const std::vector<UInt>& indirectlyConnected) const;
Expand All @@ -157,11 +165,6 @@ namespace meshkernel
/// @returns Indicates if the cleanp-up was successful or not
[[nodiscard]] bool CleanUpEdge(Mesh2D& mesh, const UInt edgeId) const;

/// @brief Find the id of the shared node for two edges.
///
/// If no node is shared then return null value
UInt FindCommonNode(const Mesh2D& mesh, const UInt edgeId1, const UInt edgeId2) const;

/// @brief Do the Casullu de-refinement
[[nodiscard]] bool DoDeRefinement(Mesh2D& mesh, const Polygons& polygon) const;

Expand Down
4 changes: 2 additions & 2 deletions libs/MeshKernel/include/MeshKernel/Operations.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -241,13 +241,13 @@ namespace meshkernel
return f1 < f2 ? x1 : x2;
}

/// @brief Get the next forward index.
/// @brief Get the next forward index, for a range in 0 .. size-1.
/// @param[in] currentIndex The current index.
/// @param[in] size The size of the vector.
/// @returns The next forward index.
[[nodiscard]] UInt NextCircularForwardIndex(UInt currentIndex, UInt size);

/// @brief Get the next backward index.
/// @brief Get the next backward index, for a range in 0 .. size-1.
/// @param[in] currentIndex The current index.
/// @param[in] size The size of the vector.
/// @returns The next backward index.
Expand Down
100 changes: 37 additions & 63 deletions libs/MeshKernel/src/CasulliDeRefinement.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,11 +79,7 @@ void meshkernel::CasulliDeRefinement::FindDirectlyConnectedCells(const Mesh2D& m
}

UInt neighbouringElementId = elementId == mesh.m_edgesFaces[edgeId][0] ? mesh.m_edgesFaces[edgeId][1] : mesh.m_edgesFaces[edgeId][0];

if (std::ranges::find(connected, neighbouringElementId) == connected.end())
{
connected.push_back(neighbouringElementId);
}
connected.push_back(neighbouringElementId);
}
}

Expand Down Expand Up @@ -290,28 +286,26 @@ void meshkernel::CasulliDeRefinement::AddElementToList(const Mesh& mesh, const s
}
}

std::vector<meshkernel::CasulliDeRefinement::ElementMask> meshkernel::CasulliDeRefinement::InitialiseElementMask(const Mesh2D& mesh,
std::vector<meshkernel::CasulliDeRefinement::ElementType> meshkernel::CasulliDeRefinement::InitialiseElementType(const Mesh2D& mesh,
const std::vector<int>& nodeTypes) const
{
UInt seedElement = FindElementSeedIndex(mesh, nodeTypes);
UInt iterationCount = 0;
UInt nMax = 1000; // fix
UInt maxIterationCount = 1000; // fix

std::vector<UInt> directlyConnected;
std::vector<UInt> indirectlyConnected;
std::vector<std::array<int, 2>> kne(nMax);
std::vector<std::array<int, 2>> kne(maximumSize);
std::vector<UInt> frontIndex;
std::vector<UInt> frontIndexCopy;

directlyConnected.reserve(nMax);
indirectlyConnected.reserve(nMax);
frontIndex.reserve(nMax);
frontIndexCopy.reserve(nMax);
directlyConnected.reserve(maximumSize);
indirectlyConnected.reserve(maximumSize);
frontIndex.reserve(maximumSize);
frontIndexCopy.reserve(maximumSize);

std::vector<ElementMask> cellMask(mesh.GetNumFaces(), ElementMask::Unassigned);
std::vector<ElementType> cellMask(mesh.GetNumFaces(), ElementType::Unassigned);

cellMask[seedElement] = ElementMask::A;
cellMask[seedElement] = ElementType::WasNodeFirst;
frontIndex.push_back(seedElement);

while (frontIndex.size() > 0 && iterationCount < maxIterationCount)
Expand All @@ -326,7 +320,7 @@ std::vector<meshkernel::CasulliDeRefinement::ElementMask> meshkernel::CasulliDeR
// get the connected cells
FindSurroundingCells(mesh, elementId, directlyConnected, indirectlyConnected, kne);

if (cellMask[elementId] == ElementMask::A)
if (cellMask[elementId] == ElementType::WasNodeFirst)
{

for (UInt j = 0; j < directlyConnected.size(); ++j)
Expand All @@ -338,10 +332,10 @@ std::vector<meshkernel::CasulliDeRefinement::ElementMask> meshkernel::CasulliDeR
continue;
}

if ((cellMask[connectedElementId] != ElementMask::A && cellMask[connectedElementId] != ElementMask::NotA) &&
(cellMask[connectedElementId] != ElementMask::B && cellMask[connectedElementId] != ElementMask::NotB))
if ((cellMask[connectedElementId] != ElementType::WasNodeFirst && cellMask[connectedElementId] != ElementType::WasNodeAfter) &&
(cellMask[connectedElementId] != ElementType::WasEdgeFirst && cellMask[connectedElementId] != ElementType::WasEdgeAfter))
{
cellMask[connectedElementId] = ElementMask::B;
cellMask[connectedElementId] = ElementType::WasEdgeFirst;
AddElementToList(mesh, frontIndex, frontIndexCopy, connectedElementId);
}
}
Expand All @@ -355,15 +349,15 @@ std::vector<meshkernel::CasulliDeRefinement::ElementMask> meshkernel::CasulliDeR
continue;
}

if (cellMask[connectedElementId] != ElementMask::C)
if (cellMask[connectedElementId] != ElementType::WasCell)
{
cellMask[connectedElementId] = ElementMask::C;
cellMask[connectedElementId] = ElementType::WasCell;
}
}

cellMask[elementId] = ElementMask::NotA;
cellMask[elementId] = ElementType::WasNodeAfter;
}
else if (cellMask[elementId] == ElementMask::B)
else if (cellMask[elementId] == ElementType::WasEdgeFirst)
{

for (UInt j = 0; j < directlyConnected.size(); ++j)
Expand All @@ -375,11 +369,11 @@ std::vector<meshkernel::CasulliDeRefinement::ElementMask> meshkernel::CasulliDeR
continue;
}

if ((cellMask[connectedElementId] != ElementMask::C) &&
(cellMask[connectedElementId] != ElementMask::A && cellMask[connectedElementId] != ElementMask::NotA) &&
(cellMask[connectedElementId] != ElementMask::B && cellMask[connectedElementId] != ElementMask::NotB))
if ((cellMask[connectedElementId] != ElementType::WasCell) &&
(cellMask[connectedElementId] != ElementType::WasNodeFirst && cellMask[connectedElementId] != ElementType::WasNodeAfter) &&
(cellMask[connectedElementId] != ElementType::WasEdgeFirst && cellMask[connectedElementId] != ElementType::WasEdgeAfter))
{
cellMask[connectedElementId] = ElementMask::A;
cellMask[connectedElementId] = ElementType::WasNodeFirst;
AddElementToList(mesh, frontIndex, frontIndexCopy, connectedElementId);
}
}
Expand All @@ -393,16 +387,16 @@ std::vector<meshkernel::CasulliDeRefinement::ElementMask> meshkernel::CasulliDeR
continue;
}

if ((cellMask[connectedElementId] != ElementMask::B && cellMask[connectedElementId] != ElementMask::NotB) &&
(cellMask[connectedElementId] != ElementMask::A && cellMask[connectedElementId] != ElementMask::NotA) &&
cellMask[connectedElementId] != ElementMask::C)
if ((cellMask[connectedElementId] != ElementType::WasEdgeFirst && cellMask[connectedElementId] != ElementType::WasEdgeAfter) &&
(cellMask[connectedElementId] != ElementType::WasNodeFirst && cellMask[connectedElementId] != ElementType::WasNodeAfter) &&
cellMask[connectedElementId] != ElementType::WasCell)
{
cellMask[connectedElementId] = ElementMask::B;
cellMask[connectedElementId] = ElementType::WasEdgeFirst;
AddElementToList(mesh, frontIndex, frontIndexCopy, connectedElementId);
}
}

cellMask[elementId] = ElementMask::NotB;
cellMask[elementId] = ElementType::WasEdgeAfter;
}
}

Expand All @@ -414,20 +408,19 @@ std::vector<meshkernel::CasulliDeRefinement::ElementMask> meshkernel::CasulliDeR

bool meshkernel::CasulliDeRefinement::DoDeRefinement(Mesh2D& mesh, const Polygons& polygon) const
{
UInt nMax = 1000;
std::vector<UInt> directlyConnected;
std::vector<UInt> indirectlyConnected;
std::vector<std::array<int, 2>> kne(nMax);
std::vector<std::array<int, 2>> kne(maximumSize);
std::vector<int> nodeTypes(ComputeNodeTypes(mesh, polygon));
directlyConnected.reserve(nMax);
indirectlyConnected.reserve(nMax);
directlyConnected.reserve(maximumSize);
indirectlyConnected.reserve(maximumSize);

std::vector<ElementMask> cellMask(InitialiseElementMask(mesh, nodeTypes));
std::vector<ElementType> cellMask(InitialiseElementType(mesh, nodeTypes));
mesh.ComputeCircumcentersMassCentersAndFaceAreas(true);

for (UInt k = 0; k < cellMask.size(); ++k)
{
if (cellMask[k] == ElementMask::NotA && mesh.m_numFacesNodes[k] > 0)
if (cellMask[k] == ElementType::WasNodeAfter && mesh.m_numFacesNodes[k] > 0)
{
FindSurroundingCells(mesh, k, directlyConnected, indirectlyConnected, kne);

Expand Down Expand Up @@ -678,7 +671,7 @@ bool meshkernel::CasulliDeRefinement::UpdateDirectlyConnectedElements(Mesh2D& me
}
else
{
// Error: No node found
throw AlgorithmError("No node found in Casulli de-refinement");
}

mesh.m_numFacesNodes[connectedElementId] = ndum;
Expand Down Expand Up @@ -748,17 +741,16 @@ bool meshkernel::CasulliDeRefinement::RemoveUnwantedBoundaryNodes(Mesh2D& mesh,

for (UInt i = 0; i < mesh.m_numFacesNodes[connectedElementId]; ++i)
{
// UInt nodeId = constants::missing::uintValue;// = mesh.m_facesNodes[connectedElementId][i];
UInt edgeId = mesh.m_facesEdges[connectedElementId][i];

if (mesh.m_edgesNumFaces[edgeId] == 1)
{
UInt im1 = RotateIndex(i, mesh.m_facesNodes[connectedElementId], false /*forward*/);
UInt im1 = NextCircularBackwardIndex(i, static_cast<UInt>(mesh.m_facesNodes[connectedElementId].size()));
UInt previousEdgeId = mesh.m_facesEdges[connectedElementId][im1];

if (mesh.m_edgesNumFaces[previousEdgeId] == 1)
{
UInt nodeId = FindCommonNode(mesh, edgeId, previousEdgeId);
UInt nodeId = mesh.FindCommonNode(edgeId, previousEdgeId);

// [sic] weird
if (nodeId == constants::missing::uintValue)
Expand Down Expand Up @@ -912,24 +904,6 @@ bool meshkernel::CasulliDeRefinement::DeleteElement(Mesh2D& mesh,
return true;
}

meshkernel::UInt meshkernel::CasulliDeRefinement::FindCommonNode(const Mesh2D& mesh, const UInt edgeId1, const UInt edgeId2) const
{
const Edge& e1 = mesh.GetEdge(edgeId1);
const Edge& e2 = mesh.GetEdge(edgeId2);

if (e1.first == e2.first || e1.first == e2.second)
{
return e1.first;
}

if (e1.second == e2.first || e1.second == e2.second)
{
return e1.second;
}

return constants::missing::uintValue;
}

bool meshkernel::CasulliDeRefinement::CleanUpEdge(Mesh2D& mesh, const UInt edgeId) const
{

Expand Down Expand Up @@ -974,7 +948,7 @@ bool meshkernel::CasulliDeRefinement::CleanUpEdge(Mesh2D& mesh, const UInt edgeI

std::vector<int> meshkernel::CasulliDeRefinement::ComputeNodeTypes(const Mesh2D& mesh, const Polygons& polygon) const
{
std::vector<int> nodeTypes(mesh.GetNumNodes());
std::vector<int> nodeTypes(mesh.GetNumNodes(), 0);

for (UInt i = 0; i < mesh.GetNumNodes(); ++i)
{
Expand All @@ -990,13 +964,13 @@ std::vector<int> meshkernel::CasulliDeRefinement::ComputeNodeTypes(const Mesh2D&
std::vector<meshkernel::Point> meshkernel::CasulliDeRefinement::ElementsToDelete(const Mesh2D& mesh, const Polygons& polygon) const
{
std::vector<int> nodeTypes(ComputeNodeTypes(mesh, polygon));
std::vector<ElementMask> cellMask(InitialiseElementMask(mesh, nodeTypes));
std::vector<ElementType> cellMask(InitialiseElementType(mesh, nodeTypes));
std::vector<Point> elementCentres;
elementCentres.reserve(cellMask.size());

for (UInt k = 0; k < cellMask.size(); ++k)
{
if (cellMask[k] == ElementMask::NotA && mesh.m_numFacesNodes[k] > 0)
if (cellMask[k] == ElementType::WasNodeAfter && mesh.m_numFacesNodes[k] > 0)
{
bool toDelete = false;

Expand Down
3 changes: 2 additions & 1 deletion libs/MeshKernel/src/MeshRefinement.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -824,7 +824,8 @@ void MeshRefinement::FindHangingNodes(UInt face)
std::fill(m_isHangingNodeCache.begin(), m_isHangingNodeCache.end(), false);
std::fill(m_isHangingEdgeCache.begin(), m_isHangingEdgeCache.end(), false);

auto kknod = numFaceNodes;
auto kknod = numFaceNodes - 1;

for (UInt n = 0; n < numFaceNodes; n++)
{
const auto edgeIndex = m_mesh.m_facesEdges[face][n];
Expand Down
25 changes: 18 additions & 7 deletions libs/MeshKernel/src/Operations.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,21 +115,32 @@ namespace meshkernel

UInt NextCircularForwardIndex(UInt currentIndex, UInt size)
{
UInt index = currentIndex + 1;
if (index >= size)
if (size == 0)
{
index = index - size;
throw ConstraintError("Invalid rotation range ");
}
return index;

if (currentIndex >= size)
{
throw ConstraintError("Index is out of range: {} not in [0 .. {}]", currentIndex, size - 1);
}

return currentIndex == size - 1 ? 0 : currentIndex + 1;
}

UInt NextCircularBackwardIndex(UInt currentIndex, UInt size)
{
if (currentIndex == 0)
if (size == 0)
{
return currentIndex + size - 1;
throw ConstraintError("Invalid rotation range ");
}
return currentIndex - 1;

if (currentIndex >= size)
{
throw ConstraintError("Index is out of range: {} not in [0 .. {}]", currentIndex, size - 1);
}

return currentIndex == 0 ? size - 1 : currentIndex - 1;
}

bool IsPointOnPole(const Point& point)
Expand Down

0 comments on commit 6ab5ae0

Please sign in to comment.