GRIDEDIT-1548 Added interpolation at mesh nodes, edge centres and fac…

…e centres for

triangulation interpolation

libs/MeshKernel/include/MeshKernel/SampleAveragingInterpolator.hpp

@@ -148,7 +148,7 @@ namespace meshkernel
         /// @brief Interpolate the sample data on the element at the interpolation point.
         double InterpolateOnElement(const UInt elementId, const Point& interpolationPoint, const std::vector<double>& sampleValues) const;
 
-        /// @brief Interpolate at the points from points found withing a polygon
+        /// @brief Interpolate at the points from points found within a polygon
         double ComputeOnPolygon(const int propertyId,
                                 std::vector<Point>& polygon,
                                 const Point& interpolationPoint,

libs/MeshKernel/include/MeshKernel/SampleTriangulationInterpolator.hpp

@@ -87,10 +87,8 @@ namespace meshkernel
         /// @brief Set sample data from std::span object
         void SetDataSpan(const int propertyId, const std::span<const double>& sampleData) override;
 
-        // DOes nothing at the moment.
-        void InterpolateSpan(const int propertyId [[maybe_unused]], const Mesh2D& mesh [[maybe_unused]], const Location location [[maybe_unused]], std::span<double>& result [[maybe_unused]]) const override
-        {
-        }
+        /// @brief Interpolate the sample data at the points for the location (nodes, edges, faces)
+        void InterpolateSpan(const int propertyId, const Mesh2D& mesh, const Location location, std::span<double>& result) const override;
 
         /// @brief Interpolate the sample data set at the interpolation nodes.
         void InterpolateSpan(const int propertyId, const std::span<const Point>& iterpolationNodes, std::span<double>& result) const override;

libs/MeshKernel/src/Mesh.cpp

@@ -921,6 +921,7 @@ void Mesh::Administrate(CompoundUndoAction* undoAction)
 
 void Mesh::AdministrateNodesEdges(CompoundUndoAction* undoAction)
 {
+    m_edgesCenters.clear();
     SetUnConnectedNodesAndEdgesToInvalid(undoAction);
 
     // return if there are no nodes or no edges

libs/MeshKernel/src/SampleTriangulationInterpolator.cpp

@@ -1,4 +1,3 @@
-#include "MeshKernel/Exceptions.hpp"
 #include "MeshKernel/SampleTriangulationInterpolator.hpp"
 
 void meshkernel::SampleTriangulationInterpolator::SetDataSpan(const int propertyId, const std::span<const double>& sampleData)
@@ -12,6 +11,33 @@ void meshkernel::SampleTriangulationInterpolator::SetDataSpan(const int property
     m_sampleData[propertyId].assign(sampleData.begin(), sampleData.end());
 }
 
+void meshkernel::SampleTriangulationInterpolator::InterpolateSpan(const int propertyId, const Mesh2D& mesh, const Location location, std::span<double>& result) const
+{
+    if (!Contains(propertyId))
+    {
+        throw ConstraintError("Sample interpolator does not contain the id: {}.", propertyId);
+    }
+
+    std::span<const Point> meshNodes;
+
+    switch (location)
+    {
+    case Location::Nodes:
+        meshNodes = std::span<const Point>(mesh.Nodes().data(), mesh.Nodes().size());
+        break;
+    case Location::Edges:
+        meshNodes = std::span<const Point>(mesh.m_edgesCenters.data(), mesh.m_edgesCenters.size());
+        break;
+    case Location::Faces:
+        meshNodes = std::span<const Point>(mesh.m_facesMassCenters.data(), mesh.m_facesMassCenters.size());
+        break;
+    default:
+        throw ConstraintError("Unknown location");
+    }
+
+    InterpolateSpan(propertyId, meshNodes, result);
+}
+
 void meshkernel::SampleTriangulationInterpolator::InterpolateSpan(const int propertyId, const std::span<const Point>& interpolationNodes, std::span<double>& result) const
 {
     if (!Contains(propertyId))

libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp

@@ -987,7 +987,7 @@ namespace meshkernelapi
         /// @brief Refine mesh using the Casulli refinement algorithm based on the depth values
         /// @param[in] meshKernelId  The id of the mesh state
         /// @param[in] polygons The polygon within which the refinement is computed.
-        /// @param[in] propertyId The identifier of the intterpoaltor ot be used
+        /// @param[in] propertyId The identifier of the interpolator to be used
         /// @param[in] meshRefinementParameters Parameters indicating how the mesh is to be refined
         /// @returns Error code
         MKERNEL_API int mkernel_mesh2d_casulli_refinement_wrt_depths(int meshKernelId,
@@ -1620,10 +1620,11 @@ namespace meshkernelapi
 
         /// @brief Sets the property data for the mesh, the sample data points do not have to match the mesh2d nodes.
         /// @param[in] projectionType The projection type used by the sample data
+        /// @param[in] interpolationType The type of interpolation required, triangulation (0) or averaging (1) (for now)
         /// @param[in] sampleData   The sample data and associated sample data points.
         /// @param[out] propertyId The id of the property
         /// @returns Error code
-        MKERNEL_API int mkernel_mesh2d_set_property(int projectionType, const GeometryList& sampleData, int& propertyId);
+        MKERNEL_API int mkernel_mesh2d_set_property(int projectionType, int interpolationType, const GeometryList& sampleData, int& propertyId);
 
         /// @brief Snaps a mesh to a land boundary.
         /// @param[in] meshKernelId The id of the mesh state

libs/MeshKernelApi/include/MeshKernelApi/PropertyCalculator.hpp

@@ -51,7 +51,7 @@ namespace meshkernelapi
         virtual void Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const = 0;
 
         /// @brief Determine the size of the vector required to store the calculated properties
-        virtual int Size(const MeshKernelState& state) const = 0;
+        virtual int Size(const MeshKernelState& state, const meshkernel::Location location) const = 0;
     };
 
     /// @brief Calculator for orthogonality of a mesh.
@@ -62,7 +62,7 @@ namespace meshkernelapi
         void Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const override;
 
         /// @brief Determine the size of the orthogonality vector required
-        int Size(const MeshKernelState& state) const override;
+        int Size(const MeshKernelState& state, const meshkernel::Location location) const override;
     };
 
     /// @brief Calculator for the edge lengths for a mesh
@@ -73,7 +73,7 @@ namespace meshkernelapi
         void Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const override;
 
         /// @brief Determine the size of the edge-length vector required
-        int Size(const MeshKernelState& state) const override;
+        int Size(const MeshKernelState& state, const meshkernel::Location location) const override;
     };
 
     /// @brief Interpolate the depths at the mesh node points.
@@ -83,6 +83,7 @@ namespace meshkernelapi
         /// @brief Constructor
         InterpolatedSamplePropertyCalculator(const GeometryList& sampleData,
                                              const meshkernel::Projection projection,
+                                             const int interpolationType,
                                              const int propertyId);
 
         /// @brief Determine is the calculator can interpolate depth values correctly
@@ -92,7 +93,7 @@ namespace meshkernelapi
         void Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const override;
 
         /// @brief Determine the size of the edge-length vector required
-        int Size(const MeshKernelState& state) const override;
+        int Size(const MeshKernelState& state, const meshkernel::Location location) const override;
 
     private:
         /// @brief Interpolator for the samples

libs/MeshKernelApi/src/MeshKernel.cpp

@@ -549,7 +549,7 @@ namespace meshkernelapi
         return lastExitCode;
     }
 
