Gridedit 2292 eliminate internal element boundaries from invalid cells polygons

cmake/compiler_config.cmake

@@ -15,6 +15,16 @@ if(APPLE)
     message(STATUS "Configuring build for macOS with ${CMAKE_CXX_COMPILER_ID} (${CMAKE_CXX_COMPILER_VERSION}).")
     # Common warning and visibility flags
     add_compile_options("-fvisibility=hidden;-Wall;-Wextra;-pedantic;-Werror;-Wno-unused-function")
+
+    cmake_host_system_information(RESULT os_release QUERY OS_RELEASE)
+
+    include(CheckCXXCompilerFlag)
+    check_cxx_compiler_flag("-Wno-character-conversion" SUPPORTS_WNO_CHARACTER_CONVERSION)
+
+    if(SUPPORTS_WNO_CHARACTER_CONVERSION)
+      add_compile_options("-Wno-character-conversion")
+    endif()
+
     # Conditionally suppress unused parameter warnings (used for platform-specific compiler issues)
     if(SUPPRESS_UNUSED_PARAMETER_WARNING)
       add_compile_options("-Wno-unused-parameter")

libs/MeshKernel/include/MeshKernel/Mesh2D.hpp

@@ -461,6 +461,18 @@ namespace meshkernel
                                        std::vector<Point>& subSequence,
                                        std::vector<Point>& illegalCells) 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 Convert all mesh boundaries to a vector of polygon nodes, including holes (copynetboundstopol)
+        ///
+        /// @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<meshkernel::Point>, std::vector<bool>> GetAllBoundaryPolygons(const std::vector<Point>& polygon);
+
         /// @brief Ensure that all polynomials are orientated in the ACW direction.
         void OrientatePolygonsAntiClockwise(std::vector<Point>& polygonNodes) const;
 

libs/MeshKernel/include/MeshKernel/Polygon.hpp

@@ -27,6 +27,7 @@
 
 #pragma once
 
+#include <span>
 #include <vector>
 
 #include "MeshKernel/BoundingBox.hpp"
@@ -54,6 +55,10 @@ namespace meshkernel
         /// @brief Default move constructor.
         Polygon(Polygon&& copy) = default;
 
+        /// @brief Constructor
+        Polygon(std::span<const Point> points,
+                Projection projection);
+
         /// @brief Constructor
         Polygon(const std::vector<Point>& points,
                 Projection projection);

libs/MeshKernel/src/Mesh2D.cpp

@@ -1511,13 +1511,20 @@ std::vector<meshkernel::UInt> Mesh2D::SortedFacesAroundNode(UInt node) const
 }
 
 std::vector<meshkernel::Point> Mesh2D::ComputeBoundaryPolygons(const std::vector<Point>& polygonNodes)
+{
+    Administrate();
+    auto [boundaryPoints, isEnclosingBoundary] = GetAllBoundaryPolygons(polygonNodes);
+    return boundaryPoints;
+}
+
+std::tuple<std::vector<meshkernel::Point>, std::vector<bool>> Mesh2D::GetAllBoundaryPolygons(const std::vector<Point>& polygonNodes)
 {
     const Polygon polygon(polygonNodes, m_projection);
 
-    // Find faces
-    Administrate();
     std::vector<bool> isVisited(GetNumEdges(), false);
     std::vector<Point> meshBoundaryPolygon;
+    std::vector<bool> isEnclosingBoundary;
+
     meshBoundaryPolygon.reserve(GetNumNodes());
 
     for (UInt e = 0; e < GetNumEdges(); e++)
@@ -1527,10 +1534,10 @@ std::vector<meshkernel::Point> Mesh2D::ComputeBoundaryPolygons(const std::vector
             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];
+        const UInt firstNodeIndex = m_edges[e].first;
+        const UInt secondNodeIndex = m_edges[e].second;
+        const Point firstNode = m_nodes[firstNodeIndex];
+        const Point secondNode = m_nodes[secondNodeIndex];
 
         bool firstNodeInPolygon = polygon.Contains(m_nodes[firstNodeIndex]);
         bool secondNodeInPolygon = polygon.Contains(m_nodes[secondNodeIndex]);
@@ -1546,9 +1553,12 @@ std::vector<meshkernel::Point> Mesh2D::ComputeBoundaryPolygons(const std::vector
             meshBoundaryPolygon.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue);
         }
 
