From f489c3fc2604933fc8429162795bdf5ae19c120e Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 9 Jan 2025 10:31:53 +0100 Subject: [PATCH] GRIDEDIT-1552 Moved orthogonalisation and smooting calculation to separate calculators --- libs/MeshKernel/CMakeLists.txt | 4 + libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 8 -- .../include/MeshKernel/MeshOrthogonality.hpp | 48 +++++++++++ .../include/MeshKernel/MeshSmoothness.hpp | 51 ++++++++++++ libs/MeshKernel/src/Mesh2D.cpp | 71 +--------------- libs/MeshKernel/src/MeshOrthogonality.cpp | 78 ++++++++++++++++++ libs/MeshKernel/src/MeshSmoothness.cpp | 81 +++++++++++++++++++ libs/MeshKernel/tests/src/Mesh2DTest.cpp | 9 ++- libs/MeshKernelApi/src/MeshKernel.cpp | 31 +++++-- libs/MeshKernelApi/src/PropertyCalculator.cpp | 21 +++-- 10 files changed, 311 insertions(+), 91 deletions(-) create mode 100644 libs/MeshKernel/include/MeshKernel/MeshOrthogonality.hpp create mode 100644 libs/MeshKernel/include/MeshKernel/MeshSmoothness.hpp create mode 100644 libs/MeshKernel/src/MeshOrthogonality.cpp create mode 100644 libs/MeshKernel/src/MeshSmoothness.cpp diff --git a/libs/MeshKernel/CMakeLists.txt b/libs/MeshKernel/CMakeLists.txt index 8de376a45..8c663d46c 100644 --- a/libs/MeshKernel/CMakeLists.txt +++ b/libs/MeshKernel/CMakeLists.txt @@ -49,6 +49,8 @@ set( ${SRC_DIR}/Mesh2DGenerateGlobal.cpp ${SRC_DIR}/Mesh2DIntersections.cpp ${SRC_DIR}/Mesh2DToCurvilinear.cpp + ${SRC_DIR}/MeshOrthogonality.cpp + ${SRC_DIR}/MeshSmoothness.cpp ${SRC_DIR}/MeshTriangulation.cpp ${SRC_DIR}/MeshRefinement.cpp ${SRC_DIR}/Network1D.cpp @@ -172,6 +174,8 @@ set( ${DOMAIN_INC_DIR}/Mesh2DToCurvilinear.hpp ${DOMAIN_INC_DIR}/MeshConversion.hpp ${DOMAIN_INC_DIR}/MeshInterpolation.hpp + ${DOMAIN_INC_DIR}/MeshOrthogonality.hpp + ${DOMAIN_INC_DIR}/MeshSmoothness.hpp ${DOMAIN_INC_DIR}/MeshTriangulation.hpp ${DOMAIN_INC_DIR}/MeshRefinement.hpp ${DOMAIN_INC_DIR}/MeshTransformation.hpp diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 2ae64fbc5..aabb5d722 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -229,14 +229,6 @@ namespace meshkernel /// @brief Computes m_nodesNodes, see class members void ComputeNodeNeighbours(); - /// @brief Get the orthogonality values, the inner product of edges and segments connecting the face circumcenters - /// @return The edge orthogonality - [[nodiscard]] std::vector GetOrthogonality() const; - - /// @brief Gets the smoothness values, ratios of the face areas - /// @return The smoothness at the edges - [[nodiscard]] std::vector GetSmoothness() const; - /// @brief Gets the aspect ratios (the ratios edges lengths to flow edges lengths) /// @param[in,out] aspectRatios The aspect ratios (passed as reference to avoid re-allocation) void ComputeAspectRatios(std::vector& aspectRatios); diff --git a/libs/MeshKernel/include/MeshKernel/MeshOrthogonality.hpp b/libs/MeshKernel/include/MeshKernel/MeshOrthogonality.hpp new file mode 100644 index 000000000..a6ff238dc --- /dev/null +++ b/libs/MeshKernel/include/MeshKernel/MeshOrthogonality.hpp @@ -0,0 +1,48 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2025. +// +// 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 "MeshKernel/Mesh2D.hpp" + +namespace meshkernel +{ + /// @brief Compute the orthogonality value for the edges + class MeshOrthogonality + { + public: + /// @brief Compute the orthogonality values returning values in a vector + std::vector Compute(const Mesh2D& mesh) const; + + /// @brief Compute the orthogonality values overwriting the values in an array + void Compute(const Mesh2D& mesh, std::span orthogonality) const; + }; + +} // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/MeshSmoothness.hpp b/libs/MeshKernel/include/MeshKernel/MeshSmoothness.hpp new file mode 100644 index 000000000..0e1d65cc7 --- /dev/null +++ b/libs/MeshKernel/include/MeshKernel/MeshSmoothness.hpp @@ -0,0 +1,51 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2025. +// +// 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 "MeshKernel/Mesh2D.hpp" + +namespace meshkernel +{ + /// @brief Compute the smoothness value for the edges + class MeshSmoothness + { + public: + /// @brief Compute the smoothness values returning values in a vector + std::vector Compute(const Mesh2D& mesh) const; + + /// @brief Compute the smoothness values overwriting the values in an array + void Compute(const Mesh2D& mesh, std::span smoothness) const; + + private: + static constexpr double m_minimumCellArea = 1e-12; ///< Minimum cell area + }; + +} // namespace meshkernel diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index d9fe7b0e8..90dc2111d 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -31,6 +31,7 @@ #include "MeshKernel/Entities.hpp" #include "MeshKernel/Exceptions.hpp" #include "MeshKernel/Mesh2DIntersections.hpp" +#include "MeshKernel/MeshOrthogonality.hpp" #include "MeshKernel/Operations.hpp" #include "MeshKernel/Polygon.hpp" #include "MeshKernel/Polygons.hpp" @@ -1255,73 +1256,6 @@ void Mesh2D::ComputeNodeNeighbours() } } -std::vector Mesh2D::GetOrthogonality() const -{ - std::vector result(GetNumEdges()); - const auto numEdges = GetNumEdges(); - for (UInt e = 0; e < numEdges; e++) - { - auto val = constants::missing::doubleValue; - const auto firstNode = m_edges[e].first; - const auto secondNode = m_edges[e].second; - const auto firstFaceIndex = m_edgesFaces[e][0]; - const auto secondFaceIndex = m_edgesFaces[e][1]; - - if (firstNode != constants::missing::uintValue && - secondNode != constants::missing::uintValue && - firstFaceIndex != constants::missing::uintValue && - secondFaceIndex != constants::missing::uintValue && !IsEdgeOnBoundary(e)) - { - val = NormalizedInnerProductTwoSegments(m_nodes[firstNode], - m_nodes[secondNode], - m_facesCircumcenters[firstFaceIndex], - m_facesCircumcenters[secondFaceIndex], - m_projection); - - if (val != constants::missing::doubleValue) - { - val = std::abs(val); - } - } - result[e] = val; - } - return result; -} - -std::vector Mesh2D::GetSmoothness() const -{ - std::vector result(GetNumEdges()); - const auto numEdges = GetNumEdges(); - for (UInt e = 0; e < numEdges; e++) - { - auto val = constants::missing::doubleValue; - const auto firstNode = m_edges[e].first; - const auto secondNode = m_edges[e].second; - const auto firstFaceIndex = m_edgesFaces[e][0]; - const auto secondFaceIndex = m_edgesFaces[e][1]; - - if (firstNode != constants::missing::uintValue && - secondNode != constants::missing::uintValue && - firstFaceIndex != constants::missing::uintValue && - secondFaceIndex != constants::missing::uintValue && !IsEdgeOnBoundary(e)) - { - const auto leftFaceArea = m_faceArea[firstFaceIndex]; - const auto rightFaceArea = m_faceArea[secondFaceIndex]; - - if (leftFaceArea > m_minimumCellArea && rightFaceArea > m_minimumCellArea) - { - val = rightFaceArea / leftFaceArea; - if (val < 1.0) - { - val = 1.0 / val; - } - } - } - result[e] = val; - } - return result; -} - void Mesh2D::ComputeAverageFlowEdgesLength(std::vector& edgesLength, std::vector& averageFlowEdgesLength) const { @@ -1789,7 +1723,8 @@ std::vector Mesh2D::FilterBasedOnMetric(Location location, std::vector result(numFaces, false); // Retrieve orthogonality values - const std::vector metricValues = GetOrthogonality(); + MeshOrthogonality meshOrthogonality; + const std::vector metricValues(meshOrthogonality.Compute(*this)); // Loop through faces and compute how many edges have the metric within the range for (UInt f = 0; f < numFaces; ++f) diff --git a/libs/MeshKernel/src/MeshOrthogonality.cpp b/libs/MeshKernel/src/MeshOrthogonality.cpp new file mode 100644 index 000000000..45c8b69b7 --- /dev/null +++ b/libs/MeshKernel/src/MeshOrthogonality.cpp @@ -0,0 +1,78 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2025. +// +// 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/MeshOrthogonality.hpp" +#include "MeshKernel/Constants.hpp" +#include "MeshKernel/Exceptions.hpp" +#include "MeshKernel/Operations.hpp" + +std::vector meshkernel::MeshOrthogonality::Compute(const Mesh2D& mesh) const +{ + std::vector orthogonality(mesh.GetNumEdges(), constants::missing::doubleValue); + Compute(mesh, orthogonality); + + return orthogonality; +} + +void meshkernel::MeshOrthogonality::Compute(const Mesh2D& mesh, std::span orthogonality) const +{ + if (orthogonality.size() != mesh.GetNumEdges()) + { + throw ConstraintError("array for orthogonality values is not the correct size"); + } + + const auto numEdges = mesh.GetNumEdges(); + + for (UInt e = 0; e < numEdges; e++) + { + auto val = constants::missing::doubleValue; + + const auto [firstNode, secondNode] = mesh.GetEdge(e); + + const auto firstFaceIndex = mesh.m_edgesFaces[e][0]; + const auto secondFaceIndex = mesh.m_edgesFaces[e][1]; + + if (firstNode != constants::missing::uintValue && + secondNode != constants::missing::uintValue && + firstFaceIndex != constants::missing::uintValue && + secondFaceIndex != constants::missing::uintValue && !mesh.IsEdgeOnBoundary(e)) + { + val = NormalizedInnerProductTwoSegments(mesh.Node(firstNode), + mesh.Node(secondNode), + mesh.m_facesCircumcenters[firstFaceIndex], + mesh.m_facesCircumcenters[secondFaceIndex], + mesh.m_projection); + + if (val != constants::missing::doubleValue) + { + val = std::abs(val); + } + } + + orthogonality[e] = val; + } +} diff --git a/libs/MeshKernel/src/MeshSmoothness.cpp b/libs/MeshKernel/src/MeshSmoothness.cpp new file mode 100644 index 000000000..02c1d35c7 --- /dev/null +++ b/libs/MeshKernel/src/MeshSmoothness.cpp @@ -0,0 +1,81 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2025. +// +// 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/MeshSmoothness.hpp" +#include "MeshKernel/Constants.hpp" +#include "MeshKernel/Exceptions.hpp" + +std::vector meshkernel::MeshSmoothness::Compute(const Mesh2D& mesh) const +{ + + std::vector smoothness(mesh.GetNumEdges(), constants::missing::doubleValue); + Compute(mesh, smoothness); + + return smoothness; +} + +void meshkernel::MeshSmoothness::Compute(const Mesh2D& mesh, std::span smoothness) const +{ + + if (smoothness.size() != mesh.GetNumEdges()) + { + throw ConstraintError("array for smoothness values is not the correct size"); + } + + const UInt numEdges = mesh.GetNumEdges(); + + for (UInt e = 0; e < numEdges; e++) + { + auto val = constants::missing::doubleValue; + + const auto [firstNode, secondNode] = mesh.GetEdge(e); + + const auto firstFaceIndex = mesh.m_edgesFaces[e][0]; + const auto secondFaceIndex = mesh.m_edgesFaces[e][1]; + + if (firstNode != constants::missing::uintValue && + secondNode != constants::missing::uintValue && + firstFaceIndex != constants::missing::uintValue && + secondFaceIndex != constants::missing::uintValue && !mesh.IsEdgeOnBoundary(e)) + { + const auto leftFaceArea = mesh.m_faceArea[firstFaceIndex]; + const auto rightFaceArea = mesh.m_faceArea[secondFaceIndex]; + + if (leftFaceArea > m_minimumCellArea && rightFaceArea > m_minimumCellArea) + { + val = rightFaceArea / leftFaceArea; + + if (val < 1.0) + { + val = 1.0 / val; + } + } + } + + smoothness[e] = val; + } +} diff --git a/libs/MeshKernel/tests/src/Mesh2DTest.cpp b/libs/MeshKernel/tests/src/Mesh2DTest.cpp index 9c40403b0..2542b0d1a 100644 --- a/libs/MeshKernel/tests/src/Mesh2DTest.cpp +++ b/libs/MeshKernel/tests/src/Mesh2DTest.cpp @@ -9,6 +9,8 @@ #include "MeshKernel/Mesh2D.hpp" #include "MeshKernel/Mesh2DIntersections.hpp" #include "MeshKernel/Mesh2DToCurvilinear.hpp" +#include "MeshKernel/MeshOrthogonality.hpp" +#include "MeshKernel/MeshSmoothness.hpp" #include "MeshKernel/Operations.hpp" #include "MeshKernel/Polygons.hpp" #include "MeshKernel/RemoveDisconnectedRegions.hpp" @@ -1417,7 +1419,8 @@ TEST(Mesh2D, GetSmoothness_OnTriangularMesh_ShouldgetSmoothnessValues) const auto mesh = ReadLegacyMesh2DFromFile(TEST_FOLDER + "/data/TestOrthogonalizationMediumTriangularGrid_net.nc"); // Execute - const auto smoothness = mesh->GetSmoothness(); + meshkernel::MeshSmoothness meshSmoothness; + const auto smoothness = meshSmoothness.Compute(*mesh); //->GetSmoothness(); // Assert const double tolerance = 1e-6; @@ -1433,7 +1436,9 @@ TEST(Mesh2D, GetOrthogonality_OnTriangularMesh_ShouldGetOrthogonalityValues) const auto mesh = ReadLegacyMesh2DFromFile(TEST_FOLDER + "/data/TestOrthogonalizationMediumTriangularGrid_net.nc"); // Execute - const auto orthogonality = mesh->GetOrthogonality(); + meshkernel::MeshOrthogonality meshOrthogonality; + const auto orthogonality = meshOrthogonality.Compute(*mesh); + // const auto orthogonality = mesh->GetOrthogonality(); // Assert const double tolerance = 1e-6; diff --git a/libs/MeshKernelApi/src/MeshKernel.cpp b/libs/MeshKernelApi/src/MeshKernel.cpp index 40c3b0a44..16fe15a87 100644 --- a/libs/MeshKernelApi/src/MeshKernel.cpp +++ b/libs/MeshKernelApi/src/MeshKernel.cpp @@ -69,7 +69,9 @@ #include #include #include +#include #include +#include #include #include #include @@ -1573,17 +1575,26 @@ namespace meshkernelapi return lastExitCode; } - const auto result = meshKernelState[meshKernelId].m_mesh2d->GetOrthogonality(); + // const auto result = meshKernelState[meshKernelId].m_mesh2d->GetOrthogonality(); - if (static_cast(geometryList.num_coordinates) != result.size()) + // if (static_cast(geometryList.num_coordinates) != result.size()) + // { + // throw meshkernel::MeshKernelError("The value array has not the same size of the result array storing the orthogonality values at the edges"); + // } + + // for (auto i = 0; i < geometryList.num_coordinates; ++i) + // { + // geometryList.values[i] = result[i]; + // } + + if (static_cast(geometryList.num_coordinates) != meshKernelState[meshKernelId].m_mesh2d->GetNumEdges()) { throw meshkernel::MeshKernelError("The value array has not the same size of the result array storing the orthogonality values at the edges"); } - for (auto i = 0; i < geometryList.num_coordinates; ++i) - { - geometryList.values[i] = result[i]; - } + std::span orthogonality(geometryList.values, geometryList.num_coordinates); + meshkernel::MeshOrthogonality meshOrthogonality; + meshOrthogonality.Compute(*meshKernelState[meshKernelId].m_mesh2d, orthogonality); } catch (...) { @@ -1697,9 +1708,13 @@ namespace meshkernelapi return lastExitCode; } - const auto result = meshKernelState[meshKernelId].m_mesh2d->GetSmoothness(); + // const auto result = meshKernelState[meshKernelId].m_mesh2d->GetSmoothness(); + + // std::ranges::copy(result, geometryList.values); - std::ranges::copy(result, geometryList.values); + std::span smoothness(geometryList.values, geometryList.num_coordinates); + meshkernel::MeshSmoothness meshSmoothness; + meshSmoothness.Compute(*meshKernelState[meshKernelId].m_mesh2d, smoothness); } catch (...) { diff --git a/libs/MeshKernelApi/src/PropertyCalculator.cpp b/libs/MeshKernelApi/src/PropertyCalculator.cpp index 8992ec80e..33b8f3dfc 100644 --- a/libs/MeshKernelApi/src/PropertyCalculator.cpp +++ b/libs/MeshKernelApi/src/PropertyCalculator.cpp @@ -1,10 +1,12 @@ #include "MeshKernelApi/PropertyCalculator.hpp" +#include "MeshKernel/MeshOrthogonality.hpp" #include "MeshKernel/SampleAveragingInterpolator.hpp" #include "MeshKernel/SampleTriangulationInterpolator.hpp" #include #include +#include bool meshkernelapi::OrthogonalityPropertyCalculator::IsValid(const MeshKernelState& state, const meshkernel::Location location) const { @@ -13,16 +15,25 @@ bool meshkernelapi::OrthogonalityPropertyCalculator::IsValid(const MeshKernelSta 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()) + if (geometryList.num_coordinates < static_cast(state.m_mesh2d->GetNumEdges())) { 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); + std::span orthogonality(geometryList.values, geometryList.num_coordinates); + meshkernel::MeshOrthogonality meshOrthogonality; + meshOrthogonality.Compute(*state.m_mesh2d, orthogonality); + + // 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