-    MKERNEL_API int mkernel_mesh2d_set_property(int projectionType, const GeometryList& sampleData, int& propertyId)
+    MKERNEL_API int mkernel_mesh2d_set_property(int projectionType, int interpolationType, const GeometryList& sampleData, int& propertyId)
     {
         lastExitCode = meshkernel::ExitCode::Success;
         propertyId = -1;
@@ -559,14 +559,16 @@ namespace meshkernelapi
             meshkernel::range_check::CheckOneOf<int>(projectionType, meshkernel::GetValidProjections(), "Projection");
             auto const projection = static_cast<meshkernel::Projection>(projectionType);
 
+            meshkernel::range_check::CheckOneOf<int>(interpolationType, {0, 1}, "InterpolationType");
+
             int localPropertyId = GeneratePropertyId();
 
             if (propertyCalculators.contains(localPropertyId))
             {
-                throw meshkernel::ConstraintError("The pproperty id already exists: id = {}.", localPropertyId);
+                throw meshkernel::ConstraintError("The property id already exists: id = {}.", localPropertyId);
             }
 
-            propertyCalculators.try_emplace(localPropertyId, std::make_unique<InterpolatedSamplePropertyCalculator>(sampleData, projection, localPropertyId));
+            propertyCalculators.try_emplace(localPropertyId, std::make_unique<InterpolatedSamplePropertyCalculator>(sampleData, projection, interpolationType, localPropertyId));
             propertyId = localPropertyId;
         }
         catch (...)
