diff --git a/libs/MeshKernel/include/MeshKernel/Entities.hpp b/libs/MeshKernel/include/MeshKernel/Entities.hpp index 5d665218c..6d539dfed 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) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh.hpp b/libs/MeshKernel/include/MeshKernel/Mesh.hpp index f330add6a..e7a1b1967 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 7551be5e6..b703c5bf0 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -244,15 +244,9 @@ 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 + /// @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 @@ -378,6 +372,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; + /// @brief Deletes the mesh faces inside a set of polygons + /// @param[in] polygon The polygon where to perform the operation + /* [[nodiscard]] */ std::unique_ptr DeleteMeshFacesInPolygons(const Polygons& polygon); + private: // orthogonalization static constexpr double m_minimumEdgeLength = 1e-4; ///< Minimum edge length @@ -435,6 +433,30 @@ 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 individual polygons + void WalkMultiBoundaryFromNode(std::vector& edgeIsVisited, + std::vector& nodeIsVisited, + UInt& currentNode, + std::vector& meshBoundaryPolygon, + 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 + /// + /// 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. @@ -461,6 +483,12 @@ namespace meshkernel std::vector& sortedNodes, 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) /// @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 e92a23c76..90bdeba3f 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -25,11 +25,14 @@ // //------------------------------------------------------------------------------ -#include "MeshKernel/Mesh2D.hpp" +#include +#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" @@ -582,6 +585,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(); @@ -1552,6 +1561,203 @@ std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector return meshBoundaryPolygon; } +std::vector Mesh2D::ComputeInnerBoundaryPolygons() const +{ + if (GetNumFaces() == 0) + { + return std::vector(); + } + + 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; + WalkMultiBoundaryFromNode(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; + WalkMultiBoundaryFromNode(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; + } + } + } + + 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 (UInt 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 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 a quadrialteral 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 +{ + // Remove outer boundary. + // It is assumed that the boundary polygon with the most points is from the outer domain boundary + + UInt outerPolygonLength = 0; + UInt outerPolygonStart = 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()) + { + outerPolygonStartIntermediate = index; + + for (UInt i = outerPolygonStartIntermediate; i < polygonNodes.size(); ++i) + { + ++index; + + if (!polygonNodes[i].IsValid()) + { + outerPolygonLengthIntermediate = i - outerPolygonStartIntermediate; + break; + } + } + + if (outerPolygonLengthIntermediate > outerPolygonLength) + { + outerPolygonLength = outerPolygonLengthIntermediate; + outerPolygonStart = outerPolygonStartIntermediate; + } + } + + if (outerPolygonLength > 0 && outerPolygonLength < polygonNodes.size()) + { + ++outerPolygonLength; + } + + std::vector innerBoundaryNodes; + + // 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() + outerPolygonStart - 1); + } + + if (outerPolygonLength > 0 && outerPolygonStart + outerPolygonLength < polygonNodes.size()) + { + innerBoundaryNodes.insert(innerBoundaryNodes.end(), polygonNodes.begin() + outerPolygonStart + outerPolygonLength, polygonNodes.end()); + } + + return innerBoundaryNodes; +} + void Mesh2D::WalkBoundaryFromNode(const Polygon& polygon, std::vector& isVisited, UInt& currentNode, @@ -1587,6 +1793,76 @@ void Mesh2D::WalkBoundaryFromNode(const Polygon& polygon, } } +void Mesh2D::WalkMultiBoundaryFromNode(std::vector& edgeIsVisited, + std::vector& nodeIsVisited, + UInt& currentNode, + std::vector& meshBoundaryPolygon, + std::vector& nodeIds, + std::vector& subSequence) const +{ + UInt e = 0; + + 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) + { + UInt i = static_cast(ii) - 1; + + if (nodeIds[i] == constants::missing::uintValue) + { + break; + } + + if (nodeIds[i] == nextNode) + { + lastIndex = i; + break; + } + } + + if (lastIndex != constants::missing::uintValue) + { + + if (!subSequence.empty()) + { + subSequence.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); + } + + for (size_t i = lastIndex; i < nodeIds.size(); ++i) + { + subSequence.emplace_back(meshBoundaryPolygon[i]); + } + + 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); + } +} + std::vector Mesh2D::GetHangingEdges() const { std::vector result; @@ -1870,6 +2146,106 @@ std::unique_ptr Mesh2D::DeleteMeshFaces(const Polygons& return deleteMeshAction; } +std::unique_ptr Mesh2D::DeleteMeshFacesInPolygons(const Polygons& polygon) +{ + if (polygon.IsEmpty()) + { + return nullptr; + } + + 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() && !polygon.IsPointInAnyPolygon(m_facesMassCenters[f])) + { + faceIndices[f] = faceIndex; + ++faceIndex; + } + } + + UpdateFaceInformation(faceIndices, *deleteMeshAction); + + return deleteMeshAction; +} + +void Mesh2D::UpdateFaceInformation(const std::vector& faceIndices, CompoundUndoAction& deleteMeshAction) +{ + 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) + { + UInt edge = m_facesEdges[faceId][e]; + bool removedFace = false; + + if (m_edgesFaces[edge][0] == faceId) + { + 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]); + removedFace = true; + } + else if (m_edgesFaces[edge][1] == faceId) + { + m_edgesFaces[edge][1] = constants::missing::uintValue; + --m_edgesNumFaces[edge]; + removedFace = true; + } + + if (removedFace && m_edgesNumFaces[edge] == 0) + { + deleteMeshAction.Add(DeleteEdge(edge)); + } + } + } + + // 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? + // 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); + m_facesEdges.erase(m_facesEdges.begin() + faceId); + m_facesMassCenters.erase(m_facesMassCenters.begin() + faceId); + m_faceArea.erase(m_faceArea.begin() + faceId); + } +} + std::unique_ptr Mesh2D::DeleteHangingEdges() { std::unique_ptr deleteAction = CompoundUndoAction::Create(); diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index 73f5725c5..05865d5bb 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -814,6 +814,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/src/Polygon.cpp b/libs/MeshKernel/src/Polygon.cpp index 1888a87b6..faadeed91 100644 --- a/libs/MeshKernel/src/Polygon.cpp +++ b/libs/MeshKernel/src/Polygon.cpp @@ -58,7 +58,16 @@ 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 << "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; + + throw ConstraintError(buffer.str()); + // throw ConstraintError("Polygon is not closed"); } if (InvalidPointCount(m_nodes) > 0) diff --git a/libs/MeshKernel/src/Polygons.cpp b/libs/MeshKernel/src/Polygons.cpp index ed58f5074..de98fa1d2 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 60039ae6c..7f1fd1fe6 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" @@ -2715,7 +2718,6 @@ TEST(MeshRefinement, RowSplittingFailureTests) // Invalid edge id EXPECT_THROW([[maybe_unused]] auto undo4 = splitMeshRow.Compute(mesh, edgeId), ConstraintError); } - TEST(MeshRefinement, MeshRefinementInsidePolygon_On21x21With50x50Samples_ShouldRefineMesh) { const meshkernel::UInt size = 21; @@ -2777,3 +2779,517 @@ TEST(MeshRefinement, MeshRefinementInsidePolygon_On21x21With50x50Samples_ShouldR ASSERT_EQ(mesh->GetNumEdges(), 1240); ASSERT_EQ(mesh->GetNumFaces(), 630); } + +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(); + + 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 boundaryNodes; + + 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::doubleValue, constants::missing::doubleValue}, + mesh.Node(node21), + mesh.Node(node22), + mesh.Node(node23), + mesh.Node(node24), + mesh.Node(node21)}; + + 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::doubleValue, constants::missing::doubleValue}, + mesh.Node(node41), + mesh.Node(node42), + mesh.Node(node43), + mesh.Node(node44), + mesh.Node(node41)}; + + 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::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()); + 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 boundaryWithMissingElements(boundaryNodes, Projection::cartesian); + + auto deleteMeshFacesUndoAction = mesh.DeleteMeshFacesInPolygons(boundaryWithMissingElements); + + // Compute interior boundary polygon points + std::vector boundaryNodes2 = mesh.ComputeInnerBoundaryPolygons(); + + UInt expectedNumberOfNodes = 26; + + 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, + 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; + 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); + + 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 boundaryNodes; + + std::vector elementNodes1{mesh.Node(node11), + mesh.Node(node12), + mesh.Node(node13), + mesh.Node(node11)}; + + 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::doubleValue, constants::missing::doubleValue}, + mesh.Node(node31), + mesh.Node(node32), + mesh.Node(node33), + mesh.Node(node34), + mesh.Node(node31)}; + + 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::doubleValue, constants::missing::doubleValue}, + mesh.Node(node51), + mesh.Node(node52), + mesh.Node(node53), + mesh.Node(node51)}; + + 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()); + 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 boundaryWithMissingElements(boundaryNodes, Projection::cartesian); + auto undoAction2 = mesh.DeleteMeshFacesInPolygons(boundaryWithMissingElements); + meshkernel::SaveVtk(mesh.Nodes(), mesh.m_facesNodes, "mesh1.vtu"); + + //-------------------------------- + + 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()); + // } + + [[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, "mesh2.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); + + 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 boundaryNodes; + + 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::doubleValue, constants::missing::doubleValue}, + mesh.Node(node21), + mesh.Node(node22), + mesh.Node(node23), + mesh.Node(node24), + mesh.Node(node21)}; + + 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::doubleValue, constants::missing::doubleValue}, + mesh.Node(node41), + mesh.Node(node42), + mesh.Node(node43), + mesh.Node(node44), + mesh.Node(node41)}; + + 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::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()); + 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, "mesh1.vtu"); +} diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp index cb41751f4..382f53db3 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp @@ -1150,6 +1150,12 @@ namespace meshkernelapi /// @returns Error code MKERNEL_API int mkernel_mesh2d_delete_hanging_edges(int meshKernelId); + /// @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 + 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 45ca15434..744c7c9ff 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 a233c3467..42cd83c2b 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 f2321c645..3bd7e7215 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 (...) { diff --git a/libs/MeshKernelApi/tests/CMakeLists.txt b/libs/MeshKernelApi/tests/CMakeLists.txt index 99f713348..ab1e61d99 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 000000000..d47772580 --- /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 ad145f84f..74f988489 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)