From a889942ae8cd3ccbd985cdf55f5453ed248ee695 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 25 Sep 2025 15:54:31 +0200 Subject: [PATCH 01/16] GRIDEDIT-1795 updated comment for edge --- libs/MeshKernel/include/MeshKernel/Entities.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libs/MeshKernel/include/MeshKernel/Entities.hpp b/libs/MeshKernel/include/MeshKernel/Entities.hpp index 5d665218c4..6d539dfedb 100644 --- a/libs/MeshKernel/include/MeshKernel/Entities.hpp +++ b/libs/MeshKernel/include/MeshKernel/Entities.hpp @@ -46,7 +46,7 @@ namespace meshkernel // t.IsValid //}; - /// @brief Describes an edge with two indices + /// @brief Describes an edge connecting two nodes. using Edge = std::pair; //// @brief Determine if the value is a valid value (true) or is the undefined value (false) From 769b62229955b770873dbccd734d7fda1a73a985 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 29 Sep 2025 17:51:00 +0200 Subject: [PATCH 02/16] GRIDEDIT-1795 New approach to illegal cells that uses mostly existing code for finding the domain boundary and deleting elements --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 2 + libs/MeshKernel/src/Mesh2D.cpp | 111 +++++++++++++++ .../tests/src/MeshRefinementTests.cpp | 132 ++++++++++++++++++ 3 files changed, 245 insertions(+) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 7551be5e6c..64003f94dd 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -378,6 +378,8 @@ namespace meshkernel /// If no such element can be found then the null value will be returned. UInt FindCommonFace(const UInt edge1, const UInt edge2) const; + [[nodiscard]] std::unique_ptr DeleteMeshFaces2(const Polygons& polygon, bool invertDeletion); + private: // orthogonalization static constexpr double m_minimumEdgeLength = 1e-4; ///< Minimum edge length diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index e92a23c76c..ec4df4e6bb 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -1870,6 +1870,117 @@ std::unique_ptr Mesh2D::DeleteMeshFaces(const Polygons& return deleteMeshAction; } +std::unique_ptr Mesh2D::DeleteMeshFaces2(const Polygons& polygon, bool invertDeletion) +{ + std::unique_ptr deleteMeshAction = CompoundUndoAction::Create(); + + Administrate(deleteMeshAction.get()); + // std::vector faceCircumcenters = algo::ComputeFaceCircumcenters(*this); + std::vector faceCircumcenters = m_facesMassCenters; + + bool isInPolygon; + UInt polygonIndex; + + for (UInt f = 0u; f < GetNumFaces(); ++f) + { + isInPolygon = false; + polygonIndex = constants::missing::uintValue; + + // const auto faceIndex = m_edgesFaces[e][f]; + + if (!faceCircumcenters[f].IsValid()) + { + continue; + } + + std::tie(isInPolygon, polygonIndex) = polygon.IsPointInPolygons(faceCircumcenters[f]); + + std::cout << "centre is in polygon: " << f << " " << std::boolalpha << isInPolygon << " " << polygonIndex << " " << faceCircumcenters[f].x << ", " << faceCircumcenters[f].y << std::endl; + + if (invertDeletion) + { + isInPolygon = !isInPolygon; + } + + if (isInPolygon) + { + + for (UInt e = 0u; e < m_facesEdges[f].size(); ++e) + { + UInt edge = m_facesEdges[f][e]; + + if (m_edgesFaces[edge][0] == f) + { + m_edgesFaces[edge][0] = constants::missing::uintValue; + --m_edgesNumFaces[edge]; + } + else if (m_edgesFaces[edge][1] == f) + { + m_edgesFaces[edge][1] = constants::missing::uintValue; + --m_edgesNumFaces[edge]; + } + } + + m_facesNodes[f].clear(); + m_numFacesNodes[f] = 0; + m_facesEdges[f].clear(); + m_facesMassCenters[f] = {constants::missing::doubleValue, constants::missing::doubleValue}; + m_faceArea[f] = constants::missing::doubleValue; + + // deleteMeshAction->Add(DeleteEdge(e)); + } + } + + // for (UInt e = 0u; e < GetNumEdges(); ++e) + // { + // for (UInt f = 0u; f < GetNumEdgesFaces(e); ++f) + // { + // isInPolygon = false; + // polygonIndex = constants::missing::uintValue; + + // const auto faceIndex = m_edgesFaces[e][f]; + // if (faceIndex == constants::missing::uintValue) + // { + // continue; + // } + + // std::tie(isInPolygon, polygonIndex) = polygon.IsPointInPolygons(faceCircumcenters[faceIndex]); + + // std::cout << "centre is in polygon: " << std::boolalpha << isInPolygon << " " << polygonIndex << " " << faceCircumcenters[faceIndex].x << ", " << faceCircumcenters[faceIndex].y << std::endl; + + // if (invertDeletion) + // { + // isInPolygon = !isInPolygon; + // } + + // if (isInPolygon) + // { + + // if (m_edgesFaces[e][0] == f) + // { + // m_edgesFaces[e][0] = constants::missing::uintValue; + // --m_edgesNumFaces[e]; + // } + // else if (m_edgesFaces[e][1] == f) + // { + // m_edgesFaces[e][1] = constants::missing::uintValue; + // --m_edgesNumFaces[e]; + // } + + // m_facesNodes[f].clear(); + // m_numFacesNodes[f] = 0; + // m_facesEdges[f].clear(); + // m_facesMassCenters[f] = {constants::missing::doubleValue, constants::missing::doubleValue}; + // m_faceArea[f] = constants::missing::doubleValue; + + // // deleteMeshAction->Add(DeleteEdge(e)); + // } + // } + // } + + return deleteMeshAction; +} + std::unique_ptr Mesh2D::DeleteHangingEdges() { std::unique_ptr deleteAction = CompoundUndoAction::Create(); diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index 4668e2f3f5..be45ab2d44 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -49,6 +49,9 @@ #include "MeshKernel/UndoActions/UndoActionStack.hpp" #include "MeshKernel/Utilities/Utilities.hpp" +#include "MeshKernel/LandBoundaries.hpp" +#include "MeshKernel/OrthogonalizationAndSmoothing.hpp" + #include "TestUtils/Definitions.hpp" #include "TestUtils/MakeCurvilinearGrids.hpp" #include "TestUtils/MakeMeshes.hpp" @@ -2676,3 +2679,132 @@ TEST(MeshRefinement, RowSplittingFailureTests) // Invalid edge id EXPECT_THROW([[maybe_unused]] auto undo4 = splitMeshRow.Compute(mesh, edgeId), ConstraintError); } + +TEST(MeshRefinement, WTF) +{ + const size_t ExpectedNumberOfPoints = 184; + const size_t ExpectedNumberOfEdges = 360; + + auto curviMesh = MakeCurvilinearGrid(0.0, 0.0, 20.0, 20.0, 11, 11); + const auto edges = curviMesh->ComputeEdges(); + const auto nodes = curviMesh->ComputeNodes(); + + Mesh2D mesh(edges, nodes, Projection::cartesian); + + const std::vector originalNodes(mesh.Nodes()); + const std::vector originalEdges(mesh.Edges()); + + std::vector patch{{45.0, 5.0}, + {155.0, 5.0}, + {155.0, 155.0}, + {45.0, 155.0}, + {45.0, 5.0}, + {constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + {65.0, 65.0}, + {115.0, 65.0}, + {115.0, 115.0}, + {65.0, 115.0}, + {65.0, 65.0}}; + + std::vector expectedPoints(ExpectedNumberOfPoints); + std::vector expectedEdgeStart(ExpectedNumberOfEdges); + std::vector expectedEdgeEnd(ExpectedNumberOfEdges); + + Polygons polygon(patch, Projection::cartesian); + + auto undoAction = CasulliRefinement::Compute(mesh, polygon); + + auto node11 = mesh.FindNodeCloseToAPoint({80.0, 0.0}, 1.0e-5); + auto node12 = mesh.FindNodeCloseToAPoint({85.0, 15.0}, 1.0e-5); + auto node13 = mesh.FindNodeCloseToAPoint({75.0, 15.0}, 1.0e-5); + + auto node21 = mesh.FindNodeCloseToAPoint({180.0, 140.0}, 1.0e-5); + auto node22 = mesh.FindNodeCloseToAPoint({200.0, 140.0}, 1.0e-5); + auto node23 = mesh.FindNodeCloseToAPoint({200.0, 160.0}, 1.0e-5); + auto node24 = mesh.FindNodeCloseToAPoint({180.0, 160.0}, 1.0e-5); + + std::cout << "nodes2: " << node21 << " " << node22 << " " << node23 << " " << node24 << std::endl; + + auto node31 = mesh.FindNodeCloseToAPoint({125.0, 125.0}, 1.0e-5); + auto node32 = mesh.FindNodeCloseToAPoint({135.0, 125.0}, 1.0e-5); + auto node33 = mesh.FindNodeCloseToAPoint({135.0, 135.0}, 1.0e-5); + auto node34 = mesh.FindNodeCloseToAPoint({125.0, 135.0}, 1.0e-5); + + std::cout << "nodes3: " << node31 << " " << node32 << " " << node33 << " " << node34 << std::endl; + + std::vector polygonNodes; + auto boundaryNodes = mesh.ComputeBoundaryPolygons(polygonNodes); + + // std::vector boundaryNodes = meshBoundaryPolygon.GatherAllEnclosureNodes(); + + std::vector elementNodes1{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node11), + mesh.Node(node12), + mesh.Node(node13), + mesh.Node(node11)}; + + std::vector elementNodes2{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node21), + mesh.Node(node22), + mesh.Node(node23), + mesh.Node(node24), + mesh.Node(node21)}; + + std::vector elementNodes3{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node31), + mesh.Node(node32), + mesh.Node(node33), + mesh.Node(node34), + mesh.Node(node31)}; + + boundaryNodes.insert(boundaryNodes.end(), elementNodes1.begin(), elementNodes1.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes2.begin(), elementNodes2.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes3.begin(), elementNodes3.end()); + + [[maybe_unused]] Polygons boundaryWithMissingElement(boundaryNodes, Projection::cartesian); + + // auto edge1 = mesh.FindEdge(node1, node2); + // auto edge2 = mesh.FindEdge(node2, node3); + // auto edge3 = mesh.FindEdge(node3, node1); + + // auto undoAction3 = ::Compute(mesh, Polygons()); + + MeshRefinementParameters refParams; + refParams.max_num_refinement_iterations = 1; + + MeshRefinement meshRefinement(mesh, Polygons(), refParams); + auto undoAction3 = meshRefinement.Compute(); + + //-------------------------------- + + const auto projectToLandBoundaryOption = LandBoundaries::ProjectToLandBoundaryOption::DoNotProjectToLandBoundary; + OrthogonalizationParameters orthogonalizationParameters; + orthogonalizationParameters.inner_iterations = 2; + orthogonalizationParameters.boundary_iterations = 25; + orthogonalizationParameters.outer_iterations = 25; + orthogonalizationParameters.orthogonalization_to_smoothing_factor = 0.975; + orthogonalizationParameters.orthogonalization_to_smoothing_factor_at_boundary = 0.975; + orthogonalizationParameters.areal_to_angle_smoothing_factor = 1.0; + + // Execute + auto orthoPolygon = std::make_unique(); + + std::vector landBoundary{}; + auto landboundaries = std::make_unique(landBoundary, mesh, *orthoPolygon); + + OrthogonalizationAndSmoothing orthogonalization(mesh, + std::move(orthoPolygon), + std::move(landboundaries), + projectToLandBoundaryOption, + orthogonalizationParameters); + + [[maybe_unused]] auto undoActionOrtho = orthogonalization.Initialize(); + orthogonalization.Compute(); + + //-------------------------------- + + auto undoAction2 = mesh.DeleteMeshFaces2(boundaryWithMissingElement, true); + + meshkernel::SaveVtk(mesh.Nodes(), mesh.m_facesNodes, "/home/wcs1/MeshKernel/MeshKernel04/build_deb/mesh1.vtu"); + // meshkernel::Print(mesh.Nodes(), mesh.Edges()); +} From df1ea3231fa0c94130a2a51837bdd7dc18600f3a Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 2 Oct 2025 18:26:40 +0200 Subject: [PATCH 03/16] GRIDEDIT-1795 Found inner boundaries for holes --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 34 +- libs/MeshKernel/src/Mesh2D.cpp | 244 ++++++++++++- libs/MeshKernel/src/MeshRefinement.cpp | 5 + .../tests/src/MeshRefinementTests.cpp | 336 +++++++++++++++++- 4 files changed, 601 insertions(+), 18 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 64003f94dd..e3ed93a8b2 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -244,15 +244,10 @@ namespace meshkernel /// @return The resulting polygon mesh boundary [[nodiscard]] std::vector ComputeBoundaryPolygons(const std::vector& polygon); - /// @brief Constructs a polygon from the meshboundary, by walking through the mesh - /// @param[in] polygon The input polygon - /// @param[in,out] isVisited the visited mesh nodes - /// @param[in,out] currentNode the current node - /// @param[out] meshBoundaryPolygon The resulting polygon points - void WalkBoundaryFromNode(const Polygon& polygon, - std::vector& isVisited, - UInt& currentNode, - std::vector& meshBoundaryPolygon) const; + /// @brief Convert all mesh boundaries to a vector of polygon nodes + /// @param[in] polygon The polygon where the operation is performed + /// @return The resulting set of polygons, describing interior mesh boundaries + std::vector ComputeInnerBoundaryPolygons() const; /// @brief Gets the hanging edges /// @return A vector with the indices of the hanging edges @@ -437,6 +432,27 @@ namespace meshkernel /// @brief Find the mesh faces that lie entirely within the polygon. std::vector FindFacesEntirelyInsidePolygon(const std::vector& isNodeInsidePolygon) const; + /// @brief Constructs a polygon from the meshboundary, by walking through the mesh + void WalkBoundaryFromNode(const Polygon& polygon, + std::vector& isVisited, + UInt& currentNode, + std::vector& meshBoundaryPolygon) const; + + /// @brief Constructs a polygon or polygons from the meshboundary, by walking through the mesh + /// + /// If there are multiple polygons connected by a single node, then these will be separated into indicidual polygons + void WalkBoundaryFromNode(std::vector& edgeIsVisited, + std::vector& nodeIsVisited, + UInt& currentNode, + std::vector& meshBoundaryPolygon, + std::vector& nodeIds, + std::vector& subsSequence) const; + + /// @brief Removes the outer domain boundary polygon from the set of polygons + /// + /// It is assumed that the outer domain polygon contains the most nodes + std::vector RemoveOuterDomainBoundaryPolygon(const std::vector& polygonNodes) const; + /// @brief Deletes the mesh faces inside a polygon /// @param[in] polygon The polygon where to perform the operation /// If this Polygons instance contains multiple polygons, the first one will be taken. diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index ec4df4e6bb..07b1d300c0 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -1485,7 +1485,7 @@ std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector const Polygon polygon(polygonNodes, m_projection); // Find faces - Administrate(); + // Administrate(); std::vector isVisited(GetNumEdges(), false); std::vector meshBoundaryPolygon; meshBoundaryPolygon.reserve(GetNumNodes()); @@ -1513,6 +1513,7 @@ std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector // Start a new polyline if (!meshBoundaryPolygon.empty()) { + // meshBoundaryPolygon.emplace_back(constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator); meshBoundaryPolygon.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); } @@ -1552,6 +1553,140 @@ std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector return meshBoundaryPolygon; } +std::vector Mesh2D::ComputeInnerBoundaryPolygons() const +{ + + std::vector meshBoundaryPolygon; + meshBoundaryPolygon.reserve(GetNumNodes()); + std::vector subSqeuence; + subSqeuence.reserve(GetNumNodes()); + + std::vector edgeIsVisited(GetNumEdges(), false); + std::vector nodeIsVisited(GetNumNodes(), false); + + std::vector nodeIds; + nodeIds.reserve(GetNumNodes()); + + for (UInt e = 0; e < GetNumEdges(); e++) + { + if (edgeIsVisited[e] || !IsEdgeOnBoundary(e)) + { + continue; + } + + const auto firstNodeIndex = m_edges[e].first; + const auto secondNodeIndex = m_edges[e].second; + const auto firstNode = m_nodes[firstNodeIndex]; + const auto secondNode = m_nodes[secondNodeIndex]; + + // Start a new polyline + if (!subSqeuence.empty()) + { + subSqeuence.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); + nodeIds.emplace_back(constants::missing::uintValue); + } + + // Put the current edge on the mesh boundary, mark it as visited + const auto startPolygonEdges = static_cast(subSqeuence.size()); + subSqeuence.emplace_back(firstNode); + subSqeuence.emplace_back(secondNode); + nodeIds.emplace_back(firstNodeIndex); + nodeIds.emplace_back(secondNodeIndex); + edgeIsVisited[e] = true; + nodeIsVisited[firstNodeIndex] = true; + nodeIsVisited[secondNodeIndex] = true; + + // walk the current mesh boundary + auto currentNode = secondNodeIndex; + WalkBoundaryFromNode(edgeIsVisited, nodeIsVisited, currentNode, subSqeuence, nodeIds, meshBoundaryPolygon); + + const auto numNodesFirstTail = static_cast(subSqeuence.size()); + + // if the boundary polygon is not closed + if (currentNode != firstNodeIndex) + { + // Now grow a polyline starting at the other side of the original link L, i.e., the second tail + currentNode = firstNodeIndex; + WalkBoundaryFromNode(edgeIsVisited, nodeIsVisited, currentNode, subSqeuence, nodeIds, meshBoundaryPolygon); + } + + // There is a nonempty second tail, so reverse the first tail, so that they connect. + if (subSqeuence.size() > numNodesFirstTail) + { + const auto start = startPolygonEdges + static_cast(std::ceil((numNodesFirstTail - startPolygonEdges + static_cast(1)) * 0.5)); + for (auto n = start; n < numNodesFirstTail; n++) + { + const auto backupPoint = subSqeuence[n]; + const auto replaceIndex = numNodesFirstTail - n + firstNodeIndex; + subSqeuence[n] = subSqeuence[replaceIndex]; + subSqeuence[replaceIndex] = backupPoint; + + const UInt backupPointIndex = nodeIds[n]; + nodeIds[n] = nodeIds[replaceIndex]; + nodeIds[replaceIndex] = backupPointIndex; + } + } + } + + return RemoveOuterDomainBoundaryPolygon(meshBoundaryPolygon); +} + +std::vector Mesh2D::RemoveOuterDomainBoundaryPolygon(const std::vector& polygonNodes) const +{ + // Remove outer boundary. + // It is assumed that the boundary polygon with the most points is from the outer domain boundary + + UInt polygonLength = 0; + UInt polygonStart = 0; + + UInt polygonLengthIntermediate = 0; + UInt polygonStartIntermediate = 0; + + UInt index = 0; + + // Find start index of outer boundary and the length of the polygon + while (index < polygonNodes.size()) + { + polygonStartIntermediate = index; + + for (UInt i = polygonStartIntermediate; i < polygonNodes.size(); ++i) + { + ++index; + + if (!polygonNodes[i].IsValid()) + { + polygonLengthIntermediate = i - polygonStartIntermediate; + break; + } + } + + if (polygonLengthIntermediate > polygonLength) + { + polygonLength = polygonLengthIntermediate; + polygonStart = polygonStartIntermediate; + } + } + + if (polygonLength > 0 && polygonLength < polygonNodes.size()) + { + ++polygonLength; + } + + std::vector innerBoundaryNodes; + + if (polygonStart > 0) + { + innerBoundaryNodes.insert(innerBoundaryNodes.begin(), polygonNodes.begin(), polygonNodes.begin() + polygonStart - 1); + } + + if (polygonLength > 0 && polygonStart + polygonLength < polygonNodes.size()) + { + innerBoundaryNodes.insert(innerBoundaryNodes.end(), polygonNodes.begin() + polygonStart + polygonLength, polygonNodes.end()); + } + + return innerBoundaryNodes; +} + void Mesh2D::WalkBoundaryFromNode(const Polygon& polygon, std::vector& isVisited, UInt& currentNode, @@ -1587,6 +1722,106 @@ void Mesh2D::WalkBoundaryFromNode(const Polygon& polygon, } } +void Mesh2D::WalkBoundaryFromNode(std::vector& edgeIsVisited, + std::vector& nodeIsVisited, + UInt& currentNode, + std::vector& meshBoundaryPolygon, + std::vector& nodeIds, + std::vector& subSequence) const +{ + std::cout << std::endl + << " Mesh2D::WalkBoundaryFromNode " << currentNode << " " << m_nodes[currentNode].x << ", " << m_nodes[currentNode].y + << " " << subSequence.size() + << std::endl; + UInt e = 0; + + // size_t start = meshBoundaryPolygon.size(); + + while (e < m_nodesNumEdges[currentNode]) + { + const auto currentEdge = m_nodesEdges[currentNode][e]; + + if (edgeIsVisited[currentEdge] || !IsEdgeOnBoundary(currentEdge)) + { + e++; + continue; + } + + UInt nextNode = OtherNodeOfEdge(m_edges[currentEdge], currentNode); + e = 0; + + if (nodeIsVisited[nextNode]) + { + UInt lastIndex = constants::missing::uintValue; + + // Find index of last time node was added + for (size_t ii = nodeIds.size(); ii >= 1; --ii) + { + size_t i = ii - 1; + + if (nodeIds[i] == constants::missing::uintValue) + { + break; + } + + if (nodeIds[i] == nextNode) + { + lastIndex = i; + break; + } + } + + std::cout << "node ids: "; + + for (size_t i = lastIndex; i < nodeIds.size(); ++i) + { + std::cout << nodeIds[i] << " "; + } + + std::cout << std::endl; + + if (lastIndex != constants::missing::uintValue) // && lastIndex != 0 && lastIndex != nodeIds.size() - 1) + { + + if (!subSequence.empty()) + { + std::cout << "adding: " << -999 << " " << -999 << ", " << -999 << std::endl; + subSequence.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); + } + + for (size_t i = lastIndex; i < nodeIds.size(); ++i) + { + std::cout << "adding: " << nodeIds[i] << " " << meshBoundaryPolygon[i].x << ", " << meshBoundaryPolygon[i].y << " " + << m_nodes[nodeIds[i]].x << ", " << m_nodes[nodeIds[i]].y + << std::endl; + // subSequence.emplace_back(m_nodes[nodeIds[i]]); + subSequence.emplace_back(meshBoundaryPolygon[i]); + } + + std::cout << "adding: " << nodeIds[lastIndex] << " " << meshBoundaryPolygon[lastIndex].x << ", " << meshBoundaryPolygon[lastIndex].y << std::endl; + std::cout << "last index " << lastIndex << " " << nodeIds[lastIndex] << " " << nodeIds.size() << std::endl; + std::cout << std::endl; + + // subSequence.emplace_back(m_nodes[nodeIds[lastIndex]]); + subSequence.emplace_back(meshBoundaryPolygon[lastIndex]); + meshBoundaryPolygon.resize(lastIndex); + nodeIds.resize(lastIndex); + } + } + + currentNode = nextNode; + meshBoundaryPolygon.emplace_back(m_nodes[currentNode]); + edgeIsVisited[currentEdge] = true; + nodeIsVisited[currentNode] = true; + nodeIds.emplace_back(currentNode); + } + + for (size_t i = 0; i < subSequence.size(); ++i) + { + std::cout << "subseq: " << subSequence[i].x << ", " << subSequence[i].y << std::endl; + } +} + std::vector Mesh2D::GetHangingEdges() const { std::vector result; @@ -1919,6 +2154,11 @@ std::unique_ptr Mesh2D::DeleteMeshFaces2(const Polygons& m_edgesFaces[edge][1] = constants::missing::uintValue; --m_edgesNumFaces[edge]; } + + if (m_edgesNumFaces[edge] == 0) + { + deleteMeshAction->Add(DeleteEdge(edge)); + } } m_facesNodes[f].clear(); @@ -1978,6 +2218,8 @@ std::unique_ptr Mesh2D::DeleteMeshFaces2(const Polygons& // } // } + // SetAdministrationRequired(true); + return deleteMeshAction; } diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index 983539c09d..9cf6d0f16f 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -807,6 +807,11 @@ std::unique_ptr MeshRefinement::RefineFacesBySplittingEd const auto numEdges = m_mesh.GetNumFaceEdges(f); + if (numEdges == 0) + { + continue; + } + m_mesh.ComputeFaceClosedPolygonWithLocalMappings(f, m_polygonNodesCache, m_localNodeIndicesCache, m_globalEdgeIndicesCache); UInt numBrotherEdges = 0; diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index be45ab2d44..59194894aa 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -2680,6 +2680,185 @@ TEST(MeshRefinement, RowSplittingFailureTests) EXPECT_THROW([[maybe_unused]] auto undo4 = splitMeshRow.Compute(mesh, edgeId), ConstraintError); } +TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygons) +{ + auto curviMesh = MakeCurvilinearGrid(0.0, 0.0, 20.0, 20.0, 11, 11); + const auto edges = curviMesh->ComputeEdges(); + const auto nodes = curviMesh->ComputeNodes(); + + Mesh2D mesh(edges, nodes, Projection::cartesian); + + const std::vector originalNodes(mesh.Nodes()); + const std::vector originalEdges(mesh.Edges()); + + std::vector patch{{45.0, 5.0}, + {155.0, 5.0}, + {155.0, 155.0}, + {45.0, 155.0}, + {45.0, 5.0}, + {constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + {65.0, 65.0}, + {115.0, 65.0}, + {115.0, 115.0}, + {65.0, 115.0}, + {65.0, 65.0}}; + + Polygons polygon(patch, Projection::cartesian); + + auto casulliUndoAction = CasulliRefinement::Compute(mesh, polygon); + + // Find nodes of elements to be deleted. + auto node11 = mesh.FindNodeCloseToAPoint({80.0, 0.0}, 1.0e-5); + auto node12 = mesh.FindNodeCloseToAPoint({85.0, 15.0}, 1.0e-5); + auto node13 = mesh.FindNodeCloseToAPoint({75.0, 15.0}, 1.0e-5); + + auto node21 = mesh.FindNodeCloseToAPoint({180.0, 140.0}, 1.0e-5); + auto node22 = mesh.FindNodeCloseToAPoint({200.0, 140.0}, 1.0e-5); + auto node23 = mesh.FindNodeCloseToAPoint({200.0, 160.0}, 1.0e-5); + auto node24 = mesh.FindNodeCloseToAPoint({180.0, 160.0}, 1.0e-5); + + auto node31 = mesh.FindNodeCloseToAPoint({125.0, 125.0}, 1.0e-5); + auto node32 = mesh.FindNodeCloseToAPoint({135.0, 125.0}, 1.0e-5); + auto node33 = mesh.FindNodeCloseToAPoint({135.0, 135.0}, 1.0e-5); + auto node34 = mesh.FindNodeCloseToAPoint({125.0, 135.0}, 1.0e-5); + + auto node41 = mesh.FindNodeCloseToAPoint({85.0, 15.0}, 1.0e-5); + auto node42 = mesh.FindNodeCloseToAPoint({95.0, 15.0}, 1.0e-5); + auto node43 = mesh.FindNodeCloseToAPoint({95.0, 25.0}, 1.0e-5); + auto node44 = mesh.FindNodeCloseToAPoint({85.0, 25.0}, 1.0e-5); + + auto node51 = mesh.FindNodeCloseToAPoint({100.0, 0.0}, 1.0e-5); + auto node52 = mesh.FindNodeCloseToAPoint({105.0, 15.0}, 1.0e-5); + auto node53 = mesh.FindNodeCloseToAPoint({95.0, 15.0}, 1.0e-5); + + auto node61 = mesh.FindNodeCloseToAPoint({120.0, 0.0}, 1.0e-5); + auto node62 = mesh.FindNodeCloseToAPoint({125.0, 15.0}, 1.0e-5); + auto node63 = mesh.FindNodeCloseToAPoint({115.0, 15.0}, 1.0e-5); + + std::vector polygonNodes; + auto boundaryNodes = mesh.ComputeBoundaryPolygons(polygonNodes); + + std::vector elementNodes1{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node11), + mesh.Node(node12), + mesh.Node(node13), + mesh.Node(node11)}; + + std::vector elementNodes2{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node21), + mesh.Node(node22), + mesh.Node(node23), + mesh.Node(node24), + mesh.Node(node21)}; + + std::vector elementNodes3{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node31), + mesh.Node(node32), + mesh.Node(node33), + mesh.Node(node34), + mesh.Node(node31)}; + + std::vector elementNodes4{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node41), + mesh.Node(node42), + mesh.Node(node43), + mesh.Node(node44), + mesh.Node(node41)}; + + std::vector elementNodes5{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node51), + mesh.Node(node52), + mesh.Node(node53), + mesh.Node(node51)}; + + std::vector elementNodes6{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node61), + mesh.Node(node62), + mesh.Node(node63), + mesh.Node(node61)}; + + // Combine all nodes to for a sequence of polygons + boundaryNodes.insert(boundaryNodes.end(), elementNodes1.begin(), elementNodes1.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes2.begin(), elementNodes2.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes3.begin(), elementNodes3.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes4.begin(), elementNodes4.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes5.begin(), elementNodes5.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes6.begin(), elementNodes6.end()); + + Polygons boundaryWithMissingElement(boundaryNodes, Projection::cartesian); + + auto deleteMeshFacesUndoAction = mesh.DeleteMeshFaces2(boundaryWithMissingElement, true); + + // Compute interior boundary polygon points + std::vector boundaryNodes2 = mesh.ComputeInnerBoundaryPolygons(); + + UInt expectedNumberOfNodes = 26; + + ASSERT_EQ(expectedNumberOfNodes, boundaryNodes2.size()); + + std::vector expectedXPoints{ + 95.0, + 100.0, + 105.0, + 95.0, + constants::missing::doubleValue, + 85.0, + 95.0, + 95.0, + 85.0, + 85.0, + constants::missing::doubleValue, + 80.0, + 85.0, + 75.0, + 80.0, + constants::missing::doubleValue, + 120.0, + 125.0, + 115.0, + 120.0, + constants::missing::doubleValue, + 125.0, + 125.0, + 135.0, + 135.0, + 125.0}; + + std::vector expectedYPoints{ + 15.0, + 0.0, + 15.0, + 15.0, + constants::missing::doubleValue, + 15.0, + 15.0, + 25.0, + 25.0, + 15.0, + constants::missing::doubleValue, + 0.0, + 15.0, + 15.0, + 0.0, + constants::missing::doubleValue, + 0.0, + 15.0, + 15.0, + 0.0, + constants::missing::doubleValue, + 125.0, + 135.0, + 135.0, + 125.0, + 125.0}; + + for (size_t i = 0; i < boundaryNodes2.size(); ++i) + { + EXPECT_EQ(expectedXPoints[i], boundaryNodes2[i].x); + EXPECT_EQ(expectedYPoints[i], boundaryNodes2[i].y); + } +} + TEST(MeshRefinement, WTF) { const size_t ExpectedNumberOfPoints = 184; @@ -2717,6 +2896,7 @@ TEST(MeshRefinement, WTF) auto node11 = mesh.FindNodeCloseToAPoint({80.0, 0.0}, 1.0e-5); auto node12 = mesh.FindNodeCloseToAPoint({85.0, 15.0}, 1.0e-5); auto node13 = mesh.FindNodeCloseToAPoint({75.0, 15.0}, 1.0e-5); + std::cout << "nodes1: " << node11 << " " << node12 << " " << node13 << std::endl; auto node21 = mesh.FindNodeCloseToAPoint({180.0, 140.0}, 1.0e-5); auto node22 = mesh.FindNodeCloseToAPoint({200.0, 140.0}, 1.0e-5); @@ -2732,6 +2912,27 @@ TEST(MeshRefinement, WTF) std::cout << "nodes3: " << node31 << " " << node32 << " " << node33 << " " << node34 << std::endl; + auto node41 = mesh.FindNodeCloseToAPoint({85.0, 15.0}, 1.0e-5); + auto node42 = mesh.FindNodeCloseToAPoint({95.0, 15.0}, 1.0e-5); + auto node43 = mesh.FindNodeCloseToAPoint({95.0, 25.0}, 1.0e-5); + auto node44 = mesh.FindNodeCloseToAPoint({85.0, 25.0}, 1.0e-5); + + std::cout << "nodes4: " << node41 << " " << node42 << " " << node43 << " " << node44 << std::endl; + + auto node51 = mesh.FindNodeCloseToAPoint({100.0, 0.0}, 1.0e-5); + auto node52 = mesh.FindNodeCloseToAPoint({105.0, 15.0}, 1.0e-5); + auto node53 = mesh.FindNodeCloseToAPoint({95.0, 15.0}, 1.0e-5); + std::cout << "nodes5: " << node51 << " " << node52 << " " << node53 << std::endl; + + auto node61 = mesh.FindNodeCloseToAPoint({120.0, 0.0}, 1.0e-5); + auto node62 = mesh.FindNodeCloseToAPoint({125.0, 15.0}, 1.0e-5); + auto node63 = mesh.FindNodeCloseToAPoint({115.0, 15.0}, 1.0e-5); + // auto node61 = mesh.FindNodeCloseToAPoint({140.0, 0.0}, 1.0e-5); + // auto node62 = mesh.FindNodeCloseToAPoint({145.0, 15.0}, 1.0e-5); + // auto node63 = mesh.FindNodeCloseToAPoint({135.0, 15.0}, 1.0e-5); + std::cout << "nodes6: " << node61 << " " << node62 << " " << node63 << std::endl; + + // mesh.Administrate(); std::vector polygonNodes; auto boundaryNodes = mesh.ComputeBoundaryPolygons(polygonNodes); @@ -2757,9 +2958,31 @@ TEST(MeshRefinement, WTF) mesh.Node(node34), mesh.Node(node31)}; - boundaryNodes.insert(boundaryNodes.end(), elementNodes1.begin(), elementNodes1.end()); - boundaryNodes.insert(boundaryNodes.end(), elementNodes2.begin(), elementNodes2.end()); - boundaryNodes.insert(boundaryNodes.end(), elementNodes3.begin(), elementNodes3.end()); + std::vector elementNodes4{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node41), + mesh.Node(node42), + mesh.Node(node43), + mesh.Node(node44), + mesh.Node(node41)}; + + std::vector elementNodes5{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node51), + mesh.Node(node52), + mesh.Node(node53), + mesh.Node(node51)}; + + std::vector elementNodes6{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node61), + mesh.Node(node62), + mesh.Node(node63), + mesh.Node(node61)}; + + // boundaryNodes.insert(boundaryNodes.end(), elementNodes1.begin(), elementNodes1.end()); + // boundaryNodes.insert(boundaryNodes.end(), elementNodes2.begin(), elementNodes2.end()); + // boundaryNodes.insert(boundaryNodes.end(), elementNodes3.begin(), elementNodes3.end()); + // boundaryNodes.insert(boundaryNodes.end(), elementNodes4.begin(), elementNodes4.end()); + // boundaryNodes.insert(boundaryNodes.end(), elementNodes5.begin(), elementNodes5.end()); + // boundaryNodes.insert(boundaryNodes.end(), elementNodes6.begin(), elementNodes6.end()); [[maybe_unused]] Polygons boundaryWithMissingElement(boundaryNodes, Projection::cartesian); @@ -2769,13 +2992,112 @@ TEST(MeshRefinement, WTF) // auto undoAction3 = ::Compute(mesh, Polygons()); + //-------------------------------- + + // const auto projectToLandBoundaryOption = LandBoundaries::ProjectToLandBoundaryOption::DoNotProjectToLandBoundary; + // OrthogonalizationParameters orthogonalizationParameters; + // orthogonalizationParameters.inner_iterations = 2; + // orthogonalizationParameters.boundary_iterations = 25; + // orthogonalizationParameters.outer_iterations = 25; + // orthogonalizationParameters.orthogonalization_to_smoothing_factor = 0.975; + // orthogonalizationParameters.orthogonalization_to_smoothing_factor_at_boundary = 0.975; + // orthogonalizationParameters.areal_to_angle_smoothing_factor = 1.0; + + // // Execute + // auto orthoPolygon = std::make_unique(); + + // std::vector landBoundary{}; + // auto landboundaries = std::make_unique(landBoundary, mesh, *orthoPolygon); + + // OrthogonalizationAndSmoothing orthogonalization(mesh, + // std::move(orthoPolygon), + // std::move(landboundaries), + // projectToLandBoundaryOption, + // orthogonalizationParameters); + + // [[maybe_unused]] auto undoActionOrtho = orthogonalization.Initialize(); + // orthogonalization.Compute(); + + //-------------------------------- + + auto undoAction2 = mesh.DeleteMeshFaces2(boundaryWithMissingElement, true); + + std::vector boundaryNodes2 = mesh.ComputeInnerBoundaryPolygons(); + + // UInt polygonLength = 0; + // UInt polygonStart = 0; + + // UInt polygonLength2 = constants::missing::uintValue; + // UInt polygonStart2 = 0; + + // UInt index = 0; + + // while (index < boundaryNodes2.size()) + // { + // polygonStart2 = index; + + // for (UInt i = polygonStart2; i < boundaryNodes2.size(); ++i) + // { + // ++index; + + // if (!boundaryNodes2[i].IsValid()) + // { + // polygonLength2 = i - polygonStart2; + // break; + // } + // } + + // if (polygonLength2 > polygonLength) + // { + // polygonLength = polygonLength2; + // polygonStart = polygonStart2; + // } + // } + + // if (polygonLength < boundaryNodes2.size()) + // { + // ++polygonLength; + // } + + // std::cout << " max polygon length " << polygonLength << " " << polygonStart << std::endl; + + // std::vector boundaryNodes3; + + // if (polygonStart > 0) + // { + // boundaryNodes3.insert(boundaryNodes3.begin(), boundaryNodes2.begin(), boundaryNodes2.begin() + polygonStart - 1); + // } + + // if (polygonStart + polygonLength < boundaryNodes2.size()) + // { + // boundaryNodes3.insert(boundaryNodes3.end(), boundaryNodes2.begin() + polygonStart + polygonLength, boundaryNodes2.end()); + // } + + int count = 0; + + for (size_t i = 0; i < boundaryNodes2.size(); ++i) + { + std::cout << "point " << count << " " << " = " << boundaryNodes2[i].x << ", " << boundaryNodes2[i].y << std::endl; + ++count; + } + + [[maybe_unused]] Polygons boundaryWithMissingElement2(boundaryNodes2, Projection::cartesian); + + // for (size_t i = 0; i < boundaryNodes3.size(); ++i) + // { + // std::cout << "point " << count << " " << " = " << boundaryNodes3[i].x << ", " << boundaryNodes3[i].y << std::endl; + // ++count; + // } + + //-------------------------------- + MeshRefinementParameters refParams; refParams.max_num_refinement_iterations = 1; MeshRefinement meshRefinement(mesh, Polygons(), refParams); - auto undoAction3 = meshRefinement.Compute(); + auto undoAction5 = meshRefinement.Compute(); - //-------------------------------- + mesh.Administrate(); const auto projectToLandBoundaryOption = LandBoundaries::ProjectToLandBoundaryOption::DoNotProjectToLandBoundary; OrthogonalizationParameters orthogonalizationParameters; @@ -2801,9 +3123,7 @@ TEST(MeshRefinement, WTF) [[maybe_unused]] auto undoActionOrtho = orthogonalization.Initialize(); orthogonalization.Compute(); - //-------------------------------- - - auto undoAction2 = mesh.DeleteMeshFaces2(boundaryWithMissingElement, true); + [[maybe_unused]] auto undoAction3 = mesh.DeleteMeshFaces2(boundaryWithMissingElement2, false); meshkernel::SaveVtk(mesh.Nodes(), mesh.m_facesNodes, "/home/wcs1/MeshKernel/MeshKernel04/build_deb/mesh1.vtu"); // meshkernel::Print(mesh.Nodes(), mesh.Edges()); From 95ddbafd7a0015fa4fd280b6a405b18d8a57fdaf Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 6 Oct 2025 14:51:17 +0200 Subject: [PATCH 04/16] GRIDEDIT-1795 Corrected orientation of all polygons, orientation should be ACW --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 21 +- libs/MeshKernel/src/Mesh2D.cpp | 215 ++++++-------- .../tests/src/MeshRefinementTests.cpp | 277 +----------------- 3 files changed, 110 insertions(+), 403 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index e3ed93a8b2..981f541054 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -373,7 +373,10 @@ namespace meshkernel /// If no such element can be found then the null value will be returned. UInt FindCommonFace(const UInt edge1, const UInt edge2) const; - [[nodiscard]] std::unique_ptr DeleteMeshFaces2(const Polygons& polygon, bool invertDeletion); + /// @brief Deletes the mesh faces inside a set of polygons + /// @param[in] polygon The polygon where to perform the operation + /// If this Polygons instance contains multiple polygons, the first one will be taken. + std::unique_ptr DeleteMeshFacesInPolygons(const Polygons& polygon); private: // orthogonalization @@ -440,13 +443,15 @@ namespace meshkernel /// @brief Constructs a polygon or polygons from the meshboundary, by walking through the mesh /// - /// If there are multiple polygons connected by a single node, then these will be separated into indicidual polygons - void WalkBoundaryFromNode(std::vector& edgeIsVisited, - std::vector& nodeIsVisited, - UInt& currentNode, - std::vector& meshBoundaryPolygon, - std::vector& nodeIds, - std::vector& subsSequence) const; + /// If there are multiple polygons connected by a single node, then these will be separated into individual polygons + void WalkMultiBoundaryFromNode(std::vector& edgeIsVisited, + std::vector& nodeIsVisited, + UInt& currentNode, + std::vector& meshBoundaryPolygon, + std::vector& nodeIds, + std::vector& subsSequence) const; + + void OrientatePolygonsAntiClockwise(std::vector& polygonNodes) const; /// @brief Removes the outer domain boundary polygon from the set of polygons /// diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 07b1d300c0..2098e60e50 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -25,11 +25,13 @@ // //------------------------------------------------------------------------------ -#include "MeshKernel/Mesh2D.hpp" +#include + #include "MeshKernel/Constants.hpp" #include "MeshKernel/Definitions.hpp" #include "MeshKernel/Entities.hpp" #include "MeshKernel/Exceptions.hpp" +#include "MeshKernel/Mesh2D.hpp" #include "MeshKernel/Mesh2DIntersections.hpp" #include "MeshKernel/MeshFaceCenters.hpp" #include "MeshKernel/MeshOrthogonality.hpp" @@ -1485,7 +1487,7 @@ std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector const Polygon polygon(polygonNodes, m_projection); // Find faces - // Administrate(); + Administrate(); std::vector isVisited(GetNumEdges(), false); std::vector meshBoundaryPolygon; meshBoundaryPolygon.reserve(GetNumNodes()); @@ -1513,7 +1515,6 @@ std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector // Start a new polyline if (!meshBoundaryPolygon.empty()) { - // meshBoundaryPolygon.emplace_back(constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator); meshBoundaryPolygon.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); } @@ -1598,7 +1599,7 @@ std::vector Mesh2D::ComputeInnerBoundaryPolygons() const // walk the current mesh boundary auto currentNode = secondNodeIndex; - WalkBoundaryFromNode(edgeIsVisited, nodeIsVisited, currentNode, subSqeuence, nodeIds, meshBoundaryPolygon); + WalkMultiBoundaryFromNode(edgeIsVisited, nodeIsVisited, currentNode, subSqeuence, nodeIds, meshBoundaryPolygon); const auto numNodesFirstTail = static_cast(subSqeuence.size()); @@ -1607,7 +1608,7 @@ std::vector Mesh2D::ComputeInnerBoundaryPolygons() const { // Now grow a polyline starting at the other side of the original link L, i.e., the second tail currentNode = firstNodeIndex; - WalkBoundaryFromNode(edgeIsVisited, nodeIsVisited, currentNode, subSqeuence, nodeIds, meshBoundaryPolygon); + WalkMultiBoundaryFromNode(edgeIsVisited, nodeIsVisited, currentNode, subSqeuence, nodeIds, meshBoundaryPolygon); } // There is a nonempty second tail, so reverse the first tail, so that they connect. @@ -1628,7 +1629,64 @@ std::vector Mesh2D::ComputeInnerBoundaryPolygons() const } } - return RemoveOuterDomainBoundaryPolygon(meshBoundaryPolygon); + std::vector innerPolygonNodes(RemoveOuterDomainBoundaryPolygon(meshBoundaryPolygon)); + OrientatePolygonsAntiClockwise(innerPolygonNodes); + return innerPolygonNodes; +} + +void Mesh2D::OrientatePolygonsAntiClockwise(std::vector& polygonNodes) const +{ + UInt polygonStart = 0; + UInt polygonLength = 0; + UInt index = 0; + + while (index < polygonNodes.size()) + { + polygonStart = index; + + for (size_t i = polygonStart; i < polygonNodes.size(); ++i) + { + ++index; + + if (!polygonNodes[i].IsValid()) + { + polygonLength = i - polygonStart; + break; + } + + if (index == polygonNodes.size()) + { + ++polygonLength; + } + } + + if (polygonLength > 0) + { + const Point inValidPoint = {constants::missing::doubleValue, constants::missing::doubleValue}; + Point zeroPoint{0.0, 0.0}; + + Point midPoint = std::accumulate(polygonNodes.begin() + polygonStart, polygonNodes.begin() + polygonStart + polygonLength - 1, zeroPoint) / static_cast(polygonLength - 1); + + if (!IsPointInPolygonNodes(midPoint, polygonNodes, m_projection, inValidPoint, polygonStart, polygonStart + polygonLength - 1)) + { + // reverse order of polygon nodes + if (polygonLength - 1 == 3) + { + // Only the second and third points need be swapped to reverse the points in the polygon + std::swap(polygonNodes[polygonStart + 1], polygonNodes[polygonStart + 2]); + } + else if (polygonLength - 1 == 4) + { + // Only the second and fourth points need be swapped to reverse the points in the polygon + std::swap(polygonNodes[polygonStart + 1], polygonNodes[polygonStart + 3]); + } + else + { + std::reverse(polygonNodes.begin() + polygonStart, polygonNodes.begin() + polygonStart + polygonLength); + } + } + } + } } std::vector Mesh2D::RemoveOuterDomainBoundaryPolygon(const std::vector& polygonNodes) const @@ -1636,52 +1694,53 @@ std::vector Mesh2D::RemoveOuterDomainBoundaryPolygon(const st // Remove outer boundary. // It is assumed that the boundary polygon with the most points is from the outer domain boundary - UInt polygonLength = 0; - UInt polygonStart = 0; + UInt outerPolygonLength = 0; + UInt outerPolygonStart = 0; - UInt polygonLengthIntermediate = 0; - UInt polygonStartIntermediate = 0; + UInt outerPolygonLengthIntermediate = 0; + UInt outerPolygonStartIntermediate = 0; UInt index = 0; // Find start index of outer boundary and the length of the polygon while (index < polygonNodes.size()) { - polygonStartIntermediate = index; + outerPolygonStartIntermediate = index; - for (UInt i = polygonStartIntermediate; i < polygonNodes.size(); ++i) + for (UInt i = outerPolygonStartIntermediate; i < polygonNodes.size(); ++i) { ++index; if (!polygonNodes[i].IsValid()) { - polygonLengthIntermediate = i - polygonStartIntermediate; + outerPolygonLengthIntermediate = i - outerPolygonStartIntermediate; break; } } - if (polygonLengthIntermediate > polygonLength) + if (outerPolygonLengthIntermediate > outerPolygonLength) { - polygonLength = polygonLengthIntermediate; - polygonStart = polygonStartIntermediate; + outerPolygonLength = outerPolygonLengthIntermediate; + outerPolygonStart = outerPolygonStartIntermediate; } } - if (polygonLength > 0 && polygonLength < polygonNodes.size()) + if (outerPolygonLength > 0 && outerPolygonLength < polygonNodes.size()) { - ++polygonLength; + ++outerPolygonLength; } std::vector innerBoundaryNodes; - if (polygonStart > 0) + // Gather all node except those from the outer domain boundary polygon (this is assumed to be the polygon with the largest number of nodes) + if (outerPolygonStart > 0) { - innerBoundaryNodes.insert(innerBoundaryNodes.begin(), polygonNodes.begin(), polygonNodes.begin() + polygonStart - 1); + innerBoundaryNodes.insert(innerBoundaryNodes.begin(), polygonNodes.begin(), polygonNodes.begin() + outerPolygonStart - 1); } - if (polygonLength > 0 && polygonStart + polygonLength < polygonNodes.size()) + if (outerPolygonLength > 0 && outerPolygonStart + outerPolygonLength < polygonNodes.size()) { - innerBoundaryNodes.insert(innerBoundaryNodes.end(), polygonNodes.begin() + polygonStart + polygonLength, polygonNodes.end()); + innerBoundaryNodes.insert(innerBoundaryNodes.end(), polygonNodes.begin() + outerPolygonStart + outerPolygonLength, polygonNodes.end()); } return innerBoundaryNodes; @@ -1722,21 +1781,15 @@ void Mesh2D::WalkBoundaryFromNode(const Polygon& polygon, } } -void Mesh2D::WalkBoundaryFromNode(std::vector& edgeIsVisited, - std::vector& nodeIsVisited, - UInt& currentNode, - std::vector& meshBoundaryPolygon, - std::vector& nodeIds, - std::vector& subSequence) const -{ - std::cout << std::endl - << " Mesh2D::WalkBoundaryFromNode " << currentNode << " " << m_nodes[currentNode].x << ", " << m_nodes[currentNode].y - << " " << subSequence.size() - << std::endl; +void Mesh2D::WalkMultiBoundaryFromNode(std::vector& edgeIsVisited, + std::vector& nodeIsVisited, + UInt& currentNode, + std::vector& meshBoundaryPolygon, + std::vector& nodeIds, + std::vector& subSequence) const +{ UInt e = 0; - // size_t start = meshBoundaryPolygon.size(); - while (e < m_nodesNumEdges[currentNode]) { const auto currentEdge = m_nodesEdges[currentNode][e]; @@ -1771,38 +1824,19 @@ void Mesh2D::WalkBoundaryFromNode(std::vector& edgeIsVisited, } } - std::cout << "node ids: "; - - for (size_t i = lastIndex; i < nodeIds.size(); ++i) - { - std::cout << nodeIds[i] << " "; - } - - std::cout << std::endl; - - if (lastIndex != constants::missing::uintValue) // && lastIndex != 0 && lastIndex != nodeIds.size() - 1) + if (lastIndex != constants::missing::uintValue) { if (!subSequence.empty()) { - std::cout << "adding: " << -999 << " " << -999 << ", " << -999 << std::endl; subSequence.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); } for (size_t i = lastIndex; i < nodeIds.size(); ++i) { - std::cout << "adding: " << nodeIds[i] << " " << meshBoundaryPolygon[i].x << ", " << meshBoundaryPolygon[i].y << " " - << m_nodes[nodeIds[i]].x << ", " << m_nodes[nodeIds[i]].y - << std::endl; - // subSequence.emplace_back(m_nodes[nodeIds[i]]); subSequence.emplace_back(meshBoundaryPolygon[i]); } - std::cout << "adding: " << nodeIds[lastIndex] << " " << meshBoundaryPolygon[lastIndex].x << ", " << meshBoundaryPolygon[lastIndex].y << std::endl; - std::cout << "last index " << lastIndex << " " << nodeIds[lastIndex] << " " << nodeIds.size() << std::endl; - std::cout << std::endl; - - // subSequence.emplace_back(m_nodes[nodeIds[lastIndex]]); subSequence.emplace_back(meshBoundaryPolygon[lastIndex]); meshBoundaryPolygon.resize(lastIndex); nodeIds.resize(lastIndex); @@ -1815,11 +1849,6 @@ void Mesh2D::WalkBoundaryFromNode(std::vector& edgeIsVisited, nodeIsVisited[currentNode] = true; nodeIds.emplace_back(currentNode); } - - for (size_t i = 0; i < subSequence.size(); ++i) - { - std::cout << "subseq: " << subSequence[i].x << ", " << subSequence[i].y << std::endl; - } } std::vector Mesh2D::GetHangingEdges() const @@ -2105,13 +2134,11 @@ std::unique_ptr Mesh2D::DeleteMeshFaces(const Polygons& return deleteMeshAction; } -std::unique_ptr Mesh2D::DeleteMeshFaces2(const Polygons& polygon, bool invertDeletion) +std::unique_ptr Mesh2D::DeleteMeshFacesInPolygons(const Polygons& polygon) { std::unique_ptr deleteMeshAction = CompoundUndoAction::Create(); Administrate(deleteMeshAction.get()); - // std::vector faceCircumcenters = algo::ComputeFaceCircumcenters(*this); - std::vector faceCircumcenters = m_facesMassCenters; bool isInPolygon; UInt polygonIndex; @@ -2123,19 +2150,12 @@ std::unique_ptr Mesh2D::DeleteMeshFaces2(const Polygons& // const auto faceIndex = m_edgesFaces[e][f]; - if (!faceCircumcenters[f].IsValid()) + if (!m_facesMassCenters[f].IsValid()) { continue; } - std::tie(isInPolygon, polygonIndex) = polygon.IsPointInPolygons(faceCircumcenters[f]); - - std::cout << "centre is in polygon: " << f << " " << std::boolalpha << isInPolygon << " " << polygonIndex << " " << faceCircumcenters[f].x << ", " << faceCircumcenters[f].y << std::endl; - - if (invertDeletion) - { - isInPolygon = !isInPolygon; - } + std::tie(isInPolygon, polygonIndex) = polygon.IsPointInPolygons(m_facesMassCenters[f]); if (isInPolygon) { @@ -2166,60 +2186,9 @@ std::unique_ptr Mesh2D::DeleteMeshFaces2(const Polygons& m_facesEdges[f].clear(); m_facesMassCenters[f] = {constants::missing::doubleValue, constants::missing::doubleValue}; m_faceArea[f] = constants::missing::doubleValue; - - // deleteMeshAction->Add(DeleteEdge(e)); } } - // for (UInt e = 0u; e < GetNumEdges(); ++e) - // { - // for (UInt f = 0u; f < GetNumEdgesFaces(e); ++f) - // { - // isInPolygon = false; - // polygonIndex = constants::missing::uintValue; - - // const auto faceIndex = m_edgesFaces[e][f]; - // if (faceIndex == constants::missing::uintValue) - // { - // continue; - // } - - // std::tie(isInPolygon, polygonIndex) = polygon.IsPointInPolygons(faceCircumcenters[faceIndex]); - - // std::cout << "centre is in polygon: " << std::boolalpha << isInPolygon << " " << polygonIndex << " " << faceCircumcenters[faceIndex].x << ", " << faceCircumcenters[faceIndex].y << std::endl; - - // if (invertDeletion) - // { - // isInPolygon = !isInPolygon; - // } - - // if (isInPolygon) - // { - - // if (m_edgesFaces[e][0] == f) - // { - // m_edgesFaces[e][0] = constants::missing::uintValue; - // --m_edgesNumFaces[e]; - // } - // else if (m_edgesFaces[e][1] == f) - // { - // m_edgesFaces[e][1] = constants::missing::uintValue; - // --m_edgesNumFaces[e]; - // } - - // m_facesNodes[f].clear(); - // m_numFacesNodes[f] = 0; - // m_facesEdges[f].clear(); - // m_facesMassCenters[f] = {constants::missing::doubleValue, constants::missing::doubleValue}; - // m_faceArea[f] = constants::missing::doubleValue; - - // // deleteMeshAction->Add(DeleteEdge(e)); - // } - // } - // } - - // SetAdministrationRequired(true); - return deleteMeshAction; } diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index 59194894aa..2b8f0d0190 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -2680,8 +2680,9 @@ TEST(MeshRefinement, RowSplittingFailureTests) EXPECT_THROW([[maybe_unused]] auto undo4 = splitMeshRow.Compute(mesh, edgeId), ConstraintError); } -TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygons) +TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygonsForSixFaces) { + auto curviMesh = MakeCurvilinearGrid(0.0, 0.0, 20.0, 20.0, 11, 11); const auto edges = curviMesh->ComputeEdges(); const auto nodes = curviMesh->ComputeNodes(); @@ -2787,7 +2788,7 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygons) Polygons boundaryWithMissingElement(boundaryNodes, Projection::cartesian); - auto deleteMeshFacesUndoAction = mesh.DeleteMeshFaces2(boundaryWithMissingElement, true); + auto deleteMeshFacesUndoAction = mesh.DeleteMeshFacesInPolygons(boundaryWithMissingElement); // Compute interior boundary polygon points std::vector boundaryNodes2 = mesh.ComputeInnerBoundaryPolygons(); @@ -2796,6 +2797,8 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygons) ASSERT_EQ(expectedNumberOfNodes, boundaryNodes2.size()); + // The edge of one of the deleted elements lies on the boundary, so will be not be part of the + // interior set of polygons std::vector expectedXPoints{ 95.0, 100.0, @@ -2858,273 +2861,3 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygons) EXPECT_EQ(expectedYPoints[i], boundaryNodes2[i].y); } } - -TEST(MeshRefinement, WTF) -{ - const size_t ExpectedNumberOfPoints = 184; - const size_t ExpectedNumberOfEdges = 360; - - auto curviMesh = MakeCurvilinearGrid(0.0, 0.0, 20.0, 20.0, 11, 11); - const auto edges = curviMesh->ComputeEdges(); - const auto nodes = curviMesh->ComputeNodes(); - - Mesh2D mesh(edges, nodes, Projection::cartesian); - - const std::vector originalNodes(mesh.Nodes()); - const std::vector originalEdges(mesh.Edges()); - - std::vector patch{{45.0, 5.0}, - {155.0, 5.0}, - {155.0, 155.0}, - {45.0, 155.0}, - {45.0, 5.0}, - {constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, - {65.0, 65.0}, - {115.0, 65.0}, - {115.0, 115.0}, - {65.0, 115.0}, - {65.0, 65.0}}; - - std::vector expectedPoints(ExpectedNumberOfPoints); - std::vector expectedEdgeStart(ExpectedNumberOfEdges); - std::vector expectedEdgeEnd(ExpectedNumberOfEdges); - - Polygons polygon(patch, Projection::cartesian); - - auto undoAction = CasulliRefinement::Compute(mesh, polygon); - - auto node11 = mesh.FindNodeCloseToAPoint({80.0, 0.0}, 1.0e-5); - auto node12 = mesh.FindNodeCloseToAPoint({85.0, 15.0}, 1.0e-5); - auto node13 = mesh.FindNodeCloseToAPoint({75.0, 15.0}, 1.0e-5); - std::cout << "nodes1: " << node11 << " " << node12 << " " << node13 << std::endl; - - auto node21 = mesh.FindNodeCloseToAPoint({180.0, 140.0}, 1.0e-5); - auto node22 = mesh.FindNodeCloseToAPoint({200.0, 140.0}, 1.0e-5); - auto node23 = mesh.FindNodeCloseToAPoint({200.0, 160.0}, 1.0e-5); - auto node24 = mesh.FindNodeCloseToAPoint({180.0, 160.0}, 1.0e-5); - - std::cout << "nodes2: " << node21 << " " << node22 << " " << node23 << " " << node24 << std::endl; - - auto node31 = mesh.FindNodeCloseToAPoint({125.0, 125.0}, 1.0e-5); - auto node32 = mesh.FindNodeCloseToAPoint({135.0, 125.0}, 1.0e-5); - auto node33 = mesh.FindNodeCloseToAPoint({135.0, 135.0}, 1.0e-5); - auto node34 = mesh.FindNodeCloseToAPoint({125.0, 135.0}, 1.0e-5); - - std::cout << "nodes3: " << node31 << " " << node32 << " " << node33 << " " << node34 << std::endl; - - auto node41 = mesh.FindNodeCloseToAPoint({85.0, 15.0}, 1.0e-5); - auto node42 = mesh.FindNodeCloseToAPoint({95.0, 15.0}, 1.0e-5); - auto node43 = mesh.FindNodeCloseToAPoint({95.0, 25.0}, 1.0e-5); - auto node44 = mesh.FindNodeCloseToAPoint({85.0, 25.0}, 1.0e-5); - - std::cout << "nodes4: " << node41 << " " << node42 << " " << node43 << " " << node44 << std::endl; - - auto node51 = mesh.FindNodeCloseToAPoint({100.0, 0.0}, 1.0e-5); - auto node52 = mesh.FindNodeCloseToAPoint({105.0, 15.0}, 1.0e-5); - auto node53 = mesh.FindNodeCloseToAPoint({95.0, 15.0}, 1.0e-5); - std::cout << "nodes5: " << node51 << " " << node52 << " " << node53 << std::endl; - - auto node61 = mesh.FindNodeCloseToAPoint({120.0, 0.0}, 1.0e-5); - auto node62 = mesh.FindNodeCloseToAPoint({125.0, 15.0}, 1.0e-5); - auto node63 = mesh.FindNodeCloseToAPoint({115.0, 15.0}, 1.0e-5); - // auto node61 = mesh.FindNodeCloseToAPoint({140.0, 0.0}, 1.0e-5); - // auto node62 = mesh.FindNodeCloseToAPoint({145.0, 15.0}, 1.0e-5); - // auto node63 = mesh.FindNodeCloseToAPoint({135.0, 15.0}, 1.0e-5); - std::cout << "nodes6: " << node61 << " " << node62 << " " << node63 << std::endl; - - // mesh.Administrate(); - std::vector polygonNodes; - auto boundaryNodes = mesh.ComputeBoundaryPolygons(polygonNodes); - - // std::vector boundaryNodes = meshBoundaryPolygon.GatherAllEnclosureNodes(); - - std::vector elementNodes1{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, - mesh.Node(node11), - mesh.Node(node12), - mesh.Node(node13), - mesh.Node(node11)}; - - std::vector elementNodes2{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, - mesh.Node(node21), - mesh.Node(node22), - mesh.Node(node23), - mesh.Node(node24), - mesh.Node(node21)}; - - std::vector elementNodes3{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, - mesh.Node(node31), - mesh.Node(node32), - mesh.Node(node33), - mesh.Node(node34), - mesh.Node(node31)}; - - std::vector elementNodes4{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, - mesh.Node(node41), - mesh.Node(node42), - mesh.Node(node43), - mesh.Node(node44), - mesh.Node(node41)}; - - std::vector elementNodes5{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, - mesh.Node(node51), - mesh.Node(node52), - mesh.Node(node53), - mesh.Node(node51)}; - - std::vector elementNodes6{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, - mesh.Node(node61), - mesh.Node(node62), - mesh.Node(node63), - mesh.Node(node61)}; - - // boundaryNodes.insert(boundaryNodes.end(), elementNodes1.begin(), elementNodes1.end()); - // boundaryNodes.insert(boundaryNodes.end(), elementNodes2.begin(), elementNodes2.end()); - // boundaryNodes.insert(boundaryNodes.end(), elementNodes3.begin(), elementNodes3.end()); - // boundaryNodes.insert(boundaryNodes.end(), elementNodes4.begin(), elementNodes4.end()); - // boundaryNodes.insert(boundaryNodes.end(), elementNodes5.begin(), elementNodes5.end()); - // boundaryNodes.insert(boundaryNodes.end(), elementNodes6.begin(), elementNodes6.end()); - - [[maybe_unused]] Polygons boundaryWithMissingElement(boundaryNodes, Projection::cartesian); - - // auto edge1 = mesh.FindEdge(node1, node2); - // auto edge2 = mesh.FindEdge(node2, node3); - // auto edge3 = mesh.FindEdge(node3, node1); - - // auto undoAction3 = ::Compute(mesh, Polygons()); - - //-------------------------------- - - // const auto projectToLandBoundaryOption = LandBoundaries::ProjectToLandBoundaryOption::DoNotProjectToLandBoundary; - // OrthogonalizationParameters orthogonalizationParameters; - // orthogonalizationParameters.inner_iterations = 2; - // orthogonalizationParameters.boundary_iterations = 25; - // orthogonalizationParameters.outer_iterations = 25; - // orthogonalizationParameters.orthogonalization_to_smoothing_factor = 0.975; - // orthogonalizationParameters.orthogonalization_to_smoothing_factor_at_boundary = 0.975; - // orthogonalizationParameters.areal_to_angle_smoothing_factor = 1.0; - - // // Execute - // auto orthoPolygon = std::make_unique(); - - // std::vector landBoundary{}; - // auto landboundaries = std::make_unique(landBoundary, mesh, *orthoPolygon); - - // OrthogonalizationAndSmoothing orthogonalization(mesh, - // std::move(orthoPolygon), - // std::move(landboundaries), - // projectToLandBoundaryOption, - // orthogonalizationParameters); - - // [[maybe_unused]] auto undoActionOrtho = orthogonalization.Initialize(); - // orthogonalization.Compute(); - - //-------------------------------- - - auto undoAction2 = mesh.DeleteMeshFaces2(boundaryWithMissingElement, true); - - std::vector boundaryNodes2 = mesh.ComputeInnerBoundaryPolygons(); - - // UInt polygonLength = 0; - // UInt polygonStart = 0; - - // UInt polygonLength2 = constants::missing::uintValue; - // UInt polygonStart2 = 0; - - // UInt index = 0; - - // while (index < boundaryNodes2.size()) - // { - // polygonStart2 = index; - - // for (UInt i = polygonStart2; i < boundaryNodes2.size(); ++i) - // { - // ++index; - - // if (!boundaryNodes2[i].IsValid()) - // { - // polygonLength2 = i - polygonStart2; - // break; - // } - // } - - // if (polygonLength2 > polygonLength) - // { - // polygonLength = polygonLength2; - // polygonStart = polygonStart2; - // } - // } - - // if (polygonLength < boundaryNodes2.size()) - // { - // ++polygonLength; - // } - - // std::cout << " max polygon length " << polygonLength << " " << polygonStart << std::endl; - - // std::vector boundaryNodes3; - - // if (polygonStart > 0) - // { - // boundaryNodes3.insert(boundaryNodes3.begin(), boundaryNodes2.begin(), boundaryNodes2.begin() + polygonStart - 1); - // } - - // if (polygonStart + polygonLength < boundaryNodes2.size()) - // { - // boundaryNodes3.insert(boundaryNodes3.end(), boundaryNodes2.begin() + polygonStart + polygonLength, boundaryNodes2.end()); - // } - - int count = 0; - - for (size_t i = 0; i < boundaryNodes2.size(); ++i) - { - std::cout << "point " << count << " " << " = " << boundaryNodes2[i].x << ", " << boundaryNodes2[i].y << std::endl; - ++count; - } - - [[maybe_unused]] Polygons boundaryWithMissingElement2(boundaryNodes2, Projection::cartesian); - - // for (size_t i = 0; i < boundaryNodes3.size(); ++i) - // { - // std::cout << "point " << count << " " << " = " << boundaryNodes3[i].x << ", " << boundaryNodes3[i].y << std::endl; - // ++count; - // } - - //-------------------------------- - - MeshRefinementParameters refParams; - refParams.max_num_refinement_iterations = 1; - - MeshRefinement meshRefinement(mesh, Polygons(), refParams); - auto undoAction5 = meshRefinement.Compute(); - - mesh.Administrate(); - - const auto projectToLandBoundaryOption = LandBoundaries::ProjectToLandBoundaryOption::DoNotProjectToLandBoundary; - OrthogonalizationParameters orthogonalizationParameters; - orthogonalizationParameters.inner_iterations = 2; - orthogonalizationParameters.boundary_iterations = 25; - orthogonalizationParameters.outer_iterations = 25; - orthogonalizationParameters.orthogonalization_to_smoothing_factor = 0.975; - orthogonalizationParameters.orthogonalization_to_smoothing_factor_at_boundary = 0.975; - orthogonalizationParameters.areal_to_angle_smoothing_factor = 1.0; - - // Execute - auto orthoPolygon = std::make_unique(); - - std::vector landBoundary{}; - auto landboundaries = std::make_unique(landBoundary, mesh, *orthoPolygon); - - OrthogonalizationAndSmoothing orthogonalization(mesh, - std::move(orthoPolygon), - std::move(landboundaries), - projectToLandBoundaryOption, - orthogonalizationParameters); - - [[maybe_unused]] auto undoActionOrtho = orthogonalization.Initialize(); - orthogonalization.Compute(); - - [[maybe_unused]] auto undoAction3 = mesh.DeleteMeshFaces2(boundaryWithMissingElement2, false); - - meshkernel::SaveVtk(mesh.Nodes(), mesh.m_facesNodes, "/home/wcs1/MeshKernel/MeshKernel04/build_deb/mesh1.vtu"); - // meshkernel::Print(mesh.Nodes(), mesh.Edges()); -} From 65a5cb25fa26749d66280b9c6f2f7f1e656075d4 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 6 Oct 2025 18:11:27 +0200 Subject: [PATCH 05/16] GRIDEDIT-1795 Passing to other computer --- libs/MeshKernel/src/Mesh2D.cpp | 41 +- .../tests/src/MeshRefinementTests.cpp | 402 ++++++++++++++++++ 2 files changed, 440 insertions(+), 3 deletions(-) diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 2098e60e50..a40c903362 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -584,6 +584,12 @@ void Mesh2D::FindFacesGivenFaceNodesMapping(const std::vector> std::vector local_node_indices; for (UInt f = 0; f < m_facesNodes.size(); ++f) { + + if (numFaceNodes[f] == 0) + { + continue; + } + local_edges.clear(); local_nodes.clear(); local_node_indices.clear(); @@ -1556,6 +1562,10 @@ std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector std::vector Mesh2D::ComputeInnerBoundaryPolygons() const { + if (GetNumFaces() == 0) + { + return std::vector(); + } std::vector meshBoundaryPolygon; meshBoundaryPolygon.reserve(GetNumNodes()); @@ -1631,6 +1641,7 @@ std::vector Mesh2D::ComputeInnerBoundaryPolygons() const std::vector innerPolygonNodes(RemoveOuterDomainBoundaryPolygon(meshBoundaryPolygon)); OrientatePolygonsAntiClockwise(innerPolygonNodes); + return innerPolygonNodes; } @@ -2136,6 +2147,11 @@ std::unique_ptr Mesh2D::DeleteMeshFaces(const Polygons& std::unique_ptr Mesh2D::DeleteMeshFacesInPolygons(const Polygons& polygon) { + if (polygon.IsEmpty()) + { + return nullptr; + } + std::unique_ptr deleteMeshAction = CompoundUndoAction::Create(); Administrate(deleteMeshAction.get()); @@ -2143,13 +2159,16 @@ std::unique_ptr Mesh2D::DeleteMeshFacesInPolygons(const bool isInPolygon; UInt polygonIndex; + std::vector faceIndicesOffset(GetNumFaces(), 0); + std::vector faceIndices(GetNumFaces(), constants::missing::uintValue); + UInt faceIndexOffset = 0; + UInt faceIndex = 0; + for (UInt f = 0u; f < GetNumFaces(); ++f) { isInPolygon = false; polygonIndex = constants::missing::uintValue; - // const auto faceIndex = m_edgesFaces[e][f]; - if (!m_facesMassCenters[f].IsValid()) { continue; @@ -2157,7 +2176,8 @@ std::unique_ptr Mesh2D::DeleteMeshFacesInPolygons(const std::tie(isInPolygon, polygonIndex) = polygon.IsPointInPolygons(m_facesMassCenters[f]); - if (isInPolygon) + // TODO I think the polygon is orientated in the wrong direction + if (!isInPolygon) { for (UInt e = 0u; e < m_facesEdges[f].size(); ++e) @@ -2181,12 +2201,27 @@ std::unique_ptr Mesh2D::DeleteMeshFacesInPolygons(const } } + ++faceIndexOffset; m_facesNodes[f].clear(); m_numFacesNodes[f] = 0; m_facesEdges[f].clear(); m_facesMassCenters[f] = {constants::missing::doubleValue, constants::missing::doubleValue}; m_faceArea[f] = constants::missing::doubleValue; } + else + { + faceIndices[f] = faceIndex; + ++faceIndex; + } + + faceIndicesOffset[f] = faceIndexOffset; + } + + std::cout << "faceIndicesOffset" << std::endl; + + for (size_t i = 0; i < faceIndicesOffset.size(); ++i) + { + std::cout << " faceIndicesOffset " << i << " = " << faceIndicesOffset[i] << " " << faceIndices[i] << std::endl; } return deleteMeshAction; diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index 2b8f0d0190..1023a686b4 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -2789,6 +2789,7 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygonsForSixFa Polygons boundaryWithMissingElement(boundaryNodes, Projection::cartesian); auto deleteMeshFacesUndoAction = mesh.DeleteMeshFacesInPolygons(boundaryWithMissingElement); + meshkernel::SaveVtk(mesh.Nodes(), mesh.m_facesNodes, "/home/wcs1/MeshKernel/MeshKernel04/build_deb/mesh1.vtu"); // Compute interior boundary polygon points std::vector boundaryNodes2 = mesh.ComputeInnerBoundaryPolygons(); @@ -2861,3 +2862,404 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygonsForSixFa EXPECT_EQ(expectedYPoints[i], boundaryNodes2[i].y); } } + +TEST(MeshRefinement, WTF) +{ + const size_t ExpectedNumberOfPoints = 184; + const size_t ExpectedNumberOfEdges = 360; + + auto curviMesh = MakeCurvilinearGrid(0.0, 0.0, 20.0, 20.0, 11, 11); + const auto edges = curviMesh->ComputeEdges(); + const auto nodes = curviMesh->ComputeNodes(); + + Mesh2D mesh(edges, nodes, Projection::cartesian); + + const std::vector originalNodes(mesh.Nodes()); + const std::vector originalEdges(mesh.Edges()); + + std::vector patch{{45.0, 5.0}, + {155.0, 5.0}, + {155.0, 155.0}, + {45.0, 155.0}, + {45.0, 5.0}, + {constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + {65.0, 65.0}, + {115.0, 65.0}, + {115.0, 115.0}, + {65.0, 115.0}, + {65.0, 65.0}}; + + std::vector expectedPoints(ExpectedNumberOfPoints); + std::vector expectedEdgeStart(ExpectedNumberOfEdges); + std::vector expectedEdgeEnd(ExpectedNumberOfEdges); + + Polygons polygon(patch, Projection::cartesian); + + auto undoAction = CasulliRefinement::Compute(mesh, polygon); + + auto node11 = mesh.FindNodeCloseToAPoint({80.0, 0.0}, 1.0e-5); + auto node12 = mesh.FindNodeCloseToAPoint({85.0, 15.0}, 1.0e-5); + auto node13 = mesh.FindNodeCloseToAPoint({75.0, 15.0}, 1.0e-5); + std::cout << "nodes1: " << node11 << " " << node12 << " " << node13 << std::endl; + + auto node21 = mesh.FindNodeCloseToAPoint({180.0, 140.0}, 1.0e-5); + auto node22 = mesh.FindNodeCloseToAPoint({200.0, 140.0}, 1.0e-5); + auto node23 = mesh.FindNodeCloseToAPoint({200.0, 160.0}, 1.0e-5); + auto node24 = mesh.FindNodeCloseToAPoint({180.0, 160.0}, 1.0e-5); + + std::cout << "nodes2: " << node21 << " " << node22 << " " << node23 << " " << node24 << std::endl; + + auto node31 = mesh.FindNodeCloseToAPoint({125.0, 125.0}, 1.0e-5); + auto node32 = mesh.FindNodeCloseToAPoint({135.0, 125.0}, 1.0e-5); + auto node33 = mesh.FindNodeCloseToAPoint({135.0, 135.0}, 1.0e-5); + auto node34 = mesh.FindNodeCloseToAPoint({125.0, 135.0}, 1.0e-5); + + std::cout << "nodes3: " << node31 << " " << node32 << " " << node33 << " " << node34 << std::endl; + + auto node41 = mesh.FindNodeCloseToAPoint({85.0, 15.0}, 1.0e-5); + auto node42 = mesh.FindNodeCloseToAPoint({95.0, 15.0}, 1.0e-5); + auto node43 = mesh.FindNodeCloseToAPoint({95.0, 25.0}, 1.0e-5); + auto node44 = mesh.FindNodeCloseToAPoint({85.0, 25.0}, 1.0e-5); + + std::cout << "nodes4: " << node41 << " " << node42 << " " << node43 << " " << node44 << std::endl; + + auto node51 = mesh.FindNodeCloseToAPoint({100.0, 0.0}, 1.0e-5); + auto node52 = mesh.FindNodeCloseToAPoint({105.0, 15.0}, 1.0e-5); + auto node53 = mesh.FindNodeCloseToAPoint({95.0, 15.0}, 1.0e-5); + std::cout << "nodes5: " << node51 << " " << node52 << " " << node53 << std::endl; + + auto node61 = mesh.FindNodeCloseToAPoint({120.0, 0.0}, 1.0e-5); + auto node62 = mesh.FindNodeCloseToAPoint({125.0, 15.0}, 1.0e-5); + auto node63 = mesh.FindNodeCloseToAPoint({115.0, 15.0}, 1.0e-5); + // auto node61 = mesh.FindNodeCloseToAPoint({140.0, 0.0}, 1.0e-5); + // auto node62 = mesh.FindNodeCloseToAPoint({145.0, 15.0}, 1.0e-5); + // auto node63 = mesh.FindNodeCloseToAPoint({135.0, 15.0}, 1.0e-5); + std::cout << "nodes6: " << node61 << " " << node62 << " " << node63 << std::endl; + + // mesh.Administrate(); + std::vector polygonNodes; + auto boundaryNodes = mesh.ComputeBoundaryPolygons(polygonNodes); + + // std::vector boundaryNodes = meshBoundaryPolygon.GatherAllEnclosureNodes(); + + std::vector elementNodes1{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node11), + mesh.Node(node12), + mesh.Node(node13), + mesh.Node(node11)}; + + std::vector elementNodes2{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node21), + mesh.Node(node22), + mesh.Node(node23), + mesh.Node(node24), + mesh.Node(node21)}; + + std::vector elementNodes3{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node31), + mesh.Node(node32), + mesh.Node(node33), + mesh.Node(node34), + mesh.Node(node31)}; + + std::vector elementNodes4{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node41), + mesh.Node(node42), + mesh.Node(node43), + mesh.Node(node44), + mesh.Node(node41)}; + + std::vector elementNodes5{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node51), + mesh.Node(node52), + mesh.Node(node53), + mesh.Node(node51)}; + + std::vector elementNodes6{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node61), + mesh.Node(node62), + mesh.Node(node63), + mesh.Node(node61)}; + + // boundaryNodes.insert(boundaryNodes.end(), elementNodes1.begin(), elementNodes1.end()); + // boundaryNodes.insert(boundaryNodes.end(), elementNodes2.begin(), elementNodes2.end()); + // boundaryNodes.insert(boundaryNodes.end(), elementNodes3.begin(), elementNodes3.end()); + // boundaryNodes.insert(boundaryNodes.end(), elementNodes4.begin(), elementNodes4.end()); + // boundaryNodes.insert(boundaryNodes.end(), elementNodes5.begin(), elementNodes5.end()); + // boundaryNodes.insert(boundaryNodes.end(), elementNodes6.begin(), elementNodes6.end()); + + [[maybe_unused]] Polygons boundaryWithMissingElement(boundaryNodes, Projection::cartesian); + + // auto edge1 = mesh.FindEdge(node1, node2); + // auto edge2 = mesh.FindEdge(node2, node3); + // auto edge3 = mesh.FindEdge(node3, node1); + + // auto undoAction3 = ::Compute(mesh, Polygons()); + + //-------------------------------- + + // const auto projectToLandBoundaryOption = LandBoundaries::ProjectToLandBoundaryOption::DoNotProjectToLandBoundary; + // OrthogonalizationParameters orthogonalizationParameters; + // orthogonalizationParameters.inner_iterations = 2; + // orthogonalizationParameters.boundary_iterations = 25; + // orthogonalizationParameters.outer_iterations = 25; + // orthogonalizationParameters.orthogonalization_to_smoothing_factor = 0.975; + // orthogonalizationParameters.orthogonalization_to_smoothing_factor_at_boundary = 0.975; + // orthogonalizationParameters.areal_to_angle_smoothing_factor = 1.0; + + // // Execute + // auto orthoPolygon = std::make_unique(); + + // std::vector landBoundary{}; + // auto landboundaries = std::make_unique(landBoundary, mesh, *orthoPolygon); + + // OrthogonalizationAndSmoothing orthogonalization(mesh, + // std::move(orthoPolygon), + // std::move(landboundaries), + // projectToLandBoundaryOption, + // orthogonalizationParameters); + + // [[maybe_unused]] auto undoActionOrtho = orthogonalization.Initialize(); + // orthogonalization.Compute(); + + //-------------------------------- + + auto undoAction2 = mesh.DeleteMeshFacesInPolygons(boundaryWithMissingElement); + + std::vector boundaryNodes2 = mesh.ComputeInnerBoundaryPolygons(); + + // UInt polygonLength = 0; + // UInt polygonStart = 0; + + // UInt polygonLength2 = constants::missing::uintValue; + // UInt polygonStart2 = 0; + + // UInt index = 0; + + // while (index < boundaryNodes2.size()) + // { + // polygonStart2 = index; + + // for (UInt i = polygonStart2; i < boundaryNodes2.size(); ++i) + // { + // ++index; + + // if (!boundaryNodes2[i].IsValid()) + // { + // polygonLength2 = i - polygonStart2; + // break; + // } + // } + + // if (polygonLength2 > polygonLength) + // { + // polygonLength = polygonLength2; + // polygonStart = polygonStart2; + // } + // } + + // if (polygonLength < boundaryNodes2.size()) + // { + // ++polygonLength; + // } + + // std::cout << " max polygon length " << polygonLength << " " << polygonStart << std::endl; + + // std::vector boundaryNodes3; + + // if (polygonStart > 0) + // { + // boundaryNodes3.insert(boundaryNodes3.begin(), boundaryNodes2.begin(), boundaryNodes2.begin() + polygonStart - 1); + // } + + // if (polygonStart + polygonLength < boundaryNodes2.size()) + // { + // boundaryNodes3.insert(boundaryNodes3.end(), boundaryNodes2.begin() + polygonStart + polygonLength, boundaryNodes2.end()); + // } + + int count = 0; + + for (size_t i = 0; i < boundaryNodes2.size(); ++i) + { + std::cout << "point " << count << " " << " = " << boundaryNodes2[i].x << ", " << boundaryNodes2[i].y << std::endl; + ++count; + } + + [[maybe_unused]] Polygons boundaryWithMissingElement2(boundaryNodes2, Projection::cartesian); + + // for (size_t i = 0; i < boundaryNodes3.size(); ++i) + // { + // std::cout << "point " << count << " " << " = " << boundaryNodes3[i].x << ", " << boundaryNodes3[i].y << std::endl; + // ++count; + // } + + //-------------------------------- + + MeshRefinementParameters refParams; + refParams.max_num_refinement_iterations = 1; + + MeshRefinement meshRefinement(mesh, Polygons(), refParams); + auto undoAction5 = meshRefinement.Compute(); + + mesh.Administrate(); + + const auto projectToLandBoundaryOption = LandBoundaries::ProjectToLandBoundaryOption::DoNotProjectToLandBoundary; + OrthogonalizationParameters orthogonalizationParameters; + orthogonalizationParameters.inner_iterations = 2; + orthogonalizationParameters.boundary_iterations = 25; + orthogonalizationParameters.outer_iterations = 25; + orthogonalizationParameters.orthogonalization_to_smoothing_factor = 0.975; + orthogonalizationParameters.orthogonalization_to_smoothing_factor_at_boundary = 0.975; + orthogonalizationParameters.areal_to_angle_smoothing_factor = 1.0; + + // Execute + auto orthoPolygon = std::make_unique(); + + std::vector landBoundary{}; + auto landboundaries = std::make_unique(landBoundary, mesh, *orthoPolygon); + + OrthogonalizationAndSmoothing orthogonalization(mesh, + std::move(orthoPolygon), + std::move(landboundaries), + projectToLandBoundaryOption, + orthogonalizationParameters); + + [[maybe_unused]] auto undoActionOrtho = orthogonalization.Initialize(); + orthogonalization.Compute(); + + [[maybe_unused]] auto undoAction3 = mesh.DeleteMeshFacesInPolygons(boundaryWithMissingElement2); + + meshkernel::SaveVtk(mesh.Nodes(), mesh.m_facesNodes, "/home/wcs1/MeshKernel/MeshKernel04/build_deb/mesh1.vtu"); + // meshkernel::Print(mesh.Nodes(), mesh.Edges()); +} + +TEST(MeshRefinement, OMG) +{ + const size_t ExpectedNumberOfPoints = 184; + const size_t ExpectedNumberOfEdges = 360; + + auto curviMesh = MakeCurvilinearGrid(0.0, 0.0, 20.0, 20.0, 11, 11); + const auto edges = curviMesh->ComputeEdges(); + const auto nodes = curviMesh->ComputeNodes(); + + Mesh2D mesh(edges, nodes, Projection::cartesian); + + const std::vector originalNodes(mesh.Nodes()); + const std::vector originalEdges(mesh.Edges()); + + std::vector patch{{45.0, 5.0}, + {155.0, 5.0}, + {155.0, 155.0}, + {45.0, 155.0}, + {45.0, 5.0}, + {constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + {65.0, 65.0}, + {115.0, 65.0}, + {115.0, 115.0}, + {65.0, 115.0}, + {65.0, 65.0}}; + + std::vector expectedPoints(ExpectedNumberOfPoints); + std::vector expectedEdgeStart(ExpectedNumberOfEdges); + std::vector expectedEdgeEnd(ExpectedNumberOfEdges); + + Polygons polygon(patch, Projection::cartesian); + + auto undoAction = CasulliRefinement::Compute(mesh, polygon); + + auto node11 = mesh.FindNodeCloseToAPoint({80.0, 0.0}, 1.0e-5); + auto node12 = mesh.FindNodeCloseToAPoint({85.0, 15.0}, 1.0e-5); + auto node13 = mesh.FindNodeCloseToAPoint({75.0, 15.0}, 1.0e-5); + std::cout << "nodes1: " << node11 << " " << node12 << " " << node13 << std::endl; + + auto node21 = mesh.FindNodeCloseToAPoint({180.0, 140.0}, 1.0e-5); + auto node22 = mesh.FindNodeCloseToAPoint({200.0, 140.0}, 1.0e-5); + auto node23 = mesh.FindNodeCloseToAPoint({200.0, 160.0}, 1.0e-5); + auto node24 = mesh.FindNodeCloseToAPoint({180.0, 160.0}, 1.0e-5); + + std::cout << "nodes2: " << node21 << " " << node22 << " " << node23 << " " << node24 << std::endl; + + auto node31 = mesh.FindNodeCloseToAPoint({125.0, 125.0}, 1.0e-5); + auto node32 = mesh.FindNodeCloseToAPoint({135.0, 125.0}, 1.0e-5); + auto node33 = mesh.FindNodeCloseToAPoint({135.0, 135.0}, 1.0e-5); + auto node34 = mesh.FindNodeCloseToAPoint({125.0, 135.0}, 1.0e-5); + + std::cout << "nodes3: " << node31 << " " << node32 << " " << node33 << " " << node34 << std::endl; + + auto node41 = mesh.FindNodeCloseToAPoint({85.0, 15.0}, 1.0e-5); + auto node42 = mesh.FindNodeCloseToAPoint({95.0, 15.0}, 1.0e-5); + auto node43 = mesh.FindNodeCloseToAPoint({95.0, 25.0}, 1.0e-5); + auto node44 = mesh.FindNodeCloseToAPoint({85.0, 25.0}, 1.0e-5); + + std::cout << "nodes4: " << node41 << " " << node42 << " " << node43 << " " << node44 << std::endl; + + auto node51 = mesh.FindNodeCloseToAPoint({100.0, 0.0}, 1.0e-5); + auto node52 = mesh.FindNodeCloseToAPoint({105.0, 15.0}, 1.0e-5); + auto node53 = mesh.FindNodeCloseToAPoint({95.0, 15.0}, 1.0e-5); + std::cout << "nodes5: " << node51 << " " << node52 << " " << node53 << std::endl; + + auto node61 = mesh.FindNodeCloseToAPoint({120.0, 0.0}, 1.0e-5); + auto node62 = mesh.FindNodeCloseToAPoint({125.0, 15.0}, 1.0e-5); + auto node63 = mesh.FindNodeCloseToAPoint({115.0, 15.0}, 1.0e-5); + std::cout << "nodes6: " << node61 << " " << node62 << " " << node63 << std::endl; + + // mesh.Administrate(); + std::vector polygonNodes; + auto boundaryNodes = mesh.ComputeBoundaryPolygons(polygonNodes); + + // std::vector boundaryNodes = meshBoundaryPolygon.GatherAllEnclosureNodes(); + + std::vector elementNodes1{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node11), + mesh.Node(node12), + mesh.Node(node13), + mesh.Node(node11)}; + + std::vector elementNodes2{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node21), + mesh.Node(node22), + mesh.Node(node23), + mesh.Node(node24), + mesh.Node(node21)}; + + std::vector elementNodes3{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node31), + mesh.Node(node32), + mesh.Node(node33), + mesh.Node(node34), + mesh.Node(node31)}; + + std::vector elementNodes4{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node41), + mesh.Node(node42), + mesh.Node(node43), + mesh.Node(node44), + mesh.Node(node41)}; + + std::vector elementNodes5{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node51), + mesh.Node(node52), + mesh.Node(node53), + mesh.Node(node51)}; + + std::vector elementNodes6{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + mesh.Node(node61), + mesh.Node(node62), + mesh.Node(node63), + mesh.Node(node61)}; + + boundaryNodes.insert(boundaryNodes.end(), elementNodes1.begin(), elementNodes1.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes2.begin(), elementNodes2.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes3.begin(), elementNodes3.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes4.begin(), elementNodes4.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes5.begin(), elementNodes5.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes6.begin(), elementNodes6.end()); + + [[maybe_unused]] Polygons boundaryWithMissingElement(boundaryNodes, Projection::cartesian); + + [[maybe_unused]] auto undoAction3 = mesh.DeleteMeshFacesInPolygons(boundaryWithMissingElement); + + Mesh2D mesh2(mesh.Edges(), mesh.Nodes(), mesh.m_facesNodes, mesh.m_numFacesNodes, mesh.m_projection); + + meshkernel::SaveVtk(mesh2.Nodes(), mesh2.m_facesNodes, "/home/wcs1/MeshKernel/MeshKernel04/build_deb/mesh1.vtu"); +} From 4c672a8594e5d3c4655cfa35dd55335038b982ec Mon Sep 17 00:00:00 2001 From: BillSenior Date: Wed, 8 Oct 2025 11:59:24 +0200 Subject: [PATCH 06/16] GRIDEDIT-1795 Removed invalid element information from mesh2d after deleting elements --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 3 + libs/MeshKernel/src/Mesh2D.cpp | 89 ++++----- libs/MeshKernel/src/Polygons.cpp | 25 +-- .../tests/src/MeshRefinementTests.cpp | 173 ++++++------------ 4 files changed, 100 insertions(+), 190 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 981f541054..a9720772ea 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -484,6 +484,9 @@ namespace meshkernel std::vector& sortedNodes, std::vector& nodalValues); + /// @brief Update information about a list of faces and face-edge information + void UpdateFaceInformation(const std::vector& faceIndices, CompoundUndoAction& deleteMeshAction); + /// @brief Checks if a triangle has an acute angle (checktriangle) /// @param[in] faceNodes The face nodes composing the triangles /// @param[in] nodes The node coordinates diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index a40c903362..6354eaf30e 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -26,6 +26,7 @@ //------------------------------------------------------------------------------ #include +#include #include "MeshKernel/Constants.hpp" #include "MeshKernel/Definitions.hpp" @@ -2154,77 +2155,63 @@ std::unique_ptr Mesh2D::DeleteMeshFacesInPolygons(const std::unique_ptr deleteMeshAction = CompoundUndoAction::Create(); - Administrate(deleteMeshAction.get()); - - bool isInPolygon; - UInt polygonIndex; - - std::vector faceIndicesOffset(GetNumFaces(), 0); - std::vector faceIndices(GetNumFaces(), constants::missing::uintValue); - UInt faceIndexOffset = 0; - UInt faceIndex = 0; + std::vector faceIndices; + faceIndices.reserve(GetNumFaces()); for (UInt f = 0u; f < GetNumFaces(); ++f) { - isInPolygon = false; - polygonIndex = constants::missing::uintValue; if (!m_facesMassCenters[f].IsValid()) { continue; } - std::tie(isInPolygon, polygonIndex) = polygon.IsPointInPolygons(m_facesMassCenters[f]); - - // TODO I think the polygon is orientated in the wrong direction - if (!isInPolygon) + if (polygon.IsPointInAnyPolygon(m_facesMassCenters[f])) { + faceIndices.push_back(f); + } + } - for (UInt e = 0u; e < m_facesEdges[f].size(); ++e) - { - UInt edge = m_facesEdges[f][e]; + UpdateFaceInformation(faceIndices, *deleteMeshAction); - if (m_edgesFaces[edge][0] == f) - { - m_edgesFaces[edge][0] = constants::missing::uintValue; - --m_edgesNumFaces[edge]; - } - else if (m_edgesFaces[edge][1] == f) - { - m_edgesFaces[edge][1] = constants::missing::uintValue; - --m_edgesNumFaces[edge]; - } + return deleteMeshAction; +} - if (m_edgesNumFaces[edge] == 0) - { - deleteMeshAction->Add(DeleteEdge(edge)); - } +void Mesh2D::UpdateFaceInformation(const std::vector& faceIndices, CompoundUndoAction& deleteMeshAction) +{ + for (UInt faceId : faceIndices) + { + for (UInt e = 0u; e < m_facesEdges[faceId].size(); ++e) + { + UInt edge = m_facesEdges[faceId][e]; + + if (m_edgesFaces[edge][0] == faceId) + { + m_edgesFaces[edge][0] = constants::missing::uintValue; + --m_edgesNumFaces[edge]; + } + else if (m_edgesFaces[edge][1] == faceId) + { + m_edgesFaces[edge][1] = constants::missing::uintValue; + --m_edgesNumFaces[edge]; } - ++faceIndexOffset; - m_facesNodes[f].clear(); - m_numFacesNodes[f] = 0; - m_facesEdges[f].clear(); - m_facesMassCenters[f] = {constants::missing::doubleValue, constants::missing::doubleValue}; - m_faceArea[f] = constants::missing::doubleValue; - } - else - { - faceIndices[f] = faceIndex; - ++faceIndex; + if (m_edgesNumFaces[edge] == 0) + { + deleteMeshAction.Add(DeleteEdge(edge)); + } } - - faceIndicesOffset[f] = faceIndexOffset; } - std::cout << "faceIndicesOffset" << std::endl; - - for (size_t i = 0; i < faceIndicesOffset.size(); ++i) + // TODO need to collect undo information? + for (UInt faceId : faceIndices | std::views::reverse) { - std::cout << " faceIndicesOffset " << i << " = " << faceIndicesOffset[i] << " " << faceIndices[i] << std::endl; + m_facesNodes.erase(m_facesNodes.begin() + faceId); + m_numFacesNodes.erase(m_numFacesNodes.begin() + faceId); + m_facesEdges.erase(m_facesEdges.begin() + faceId); + m_facesMassCenters.erase(m_facesMassCenters.begin() + faceId); + m_faceArea.erase(m_faceArea.begin() + faceId); } - - return deleteMeshAction; } std::unique_ptr Mesh2D::DeleteHangingEdges() diff --git a/libs/MeshKernel/src/Polygons.cpp b/libs/MeshKernel/src/Polygons.cpp index ed58f50749..de98fa1d2b 100644 --- a/libs/MeshKernel/src/Polygons.cpp +++ b/libs/MeshKernel/src/Polygons.cpp @@ -309,29 +309,8 @@ std::tuple Polygons::IsPointInPolygons(const Point& poin bool Polygons::IsPointInAnyPolygon(const Point& point) const { - // empty polygon means everything is included - if (IsEmpty()) - { - return true; - } - - for (UInt i = 0; i < m_enclosures.size(); ++i) - { - PolygonalEnclosure::Region containingRegion = m_enclosures[i].ContainsRegion(point); - - if (containingRegion == PolygonalEnclosure::Region::Exterior) - { - // Point was found in - return true; - } - else if (containingRegion == PolygonalEnclosure::Region::Interior) - { - // Point can be found in an hole in the polygon - break; - } - } - - return false; + auto [isInPolygon, polygonIndex] = IsPointInPolygons(point); + return isInPolygon; } std::vector Polygons::PointsInPolygons(const std::vector& points) const diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index 65eb202c67..71384b8af7 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -2836,43 +2836,42 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygonsForSixFa auto node62 = mesh.FindNodeCloseToAPoint({125.0, 15.0}, 1.0e-5); auto node63 = mesh.FindNodeCloseToAPoint({115.0, 15.0}, 1.0e-5); - std::vector polygonNodes; - auto boundaryNodes = mesh.ComputeBoundaryPolygons(polygonNodes); + std::vector boundaryNodes; - std::vector elementNodes1{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes1{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node11), mesh.Node(node12), mesh.Node(node13), mesh.Node(node11)}; - std::vector elementNodes2{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes2{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node21), mesh.Node(node22), mesh.Node(node23), mesh.Node(node24), mesh.Node(node21)}; - std::vector elementNodes3{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes3{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node31), mesh.Node(node32), mesh.Node(node33), mesh.Node(node34), mesh.Node(node31)}; - std::vector elementNodes4{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes4{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node41), mesh.Node(node42), mesh.Node(node43), mesh.Node(node44), mesh.Node(node41)}; - std::vector elementNodes5{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes5{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node51), mesh.Node(node52), mesh.Node(node53), mesh.Node(node51)}; - std::vector elementNodes6{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes6{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node61), mesh.Node(node62), mesh.Node(node63), @@ -2886,10 +2885,9 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygonsForSixFa boundaryNodes.insert(boundaryNodes.end(), elementNodes5.begin(), elementNodes5.end()); boundaryNodes.insert(boundaryNodes.end(), elementNodes6.begin(), elementNodes6.end()); - Polygons boundaryWithMissingElement(boundaryNodes, Projection::cartesian); + Polygons boundaryWithMissingElements(boundaryNodes, Projection::cartesian); - auto deleteMeshFacesUndoAction = mesh.DeleteMeshFacesInPolygons(boundaryWithMissingElement); - meshkernel::SaveVtk(mesh.Nodes(), mesh.m_facesNodes, "/home/wcs1/MeshKernel/MeshKernel04/build_deb/mesh1.vtu"); + auto deleteMeshFacesUndoAction = mesh.DeleteMeshFacesInPolygons(boundaryWithMissingElements); // Compute interior boundary polygon points std::vector boundaryNodes2 = mesh.ComputeInnerBoundaryPolygons(); @@ -3000,132 +2998,85 @@ TEST(MeshRefinement, WTF) auto node11 = mesh.FindNodeCloseToAPoint({80.0, 0.0}, 1.0e-5); auto node12 = mesh.FindNodeCloseToAPoint({85.0, 15.0}, 1.0e-5); auto node13 = mesh.FindNodeCloseToAPoint({75.0, 15.0}, 1.0e-5); - std::cout << "nodes1: " << node11 << " " << node12 << " " << node13 << std::endl; auto node21 = mesh.FindNodeCloseToAPoint({180.0, 140.0}, 1.0e-5); auto node22 = mesh.FindNodeCloseToAPoint({200.0, 140.0}, 1.0e-5); auto node23 = mesh.FindNodeCloseToAPoint({200.0, 160.0}, 1.0e-5); auto node24 = mesh.FindNodeCloseToAPoint({180.0, 160.0}, 1.0e-5); - std::cout << "nodes2: " << node21 << " " << node22 << " " << node23 << " " << node24 << std::endl; - auto node31 = mesh.FindNodeCloseToAPoint({125.0, 125.0}, 1.0e-5); auto node32 = mesh.FindNodeCloseToAPoint({135.0, 125.0}, 1.0e-5); auto node33 = mesh.FindNodeCloseToAPoint({135.0, 135.0}, 1.0e-5); auto node34 = mesh.FindNodeCloseToAPoint({125.0, 135.0}, 1.0e-5); - std::cout << "nodes3: " << node31 << " " << node32 << " " << node33 << " " << node34 << std::endl; - auto node41 = mesh.FindNodeCloseToAPoint({85.0, 15.0}, 1.0e-5); auto node42 = mesh.FindNodeCloseToAPoint({95.0, 15.0}, 1.0e-5); auto node43 = mesh.FindNodeCloseToAPoint({95.0, 25.0}, 1.0e-5); auto node44 = mesh.FindNodeCloseToAPoint({85.0, 25.0}, 1.0e-5); - std::cout << "nodes4: " << node41 << " " << node42 << " " << node43 << " " << node44 << std::endl; - auto node51 = mesh.FindNodeCloseToAPoint({100.0, 0.0}, 1.0e-5); auto node52 = mesh.FindNodeCloseToAPoint({105.0, 15.0}, 1.0e-5); auto node53 = mesh.FindNodeCloseToAPoint({95.0, 15.0}, 1.0e-5); - std::cout << "nodes5: " << node51 << " " << node52 << " " << node53 << std::endl; auto node61 = mesh.FindNodeCloseToAPoint({120.0, 0.0}, 1.0e-5); auto node62 = mesh.FindNodeCloseToAPoint({125.0, 15.0}, 1.0e-5); auto node63 = mesh.FindNodeCloseToAPoint({115.0, 15.0}, 1.0e-5); - // auto node61 = mesh.FindNodeCloseToAPoint({140.0, 0.0}, 1.0e-5); - // auto node62 = mesh.FindNodeCloseToAPoint({145.0, 15.0}, 1.0e-5); - // auto node63 = mesh.FindNodeCloseToAPoint({135.0, 15.0}, 1.0e-5); - std::cout << "nodes6: " << node61 << " " << node62 << " " << node63 << std::endl; - // mesh.Administrate(); - std::vector polygonNodes; - auto boundaryNodes = mesh.ComputeBoundaryPolygons(polygonNodes); + std::vector boundaryNodes; - // std::vector boundaryNodes = meshBoundaryPolygon.GatherAllEnclosureNodes(); - - std::vector elementNodes1{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes1{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node11), mesh.Node(node12), mesh.Node(node13), mesh.Node(node11)}; - std::vector elementNodes2{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes2{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node21), mesh.Node(node22), mesh.Node(node23), mesh.Node(node24), mesh.Node(node21)}; - std::vector elementNodes3{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes3{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node31), mesh.Node(node32), mesh.Node(node33), mesh.Node(node34), mesh.Node(node31)}; - std::vector elementNodes4{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes4{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node41), mesh.Node(node42), mesh.Node(node43), mesh.Node(node44), mesh.Node(node41)}; - std::vector elementNodes5{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes5{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node51), mesh.Node(node52), mesh.Node(node53), mesh.Node(node51)}; - std::vector elementNodes6{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes6{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node61), mesh.Node(node62), mesh.Node(node63), mesh.Node(node61)}; - // boundaryNodes.insert(boundaryNodes.end(), elementNodes1.begin(), elementNodes1.end()); - // boundaryNodes.insert(boundaryNodes.end(), elementNodes2.begin(), elementNodes2.end()); - // boundaryNodes.insert(boundaryNodes.end(), elementNodes3.begin(), elementNodes3.end()); - // boundaryNodes.insert(boundaryNodes.end(), elementNodes4.begin(), elementNodes4.end()); - // boundaryNodes.insert(boundaryNodes.end(), elementNodes5.begin(), elementNodes5.end()); - // boundaryNodes.insert(boundaryNodes.end(), elementNodes6.begin(), elementNodes6.end()); - - [[maybe_unused]] Polygons boundaryWithMissingElement(boundaryNodes, Projection::cartesian); - - // auto edge1 = mesh.FindEdge(node1, node2); - // auto edge2 = mesh.FindEdge(node2, node3); - // auto edge3 = mesh.FindEdge(node3, node1); - - // auto undoAction3 = ::Compute(mesh, Polygons()); - - //-------------------------------- - - // const auto projectToLandBoundaryOption = LandBoundaries::ProjectToLandBoundaryOption::DoNotProjectToLandBoundary; - // OrthogonalizationParameters orthogonalizationParameters; - // orthogonalizationParameters.inner_iterations = 2; - // orthogonalizationParameters.boundary_iterations = 25; - // orthogonalizationParameters.outer_iterations = 25; - // orthogonalizationParameters.orthogonalization_to_smoothing_factor = 0.975; - // orthogonalizationParameters.orthogonalization_to_smoothing_factor_at_boundary = 0.975; - // orthogonalizationParameters.areal_to_angle_smoothing_factor = 1.0; - - // // Execute - // auto orthoPolygon = std::make_unique(); - - // std::vector landBoundary{}; - // auto landboundaries = std::make_unique(landBoundary, mesh, *orthoPolygon); - - // OrthogonalizationAndSmoothing orthogonalization(mesh, - // std::move(orthoPolygon), - // std::move(landboundaries), - // projectToLandBoundaryOption, - // orthogonalizationParameters); + // Combine all nodes to for a sequence of polygons + boundaryNodes.insert(boundaryNodes.end(), elementNodes1.begin(), elementNodes1.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes2.begin(), elementNodes2.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes3.begin(), elementNodes3.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes4.begin(), elementNodes4.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes5.begin(), elementNodes5.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes6.begin(), elementNodes6.end()); - // [[maybe_unused]] auto undoActionOrtho = orthogonalization.Initialize(); - // orthogonalization.Compute(); + [[maybe_unused]] Polygons boundaryWithMissingElements(boundaryNodes, Projection::cartesian); + auto undoAction2 = mesh.DeleteMeshFacesInPolygons(boundaryWithMissingElements); + meshkernel::SaveVtk(mesh.Nodes(), mesh.m_facesNodes, "mesh1.vtu"); //-------------------------------- - auto undoAction2 = mesh.DeleteMeshFacesInPolygons(boundaryWithMissingElement); - std::vector boundaryNodes2 = mesh.ComputeInnerBoundaryPolygons(); // UInt polygonLength = 0; @@ -3179,6 +3130,8 @@ TEST(MeshRefinement, WTF) int count = 0; + std::cout << "boundaryNodes2.size " << boundaryNodes2.size() << std::endl; + for (size_t i = 0; i < boundaryNodes2.size(); ++i) { std::cout << "point " << count << " " << " = " << boundaryNodes2[i].x << ", " << boundaryNodes2[i].y << std::endl; @@ -3203,33 +3156,33 @@ TEST(MeshRefinement, WTF) mesh.Administrate(); - const auto projectToLandBoundaryOption = LandBoundaries::ProjectToLandBoundaryOption::DoNotProjectToLandBoundary; - OrthogonalizationParameters orthogonalizationParameters; - orthogonalizationParameters.inner_iterations = 2; - orthogonalizationParameters.boundary_iterations = 25; - orthogonalizationParameters.outer_iterations = 25; - orthogonalizationParameters.orthogonalization_to_smoothing_factor = 0.975; - orthogonalizationParameters.orthogonalization_to_smoothing_factor_at_boundary = 0.975; - orthogonalizationParameters.areal_to_angle_smoothing_factor = 1.0; + // const auto projectToLandBoundaryOption = LandBoundaries::ProjectToLandBoundaryOption::DoNotProjectToLandBoundary; + // OrthogonalizationParameters orthogonalizationParameters; + // orthogonalizationParameters.inner_iterations = 2; + // orthogonalizationParameters.boundary_iterations = 25; + // orthogonalizationParameters.outer_iterations = 25; + // orthogonalizationParameters.orthogonalization_to_smoothing_factor = 0.975; + // orthogonalizationParameters.orthogonalization_to_smoothing_factor_at_boundary = 0.975; + // orthogonalizationParameters.areal_to_angle_smoothing_factor = 1.0; - // Execute - auto orthoPolygon = std::make_unique(); + // // Execute + // auto orthoPolygon = std::make_unique(); - std::vector landBoundary{}; - auto landboundaries = std::make_unique(landBoundary, mesh, *orthoPolygon); + // std::vector landBoundary{}; + // auto landboundaries = std::make_unique(landBoundary, mesh, *orthoPolygon); - OrthogonalizationAndSmoothing orthogonalization(mesh, - std::move(orthoPolygon), - std::move(landboundaries), - projectToLandBoundaryOption, - orthogonalizationParameters); + // OrthogonalizationAndSmoothing orthogonalization(mesh, + // std::move(orthoPolygon), + // std::move(landboundaries), + // projectToLandBoundaryOption, + // orthogonalizationParameters); - [[maybe_unused]] auto undoActionOrtho = orthogonalization.Initialize(); - orthogonalization.Compute(); + // [[maybe_unused]] auto undoActionOrtho = orthogonalization.Initialize(); + // orthogonalization.Compute(); [[maybe_unused]] auto undoAction3 = mesh.DeleteMeshFacesInPolygons(boundaryWithMissingElement2); - meshkernel::SaveVtk(mesh.Nodes(), mesh.m_facesNodes, "/home/wcs1/MeshKernel/MeshKernel04/build_deb/mesh1.vtu"); + meshkernel::SaveVtk(mesh.Nodes(), mesh.m_facesNodes, "mesh2.vtu"); // meshkernel::Print(mesh.Nodes(), mesh.Edges()); } @@ -3270,84 +3223,72 @@ TEST(MeshRefinement, OMG) auto node11 = mesh.FindNodeCloseToAPoint({80.0, 0.0}, 1.0e-5); auto node12 = mesh.FindNodeCloseToAPoint({85.0, 15.0}, 1.0e-5); auto node13 = mesh.FindNodeCloseToAPoint({75.0, 15.0}, 1.0e-5); - std::cout << "nodes1: " << node11 << " " << node12 << " " << node13 << std::endl; auto node21 = mesh.FindNodeCloseToAPoint({180.0, 140.0}, 1.0e-5); auto node22 = mesh.FindNodeCloseToAPoint({200.0, 140.0}, 1.0e-5); auto node23 = mesh.FindNodeCloseToAPoint({200.0, 160.0}, 1.0e-5); auto node24 = mesh.FindNodeCloseToAPoint({180.0, 160.0}, 1.0e-5); - std::cout << "nodes2: " << node21 << " " << node22 << " " << node23 << " " << node24 << std::endl; - auto node31 = mesh.FindNodeCloseToAPoint({125.0, 125.0}, 1.0e-5); auto node32 = mesh.FindNodeCloseToAPoint({135.0, 125.0}, 1.0e-5); auto node33 = mesh.FindNodeCloseToAPoint({135.0, 135.0}, 1.0e-5); auto node34 = mesh.FindNodeCloseToAPoint({125.0, 135.0}, 1.0e-5); - std::cout << "nodes3: " << node31 << " " << node32 << " " << node33 << " " << node34 << std::endl; - auto node41 = mesh.FindNodeCloseToAPoint({85.0, 15.0}, 1.0e-5); auto node42 = mesh.FindNodeCloseToAPoint({95.0, 15.0}, 1.0e-5); auto node43 = mesh.FindNodeCloseToAPoint({95.0, 25.0}, 1.0e-5); auto node44 = mesh.FindNodeCloseToAPoint({85.0, 25.0}, 1.0e-5); - std::cout << "nodes4: " << node41 << " " << node42 << " " << node43 << " " << node44 << std::endl; - auto node51 = mesh.FindNodeCloseToAPoint({100.0, 0.0}, 1.0e-5); auto node52 = mesh.FindNodeCloseToAPoint({105.0, 15.0}, 1.0e-5); auto node53 = mesh.FindNodeCloseToAPoint({95.0, 15.0}, 1.0e-5); - std::cout << "nodes5: " << node51 << " " << node52 << " " << node53 << std::endl; auto node61 = mesh.FindNodeCloseToAPoint({120.0, 0.0}, 1.0e-5); auto node62 = mesh.FindNodeCloseToAPoint({125.0, 15.0}, 1.0e-5); auto node63 = mesh.FindNodeCloseToAPoint({115.0, 15.0}, 1.0e-5); - std::cout << "nodes6: " << node61 << " " << node62 << " " << node63 << std::endl; - // mesh.Administrate(); - std::vector polygonNodes; - auto boundaryNodes = mesh.ComputeBoundaryPolygons(polygonNodes); + std::vector boundaryNodes; - // std::vector boundaryNodes = meshBoundaryPolygon.GatherAllEnclosureNodes(); - - std::vector elementNodes1{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes1{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node11), mesh.Node(node12), mesh.Node(node13), mesh.Node(node11)}; - std::vector elementNodes2{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes2{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node21), mesh.Node(node22), mesh.Node(node23), mesh.Node(node24), mesh.Node(node21)}; - std::vector elementNodes3{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes3{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node31), mesh.Node(node32), mesh.Node(node33), mesh.Node(node34), mesh.Node(node31)}; - std::vector elementNodes4{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes4{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node41), mesh.Node(node42), mesh.Node(node43), mesh.Node(node44), mesh.Node(node41)}; - std::vector elementNodes5{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes5{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node51), mesh.Node(node52), mesh.Node(node53), mesh.Node(node51)}; - std::vector elementNodes6{{constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + std::vector elementNodes6{{constants::missing::doubleValue, constants::missing::doubleValue}, mesh.Node(node61), mesh.Node(node62), mesh.Node(node63), mesh.Node(node61)}; + // Combine all nodes to for a sequence of polygons boundaryNodes.insert(boundaryNodes.end(), elementNodes1.begin(), elementNodes1.end()); boundaryNodes.insert(boundaryNodes.end(), elementNodes2.begin(), elementNodes2.end()); boundaryNodes.insert(boundaryNodes.end(), elementNodes3.begin(), elementNodes3.end()); @@ -3361,5 +3302,5 @@ TEST(MeshRefinement, OMG) Mesh2D mesh2(mesh.Edges(), mesh.Nodes(), mesh.m_facesNodes, mesh.m_numFacesNodes, mesh.m_projection); - meshkernel::SaveVtk(mesh2.Nodes(), mesh2.m_facesNodes, "/home/wcs1/MeshKernel/MeshKernel04/build_deb/mesh1.vtu"); + meshkernel::SaveVtk(mesh2.Nodes(), mesh2.m_facesNodes, "mesh1.vtu"); } From 8e001fdd8df63cd231111a7845584d405c5cc748 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Wed, 8 Oct 2025 15:49:44 +0200 Subject: [PATCH 07/16] GRIDEDIT-1795 Update edge-face connectivity info when one of more faces has been deleted --- libs/MeshKernel/include/MeshKernel/Mesh.hpp | 6 +-- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 6 ++- libs/MeshKernel/src/Mesh2D.cpp | 45 ++++++++++++++++--- .../tests/src/MeshRefinementTests.cpp | 3 +- 4 files changed, 47 insertions(+), 13 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh.hpp b/libs/MeshKernel/include/MeshKernel/Mesh.hpp index f330add6a6..e7a1b19675 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh.hpp @@ -319,6 +319,9 @@ namespace meshkernel /// @brief Perform node and edges administration void AdministrateNodesEdges(CompoundUndoAction* undoAction = nullptr); + /// @brief Determine if a administration is required + bool AdministrationRequired() const; + /// @brief Sort mesh edges around a node in counterclockwise order (Sort_links_ccw) /// @param[in] startNode The first node index where to perform edge sorting. /// @param[in] endNode The last node index where to perform edge sorting. @@ -477,9 +480,6 @@ namespace meshkernel Projection m_projection; ///< The projection used protected: - /// @brief Determine if a administration is required - bool AdministrationRequired() const; - /// @brief Determine if a administration is required void SetAdministrationRequired(const bool value); diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index a9720772ea..0536529b53 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -375,8 +375,7 @@ namespace meshkernel /// @brief Deletes the mesh faces inside a set of polygons /// @param[in] polygon The polygon where to perform the operation - /// If this Polygons instance contains multiple polygons, the first one will be taken. - std::unique_ptr DeleteMeshFacesInPolygons(const Polygons& polygon); + /* [[nodiscard]] */ std::unique_ptr DeleteMeshFacesInPolygons(const Polygons& polygon); private: // orthogonalization @@ -485,6 +484,9 @@ namespace meshkernel std::vector& nodalValues); /// @brief Update information about a list of faces and face-edge information + /// + /// @param [in] faceIndices Mapping between old and new face-ids, the invalid value implies a deleted element + /// @param [in,out] deleteMeshAction Gather any additional undo information void UpdateFaceInformation(const std::vector& faceIndices, CompoundUndoAction& deleteMeshAction); /// @brief Checks if a triangle has an acute angle (checktriangle) diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 6354eaf30e..2020fd30e5 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -2154,9 +2154,10 @@ std::unique_ptr Mesh2D::DeleteMeshFacesInPolygons(const } std::unique_ptr deleteMeshAction = CompoundUndoAction::Create(); + ComputeFaceAreaAndMassCenters(true); - std::vector faceIndices; - faceIndices.reserve(GetNumFaces()); + std::vector faceIndices(GetNumFaces(), constants::missing::uintValue); + UInt faceIndex = 0; for (UInt f = 0u; f < GetNumFaces(); ++f) { @@ -2166,9 +2167,10 @@ std::unique_ptr Mesh2D::DeleteMeshFacesInPolygons(const continue; } - if (polygon.IsPointInAnyPolygon(m_facesMassCenters[f])) + if (!polygon.IsPointInAnyPolygon(m_facesMassCenters[f])) { - faceIndices.push_back(f); + faceIndices[f] = faceIndex; + ++faceIndex; } } @@ -2179,7 +2181,20 @@ std::unique_ptr Mesh2D::DeleteMeshFacesInPolygons(const void Mesh2D::UpdateFaceInformation(const std::vector& faceIndices, CompoundUndoAction& deleteMeshAction) { - for (UInt faceId : faceIndices) + std::vector facesToDelete; + facesToDelete.reserve(GetNumFaces()); + + // Collect face-ids of faces that are marked for deletion + for (UInt i = 0; i < faceIndices.size(); ++i) + { + if (faceIndices[i] == constants::missing::uintValue) + { + facesToDelete.push_back(i); + } + } + + // Remove deleted faces from edge-face connectivity + for (UInt faceId : facesToDelete) { for (UInt e = 0u; e < m_facesEdges[faceId].size(); ++e) { @@ -2189,6 +2204,9 @@ void Mesh2D::UpdateFaceInformation(const std::vector& faceIndices, Compoun { m_edgesFaces[edge][0] = constants::missing::uintValue; --m_edgesNumFaces[edge]; + // Ensure the first connected face is on the 0th side + // If the 1st side is already the invalid value then this operation does not change anything. + std::swap(m_edgesFaces[edge][0], m_edgesFaces[edge][1]); } else if (m_edgesFaces[edge][1] == faceId) { @@ -2203,8 +2221,23 @@ void Mesh2D::UpdateFaceInformation(const std::vector& faceIndices, Compoun } } + // Renumber existing edge-face ids to match new face-ids + for (size_t edge = 0; edge < m_edgesFaces.size(); ++edge) + { + if (m_edgesFaces[edge][0] != constants::missing::uintValue) + { + m_edgesFaces[edge][0] = faceIndices[m_edgesFaces[edge][0]]; + } + + if (m_edgesFaces[edge][1] != constants::missing::uintValue) + { + m_edgesFaces[edge][1] = faceIndices[m_edgesFaces[edge][1]]; + } + } + // TODO need to collect undo information? - for (UInt faceId : faceIndices | std::views::reverse) + // Shift face connectivity in arrays where deleted faces have been removed from arrays + for (UInt faceId : facesToDelete | std::views::reverse) { m_facesNodes.erase(m_facesNodes.begin() + faceId); m_numFacesNodes.erase(m_numFacesNodes.begin() + faceId); diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index 71384b8af7..315c405fb5 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -3024,8 +3024,7 @@ TEST(MeshRefinement, WTF) std::vector boundaryNodes; - std::vector elementNodes1{{constants::missing::doubleValue, constants::missing::doubleValue}, - mesh.Node(node11), + std::vector elementNodes1{mesh.Node(node11), mesh.Node(node12), mesh.Node(node13), mesh.Node(node11)}; From dd21ce2cb0b6c0f94e7c706b201ab9940c76feb5 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 9 Oct 2025 09:35:42 +0200 Subject: [PATCH 08/16] GRIDEDIT-1795 Small refactoring and added comments --- libs/MeshKernel/src/Mesh2D.cpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 2020fd30e5..ffb9ed1d4d 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -2156,18 +2156,14 @@ std::unique_ptr Mesh2D::DeleteMeshFacesInPolygons(const std::unique_ptr deleteMeshAction = CompoundUndoAction::Create(); ComputeFaceAreaAndMassCenters(true); + // A mapping between the old face-id and the new face-id. + // Any deleted elements will be the invalid uint value. std::vector faceIndices(GetNumFaces(), constants::missing::uintValue); UInt faceIndex = 0; for (UInt f = 0u; f < GetNumFaces(); ++f) { - - if (!m_facesMassCenters[f].IsValid()) - { - continue; - } - - if (!polygon.IsPointInAnyPolygon(m_facesMassCenters[f])) + if (m_facesMassCenters[f].IsValid() && !polygon.IsPointInAnyPolygon(m_facesMassCenters[f])) { faceIndices[f] = faceIndex; ++faceIndex; @@ -2199,6 +2195,7 @@ void Mesh2D::UpdateFaceInformation(const std::vector& faceIndices, Compoun for (UInt e = 0u; e < m_facesEdges[faceId].size(); ++e) { UInt edge = m_facesEdges[faceId][e]; + bool removedFace = false; if (m_edgesFaces[edge][0] == faceId) { @@ -2207,14 +2204,16 @@ void Mesh2D::UpdateFaceInformation(const std::vector& faceIndices, Compoun // Ensure the first connected face is on the 0th side // If the 1st side is already the invalid value then this operation does not change anything. std::swap(m_edgesFaces[edge][0], m_edgesFaces[edge][1]); + removedFace = true; } else if (m_edgesFaces[edge][1] == faceId) { m_edgesFaces[edge][1] = constants::missing::uintValue; --m_edgesNumFaces[edge]; + removedFace = true; } - if (m_edgesNumFaces[edge] == 0) + if (removedFace && m_edgesNumFaces[edge] == 0) { deleteMeshAction.Add(DeleteEdge(edge)); } From 61a76bcbd1054b9371517b136ec31d58b1915cab Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 9 Oct 2025 09:36:29 +0200 Subject: [PATCH 09/16] GRIDEDIT-1795 Apply invalid cells mask after steps that will call an administrate or change the mesh --- .../include/MeshKernelApi/MeshKernel.hpp | 5 + .../include/MeshKernelApi/State.hpp | 12 +- libs/MeshKernelApi/src/MKStateUndoAction.cpp | 2 + libs/MeshKernelApi/src/MeshKernel.cpp | 228 +++++++++++++++++- 4 files changed, 236 insertions(+), 11 deletions(-) diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp index cb41751f4d..3c8cf1fe5b 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp @@ -1150,6 +1150,11 @@ namespace meshkernelapi /// @returns Error code MKERNEL_API int mkernel_mesh2d_delete_hanging_edges(int meshKernelId); + /// @brief Deletes + /// @param[in] meshKernelId The id of the mesh state + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_delete_faces_in_polygons(int meshKernelId, const GeometryList& polygon); + /// @brief Deletes a mesh2d node /// @param[in] meshKernelId The id of the mesh state /// @param[in] nodeIndex The index of the node to delete diff --git a/libs/MeshKernelApi/include/MeshKernelApi/State.hpp b/libs/MeshKernelApi/include/MeshKernelApi/State.hpp index 45ca154349..744c7c9fff 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/State.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/State.hpp @@ -29,6 +29,7 @@ #include #include +#include #include "MeshKernel/Contacts.hpp" #include "MeshKernel/CurvilinearGrid/CurvilinearGridFromSplines.hpp" @@ -72,11 +73,12 @@ namespace meshkernelapi } // Geometrical entities instances - std::shared_ptr m_mesh1d; ///< Shared pointer to meshkernel::Mesh1D instance - std::shared_ptr m_network1d; ///< Shared pointer to meshkernel::Network1D instance - std::shared_ptr m_mesh2d; ///< Shared pointer to meshkernel::Mesh2D instance - std::shared_ptr m_contacts; ///< Shared pointer to meshkernel::Contacts instance - std::shared_ptr m_curvilinearGrid; ///< Shared pointer to meshkernel::CurvilinearGrid instance + std::shared_ptr m_mesh1d; ///< Shared pointer to meshkernel::Mesh1D instance + std::shared_ptr m_network1d; ///< Shared pointer to meshkernel::Network1D instance + std::shared_ptr m_mesh2d; ///< Shared pointer to meshkernel::Mesh2D instance + std::shared_ptr m_contacts; ///< Shared pointer to meshkernel::Contacts instance + std::shared_ptr m_curvilinearGrid; ///< Shared pointer to meshkernel::CurvilinearGrid instance + std::shared_ptr> m_invalidCellPolygons; ///< Set of polygon points for invalid regions or holes in mesh // Algorithms instances (interactivity) std::shared_ptr m_meshOrthogonalization; ///< Shared pointer to meshkernel::OrthogonalizationAndSmoothing instance diff --git a/libs/MeshKernelApi/src/MKStateUndoAction.cpp b/libs/MeshKernelApi/src/MKStateUndoAction.cpp index a233c3467f..42cd83c2b9 100644 --- a/libs/MeshKernelApi/src/MKStateUndoAction.cpp +++ b/libs/MeshKernelApi/src/MKStateUndoAction.cpp @@ -43,6 +43,7 @@ meshkernelapi::MKStateUndoAction::MKStateUndoAction(MeshKernelState& mkState) : m_mkState.m_curvilinearGridFromSplines = mkState.m_curvilinearGridFromSplines; m_mkState.m_curvilinearGridLineShift = mkState.m_curvilinearGridLineShift; m_mkState.m_projection = mkState.m_projection; + m_mkState.m_invalidCellPolygons = mkState.m_invalidCellPolygons; } void meshkernelapi::MKStateUndoAction::SwapContents() @@ -52,6 +53,7 @@ void meshkernelapi::MKStateUndoAction::SwapContents() std::swap(m_mkState.m_network1d, m_mkStateReference.m_network1d); std::swap(m_mkState.m_contacts, m_mkStateReference.m_contacts); std::swap(m_mkState.m_curvilinearGrid, m_mkStateReference.m_curvilinearGrid); + std::swap(m_mkState.m_invalidCellPolygons, m_mkStateReference.m_invalidCellPolygons); } void meshkernelapi::MKStateUndoAction::DoCommit() diff --git a/libs/MeshKernelApi/src/MeshKernel.cpp b/libs/MeshKernelApi/src/MeshKernel.cpp index f2321c6454..3bd7e7215b 100644 --- a/libs/MeshKernelApi/src/MeshKernel.cpp +++ b/libs/MeshKernelApi/src/MeshKernel.cpp @@ -88,6 +88,7 @@ #include #include #include +#include #include "MeshKernelApi/ApiCache/CachedPointValues.hpp" #include "MeshKernelApi/ApiCache/CurvilinearBoundariesAsPolygonCache.hpp" @@ -203,6 +204,7 @@ namespace meshkernelapi auto const projection = static_cast(projectionType); meshKernelState.insert({meshKernelId, MeshKernelState(projection)}); meshKernelState[meshKernelId].m_propertyCalculators = allocateDefaultPropertyCalculators(); + meshKernelState[meshKernelId].m_invalidCellPolygons = std::make_shared>(); } catch (...) { @@ -458,7 +460,18 @@ namespace meshkernelapi const auto deletionOptionEnum = static_cast(deletionOption); meshKernelUndoStack.Add(meshKernelState[meshKernelId].m_mesh2d->DeleteMesh(meshKernelPolygon, deletionOptionEnum, invertDeletionBool), meshKernelId); + + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } + + // Invalid cells may need to be updated after deleting + *meshKernelState[meshKernelId].m_invalidCellPolygons = meshKernelState[meshKernelId].m_mesh2d->ComputeInnerBoundaryPolygons(); } + catch (...) { lastExitCode = HandleException(); @@ -515,6 +528,7 @@ namespace meshkernelapi meshKernelState[meshKernelId].m_projection); } + *meshKernelState[meshKernelId].m_invalidCellPolygons = meshKernelState[meshKernelId].m_mesh2d->ComputeInnerBoundaryPolygons(); meshKernelUndoStack.Add(std::move(undoAction), meshKernelId); } catch (...) @@ -663,6 +677,13 @@ namespace meshkernelapi meshkernel::SplitRowColumnOfMesh splitAlongRow; meshKernelUndoStack.Add(splitAlongRow.Compute(*meshKernelState[meshKernelId].m_mesh2d, edgeId), meshKernelId); + + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } } catch (...) { @@ -881,7 +902,20 @@ namespace meshkernelapi throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); } - meshKernelState[meshKernelId].m_mesh2d->Administrate(); + bool doAdministration = meshKernelState[meshKernelId].m_mesh2d->AdministrationRequired(); + + if (doAdministration) + { + meshKernelState[meshKernelId].m_mesh2d->Administrate(); + } + + if (doAdministration && meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } + SetMesh2dApiDimensions(*meshKernelState[meshKernelId].m_mesh2d, mesh2d); } catch (...) @@ -1298,6 +1332,13 @@ namespace meshkernelapi auto undoAction = MKStateUndoAction::Create(meshKernelState[meshKernelId]); meshkernel::Mesh2DToCurvilinear mesh2DToCurvilinear(*meshKernelState[meshKernelId].m_mesh2d); + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } + meshKernelState[meshKernelId].m_curvilinearGrid = mesh2DToCurvilinear.Compute({xPointCoordinate, yPointCoordinate}); meshKernelUndoStack.Add(std::move(undoAction), meshKernelId); @@ -1325,7 +1366,20 @@ namespace meshkernelapi throw meshkernel::MeshKernelError("Polygon Hanging edge has already been cached. Cached values will be delelted."); } - meshKernelState[meshKernelId].m_mesh2d->Administrate(); + bool doAdministration = meshKernelState[meshKernelId].m_mesh2d->AdministrationRequired(); + + if (doAdministration) + { + meshKernelState[meshKernelId].m_mesh2d->Administrate(); + } + + if (doAdministration && meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } + const auto hangingEdges = meshKernelState[meshKernelId].m_mesh2d->GetHangingEdges(); meshKernelState[meshKernelId].m_hangingEdgeCache = std::make_shared(hangingEdges); numHangingEdges = meshKernelState[meshKernelId].m_hangingEdgeCache->Size(); @@ -1381,6 +1435,29 @@ namespace meshkernelapi return lastExitCode; } + MKERNEL_API int mkernel_mesh2d_delete_faces_in_polygons(int meshKernelId, const GeometryList& polygon) + { + lastExitCode = meshkernel::ExitCode::Success; + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + const std::vector polygonPoints = ConvertGeometryListToPointVector(polygon); + const meshkernel::Polygons invalidCellsPolygon(polygonPoints, meshKernelState[meshKernelId].m_mesh2d->m_projection); + + meshKernelUndoStack.Add(meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidCellsPolygon), meshKernelId); + *meshKernelState[meshKernelId].m_invalidCellPolygons = meshKernelState[meshKernelId].m_mesh2d->ComputeInnerBoundaryPolygons(); + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + MKERNEL_API int mkernel_mesh2d_compute_orthogonalization(int meshKernelId, int projectToLandBoundaryOption, const meshkernel::OrthogonalizationParameters& orthogonalizationParameters, @@ -1416,6 +1493,13 @@ namespace meshkernelapi orthogonalizationParameters); meshKernelUndoStack.Add(ortogonalization.Initialize(), meshKernelId); ortogonalization.Compute(); + + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } } catch (...) { @@ -2046,6 +2130,14 @@ namespace meshkernelapi const meshkernel::Polygons polygons(boundaryPolygonPoints, meshKernelState[meshKernelId].m_projection); meshkernel::Mesh2DIntersections mesh2DIntersections(*meshKernelState[meshKernelId].m_mesh2d); + + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } + mesh2DIntersections.Compute(polygons); auto edgeIntersections = mesh2DIntersections.EdgeIntersections(); auto faceIntersections = mesh2DIntersections.FaceIntersections(); @@ -2635,7 +2727,21 @@ namespace meshkernelapi { throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); } - meshKernelState[meshKernelId].m_mesh2d->Administrate(); + + bool doAdministration = meshKernelState[meshKernelId].m_mesh2d->AdministrationRequired(); + + if (doAdministration) + { + meshKernelState[meshKernelId].m_mesh2d->Administrate(); + } + + if (doAdministration && meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } + const auto numFaces = meshKernelState[meshKernelId].m_mesh2d->GetNumFaces(); std::vector validFace(numFaces, false); for (meshkernel::UInt f = 0; f < numFaces; ++f) @@ -2665,7 +2771,21 @@ namespace meshkernelapi { throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); } - meshKernelState[meshKernelId].m_mesh2d->Administrate(); + + bool doAdministration = meshKernelState[meshKernelId].m_mesh2d->AdministrationRequired(); + + if (doAdministration) + { + meshKernelState[meshKernelId].m_mesh2d->Administrate(); + } + + if (doAdministration && meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } + const auto numFaces = meshKernelState[meshKernelId].m_mesh2d->GetNumFaces(); int numMatchingFaces = 0; for (meshkernel::UInt f = 0; f < numFaces; ++f) @@ -2899,6 +3019,13 @@ namespace meshkernelapi std::move(averaging), meshRefinementParameters); meshKernelUndoStack.Add(meshRefinement.Compute(), meshKernelId); + + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } } catch (...) { @@ -2957,6 +3084,13 @@ namespace meshkernelapi std::move(averaging), meshRefinementParameters); meshKernelUndoStack.Add(meshRefinement.Compute(), meshKernelId); + + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } } catch (...) { @@ -2995,6 +3129,13 @@ namespace meshkernelapi meshRefinementParameters, useNodalRefinement); meshKernelUndoStack.Add(meshRefinement.Compute(), meshKernelId); + + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } } catch (...) { @@ -3025,6 +3166,16 @@ namespace meshkernelapi meshkernel::MeshRefinement meshRefinement(*meshKernelState[meshKernelId].m_mesh2d, polygon, meshRefinementParameters); meshKernelUndoStack.Add(meshRefinement.Compute(), meshKernelId); + + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } + + meshkernel::SaveVtk(meshKernelState[meshKernelId].m_mesh2d->Nodes(), meshKernelState[meshKernelId].m_mesh2d->m_facesNodes, "mesh2.vtu"); + meshkernel::Print(meshKernelState[meshKernelId].m_mesh2d->Nodes(), meshKernelState[meshKernelId].m_mesh2d->Edges()); } catch (...) { @@ -3049,6 +3200,13 @@ namespace meshkernelapi meshkernel::RemoveDisconnectedRegions removeDisconnectedRegions; meshKernelUndoStack.Add(removeDisconnectedRegions.Compute(*meshKernelState[meshKernelId].m_mesh2d), meshKernelId); + + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } } catch (...) { @@ -3220,6 +3378,13 @@ namespace meshkernelapi } meshKernelUndoStack.Add(meshkernel::CasulliRefinement::Compute(*meshKernelState[meshKernelId].m_mesh2d), meshKernelId); + + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } } catch (...) { @@ -3247,6 +3412,13 @@ namespace meshkernelapi meshKernelUndoStack.Add(meshkernel::CasulliRefinement::Compute(*meshKernelState[meshKernelId].m_mesh2d, meshKernelPolygons), meshKernelId); + + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } } catch (...) { @@ -3293,6 +3465,13 @@ namespace meshkernelapi meshRefinementParameters, minimumRefinementDepth); + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } + meshKernelUndoStack.Add(std::move(undoAction), meshKernelId); } } @@ -3388,6 +3567,13 @@ namespace meshkernelapi } meshKernelUndoStack.Add(meshkernel::CasulliDeRefinement::Compute(*meshKernelState[meshKernelId].m_mesh2d), meshKernelId); + + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } } catch (...) { @@ -3414,6 +3600,13 @@ namespace meshkernelapi meshKernelUndoStack.Add(meshkernel::CasulliDeRefinement::Compute(*meshKernelState[meshKernelId].m_mesh2d, meshKernelPolygons), meshKernelId); + + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } } catch (...) { @@ -4052,7 +4245,23 @@ namespace meshkernelapi meshKernelState[meshKernelId].m_mesh2d->SetNodes(mergedMeshes->Nodes()); meshKernelState[meshKernelId].m_mesh2d->SetEdges(mergedMeshes->Edges()); - meshKernelState[meshKernelId].m_mesh2d->Administrate(); + + bool doAdministration = meshKernelState[meshKernelId].m_mesh2d->AdministrationRequired(); + + if (doAdministration) + { + meshKernelState[meshKernelId].m_mesh2d->Administrate(); + } + + if (doAdministration && meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } + + // Update the missing elements polygons + *meshKernelState[meshKernelId].m_invalidCellPolygons = meshKernelState[meshKernelId].m_mesh2d->ComputeInnerBoundaryPolygons(); meshKernelUndoStack.Add(std::move(undoAction), meshKernelId); } catch (...) @@ -5164,7 +5373,7 @@ namespace meshkernelapi // The undo action for conversion of clg to m2d is made in two steps std::unique_ptr undoAction = meshkernel::CompoundUndoAction::Create(); - // 1. Keep the mesh kernel state to be able to restore (manily) the curvilinear grid and (secondly) other members. + // 1. Keep the mesh kernel state to be able to restore (firstly) the curvilinear grid and (secondly) other members. undoAction->Add(MKStateUndoAction::Create(meshKernelState[meshKernelId])); // 2. Keep track of the undo required to restore the mesh2d to its pre-converted state. @@ -5173,6 +5382,13 @@ namespace meshkernelapi // curvilinear grid must be reset to an empty curvilinear grid meshKernelState[meshKernelId].m_curvilinearGrid = std::make_unique(); meshKernelUndoStack.Add(std::move(undoAction), meshKernelId); + + if (meshKernelState[meshKernelId].m_invalidCellPolygons->size() > 0) + { + meshkernel::Polygons invalidElementPolygon(*meshKernelState[meshKernelId].m_invalidCellPolygons, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + meshKernelState[meshKernelId].m_mesh2d->DeleteMeshFacesInPolygons(invalidElementPolygon); + } } catch (...) { From 5ac8946acb83b53199ec01e0e66085c65589c401 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 9 Oct 2025 09:39:43 +0200 Subject: [PATCH 10/16] GRIDEDIT-1795 Added a new file containing tests for the keeping of invalid cells upto date --- libs/MeshKernelApi/tests/CMakeLists.txt | 1 + .../tests/src/InvalidCellsPolygonsTests.cpp | 182 ++++++++++++++++++ .../tests/src/Mesh2DRefinmentTests.cpp | 2 - 3 files changed, 183 insertions(+), 2 deletions(-) create mode 100644 libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp diff --git a/libs/MeshKernelApi/tests/CMakeLists.txt b/libs/MeshKernelApi/tests/CMakeLists.txt index 99f713348c..ab1e61d99b 100644 --- a/libs/MeshKernelApi/tests/CMakeLists.txt +++ b/libs/MeshKernelApi/tests/CMakeLists.txt @@ -17,6 +17,7 @@ set(SRC_LIST ${SRC_DIR}/CurvilinearGridTests.cpp ${SRC_DIR}/CurvilinearGridUndoTests.cpp ${SRC_DIR}/ErrorHandlingTests.cpp + ${SRC_DIR}/InvalidCellsPolygonsTests.cpp ${SRC_DIR}/LandBoundaryTests.cpp ${SRC_DIR}/Mesh1DTests.cpp ${SRC_DIR}/Mesh2DCasulliRefinmentTests.cpp diff --git a/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp b/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp new file mode 100644 index 0000000000..d477725807 --- /dev/null +++ b/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp @@ -0,0 +1,182 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2025. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation version 3. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include "CartesianApiTestFixture.hpp" +#include "MeshKernel/AveragingInterpolation.hpp" +#include "MeshKernel/MeshRefinement.hpp" +#include "SampleGenerator.hpp" + +#include "MeshKernel/Utilities/Utilities.hpp" + +TEST(InvalidCellsPolygonsTests, OMG) +{ + // Prepare + int meshKernelId; + constexpr int isSpherical = 0; + meshkernelapi::mkernel_allocate_state(isSpherical, meshKernelId); + + meshkernel::MakeGridParameters gridParameters; + gridParameters.num_columns = 10; + gridParameters.num_rows = 10; + gridParameters.block_size_x = 20.0; + gridParameters.block_size_y = 20.0; + gridParameters.origin_x = 0.0; + gridParameters.origin_y = 0.0; + gridParameters.angle = 0.0; + + auto errorCode = meshkernelapi::mkernel_mesh2d_make_rectangular_mesh(meshKernelId, gridParameters); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + //-------------------------------- + + std::vector patchX{45.0, 155.0, 155.0, 45.0, 45.0, meshkernel::constants::missing::innerOuterSeparator, 65.0, 115.0, 115.0, 65.0, 65.0}; + std::vector patchY{5.0, 5.0, 155.0, 155.0, 5.0, meshkernel::constants::missing::innerOuterSeparator, 65.0, 65.0, 115.0, 115.0, 65.0}; + + meshkernelapi::GeometryList refinementPolygon; + refinementPolygon.num_coordinates = static_cast(patchX.size()); + refinementPolygon.geometry_separator = meshkernel::constants::missing::doubleValue; + refinementPolygon.coordinates_x = patchX.data(); + refinementPolygon.coordinates_y = patchY.data(); + + errorCode = meshkernelapi::mkernel_mesh2d_casulli_refinement_on_polygon(meshKernelId, refinementPolygon); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + //-------------------------------- + + std::vector invalidCellsX{80.0, 85.0, 75.0, 80.0, meshkernel::constants::missing::doubleValue, + 180.0, 200.0, 200.0, 180.0, 180.0, meshkernel::constants::missing::doubleValue, + 125.0, 135.0, 135.0, 125.0, 125.0, meshkernel::constants::missing::doubleValue, + 85.0, 95.0, 95.0, 85.0, 85.0, meshkernel::constants::missing::doubleValue, + 100.0, 105.0, 95.0, 100.0, meshkernel::constants::missing::doubleValue, + 120.0, 125.0, 115.0, 120.0}; + + std::vector invalidCellsY{0.0, 15.0, 15.0, 0.0, meshkernel::constants::missing::doubleValue, + 140.0, 140.0, 160.0, 160.0, 140.0, meshkernel::constants::missing::doubleValue, + 125.0, 125.0, 135.0, 135.0, 125.0, meshkernel::constants::missing::doubleValue, + 15.0, 15.0, 25.0, 25.0, 15.0, meshkernel::constants::missing::doubleValue, + 0.0, 15.0, 15.0, 0.0, meshkernel::constants::missing::doubleValue, + 0.0, 15.0, 15.0, 0.0}; + + meshkernelapi::GeometryList invalidCells; + invalidCells.num_coordinates = static_cast(invalidCellsX.size()); + invalidCells.geometry_separator = meshkernel::constants::missing::doubleValue; + invalidCells.coordinates_x = invalidCellsX.data(); + invalidCells.coordinates_y = invalidCellsY.data(); + + errorCode = meshkernelapi::mkernel_mesh2d_delete_faces_in_polygons(meshKernelId, invalidCells); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + meshkernelapi::GeometryList geometryListIn; + geometryListIn.geometry_separator = meshkernel::constants::missing::doubleValue; + // geometryListIn.coordinates_x = xCoordinatesIn.data(); + // geometryListIn.coordinates_y = yCoordinatesIn.data(); + // geometryListIn.values = valuesIn.data(); + geometryListIn.num_coordinates = 0; // static_cast(xCoordinatesIn.size()); + + meshkernel::MeshRefinementParameters meshRefinementParameters; + meshRefinementParameters.max_num_refinement_iterations = 2; + // meshRefinementParameters.min_edge_size = 200000.0; + + const int isTriangulationRequired = 1; + const int projectToLandBoundaryOption = 1; + meshkernelapi::GeometryList selectingPolygon{}; + meshkernelapi::GeometryList landBoundaries{}; + errorCode = mkernel_mesh2d_flip_edges(meshKernelId, + isTriangulationRequired, + projectToLandBoundaryOption, + selectingPolygon, + landBoundaries); + + // Execute + errorCode = mkernel_mesh2d_refine_based_on_polygon(meshKernelId, geometryListIn, meshRefinementParameters); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + // errorCode = meshkernelapi::mkernel_mesh2d_casulli_refinement(meshKernelId); + // ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + // // Get the mesh dimensions + // meshkernelapi::Mesh2D mesh2d{}; + // errorCode = meshkernelapi::mkernel_mesh2d_get_dimensions(meshKernelId, mesh2d); + // ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + // ASSERT_EQ(mesh2d.num_nodes, static_cast(expectedX.size())); + // ASSERT_EQ(2 * mesh2d.num_edges, static_cast(expectedEdgesNodes.size())); + + // std::vector edge_faces(mesh2d.num_edges * 2); + // std::vector edge_nodes(mesh2d.num_edges * 2); + // std::vector face_nodes(mesh2d.num_face_nodes); + // std::vector face_edges(mesh2d.num_face_nodes); + // std::vector nodes_per_face(mesh2d.num_faces); + // std::vector node_x(mesh2d.num_nodes); + // std::vector node_y(mesh2d.num_nodes); + // std::vector edge_x(mesh2d.num_edges); + // std::vector edge_y(mesh2d.num_edges); + // std::vector face_x(mesh2d.num_faces); + // std::vector face_y(mesh2d.num_faces); + + // mesh2d.edge_faces = edge_faces.data(); + // mesh2d.edge_nodes = edge_nodes.data(); + // mesh2d.face_nodes = face_nodes.data(); + // mesh2d.face_edges = face_edges.data(); + // mesh2d.nodes_per_face = nodes_per_face.data(); + // mesh2d.node_x = node_x.data(); + // mesh2d.node_y = node_y.data(); + // mesh2d.edge_x = edge_x.data(); + // mesh2d.edge_y = edge_y.data(); + // mesh2d.face_x = face_x.data(); + // mesh2d.face_y = face_y.data(); + // errorCode = mkernel_mesh2d_get_data(meshKernelId, mesh2d); + // ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + // const double tolerance = 1e-12; + + // for (int i = 0; i < mesh2d.num_nodes; ++i) + // { + // EXPECT_NEAR(expectedX[i], node_x[i], tolerance); + // EXPECT_NEAR(expectedY[i], node_y[i], tolerance); + // } + + // for (int i = 0; i < 2 * mesh2d.num_edges; ++i) + // { + // EXPECT_EQ(expectedEdgesNodes[i], edge_nodes[i]); + // } +} diff --git a/libs/MeshKernelApi/tests/src/Mesh2DRefinmentTests.cpp b/libs/MeshKernelApi/tests/src/Mesh2DRefinmentTests.cpp index ad145f84fe..74f988489b 100644 --- a/libs/MeshKernelApi/tests/src/Mesh2DRefinmentTests.cpp +++ b/libs/MeshKernelApi/tests/src/Mesh2DRefinmentTests.cpp @@ -45,8 +45,6 @@ #include "MeshKernel/MeshRefinement.hpp" #include "SampleGenerator.hpp" -#include "MeshKernel/Utilities/Utilities.hpp" - namespace fs = std::filesystem; TEST_F(CartesianApiTestFixture, RefineAPolygonThroughApi) From 86d1ca7654ed4408ed78f66466e070789ff27cd5 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 9 Oct 2025 09:53:28 +0200 Subject: [PATCH 11/16] GRIDEDIT-1795 Merged master into the branch --- libs/MeshKernel/tests/src/MeshRefinementTests.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index 315c405fb5..44a6a7871b 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -3133,7 +3133,8 @@ TEST(MeshRefinement, WTF) for (size_t i = 0; i < boundaryNodes2.size(); ++i) { - std::cout << "point " << count << " " << " = " << boundaryNodes2[i].x << ", " << boundaryNodes2[i].y << std::endl; + std::cout << "point " << count << " " << " = " + << boundaryNodes2[i].x << ", " << boundaryNodes2[i].y << std::endl; ++count; } From 86da75a1f055ec3cb2268a798758b566b836d10f Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 9 Oct 2025 10:02:49 +0200 Subject: [PATCH 12/16] GRIDEDIT-1795 Fix windows build, update comments and fix clang formatting warning --- libs/MeshKernel/src/Mesh2D.cpp | 8 ++++---- libs/MeshKernel/tests/src/MeshRefinementTests.cpp | 11 ----------- .../include/MeshKernelApi/MeshKernel.hpp | 1 + 3 files changed, 5 insertions(+), 15 deletions(-) diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index ffb9ed1d4d..90bdeba3f5 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -1656,7 +1656,7 @@ void Mesh2D::OrientatePolygonsAntiClockwise(std::vector& polygonNodes) co { polygonStart = index; - for (size_t i = polygonStart; i < polygonNodes.size(); ++i) + for (UInt i = polygonStart; i < polygonNodes.size(); ++i) { ++index; @@ -1684,12 +1684,12 @@ void Mesh2D::OrientatePolygonsAntiClockwise(std::vector& polygonNodes) co // reverse order of polygon nodes if (polygonLength - 1 == 3) { - // Only the second and third points need be swapped to reverse the points in the polygon + // Only the second and third points need be swapped to reverse the points in a triangle polygon std::swap(polygonNodes[polygonStart + 1], polygonNodes[polygonStart + 2]); } else if (polygonLength - 1 == 4) { - // Only the second and fourth points need be swapped to reverse the points in the polygon + // Only the second and fourth points need be swapped to reverse the points in a quadrialteral polygon std::swap(polygonNodes[polygonStart + 1], polygonNodes[polygonStart + 3]); } else @@ -1822,7 +1822,7 @@ void Mesh2D::WalkMultiBoundaryFromNode(std::vector& edgeIsVisited, // Find index of last time node was added for (size_t ii = nodeIds.size(); ii >= 1; --ii) { - size_t i = ii - 1; + UInt i = static_cast(ii) - 1; if (nodeIds[i] == constants::missing::uintValue) { diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index 44a6a7871b..7f1fd1fe6a 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -3127,17 +3127,6 @@ TEST(MeshRefinement, WTF) // boundaryNodes3.insert(boundaryNodes3.end(), boundaryNodes2.begin() + polygonStart + polygonLength, boundaryNodes2.end()); // } - int count = 0; - - std::cout << "boundaryNodes2.size " << boundaryNodes2.size() << std::endl; - - for (size_t i = 0; i < boundaryNodes2.size(); ++i) - { - std::cout << "point " << count << " " << " = " - << boundaryNodes2[i].x << ", " << boundaryNodes2[i].y << std::endl; - ++count; - } - [[maybe_unused]] Polygons boundaryWithMissingElement2(boundaryNodes2, Projection::cartesian); // for (size_t i = 0; i < boundaryNodes3.size(); ++i) diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp index 3c8cf1fe5b..495ac9e609 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp @@ -1152,6 +1152,7 @@ namespace meshkernelapi /// @brief Deletes /// @param[in] meshKernelId The id of the mesh state + /// @param[in] polygon The polygon regions where the faces are to be deleted /// @returns Error code MKERNEL_API int mkernel_mesh2d_delete_faces_in_polygons(int meshKernelId, const GeometryList& polygon); From 5f14efd4d11e9a22afe56820597270bc87d50aa2 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 9 Oct 2025 10:21:15 +0200 Subject: [PATCH 13/16] GRIDEDIT-1795 Updated comments --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 1 - libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 0536529b53..d549dea139 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -245,7 +245,6 @@ namespace meshkernel [[nodiscard]] std::vector ComputeBoundaryPolygons(const std::vector& polygon); /// @brief Convert all mesh boundaries to a vector of polygon nodes - /// @param[in] polygon The polygon where the operation is performed /// @return The resulting set of polygons, describing interior mesh boundaries std::vector ComputeInnerBoundaryPolygons() const; diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp index 495ac9e609..382f53db3b 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp @@ -1150,7 +1150,7 @@ namespace meshkernelapi /// @returns Error code MKERNEL_API int mkernel_mesh2d_delete_hanging_edges(int meshKernelId); - /// @brief Deletes + /// @brief Deletes the mesh faces inside a set of polygons /// @param[in] meshKernelId The id of the mesh state /// @param[in] polygon The polygon regions where the faces are to be deleted /// @returns Error code From 3ef4478ed4296a9b8e80310905cd543242d99b1d Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 9 Oct 2025 11:23:36 +0200 Subject: [PATCH 14/16] GRIDEDIT-1795 Updated comments --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index d549dea139..b703c5bf07 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -449,6 +449,7 @@ namespace meshkernel std::vector& nodeIds, std::vector& subsSequence) const; + /// @brief Ensure that all polynomials are orientated in the ACW direction. void OrientatePolygonsAntiClockwise(std::vector& polygonNodes) const; /// @brief Removes the outer domain boundary polygon from the set of polygons From 7d80ffe5802eaf2e0ba72438d2b779d13b10da50 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 9 Oct 2025 16:52:50 +0200 Subject: [PATCH 15/16] GRIDEDIT-1795 Added extra information to polygon not closed error --- libs/MeshKernel/src/Polygon.cpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/libs/MeshKernel/src/Polygon.cpp b/libs/MeshKernel/src/Polygon.cpp index 1888a87b62..cb9aa80cfb 100644 --- a/libs/MeshKernel/src/Polygon.cpp +++ b/libs/MeshKernel/src/Polygon.cpp @@ -58,7 +58,15 @@ void meshkernel::Polygon::Initialise() if (m_nodes.size() > 0 && m_nodes[0] != m_nodes[m_nodes.size() - 1]) { - throw ConstraintError("Polygon is not closed"); + std::stringstream buffer; + + buffer << "Polygon is not closed" << std::endl; + buffer.precision(15); + buffer << "first point " << m_nodes[0].x << ", " << m_nodes[0].y << std::endl; + buffer << "last point " << m_nodes[m_nodes.size() - 1].x << ", " << m_nodes[m_nodes.size() - 1].y << std::endl; + + throw ConstraintError(buffer.str()); + // throw ConstraintError("Polygon is not closed"); } if (InvalidPointCount(m_nodes) > 0) From 94b84ffcd72f40ba84c7762c25423e8a7213be27 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 9 Oct 2025 16:53:46 +0200 Subject: [PATCH 16/16] GRIDEDIT-1795 Added extra information to polygon not closed error --- libs/MeshKernel/src/Polygon.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/libs/MeshKernel/src/Polygon.cpp b/libs/MeshKernel/src/Polygon.cpp index cb9aa80cfb..faadeed91a 100644 --- a/libs/MeshKernel/src/Polygon.cpp +++ b/libs/MeshKernel/src/Polygon.cpp @@ -62,6 +62,7 @@ void meshkernel::Polygon::Initialise() buffer << "Polygon is not closed" << std::endl; buffer.precision(15); + buffer << "size " << m_nodes.size() << std::endl; buffer << "first point " << m_nodes[0].x << ", " << m_nodes[0].y << std::endl; buffer << "last point " << m_nodes[m_nodes.size() - 1].x << ", " << m_nodes[m_nodes.size() - 1].y << std::endl;