+        const size_t boundaryPolygonStartIndex = meshBoundaryPolygon.size();
+
         // 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
@@ -1571,8 +1581,17 @@ std::vector<meshkernel::Point> Mesh2D::ComputeBoundaryPolygons(const std::vector
             std::reverse(meshBoundaryPolygon.begin() + numNodesFirstTail, meshBoundaryPolygon.end());
             meshBoundaryPolygon.push_back(meshBoundaryPolygon.front());
         }
+
+        const size_t boundaryPolygonEndIndex = meshBoundaryPolygon.size();
+
+        std::span<const Point> currentPolygonSpan(std::span<const Point>(meshBoundaryPolygon.data() + boundaryPolygonStartIndex,
+                                                                         meshBoundaryPolygon.data() + boundaryPolygonEndIndex));
+        Polygon currentPolygon(currentPolygonSpan, m_projection);
+
+        isEnclosingBoundary.push_back(currentPolygon.Contains(m_facesMassCenters[m_edgesFaces[e][0]]));
     }
-    return meshBoundaryPolygon;
+
+    return {meshBoundaryPolygon, isEnclosingBoundary};
 }
 
 std::vector<meshkernel::Point> Mesh2D::ComputeInnerBoundaryPolygons() const
@@ -2183,6 +2202,45 @@ std::unique_ptr<meshkernel::UndoAction> Mesh2D::DeleteMeshFacesInPolygon(const P
     return UpdateFaceInformation(faceIndices, appendDeletedFaces);
 }
 
+void Mesh2D::ReconstructInvalidCellsPolygon()
+{
+    auto [boundaryPoints, isEnclosingBoundary] = GetAllBoundaryPolygons(std::vector<meshkernel::Point>());
+
+    Polygons polygons(boundaryPoints, m_projection);
+
+    if (polygons.GetNumPolygons() <= 1)
+    {
+        // There are no interior boundary polygons
+        return;
+    }
+
+    std::vector<Point> innerBoundaryPolygons;
+    innerBoundaryPolygons.reserve(m_invalidCellPolygons.size());
+    bool firstElement = true;
+
+    for (UInt p = 0; p < polygons.GetNumPolygons(); ++p)
+    {
+        if (isEnclosingBoundary[p])
+        {
+            continue;
+        }
+
+        if (!firstElement)
+        {
+            innerBoundaryPolygons.push_back(Point(constants::missing::doubleValue, constants::missing::doubleValue));
+        }
+        else
+        {
+            firstElement = false;
+        }
+
+        const std::vector<Point>& polygonPoints(polygons.Enclosure(p).Outer().Nodes());
+        innerBoundaryPolygons.insert(innerBoundaryPolygons.end(), polygonPoints.begin(), polygonPoints.end());
+    }
+
+    m_invalidCellPolygons = std::move(innerBoundaryPolygons);
+}
+
 std::unique_ptr<meshkernel::UndoAction> Mesh2D::UpdateFaceInformation(const std::vector<UInt>& faceIndices, const bool appendDeletedFaces)
 {
     std::vector<UInt> facesToDelete;
@@ -2273,6 +2331,8 @@ std::unique_ptr<meshkernel::UndoAction> Mesh2D::UpdateFaceInformation(const std:
         m_faceArea.erase(m_faceArea.begin() + faceId);
     }
 
+    ReconstructInvalidCellsPolygon();
+
     return deleteMeshAction;
 }
 

libs/MeshKernel/src/Polygon.cpp