@@ -1610,11 +1612,11 @@ namespace meshkernelapi
                 throw meshkernel::MeshKernelError("The property calculator does not exist.");
             }
 
-            if (geometryList.num_coordinates < propertyCalculators[propertyValue]->Size(meshKernelState.at(meshKernelId)))
+            if (geometryList.num_coordinates < propertyCalculators[propertyValue]->Size(meshKernelState.at(meshKernelId), meshkernel::Location::Edges))
             {
                 throw meshkernel::ConstraintError("Array size too small to store property values {} < {}.",
                                                   geometryList.num_coordinates,
-                                                  propertyCalculators[propertyValue]->Size(meshKernelState.at(meshKernelId)));
+                                                  propertyCalculators[propertyValue]->Size(meshKernelState.at(meshKernelId), meshkernel::Location::Edges));
             }
 
             if (geometryList.values == nullptr)
@@ -1624,6 +1626,11 @@ namespace meshkernelapi
 
             if (propertyCalculators[propertyValue]->IsValid(meshKernelState[meshKernelId]))
             {
+                if (meshKernelState[meshKernelId].m_mesh2d->m_edgesCenters.size() == 0)
+                {
+                    meshKernelState[meshKernelId].m_mesh2d->ComputeEdgesCenters();
+                }
+
                 propertyCalculators[propertyValue]->Calculate(meshKernelState[meshKernelId], meshkernel::Location::Edges, geometryList);
             }
             else
@@ -1657,7 +1664,7 @@ namespace meshkernelapi
 
             if (propertyCalculators.contains(propertyValue) && propertyCalculators[propertyValue] != nullptr)
             {
-                dimension = propertyCalculators[propertyValue]->Size(meshKernelState[meshKernelId]);
+                dimension = propertyCalculators[propertyValue]->Size(meshKernelState[meshKernelId], meshkernel::Location::Edges);
             }
             else
             {
@@ -3201,7 +3208,7 @@ namespace meshkernelapi
             std::vector<double> depthValues(meshKernelState[meshKernelId].m_mesh2d->GetNumEdges());
             GeometryList depthGeomValues;
             depthGeomValues.num_coordinates = static_cast<int>(meshKernelState[meshKernelId].m_mesh2d->GetNumEdges());
-            depthGeomValues.values = depthArray.data();
+            depthGeomValues.values = depthValues.data();
 
             // lastExitCode will be set here if there is an exception when computing the property.
             int error = mkernel_mesh2d_get_property(meshKernelId, propertyId, depthGeomValues);

libs/MeshKernelApi/src/PropertyCalculator.cpp

@@ -11,26 +11,26 @@ bool meshkernelapi::PropertyCalculator::IsValid(const MeshKernelState& state) co
     return state.m_mesh2d != nullptr && state.m_mesh2d->GetNumNodes() > 0;
 }
 
