From 4b02a5d7437163875a0aeac8f322be6c1df83b88 Mon Sep 17 00:00:00 2001
From: Luca Carniato <luca.carniato@deltares.nl>
Date: Mon, 23 Oct 2023 16:50:50 +0200
Subject: [PATCH] bug fix, add unit test

---
 .../MeshKernel/include/MeshKernel/Polygon.hpp |  4 +--
 libs/MeshKernel/src/Polygon.cpp               | 16 ++++--------
 libs/MeshKernel/tests/src/Mesh2DTest.cpp      | 26 +++++++++++++++++++
 3 files changed, 33 insertions(+), 13 deletions(-)

diff --git a/libs/MeshKernel/include/MeshKernel/Polygon.hpp b/libs/MeshKernel/include/MeshKernel/Polygon.hpp
index df6a9c6d6..89c546d97 100644
--- a/libs/MeshKernel/include/MeshKernel/Polygon.hpp
+++ b/libs/MeshKernel/include/MeshKernel/Polygon.hpp
@@ -135,8 +135,8 @@ namespace meshkernel
 
         /// @brief Determine if the polygon contains the point for Cartesian coordinate system
         ///
-        /// Also for spherical coordinates
-        bool ContainsCartesian(const Point& point) const;
+        /// For cartesian and spherical coordinates
+        bool ContainsCartesianOrSpherical(const Point& point) const;
 
         /// @brief Determine if the polygon contains the point for accurate spherical coordinate system
         bool ContainsSphericalAccurate(const Point& point) const;
diff --git a/libs/MeshKernel/src/Polygon.cpp b/libs/MeshKernel/src/Polygon.cpp
index 8991516e8..43f05856c 100644
--- a/libs/MeshKernel/src/Polygon.cpp
+++ b/libs/MeshKernel/src/Polygon.cpp
@@ -78,16 +78,13 @@ void meshkernel::Polygon::Reset(const std::vector<Point>& points,
     Initialise();
 }
 
-bool meshkernel::Polygon::ContainsCartesian(const Point& point) const
+bool meshkernel::Polygon::ContainsCartesianOrSpherical(const Point& point) const
 {
     int windingNumber = 0;
 
     for (size_t n = 0; n < m_nodes.size() - 1; n++)
     {
-        // TODO always Cartesian
-        // So Dx and Dy can be simplified (no branching)
-        // Then for 2 or more points, return multiple cross product values
-        const auto crossProductValue = crossProduct(m_nodes[n], m_nodes[n + 1], m_nodes[n], point, Projection::cartesian);
+        const auto crossProductValue = crossProduct(m_nodes[n], m_nodes[n + 1], m_nodes[n], point, m_projection);
 
         if (IsEqual(crossProductValue, 0.0))
         {
@@ -206,13 +203,10 @@ bool meshkernel::Polygon::Contains(const Point& pnt) const
 
     if (m_projection == Projection::cartesian || m_projection == Projection::spherical)
     {
-        return ContainsCartesian(pnt);
-    }
-    else
-    {
-        // projection = Projection::sphericalAccurate
-        return ContainsSphericalAccurate(pnt);
+        return ContainsCartesianOrSpherical(pnt);
     }
+    // projection = Projection::sphericalAccurate
+    return ContainsSphericalAccurate(pnt);
 }
 
 // TODO does this need start and end points, probably not all the polygon is to be snapped
diff --git a/libs/MeshKernel/tests/src/Mesh2DTest.cpp b/libs/MeshKernel/tests/src/Mesh2DTest.cpp
index 92d18d417..b54483cd8 100644
--- a/libs/MeshKernel/tests/src/Mesh2DTest.cpp
+++ b/libs/MeshKernel/tests/src/Mesh2DTest.cpp
@@ -986,3 +986,29 @@ TEST(Mesh2D, DeleteMesh_WhenFacesAreIntersected_ShouldNotDeleteFaces)
     // Assert
     EXPECT_EQ(mesh->GetNumFaces(), 9);
 }
+
+TEST(Mesh2D, DeleteMesh_WhenFacesAreIntersectedSpherical_ShouldNotDeleteFaces)
+{
+    // Prepare
+    const auto mesh = MakeRectangularMeshForTesting(4, 4, 3, 3, meshkernel::Projection::spherical, meshkernel::Point{0, 0});
+
+    // a polygon including all nodes of a face, but also intersecting
+    std::vector<meshkernel::Point> polygonNodes{
+        {1.87622950819672, -0.299180327868853},
+        {1.86885245901639, 0.187704918032786},
+        {3.27049180327869, 0.195081967213114},
+        {3.27049180327869, 0.320491803278688},
+        {1.87622950819672, 0.320491803278688},
+        {1.86147540983607, 1.16147540983606},
+        {3.54344262295082, 1.18360655737705},
+        {3.55081967213115, -0.358196721311476},
+        {1.87622950819672, -0.299180327868853}};
+
+    auto polygon = meshkernel::Polygons(polygonNodes, meshkernel::Projection::spherical);
+
+    // Execute
+    mesh->DeleteMesh(polygon, 0, false);
+
+    // Assert
+    EXPECT_EQ(mesh->GetNumFaces(), 9);
+}