@@ -42,6 +42,12 @@ meshkernel::Polygon::Polygon(const std::vector<Point>& points,
     Initialise();
 }
 
+meshkernel::Polygon::Polygon(std::span<const Point> points,
+                             Projection projection) : m_nodes(points.begin(), points.end()), m_projection(projection)
+{
+    Initialise();
+}
+
 meshkernel::Polygon::Polygon(std::vector<Point>&& points,
                              Projection projection) : m_nodes(points), m_projection(projection)
 {

libs/MeshKernel/src/Utilities/Utilities.cpp

@@ -103,9 +103,6 @@ void meshkernel::SaveVtk(const std::vector<Point>& nodes, const std::vector<std:
 
     std::string meshType = "UnstructuredGrid";
     std::string versionNumber = "1.0";
-    std::array<UInt, 15> unsavedElements;
-
-    unsavedElements.fill(0);
 
     UInt numberOfElements = 0;
 
@@ -115,16 +112,6 @@ void meshkernel::SaveVtk(const std::vector<Point>& nodes, const std::vector<std:
         {
             ++numberOfElements;
         }
-
-        if (faces[i].size() < 15)
-        {
-            ++unsavedElements[faces[i].size()];
-        }
-    }
-
-    for (size_t i = 0; i < unsavedElements.size(); ++i)
-    {
-        std::cout << "elements " << i << " = " << unsavedElements[i] << std::endl;
     }
 
     vtkFile << "<VTKFile type=\""
@@ -258,9 +245,6 @@ void meshkernel::SaveVtk(const std::vector<double>& nodesX, const std::vector<do
 
     std::string meshType = "UnstructuredGrid";
     std::string versionNumber = "1.0";
-    std::array<UInt, 15> unsavedElements;
-
-    unsavedElements.fill(0);
 
     UInt numberOfElements = 0;
 
@@ -270,16 +254,6 @@ void meshkernel::SaveVtk(const std::vector<double>& nodesX, const std::vector<do
         {
             ++numberOfElements;
         }
-
-        if (facesNumNodes[i] < 15)
-        {
-            ++unsavedElements[facesNumNodes[i]];
-        }
-    }
-
-    for (size_t i = 0; i < unsavedElements.size(); ++i)
-    {
-        std::cout << "elements " << i << " = " << unsavedElements[i] << std::endl;
     }
 
     vtkFile << "<VTKFile type=\""

libs/MeshKernel/tests/src/MeshRefinementTests.cpp

@@ -3082,67 +3082,23 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon
     std::vector<Point> innerBoundaryPoints = mesh2.GetInnerBoundaryPolygons();
 
     // The expected number of points, should not include any land boundary points
-    UInt expectedNumberOfNodes = 26;
+    constexpr UInt expectedNumberOfNodes = 22;
 
     ASSERT_EQ(expectedNumberOfNodes, innerBoundaryPoints.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<double> 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<double> 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<double> expectedXPoints{80.0, 85.0, 95.0, 100.0, 105.0, 95.0, 95.0, 85.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<double> expectedYPoints{0.0, 15.0, 15.0, 0.0, 15.0, 15.0, 25.0, 25.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 < expectedXPoints.size(); ++i)
     {

libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp

@@ -191,15 +191,25 @@ TEST(InvalidCellsPolygonsTests, MeshHolesAreMainainedAfterRefinement)
     //-----------------------
     // Check the inner boundary polygon are correct
 
-    std::vector<double> 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<double> 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<double> expectedInnerX{80.0, 85.0, 95.0, 100.0, 105.0, 95.0, 95.0, 85.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, 125.0, 135.0, 135.0, 125.0};
+
+    std::vector<double> expectedInnerY{0.0, 15.0, 15.0, 0.0, 15.0, 15.0, 25.0, 25.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, 135.0, 135.0, 125.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, 22);
 
     innerPolygon.num_coordinates = innerPolygonSize;
     std::vector<double> xInner(innerPolygon.num_coordinates);