-void meshkernelapi::OrthogonalityPropertyCalculator::Calculate(const MeshKernelState& state, const meshkernel::Location location [[maybe_unused]], const GeometryList& geometryList) const
+void meshkernelapi::OrthogonalityPropertyCalculator::Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const
 {
 
     std::vector<double> values = state.m_mesh2d->GetOrthogonality();
 
     if (static_cast<size_t>(geometryList.num_coordinates) < values.size())
     {
         throw meshkernel::ConstraintError("GeometryList with wrong dimensions, {} must be greater than or equal to {}",
-                                          geometryList.num_coordinates, Size(state));
+                                          geometryList.num_coordinates, Size(state, location));
     }
 
     std::ranges::copy(values, geometryList.values);
 }
 
-int meshkernelapi::OrthogonalityPropertyCalculator::Size(const MeshKernelState& state) const
+int meshkernelapi::OrthogonalityPropertyCalculator::Size(const MeshKernelState& state, const meshkernel::Location location [[maybe_unused]]) const
 {
     return static_cast<int>(state.m_mesh2d->GetNumEdges());
 }
 
-void meshkernelapi::EdgeLengthPropertyCalculator::Calculate(const MeshKernelState& state, const meshkernel::Location location [[maybe_unused]], const GeometryList& geometryList) const
+void meshkernelapi::EdgeLengthPropertyCalculator::Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const
 {
 
     state.m_mesh2d->ComputeEdgesLengths();
@@ -39,25 +39,37 @@ void meshkernelapi::EdgeLengthPropertyCalculator::Calculate(const MeshKernelStat
     if (static_cast<size_t>(geometryList.num_coordinates) < values.size())
     {
         throw meshkernel::ConstraintError("GeometryList with wrong dimensions, {} must be greater than or equal to {}",
-                                          geometryList.num_coordinates, Size(state));
+                                          geometryList.num_coordinates, Size(state, location));
     }
 
     std::ranges::copy(values, geometryList.values);
 }
 
-int meshkernelapi::EdgeLengthPropertyCalculator::Size(const MeshKernelState& state) const
+int meshkernelapi::EdgeLengthPropertyCalculator::Size(const MeshKernelState& state, const meshkernel::Location location [[maybe_unused]]) const
 {
     return static_cast<int>(state.m_mesh2d->GetNumEdges());
 }
 
