GRIDEDIT-1548 Passing to other computer

libs/MeshKernel/include/MeshKernel/MeshTriangulation.hpp

@@ -86,6 +86,11 @@ namespace meshkernel
         /// @brief Get the node id's of the element
         std::array<UInt, 3> GetNodeIds(const UInt faceId) const;
 
+        /// @brief Get the edge id's of the element
+        std::array<UInt, 3> GetEdgeIds(const UInt faceId) const;
+
+        const std::array<UInt, 2>& GetFaceIds (const UInt edgeId) const;
+
         /// @brief Find the nearest face to the point
         UInt FindNearestFace(const Point& pnt) const;
 
@@ -104,6 +109,7 @@ namespace meshkernel
         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;       ///< Face edges flat array passed to the triangulation library
         std::vector<Point> m_elementCentres; ///< Array of the centres of the elements
 
         UInt m_numEdges{0}; ///< number of triangulated edges
@@ -220,3 +226,18 @@ inline meshkernel::Edge meshkernel::MeshTriangulation::GetEdge(const UInt edgeId
 
     return {m_edgeNodes[2 * edgeId], m_edgeNodes[2 * edgeId + 1]};
 }
+
+inline const std::array<meshkernel::UInt, 2>& meshkernel::MeshTriangulation::GetFaceIds (const UInt edgeId) const
+{
+    if (edgeId == constants::missing::uintValue)
+    {
+        throw ConstraintError("Invalid edge id");
+    }
+
+    if (edgeId >= m_numEdges)
+    {
+        throw ConstraintError("edge id out of range: {} >= {}", edgeId, m_numEdges);
+    }
+
+    return m_edgesFaces [edgeId];
+}

libs/MeshKernel/include/MeshKernel/Operations.hpp

@@ -749,9 +749,24 @@ template <class PointVector>
                     --windingNumber; // have  a valid down intersect
                 }
             }
+
+            std::cout << " crossProductValue " << crossProductValue << "  "
+                      << triangleNodes[n].y << "  "
+                      << triangleNodes[endIndex].y << "  " << point.y << "  "
+                      << "  " << windingNumber << std::endl;
+
         }
 
         isInTriangle = windingNumber == 0 ? false : true;
+        std::cout << " isInTriangle "<< std::boolalpha << isInTriangle << std::endl;
+
+        if (7.11e5 <= point.x && point.x <= 7.125e5 && -5.598e6 <= point.y && point.y <= -5.597e6) {
+            std::cout << "triangle points: " << point.x << ", " << point.y << std::endl;
+            std::cout << triangleNodes[0].x << ", " << triangleNodes[0].y << std::endl;
+            std::cout << triangleNodes[1].x << ", " << triangleNodes[1].y << std::endl;
+            std::cout << triangleNodes[2].x << ", " << triangleNodes[2].y << std::endl;
+        }
+
     }
 
     if (projection == Projection::sphericalAccurate)

libs/MeshKernel/include/MeshKernel/SampleInterpolator.hpp

@@ -77,6 +77,12 @@ namespace meshkernel
         template <meshkernel::ValidConstPointArray PointVectorType, meshkernel::ValidConstDoubleArray ScalarVectorType>
         void Interpolate(const int sampleId, const PointVectorType& iterpolationNodes, ScalarVectorType& result) const;
 
+        /// @brief Interpolate the sample data set at a single interpolation point.
+        ///
+        /// If interpolation at multiple points is required then better performance
+        /// can be obtained using the Interpolate function above.
+        double Interpolate (const int sampleId, const Point& evaluationPoint) const;
+
         /// @brief Determine if the SampleInterpolator already has this sample set.
         bool Contains(const int sampleId) const;
 

libs/MeshKernel/include/MeshKernel/Utilities/RTree.hpp

@@ -137,6 +137,11 @@ namespace meshkernel
         /// @param[in] node The node
         void SearchNearestPoint(Point const& node) override;
 
+        /// @brief Gets the nearest of all nodes
+        /// @param[in] node The node
+        /// @param[in] count The number of nearest points to find (must be strictly greater than 0)
+        void SearchNearestPointCount(Point const& node, UInt count) override;
+
         /// @brief Deletes a node
         /// @param[in] position The index of the point to remove in m_points
         void DeleteNode(UInt position) override;
@@ -263,6 +268,27 @@ namespace meshkernel
 
     template <typename projection>
     void RTree<projection>::SearchNearestPoint(Point const& node)
