GRIDEDIT-1548 Fixed a couple of bugs, e.g. not initialising the node …

…tree in that interpolator
Deltares · Dec 9, 2024 · ad44119 · ad44119
1 parent 8a53c85
commit ad44119
Show file tree

Hide file tree

Showing 7 changed files with 54 additions and 15 deletions.
diff --git a/libs/MeshKernel/include/MeshKernel/MeshTriangulation.hpp b/libs/MeshKernel/include/MeshKernel/MeshTriangulation.hpp
@@ -48,13 +48,21 @@ namespace meshkernel
     class BoundedArray
     {
     public:
+        BoundedArray() : m_size(0) {}
+
         UInt size() const
         {
             return m_size;
         }
 
         void push_back(const UInt index)
         {
+            if (m_size == Dimension - 1)
+            {
+                throw ConstraintError("array already at size limit: {}", Dimension);
+            }
+
+            // m_indices.push_back(index);
             m_indices[m_size] = index;
             ++m_size;
         }
@@ -71,6 +79,16 @@ namespace meshkernel
 
         // begin and end (probably const only needed) for stl algos
 
+        // auto begin() const
+        // {
+        //     return m_indices.begin();
+        // }
+
+        // auto end() const
+        // {
+        //     return m_indices.begin() + m_size;
+        // }
+
         std::array<UInt, Dimension>::const_iterator begin() const
         {
             return m_indices.begin();
@@ -87,6 +105,7 @@ namespace meshkernel
         }
 
     private:
+        // std::vector<UInt> m_indices;
         std::array<UInt, Dimension> m_indices;
         UInt m_size = 0;
     };
@@ -154,12 +173,13 @@ namespace meshkernel
         void Compute(const std::span<const double>& xNodes,
                      const std::span<const double>& yNodes);
 
-        std::vector<Point> m_nodes;                                          ///< x-node values
-        std::vector<UInt> m_faceNodes;                                       ///< Face nodes flat array passed to the triangulation library
-        std::vector<UInt> m_edgeNodes;                                       ///< Edge nodes flat array passed to the triangulation library
-        std::vector<UInt> m_faceEdges;                                       ///< Face edges flat array passed to the triangulation library
-        std::vector<std::array<UInt, 2>> m_edgesFaces;                       ///< edge-face connectivity, generated from triangulation data
-        std::vector<BoundedArray<MaximumNumberOfEdgesPerNode>> m_nodesEdges; ///< node-edge connectivity, generated from triangulation data
+        std::vector<Point> m_nodes;                    ///< x-node values
+        std::vector<UInt> m_faceNodes;                 ///< Face nodes flat array passed to the triangulation library
+        std::vector<UInt> m_edgeNodes;                 ///< Edge nodes flat array passed to the triangulation library
+        std::vector<UInt> m_faceEdges;                 ///< Face edges flat array passed to the triangulation library
+        std::vector<std::array<UInt, 2>> m_edgesFaces; ///< edge-face connectivity, generated from triangulation data
+        std::vector<std::vector<UInt>> m_nodesEdges;   ///< node-edge connectivity, generated from triangulation data
+        // std::vector<BoundedArray<2 * MaximumNumberOfEdgesPerNode>> m_nodesEdges; ///< node-edge connectivity, generated from triangulation data
 
         std::vector<Point> m_elementCentres; ///< Array of the centres of the elements
 

diff --git a/libs/MeshKernel/include/MeshKernel/SampleInterpolator.hpp b/libs/MeshKernel/include/MeshKernel/SampleInterpolator.hpp
@@ -197,7 +197,11 @@ namespace meshkernel
               m_interpolationParameters(interpolationParameters),
               m_strategy(averaging::AveragingStrategyFactory::GetAveragingStrategy(interpolationParameters.m_method,
                                                                                    interpolationParameters.m_minimumNumberOfSamples,
-                                                                                   projection)) {}
+                                                                                   projection)),
+              m_nodeRTree(RTreeFactory::Create(projection))
+        {
+            m_nodeRTree->BuildTree(m_samplePoints);
+        }
 
         /// @brief Constructor.
         ///
