diff --git a/libs/MeshKernel/CMakeLists.txt b/libs/MeshKernel/CMakeLists.txt index 2fb3902ad..e2dd97211 100644 --- a/libs/MeshKernel/CMakeLists.txt +++ b/libs/MeshKernel/CMakeLists.txt @@ -45,6 +45,7 @@ set( ${SRC_DIR}/LandBoundaries.cpp ${SRC_DIR}/LandBoundary.cpp ${SRC_DIR}/Mesh.cpp + ${SRC_DIR}/MeshBoundaryExtractor.cpp ${SRC_DIR}/MeshEdgeLength.cpp ${SRC_DIR}/Mesh1D.cpp ${SRC_DIR}/Mesh2D.cpp @@ -181,6 +182,7 @@ set( ${DOMAIN_INC_DIR}/LandBoundaries.hpp ${DOMAIN_INC_DIR}/LandBoundary.hpp ${DOMAIN_INC_DIR}/Mesh.hpp + ${DOMAIN_INC_DIR}/MeshBoundaryExtractor.hpp ${DOMAIN_INC_DIR}/MeshEdgeLength.hpp ${DOMAIN_INC_DIR}/Mesh1D.hpp ${DOMAIN_INC_DIR}/Mesh2D.hpp diff --git a/libs/MeshKernel/include/MeshKernel/Definitions.hpp b/libs/MeshKernel/include/MeshKernel/Definitions.hpp index c02d6c250..8cc1ba3c3 100644 --- a/libs/MeshKernel/include/MeshKernel/Definitions.hpp +++ b/libs/MeshKernel/include/MeshKernel/Definitions.hpp @@ -174,4 +174,12 @@ namespace meshkernel AllNetlinksLoop }; + /// \brief Boundary selection indicator + enum class BoundarySelection + { + ExteriorOnly, ///< Exterior boundaries only + InteriorOnly, ///< Interior boundaries only + All ///< All boundaries, exterior and interior + }; + } // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index f1f163975..ae7250e58 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -254,10 +254,6 @@ namespace meshkernel /// @return The resulting polygon mesh boundary [[nodiscard]] std::vector ComputeBoundaryPolygons(const std::vector& polygon); - /// @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 [[nodiscard]] std::vector GetHangingEdges() const; @@ -444,25 +440,17 @@ 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 Reconstruct the invalid cell polygons + /// + /// When constructing the invalid cell polygons, they can be computed with many smaller polygons. + /// If these smaller polygons form a single patch on the domain, then they need to be combined + void ReconstructInvalidCellsPolygon(); - /// @brief Constructs a polygon or polygons from the meshboundary, by walking through the mesh + /// @brief Convert all mesh boundaries to a vector of polygon nodes, including holes (copynetboundstopol) /// - /// 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& subSequence, - std::vector& illegalCells) const; - - /// @brief Ensure that all polynomials are orientated in the ACW direction. - void OrientatePolygonsAntiClockwise(std::vector& polygonNodes) const; + /// @return a sequence of boundary points, which may be separated by the invalid point, and a matching sequence of Boolean values indicating which + /// polygonal sub-sequence forms a external boundary + [[nodiscard]] std::tuple, std::vector> GetAllBoundaryPolygons(const std::vector& polygon); /// @brief Removes the outer domain boundary polygon from the set of polygons /// diff --git a/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp new file mode 100644 index 000000000..41fc0c4a4 --- /dev/null +++ b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp @@ -0,0 +1,126 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2026. +// +// 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. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include +#include + +#include "MeshKernel/Definitions.hpp" +#include "MeshKernel/Mesh2D.hpp" +#include "MeshKernel/Point.hpp" +#include "MeshKernel/Polygon.hpp" + +namespace meshkernel +{ + + /// @brief Extract the boundary polygon from the mesh + class MeshBoundaryExtractor + { + public: + /// @brief Extract all boundaries as a single sequence of points, separated by an invalid point + static std::vector ExtractConcatenated(const Mesh2D& mesh, BoundarySelection boundaryType = BoundarySelection::All); + + /// @brief Extract all boundaries keeping them separated + /// + /// The result consists of an array of each of the boundary polygons + /// Additionally, an array indicating if the boundary polygon is a exterior boundary or not. + /// True => is-exterior, False => otherwise + static std::tuple>, std::vector> Extract(const Mesh2D& mesh); + + /// @brief Extract all boundaries contained within a constraining polygon + /// + /// The result consists of an array of each of the boundary polygons + /// Additionally, an array indicating if the boundary polygon is a exterior boundary or not. + /// True => is-exterior, False => otherwise + static std::tuple>, std::vector> Extract(const Mesh2D& mesh, const Polygon& polygon); + + private: + /// @brief The minimum number of points in a polygon, excluding the closing point + static constexpr size_t MinimumNumberOfPoints = 3; + + /// @brief Temporary struct, used when computing the boundaries + struct BoundaryEdge + { + UInt edgeId; ///< Id of edge + UInt neighbourNode; ///< Id of node at opposite end of edge + UInt leftFace; ///< Store face mapping on the edge for easy retrieval during polygon trace + double angle; ///< Angle of the edge pointing away from the pivot node + }; + + /// @brief Ensure the angle lies between 0 .. 2 pi. + static double NormalizeAngle(double angle); + + /// @brief Find the edge that has the smallest angle [0 .. 2pi), to the incident edge + static UInt FindEdgeWithMinumumAngle(const std::vector& boundaryEdges, + const std::vector& edgeVisited, + const double incomingAngle); + + /// @brief Construct mapping from node-id to all impinging edges for boundary all edges + static void FindAllBoundarEdges(const std::vector& nodes, + const std::vector& edges, + const std::vector>& edgesFaces, + std::unordered_map>& boundaryAdjacency); + + /// @brief Append the boundary polygon to the set of all boundary polygon + /// + /// The boundary polygons may be reversed if found to be in clockwise direction and clipped to be + /// inside a constraining polygon + static void Append(const Polygon& polygon, + const Point& centre, + const Projection projection, + std::vector& boundaryPolygon, + std::vector& isExterior, + std::vector>& separatedBoundaryPolygons); + + /// @brief Clip a boundary polygon node sequence to be contains within a constraining polygon + /// + /// Edges that cross the constraining polygons are also included, i.e. edges that have 1 one node + /// contained with the polygons and another not. + static void ClipToConstrainingPolygon(const Polygon& polygon, std::vector& nodes); + + /// @brief Find boundary loops + /// + /// Any boundary loops found may need to be processed further as they may themselves contain sub-loops + static void FindBoundaryPolygons(const std::vector& nodes, + const std::vector& edges, + const std::vector>& edgesFaces, + std::vector>& allPolygons, + std::vector>& allTouchedFaces); + + /// @brief Separate polygons that contains multiple sub-polygons and determine externality + /// + /// allBoundaryPolygons is not const because it may be updated. + /// @note allBoundaryPolygons should not be accessed after calling this function + static std::tuple>, std::vector> + SeparateAndDetermineExternality(const Mesh2D& mesh, + const Polygon& polygon, + std::vector>& allBoundaryPolygons, + const std::vector>& allTouchedFaces); + }; + +} // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/Operations.hpp b/libs/MeshKernel/include/MeshKernel/Operations.hpp index 8b9f5d07e..57ffdda79 100644 --- a/libs/MeshKernel/include/MeshKernel/Operations.hpp +++ b/libs/MeshKernel/include/MeshKernel/Operations.hpp @@ -298,6 +298,12 @@ namespace meshkernel const std::span triangleNodes, const Projection& projection); + /// \brief Compute the area and the centre of mass of a polygonal area + /// + /// The area of the polygon will be positive is the boundary points traverse the boundary in a anti clockwise direction, + /// If the are is negative then the points traverse the boundary in a clockwise direction. + std::tuple ComputePolygonAreaAndCentre(const std::vector& boundaryPolygon, const Projection& projection); + /// @brief Computes three base components void ComputeThreeBaseComponents(const Point& point, std::array& exxp, std::array& eyyp, std::array& ezzp); @@ -621,4 +627,42 @@ namespace meshkernel } } + /// @brief Determine is a polygon is comprised of multiple polygons intersecting at the control points only + /// + /// @note This currently works only for polygons that intersect at the control points + bool IsMultiPolygon(std::span boundary); + + /// @brief Split a single polygon line comprised of multiple polygons into separate polygons. + /// + /// @note This currently works only for polygons that intersect at the control points + std::tuple>, std::vector> SplitMultiplePolygons(std::span boundary, std::span elementIds); + + /// @brief Concaenate vectors of vectors of points to a vector of points separated by the invalid point + template + std::vector ConcatenatePointVectors(const std::vector>& pointVectors, Predicate predicate) + { + + std::vector combinedPoints; + bool isFirst = true; + + for (size_t count = 0; const auto& points : pointVectors) + { + + if (predicate(count)) + { + if (!isFirst) + { + combinedPoints.push_back({constants::missing::doubleValue, constants::missing::doubleValue}); + } + + combinedPoints.insert(combinedPoints.end(), points.begin(), points.end()); + isFirst = false; + } + + ++count; + } + + return combinedPoints; + } + } // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/Point.hpp b/libs/MeshKernel/include/MeshKernel/Point.hpp index e40775df0..824c38401 100644 --- a/libs/MeshKernel/include/MeshKernel/Point.hpp +++ b/libs/MeshKernel/include/MeshKernel/Point.hpp @@ -125,6 +125,14 @@ namespace meshkernel return !isInvalid; } + + /// @brief Required by std::set + bool operator<(const Point& other) const + { + if (x != other.x) + return x < other.x; + return y < other.y; + } }; /// @brief Compute the dot product of a point with itself. diff --git a/libs/MeshKernel/include/MeshKernel/Polygon.hpp b/libs/MeshKernel/include/MeshKernel/Polygon.hpp index 2c390c295..487d95ab2 100644 --- a/libs/MeshKernel/include/MeshKernel/Polygon.hpp +++ b/libs/MeshKernel/include/MeshKernel/Polygon.hpp @@ -27,6 +27,7 @@ #pragma once +#include #include #include "MeshKernel/BoundingBox.hpp" @@ -54,6 +55,10 @@ namespace meshkernel /// @brief Default move constructor. Polygon(Polygon&& copy) = default; + /// @brief Constructor + Polygon(std::span points, + Projection projection); + /// @brief Constructor Polygon(const std::vector& points, Projection projection); @@ -73,6 +78,9 @@ namespace meshkernel /// @brief Move assignment operator Polygon& operator=(Polygon&& copy); + /// @brief Determine if the polygon is empty, i.e. has no control points. + bool IsEmpty() const; + /// @brief Return the number of points in the polygon UInt Size() const; @@ -218,6 +226,11 @@ inline meshkernel::UInt meshkernel::Polygon::Size() const return static_cast(m_nodes.size()); } +inline bool meshkernel::Polygon::IsEmpty() const +{ + return m_nodes.empty(); +} + inline const std::vector& meshkernel::Polygon::Nodes() const { return m_nodes; diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index e5fa3f6c7..1621e73f7 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -33,6 +33,7 @@ #include "MeshKernel/Exceptions.hpp" #include "MeshKernel/Mesh2D.hpp" #include "MeshKernel/Mesh2DIntersections.hpp" +#include "MeshKernel/MeshBoundaryExtractor.hpp" #include "MeshKernel/MeshFaceCenters.hpp" #include "MeshKernel/MeshOrthogonality.hpp" #include "MeshKernel/Operations.hpp" @@ -217,7 +218,8 @@ void Mesh2D::DoAdministrationGivenFaceNodesMapping(const std::vector Mesh2D::SortedFacesAroundNode(UInt node) const std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector& polygonNodes) { - const Polygon polygon(polygonNodes, m_projection); - - // Find faces Administrate(); - std::vector isVisited(GetNumEdges(), false); - std::vector meshBoundaryPolygon; - meshBoundaryPolygon.reserve(GetNumNodes()); - - for (UInt e = 0; e < GetNumEdges(); e++) - { - if (isVisited[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]; - - bool firstNodeInPolygon = polygon.Contains(m_nodes[firstNodeIndex]); - bool secondNodeInPolygon = polygon.Contains(m_nodes[secondNodeIndex]); - - if (!firstNodeInPolygon && !secondNodeInPolygon) - { - continue; - } - - // Start a new polyline - if (!meshBoundaryPolygon.empty()) - { - meshBoundaryPolygon.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); - } - - // Put the current edge on the mesh boundary, mark it as visited - meshBoundaryPolygon.emplace_back(firstNode); - meshBoundaryPolygon.emplace_back(secondNode); - isVisited[e] = true; - - // walk the current mesh boundary - auto currentNode = secondNodeIndex; - WalkBoundaryFromNode(polygon, isVisited, currentNode, meshBoundaryPolygon); - - const auto numNodesFirstTail = static_cast(meshBoundaryPolygon.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(polygon, isVisited, currentNode, meshBoundaryPolygon); - } - - // There is a nonempty second tail: reverse the second tail so that the tails connect and close the polygon. - if (meshBoundaryPolygon.size() > numNodesFirstTail) - { - std::reverse(meshBoundaryPolygon.begin() + numNodesFirstTail, meshBoundaryPolygon.end()); - meshBoundaryPolygon.push_back(meshBoundaryPolygon.front()); - } - } - return meshBoundaryPolygon; + auto [boundaryPoints, isEnclosingBoundary] = GetAllBoundaryPolygons(polygonNodes); + return boundaryPoints; } -std::vector Mesh2D::ComputeInnerBoundaryPolygons() const +std::tuple, std::vector> Mesh2D::GetAllBoundaryPolygons(const std::vector& polygonNodes) { - if (GetNumFaces() == 0) - { - return std::vector(); - } - - std::vector illegalCells; - illegalCells.reserve(GetNumNodes()); - std::vector meshBoundaryPolygon; - meshBoundaryPolygon.reserve(GetNumNodes()); - std::vector subSequence; - subSequence.reserve(GetNumNodes()); + const Polygon polygon(polygonNodes, m_projection); - std::vector edgeIsVisited(GetNumEdges(), false); - std::vector nodeIsVisited(GetNumNodes(), false); + MeshBoundaryExtractor meshBoundaryExtractor; - std::vector nodeIds; - nodeIds.reserve(GetNumNodes()); + // The use of std::tie instead of a structured binding is due to limitations in the macos compiler + // The compiler error "error: capturing a structured binding is not yet supported in OpenMP" + // + // auto [boundaryPoints, isEnclosingBoundary] = meshBoundaryExtractor.Extract(*this, polygon); + // Replace the 3 lines below with the line above + std::vector> boundaryPoints; + std::vector isEnclosingBoundary; + std::tie(boundaryPoints, isEnclosingBoundary) = meshBoundaryExtractor.Extract(*this, polygon); - for (UInt e = 0; e < GetNumEdges(); e++) + if (polygon.IsEmpty()) { - 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 (!subSequence.empty()) - { - subSequence.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(subSequence.size()); - subSequence.emplace_back(firstNode); - subSequence.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, subSequence, nodeIds, meshBoundaryPolygon, illegalCells); - - const auto numNodesFirstTail = static_cast(subSequence.size()); - - // if the boundary polygon is not closed - if (currentNode != firstNodeIndex) + auto alwaysTrue = [](size_t idx [[maybe_unused]]) { - // 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, subSequence, nodeIds, meshBoundaryPolygon, illegalCells); - } - - // There is a nonempty second tail, so reverse the first tail, so that they connect. - if (subSequence.size() > numNodesFirstTail) - { - const auto start = startPolygonEdges + static_cast(std::ceil((numNodesFirstTail - startPolygonEdges + static_cast(1)) * 0.5)); + return true; + }; - for (auto n = start; n < numNodesFirstTail; n++) - { - const auto backupPoint = subSequence[n]; - const auto replaceIndex = numNodesFirstTail - n + firstNodeIndex; - subSequence[n] = subSequence[replaceIndex]; - subSequence[replaceIndex] = backupPoint; - - const UInt backupPointIndex = nodeIds[n]; - nodeIds[n] = nodeIds[replaceIndex]; - nodeIds[replaceIndex] = backupPointIndex; - } - } + return {ConcatenatePointVectors(boundaryPoints, alwaysTrue), isEnclosingBoundary}; } - OrientatePolygonsAntiClockwise(illegalCells); + std::vector containedIsEnclosingBoundary; - return illegalCells; -} - -void Mesh2D::OrientatePolygonsAntiClockwise(std::vector& polygonNodes) const -{ - UInt polygonStart = 0; - UInt polygonLength = 0; - UInt index = 0; - - while (index < polygonNodes.size()) + // NOTE: addNonEmpty lambda also updates containedIsEnclosingBoundary + auto addNonEmpty = [&boundaryPoints, &isEnclosingBoundary, &containedIsEnclosingBoundary](size_t idx) mutable { - polygonStart = index; - polygonLength = 0; - - for (UInt i = polygonStart; i < polygonNodes.size(); ++i) + if (boundaryPoints[idx].size() > 0) { - ++index; - - if (!polygonNodes[i].IsValid()) - { - polygonLength = i - polygonStart; - break; - } - - if (index == polygonNodes.size()) - { - ++polygonLength; - } + containedIsEnclosingBoundary.push_back(isEnclosingBoundary[idx]); + return true; } - if (polygonLength > 0) - { - const Point inValidPoint = {constants::missing::doubleValue, constants::missing::doubleValue}; - Point zeroPoint{0.0, 0.0}; + return false; + }; - Point midPoint = std::accumulate(polygonNodes.begin() + polygonStart, polygonNodes.begin() + polygonStart + polygonLength - 1, zeroPoint) / static_cast(polygonLength - 1); + std::vector containedBoundaryPoints = ConcatenatePointVectors(boundaryPoints, addNonEmpty); - if (!IsPointInPolygonNodes(midPoint, polygonNodes, m_projection, inValidPoint, polygonStart, polygonStart + polygonLength)) - { - // 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 quadrilateral polygon - std::swap(polygonNodes[polygonStart + 1], polygonNodes[polygonStart + 3]); - } - else - { - std::reverse(polygonNodes.begin() + polygonStart, polygonNodes.begin() + polygonStart + polygonLength); - } - } - } - } + return {containedBoundaryPoints, containedIsEnclosingBoundary}; } std::vector Mesh2D::RemoveOuterDomainBoundaryPolygon(const std::vector& polygonNodes) const @@ -1775,130 +1621,6 @@ std::vector Mesh2D::RemoveOuterDomainBoundaryPolygon(const st return innerBoundaryNodes; } -void Mesh2D::WalkBoundaryFromNode(const Polygon& polygon, - std::vector& isVisited, - UInt& currentNode, - std::vector& meshBoundaryPolygon) const -{ - UInt e = 0; - bool currentNodeInPolygon = false; - while (e < m_nodesNumEdges[currentNode]) - { - if (!currentNodeInPolygon) - { - currentNodeInPolygon = polygon.Contains(m_nodes[currentNode]); - } - - if (!currentNodeInPolygon) - { - break; - } - - const auto currentEdge = m_nodesEdges[currentNode][e]; - if (isVisited[currentEdge] || !IsEdgeOnBoundary(currentEdge)) - { - e++; - continue; - } - - currentNode = OtherNodeOfEdge(m_edges[currentEdge], currentNode); - e = 0; - currentNodeInPolygon = false; - - meshBoundaryPolygon.emplace_back(m_nodes[currentNode]); - isVisited[currentEdge] = true; - } -} - -void Mesh2D::WalkMultiBoundaryFromNode(std::vector& edgeIsVisited, - std::vector& nodeIsVisited, - UInt& currentNode, - std::vector& subSequence, - std::vector& nodeIds, - std::vector& meshBoundaryPolygon, - std::vector& illegalCells) 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) - { - size_t start = meshBoundaryPolygon.size(); - - if (!meshBoundaryPolygon.empty()) - { - meshBoundaryPolygon.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); - ++start; - } - - meshBoundaryPolygon.insert(meshBoundaryPolygon.end(), subSequence.begin() + lastIndex, subSequence.end()); - meshBoundaryPolygon.emplace_back(subSequence[lastIndex]); - - // the points making up the last found polygon - std::span currentPolygon(meshBoundaryPolygon.data() + start, meshBoundaryPolygon.data() + meshBoundaryPolygon.size()); - // Since the edge lies on a boundary, there will be only 1 attached element. - // This element will be in the 0th position - UInt connectedFace = m_edgesFaces[currentEdge][0]; - - if (!IsPointInPolygonNodes(m_facesMassCenters[connectedFace], currentPolygon, m_projection)) - { - // If the centre of this element does not lie within the polygon, then the polygon defines a hole in the mesh. - - if (!illegalCells.empty()) - { - // If illegal cells array is not empty then add the polygon separator - illegalCells.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); - } - - illegalCells.insert(illegalCells.end(), currentPolygon.begin(), currentPolygon.end()); - } - - subSequence.resize(lastIndex); - nodeIds.resize(lastIndex); - } - } - - currentNode = nextNode; - subSequence.emplace_back(m_nodes[currentNode]); - edgeIsVisited[currentEdge] = true; - nodeIsVisited[currentNode] = true; - nodeIds.emplace_back(currentNode); - } -} - std::vector Mesh2D::GetHangingEdges() const { std::vector result; @@ -2183,6 +1905,32 @@ std::unique_ptr Mesh2D::DeleteMeshFacesInPolygon(const P return UpdateFaceInformation(faceIndices, appendDeletedFaces); } +void Mesh2D::ReconstructInvalidCellsPolygon() +{ + MeshBoundaryExtractor meshBoundaryExtractor; + + // The use of std::tie instead of a structured binding is due to limitations in the macos compiler + // The compiler error "error: capturing a structured binding is not yet supported in OpenMP" + // + // auto [boundaryPoints, isEnclosingBoundary] = meshBoundaryExtractor.Extract(*this); + // Replace the 3 lines below with the line above + std::vector> boundaryPoints; + std::vector isEnclosingBoundary; + std::tie(boundaryPoints, isEnclosingBoundary) = meshBoundaryExtractor.Extract(*this); + + if (boundaryPoints.size() <= 1) + { + // There are no interior boundary polygons + m_invalidCellPolygons.clear(); + return; + } + + auto isInteriorBoundary = [&isEnclosingBoundary](size_t i) + { return !isEnclosingBoundary[i]; }; + + m_invalidCellPolygons = ConcatenatePointVectors(boundaryPoints, isInteriorBoundary); +} + std::unique_ptr Mesh2D::UpdateFaceInformation(const std::vector& faceIndices, const bool appendDeletedFaces) { std::vector facesToDelete; @@ -2273,6 +2021,8 @@ std::unique_ptr Mesh2D::UpdateFaceInformation(const std: m_faceArea.erase(m_faceArea.begin() + faceId); } + ReconstructInvalidCellsPolygon(); + return deleteMeshAction; } diff --git a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp new file mode 100644 index 000000000..ccd28fbe8 --- /dev/null +++ b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp @@ -0,0 +1,352 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2026. +// +// 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 "MeshKernel/MeshBoundaryExtractor.hpp" + +#include +#include +#include + +#include "MeshKernel/Operations.hpp" + +std::vector meshkernel::MeshBoundaryExtractor::ExtractConcatenated(const Mesh2D& mesh, BoundarySelection boundaryType) +{ + + // The use of std::tie instead of a structured binding is due to limitations in the macos compiler + // The compiler error "error: capturing a structured binding is not yet supported in OpenMP" + // + // auto [boundarySequences, isExterior] = Extract(mesh); + // Replace the 3 lines below with the line above + std::vector> boundarySequences; + std::vector isExterior; + std::tie(boundarySequences, isExterior) = Extract(mesh); + + std::vector allPoints; + + auto isBoundarySelection = [boundaryType, &isExterior](size_t idx) + { + using enum BoundarySelection; + + return (boundaryType == All) || + (boundaryType == ExteriorOnly && isExterior[idx]) || + (boundaryType == InteriorOnly && !isExterior[idx]); + }; + + return ConcatenatePointVectors(boundarySequences, isBoundarySelection); +} + +std::tuple>, std::vector> meshkernel::MeshBoundaryExtractor::Extract(const Mesh2D& mesh) +{ + const Polygon emptyPolygon; + return Extract(mesh, emptyPolygon); +} + +std::tuple>, std::vector> meshkernel::MeshBoundaryExtractor::Extract(const Mesh2D& mesh, const Polygon& polygon) +{ + + std::vector> allBoundaryPolygons; + std::vector> allTouchedFaces; + + const std::vector& meshNodes(mesh.Nodes()); + + FindBoundaryPolygons(meshNodes, mesh.Edges(), mesh.m_edgesFaces, allBoundaryPolygons, allTouchedFaces); + return SeparateAndDetermineExternality(mesh, polygon, allBoundaryPolygons, allTouchedFaces); +} + +void meshkernel::MeshBoundaryExtractor::ClipToConstrainingPolygon(const Polygon& polygon, std::vector& nodes) +{ + if (polygon.IsEmpty() || nodes.size() <= MinimumNumberOfPoints) + { + return; + } + + const size_t uniquePointCount = nodes.size() - 1; + std::vector clippedPolygon; + std::vector nodeIsContained(uniquePointCount); + UInt numberOfNodesContained = 0; + + for (size_t i = 0; i < uniquePointCount; ++i) + { + nodeIsContained[i] = polygon.Contains(nodes[i]); + numberOfNodesContained += nodeIsContained[i] ? 1 : 0; + } + + if (numberOfNodesContained < MinimumNumberOfPoints) + { + // Cannot make a reasonable polygon with only 2 nodes (and the closing node htat will be added later, making 3 nodes) + nodes.clear(); + return; + } + + clippedPolygon.reserve(nodes.size()); + + // At this point the boundary polygon is closed, so we have to ignore the last point on the sequence. + for (size_t i = 0; i < uniquePointCount; ++i) + { + size_t nextNodeId = (i + 1) % uniquePointCount; + size_t previousNodeId = (i + uniquePointCount - 1) % uniquePointCount; + + if (nodeIsContained[i] || nodeIsContained[nextNodeId] || nodeIsContained[previousNodeId]) + { + clippedPolygon.push_back(nodes[i]); + } + } + + // Now re-close the polygon + clippedPolygon.push_back(clippedPolygon[0]); + + nodes = std::move(clippedPolygon); +} + +void meshkernel::MeshBoundaryExtractor::Append(const Polygon& polygon, + const Point& centre, + const Projection projection, + std::vector& boundaryPolygon, + std::vector& isExterior, + std::vector>& separatedBoundaryPolygons) +{ + + if (auto [area, centreOfMass] = ComputePolygonAreaAndCentre(boundaryPolygon, projection); area < 0.0) + { + std::ranges::reverse(boundaryPolygon); + } + + const bool isExteriorBoundary = IsPointInPolygonNodes(centre, boundaryPolygon, projection); + + // Clipping to boundary polygon must be done after determining if the boundary-points form an exterior or interior boundary. + // Because the element centre (taken from the first element the boundary-polygon touches) may no longer be inside the clipped + // boundary-polygon + ClipToConstrainingPolygon(polygon, boundaryPolygon); + + if (boundaryPolygon.size() > MinimumNumberOfPoints) + { + // Only add the polygon and the exterior/interior indicator if the polygon has a sufficient number of points + isExterior.push_back(isExteriorBoundary); + separatedBoundaryPolygons.push_back(std::move(boundaryPolygon)); + } +} + +std::tuple>, std::vector> +meshkernel::MeshBoundaryExtractor::SeparateAndDetermineExternality(const Mesh2D& mesh, + const Polygon& polygon, + std::vector>& allBoundaryPolygons, + const std::vector>& allTouchedFaces) +{ + std::vector> separatedBoundaryPolygons; + std::vector isExterior; + + for (size_t i = 0; i < allBoundaryPolygons.size(); ++i) + { + + if (IsMultiPolygon(allBoundaryPolygons[i])) + { + // It can be that some of the boundaries that are found are composed of multiple sub-boundaries. + // Some of the sub-boundaries may be combined in a non conforming way. + // In either case, the boundaries are separated into distinct boundary polygons. + auto [individualBoundaryPolygons, firstElement] = SplitMultiplePolygons(allBoundaryPolygons[i], allTouchedFaces[i]); + + for (size_t j = 0; j < individualBoundaryPolygons.size(); ++j) + { + const Point centre = mesh.m_facesMassCenters[firstElement[j]]; + + Append(polygon, centre, mesh.m_projection, individualBoundaryPolygons[j], isExterior, separatedBoundaryPolygons); + } + } + else + { + const Point centre = mesh.m_facesMassCenters[allTouchedFaces[i][0]]; + + Append(polygon, centre, mesh.m_projection, allBoundaryPolygons[i], isExterior, separatedBoundaryPolygons); + } + } + + return {separatedBoundaryPolygons, isExterior}; +} + +double meshkernel::MeshBoundaryExtractor::NormalizeAngle(double angle) +{ + while (angle < 0) + { + angle += 2.0 * std::numbers::pi; + } + + while (angle >= 2.0 * std::numbers::pi) + { + angle -= 2.0 * std::numbers::pi; + } + + return angle; +} + +void meshkernel::MeshBoundaryExtractor::FindAllBoundarEdges(const std::vector& nodes, + const std::vector& edges, + const std::vector>& edgesFaces, + std::unordered_map>& boundaryAdjacency) +{ + + // Collect all boundary edges and compute the angle + for (UInt count = 0; count < edges.size(); ++count) + { + + if (edgesFaces[count][1] != constants::missing::uintValue || !IsValidEdge(edges[count])) + { + // Edge is either invalid or not on boundaryy + continue; + } + + double dx = nodes[edges[count].second].x - nodes[edges[count].first].x; + double dy = nodes[edges[count].second].y - nodes[edges[count].first].y; + + // Compute angle for "crossong/pinch point" sorting + // Pinch points are when the boundary polygon intersects itself, the points will have identical values. + double angleStartToEnd = NormalizeAngle(std::atan2(dy, dx)); + double angleEndToStart = NormalizeAngle(std::atan2(-dy, -dx)); + + // Note: edgesFaces[count][0] is the internal mesh face valid for both directions of boundary traversal + boundaryAdjacency[edges[count].first].push_back({count, edges[count].second, edgesFaces[count][0], angleStartToEnd}); + boundaryAdjacency[edges[count].second].push_back({count, edges[count].first, edgesFaces[count][0], angleEndToStart}); + } + + auto boundaryAngleLessThan = [](const BoundaryEdge& a, const BoundaryEdge& b) + { return a.angle < b.angle; }; + + // Sort outgoing edges anti-clockwise + for (auto& [node_id, connected_edges] : boundaryAdjacency) + { + std::sort(connected_edges.begin(), connected_edges.end(), boundaryAngleLessThan); + } +} + +meshkernel::UInt meshkernel::MeshBoundaryExtractor::FindEdgeWithMinumumAngle(const std::vector& boundaryEdges, + const std::vector& edgeVisited, + const double incomingAngle) +{ + + // Find the smallest visited angle. + // All edge angles must be in interval [0, 2pi], initialise with number greater than 2pi + double deltaAngle = 3.0 * std::numbers::pi; + UInt edgeIndex = constants::missing::uintValue; + + // Find the boundary edge that tracks closest clockwise to keep empty space on the right + for (UInt e = 0; e < boundaryEdges.size(); ++e) + { + if (edgeVisited[boundaryEdges[e].edgeId]) + { + continue; + } + + if (double delta = NormalizeAngle(incomingAngle - boundaryEdges[e].angle); delta < deltaAngle) + { + deltaAngle = delta; + edgeIndex = e; + } + } + + return edgeIndex; +} + +void meshkernel::MeshBoundaryExtractor::FindBoundaryPolygons(const std::vector& nodes, + const std::vector& edges, + const std::vector>& edgesFaces, + std::vector>& allPolygons, + std::vector>& allTouchedFaces) +{ + allPolygons.clear(); + allTouchedFaces.clear(); + + // Mapping from mesh node-id to a sequence of boundary edges, the boundary edges will be sorted by angle + std::unordered_map> boundaryAdjacency; + + FindAllBoundarEdges(nodes, edges, edgesFaces, boundaryAdjacency); + + std::vector edgeVisited(edges.size(), false); + + // Trace boundary polygons + for (UInt count = 0; count < edges.size(); ++count) + { + + if (!IsValidEdge(edges[count])) + { + continue; + } + + if (edgesFaces[count][1] != constants::missing::uintValue || edgeVisited[count]) + { + continue; + } + + std::vector currentNodes; + std::vector currentFaces; + + UInt prevNodeIndex = edges[count].first; + UInt currentNodeIndex = edges[count].second; + + currentNodes.push_back(nodes[prevNodeIndex]); + currentFaces.push_back(edgesFaces[count][0]); // First face touched by the initial edge + edgeVisited[count] = true; // Mark edge as having been visited + + double dx = nodes[currentNodeIndex].x - nodes[prevNodeIndex].x; + double dy = nodes[currentNodeIndex].y - nodes[prevNodeIndex].y; + double incomingAngle = NormalizeAngle(std::atan2(dy, dx)); + + // Polygon until we find the start node, making a closed boundary polygon + while (currentNodeIndex != edges[count].first) + { + currentNodes.push_back(nodes[currentNodeIndex]); + + const std::vector& boundaryEdges = boundaryAdjacency[currentNodeIndex]; + + UInt edgeIndex = FindEdgeWithMinumumAngle(boundaryEdges, edgeVisited, incomingAngle); + + // No unvisited edges were found + // So eigher the boundary polygon was completed or a dead-end reached. + if (edgeIndex == constants::missing::uintValue) + { + break; + } + + const BoundaryEdge& chosenEdge = boundaryEdges[edgeIndex]; + edgeVisited[chosenEdge.edgeId] = true; + + currentFaces.push_back(chosenEdge.leftFace); + + prevNodeIndex = currentNodeIndex; + currentNodeIndex = chosenEdge.neighbourNode; + incomingAngle = chosenEdge.angle; + } + + if (currentNodes.size() >= MinimumNumberOfPoints) + { + // Close the polygon + // The polygon needs to be closed here. If it is composed of multiple sub-polygons then closing here + // makes separating them later a much easier task + currentNodes.push_back(currentNodes.front()); + allPolygons.emplace_back(std::move(currentNodes)); + allTouchedFaces.push_back(std::move(currentFaces)); + } + } +} diff --git a/libs/MeshKernel/src/Operations.cpp b/libs/MeshKernel/src/Operations.cpp index 1ede27d7c..e9c0f5044 100644 --- a/libs/MeshKernel/src/Operations.cpp +++ b/libs/MeshKernel/src/Operations.cpp @@ -25,6 +25,7 @@ // //------------------------------------------------------------------------------ +#include #include #include "MeshKernel/Cartesian3DPoint.hpp" @@ -427,6 +428,110 @@ namespace meshkernel } } + std::tuple ComputePolygonAreaAndCentre(const std::vector& polygon, const Projection& projection) + { + + if (polygon.size() < constants::geometric::numNodesInTriangle) + { + throw std::invalid_argument("FaceAreaAndCenterOfMass: The polygon has less than 3 unique nodes."); + } + + Point centreOfMass(0.0, 0.0); + double area = 0.0; + + const double minArea = 1e-8; + const Point reference = ReferencePoint(polygon, projection); + const auto numberOfPointsOpenedPolygon = static_cast(polygon.size()) - 1; + + if (numberOfPointsOpenedPolygon == constants::geometric::numNodesInTriangle) + { + Vector delta1 = GetDelta(reference, polygon[0], projection); + Vector delta2 = GetDelta(reference, polygon[1], projection); + Vector delta3 = GetDelta(reference, polygon[2], projection); + + Vector middle1 = 0.5 * (delta1 + delta2); + Vector middle2 = 0.5 * (delta2 + delta3); + Vector middle3 = 0.5 * (delta3 + delta1); + + delta1 = GetDelta(polygon[0], polygon[1], projection); + delta2 = GetDelta(polygon[1], polygon[2], projection); + delta3 = GetDelta(polygon[2], polygon[0], projection); + + double xds1 = delta1.y() * middle1.x() - delta1.x() * middle1.y(); + double xds2 = delta2.y() * middle2.x() - delta2.x() * middle2.y(); + double xds3 = delta3.y() * middle3.x() - delta3.x() * middle3.y(); + + area = 0.5 * (xds1 + xds2 + xds3); + + centreOfMass += xds1 * middle1; + centreOfMass += xds2 * middle2; + centreOfMass += xds3 * middle3; + } + else if (numberOfPointsOpenedPolygon == constants::geometric::numNodesInQuadrilateral) + { + Vector delta1 = GetDelta(reference, polygon[0], projection); + Vector delta2 = GetDelta(reference, polygon[1], projection); + Vector delta3 = GetDelta(reference, polygon[2], projection); + Vector delta4 = GetDelta(reference, polygon[3], projection); + + Vector middle1 = 0.5 * (delta1 + delta2); + Vector middle2 = 0.5 * (delta2 + delta3); + Vector middle3 = 0.5 * (delta3 + delta4); + Vector middle4 = 0.5 * (delta4 + delta1); + + delta1 = GetDelta(polygon[0], polygon[1], projection); + delta2 = GetDelta(polygon[1], polygon[2], projection); + delta3 = GetDelta(polygon[2], polygon[3], projection); + delta4 = GetDelta(polygon[3], polygon[0], projection); + + double xds1 = delta1.y() * middle1.x() - delta1.x() * middle1.y(); + double xds2 = delta2.y() * middle2.x() - delta2.x() * middle2.y(); + double xds3 = delta3.y() * middle3.x() - delta3.x() * middle3.y(); + double xds4 = delta4.y() * middle4.x() - delta4.x() * middle4.y(); + + area = 0.5 * (xds1 + xds2 + xds3 + xds4); + + centreOfMass += xds1 * middle1; + centreOfMass += xds2 * middle2; + centreOfMass += xds3 * middle3; + centreOfMass += xds4 * middle4; + } + else + { + + for (UInt n = 0; n < numberOfPointsOpenedPolygon; ++n) + { + const auto nextNode = NextCircularForwardIndex(n, numberOfPointsOpenedPolygon); + + Vector delta = GetDelta(reference, polygon[n], projection); + Vector deltaNext = GetDelta(reference, polygon[nextNode], projection); + Vector middle = 0.5 * (delta + deltaNext); + delta = GetDelta(polygon[n], polygon[nextNode], projection); + + // Rotate by 3pi/2 + Vector normal(delta.y(), -delta.x()); + double xds = dot(normal, middle); + area += 0.5 * xds; + + centreOfMass += xds * middle; + } + } + + area = std::abs(area) < minArea ? minArea : area; + centreOfMass *= 1.0 / (3.0 * area); + + // TODO SHould this also apply to spheciral accurate? + if (projection == Projection::spherical) + { + centreOfMass.y /= (constants::geometric::earth_radius * constants::conversion::degToRad); + centreOfMass.x /= (constants::geometric::earth_radius * constants::conversion::degToRad * std::cos((centreOfMass.y + reference.y) * constants::conversion::degToRad)); + } + + centreOfMass += reference; + + return {area, centreOfMass}; + } + void ComputeThreeBaseComponents(const Point& point, std::array& exxp, std::array& eyyp, std::array& ezzp) { const double phi0 = point.y * constants::conversion::degToRad; @@ -1723,4 +1828,81 @@ namespace meshkernel return (matCoefficients[0] * x[0] + matCoefficients[1] * x[1]) * y[0] + (matCoefficients[2] * x[0] + matCoefficients[3] * x[1]) * y[1]; } + bool IsMultiPolygon(std::span boundary) + { + + if (boundary.size() < 4) + { + return false; + } + + std::set uniquePoints; + + // Since polygons are always closed (first = last) then start from 1 after the first + for (size_t i = 1; i < boundary.size(); ++i) + { + if (!uniquePoints.insert(boundary[i]).second) + { + return true; + } + } + + return false; + } + + std::tuple>, std::vector> SplitMultiplePolygons(std::span boundaryPoints, std::span elementIds) + { + std::vector> completedPolygons; + std::vector firstElementIds; + + std::vector> pointFaceStack; + // Maps a point to its current index in the pointFaceStack + std::map activePoints; + + for (size_t i = 0; i < boundaryPoints.size(); ++i) + { + const Point& currentPoint = boundaryPoints[i]; + + if (!currentPoint.IsValid()) + { + continue; + } + + // Check if this node closes a loop with a previously visited point + if (auto it = activePoints.find(currentPoint); it != activePoints.end()) + { + size_t loopStartIndex = it->second; + + std::vector subPolygon; + subPolygon.reserve((pointFaceStack.size() - loopStartIndex) + 1); + + int firstEdgeId = pointFaceStack[loopStartIndex].second; + + for (size_t j = loopStartIndex; j < pointFaceStack.size(); ++j) + { + subPolygon.push_back(pointFaceStack[j].first); + activePoints.erase(pointFaceStack[j].first); + } + + // close the polygon + if (!subPolygon.empty()) + { + subPolygon.push_back(subPolygon.front()); + } + + completedPolygons.push_back(subPolygon); + firstElementIds.push_back(firstEdgeId); + + pointFaceStack.resize(loopStartIndex); + } + + int currentEdge = (i < elementIds.size()) ? elementIds[i] : -1; + + activePoints[currentPoint] = pointFaceStack.size(); + pointFaceStack.emplace_back(currentPoint, currentEdge); + } + + return {completedPolygons, firstElementIds}; + } + } // namespace meshkernel diff --git a/libs/MeshKernel/src/Polygon.cpp b/libs/MeshKernel/src/Polygon.cpp index f72fd156a..c9f045dbd 100644 --- a/libs/MeshKernel/src/Polygon.cpp +++ b/libs/MeshKernel/src/Polygon.cpp @@ -42,6 +42,12 @@ meshkernel::Polygon::Polygon(const std::vector& points, Initialise(); } +meshkernel::Polygon::Polygon(std::span points, + Projection projection) : m_nodes(points.begin(), points.end()), m_projection(projection) +{ + Initialise(); +} + meshkernel::Polygon::Polygon(std::vector&& points, Projection projection) : m_nodes(points), m_projection(projection) { @@ -748,100 +754,10 @@ std::tuple meshkernel throw std::invalid_argument("FaceAreaAndCenterOfMass: The polygon has less than 3 unique nodes."); } - Point centreOfMass(0.0, 0.0); - double area = 0.0; - - const double minArea = 1e-8; - const Point reference = ReferencePoint(polygon, projection); - const auto numberOfPointsOpenedPolygon = static_cast(polygon.size()) - 1; - - if (numberOfPointsOpenedPolygon == constants::geometric::numNodesInTriangle) - { - Vector delta1 = GetDelta(reference, polygon[0], projection); - Vector delta2 = GetDelta(reference, polygon[1], projection); - Vector delta3 = GetDelta(reference, polygon[2], projection); - - Vector middle1 = 0.5 * (delta1 + delta2); - Vector middle2 = 0.5 * (delta2 + delta3); - Vector middle3 = 0.5 * (delta3 + delta1); - - delta1 = GetDelta(polygon[0], polygon[1], projection); - delta2 = GetDelta(polygon[1], polygon[2], projection); - delta3 = GetDelta(polygon[2], polygon[0], projection); - - double xds1 = delta1.y() * middle1.x() - delta1.x() * middle1.y(); - double xds2 = delta2.y() * middle2.x() - delta2.x() * middle2.y(); - double xds3 = delta3.y() * middle3.x() - delta3.x() * middle3.y(); - - area = 0.5 * (xds1 + xds2 + xds3); - - centreOfMass += xds1 * middle1; - centreOfMass += xds2 * middle2; - centreOfMass += xds3 * middle3; - } - else if (numberOfPointsOpenedPolygon == constants::geometric::numNodesInQuadrilateral) - { - Vector delta1 = GetDelta(reference, polygon[0], projection); - Vector delta2 = GetDelta(reference, polygon[1], projection); - Vector delta3 = GetDelta(reference, polygon[2], projection); - Vector delta4 = GetDelta(reference, polygon[3], projection); - - Vector middle1 = 0.5 * (delta1 + delta2); - Vector middle2 = 0.5 * (delta2 + delta3); - Vector middle3 = 0.5 * (delta3 + delta4); - Vector middle4 = 0.5 * (delta4 + delta1); - - delta1 = GetDelta(polygon[0], polygon[1], projection); - delta2 = GetDelta(polygon[1], polygon[2], projection); - delta3 = GetDelta(polygon[2], polygon[3], projection); - delta4 = GetDelta(polygon[3], polygon[0], projection); - - double xds1 = delta1.y() * middle1.x() - delta1.x() * middle1.y(); - double xds2 = delta2.y() * middle2.x() - delta2.x() * middle2.y(); - double xds3 = delta3.y() * middle3.x() - delta3.x() * middle3.y(); - double xds4 = delta4.y() * middle4.x() - delta4.x() * middle4.y(); - - area = 0.5 * (xds1 + xds2 + xds3 + xds4); - - centreOfMass += xds1 * middle1; - centreOfMass += xds2 * middle2; - centreOfMass += xds3 * middle3; - centreOfMass += xds4 * middle4; - } - else - { - - for (UInt n = 0; n < numberOfPointsOpenedPolygon; ++n) - { - const auto nextNode = NextCircularForwardIndex(n, numberOfPointsOpenedPolygon); - - Vector delta = GetDelta(reference, polygon[n], projection); - Vector deltaNext = GetDelta(reference, polygon[nextNode], projection); - Vector middle = 0.5 * (delta + deltaNext); - delta = GetDelta(polygon[n], polygon[nextNode], projection); - - // Rotate by 3pi/2 - Vector normal(delta.y(), -delta.x()); - double xds = dot(normal, middle); - area += 0.5 * xds; - - centreOfMass += xds * middle; - } - } + auto [area, centreOfMass] = ComputePolygonAreaAndCentre(polygon, projection); TraversalDirection direction = area > 0.0 ? TraversalDirection::AntiClockwise : TraversalDirection::Clockwise; - area = std::abs(area) < minArea ? minArea : area; - centreOfMass *= 1.0 / (3.0 * area); - - // TODO SHould this also apply to spheciral accurate? - if (projection == Projection::spherical) - { - centreOfMass.y /= (constants::geometric::earth_radius * constants::conversion::degToRad); - centreOfMass.x /= (constants::geometric::earth_radius * constants::conversion::degToRad * std::cos((centreOfMass.y + reference.y) * constants::conversion::degToRad)); - } - - centreOfMass += reference; return {std::abs(area), centreOfMass, direction}; } diff --git a/libs/MeshKernel/src/Utilities/Utilities.cpp b/libs/MeshKernel/src/Utilities/Utilities.cpp index 28cdec407..32151435e 100644 --- a/libs/MeshKernel/src/Utilities/Utilities.cpp +++ b/libs/MeshKernel/src/Utilities/Utilities.cpp @@ -103,9 +103,6 @@ void meshkernel::SaveVtk(const std::vector& nodes, const std::vector unsavedElements; - - unsavedElements.fill(0); UInt numberOfElements = 0; @@ -115,16 +112,6 @@ void meshkernel::SaveVtk(const std::vector& nodes, const std::vector& nodesX, const std::vector unsavedElements; - - unsavedElements.fill(0); UInt numberOfElements = 0; @@ -270,16 +254,6 @@ void meshkernel::SaveVtk(const std::vector& nodesX, const std::vector polygonNodes; + polygonNodes.push_back({-5.0, 35.0}); polygonNodes.push_back({15.0, 35.0}); polygonNodes.push_back({15.0, -5.0}); @@ -305,31 +306,103 @@ TEST(Mesh2D, MeshBoundaryToPolygonWithSelection) polygonNodes.push_back({-5.0, 35.0}); // 2 Execution - const auto meshBoundaryPolygon = mesh.ComputeBoundaryPolygons(polygonNodes); + auto meshBoundaryPolygon = mesh.ComputeBoundaryPolygons(polygonNodes); // 3 Validation const double tolerance = 1e-5; ASSERT_EQ(9, meshBoundaryPolygon.size()); ASSERT_NEAR(0.0, meshBoundaryPolygon[0].x, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[0].y, tolerance); - ASSERT_NEAR(0.0, meshBoundaryPolygon[1].x, tolerance); - ASSERT_NEAR(10.0, meshBoundaryPolygon[1].y, tolerance); - ASSERT_NEAR(0.0, meshBoundaryPolygon[2].x, tolerance); - ASSERT_NEAR(20.0, meshBoundaryPolygon[2].y, tolerance); - ASSERT_NEAR(0.0, meshBoundaryPolygon[3].x, tolerance); + + ASSERT_NEAR(10.0, meshBoundaryPolygon[1].x, tolerance); + ASSERT_NEAR(0.0, meshBoundaryPolygon[1].y, tolerance); + + ASSERT_NEAR(20.0, meshBoundaryPolygon[2].x, tolerance); + ASSERT_NEAR(0.0, meshBoundaryPolygon[2].y, tolerance); + + ASSERT_NEAR(20.0, meshBoundaryPolygon[3].x, tolerance); ASSERT_NEAR(30.0, meshBoundaryPolygon[3].y, tolerance); + ASSERT_NEAR(10.0, meshBoundaryPolygon[4].x, tolerance); ASSERT_NEAR(30.0, meshBoundaryPolygon[4].y, tolerance); - ASSERT_NEAR(20.0, meshBoundaryPolygon[5].x, tolerance); + + ASSERT_NEAR(0.0, meshBoundaryPolygon[5].x, tolerance); ASSERT_NEAR(30.0, meshBoundaryPolygon[5].y, tolerance); - ASSERT_NEAR(20.0, meshBoundaryPolygon[6].x, tolerance); - ASSERT_NEAR(00.0, meshBoundaryPolygon[6].y, tolerance); - ASSERT_NEAR(10.0, meshBoundaryPolygon[7].x, tolerance); - ASSERT_NEAR(00.0, meshBoundaryPolygon[7].y, tolerance); + + ASSERT_NEAR(0.0, meshBoundaryPolygon[6].x, tolerance); + ASSERT_NEAR(20.0, meshBoundaryPolygon[6].y, tolerance); + + ASSERT_NEAR(0.0, meshBoundaryPolygon[7].x, tolerance); + ASSERT_NEAR(10.0, meshBoundaryPolygon[7].y, tolerance); + ASSERT_NEAR(0.0, meshBoundaryPolygon[8].x, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[8].y, tolerance); } +TEST(Mesh2D, MeshBoundaryToPolygonWithAnotherSelection) +{ + // 1 Setup + auto mesh = MakeRectangularMeshForTesting(7, 7, 5.0, meshkernel::Projection::cartesian); + + std::vector polygonNodes; + + polygonNodes.push_back({-5.0, 35.0}); + polygonNodes.push_back({15.0, 35.0}); + polygonNodes.push_back({15.0, -5.0}); + polygonNodes.push_back({-5.0, -5.0}); + // Add small kink in the polygon, it should miss out a small section of the boundary + polygonNodes.push_back({-1.0, 7.5}); + polygonNodes.push_back({10.0, 12.0}); + polygonNodes.push_back({-1.0, 22.0}); + polygonNodes.push_back({-5.0, 35.0}); + + // 2 Execution + auto meshBoundaryPolygon = mesh->ComputeBoundaryPolygons(polygonNodes); + + // Note missing point at (0.0, 15.0) + std::vector expectedXs{0.0, 5.0, 10.0, 15.0, 20.0, 20.0, 15.0, 10.0, 5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + std::vector expectedYs{0.0, 0.0, 0.0, 0.0, 0.0, 30.0, 30.0, 30.0, 30.0, 30.0, 25.0, 20.0, 10.0, 5.0, 0.0}; + + // 3 Validation + const double tolerance = 1e-5; + ASSERT_EQ(15, meshBoundaryPolygon.size()); + + for (size_t i = 0; i < meshBoundaryPolygon.size(); ++i) + { + EXPECT_NEAR(expectedXs[i], meshBoundaryPolygon[i].x, tolerance); + EXPECT_NEAR(expectedYs[i], meshBoundaryPolygon[i].y, tolerance); + } +} + +TEST(Mesh2D, MeshBoundaryToPolygonWithSelectionExcludingFirstFoundPoint) +{ + // 1 Setup + auto mesh = MakeRectangularMeshForTesting(7, 7, 5.0, meshkernel::Projection::cartesian); + + // The way that the mesh is constructed means that when the full boundary polygon is computed + // the origin (0.0) and connected edges will be part of this polygon. This point is also used + // to close the polygon, since it is the first in the seuqnece. + // However this clipping polygons defined below, this point and its neighbouring points will outside of this clipping polygon + // The test is to check that this is handled correctly and that the boudnary polygon is then closed correctly. + std::vector polygonNodes{{7.5, -1.0}, {18.0, -1.0}, {18.0, 31.0}, {-1.0, 31.0}, {-1.0, 14.0}, {7.5, -1.0}}; + + // 2 Execution + auto meshBoundaryPolygon = mesh->ComputeBoundaryPolygons(polygonNodes); + + std::vector expectedXs{5.0, 10.0, 15.0, 20.0, 20.0, 15.0, 10.0, 5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.0}; + std::vector expectedYs{0.0, 0.0, 0.0, 0.0, 30.0, 30.0, 30.0, 30.0, 30.0, 25.0, 20.0, 15.0, 10.0, 0.0}; + + // 3 Validation + const double tolerance = 1e-5; + ASSERT_EQ(14, meshBoundaryPolygon.size()); + + for (size_t i = 0; i < meshBoundaryPolygon.size(); ++i) + { + EXPECT_NEAR(expectedXs[i], meshBoundaryPolygon[i].x, tolerance); + EXPECT_NEAR(expectedYs[i], meshBoundaryPolygon[i].y, tolerance); + } +} + TEST(Mesh2D, HangingEdge) { // 1 Setup diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index e22d177eb..a92a046a4 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -36,6 +36,7 @@ #include "MeshKernel/CasulliDeRefinement.hpp" #include "MeshKernel/CasulliRefinement.hpp" #include "MeshKernel/Mesh2D.hpp" +#include "MeshKernel/MeshBoundaryExtractor.hpp" #include "MeshKernel/MeshEdgeLength.hpp" #include "MeshKernel/MeshFaceCenters.hpp" #include "MeshKernel/MeshRefinement.hpp" @@ -2890,7 +2891,7 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygonsForSixFa auto deleteMeshFacesUndoAction = mesh.DeleteMeshFacesInPolygon(boundaryWithMissingElements); // Compute interior boundary polygon points - std::vector boundaryNodes2 = mesh.ComputeInnerBoundaryPolygons(); + std::vector boundaryNodes2 = MeshBoundaryExtractor::ExtractConcatenated(mesh, BoundarySelection::InteriorOnly); // The expected number of points include the land boundary points UInt expectedNumberOfNodes = 26; @@ -2901,43 +2902,43 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygonsForSixFa // interior set of polygons std::vector expectedXPoints{ 95.0, - 105.0, 100.0, + 105.0, 95.0, constants::missing::doubleValue, 85.0, - 85.0, 95.0, 95.0, 85.0, + 85.0, constants::missing::doubleValue, 80.0, - 75.0, 85.0, + 75.0, 80.0, constants::missing::doubleValue, 120.0, - 115.0, 125.0, + 115.0, 120.0, constants::missing::doubleValue, 125.0, - 125.0, 135.0, 135.0, + 125.0, 125.0}; std::vector expectedYPoints{ - 15.0, 15.0, 0.0, 15.0, + 15.0, constants::missing::doubleValue, 15.0, + 15.0, 25.0, 25.0, 15.0, - 15.0, constants::missing::doubleValue, 0.0, 15.0, @@ -2950,9 +2951,9 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygonsForSixFa 0.0, constants::missing::doubleValue, 125.0, + 125.0, 135.0, 135.0, - 125.0, 125.0}; for (size_t i = 0; i < expectedXPoints.size(); ++i) @@ -3014,10 +3015,6 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon 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}, @@ -3053,19 +3050,12 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon 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 form 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); @@ -3078,75 +3068,59 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon // This should not fill in the holes in the mesh mesh2.Administrate(); - // Get interior boundary polygon points - std::vector innerBoundaryPoints = mesh2.GetInnerBoundaryPolygons(); + // After having computed an administrate, the deleted elememnts will have ebeen found and removed again + // by the illegal cells polygons. + + 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 elementNodes6{{constants::missing::doubleValue, constants::missing::doubleValue}, + mesh.Node(node61), + mesh.Node(node62), + mesh.Node(node63), + mesh.Node(node61)}; + + // Now delete another cell and compute another administrate. + // The original deletd cells should remain deleted, and with the addiitonal deleted cell. + Polygons boundaryWithMissingElements2(elementNodes6, Projection::cartesian); + auto deleteMeshFacesUndoAction2 = mesh2.DeleteMeshFacesInPolygon(boundaryWithMissingElements2); + mesh2.Administrate(); + + MeshBoundaryExtractor extractor; + + auto interiorBoundaryPoints = extractor.ExtractConcatenated(mesh2, meshkernel::BoundarySelection::InteriorOnly); // The expected number of points, should not include any land boundary points - UInt expectedNumberOfNodes = 26; + constexpr UInt expectedNumberOfNodes = 26; - ASSERT_EQ(expectedNumberOfNodes, innerBoundaryPoints.size()); + ASSERT_EQ(expectedNumberOfNodes, interiorBoundaryPoints.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, - 105.0, - 100.0, - 95.0, - constants::missing::doubleValue, - 85.0, - 85.0, - 95.0, - 95.0, - 85.0, - constants::missing::doubleValue, - 80.0, - 75.0, - 85.0, - 80.0, - constants::missing::doubleValue, - 120.0, - 115.0, - 125.0, - 120.0, - constants::missing::doubleValue, - 125.0, - 125.0, - 135.0, - 135.0, - 125.0}; - - std::vector expectedYPoints{ - 15.0, - 15.0, - 0.0, - 15.0, - constants::missing::doubleValue, - 15.0, - 25.0, - 25.0, - 15.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}; + 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, 135.0, 135.0, 125.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, 125.0, 135.0, 135.0, 125.0}; for (size_t i = 0; i < expectedXPoints.size(); ++i) { - EXPECT_EQ(expectedXPoints[i], innerBoundaryPoints[i].x); - EXPECT_EQ(expectedYPoints[i], innerBoundaryPoints[i].y); + EXPECT_EQ(expectedXPoints[i], interiorBoundaryPoints[i].x); + EXPECT_EQ(expectedYPoints[i], interiorBoundaryPoints[i].y); } } diff --git a/libs/MeshKernel/tests/src/MeshTests.cpp b/libs/MeshKernel/tests/src/MeshTests.cpp index d89e2e041..aef9995a0 100644 --- a/libs/MeshKernel/tests/src/MeshTests.cpp +++ b/libs/MeshKernel/tests/src/MeshTests.cpp @@ -271,9 +271,9 @@ TEST(Mesh, MeshBoundaryToPolygon) ASSERT_NEAR(0.0, meshBoundaryPolygon[4].x, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[0].y, tolerance); - ASSERT_NEAR(5.0, meshBoundaryPolygon[1].y, tolerance); + ASSERT_NEAR(-5.0, meshBoundaryPolygon[1].y, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[2].y, tolerance); - ASSERT_NEAR(-5.0, meshBoundaryPolygon[3].y, tolerance); + ASSERT_NEAR(5.0, meshBoundaryPolygon[3].y, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[4].y, tolerance); } diff --git a/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp b/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp index 33c5987c4..5e799a062 100644 --- a/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp +++ b/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp @@ -82,18 +82,28 @@ TEST(InvalidCellsPolygonsTests, MeshHolesAreMainainedAfterRefinement) //-------------------------------- - 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, + 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, + 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; @@ -191,15 +201,33 @@ TEST(InvalidCellsPolygonsTests, MeshHolesAreMainainedAfterRefinement) //----------------------- // Check the inner boundary polygon are correct - std::vector expectedInnerX{80.0, 85.0, 75.0, 80.0, -999.0, 100.0, 105.0, 95.0, 100.0, -999.0, 120.0, 125.0, 115.0, 120.0, -999.0, 180.0, 200.0, 200.0, 180.0, 180.0, -999.0, 95.0, 95.0, 85.0, 85.0, 95.0, -999.0, 125.0, 135.0, 135.0, 125.0, 125.0}; - std::vector expectedInnerY{0.0, 15.0, 15.0, 0.0, -999.0, 0.0, 15.0, 15.0, 0.0, -999.0, 0.0, 15.0, 15.0, 0.0, -999.0, 140.0, 140.0, 160.0, 160.0, 140.0, -999.0, 15.0, 25.0, 25.0, 15.0, 15.0, -999.0, 125.0, 125.0, 135.0, 135.0, 125.0}; + std::vector expectedInnerX{95.0, 100.0, 105.0, 95.0, + meshkernel::constants::missing::doubleValue, + 85.0, 95.0, 95.0, 85.0, 85.0, + meshkernel::constants::missing::doubleValue, + 80.0, 85.0, 75.0, 80.0, + meshkernel::constants::missing::doubleValue, + 120.0, 125.0, 115.0, 120.0, + meshkernel::constants::missing::doubleValue, + 125.0, 135.0, 135.0, 125.0, 125.0}; + + std::vector expectedInnerY{15.0, 0.0, 15.0, 15.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, + meshkernel::constants::missing::doubleValue, + 125.0, 125.0, 135.0, 135.0, 125.0}; int innerPolygonSize = 0; meshkernelapi::GeometryList innerPolygon; errorCode = meshkernelapi::mkernel_mesh2d_get_mesh_inner_boundaries_as_polygons_dimension(meshKernelId, innerPolygonSize); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); - ASSERT_EQ(innerPolygonSize, 32); + + ASSERT_EQ(innerPolygonSize, 26); innerPolygon.num_coordinates = innerPolygonSize; std::vector xInner(innerPolygon.num_coordinates);