-meshkernelapi::InterpolatedSamplePropertyCalculator::InterpolatedSamplePropertyCalculator(const GeometryList& sampleData, const meshkernel::Projection projection, const int propertyId)
+meshkernelapi::InterpolatedSamplePropertyCalculator::InterpolatedSamplePropertyCalculator(const GeometryList& sampleData,
+                                                                                          const meshkernel::Projection projection,
+                                                                                          const int interpolationType,
+                                                                                          const int propertyId)
     : m_projection(projection),
       m_propertyId(propertyId)
 {
     std::span<const double> xNodes(sampleData.coordinates_x, sampleData.num_coordinates);
     std::span<const double> yNodes(sampleData.coordinates_y, sampleData.num_coordinates);
 
-    m_sampleInterpolator = std::make_unique<meshkernel::SampleTriangulationInterpolator>(xNodes, yNodes, m_projection);
+    if (interpolationType == 0)
+    {
+        m_sampleInterpolator = std::make_unique<meshkernel::SampleTriangulationInterpolator>(xNodes, yNodes, m_projection);
+    }
+    else if (interpolationType == 1)
+    {
+        // Need to pass from api.
+        meshkernel::InterpolationParameters interpolationParameters{};
+        m_sampleInterpolator = std::make_unique<meshkernel::SampleAveragingInterpolator>(xNodes, yNodes, m_projection, interpolationParameters);
+    }
 
     std::span<const double> dataSamples(sampleData.values, sampleData.num_coordinates);
     m_sampleInterpolator->SetData(m_propertyId, dataSamples);
@@ -77,7 +89,19 @@ void meshkernelapi::InterpolatedSamplePropertyCalculator::Calculate(const MeshKe
     m_sampleInterpolator->Interpolate(m_propertyId, *state.m_mesh2d, location, interpolatedSampleData);
 }
 
-int meshkernelapi::InterpolatedSamplePropertyCalculator::Size(const MeshKernelState& state) const
+int meshkernelapi::InterpolatedSamplePropertyCalculator::Size(const MeshKernelState& state, const meshkernel::Location location) const
 {
-    return static_cast<int>(state.m_mesh2d->GetNumNodes());
+    using enum meshkernel::Location;
+
+    switch (location)
+    {
+    case Nodes:
+        return static_cast<int>(state.m_mesh2d->GetNumNodes());
+    case Edges:
+        return static_cast<int>(state.m_mesh2d->GetNumEdges());
+    case Faces:
+        return static_cast<int>(state.m_mesh2d->GetNumFaces());
+    default:
+        return -1;
+    }
 }

libs/MeshKernelApi/tests/src/MeshPropertyTests.cpp

@@ -20,7 +20,7 @@ TEST(MeshPropertyTests, BathymetryTest)
 
     const int numXCoords = 4;
     const int numYCoords = 4;
-    const int numberOfCoordinates = numXCoords * numYCoords;
+    const int numberOfEdges = (numXCoords - 1) * numYCoords + numXCoords * (numYCoords - 1);
 
     const int numXBathyCoords = 5;
     const int numYBathyCoords = 5;
@@ -110,27 +110,31 @@ TEST(MeshPropertyTests, BathymetryTest)
     sampleData.coordinates_x = bathymetryXNodes.data();
     sampleData.coordinates_y = bathymetryYNodes.data();
 
-    errorCode = mkapi::mkernel_mesh2d_set_property(projectionType, sampleData, bathymetryPropertyId);
+    errorCode = mkapi::mkernel_mesh2d_set_property(projectionType, 0 /*use interpolation based on triangulation*/, sampleData, bathymetryPropertyId);
     ASSERT_EQ(mk::ExitCode::Success, errorCode);
 
     const double tolerance = 1.0e-13;
 
     std::vector<double> expectedInterpolatedData{0.0, 1.0, 2.0, 3.0,
                                                  0.0, 1.0, 2.0, 3.0,
                                                  0.0, 1.0, 2.0, 3.0,
-                                                 0.0, 1.0, 2.0, 3.0};
+                                                 0.5, 1.5, 2.5, 0.5,
+                                                 1.5, 2.5, 0.5, 1.5,
+                                                 2.5, 0.5, 1.5, 2.5};
 
     int sampleDataSize = -1;
     errorCode = mkapi::mkernel_mesh2d_get_property_dimension(meshKernelId, bathymetryPropertyId, sampleDataSize);
+
     ASSERT_EQ(mk::ExitCode::Success, errorCode);
-    ASSERT_EQ(sampleDataSize, numberOfCoordinates);
+    ASSERT_EQ(sampleDataSize, numberOfEdges);
 
     mkapi::GeometryList propertyData{};
-    std::vector<double> retrievedPropertyData(numberOfCoordinates, -1.0);
-    propertyData.num_coordinates = numberOfCoordinates;
+    std::vector<double> retrievedPropertyData(numberOfEdges, -1.0);
+    propertyData.num_coordinates = numberOfEdges;
     propertyData.values = retrievedPropertyData.data();
 
     errorCode = mkapi::mkernel_mesh2d_get_property(meshKernelId, bathymetryPropertyId, propertyData);
+
     ASSERT_EQ(mk::ExitCode::Success, errorCode);
 
     for (size_t i = 0; i < retrievedPropertyData.size(); ++i)
@@ -250,7 +254,7 @@ TEST(MeshPropertyTests, PropertyFailureTest)
     ASSERT_EQ(mk::ExitCode::Success, errorCode);
     EXPECT_FALSE(hasBathymetryData);
 
-    errorCode = mkapi::mkernel_mesh2d_set_property(projectionType, sampleData, bathymetryPropertyId);
+    errorCode = mkapi::mkernel_mesh2d_set_property(projectionType, 0 /*use interpolation based on triangulation*/, sampleData, bathymetryPropertyId);
     ASSERT_EQ(mk::ExitCode::Success, errorCode);
 
     errorCode = mkapi::mkernel_mesh2d_is_valid_property(meshKernelId, bathymetryPropertyId, hasBathymetryData);