+    {
+        SearchNearestPointCount (node, 1);
+        // if (Empty())
+        // {
+        //     throw AlgorithmError("RTree is empty, search cannot be performed");
+        // }
+
+        // m_queryCache.reserve(m_queryVectorCapacity);
+        // m_queryCache.clear();
+        // const Point2D nodeSought = Point2D(node.x, node.y);
+        // m_rtree2D.query(bgi::nearest(nodeSought, 1), std::back_inserter(m_queryCache));
+
+        // if (!m_queryCache.empty())
+        // {
+        //     m_queryIndices.clear();
+        //     m_queryIndices.emplace_back(m_queryCache[0].second);
+        // }
+    }
+
+    template <typename projection>
+    void RTree<projection>::SearchNearestPointCount(Point const& node, UInt count)
     {
         if (Empty())
         {
@@ -272,7 +298,7 @@ namespace meshkernel
         m_queryCache.reserve(m_queryVectorCapacity);
         m_queryCache.clear();
         const Point2D nodeSought = Point2D(node.x, node.y);
-        m_rtree2D.query(bgi::nearest(nodeSought, 1), std::back_inserter(m_queryCache));
+        m_rtree2D.query(bgi::nearest(nodeSought, count), std::back_inserter(m_queryCache));
 
         if (!m_queryCache.empty())
         {

libs/MeshKernel/include/MeshKernel/Utilities/RTreeBase.hpp

@@ -71,6 +71,11 @@ namespace meshkernel
         /// @param[in] node The node
         virtual void SearchNearestPoint(Point const& node) = 0;
 
+        /// @brief Gets the nearest of all nodes
+        /// @param[in] node The node
+        /// @param[in] count The number of nearest points to find (must be strictly greater than 0)
+        virtual void SearchNearestPointCount(Point const& node, UInt count) = 0;
+
         /// @brief Deletes a node
         /// @param[in] position The index of the point to remove in m_points
         virtual void DeleteNode(UInt position) = 0;

libs/MeshKernel/src/CasulliRefinement.cpp

@@ -242,14 +242,19 @@ void meshkernel::CasulliRefinement::RefineNodeMaskBasedOnDepths(const Mesh2D& me
         for (size_t j = 0; j < mesh.m_nodesEdges[i].size(); ++j)
         {
             UInt edgeId = mesh.m_nodesEdges[i][j];
+            double depth = interpolatedDepth[edgeId];
+
+            if (depth == constants::missing::doubleValue) {
+                continue;
+            }
 
             // If the current edge is already shorter than the minimum edge size then do not refine further
             if (mesh.m_edgeLengths[edgeId] <= refinementParameters.min_edge_size)
             {
                 continue;
             }
 
-            const double celerity = constants::physical::sqrt_gravity * std::sqrt(std::abs(interpolatedDepth[edgeId]));
+            const double celerity = constants::physical::sqrt_gravity * std::sqrt(std::abs(depth));
             const double waveCourant = celerity * maxDtCourant / mesh.m_edgeLengths[edgeId];
 
             std::cout << waveCourant << "  ";

libs/MeshKernel/src/MeshTriangulation.cpp

@@ -148,13 +148,107 @@ void meshkernel::MeshTriangulation::Compute(const std::span<const double>& xNode
         m_elementCentreRTree = RTreeFactory::Create(m_projection);
         m_elementCentreRTree->BuildTree(m_elementCentres);
     }
+
+
+    m_edgesFaces.resize (m_numEdges, {constants::missing::uintValue, constants::missing::uintValue});
+    size_t count = 0;
+
+    for (UInt i = 0; i < m_numFaces; ++i)
+    {
+
+        for (UInt j = 0; j < constants::geometric::numNodesInTriangle; ++j)
+        {
+            UInt edge = m_faceEdges [count];
+            ++count;
+
+            if (m_edgesFaces [edge][0] == constants::missing::uintValue)
+            {
+                m_edgesFaces [edge][0] = i;
+            }
+            else
+            {
+                m_edgesFaces [edge][1] = i;
+            }
+
+        }
+
+    }
+
 }
 
 meshkernel::UInt meshkernel::MeshTriangulation::FindNearestFace(const Point& pnt) const
 {
     m_elementCentreRTree->SearchNearestPoint(pnt);
+    // return m_elementCentreRTree->HasQueryResults() ? m_elementCentreRTree->GetQueryResult(0) : constants::missing::uintValue;
+
+    if (m_elementCentreRTree->HasQueryResults())
+    {
+        UInt faceId = m_elementCentreRTree->GetQueryResult(0);
+
+        if (PointIsInElement (pnt, faceId))
+        {
+            return faceId;
+        }
+        else
+        {
+            const auto edgeIds = GetEdgeIds (faceId);
+
+            // TODO collect the 3 neighbour elements ids, so that later
+            // if we have to check the elements that are connected to a
+            // node then we do not have to check those elements we have
+            // checked already
+
+            for (UInt i = 0; i < edgeIds.size (); ++i)
+            {
+                const auto [neighbour1, neighbour2] = GetFaceIds (edgeIds [i]);
+
+                if (neighbour1 != faceId && neighbour1 != constants::missing::uintValue)
+                {
+                    if (PointIsInElement (pnt, neighbour1))
+                    {
+                        return neighbour1;
+                    }
+
+                }
+                else if (neighbour2 != faceId && neighbour2 != constants::missing::uintValue)
+                {
+                    if (PointIsInElement (pnt, neighbour2))
+                    {
+                        return neighbour2;
+                    }
+
+                }
+
+            }
+
+        }
+    }
+
+    // else
+    {
+        return constants::missing::uintValue;
+    }
+
+    // if (m_elementCentreRTree->HasQueryResults())
+    // {
+
+    //     for (UInt i = 0; i < m_elementCentreRTree->GetQueryResultSize(); ++i )
+    //     {
+    //         UInt id = m_elementCentreRTree->GetQueryResult(i);
+
+    //         if (id < NumberOfFaces() && PointIsInElement (pnt, id))
+    //         {
+    //             return id;
+    //         }
+
+    //     }
+    // }
+
+    // // else
+    // {
+    //     return constants::missing::uintValue;
+    // }
 
-    return m_elementCentreRTree->HasQueryResults() ? m_elementCentreRTree->GetQueryResult(0) : constants::missing::uintValue;
 }
 
 std::array<meshkernel::Point, 3> meshkernel::MeshTriangulation::GetNodes(const UInt faceId) const
@@ -191,6 +285,23 @@ std::array<meshkernel::UInt, 3> meshkernel::MeshTriangulation::GetNodeIds(const
     return {m_faceNodes[faceIndexStart], m_faceNodes[faceIndexStart + 1], m_faceNodes[faceIndexStart + 2]};
 }
 
+std::array<meshkernel::UInt, 3> meshkernel::MeshTriangulation::GetEdgeIds(const UInt faceId) const
+{
+    if (faceId == constants::missing::uintValue)
+    {
+        throw ConstraintError("Invalid face id");
+    }
+
+    if (faceId >= m_numFaces)
+    {
+        throw ConstraintError("face id out of range: {} >= {}", faceId, m_numFaces);
+    }
+
+    UInt faceIndexStart = 3 * faceId;
+
+    return {m_faceEdges[faceIndexStart], m_faceEdges[faceIndexStart + 1], m_faceEdges[faceIndexStart + 2]};
+}
+
 bool meshkernel::MeshTriangulation::PointIsInElement(const Point& pnt, const UInt faceId) const
 {
     if (faceId == constants::missing::uintValue)

libs/MeshKernel/src/SampleInterpolator.cpp

@@ -38,6 +38,7 @@ double meshkernel::SampleInterpolator::InterpolateOnElement(const UInt elementId
         sampleValues[id2] == constants::missing::doubleValue ||
         sampleValues[id3] == constants::missing::doubleValue) [[unlikely]]
     {
+        std::cout << " invalid sample " << std::endl;
         return result;
     }
 
@@ -56,6 +57,7 @@ double meshkernel::SampleInterpolator::InterpolateOnElement(const UInt elementId
 
     if (std::abs(det) < 1e-12) [[unlikely]]
     {
+        std::cout << " det = " << det << std::endl;
         return result;
     }
 
@@ -66,3 +68,39 @@ double meshkernel::SampleInterpolator::InterpolateOnElement(const UInt elementId
 
     return result;
 }
+
+double meshkernel::SampleInterpolator::Interpolate (const int sampleId, const Point& evaluationPoint) const
+{
+
+    if (!Contains(sampleId))
+    {
+        throw ConstraintError("Sample interpolator does not contain the id: {}.", sampleId);
+    }
+
+    double result = constants::missing::doubleValue;
+
+    if (!evaluationPoint.IsValid())
+    {
+        std::cout << " invalid point " << std::endl;
+        return result;
+    }
+
+    const UInt elementId = m_triangulation.FindNearestFace(evaluationPoint);
+
+    if (elementId == constants::missing::uintValue)
+    {
+        std::cout << " element invalid " << std::endl;
+        return result;
+    }
+
+    if (m_triangulation.PointIsInElement(evaluationPoint, elementId))
+    {
+        const std::vector<double>& propertyValues = m_sampleData.at(sampleId);
+        result = InterpolateOnElement(elementId, evaluationPoint, propertyValues);
+    } else
+    {
+        std::cout << " point not in element " << std::endl;
+    }
+
+    return result;
+}