@@ -211,7 +215,11 @@ namespace meshkernel
               m_interpolationParameters(interpolationParameters),
               m_strategy(averaging::AveragingStrategyFactory::GetAveragingStrategy(interpolationParameters.m_method,
                                                                                    interpolationParameters.m_minimumNumberOfSamples,
-                                                                                   projection)) {}
+                                                                                   projection)),
+              m_nodeRTree(RTreeFactory::Create(projection))
+        {
+            m_nodeRTree->BuildTree(m_samplePoints);
+        }
 
         /// @brief Get the number of nodes of size of the sample data.
         UInt Size() const override;
@@ -277,7 +285,6 @@ namespace meshkernel
                                                        const Projection projection,
                                                        std::vector<Sample>& sampleCache) const;
 
-
         std::vector<Point> m_samplePoints;
 
         // SHould use the m_projection from the mesh in the interpolate function or
@@ -293,7 +300,6 @@ namespace meshkernel
 
         ///< RTree of mesh nodes
         std::unique_ptr<RTreeBase> m_nodeRTree;
-
     };
 
 } // namespace meshkernel

diff --git a/libs/MeshKernel/src/CasulliRefinement.cpp b/libs/MeshKernel/src/CasulliRefinement.cpp
@@ -802,6 +802,13 @@ void meshkernel::CasulliRefinement::ComputeNewNodes(Mesh2D& mesh, std::vector<Ed
 
 void meshkernel::CasulliRefinement::StoreNewNode(const Mesh2D& mesh, const UInt nodeId, const UInt edge1Index, const UInt edge2Index, const UInt newNodeId, std::vector<EdgeNodes>& newNodes)
 {
+
+    if (newNodeId == constants::missing::uintValue)
+    {
+        // TODO is this ok to do?
+        return;
+    }
+
     UInt edgeId1 = edge1Index;
     UInt edgeId2 = edge2Index;
 

diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp
@@ -1505,9 +1505,15 @@ void Mesh2D::MakeDualFace(UInt node, double enlargementFactor, std::vector<Point
 {
     const auto sortedFacesIndices = SortedFacesAroundNode(node);
     const auto numEdges = m_nodesNumEdges[node];
+
     dualFace.reserve(m_maximumNumberOfEdgesPerNode);
     dualFace.clear();
 
+    if (sortedFacesIndices.size() == 0)
+    {
+        return;
+    }
+
     for (UInt e = 0; e < numEdges; ++e)
     {
         const auto edgeIndex = m_nodesEdges[node][e];

diff --git a/libs/MeshKernel/src/MeshTriangulation.cpp b/libs/MeshKernel/src/MeshTriangulation.cpp
@@ -201,7 +201,7 @@ meshkernel::UInt meshkernel::MeshTriangulation::FindNearestFace(const Point& pnt
         {
             const auto edgeIds = GetEdgeIds(faceId);
 
-            BoundedArray<3 * MaximumNumberOfEdgesPerNode> elementsChecked;
+            BoundedArray<4 * MaximumNumberOfEdgesPerNode> elementsChecked;
             elementsChecked.push_back(faceId);
 
             for (UInt i = 0; i < edgeIds.size(); ++i)
@@ -269,6 +269,7 @@ meshkernel::UInt meshkernel::MeshTriangulation::FindNearestFace(const Point& pnt
                     }
                 }
             }
+            std::cout << "elementsChecked.size " << elementsChecked.size() << std::endl;
         }
     }
 

diff --git a/libs/MeshKernel/src/SampleInterpolator.cpp b/libs/MeshKernel/src/SampleInterpolator.cpp
@@ -177,7 +177,6 @@ double meshkernel::SampleAveragingInterpolator::GetSearchRadiusSquared(const std
     for (const auto& value : searchPolygon)
     {
         auto const squaredDistance = ComputeSquaredDistance(interpolationPoint, value, projection);
-        std::cout << "search radius: " << value.x << ", " << value.y << "  " << squaredDistance << std::endl;
         result = std::max(result, squaredDistance);
     }
 

diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp
@@ -2907,7 +2907,7 @@ TEST(MeshRefinement, CasulliRefinementBasedOnDepthReal)
 
     mk::Print(mesh.Nodes(), mesh.Edges());
 
-    // std::cout << "min max x " << minX << "  " << maxX << std::endl;
-    // std::cout << "min max y " << minY << "  " << maxY << std::endl;
-    // std::cout << "min max d " << minD << "  " << maxD << std::endl;
+    std::cout << "min max x " << minX << "  " << maxX << std::endl;
+    std::cout << "min max y " << minY << "  " << maxY << std::endl;
+    std::cout << "min max d " << minD << "  " << maxD << std::endl;
 }