From 8b0aef28283381533ede4275810febbf9e93130a Mon Sep 17 00:00:00 2001 From: BillSenior <89970704+BillSenior@users.noreply.github.com> Date: Thu, 9 Jan 2025 18:12:05 +0100 Subject: [PATCH] not casulli refinement depths (#393 | gridedit 1548) --- .../casulli_refinement_patch_with_hole.txt | 734 +++++++++--------- libs/MeshKernel/CMakeLists.txt | 11 +- .../include/MeshKernel/CasulliRefinement.hpp | 6 + .../include/MeshKernel/Definitions.hpp | 1 + libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 2 +- .../include/MeshKernel/MeshTriangulation.hpp | 250 ++++++ .../include/MeshKernel/Operations.hpp | 32 +- .../SampleAveragingInterpolator.hpp | 175 +++++ .../include/MeshKernel/SampleInterpolator.hpp | 93 +++ .../SampleTriangulationInterpolator.hpp | 92 +++ .../MeshKernel/TriangulationWrapper.hpp | 16 - .../MeshKernel/Utilities/Utilities.hpp | 7 + .../MeshKernel/src/AveragingInterpolation.cpp | 7 +- libs/MeshKernel/src/CasulliRefinement.cpp | 37 +- libs/MeshKernel/src/Mesh.cpp | 1 + libs/MeshKernel/src/Mesh2D.cpp | 21 +- libs/MeshKernel/src/MeshTriangulation.cpp | 393 ++++++++++ libs/MeshKernel/src/Operations.cpp | 163 +++- .../src/SampleAveragingInterpolator.cpp | 367 +++++++++ libs/MeshKernel/src/SampleInterpolator.cpp | 51 ++ .../src/SampleTriangulationInterpolator.cpp | 136 ++++ libs/MeshKernel/src/Utilities/Utilities.cpp | 142 ++++ libs/MeshKernel/tests/CMakeLists.txt | 2 + .../tests/src/MeshPropertyTests.cpp | 156 ++++ .../tests/src/MeshRefinementTests.cpp | 59 +- .../tests/src/SampleInterpolationTests.cpp | 240 ++++++ libs/MeshKernelApi/CMakeLists.txt | 2 + .../include/MeshKernelApi/MeshKernel.hpp | 28 +- .../MeshKernelApi/PropertyCalculator.hpp | 122 +++ libs/MeshKernelApi/src/MeshKernel.cpp | 174 ++++- libs/MeshKernelApi/src/PropertyCalculator.cpp | 112 +++ libs/MeshKernelApi/tests/CMakeLists.txt | 1 + libs/MeshKernelApi/tests/src/Mesh2DTests.cpp | 3 +- .../tests/src/MeshPropertyTests.cpp | 308 ++++++++ 34 files changed, 3451 insertions(+), 493 deletions(-) create mode 100644 libs/MeshKernel/include/MeshKernel/MeshTriangulation.hpp create mode 100644 libs/MeshKernel/include/MeshKernel/SampleAveragingInterpolator.hpp create mode 100644 libs/MeshKernel/include/MeshKernel/SampleInterpolator.hpp create mode 100644 libs/MeshKernel/include/MeshKernel/SampleTriangulationInterpolator.hpp create mode 100644 libs/MeshKernel/src/MeshTriangulation.cpp create mode 100644 libs/MeshKernel/src/SampleAveragingInterpolator.cpp create mode 100644 libs/MeshKernel/src/SampleInterpolator.cpp create mode 100644 libs/MeshKernel/src/SampleTriangulationInterpolator.cpp create mode 100644 libs/MeshKernel/tests/src/MeshPropertyTests.cpp create mode 100644 libs/MeshKernel/tests/src/SampleInterpolationTests.cpp create mode 100644 libs/MeshKernelApi/include/MeshKernelApi/PropertyCalculator.hpp create mode 100644 libs/MeshKernelApi/src/PropertyCalculator.cpp create mode 100644 libs/MeshKernelApi/tests/src/MeshPropertyTests.cpp diff --git a/data/test/data/CasulliRefinement/casulli_refinement_patch_with_hole.txt b/data/test/data/CasulliRefinement/casulli_refinement_patch_with_hole.txt index 6f67d9a87..18c745b28 100644 --- a/data/test/data/CasulliRefinement/casulli_refinement_patch_with_hole.txt +++ b/data/test/data/CasulliRefinement/casulli_refinement_patch_with_hole.txt @@ -397,12 +397,18 @@ 33 34 35 +36 +37 +38 +39 +40 41 42 43 44 45 46 +47 48 49 52 @@ -411,43 +417,37 @@ 55 56 57 +58 +59 +60 +61 +62 63 64 65 66 67 68 +69 +70 +71 +72 +73 74 75 76 77 78 79 +80 +81 +82 +83 +84 85 86 87 88 -89 -90 -91 -92 -93 -94 -95 -96 -97 -98 -99 -100 -101 -102 -103 -104 -105 -106 -107 -108 -109 0 1 2 @@ -480,27 +480,46 @@ 31 33 34 -41 +36 +37 +39 +40 42 44 45 +47 48 +50 52 53 55 56 +58 59 -63 +61 +62 64 -66 +65 67 +68 +69 +70 +71 +72 +73 74 75 -77 +76 78 +79 +80 +81 +82 +83 +84 85 86 -88 +87 89 90 91 @@ -510,222 +529,203 @@ 95 96 97 -99 -100 -101 +98 +25 +25 102 -103 +26 +26 104 -105 +27 +27 106 -107 +28 +28 108 +29 +29 +109 +112 110 -111 112 -113 114 115 +113 +115 +117 116 117 +120 118 -119 -25 -25 +120 123 -26 -26 +124 +121 +124 125 -27 -27 -127 -28 -28 +128 +126 +128 129 -29 -29 +132 130 -133 -131 -133 +132 135 136 -134 +133 136 -138 137 +140 138 -141 -139 -141 -144 -145 +140 142 -145 -146 -149 +50 +50 +144 +51 +51 147 -149 -150 -153 +148 +145 +148 151 +152 +149 +152 153 156 -157 154 -157 -158 -161 +156 159 -161 +160 +157 +160 163 -59 -59 +164 +161 +164 +167 +168 165 -60 -60 168 -169 -166 +171 +172 169 172 173 -170 -173 +175 174 -177 175 177 -180 -181 +176 +177 +179 178 +179 +181 +180 181 -184 -185 +183 182 -185 -188 -189 -186 -189 -192 -193 -190 -193 -194 -196 -195 -196 -198 -197 -198 -200 -199 -200 -202 -201 -202 -204 -203 -204 +183 35 35 -121 +100 +102 +112 +102 +101 +104 +115 +104 +103 +106 +117 +106 +105 +108 +120 +108 +107 +109 +124 +109 +41 +41 +111 +114 +128 +114 +43 +43 +119 123 -133 +132 123 122 125 136 125 -124 -127 -138 +49 +49 127 -126 129 -141 +140 129 -128 -130 -145 -130 -46 -46 -132 +51 +51 +131 135 -149 +148 135 -49 -49 -140 +134 +137 +152 +137 +57 +57 +139 +142 +156 +142 +141 144 -153 +160 144 143 +147 +164 +147 146 -157 -146 -57 -57 -148 -150 -161 +151 +168 +151 150 -60 -60 -152 -156 -169 -156 +153 +172 +153 +63 +63 155 +159 +175 +159 158 -173 -158 -68 -68 -160 163 177 163 162 -165 -181 -165 -164 -168 -185 -168 167 -172 -189 -172 +179 +167 +166 171 -174 -193 -174 -79 -79 -176 -180 -196 -180 -179 -184 -198 -184 +181 +171 +170 +173 183 -188 -200 -188 -187 -192 -202 -192 -191 -194 -204 -194 +173 11 12 13 @@ -751,63 +751,63 @@ 33 34 35 +36 +37 +38 +39 +40 41 -42 -43 44 45 46 +47 +48 +49 +50 +51 52 53 54 55 56 57 +58 59 60 +61 +62 63 64 65 66 67 68 -74 +69 75 76 77 78 79 +80 +81 +82 +83 +84 85 86 87 88 89 90 +91 +92 +93 +94 +95 96 97 98 99 -100 -101 -102 -103 -104 -105 -106 -107 -108 -109 -110 -111 -112 -113 -114 -115 -116 -117 -118 -119 -120 1 2 3 @@ -840,27 +840,46 @@ 32 34 35 -42 +37 +38 +40 +41 43 45 46 +48 49 +51 53 54 56 57 +59 60 -64 +62 +63 65 -67 +66 68 +69 +70 +71 +72 +73 +74 75 76 -78 +77 79 +80 +81 +82 +83 +84 +85 86 87 -89 +88 90 91 92 @@ -870,219 +889,200 @@ 96 97 98 -100 -101 +99 102 -103 +100 +100 104 -105 +101 +101 106 -107 +103 +103 108 +105 +105 109 +107 +107 +114 111 -112 +110 +111 +42 +42 113 -114 -115 +43 +43 116 -117 +123 +119 118 119 -120 -123 -121 -121 125 122 +121 122 -127 -124 -124 129 +127 126 -126 -130 -128 -128 +127 135 -132 131 -132 -48 -48 -134 -49 -49 +130 +131 137 -144 -140 +134 +133 +134 +142 139 -140 -146 +138 +139 +144 +141 +141 +147 143 -142 143 -150 -148 -147 -148 -156 -152 151 -152 -158 +146 +145 +146 +153 +150 +149 +150 +159 155 154 155 163 -160 -159 -160 -165 +158 +157 +158 +167 162 +161 162 -168 -164 -164 -172 -167 -166 -167 -174 171 +166 +165 +166 +173 170 -171 -180 -176 -175 +169 +170 +70 +70 +174 +71 +71 176 -184 -179 +72 +72 178 -179 -188 -183 +73 +73 +180 +74 +74 182 -183 -192 -187 -186 -187 -194 -191 -190 -191 -91 -91 -195 -92 -92 -197 -93 -93 -199 -94 -94 -201 -95 -95 -203 +100 +110 +110 +101 +113 +112 +113 +103 +116 +115 +116 +105 +118 +117 +118 +107 121 -131 -131 +120 +121 +36 +36 +124 +111 +126 +126 +42 +42 +128 +119 +130 +130 122 -134 133 -134 -124 -137 +132 +133 +44 +44 136 -137 -126 -139 +127 138 -139 -128 -142 -141 -142 -41 -41 +138 +50 +50 +140 +131 145 -132 -147 -147 -48 -48 +145 +134 +149 +148 149 -140 -151 -151 -143 -154 -153 -154 52 52 +152 +139 +154 +154 +141 157 -148 -159 -159 -59 -59 +156 +157 +143 161 -152 -166 -166 -155 -170 -169 -170 -63 -63 -173 160 +161 +146 +165 +164 +165 +150 +169 +168 +169 +58 +58 +172 +155 +174 +174 +158 +176 175 -175 +176 162 178 177 178 -164 +166 +180 +179 +180 +170 182 181 182 -167 -186 -185 -186 -171 -190 -189 -190 -74 -74 -193 -176 -195 -195 -179 -197 -196 -197 +64 +64 183 -199 -198 -199 -187 -201 -200 -201 -191 -203 -202 -203 -85 -85 -204 diff --git a/libs/MeshKernel/CMakeLists.txt b/libs/MeshKernel/CMakeLists.txt index 1ef6ebad0..8de376a45 100644 --- a/libs/MeshKernel/CMakeLists.txt +++ b/libs/MeshKernel/CMakeLists.txt @@ -49,6 +49,7 @@ set( ${SRC_DIR}/Mesh2DGenerateGlobal.cpp ${SRC_DIR}/Mesh2DIntersections.cpp ${SRC_DIR}/Mesh2DToCurvilinear.cpp + ${SRC_DIR}/MeshTriangulation.cpp ${SRC_DIR}/MeshRefinement.cpp ${SRC_DIR}/Network1D.cpp ${SRC_DIR}/Operations.cpp @@ -59,6 +60,9 @@ set( ${SRC_DIR}/Polygons.cpp ${SRC_DIR}/RemoveDisconnectedRegions.cpp ${SRC_DIR}/SamplesHessianCalculator.cpp + ${SRC_DIR}/SampleAveragingInterpolator.cpp + ${SRC_DIR}/SampleInterpolator.cpp + ${SRC_DIR}/SampleTriangulationInterpolator.cpp ${SRC_DIR}/Smoother.cpp ${SRC_DIR}/SplineAlgorithms.cpp ${SRC_DIR}/Splines.cpp @@ -168,6 +172,7 @@ set( ${DOMAIN_INC_DIR}/Mesh2DToCurvilinear.hpp ${DOMAIN_INC_DIR}/MeshConversion.hpp ${DOMAIN_INC_DIR}/MeshInterpolation.hpp + ${DOMAIN_INC_DIR}/MeshTriangulation.hpp ${DOMAIN_INC_DIR}/MeshRefinement.hpp ${DOMAIN_INC_DIR}/MeshTransformation.hpp ${DOMAIN_INC_DIR}/Network1D.hpp @@ -183,6 +188,9 @@ set( ${DOMAIN_INC_DIR}/RangeCheck.hpp ${DOMAIN_INC_DIR}/RemoveDisconnectedRegions.hpp ${DOMAIN_INC_DIR}/SamplesHessianCalculator.hpp + ${DOMAIN_INC_DIR}/SampleAveragingInterpolator.hpp + ${DOMAIN_INC_DIR}/SampleInterpolator.hpp + ${DOMAIN_INC_DIR}/SampleTriangulationInterpolator.hpp ${DOMAIN_INC_DIR}/Smoother.hpp ${DOMAIN_INC_DIR}/SplineAlgorithms.hpp ${DOMAIN_INC_DIR}/Splines.hpp @@ -285,7 +293,8 @@ target_sources( ${AVERAGING_STRATEGIES_SRC_LIST} ${CURVILINEAR_GRID_SRC_LIST} ${CURVILINEAR_GRID_UNDO_ACTION_SRC_LIST} - PUBLIC + ${UTILITIES_SRC_LIST} + PUBLIC FILE_SET HEADERS BASE_DIRS ${INC_DIR} diff --git a/libs/MeshKernel/include/MeshKernel/CasulliRefinement.hpp b/libs/MeshKernel/include/MeshKernel/CasulliRefinement.hpp index d85c0dd3f..f96e03529 100644 --- a/libs/MeshKernel/include/MeshKernel/CasulliRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/CasulliRefinement.hpp @@ -32,6 +32,7 @@ #include #include "MeshKernel/Mesh2D.hpp" +#include "MeshKernel/Parameters.hpp" #include "MeshKernel/Polygons.hpp" #include "MeshKernel/UndoActions/UndoAction.hpp" @@ -94,6 +95,11 @@ namespace meshkernel /// @param [in, out] nodeMask Node mask information static void InitialiseFaceNodes(const Mesh2D& mesh, std::vector& nodeMask); + /// @brief Set the node mask based on point contained in a polygon + static void RegisterNodesInsidePolygon(const Mesh2D& mesh, + const Polygons& polygon, + std::vector& nodeMask); + /// @brief Initialise the node mask array. /// /// @param [in] mesh Mesh used to initialise the node mask diff --git a/libs/MeshKernel/include/MeshKernel/Definitions.hpp b/libs/MeshKernel/include/MeshKernel/Definitions.hpp index ebca71a72..a521e4ee6 100644 --- a/libs/MeshKernel/include/MeshKernel/Definitions.hpp +++ b/libs/MeshKernel/include/MeshKernel/Definitions.hpp @@ -27,6 +27,7 @@ #pragma once +#include #include #include #include diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 70a4a6ae9..3ee122688 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -254,7 +254,7 @@ namespace meshkernel /// @param[in] node The node index /// @param[in] enlargementFactor The factor by which the dual face is enlarged /// @param[out] dualFace The dual face to be calculated - void MakeDualFace(UInt node, double enlargementFactor, std::vector& dualFace); + void MakeDualFace(UInt node, double enlargementFactor, std::vector& dualFace) const; /// @brief Sorts the faces around a node, sorted in counter clock wise order /// @param[in] node The node index diff --git a/libs/MeshKernel/include/MeshKernel/MeshTriangulation.hpp b/libs/MeshKernel/include/MeshKernel/MeshTriangulation.hpp new file mode 100644 index 000000000..0b899f1f1 --- /dev/null +++ b/libs/MeshKernel/include/MeshKernel/MeshTriangulation.hpp @@ -0,0 +1,250 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation version 3. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include +#include +#include +#include + +#include "MeshKernel/Constants.hpp" +#include "MeshKernel/Definitions.hpp" +#include "MeshKernel/Entities.hpp" +#include "MeshKernel/Exceptions.hpp" +#include "MeshKernel/Operations.hpp" +#include "MeshKernel/Point.hpp" +#include "MeshKernel/Utilities/RTreeFactory.hpp" + +namespace meshkernel +{ + + /// @brief A simple bounded stack + template + class BoundedStack + { + public: + /// @brief Number of elements in the array + UInt size() const + { + return m_size; + } + + /// @brief Add an element to the end of the array + void push_back(const UInt index) + { + if (m_size == Dimension - 1) + { + throw ConstraintError("array already at size limit: {}", Dimension); + } + + m_indices[m_size] = index; + ++m_size; + } + + /// @brief Get the element at the position + UInt operator[](const UInt index) const + { + return m_indices[index]; + } + + /// @brief Get the element at the position + UInt& operator[](const UInt index) + { + return m_indices[index]; + } + + /// @brief The iterator at the start of the array + std::array::const_iterator begin() const + { + return m_indices.begin(); + } + + /// @brief The iterator at the end of the array + std::array::const_iterator end() const + { + return m_indices.begin() + m_size; + } + + /// @brief Does the array contain the element value or not. + bool contains(const UInt index) const + { + return std::find(begin(), end(), index) != end(); + } + + private: + /// @brief stack based array containing the values. + std::array m_indices; + + /// @brief The current number of elements in the array + UInt m_size = 0; + }; + + /// @brief Contains a mesh triangulated from a set of points. + /// + /// Contains the original set of nodes, the edges connecting nodes + /// and the set of elements/faces making up the triangulation. + class MeshTriangulation + { + public: + /// @brief Constructor with array of points + MeshTriangulation(const std::span nodes, + const Projection projection); + + /// @brief Constructor with separate arrays of x- and y-coordinates + MeshTriangulation(const std::span xNodes, + const std::span yNodes, + const Projection projection); + + /// @brief Get the projection type used in the triangulation + Projection GetProjection() const; + + /// @brief Get the number of nodes in the triangulation + UInt NumberOfNodes() const; + + /// @brief Get the number of edges in the triangulation + UInt NumberOfEdges() const; + + /// @brief Get the number of faces/elements in the triangulation + UInt NumberOfFaces() const; + + /// @brief Get node as the position + Point GetNode(const UInt nodeId) const; + + /// @brief Get edge as the position + Edge GetEdge(const UInt nodeId) const; + + /// @brief Get the node values of the element + std::array GetNodes(const UInt faceId) const; + + /// @brief Get the node id's of the element + std::array GetNodeIds(const UInt faceId) const; + + /// @brief Get the edge id's of the element + std::array GetEdgeIds(const UInt faceId) const; + + /// @brief Get the id's of faces either side of the edge. + /// + /// May return invalid identifier in one or both values + const std::array& GetFaceIds(const UInt edgeId) const; + + /// @brief Find the nearest face to the point + UInt FindNearestFace(const Point& pnt) const; + + /// @brief Determine if the point lies within the element + bool PointIsInElement(const Point& pnt, const UInt faceId) const; + + private: + static constexpr UInt MaximumNumberOfEdgesPerNode = 16; ///< Maximum number of edges per node + + /// @brief Compute the triangulation. + void Compute(const std::span& xNodes, + const std::span& yNodes); + + std::vector m_nodes; ///< x-node values + std::vector m_faceNodes; ///< Face nodes flat array passed to the triangulation library + std::vector m_edgeNodes; ///< Edge nodes flat array passed to the triangulation library + std::vector m_faceEdges; ///< Face edges flat array passed to the triangulation library + std::vector> m_edgesFaces; ///< edge-face connectivity, generated from triangulation data + std::vector> m_nodesEdges; ///< node-edge connectivity, generated from triangulation data + + std::vector m_elementCentres; ///< Array of the centres of the elements + + UInt m_numEdges{0}; ///< number of triangulated edges + UInt m_numFaces{0}; ///< number of triangulated faces + + Projection m_projection = Projection::cartesian; ///< The projection used + std::unique_ptr m_elementCentreRTree; ///< RTree of element centres + std::unique_ptr m_nodeRTree; ///< RTree of mesh nods + }; + +} // namespace meshkernel + +inline meshkernel::Projection meshkernel::MeshTriangulation::GetProjection() const +{ + return m_projection; +} + +inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfNodes() const +{ + return static_cast(m_nodes.size()); +} + +inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfEdges() const +{ + return m_numEdges; +} + +inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfFaces() const +{ + return m_numFaces; +} + +inline meshkernel::Point meshkernel::MeshTriangulation::GetNode(const UInt nodeId) const +{ + if (nodeId == constants::missing::uintValue) + { + throw ConstraintError("Invalid node id"); + } + + if (nodeId >= NumberOfNodes()) + { + throw ConstraintError("node id out of range: {} >= {}", nodeId, NumberOfNodes()); + } + + return m_nodes[nodeId]; +} + +inline meshkernel::Edge meshkernel::MeshTriangulation::GetEdge(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_edgeNodes[2 * edgeId], m_edgeNodes[2 * edgeId + 1]}; +} + +inline const std::array& 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]; +} diff --git a/libs/MeshKernel/include/MeshKernel/Operations.hpp b/libs/MeshKernel/include/MeshKernel/Operations.hpp index 0a86b21a1..a34bf5296 100644 --- a/libs/MeshKernel/include/MeshKernel/Operations.hpp +++ b/libs/MeshKernel/include/MeshKernel/Operations.hpp @@ -29,6 +29,8 @@ #include #include +#include +#include #include "MeshKernel/BoundingBox.hpp" #include "MeshKernel/Constants.hpp" @@ -300,7 +302,7 @@ namespace meshkernel /// @param[in] endNode The end index in polygonNodes /// @param[in] projection The coordinate system projection. /// @param[in] polygonCenter A coordinate needed in case of sphericalAccurate projection - /// @returns If point is inside the designed polygon + /// @returns If point is inside the designated polygon [[nodiscard]] bool IsPointInPolygonNodes(const Point& point, const std::vector& polygonNodes, const Projection& projection, @@ -308,6 +310,32 @@ namespace meshkernel UInt startNode = constants::missing::uintValue, UInt endNode = constants::missing::uintValue); + /// @brief Checks if a point is in polygonNodes using the winding number method + /// @param[in] point The point to check + /// @param[in] polygonNodes A series of closed polygons + /// @param[in] projection The coordinate system projection. + /// @param[in] boundingBox The bounding box of the polygon nodes + /// @param[in] startNode The start index in polygonNodes + /// @param[in] endNode The end index in polygonNodes + /// @param[in] polygonCenter A coordinate needed in case of sphericalAccurate projection + /// @returns If point is inside the designated polygon + [[nodiscard]] bool IsPointInPolygonNodes(const Point& point, + const std::vector& polygonNodes, + const Projection& projection, + const BoundingBox& boundingBox, + Point polygonCenter = {constants::missing::doubleValue, constants::missing::doubleValue}, + UInt startNode = constants::missing::uintValue, + UInt endNode = constants::missing::uintValue); + + /// @brief Checks if a point is in triangle using the winding number method + /// @param[in] point The point to check + /// @param[in] triangleNodes A series of node forming an open triangle + /// @param[in] projection The coordinate system projection. + /// @returns If point is inside the designated triangle + [[nodiscard]] bool IsPointInTriangle(const Point& point, + const std::span triangleNodes, + const Projection& projection); + /// @brief Computes three base components void ComputeThreeBaseComponents(const Point& point, std::array& exxp, std::array& eyyp, std::array& ezzp); @@ -566,7 +594,7 @@ namespace meshkernel /// @param[in] points The point series. /// @param[in] projection The projection to use. /// @return The average coordinate. - [[nodiscard]] Point ComputeAverageCoordinate(const std::vector& points, const Projection& projection); + [[nodiscard]] Point ComputeAverageCoordinate(const std::span points, const Projection& projection); /// @brief Cartesian projection of a point on a segment defined by other two points /// @param firstNode The first node of the segment diff --git a/libs/MeshKernel/include/MeshKernel/SampleAveragingInterpolator.hpp b/libs/MeshKernel/include/MeshKernel/SampleAveragingInterpolator.hpp new file mode 100644 index 000000000..9638291d1 --- /dev/null +++ b/libs/MeshKernel/include/MeshKernel/SampleAveragingInterpolator.hpp @@ -0,0 +1,175 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation version 3. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include +#include +#include +#include + +#include "MeshKernel/AveragingInterpolation.hpp" // Only for the enum Method. should move to Definitions.hpp +#include "MeshKernel/AveragingStrategies/AveragingStrategy.hpp" +#include "MeshKernel/AveragingStrategies/AveragingStrategyFactory.hpp" +#include "MeshKernel/Constants.hpp" +#include "MeshKernel/Definitions.hpp" +#include "MeshKernel/Entities.hpp" +#include "MeshKernel/Exceptions.hpp" +#include "MeshKernel/Mesh2D.hpp" +#include "MeshKernel/MeshTriangulation.hpp" +#include "MeshKernel/Operations.hpp" +#include "MeshKernel/Point.hpp" +#include "MeshKernel/SampleInterpolator.hpp" +#include "MeshKernel/Utilities/RTreeFactory.hpp" + +namespace meshkernel +{ + + /// @brief Parameters used by the averaging interpolation + struct InterpolationParameters + { + /// @brief Which averaging method should be used + AveragingInterpolation::Method m_method = AveragingInterpolation::Method::SimpleAveraging; + + /// @brief The absolute search radius + double m_absoluteSearchRadius = 100.0; + + /// @brief The relative search radius + double m_relativeSearchRadius = 1.0; + + /// @brief If no point is found in polygon then just used the closest point + bool m_useClosestIfNoneFound = true; + + /// @brief The minimum number of samples for several averaging methods. + UInt m_minimumNumberOfSamples = 10; + }; + + /// @brief Interpolator for sample data using an averaging scheme + class SampleAveragingInterpolator : public SampleInterpolator + { + public: + /// @brief Constructor. + SampleAveragingInterpolator(const std::span xNodes, + const std::span yNodes, + const Projection projection, + const InterpolationParameters& interpolationParameters); + + /// @brief Constructor. + SampleAveragingInterpolator(const std::span nodes, + const Projection projection, + const InterpolationParameters& interpolationParameters); + + /// @brief Get the number of nodes of size of the sample data. + UInt Size() const override; + + /// @brief Interpolate the sample data set at the interpolation nodes. + void Interpolate(const int propertyId, const std::span iterpolationNodes, std::span result) const override; + + /// @brief Interpolate the sample data set at the locationd defined. + void Interpolate(const int propertyId, const Mesh2D& mesh, const Location location, std::span result) const override; + + /// @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 InterpolateValue(const int propertyId, const Point& evaluationPoint) const override; + + private: + static constexpr UInt MaximumNumberOfEdgesPerNode = 16; ///< Maximum number of edges per node + + /// @brief Combine x- and y-coordinate arrays to a point array + static std::vector CombineCoordinates(const std::span xNodes, const std::span yNodes); + + /// @brief Interpolate the sample data on the element at the interpolation point. + double InterpolateOnElement(const UInt elementId, const Point& interpolationPoint, const std::vector& sampleValues) const; + + /// @brief Interpolate at the points from points found within a polygon + double ComputeOnPolygon(const int propertyId, + std::vector& polygon, + const Point& interpolationPoint, + const Projection projection, + std::vector& sampleCache) const; + + /// @brief Compute the search radius + double GetSearchRadiusSquared(const std::vector& searchPolygon, + const Point& interpolationPoint, + const Projection projection) const; + + /// @brief Compute the search polygon + void GenerateSearchPolygon(const double relativeSearchRadius, + const Point& interpolationPoint, + std::vector& polygon, + const Projection projection) const; + + /// @brief Get the sample value at the index from the rtree result cache + double GetSampleValueFromRTree(const int propertyId, const UInt index) const; + + /// @brief Compute the average from the neighbours. + double ComputeInterpolationResultFromNeighbors(const int propertyId, + const Point& interpolationPoint, + const std::vector& searchPolygon, + const Projection projection, + std::vector& sampleCache) const; + + /// @brief Interpolate at the mesh nodes + void InterpolateAtNodes(const int propertyId, const Mesh2D& mesh, std::span& result) const; + + /// @brief Interpolate at edge centres by averaging the node values + void InterpolateAtEdgeCentres(const Mesh2D& mesh, + const std::span& nodeResult, + std::span& result) const; + + /// @brief Interpolate at the face centres + void InterpolateAtFaces(const int propertyId, const Mesh2D& mesh, std::span& result) const; + + /// @brief The sample points + std::vector m_samplePoints; + + // Should use the m_projection from the mesh in the interpolate function or + // needs to be passed to the interpolate function in the no mesh version + /// @brief The projection used + Projection m_projection = Projection::cartesian; + + /// @brief Interpolation parameter + InterpolationParameters m_interpolationParameters; + + /// @brief Averaging strategy + std::unique_ptr m_strategy; + + /// @brief Map from sample id (int) to sample data. + std::map> m_sampleData; + + /// @brief RTree of mesh nodes + std::unique_ptr m_nodeRTree; + }; + +} // namespace meshkernel + +inline meshkernel::UInt meshkernel::SampleAveragingInterpolator::Size() const +{ + return static_cast(m_samplePoints.size()); +} diff --git a/libs/MeshKernel/include/MeshKernel/SampleInterpolator.hpp b/libs/MeshKernel/include/MeshKernel/SampleInterpolator.hpp new file mode 100644 index 000000000..502f71adc --- /dev/null +++ b/libs/MeshKernel/include/MeshKernel/SampleInterpolator.hpp @@ -0,0 +1,93 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation version 3. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include +#include +#include +#include + +#include "MeshKernel/AveragingInterpolation.hpp" +#include "MeshKernel/AveragingStrategies/AveragingStrategy.hpp" +#include "MeshKernel/AveragingStrategies/AveragingStrategyFactory.hpp" +#include "MeshKernel/Constants.hpp" +#include "MeshKernel/Definitions.hpp" +#include "MeshKernel/Entities.hpp" +#include "MeshKernel/Exceptions.hpp" +#include "MeshKernel/Mesh2D.hpp" +#include "MeshKernel/MeshTriangulation.hpp" +#include "MeshKernel/Operations.hpp" +#include "MeshKernel/Point.hpp" +#include "MeshKernel/Utilities/RTreeFactory.hpp" + +namespace meshkernel +{ + + /// @brief Interface for sample interpolation + class SampleInterpolator + { + public: + /// @brief Destructor + virtual ~SampleInterpolator() = default; + + /// @brief Set sample data + void SetData(const int propertyId, const std::span sampleData); + + /// @brief Get the number of sample points + virtual UInt Size() const = 0; + + /// @brief Determine if the SampleInterpolator already has this sample set. + bool Contains(const int propertyId) const; + + /// @brief Interpolate the sample data set at the interpolation nodes. + virtual void Interpolate(const int propertyId, const std::span iterpolationNodes, std::span result) const = 0; + + /// @brief Interpolate the sample data. + virtual void Interpolate(const int propertyId, const Mesh2D& mesh, const Location location, std::span result) const = 0; + + /// @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. + virtual double InterpolateValue(const int propertyId, const Point& evaluationPoint) const = 0; + + protected: + /// @brief Get the sample property data for the id. + const std::vector& GetSampleData(const int propertyId) const; + + private: + /// @brief Map from sample id (int) to sample data. + std::map> m_sampleData; + }; + +} // namespace meshkernel + +inline bool meshkernel::SampleInterpolator::Contains(const int propertyId) const +{ + return m_sampleData.contains(propertyId); +} diff --git a/libs/MeshKernel/include/MeshKernel/SampleTriangulationInterpolator.hpp b/libs/MeshKernel/include/MeshKernel/SampleTriangulationInterpolator.hpp new file mode 100644 index 000000000..5b931e755 --- /dev/null +++ b/libs/MeshKernel/include/MeshKernel/SampleTriangulationInterpolator.hpp @@ -0,0 +1,92 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation version 3. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include +#include +#include + +#include "MeshKernel/Definitions.hpp" +#include "MeshKernel/Entities.hpp" +#include "MeshKernel/Exceptions.hpp" +#include "MeshKernel/Mesh2D.hpp" +#include "MeshKernel/MeshTriangulation.hpp" +#include "MeshKernel/Operations.hpp" +#include "MeshKernel/Point.hpp" +#include "MeshKernel/SampleInterpolator.hpp" +#include "MeshKernel/Utilities/RTreeFactory.hpp" + +namespace meshkernel +{ + + /// @brief Interpolator for sample data on a triangulated grid. + /// + /// The triangulation does not have to match any mesh. + class SampleTriangulationInterpolator : public SampleInterpolator + { + public: + /// @brief Constructor. + SampleTriangulationInterpolator(const std::span xNodes, + const std::span yNodes, + const Projection projection) + : m_triangulation(xNodes, yNodes, projection) {} + + /// @brief Constructor. + SampleTriangulationInterpolator(const std::span nodes, + const Projection projection) + : m_triangulation(nodes, projection) {} + + /// @brief Get the number of nodes of size of the sample data. + UInt Size() const override; + + /// @brief Interpolate the sample data at the points for the location (nodes, edges, faces) + void Interpolate(const int propertyId, const Mesh2D& mesh, const Location location, std::span result) const override; + + /// @brief Interpolate the sample data set at the interpolation nodes. + void Interpolate(const int propertyId, const std::span iterpolationNodes, std::span result) const override; + + /// @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 InterpolateValue(const int propertyId, const Point& evaluationPoint) const override; + + private: + /// @brief Interpolate the sample data on the element at the interpolation point. + double InterpolateOnElement(const UInt elementId, const Point& interpolationPoint, const std::vector& sampleValues) const; + + /// @brief Triangulation of the sample points + MeshTriangulation m_triangulation; + }; + +} // namespace meshkernel + +inline meshkernel::UInt meshkernel::SampleTriangulationInterpolator::Size() const +{ + return m_triangulation.NumberOfNodes(); +} diff --git a/libs/MeshKernel/include/MeshKernel/TriangulationWrapper.hpp b/libs/MeshKernel/include/MeshKernel/TriangulationWrapper.hpp index 16a7809f6..0e4729908 100644 --- a/libs/MeshKernel/include/MeshKernel/TriangulationWrapper.hpp +++ b/libs/MeshKernel/include/MeshKernel/TriangulationWrapper.hpp @@ -226,22 +226,6 @@ namespace meshkernel return m_edgesFaces[edgeIndex][faceIndex]; } - /// @brief Retrieves the x coordinate of a triangulated node - /// @param nodeIndex The index of the node to retrieve - /// @return const reference to the x coordinate - [[nodiscard]] double GetXCoord(const UInt nodeIndex) const - { - return m_xCoordFlat[nodeIndex]; - } - - /// @brief Retrieves the y coordinate of a triangulated node - /// @param nodeIndex The index of the node to retrieve - /// @return const reference to the y coordinate - [[nodiscard]] double GetYCoord(const UInt nodeIndex) const - { - return m_yCoordFlat[nodeIndex]; - } - /// @brief Retrieves the (x,y) coordinate of a triangulated node /// @param nodeIndex The index of the node to retrieve /// @return Point diff --git a/libs/MeshKernel/include/MeshKernel/Utilities/Utilities.hpp b/libs/MeshKernel/include/MeshKernel/Utilities/Utilities.hpp index 2575c534f..678fd5f2c 100644 --- a/libs/MeshKernel/include/MeshKernel/Utilities/Utilities.hpp +++ b/libs/MeshKernel/include/MeshKernel/Utilities/Utilities.hpp @@ -28,9 +28,11 @@ #pragma once #include +#include #include #include "MeshKernel/Entities.hpp" +#include "MeshKernel/Mesh2D.hpp" #include "MeshKernel/Point.hpp" namespace meshkernel @@ -48,4 +50,9 @@ namespace meshkernel const std::vector& yNodes, const std::vector& edges, std::ostream& out = std::cout); + /// @brief Save the mesh data in a vtk file format + /// + /// @note saves only triangle and quadrilateral elements. + void SaveVtk(const std::vector& nodes, const std::vector>& faces, const std::string& fileName); + } // namespace meshkernel diff --git a/libs/MeshKernel/src/AveragingInterpolation.cpp b/libs/MeshKernel/src/AveragingInterpolation.cpp index 570ed540a..b9ccf2e5d 100644 --- a/libs/MeshKernel/src/AveragingInterpolation.cpp +++ b/libs/MeshKernel/src/AveragingInterpolation.cpp @@ -196,19 +196,21 @@ double AveragingInterpolation::ComputeInterpolationResultFromNeighbors(const Poi { m_interpolationSampleCache.clear(); + BoundingBox boundingBox(searchPolygon); + for (UInt i = 0; i < m_samplesRtree->GetQueryResultSize(); ++i) { auto const sampleIndex = m_samplesRtree->GetQueryResult(i); auto const sampleValue = m_samples[sampleIndex].value; - if (sampleValue <= constants::missing::doubleValue) + if (sampleValue == constants::missing::doubleValue) { continue; } Point samplePoint{m_samples[sampleIndex].x, m_samples[sampleIndex].y}; - if (IsPointInPolygonNodes(samplePoint, searchPolygon, m_mesh.m_projection)) + if (IsPointInPolygonNodes(samplePoint, searchPolygon, m_mesh.m_projection, boundingBox)) { m_interpolationSampleCache.emplace_back(samplePoint.x, samplePoint.y, sampleValue); } @@ -243,7 +245,6 @@ double AveragingInterpolation::ComputeOnPolygon(const std::vector& polygo } if (m_samplesRtree->HasQueryResults()) { - return ComputeInterpolationResultFromNeighbors(interpolationPoint, searchPolygon); } diff --git a/libs/MeshKernel/src/CasulliRefinement.cpp b/libs/MeshKernel/src/CasulliRefinement.cpp index c2c1454e3..fd2535aca 100644 --- a/libs/MeshKernel/src/CasulliRefinement.cpp +++ b/libs/MeshKernel/src/CasulliRefinement.cpp @@ -180,13 +180,10 @@ void meshkernel::CasulliRefinement::InitialiseFaceNodes(const Mesh2D& mesh, std: } } -std::vector meshkernel::CasulliRefinement::InitialiseNodeMask(const Mesh2D& mesh, const Polygons& polygon) +void meshkernel::CasulliRefinement::RegisterNodesInsidePolygon(const Mesh2D& mesh, + const Polygons& polygon, + std::vector& nodeMask) { - std::vector nodeMask(10 * mesh.GetNumNodes(), NodeMask::Unassigned); - - // Find nodes that are inside the polygon. - // If the polygon is empty then all nodes will be taken into account. - for (UInt i = 0; i < mesh.GetNumNodes(); ++i) { auto [containsPoint, pointIndex] = polygon.IsPointInPolygons(mesh.Node(i)); @@ -196,7 +193,16 @@ std::vector meshkernel::CasulliRefineme nodeMask[i] = NodeMask::RegisteredNode; } } +} + +std::vector meshkernel::CasulliRefinement::InitialiseNodeMask(const Mesh2D& mesh, const Polygons& polygon) +{ + std::vector nodeMask(10 * mesh.GetNumNodes(), NodeMask::Unassigned); + + // Find nodes that are inside the polygon. + // If the polygon is empty then all nodes will be taken into account. + RegisterNodesInsidePolygon(mesh, polygon, nodeMask); InitialiseBoundaryNodes(mesh, nodeMask); InitialiseCornerNodes(mesh, nodeMask); InitialiseFaceNodes(mesh, nodeMask); @@ -216,6 +222,7 @@ void meshkernel::CasulliRefinement::Administrate(Mesh2D& mesh, const UInt numNod } } + mesh.DeleteInvalidNodesAndEdges(); mesh.Administrate(); } @@ -406,6 +413,11 @@ void meshkernel::CasulliRefinement::ConnectEdges(Mesh2D& mesh, const UInt curren { UInt edgeId = mesh.m_nodesEdges[currentNode][j]; + if (mesh.m_edgesNumFaces[edgeId] == 0) + { + continue; + } + if (mesh.m_edgesNumFaces[edgeId] == 1) { if (edgeCount >= newEdges.size()) @@ -539,6 +551,11 @@ void meshkernel::CasulliRefinement::ConnectNewNodes(Mesh2D& mesh, const std::vec { const UInt edgeId = mesh.m_nodesEdges[i][j]; + if (mesh.m_edgesNumFaces[edgeId] == 0) + { + continue; + } + if (mesh.GetEdge(edgeId).first == i) { // With passing false for the collectUndo, ConnectNodes should return a null pointer @@ -627,7 +644,7 @@ void meshkernel::CasulliRefinement::ComputeNewEdgeNodes(Mesh2D& mesh, const UInt const UInt node1 = mesh.GetEdge(i).first; const UInt node2 = mesh.GetEdge(i).second; - if (node1 == constants::missing::uintValue && node2 == constants::missing::uintValue) + if (node1 == constants::missing::uintValue || node2 == constants::missing::uintValue) { continue; } @@ -675,6 +692,12 @@ void meshkernel::CasulliRefinement::ComputeNewNodes(Mesh2D& mesh, std::vector& newNodes) { + + if (newNodeId == constants::missing::uintValue) + { + return; + } + UInt edgeId1 = edge1Index; UInt edgeId2 = edge2Index; diff --git a/libs/MeshKernel/src/Mesh.cpp b/libs/MeshKernel/src/Mesh.cpp index 4b255ef0d..a6ae68ed8 100644 --- a/libs/MeshKernel/src/Mesh.cpp +++ b/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 diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 1a84f1a25..99cb616e0 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -1506,16 +1506,30 @@ std::unique_ptr Mesh2D::TriangulateFaces() return triangulationAction; } -void Mesh2D::MakeDualFace(UInt node, double enlargementFactor, std::vector& dualFace) +void Mesh2D::MakeDualFace(UInt node, double enlargementFactor, std::vector& dualFace) const { const auto sortedFacesIndices = SortedFacesAroundNode(node); const auto numEdges = m_nodesNumEdges[node]; + dualFace.reserve(m_maximumNumberOfEdgesPerNode); dualFace.clear(); + if (sortedFacesIndices.empty()) + { + return; + } + + UInt sortedFacesCount = 0; + for (UInt e = 0; e < numEdges; ++e) { const auto edgeIndex = m_nodesEdges[node][e]; + + if (m_edgesNumFaces[edgeIndex] == 0) + { + continue; + } + auto edgeCenter = m_edgesCenters[edgeIndex]; if (m_projection == Projection::spherical) @@ -1539,7 +1553,9 @@ void Mesh2D::MakeDualFace(UInt node, double enlargementFactor, std::vector. +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#include "MeshKernel/MeshTriangulation.hpp" +#include "MeshKernel/Constants.hpp" +#include "MeshKernel/Utilities/Utilities.hpp" + +#include +#include + +extern "C" +{ + /// @brief Function of the Triangle library + /// + /// \see https://www.cs.cmu.edu/~quake/triangle.html + void Triangulation(int jatri, + double const* const xs, + double const* const ys, + int ns, + int* const indx, + int* const numtri, + int* const edgeidx, + int* const numedge, + int* const triedge, + double* const xs3, + double* const ys3, + int* const ns3, + double trisize); +} + +meshkernel::MeshTriangulation::MeshTriangulation(const std::span nodes, + const Projection projection) + : m_nodes(nodes.begin(), nodes.end()), + m_projection(projection) +{ + if (nodes.size() < constants::geometric::numNodesInTriangle) + { + throw ConstraintError("Not enough nodes for a single triangle: {}", nodes.size()); + } + + std::vector xNodes(nodes.size()); + std::vector yNodes(nodes.size()); + + std::ranges::transform(nodes, xNodes.begin(), [](const Point& p) + { return p.x; }); + std::ranges::transform(nodes, yNodes.begin(), [](const Point& p) + { return p.y; }); + + Compute(xNodes, yNodes); +} + +meshkernel::MeshTriangulation::MeshTriangulation(const std::span xNodes, + const std::span yNodes, + const Projection projection) + : m_nodes(xNodes.size()), + m_projection(projection) +{ + if (xNodes.size() < constants::geometric::numNodesInTriangle) + { + throw ConstraintError("Not enough nodes for a single triangle: {}", xNodes.size()); + } + + if (xNodes.size() != yNodes.size()) + { + throw ConstraintError("Size mismatch between x- and y-node values: {} /= {}", + xNodes.size(), yNodes.size()); + } + + for (size_t i = 0; i < xNodes.size(); ++i) + { + m_nodes[i] = Point{xNodes[i], yNodes[i]}; + } + + Compute(xNodes, yNodes); +} + +void meshkernel::MeshTriangulation::Compute(const std::span& xNodes, + const std::span& yNodes) +{ + + if (xNodes.size() != yNodes.size()) + { + throw ConstraintError("x- and y-points are not the same size: {} /= {}", xNodes.size(), yNodes.size()); + } + + if (xNodes.size() < 3) + { + throw ConstraintError("There are not enough points in mesh for a triangulation: {}", xNodes.size()); + } + + m_numEdges = 0; + m_numFaces = 0; + + UInt estimatedNumberOfTriangles = 0; + + int numInputNodes = static_cast(xNodes.size()); + + if (estimatedNumberOfTriangles == 0) + { + estimatedNumberOfTriangles = static_cast(xNodes.size()) * 6 + 10; + } + + double averageTriangleArea = 0.0; + + int intTriangulationOption = 3; + int numFaces = -1; + int numEdges = 0; + + std::vector faceNodes; + std::vector edgeNodes; + std::vector faceEdges; + + // If the number of estimated triangles is not sufficient, triangulation must be repeated + while (numFaces < 0) + { + numFaces = static_cast(estimatedNumberOfTriangles); + + faceNodes.resize(estimatedNumberOfTriangles * 3); + std::ranges::fill(faceNodes, 0); + + edgeNodes.resize(estimatedNumberOfTriangles * 2); + std::ranges::fill(edgeNodes, 0); + + faceEdges.resize(estimatedNumberOfTriangles * 3); + std::ranges::fill(faceEdges, 0); + + Triangulation(intTriangulationOption, + xNodes.data(), + yNodes.data(), + numInputNodes, + faceNodes.data(), // INDX + &numFaces, + edgeNodes.data(), // EDGEINDX + &numEdges, + faceEdges.data(), // TRIEDGE + nullptr, nullptr, nullptr, + averageTriangleArea); + + if (estimatedNumberOfTriangles > 0) + { + estimatedNumberOfTriangles = static_cast(-numFaces); + } + } + + m_numFaces = static_cast(numFaces); + m_numEdges = static_cast(numEdges); + + m_faceNodes.resize(3 * m_numFaces); + m_edgeNodes.resize(2 * m_numEdges); + m_faceEdges.resize(3 * m_numFaces); + + std::transform(faceNodes.begin(), faceNodes.begin() + 3 * m_numFaces, m_faceNodes.begin(), [](const int val) + { return static_cast(val - 1); }); + + std::transform(faceEdges.begin(), faceEdges.begin() + 3 * m_numFaces, m_faceEdges.begin(), [](const int val) + { return static_cast(val - 1); }); + + std::transform(edgeNodes.begin(), edgeNodes.begin() + 2 * m_numEdges, m_edgeNodes.begin(), [](const int val) + { return static_cast(val - 1); }); + + m_elementCentres.resize(m_numFaces); + + for (UInt i = 0; i < m_numFaces; ++i) + { + m_elementCentres[i] = ComputeAverageCoordinate(GetNodes(i), m_projection); + } + + if (m_numFaces > 0) + { + m_elementCentreRTree = RTreeFactory::Create(m_projection); + m_elementCentreRTree->BuildTree(m_elementCentres); + + m_nodeRTree = RTreeFactory::Create(m_projection); + m_nodeRTree->BuildTree(m_nodes); + } + + 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; + } + } + } + + m_nodesEdges.resize(m_nodes.size()); + + for (UInt e = 0; e < m_numEdges; ++e) + { + auto [edgeFirst, edgeSecond] = GetEdge(e); + + m_nodesEdges[edgeFirst].push_back(e); + m_nodesEdges[edgeSecond].push_back(e); + } +} + +meshkernel::UInt meshkernel::MeshTriangulation::FindNearestFace(const Point& pnt) const +{ + m_elementCentreRTree->SearchNearestPoint(pnt); + + if (!m_elementCentreRTree->HasQueryResults()) + { + return constants::missing::uintValue; + } + + UInt faceId = m_elementCentreRTree->GetQueryResult(0); + + if (PointIsInElement(pnt, faceId)) + { + return faceId; + } + + const auto edgeIds = GetEdgeIds(faceId); + + BoundedStack<4 * MaximumNumberOfEdgesPerNode> elementsChecked; + elementsChecked.push_back(faceId); + + 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; + } + + elementsChecked.push_back(neighbour1); + } + + if (neighbour2 != faceId && neighbour2 != constants::missing::uintValue) + { + if (PointIsInElement(pnt, neighbour2)) + { + return neighbour2; + } + + elementsChecked.push_back(neighbour2); + } + } + + // Point not in direct neighbours of the element + + m_nodeRTree->SearchNearestPoint(pnt); + + if (m_nodeRTree->HasQueryResults()) + { + UInt nodeId = m_nodeRTree->GetQueryResult(0); + + if (nodeId == constants::missing::uintValue) + { + return constants::missing::uintValue; + } + + for (UInt n = 0; n < m_nodesEdges[nodeId].size(); ++n) + { + const auto faces = GetFaceIds(m_nodesEdges[nodeId][n]); + + if (faces[0] != constants::missing::uintValue && !elementsChecked.contains(faces[0])) + { + + if (PointIsInElement(pnt, faces[0])) + { + return faces[0]; + } + + elementsChecked.push_back(faces[0]); + } + + if (faces[1] != constants::missing::uintValue && !elementsChecked.contains(faces[1])) + { + + if (PointIsInElement(pnt, faces[1])) + { + return faces[1]; + } + + elementsChecked.push_back(faces[1]); + } + } + } + + return constants::missing::uintValue; +} + +std::array meshkernel::MeshTriangulation::GetNodes(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 {GetNode(m_faceNodes[faceIndexStart]), GetNode(m_faceNodes[faceIndexStart + 1]), GetNode(m_faceNodes[faceIndexStart + 2])}; +} + +std::array meshkernel::MeshTriangulation::GetNodeIds(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_faceNodes[faceIndexStart], m_faceNodes[faceIndexStart + 1], m_faceNodes[faceIndexStart + 2]}; +} + +std::array 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) + { + throw ConstraintError("Invalid face id"); + } + + if (faceId >= m_numFaces) + { + throw ConstraintError("face id out of range: {} >= {}", faceId, m_numFaces); + } + + if (!pnt.IsValid()) + { + throw ConstraintError("Point is not valid"); + } + + return IsPointInTriangle(pnt, GetNodes(faceId), m_projection); +} diff --git a/libs/MeshKernel/src/Operations.cpp b/libs/MeshKernel/src/Operations.cpp index 3e3bc6660..95f181be0 100644 --- a/libs/MeshKernel/src/Operations.cpp +++ b/libs/MeshKernel/src/Operations.cpp @@ -286,6 +286,74 @@ namespace meshkernel return dx1 * dy2 - dy1 * dx2; } + bool IsPointInTriangle(const Point& point, + const std::span triangleNodes, + const Projection& projection) + { + if (triangleNodes.empty()) + { + return true; + } + + bool isInTriangle = false; + + if (projection == Projection::cartesian || projection == Projection::spherical) + { + int windingNumber = 0; + + for (UInt n = 0; n < constants::geometric::numNodesInTriangle; n++) + { + UInt endIndex = n == 2 ? 0 : n + 1; + + const auto crossProductValue = crossProduct(triangleNodes[n], triangleNodes[endIndex], triangleNodes[n], point, Projection::cartesian); + + if (IsEqual(crossProductValue, 0.0)) + { + // point on the line + return true; + } + + if (triangleNodes[n].y <= point.y) // an upward crossing + { + if (triangleNodes[endIndex].y > point.y && crossProductValue > 0.0) + { + ++windingNumber; // have a valid up intersect + } + } + else + { + if (triangleNodes[endIndex].y <= point.y && crossProductValue < 0.0) // a downward crossing + { + --windingNumber; // have a valid down intersect + } + } + } + + isInTriangle = windingNumber == 0 ? false : true; + } + + if (projection == Projection::sphericalAccurate) + { + std::vector triangleCopy; + triangleCopy.reserve(constants::geometric::numNodesInTriangle + 1); + Point centre(0.0, 0.0); + + for (UInt i = 0; i < constants::geometric::numNodesInTriangle; ++i) + { + centre += triangleNodes[i]; + triangleCopy.emplace_back(triangleNodes[i]); + } + + triangleCopy.emplace_back(triangleNodes[0]); + + centre *= 1.0 / 3.0; + + isInTriangle = IsPointInPolygonNodes(point, triangleCopy, projection, centre); + } + + return isInTriangle; + } + bool IsPointInPolygonNodes(const Point& point, const std::vector& polygonNodes, const Projection& projection, @@ -335,6 +403,56 @@ namespace meshkernel } } + bool IsPointInPolygonNodes(const Point& point, + const std::vector& polygonNodes, + const Projection& projection, + const BoundingBox& boundingBox, + Point polygonCenter, + UInt startNode, + UInt endNode) + { + + if (polygonNodes.empty()) + { + return true; + } + + if (startNode == constants::missing::uintValue && endNode == constants::missing::uintValue) + { + startNode = 0; + endNode = static_cast(polygonNodes.size()) - 1; // closed polygon + } + + if (endNode <= startNode) + { + return true; + } + + if (const auto currentPolygonSize = endNode - startNode + 1; currentPolygonSize < constants::geometric::numNodesInTriangle || polygonNodes.size() < currentPolygonSize) + { + return false; + } + + if (polygonNodes[startNode] != polygonNodes[endNode]) + { + return false; + } + + if (!boundingBox.Contains(point)) + { + return false; + } + + if (projection == Projection::cartesian || projection == Projection::spherical) + { + return impl::IsPointInPolygonNodesCartesian(point, polygonNodes, startNode, endNode); + } + else + { + return impl::IsPointInPolygonNodesSphericalAccurate(point, polygonNodes, polygonCenter, startNode, endNode); + } + } + void ComputeThreeBaseComponents(const Point& point, std::array& exxp, std::array& eyyp, std::array& ezzp) { const double phi0 = point.y * constants::conversion::degToRad; @@ -1528,40 +1646,47 @@ namespace meshkernel return result; } - Point ComputeAverageCoordinate(const std::vector& points, const Projection& projection) + meshkernel::Point ComputeAverageCoordinate(const std::span points, const Projection& projection) { - std::vector validPoints; - validPoints.reserve(points.size()); - for (const auto& point : points) - { - if (!point.IsValid()) - { - continue; - } - validPoints.emplace_back(point); - } + size_t validCount = std::ranges::count_if(points, [](const Point& p) + { return p.IsValid(); }); if (projection == Projection::sphericalAccurate) { + UInt firstValidPoint = 0; + + if (validCount != points.size()) + { + auto iterator = std::ranges::find_if(points, [](const Point& p) + { return p.IsValid(); }); + firstValidPoint = static_cast(iterator - points.begin()); + } + Cartesian3DPoint averagePoint3D{0.0, 0.0, 0.0}; - for (const auto& point : validPoints) + for (const auto& point : points) { + if (!point.IsValid()) + { + continue; + } + const Cartesian3DPoint point3D{SphericalToCartesian3D(point)}; averagePoint3D.x += point3D.x; averagePoint3D.y += point3D.y; averagePoint3D.z += point3D.z; } - averagePoint3D.x = averagePoint3D.x / static_cast(validPoints.size()); - averagePoint3D.y = averagePoint3D.y / static_cast(validPoints.size()); - averagePoint3D.z = averagePoint3D.z / static_cast(validPoints.size()); + averagePoint3D.x = averagePoint3D.x / static_cast(validCount); + averagePoint3D.y = averagePoint3D.y / static_cast(validCount); + averagePoint3D.z = averagePoint3D.z / static_cast(validCount); - return Cartesian3DToSpherical(averagePoint3D, points[0].x); + return Cartesian3DToSpherical(averagePoint3D, points[firstValidPoint].x); } - auto result = std::accumulate(validPoints.begin(), validPoints.end(), Point{0.0, 0.0}); - result.x = result.x / static_cast(validPoints.size()); - result.y = result.y / static_cast(validPoints.size()); + auto result = std::accumulate(points.begin(), points.end(), Point{0.0, 0.0}, [](const Point& sum, const Point& current) + { return current.IsValid() ? sum + current : sum; }); + result.x = result.x / static_cast(validCount); + result.y = result.y / static_cast(validCount); return result; } diff --git a/libs/MeshKernel/src/SampleAveragingInterpolator.cpp b/libs/MeshKernel/src/SampleAveragingInterpolator.cpp new file mode 100644 index 000000000..2fd911e39 --- /dev/null +++ b/libs/MeshKernel/src/SampleAveragingInterpolator.cpp @@ -0,0 +1,367 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation version 3. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#include "MeshKernel/SampleAveragingInterpolator.hpp" +#include "MeshKernel/Operations.hpp" + +std::vector meshkernel::SampleAveragingInterpolator::CombineCoordinates(const std::span xNodes, + const std::span yNodes) +{ + std::vector result(xNodes.size()); + + for (size_t i = 0; i < xNodes.size(); ++i) + { + result[i].x = xNodes[i]; + result[i].y = yNodes[i]; + } + + return result; +} + +meshkernel::SampleAveragingInterpolator::SampleAveragingInterpolator(const std::span xNodes, + const std::span yNodes, + const Projection projection, + const InterpolationParameters& interpolationParameters) + : m_samplePoints(CombineCoordinates(xNodes, yNodes)), + m_projection(projection), + m_interpolationParameters(interpolationParameters), + m_strategy(averaging::AveragingStrategyFactory::GetAveragingStrategy(interpolationParameters.m_method, + interpolationParameters.m_minimumNumberOfSamples, + projection)), + m_nodeRTree(RTreeFactory::Create(projection)) +{ + m_nodeRTree->BuildTree(m_samplePoints); +} + +meshkernel::SampleAveragingInterpolator::SampleAveragingInterpolator(const std::span nodes, + const Projection projection, + const InterpolationParameters& interpolationParameters) + : m_samplePoints(nodes.begin(), nodes.end()), + m_projection(projection), + m_interpolationParameters(interpolationParameters), + m_strategy(averaging::AveragingStrategyFactory::GetAveragingStrategy(interpolationParameters.m_method, + interpolationParameters.m_minimumNumberOfSamples, + projection)), + m_nodeRTree(RTreeFactory::Create(projection)) +{ + m_nodeRTree->BuildTree(m_samplePoints); +} + +void meshkernel::SampleAveragingInterpolator::Interpolate(const int propertyId, const std::span interpolationNodes, std::span result) const +{ + const std::vector& propertyValues = GetSampleData(propertyId); + std::vector sampleCache; + sampleCache.reserve(100); + + double searchRadiusSquared = m_interpolationParameters.m_absoluteSearchRadius * m_interpolationParameters.m_absoluteSearchRadius; + + for (size_t i = 0; i < interpolationNodes.size(); ++i) + { + m_nodeRTree->SearchPoints(interpolationNodes[i], searchRadiusSquared); + + double resultValue = constants::missing::doubleValue; + + if (!m_nodeRTree->HasQueryResults()) + { + resultValue = constants::missing::doubleValue; + } + else + { + sampleCache.clear(); + + for (UInt j = 0; j < m_nodeRTree->GetQueryResultSize(); ++j) + { + auto const sampleIndex = m_nodeRTree->GetQueryResult(j); + auto const sampleValue = propertyValues[sampleIndex]; + + if (sampleValue == constants::missing::doubleValue) + { + continue; + } + + const Point& samplePoint(m_samplePoints[sampleIndex]); + + sampleCache.emplace_back(samplePoint.x, samplePoint.y, sampleValue); + } + + resultValue = m_strategy->Calculate(interpolationNodes[i], sampleCache); + } + + result[i] = resultValue; + } +} + +double meshkernel::SampleAveragingInterpolator::GetSearchRadiusSquared(const std::vector& searchPolygon, + const Point& interpolationPoint, + const Projection projection) const +{ + double result = std::numeric_limits::lowest(); + + for (const auto& value : searchPolygon) + { + auto const squaredDistance = ComputeSquaredDistance(interpolationPoint, value, projection); + result = std::max(result, squaredDistance); + } + + return result; +} + +void meshkernel::SampleAveragingInterpolator::GenerateSearchPolygon(const double relativeSearchRadius, + const Point& interpolationPoint, + std::vector& polygon, + const Projection projection) const +{ + std::ranges::transform(std::begin(polygon), + std::end(polygon), + std::begin(polygon), + [&](const Point& p) + { return p * relativeSearchRadius + interpolationPoint * (1.0 - relativeSearchRadius); }); + + if (projection == Projection::spherical) + { + const auto boundingBox = BoundingBox(polygon); + const auto lowerLeft = boundingBox.lowerLeft(); + const auto upperRight = boundingBox.upperRight(); + + if (upperRight.x - lowerLeft.x <= 180.0) + { + return; + } + + auto const x_mean = 0.5 * (upperRight.x + lowerLeft.x); + + for (auto& value : polygon) + { + if (value.x < x_mean) + { + value.x = value.x + 360.0; + } + } + } +} + +double meshkernel::SampleAveragingInterpolator::GetSampleValueFromRTree(const int propertyId, const UInt index) const +{ + UInt sampleIndex = m_nodeRTree->GetQueryResult(index); + return GetSampleData(propertyId)[sampleIndex]; +} + +double meshkernel::SampleAveragingInterpolator::ComputeInterpolationResultFromNeighbors(const int propertyId, + const Point& interpolationPoint, + const std::vector& searchPolygon, + const Projection projection, + std::vector& sampleCache) const +{ + sampleCache.clear(); + + const std::vector& propertyData(GetSampleData(propertyId)); + + Point polygonCentre{constants::missing::doubleValue, constants::missing::doubleValue}; + + if (projection == Projection::sphericalAccurate) + { + polygonCentre = ComputeAverageCoordinate(searchPolygon, projection); + } + + BoundingBox boundingBox(searchPolygon); + + for (UInt i = 0; i < m_nodeRTree->GetQueryResultSize(); ++i) + { + auto const sampleIndex = m_nodeRTree->GetQueryResult(i); + auto const sampleValue = propertyData[sampleIndex]; + + if (sampleValue == constants::missing::doubleValue) + { + continue; + } + + const Point& samplePoint(m_samplePoints[sampleIndex]); + + if (IsPointInPolygonNodes(samplePoint, searchPolygon, projection, boundingBox, polygonCentre)) + { + sampleCache.emplace_back(samplePoint.x, samplePoint.y, sampleValue); + } + } + + return m_strategy->Calculate(interpolationPoint, sampleCache); +} + +double meshkernel::SampleAveragingInterpolator::ComputeOnPolygon(const int propertyId, + std::vector& polygon, + const Point& interpolationPoint, + const Projection projection, + std::vector& sampleCache) const +{ + + if (!interpolationPoint.IsValid()) + { + throw ConstraintError("Invalid interpolation point"); + } + + GenerateSearchPolygon(m_interpolationParameters.m_relativeSearchRadius, interpolationPoint, polygon, projection); + + const double searchRadiusSquared = GetSearchRadiusSquared(polygon, interpolationPoint, projection); + + if (searchRadiusSquared <= 0.0) + { + throw ConstraintError("Search radius: {} <= 0", searchRadiusSquared); + } + + m_nodeRTree->SearchPoints(interpolationPoint, searchRadiusSquared); + + if (!m_nodeRTree->HasQueryResults() && m_interpolationParameters.m_useClosestIfNoneFound) + { + m_nodeRTree->SearchNearestPoint(interpolationPoint); + return m_nodeRTree->HasQueryResults() ? GetSampleValueFromRTree(propertyId, 0) : constants::missing::doubleValue; + } + else if (m_nodeRTree->HasQueryResults()) + { + return ComputeInterpolationResultFromNeighbors(propertyId, interpolationPoint, polygon, projection, sampleCache); + } + + return constants::missing::doubleValue; +} + +void meshkernel::SampleAveragingInterpolator::InterpolateAtNodes(const int propertyId, const Mesh2D& mesh, std::span& result) const +{ + std::vector dualFacePolygon; + dualFacePolygon.reserve(MaximumNumberOfEdgesPerNode); + std::vector sampleCache; + sampleCache.reserve(100); + + for (UInt n = 0; n < mesh.GetNumNodes(); ++n) + { + mesh.MakeDualFace(n, m_interpolationParameters.m_relativeSearchRadius, dualFacePolygon); + + double resultValue = constants::missing::doubleValue; + + if (!dualFacePolygon.empty()) + { + resultValue = ComputeOnPolygon(propertyId, + dualFacePolygon, + mesh.Node(n), + mesh.m_projection, + sampleCache); + } + + result[n] = resultValue; + } +} + +void meshkernel::SampleAveragingInterpolator::InterpolateAtEdgeCentres(const Mesh2D& mesh, + const std::span& nodeResult, + std::span& result) const +{ + std::ranges::fill(result, constants::missing::doubleValue); + + for (UInt e = 0; e < mesh.GetNumEdges(); ++e) + { + const auto& [first, second] = mesh.GetEdge(e); + + if (first == constants::missing::uintValue || second == constants::missing::uintValue) + { + continue; + } + + const double& firstValue = nodeResult[first]; + const double& secondValue = nodeResult[second]; + + if (firstValue != constants::missing::doubleValue && secondValue != constants::missing::doubleValue) + { + result[e] = 0.5 * (firstValue + secondValue); + } + } +} + +void meshkernel::SampleAveragingInterpolator::InterpolateAtFaces(const int propertyId, const Mesh2D& mesh, std::span& result) const +{ + std::vector polygonNodesCache(MaximumNumberOfEdgesPerNode + 1); + std::ranges::fill(result, constants::missing::doubleValue); + std::vector sampleCache; + sampleCache.reserve(100); + + for (UInt f = 0; f < mesh.GetNumFaces(); ++f) + { + polygonNodesCache.clear(); + + for (UInt n = 0; n < mesh.GetNumFaceEdges(f); ++n) + { + polygonNodesCache.emplace_back(mesh.m_facesMassCenters[f] + (mesh.Node(mesh.m_facesNodes[f][n]) - mesh.m_facesMassCenters[f]) * m_interpolationParameters.m_relativeSearchRadius); + } + + // Close the polygon + polygonNodesCache.emplace_back(polygonNodesCache[0]); + + sampleCache.clear(); + result[f] = ComputeOnPolygon(propertyId, + polygonNodesCache, + mesh.m_facesMassCenters[f], + mesh.m_projection, + sampleCache); + } +} + +void meshkernel::SampleAveragingInterpolator::Interpolate(const int propertyId, const Mesh2D& mesh, const Location location, std::span result) const +{ + std::ranges::fill(result, constants::missing::doubleValue); + + // Only allocate intermediateNodeResult if needed, i.e. when location is edges + std::vector intermediateNodeResult; + std::span nodeResult; + + using enum Location; + + if (location == Nodes) + { + nodeResult = std::span(result.data(), result.size()); + } + else if (location == Edges) + { + intermediateNodeResult.resize(mesh.GetNumNodes(), 0.0); + nodeResult = std::span(intermediateNodeResult.data(), intermediateNodeResult.size()); + } + + if (location == Nodes || location == Edges) + { + InterpolateAtNodes(propertyId, mesh, nodeResult); + } + + if (location == Edges) + { + InterpolateAtEdgeCentres(mesh, nodeResult, result); + } + + if (location == Faces) + { + InterpolateAtFaces(propertyId, mesh, result); + } +} + +double meshkernel::SampleAveragingInterpolator::InterpolateValue(const int propertyId [[maybe_unused]], const Point& evaluationPoint [[maybe_unused]]) const +{ + return constants::missing::doubleValue; +} diff --git a/libs/MeshKernel/src/SampleInterpolator.cpp b/libs/MeshKernel/src/SampleInterpolator.cpp new file mode 100644 index 000000000..0464c50e7 --- /dev/null +++ b/libs/MeshKernel/src/SampleInterpolator.cpp @@ -0,0 +1,51 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation version 3. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#include "MeshKernel/SampleInterpolator.hpp" +#include "MeshKernel/Exceptions.hpp" + +void meshkernel::SampleInterpolator::SetData(const int propertyId, const std::span sampleData) +{ + if (Size() != sampleData.size()) + { + throw ConstraintError("The sample data array does not have the same number of values as the sample point set: {} /= {}", + Size(), sampleData.size()); + } + + m_sampleData[propertyId].assign(sampleData.begin(), sampleData.end()); +} + +const std::vector& meshkernel::SampleInterpolator::GetSampleData(const int propertyId) const +{ + + if (!Contains(propertyId)) + { + throw ConstraintError("The sample data set for property id {}, has not been defined", propertyId); + } + + return m_sampleData.at(propertyId); +} diff --git a/libs/MeshKernel/src/SampleTriangulationInterpolator.cpp b/libs/MeshKernel/src/SampleTriangulationInterpolator.cpp new file mode 100644 index 000000000..40b005267 --- /dev/null +++ b/libs/MeshKernel/src/SampleTriangulationInterpolator.cpp @@ -0,0 +1,136 @@ +#include "MeshKernel/SampleTriangulationInterpolator.hpp" + +void meshkernel::SampleTriangulationInterpolator::Interpolate(const int propertyId, const Mesh2D& mesh, const Location location, std::span result) const +{ + if (!Contains(propertyId)) + { + throw ConstraintError("Sample interpolator does not contain the id: {}.", propertyId); + } + + std::span meshNodes; + + switch (location) + { + case Location::Nodes: + meshNodes = std::span(mesh.Nodes()); + break; + case Location::Edges: + meshNodes = std::span(mesh.m_edgesCenters); + break; + case Location::Faces: + meshNodes = std::span(mesh.m_facesMassCenters); + break; + default: + throw ConstraintError("Unknown location"); + } + + Interpolate(propertyId, meshNodes, result); +} + +void meshkernel::SampleTriangulationInterpolator::Interpolate(const int propertyId, const std::span interpolationNodes, std::span result) const +{ + if (!Contains(propertyId)) + { + throw ConstraintError("Sample interpolator does not contain the id: {}.", propertyId); + } + + if (interpolationNodes.size() != result.size()) + { + throw ConstraintError("The arrays for interpolation nodes and the results are different sizes: {} /= {}", + interpolationNodes.size(), result.size()); + } + + const std::vector& propertyValues = GetSampleData(propertyId); + + for (size_t i = 0; i < interpolationNodes.size(); ++i) + { + result[i] = constants::missing::doubleValue; + + if (!interpolationNodes[i].IsValid()) + { + continue; + } + + const UInt elementId = m_triangulation.FindNearestFace(interpolationNodes[i]); + + if (elementId == constants::missing::uintValue) + { + continue; + } + + if (m_triangulation.PointIsInElement(interpolationNodes[i], elementId)) + { + result[i] = InterpolateOnElement(elementId, interpolationNodes[i], propertyValues); + } + } +} + +double meshkernel::SampleTriangulationInterpolator::InterpolateOnElement(const UInt elementId, const Point& interpolationPoint, const std::vector& sampleValues) const +{ + double result = constants::missing::doubleValue; + + auto [id1, id2, id3] = m_triangulation.GetNodeIds(elementId); + + if (sampleValues[id1] == constants::missing::doubleValue || + sampleValues[id2] == constants::missing::doubleValue || + sampleValues[id3] == constants::missing::doubleValue) [[unlikely]] + { + return result; + } + + auto [p1, p2, p3] = m_triangulation.GetNodes(elementId); + + const double a11 = GetDx(p1, p2, m_triangulation.GetProjection()); + const double a21 = GetDy(p1, p2, m_triangulation.GetProjection()); + + const double a12 = GetDx(p1, p3, m_triangulation.GetProjection()); + const double a22 = GetDy(p1, p3, m_triangulation.GetProjection()); + + const double b1 = GetDx(p1, interpolationPoint, m_triangulation.GetProjection()); + const double b2 = GetDy(p1, interpolationPoint, m_triangulation.GetProjection()); + + const double det = a11 * a22 - a12 * a21; + + if (std::abs(det) < 1e-12) [[unlikely]] + { + return result; + } + + const double rlam = (a22 * b1 - a12 * b2) / det; + const double rmhu = (a11 * b2 - a21 * b1) / det; + + result = sampleValues[id1] + rlam * (sampleValues[id2] - sampleValues[id1]) + rmhu * (sampleValues[id3] - sampleValues[id1]); + + return result; +} + +double meshkernel::SampleTriangulationInterpolator::InterpolateValue(const int propertyId, const Point& evaluationPoint) const +{ + + if (!Contains(propertyId)) + { + throw ConstraintError("Sample interpolator does not contain the id: {}.", propertyId); + } + + double result = constants::missing::doubleValue; + + if (!evaluationPoint.IsValid()) + { + return result; + } + + const UInt elementId = m_triangulation.FindNearestFace(evaluationPoint); + + if (elementId == constants::missing::uintValue) + { + return result; + } + + if (m_triangulation.PointIsInElement(evaluationPoint, elementId)) + { + const std::vector& propertyValues = GetSampleData(propertyId); + result = InterpolateOnElement(elementId, evaluationPoint, propertyValues); + } + + return result; +} diff --git a/libs/MeshKernel/src/Utilities/Utilities.cpp b/libs/MeshKernel/src/Utilities/Utilities.cpp index 5d2b34ee9..b5ff32b7e 100644 --- a/libs/MeshKernel/src/Utilities/Utilities.cpp +++ b/libs/MeshKernel/src/Utilities/Utilities.cpp @@ -28,6 +28,9 @@ #include "MeshKernel/Utilities/Utilities.hpp" #include "MeshKernel/Exceptions.hpp" +#include +#include + void meshkernel::Print(const std::vector& nodes, const std::vector& edges, std::ostream& out) { out << "nullId = " << constants::missing::uintValue << ";" << std::endl; @@ -90,3 +93,142 @@ void meshkernel::Print(const std::vector& xNodes, out << "edges ( " << i + 1 << ", 2 ) = " << edges[2 * i + 1] + 1 << ";" << std::endl; } } + +void meshkernel::SaveVtk(const std::vector& nodes, const std::vector>& faces, const std::string& fileName) +{ + + std::ofstream vtkFile(fileName.c_str()); + + vtkFile.precision(18); + + std::string meshType = "UnstructuredGrid"; + std::string versionNumber = "1.0"; + + UInt numberOfElements = 0; + + for (size_t i = 0; i < faces.size(); ++i) + { + if (faces[i].size() == 3 || faces[i].size() == 4) + { + ++numberOfElements; + } + } + + vtkFile << "" + << std::endl; + vtkFile << " " << std::endl; + + vtkFile << " " + << std::endl; + + vtkFile << " " << std::endl; + vtkFile << " " << std::endl; + + for (size_t i = 0; i < nodes.size(); ++i) + { + vtkFile << std::setw(8) << " " << nodes[i].x << " " << std::setw(8) << nodes[i].y << " " << 0.0 << std::endl; + } + + vtkFile << " " << std::endl; + vtkFile << " " << std::endl; + + //-------------------------------- + vtkFile << " " << std::endl; + vtkFile << " " + << std::endl; + + for (size_t i = 0; i < faces.size(); ++i) + { + + if (faces[i].size() == 3 || faces[i].size() == 4) + { + for (size_t j = 0; j < faces[i].size(); ++j) + { + vtkFile << " " << std::setw(8) << faces[i][j]; + } + + vtkFile << std::endl; + } + } + + vtkFile << std::endl; + vtkFile << " " << std::endl; + + //-------------------------------- + vtkFile << " " + << std::endl; + + size_t sum = 0; + + for (size_t i = 0; i < faces.size(); ++i) + { + if (faces[i].size() == 3 || faces[i].size() == 4) + { + sum += faces[i].size(); + vtkFile << " " << std::setw(8) << sum; + } + } + + vtkFile << std::endl; + vtkFile << " " << std::endl; + + //-------------------------------- + vtkFile << " " + << std::endl; + + for (size_t i = 0; i < faces.size(); ++i) + { + UInt type = 0; + + if (faces[i].size() == 3) + { + type = 69; + } + else if (faces[i].size() == 4) + { + type = 70; + } + + if (faces[i].size() == 3 || faces[i].size() == 4) + { + vtkFile << " " << std::setw(8) << type; + } + } + + vtkFile << std::endl; + vtkFile << " " << std::endl; + vtkFile << " " << std::endl; + + vtkFile << " " << std::endl; + vtkFile << " " << std::endl; + vtkFile << "" << std::endl; + + vtkFile.close(); +} diff --git a/libs/MeshKernel/tests/CMakeLists.txt b/libs/MeshKernel/tests/CMakeLists.txt index 112150fdb..81cca6b35 100644 --- a/libs/MeshKernel/tests/CMakeLists.txt +++ b/libs/MeshKernel/tests/CMakeLists.txt @@ -43,6 +43,7 @@ set( ${SRC_DIR}/Mesh2DTest.cpp ${SRC_DIR}/MeshConversionTests.cpp ${SRC_DIR}/MeshConversionTests.cpp + ${SRC_DIR}/MeshPropertyTests.cpp ${SRC_DIR}/MeshRefinementTests.cpp ${SRC_DIR}/MeshTests.cpp ${SRC_DIR}/MeshTransformationTest.cpp @@ -53,6 +54,7 @@ set( ${SRC_DIR}/PolygonalEnclosureTests.cpp ${SRC_DIR}/PolygonsTests.cpp ${SRC_DIR}/RangeCheckTests.cpp + ${SRC_DIR}/SampleInterpolationTests.cpp ${SRC_DIR}/SpatialTreesTests.cpp ${SRC_DIR}/SplineAlgorithmTests.cpp ${SRC_DIR}/SplineTests.cpp diff --git a/libs/MeshKernel/tests/src/MeshPropertyTests.cpp b/libs/MeshKernel/tests/src/MeshPropertyTests.cpp new file mode 100644 index 000000000..6ebb10e09 --- /dev/null +++ b/libs/MeshKernel/tests/src/MeshPropertyTests.cpp @@ -0,0 +1,156 @@ +#include +#include +#include +#include +#include + +#include "MeshKernel/Constants.hpp" +#include "MeshKernel/Entities.hpp" +#include "MeshKernel/MeshTriangulation.hpp" +#include "MeshKernel/SampleTriangulationInterpolator.hpp" + +namespace mk = meshkernel; + +TEST(MeshPropertyTests, TriangulationTest) +{ + std::vector xValues{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}; + std::vector yValues{0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0}; + + std::vector nodes{{0.25, 0.25}, {1.25, 0.25}, {2.25, 0.25}, {0.25, 1.25}, {1.25, 1.25}, {2.25, 1.25}, {0.25, 2.25}, {1.25, 2.25}, {2.25, 2.25}}; + std::vector interpolated(nodes.size(), 0.0); + + mk::MeshTriangulation triangulation(xValues, yValues, mk::Projection::cartesian); + + EXPECT_EQ(triangulation.GetProjection(), mk::Projection::cartesian); + ASSERT_EQ(triangulation.NumberOfNodes(), 16); + ASSERT_EQ(triangulation.NumberOfEdges(), 33); + ASSERT_EQ(triangulation.NumberOfFaces(), 18); + + std::vector expectedEdgeIds{{4, 0}, {0, 1}, {1, 4}, {1, 5}, {5, 4}, {1, 2}, {2, 5}, {5, 9}, {9, 4}, {12, 8}, {8, 9}, {9, 12}, {9, 13}, {13, 12}, {9, 10}, {10, 13}, {8, 4}, {5, 6}, {6, 9}, {2, 6}, {2, 3}, {3, 6}, {3, 7}, {7, 6}, {7, 11}, {11, 6}, {10, 14}, {14, 13}, {10, 11}, {11, 14}, {11, 15}, {15, 14}, {10, 6}}; + std::vector> expectedTriangleNodeIds{{4, 0, 1}, {1, 5, 4}, {5, 1, 2}, {4, 5, 9}, {12, 8, 9}, {9, 13, 12}, {13, 9, 10}, {8, 4, 9}, {5, 6, 9}, {2, 6, 5}, {6, 2, 3}, {3, 7, 6}, {6, 7, 11}, {10, 14, 13}, {14, 10, 11}, {11, 15, 14}, {10, 6, 11}, {6, 10, 9}}; + + for (mk::UInt i = 0; i < triangulation.NumberOfEdges(); ++i) + { + auto [edgeStart, edgeEnd] = triangulation.GetEdge(i); + EXPECT_EQ(edgeStart, expectedEdgeIds[i].first); + EXPECT_EQ(edgeEnd, expectedEdgeIds[i].second); + } + + for (mk::UInt i = 0; i < triangulation.NumberOfFaces(); ++i) + { + auto [id1, id2, id3] = triangulation.GetNodeIds(i); + EXPECT_EQ(id1, expectedTriangleNodeIds[i][0]); + EXPECT_EQ(id2, expectedTriangleNodeIds[i][1]); + EXPECT_EQ(id3, expectedTriangleNodeIds[i][2]); + } + + for (mk::UInt i = 0; i < triangulation.NumberOfFaces(); ++i) + { + auto [p1, p2, p3] = triangulation.GetNodes(i); + EXPECT_EQ(p1.x, triangulation.GetNode(expectedTriangleNodeIds[i][0]).x); + EXPECT_EQ(p1.y, triangulation.GetNode(expectedTriangleNodeIds[i][0]).y); + + EXPECT_EQ(p2.x, triangulation.GetNode(expectedTriangleNodeIds[i][1]).x); + EXPECT_EQ(p2.y, triangulation.GetNode(expectedTriangleNodeIds[i][1]).y); + + EXPECT_EQ(p3.x, triangulation.GetNode(expectedTriangleNodeIds[i][2]).x); + EXPECT_EQ(p3.y, triangulation.GetNode(expectedTriangleNodeIds[i][2]).y); + } + + // Origin of mesh is {0, 0}, with delta {1,1} and there are 3x3 nodes + + std::vector pointsToCheck{{-1.0, -1.0}, {0.25, 0.25}, {0.25, 2.25}, {3.5, 0.25}, {4.0, 4.0}, {2.5, 2.5}}; + std::vector expectedElementId{mk::constants::missing::uintValue, + 0, 4, mk::constants::missing::uintValue, mk::constants::missing::uintValue, 15}; + std::vector expectedIsInElement{false, true, true, false, false, true}; + + for (size_t i = 0; i < pointsToCheck.size(); ++i) + { + mk::UInt faceId = triangulation.FindNearestFace(pointsToCheck[i]); + EXPECT_EQ(faceId, expectedElementId[i]); + + if (faceId != mk::constants::missing::uintValue) + { + bool isInFace = triangulation.PointIsInElement(pointsToCheck[i], faceId); + EXPECT_EQ(isInFace, expectedIsInElement[i]); + } + } +} + +TEST(MeshPropertyTests, TriangulationFailureTest) +{ + // Not enough points + std::vector twoPointsX{0.0, 1.0}; + std::vector twoPointsY{0.0, 1.0}; + EXPECT_THROW(mk::MeshTriangulation triangulation(twoPointsX, twoPointsY, mk::Projection::cartesian), mk::ConstraintError); + // x- and y-points not same size + std::vector threePointsX{0.0, 1.0, 2.0}; + EXPECT_THROW(mk::MeshTriangulation triangulation(threePointsX, twoPointsY, mk::Projection::cartesian), mk::ConstraintError); + + std::vector xValues{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}; + std::vector yValues{0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0}; + mk::MeshTriangulation triangulation(xValues, yValues, mk::Projection::cartesian); + + // Invalid point passed to PointIsInElement + EXPECT_THROW(triangulation.PointIsInElement({mk::constants::missing::doubleValue, mk::constants::missing::doubleValue}, 1), mk::ConstraintError); + // Face id is null + EXPECT_THROW(triangulation.PointIsInElement({0.25, 0.25}, mk::constants::missing::uintValue), mk::ConstraintError); + // Face id is out of range + EXPECT_THROW(triangulation.PointIsInElement({0.25, 0.25}, 100), mk::ConstraintError); +} + +TEST(MeshPropertyTests, SimpleSampleInterpolationTest) +{ + const double tolerance = 1.0e-13; + + const int xDepth = 1; + const int yDepth = 2; + + std::vector xValues{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}; + std::vector yValues{0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0}; + + std::vector nodes{{0.25, 0.25}, {1.25, 0.25}, {2.25, 0.25}, {0.25, 1.25}, {1.25, 1.25}, {2.25, 1.25}, {0.25, 2.25}, {1.25, 2.25}, {2.25, 2.25}}; + std::vector interpolated(nodes.size(), 0.0); + + mk::SampleTriangulationInterpolator properties(xValues, yValues, mk::Projection::cartesian); + + std::vector expectedXDepths{0.25, 1.25, 2.25, 0.25, 1.25, 2.25, 0.25, 1.25, 2.25}; + std::vector expectedYDepths{0.25, 0.25, 0.25, 1.25, 1.25, 1.25, 2.25, 2.25, 2.25}; + + properties.SetData(xDepth, std::vector{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}); + + properties.SetData(yDepth, std::vector{0.0, 0.0, 0.0, 0.0, + 1.0, 1.0, 1.0, 1.0, + 2.0, 2.0, 2.0, 2.0, + 3.0, 3.0, 3.0, 3.0}); + + properties.Interpolate(xDepth, nodes, interpolated); + + for (size_t i = 0; i < interpolated.size(); ++i) + { + EXPECT_NEAR(interpolated[i], expectedXDepths[i], tolerance); + } + + std::span interpolatedData(interpolated.data(), interpolated.size()); + + properties.Interpolate(yDepth, nodes, interpolatedData); + + for (size_t i = 0; i < interpolated.size(); ++i) + { + EXPECT_NEAR(interpolatedData[i], expectedYDepths[i], tolerance); + } +} + +TEST(MeshPropertyTests, AveragePointTest) +{ + const double tolerance = 1.0e-13; + + std::vector pnts{{1.0, 2.0}, {-999.0, -999.0}, {1.0, 2.0}, {1.0, 2.0}}; + mk::Point average = mk::ComputeAverageCoordinate(pnts, mk::Projection::sphericalAccurate); + + EXPECT_NEAR(average.x, 1.0, tolerance); + EXPECT_NEAR(average.y, 2.0, tolerance); +} diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index afb2da8a4..012b23c60 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -25,6 +25,9 @@ // //------------------------------------------------------------------------------ +#include +#include + #include #include @@ -34,21 +37,25 @@ #include "MeshKernel/CasulliRefinement.hpp" #include "MeshKernel/Mesh2D.hpp" #include "MeshKernel/MeshRefinement.hpp" +#include "MeshKernel/Operations.hpp" #include "MeshKernel/Parameters.hpp" #include "MeshKernel/Polygons.hpp" +#include "MeshKernel/RemoveDisconnectedRegions.hpp" +#include "MeshKernel/SampleAveragingInterpolator.hpp" +#include "MeshKernel/SampleInterpolator.hpp" +#include "MeshKernel/SampleTriangulationInterpolator.hpp" #include "MeshKernel/SamplesHessianCalculator.hpp" #include "MeshKernel/SplitRowColumnOfMesh.hpp" #include "MeshKernel/UndoActions/UndoActionStack.hpp" +#include "MeshKernel/Utilities/Utilities.hpp" + #include "TestUtils/Definitions.hpp" +#include "TestUtils/MakeCurvilinearGrids.hpp" #include "TestUtils/MakeMeshes.hpp" #include "TestUtils/MeshReaders.hpp" #include "TestUtils/SampleFileReader.hpp" #include "TestUtils/SampleGenerator.hpp" -#include - -#include - using namespace meshkernel; namespace mk = meshkernel; @@ -1674,31 +1681,25 @@ TEST(MeshRefinement, CasulliRefinement) {12.5, 20}, {17.5, 20}}; - std::vector expectedEdgesStart = {9, 25, 9, 12, 13, - 10, 13, 16, 27, 14, - 27, 28, 17, 29, 17, - 20, 21, 18, 21, 24, - 31, 22, 31, 32, 33, - 9, 33, 34, 35, 13, - 35, 36, 12, 17, 12, - 11, 16, 21, 16, 15, - 20, 37, 20, 19, 24, - 39, 24, 23, 0, 0, - 34, 2, 2, 26, 28, - 6, 6, 39, 8, 8}; - - std::vector expectedEdgesEnd = {12, 26, 25, 26, 16, - 11, 10, 11, 28, 15, - 14, 15, 20, 30, 29, - 30, 24, 19, 18, 19, - 32, 23, 22, 23, 34, - 10, 9, 10, 36, 14, - 13, 14, 11, 18, 17, - 18, 15, 22, 21, 22, - 19, 38, 37, 38, 23, - 40, 39, 40, 25, 33, - 35, 27, 36, 29, 31, - 30, 37, 38, 32, 40}; + std::vector expectedEdgesStart = {4, 20, 4, 7, 8, 5, 8, + 11, 22, 9, 22, 23, 12, + 24, 12, 15, 16, 13, 16, + 19, 26, 17, 26, 27, 28, + 4, 28, 29, 30, 8, 30, 31, + 7, 12, 7, 6, 11, 16, 11, + 10, 15, 32, 15, 14, 19, + 34, 19, 18, 0, 0, 29, 1, + 1, 21, 23, 2, 2, 34, 3, 3}; + + std::vector expectedEdgesEnd = {7, 21, 20, 21, 11, 6, 5, + 6, 23, 10, 9, 10, 15, 25, + 24, 25, 19, 14, 13, 14, 27, + 18, 17, 18, 29, 5, 4, 5, 31, + 9, 8, 9, 6, 13, 12, 13, 10, + 17, 16, 17, 14, 33, 32, 33, + 18, 35, 34, 35, 20, 28, 30, + 22, 31, 24, 26, 25, 32, 33, + 27, 35}; auto undoAction = CasulliRefinement::Compute(mesh); diff --git a/libs/MeshKernel/tests/src/SampleInterpolationTests.cpp b/libs/MeshKernel/tests/src/SampleInterpolationTests.cpp new file mode 100644 index 000000000..ed1f4d914 --- /dev/null +++ b/libs/MeshKernel/tests/src/SampleInterpolationTests.cpp @@ -0,0 +1,240 @@ +#include +#include + +#include + +#include "MeshKernel/SampleAveragingInterpolator.hpp" +#include "MeshKernel/SampleTriangulationInterpolator.hpp" +#include "TestUtils/Definitions.hpp" +#include "TestUtils/MakeMeshes.hpp" +#include "TestUtils/SampleFileReader.hpp" + +namespace mk = meshkernel; + +TEST(SampleInterpolationTests, AveragingInterpolationWithPoints) +{ + const mk::UInt numberOfPointsX = 11; + const mk::UInt numberOfPointsY = 11; + const mk::UInt numberOfPoints = numberOfPointsX * numberOfPointsY; + + std::vector xPoints(numberOfPoints); + std::vector yPoints(numberOfPoints); + std::vector data(numberOfPoints); + + // Generate sample data points + double delta = 1000.0; + double x = 0.0; + double y = 0.0; + size_t count = 0; + + for (size_t i = 0; i < numberOfPointsY; ++i) + { + x = 0.0; + + for (size_t j = 0; j < numberOfPointsX; ++j) + { + xPoints[count] = x; + yPoints[count] = y; + data[count] = x; + x += delta; + ++count; + } + + y += delta; + } + + mk::InterpolationParameters params{.m_absoluteSearchRadius = 3.0 * delta}; + mk::SampleAveragingInterpolator interpolator(xPoints, yPoints, mk::Projection::cartesian, params); + + ASSERT_EQ(static_cast(interpolator.Size()), xPoints.size()); + + int propertyId = 1; + interpolator.SetData(propertyId, data); + + // Execute + ASSERT_EQ(interpolator.Size(), numberOfPoints); + + std::vector interpolationPoints{{0.5 * delta, 0.5 * delta}, {9.5 * delta, 9.5 * delta}}; + const double initialValue = -1.0e20; + std::vector interpolationResult(interpolationPoints.size(), initialValue); + std::vector expectedResult{1400.0, 8600.0}; + + interpolator.Interpolate(propertyId, interpolationPoints, interpolationResult); + + const double tolerance = 1.0e-8; + + for (size_t i = 0; i < expectedResult.size(); ++i) + { + EXPECT_NEAR(expectedResult[i], interpolationResult[i], tolerance); + } +} + +TEST(SampleInterpolationTests, AveragingInterpolationWithMesh) +{ + const mk::UInt numberOfPointsX = 11; + const mk::UInt numberOfPointsY = 11; + const mk::UInt numberOfPoints = numberOfPointsX * numberOfPointsY; + + std::vector samplePoints(numberOfPoints); + std::vector data(numberOfPoints); + + // Generate sample data points + double delta = 1000.0; + double x = 0.0; + double y = 0.0; + size_t count = 0; + + for (size_t i = 0; i < numberOfPointsY; ++i) + { + x = 0.0; + + for (size_t j = 0; j < numberOfPointsX; ++j) + { + samplePoints[i].x = x; + samplePoints[i].y = y; + data[count] = x; + x += delta; + ++count; + } + + y += delta; + } + + mk::InterpolationParameters params{.m_absoluteSearchRadius = 3.0 * delta, .m_minimumNumberOfSamples = 1}; + mk::SampleAveragingInterpolator interpolator(samplePoints, mk::Projection::cartesian, params); + + ASSERT_EQ(static_cast(interpolator.Size()), samplePoints.size()); + + //-------------------------------- + + const mk::UInt meshPointsX = 5; + const mk::UInt meshPointsY = 5; + + const auto mesh = MakeRectangularMeshForTesting(meshPointsX, + meshPointsY, + 8.0 * delta, + 8.0 * delta, + mk::Projection::cartesian, + {0.5 * delta, 0.5 * delta}); + + mesh->ComputeEdgesCenters(); + + //-------------------------------- + + int propertyId = 1; + interpolator.SetData(propertyId, data); + + // Execute + ASSERT_EQ(interpolator.Size(), numberOfPoints); + + // Test node data + + const double initialValue = -1.0e20; + std::vector interpolationNodeResult(mesh->GetNumNodes(), initialValue); + std::vector expectedNodeResult{0.0, 2000.0, 4000.0, 6000.0, 8000.0, + 0.0, 2000.0, 4000.0, 6000.0, 8000.0, + 0.0, 2000.0, 4000.0, 6000.0, 8000.0, + 0.0, 2000.0, 4000.0, 6000.0, 8000.0, + 0.0, 2000.0, 4000.0, 6000.0, 8000.0}; + + interpolator.Interpolate(propertyId, *mesh, mk::Location::Nodes, interpolationNodeResult); + + const double tolerance = 1.0e-8; + + for (size_t i = 0; i < interpolationNodeResult.size(); ++i) + { + EXPECT_NEAR(expectedNodeResult[i], interpolationNodeResult[i], tolerance); + } + + // Test edge data + + std::vector interpolationEdgeResult(mesh->GetNumEdges(), initialValue); + std::vector expectedEdgeResult{0.0, 2000.0, 4000.0, 6000.0, 8000.0, + 0.0, 2000.0, 4000.0, 6000.0, 8000.0, + 0.0, 2000.0, 4000.0, 6000.0, 8000.0, + 0.0, 2000.0, 4000.0, 6000.0, 8000.0, + 1000.0, 3000.0, 5000.0, 7000.0, + 1000.0, 3000.0, 5000.0, 7000.0, + 1000.0, 3000.0, 5000.0, 7000.0, + 1000.0, 3000.0, 5000.0, 7000.0, + 1000.0, 3000.0, 5000.0, 7000.0}; + + interpolator.Interpolate(propertyId, *mesh, mk::Location::Edges, interpolationEdgeResult); + + for (size_t i = 0; i < interpolationEdgeResult.size(); ++i) + { + EXPECT_NEAR(expectedEdgeResult[i], interpolationEdgeResult[i], tolerance); + } + + // Test face data + + std::vector interpolationFaceResult(mesh->GetNumFaces(), initialValue); + std::vector expectedFaceResult{1000.0, 3000.0, 5000.0, 7000.0, + 1000.0, 3000.0, 5000.0, 7000.0, + 1000.0, 3000.0, 5000.0, 7000.0, + 1000.0, 3000.0, 5000.0, 7000.0}; + + interpolator.Interpolate(propertyId, *mesh, mk::Location::Faces, interpolationFaceResult); + + for (size_t i = 0; i < interpolationFaceResult.size(); ++i) + { + EXPECT_NEAR(expectedFaceResult[i], interpolationFaceResult[i], tolerance); + } +} + +TEST(SampleInterpolationTests, TriangulationInterpolationWithPoints) +{ + const mk::UInt numberOfPointsX = 11; + const mk::UInt numberOfPointsY = 11; + const mk::UInt numberOfPoints = numberOfPointsX * numberOfPointsY; + + std::vector xPoints(numberOfPoints); + std::vector yPoints(numberOfPoints); + std::vector data(numberOfPoints); + + // Generate sample data points + double delta = 1000.0; + double x = 0.0; + double y = 0.0; + size_t count = 0; + + for (size_t i = 0; i < numberOfPointsY; ++i) + { + x = 0.0; + + for (size_t j = 0; j < numberOfPointsX; ++j) + { + xPoints[count] = x; + yPoints[count] = y; + data[count] = x; + x += delta; + ++count; + } + + y += delta; + } + + mk::SampleTriangulationInterpolator interpolator(xPoints, yPoints, mk::Projection::cartesian); + + ASSERT_EQ(static_cast(interpolator.Size()), xPoints.size()); + + int propertyId = 1; + interpolator.SetData(propertyId, data); + + // Execute + ASSERT_EQ(interpolator.Size(), numberOfPoints); + + std::vector interpolationPoints{{0.5 * delta, 0.5 * delta}, {9.5 * delta, 9.5 * delta}, {11.0 * delta, 11.0 * delta}}; + const double initialValue = -1.0e20; + std::vector interpolationResult(interpolationPoints.size(), initialValue); + std::vector expectedResult{500.0, 9500.0, mk::constants::missing::doubleValue}; + + interpolator.Interpolate(propertyId, interpolationPoints, interpolationResult); + + const double tolerance = 1.0e-8; + + for (size_t i = 0; i < interpolationResult.size(); ++i) + { + EXPECT_NEAR(expectedResult[i], interpolationResult[i], tolerance); + } +} diff --git a/libs/MeshKernelApi/CMakeLists.txt b/libs/MeshKernelApi/CMakeLists.txt index 9a605c463..e18ab5b76 100644 --- a/libs/MeshKernelApi/CMakeLists.txt +++ b/libs/MeshKernelApi/CMakeLists.txt @@ -25,6 +25,7 @@ set(VERSION_INC_DIR ${CMAKE_SOURCE_DIR}/tools) set(SRC_LIST ${SRC_DIR}/MKStateUndoAction.cpp ${SRC_DIR}/MeshKernel.cpp + ${SRC_DIR}/PropertyCalculator.cpp ) set(CACHE_SRC_LIST @@ -52,6 +53,7 @@ set( ${DOMAIN_INC_DIR}/Mesh1D.hpp ${DOMAIN_INC_DIR}/Mesh2D.hpp ${DOMAIN_INC_DIR}/MeshKernel.hpp + ${DOMAIN_INC_DIR}/PropertyCalculator.hpp ${DOMAIN_INC_DIR}/State.hpp ${DOMAIN_INC_DIR}/Utils.hpp ${VERSION_INC_DIR}/Version/Version.hpp diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp index c2b65ec6b..18517f5aa 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp @@ -685,6 +685,11 @@ namespace meshkernelapi double regionControlPointX = mkernel_get_separator(), double regionControlPointY = mkernel_get_separator()); + /// @brief Deallocate property calculator + /// @param[in] propertyId The id of the property + /// @returns Error code + MKERNEL_API int mkernel_deallocate_property(int propertyId); + /// @brief Deallocate mesh state /// @param[in] meshKernelId The id of the mesh state /// @returns Error code @@ -1344,15 +1349,16 @@ namespace meshkernelapi /// /// @param[in] meshKernelId The id of the mesh state /// @param[in] propertyValue The value representing the specific property - /// @param[in] geometrylist A reference to a GeometryList object that will be populated with the values of the requested property + /// @param[in] locationId The location (nodes, edge centres or face centres) at which the samples should interpolated. + /// @param[in,out] geometrylist A reference to a GeometryList object that will be populated with the values of the requested property /// @returns Error code - MKERNEL_API int mkernel_mesh2d_get_property(int meshKernelId, int propertyValue, const GeometryList& geometrylist); + MKERNEL_API int mkernel_mesh2d_get_property(int meshKernelId, int propertyValue, int locationId, const GeometryList& geometrylist); /// @brief The dimension of a specified property of a 2D mesh. /// /// @param[in] meshKernelId The id of the mesh state /// @param[in] propertyValue The value representing the specific property - /// @param[in] dimension The dimension of the specified property + /// @param[out] dimension The dimension of the specified property /// @returns Error code MKERNEL_API int mkernel_mesh2d_get_property_dimension(int meshKernelId, int propertyValue, int& dimension); @@ -1445,6 +1451,14 @@ namespace meshkernelapi int* faceNumEdges, int* faceEdgeIndex); + /// @brief Determine if the property data for the mesh can be computed + /// @param[in] meshKernelId The id of the mesh state + /// @param[in] propertyId The id of the property + /// @param[in] locationId The location (nodes, edge centres or face centres) at which the samples should interpolated. + /// @param[out] propertyIsAvailable Indicate (true or false) if the property can be calculated. + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_is_valid_property(int meshKernelId, const int propertyId, const int locationId, bool& propertyIsAvailable); + /// @brief Compute the global mesh with a given number of points along the longitude and latitude directions. /// @param[in] meshKernelId The id of the mesh state /// @param [in] numLongitudeNodes The number of points along the longitude. @@ -1595,6 +1609,14 @@ namespace meshkernelapi /// @returns Error code MKERNEL_API int mkernel_mesh2d_set(int meshKernelId, const Mesh2D& mesh2d); + /// @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, int interpolationType, const GeometryList& sampleData, int& propertyId); + /// @brief Snaps a mesh to a land boundary. /// @param[in] meshKernelId The id of the mesh state /// @param[in] selectingPolygon The polygon where to perform the snapping diff --git a/libs/MeshKernelApi/include/MeshKernelApi/PropertyCalculator.hpp b/libs/MeshKernelApi/include/MeshKernelApi/PropertyCalculator.hpp new file mode 100644 index 000000000..5ba7f90fb --- /dev/null +++ b/libs/MeshKernelApi/include/MeshKernelApi/PropertyCalculator.hpp @@ -0,0 +1,122 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation version 3. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include "MeshKernel/Definitions.hpp" +#include "MeshKernel/SampleInterpolator.hpp" + +#include "MeshKernelApi/GeometryList.hpp" +#include "MeshKernelApi/State.hpp" + +namespace meshkernelapi +{ + + /// @brief Base class for calculating properties for a mesh + class PropertyCalculator + { + public: + /// @brief Destructor + virtual ~PropertyCalculator() = default; + + /// @brief Determine is the calculator can compute the desired results correctly. + virtual bool IsValid(const MeshKernelState& state, const meshkernel::Location location) const = 0; + + /// @brief Calculate the property + 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 meshkernel::Location location) const = 0; + }; + + /// @brief Calculator for orthogonality of a mesh. + class OrthogonalityPropertyCalculator : public PropertyCalculator + { + public: + /// @brief Determine is the calculator can compute the desired results correctly. + /// + /// This has a default of checking that the mesh2d is valid and the location is at edges + virtual bool IsValid(const MeshKernelState& state, const meshkernel::Location location) const override; + + /// @brief Calculate the orthogonality for a mesh + 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 meshkernel::Location location) const override; + }; + + /// @brief Calculator for the edge lengths for a mesh + class EdgeLengthPropertyCalculator : public PropertyCalculator + { + public: + /// @brief Determine is the calculator can compute the desired results correctly. + /// + /// This has a default of checking that the mesh2d is valid and the location is at edges + virtual bool IsValid(const MeshKernelState& state, const meshkernel::Location location) const override; + + /// @brief Calculate the edge-length for a mesh + /// + /// \note This calculator is for mesh edges only + 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 meshkernel::Location location) const override; + }; + + /// @brief Interpolate the depths at the mesh node points. + class InterpolatedSamplePropertyCalculator : public PropertyCalculator + { + public: + /// @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 + bool IsValid(const MeshKernelState& state, const meshkernel::Location location) const override; + + /// @brief Calculate the edge-length for a mesh + /// + /// \note This calculator is for mesh edges only + 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 meshkernel::Location location) const override; + + private: + /// @brief Interpolator for the samples + std::unique_ptr m_sampleInterpolator; + + /// @brief Projection sued for sample data. + meshkernel::Projection m_projection; + + /// @brief Property id. + int m_propertyId = -1; + }; + +} // namespace meshkernelapi diff --git a/libs/MeshKernelApi/src/MeshKernel.cpp b/libs/MeshKernelApi/src/MeshKernel.cpp index af5c94673..f1c98cbc2 100644 --- a/libs/MeshKernelApi/src/MeshKernel.cpp +++ b/libs/MeshKernelApi/src/MeshKernel.cpp @@ -99,12 +99,14 @@ #include "MeshKernelApi/ApiCache/SmallFlowEdgeCentreCache.hpp" #include "MeshKernelApi/MKStateUndoAction.hpp" #include "MeshKernelApi/MeshKernel.hpp" +#include "MeshKernelApi/PropertyCalculator.hpp" #include "MeshKernelApi/State.hpp" #include "MeshKernelApi/Utils.hpp" #include #include +#include #include #include @@ -125,6 +127,31 @@ namespace meshkernelapi /// @brief Stack of undo actions static meshkernel::UndoActionStack meshKernelUndoStack; + int GeneratePropertyId() + { + // The current property id, initialised with a value equal to the last enum in Mesh2D:::Property enum values + static int currentPropertyId = static_cast(meshkernel::Mesh2D::Property::EdgeLength); + + // Increment and return the current property id value. + return ++currentPropertyId; + } + + std::map> allocatePropertyCalculators() + { + std::map> propertyMap; + + int propertyId = static_cast(meshkernel::Mesh2D::Property::Orthogonality); + propertyMap.emplace(propertyId, std::make_unique()); + + propertyId = static_cast(meshkernel::Mesh2D::Property::EdgeLength); + propertyMap.emplace(propertyId, std::make_unique()); + + return propertyMap; + } + + /// @brief Map of property calculators, from an property identifier to the calculator. + static std::map> propertyCalculators = allocatePropertyCalculators(); + static meshkernel::ExitCode HandleException(std::exception_ptr exception_ptr = std::current_exception()) { try @@ -484,6 +511,74 @@ namespace meshkernelapi return lastExitCode; } + MKERNEL_API int mkernel_deallocate_property(int propertyId) + { + lastExitCode = meshkernel::ExitCode::Success; + try + { + + if (!propertyCalculators.contains(propertyId) || propertyCalculators[propertyId] == nullptr) + { + throw meshkernel::MeshKernelError("The property id does not exist: {}.", propertyId); + } + + propertyCalculators.erase(propertyId); + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + + MKERNEL_API int mkernel_mesh2d_is_valid_property(int meshKernelId, const int propertyId, const int locationId, bool& propertyIsAvailable) + { + lastExitCode = meshkernel::ExitCode::Success; + propertyIsAvailable = false; + try + { + const meshkernel::Location location = static_cast(locationId); + propertyIsAvailable = meshKernelState.contains(meshKernelId) && + propertyCalculators.contains(propertyId) && + propertyCalculators[propertyId] != nullptr && + propertyCalculators[propertyId]->IsValid(meshKernelState.at(meshKernelId), location); + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + + MKERNEL_API int mkernel_mesh2d_set_property(int projectionType, int interpolationType, const GeometryList& sampleData, int& propertyId) + { + lastExitCode = meshkernel::ExitCode::Success; + propertyId = -1; + + try + { + meshkernel::range_check::CheckOneOf(projectionType, meshkernel::GetValidProjections(), "Projection"); + auto const projection = static_cast(projectionType); + + meshkernel::range_check::CheckOneOf(interpolationType, {0, 1}, "InterpolationType"); + + int localPropertyId = GeneratePropertyId(); + + if (propertyCalculators.contains(localPropertyId)) + { + throw meshkernel::ConstraintError("The property id already exists: id = {}.", localPropertyId); + } + + propertyCalculators.try_emplace(localPropertyId, std::make_unique(sampleData, projection, interpolationType, localPropertyId)); + propertyId = localPropertyId; + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + MKERNEL_API int mkernel_mesh2d_snap_to_landboundary(int meshKernelId, const GeometryList& selectingPolygon, const GeometryList& landBoundaries) { lastExitCode = meshkernel::ExitCode::Success; @@ -1497,9 +1592,10 @@ namespace meshkernelapi return lastExitCode; } - MKERNEL_API int mkernel_mesh2d_get_property(int meshKernelId, int propertyValue, const GeometryList& geometryList) + MKERNEL_API int mkernel_mesh2d_get_property(int meshKernelId, int propertyValue, int locationId, const GeometryList& geometryList) { lastExitCode = meshkernel::ExitCode::Success; + try { if (!meshKernelState.contains(meshKernelId)) @@ -1507,39 +1603,43 @@ namespace meshkernelapi throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); } - const auto& mesh2d = meshKernelState.at(meshKernelId).m_mesh2d; - - if (!mesh2d || mesh2d->GetNumNodes() <= 0) + if (const auto& mesh2d = meshKernelState.at(meshKernelId).m_mesh2d; mesh2d == nullptr || mesh2d->GetNumNodes() <= 0) { return lastExitCode; } - const auto propertyValueEnum = static_cast(propertyValue); - switch (propertyValueEnum) + if (!propertyCalculators.contains(propertyValue) || propertyCalculators[propertyValue] == nullptr) { - case meshkernel::Mesh2D::Property::Orthogonality: + throw meshkernel::MeshKernelError("The property calculator does not exist."); + } + + const meshkernel::Location location = static_cast(locationId); + + if (geometryList.num_coordinates < propertyCalculators[propertyValue]->Size(meshKernelState.at(meshKernelId), location)) { - std::vector values = mesh2d->GetOrthogonality(); - if (static_cast(geometryList.num_coordinates) < values.size()) - { - throw meshkernel::MeshKernelError("GeometryList with wrong dimensions"); - } - std::copy(values.begin(), values.end(), geometryList.values); + throw meshkernel::ConstraintError("Array size too small to store property values {} < {}.", + geometryList.num_coordinates, + propertyCalculators[propertyValue]->Size(meshKernelState.at(meshKernelId), location)); + } + + if (geometryList.values == nullptr) + { + throw meshkernel::ConstraintError("The property values are null."); } - break; - case meshkernel::Mesh2D::Property::EdgeLength: + + if (propertyCalculators[propertyValue]->IsValid(meshKernelState[meshKernelId], location)) { - mesh2d->ComputeEdgesLengths(); - std::vector values = mesh2d->m_edgeLengths; - if (static_cast(geometryList.num_coordinates) < values.size()) + + if (location == meshkernel::Location::Edges && meshKernelState[meshKernelId].m_mesh2d->m_edgesCenters.empty()) { - throw meshkernel::MeshKernelError("GeometryList with wrong dimensions"); + meshKernelState[meshKernelId].m_mesh2d->ComputeEdgesCenters(); } - std::copy(values.begin(), values.end(), geometryList.values); + + propertyCalculators[propertyValue]->Calculate(meshKernelState[meshKernelId], location, geometryList); } - break; - default: - throw meshkernel::MeshKernelError("Property not supported"); + else + { + throw meshkernel::MeshKernelError("Property not supported at this location"); } } catch (...) @@ -1552,6 +1652,8 @@ namespace meshkernelapi MKERNEL_API int mkernel_mesh2d_get_property_dimension(int meshKernelId, int propertyValue, int& dimension) { lastExitCode = meshkernel::ExitCode::Success; + dimension = -1; + try { if (!meshKernelState.contains(meshKernelId)) @@ -1559,26 +1661,17 @@ namespace meshkernelapi throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); } - const auto& mesh2d = meshKernelState.at(meshKernelId).m_mesh2d; - - if (!mesh2d || mesh2d->GetNumNodes() <= 0) + if (const auto& mesh2d = meshKernelState.at(meshKernelId).m_mesh2d; !mesh2d || mesh2d->GetNumNodes() <= 0) { return lastExitCode; } - const auto propertyValueEnum = static_cast(propertyValue); - dimension = -1; - switch (propertyValueEnum) + if (propertyCalculators.contains(propertyValue) && propertyCalculators[propertyValue] != nullptr) + { + dimension = propertyCalculators[propertyValue]->Size(meshKernelState[meshKernelId], meshkernel::Location::Edges); + } + else { - case meshkernel::Mesh2D::Property::Orthogonality: - dimension = static_cast(mesh2d->GetOrthogonality().size()); - break; - - case meshkernel::Mesh2D::Property::EdgeLength: - mesh2d->ComputeEdgesLengths(); - dimension = static_cast(mesh2d->m_edgeLengths.size()); - break; - default: throw meshkernel::MeshKernelError("Property not supported"); } } @@ -1606,10 +1699,7 @@ namespace meshkernelapi const auto result = meshKernelState[meshKernelId].m_mesh2d->GetSmoothness(); - for (auto i = 0; i < geometryList.num_coordinates; ++i) - { - geometryList.values[i] = result[i]; - } + std::ranges::copy(result, geometryList.values); } catch (...) { diff --git a/libs/MeshKernelApi/src/PropertyCalculator.cpp b/libs/MeshKernelApi/src/PropertyCalculator.cpp new file mode 100644 index 000000000..8992ec80e --- /dev/null +++ b/libs/MeshKernelApi/src/PropertyCalculator.cpp @@ -0,0 +1,112 @@ +#include "MeshKernelApi/PropertyCalculator.hpp" + +#include "MeshKernel/SampleAveragingInterpolator.hpp" +#include "MeshKernel/SampleTriangulationInterpolator.hpp" + +#include +#include + +bool meshkernelapi::OrthogonalityPropertyCalculator::IsValid(const MeshKernelState& state, const meshkernel::Location location) const +{ + return state.m_mesh2d != nullptr && state.m_mesh2d->GetNumNodes() > 0 && location == meshkernel::Location::Edges; +} + +void meshkernelapi::OrthogonalityPropertyCalculator::Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const +{ + + std::vector values = state.m_mesh2d->GetOrthogonality(); + + if (static_cast(geometryList.num_coordinates) < values.size()) + { + throw meshkernel::ConstraintError("GeometryList with wrong dimensions, {} must be greater than or equal to {}", + geometryList.num_coordinates, Size(state, location)); + } + + std::ranges::copy(values, geometryList.values); +} + +int meshkernelapi::OrthogonalityPropertyCalculator::Size(const MeshKernelState& state, const meshkernel::Location location [[maybe_unused]]) const +{ + return static_cast(state.m_mesh2d->GetNumEdges()); +} + +bool meshkernelapi::EdgeLengthPropertyCalculator::IsValid(const MeshKernelState& state, const meshkernel::Location location) const +{ + return state.m_mesh2d != nullptr && state.m_mesh2d->GetNumNodes() > 0 && location == meshkernel::Location::Edges; +} + +void meshkernelapi::EdgeLengthPropertyCalculator::Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const +{ + + state.m_mesh2d->ComputeEdgesLengths(); + std::vector values = state.m_mesh2d->m_edgeLengths; + + if (static_cast(geometryList.num_coordinates) < values.size()) + { + throw meshkernel::ConstraintError("GeometryList with wrong dimensions, {} must be greater than or equal to {}", + geometryList.num_coordinates, Size(state, location)); + } + + std::ranges::copy(values, geometryList.values); +} + +int meshkernelapi::EdgeLengthPropertyCalculator::Size(const MeshKernelState& state, const meshkernel::Location location [[maybe_unused]]) const +{ + return static_cast(state.m_mesh2d->GetNumEdges()); +} + +meshkernelapi::InterpolatedSamplePropertyCalculator::InterpolatedSamplePropertyCalculator(const GeometryList& sampleData, + const meshkernel::Projection projection, + const int interpolationType, + const int propertyId) + : m_projection(projection), + m_propertyId(propertyId) +{ + std::span xNodes(sampleData.coordinates_x, sampleData.num_coordinates); + std::span yNodes(sampleData.coordinates_y, sampleData.num_coordinates); + + if (interpolationType == 0) + { + m_sampleInterpolator = std::make_unique(xNodes, yNodes, m_projection); + } + else if (interpolationType == 1) + { + // Need to pass from api. + meshkernel::InterpolationParameters interpolationParameters{}; + m_sampleInterpolator = std::make_unique(xNodes, yNodes, m_projection, interpolationParameters); + } + + std::span dataSamples(sampleData.values, sampleData.num_coordinates); + m_sampleInterpolator->SetData(m_propertyId, dataSamples); +} + +bool meshkernelapi::InterpolatedSamplePropertyCalculator::IsValid(const MeshKernelState& state, const meshkernel::Location location [[maybe_unused]]) const +{ + return state.m_mesh2d != nullptr && + state.m_mesh2d->GetNumNodes() > 0 && + m_sampleInterpolator->Contains(m_propertyId) && + m_projection == state.m_projection; +} + +void meshkernelapi::InterpolatedSamplePropertyCalculator::Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const +{ + std::span interpolatedSampleData(geometryList.values, geometryList.num_coordinates); + m_sampleInterpolator->Interpolate(m_propertyId, *state.m_mesh2d, location, interpolatedSampleData); +} + +int meshkernelapi::InterpolatedSamplePropertyCalculator::Size(const MeshKernelState& state, const meshkernel::Location location) const +{ + using enum meshkernel::Location; + + switch (location) + { + case Nodes: + return static_cast(state.m_mesh2d->GetNumNodes()); + case Edges: + return static_cast(state.m_mesh2d->GetNumEdges()); + case Faces: + return static_cast(state.m_mesh2d->GetNumFaces()); + default: + return -1; + } +} diff --git a/libs/MeshKernelApi/tests/CMakeLists.txt b/libs/MeshKernelApi/tests/CMakeLists.txt index 27ea964bc..4a7d27a42 100644 --- a/libs/MeshKernelApi/tests/CMakeLists.txt +++ b/libs/MeshKernelApi/tests/CMakeLists.txt @@ -22,6 +22,7 @@ set(SRC_LIST ${SRC_DIR}/Mesh2DCasulliRefinmentTests.cpp ${SRC_DIR}/Mesh2DRefinmentTests.cpp ${SRC_DIR}/Mesh2DTests.cpp + ${SRC_DIR}/MeshPropertyTests.cpp ${SRC_DIR}/PolygonTests.cpp ${SRC_DIR}/SmoothnessOrthogonalisationTests.cpp ${SRC_DIR}/TestMeshGeneration.cpp diff --git a/libs/MeshKernelApi/tests/src/Mesh2DTests.cpp b/libs/MeshKernelApi/tests/src/Mesh2DTests.cpp index 8eec93f5b..163cff06c 100644 --- a/libs/MeshKernelApi/tests/src/Mesh2DTests.cpp +++ b/libs/MeshKernelApi/tests/src/Mesh2DTests.cpp @@ -174,12 +174,13 @@ TEST(Mesh2DTests, Mesh2DGetPropertyTest) ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); // Execute + int locationId = static_cast(meshkernel::Location::Edges); meshkernelapi::GeometryList propertyvalues{}; propertyvalues.num_coordinates = geometryListDimension; propertyvalues.geometry_separator = meshkernel::constants::missing::doubleValue; std::vector values(geometryListDimension); propertyvalues.values = values.data(); - errorCode = mkernel_mesh2d_get_property(meshKernelId, 0, propertyvalues); + errorCode = mkernel_mesh2d_get_property(meshKernelId, 0, locationId, propertyvalues); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); // Assert diff --git a/libs/MeshKernelApi/tests/src/MeshPropertyTests.cpp b/libs/MeshKernelApi/tests/src/MeshPropertyTests.cpp new file mode 100644 index 000000000..16c33a28f --- /dev/null +++ b/libs/MeshKernelApi/tests/src/MeshPropertyTests.cpp @@ -0,0 +1,308 @@ +#include +#include +#include +#include + +#include "CartesianApiTestFixture.hpp" +#include "MeshKernel/Parameters.hpp" +#include "MeshKernelApi/BoundingBox.hpp" +#include "MeshKernelApi/Mesh2D.hpp" +#include "MeshKernelApi/MeshKernel.hpp" +#include "TestUtils/MakeMeshes.hpp" + +// namespace aliases +namespace mk = meshkernel; +namespace mkapi = meshkernelapi; + +TEST(MeshPropertyTests, BathymetryTest) +{ + const int projectionType = 0; + + const int numXCoords = 4; + const int numYCoords = 4; + const int numberOfEdges = (numXCoords - 1) * numYCoords + numXCoords * (numYCoords - 1); + + const int numXBathyCoords = 5; + const int numYBathyCoords = 5; + const int numberOfBathyCoordinates = numXBathyCoords * numYBathyCoords; + + const double origin = 0.0; + const double delta = 1.0; + + int bathymetryPropertyId = -1; + + int meshKernelId = meshkernel::constants::missing::intValue; + int errorCode; + + // Clear the meshkernel state and undo stack before starting the test. + errorCode = mkapi::mkernel_clear_state(); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + + errorCode = mkapi::mkernel_allocate_state(0, meshKernelId); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + + meshkernel::MakeGridParameters makeGridParameters; + + // num_columns and num_rows indicate number of elements in each direction + makeGridParameters.num_columns = numXCoords - 1; + makeGridParameters.num_rows = numYCoords - 1; + makeGridParameters.angle = 0.0; + makeGridParameters.origin_x = origin; + makeGridParameters.origin_y = origin; + makeGridParameters.block_size_x = delta; + makeGridParameters.block_size_y = delta; + + // Generate curvilinear grid. + errorCode = mkapi::mkernel_curvilinear_compute_rectangular_grid(meshKernelId, makeGridParameters); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + + // Convert the curvilinear grid to an unstructured grid. + errorCode = mkapi::mkernel_curvilinear_convert_to_mesh2d(meshKernelId); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + + //-------------------------------- + + mkapi::Mesh2D meshData{}; + mkapi::mkernel_mesh2d_get_dimensions(meshKernelId, meshData); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + + std::vector xs(meshData.num_nodes); + std::vector ys(meshData.num_nodes); + std::vector edges(meshData.num_edges * 2); + + meshData.node_x = xs.data(); + meshData.node_y = ys.data(); + meshData.edge_nodes = edges.data(); + + mkapi::mkernel_mesh2d_get_node_edge_data(meshKernelId, meshData); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + + //-------------------------------- + + std::vector bathymetryData(numberOfBathyCoordinates, 1.0); + std::vector bathymetryXNodes(numberOfBathyCoordinates); + std::vector bathymetryYNodes(numberOfBathyCoordinates); + + double bathyY = origin - 0.5 * delta; + size_t count = 0; + + for (size_t i = 0; i < numXBathyCoords; ++i) + { + double bathyX = origin - 0.5 * delta; + + for (size_t j = 0; j < numYBathyCoords; ++j) + { + bathymetryXNodes[count] = bathyX; + bathymetryYNodes[count] = bathyY; + bathymetryData[count] = bathyX; + + bathyX += delta; + ++count; + } + + bathyY += delta; + } + + mkapi::GeometryList sampleData{}; + + sampleData.num_coordinates = numberOfBathyCoordinates; + sampleData.values = bathymetryData.data(); + sampleData.coordinates_x = bathymetryXNodes.data(); + sampleData.coordinates_y = bathymetryYNodes.data(); + + 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 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.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, numberOfEdges); + + auto locationId = static_cast(meshkernel::Location::Edges); + + mkapi::GeometryList propertyData{}; + std::vector retrievedPropertyData(numberOfEdges, -1.0); + propertyData.num_coordinates = numberOfEdges; + propertyData.values = retrievedPropertyData.data(); + + errorCode = mkapi::mkernel_mesh2d_get_property(meshKernelId, bathymetryPropertyId, locationId, propertyData); + + ASSERT_EQ(mk::ExitCode::Success, errorCode); + + for (size_t i = 0; i < retrievedPropertyData.size(); ++i) + { + EXPECT_NEAR(retrievedPropertyData[i], expectedInterpolatedData[i], tolerance); + } + + //-------------------------------- + // Remove kernel and property id from state. + + errorCode = mkapi::mkernel_expunge_state(meshKernelId); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + + errorCode = mkapi::mkernel_deallocate_property(bathymetryPropertyId); + ASSERT_EQ(mk::ExitCode::Success, errorCode); +} + +TEST(MeshPropertyTests, PropertyFailureTest) +{ + const int projectionType = 0; + + const int numXCoords = 4; + const int numYCoords = 4; + const int numberOfCoordinates = numXCoords * numYCoords; + + const int numXBathyCoords = 5; + const int numYBathyCoords = 5; + const int numberOfBathyCoordinates = numXBathyCoords * numYBathyCoords; + + const double origin = 0.0; + const double delta = 1.0; + + int bathymetryPropertyId = -1; + + int meshKernelId = meshkernel::constants::missing::intValue; + int errorCode; + + // Clear the meshkernel state and undo stack before starting the test. + errorCode = mkapi::mkernel_clear_state(); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + + errorCode = mkapi::mkernel_allocate_state(0, meshKernelId); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + + meshkernel::MakeGridParameters makeGridParameters; + + // num_columns and num_rows indicate number of elements in each direction + makeGridParameters.num_columns = numXCoords - 1; + makeGridParameters.num_rows = numYCoords - 1; + makeGridParameters.angle = 0.0; + makeGridParameters.origin_x = origin; + makeGridParameters.origin_y = origin; + makeGridParameters.block_size_x = delta; + makeGridParameters.block_size_y = delta; + + // Generate curvilinear grid. + errorCode = mkapi::mkernel_curvilinear_compute_rectangular_grid(meshKernelId, makeGridParameters); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + + // Convert the curvilinear grid to an unstructured grid. + errorCode = mkapi::mkernel_curvilinear_convert_to_mesh2d(meshKernelId); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + + //-------------------------------- + + mkapi::Mesh2D meshData{}; + mkapi::mkernel_mesh2d_get_dimensions(meshKernelId, meshData); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + + std::vector xs(meshData.num_nodes); + std::vector ys(meshData.num_nodes); + std::vector edges(meshData.num_edges * 2); + + meshData.node_x = xs.data(); + meshData.node_y = ys.data(); + meshData.edge_nodes = edges.data(); + + mkapi::mkernel_mesh2d_get_node_edge_data(meshKernelId, meshData); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + + //-------------------------------- + + std::vector bathymetryData(numberOfBathyCoordinates, 1.0); + std::vector bathymetryXNodes(numberOfBathyCoordinates); + std::vector bathymetryYNodes(numberOfBathyCoordinates); + + double bathyY = origin - 0.5 * delta; + size_t count = 0; + + for (size_t i = 0; i < numXBathyCoords; ++i) + { + double bathyX = origin - 0.5 * delta; + + for (size_t j = 0; j < numYBathyCoords; ++j) + { + bathymetryXNodes[count] = bathyX; + bathymetryYNodes[count] = bathyY; + bathymetryData[count] = bathyX; + + bathyX += delta; + ++count; + } + + bathyY += delta; + } + + mkapi::GeometryList sampleData{}; + + sampleData.num_coordinates = numberOfBathyCoordinates; + sampleData.values = bathymetryData.data(); + sampleData.coordinates_x = bathymetryXNodes.data(); + sampleData.coordinates_y = bathymetryYNodes.data(); + + auto locationId = static_cast(meshkernel::Location::Edges); + bool hasBathymetryData = false; + + errorCode = mkapi::mkernel_mesh2d_is_valid_property(meshKernelId, bathymetryPropertyId, locationId, hasBathymetryData); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + EXPECT_FALSE(hasBathymetryData); + + 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, locationId, hasBathymetryData); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + EXPECT_TRUE(hasBathymetryData); + + //-------------------------------- + // Start test + + int sampleDataSize = -1; + + // Expected test failure due to invalid meshkernel id + errorCode = mkapi::mkernel_mesh2d_get_property_dimension(meshKernelId + 100, bathymetryPropertyId, sampleDataSize); + ASSERT_EQ(mk::ExitCode::MeshKernelErrorCode, errorCode); + + mkapi::GeometryList propertyData{}; + std::vector retrievedPropertyData(numberOfCoordinates, -1.0); + propertyData.num_coordinates = numberOfCoordinates; + + // Expected test failure due to values set to null + propertyData.values = nullptr; + errorCode = mkapi::mkernel_mesh2d_get_property(meshKernelId, bathymetryPropertyId, locationId, propertyData); + ASSERT_EQ(mk::ExitCode::ConstraintErrorCode, errorCode); + + propertyData.values = retrievedPropertyData.data(); + + // Expected test failure due to invalid property id + errorCode = mkapi::mkernel_mesh2d_get_property(meshKernelId, bathymetryPropertyId + 100, locationId, propertyData); + ASSERT_EQ(mk::ExitCode::MeshKernelErrorCode, errorCode); + + // Expected test failure due to incorrect size + propertyData.num_coordinates = numberOfCoordinates - 1; + errorCode = mkapi::mkernel_mesh2d_get_property(meshKernelId, bathymetryPropertyId, locationId, propertyData); + ASSERT_EQ(mk::ExitCode::ConstraintErrorCode, errorCode); + + //-------------------------------- + // Remove kernel and property id from state. + + errorCode = mkapi::mkernel_deallocate_property(bathymetryPropertyId); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + + errorCode = mkapi::mkernel_mesh2d_is_valid_property(meshKernelId, bathymetryPropertyId, locationId, hasBathymetryData); + ASSERT_EQ(mk::ExitCode::Success, errorCode); + EXPECT_FALSE(hasBathymetryData); + + errorCode = mkapi::mkernel_expunge_state(meshKernelId); + ASSERT_EQ(mk::ExitCode::Success, errorCode); +}