From b914cbc9e0c7f48c3e0003dd7b891f28368b60bd Mon Sep 17 00:00:00 2001 From: BillSenior Date: Fri, 24 Nov 2023 18:15:31 +0100 Subject: [PATCH 01/15] GRIDEDIT-785 Small refactoring --- .../include/MeshKernel/MeshRefinement.hpp | 66 ++++++++- libs/MeshKernel/src/MeshRefinement.cpp | 139 ++++++++++++------ 2 files changed, 160 insertions(+), 45 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp index 38604049a..916aaee0f 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp @@ -37,6 +37,59 @@ namespace meshkernel // Forward declarations class Mesh2D; +#if 0 + + class ComputeRefinement + { + public: + virtual ~ComputeRefinement() = default; + + virtual void compute() const = 0; + + private: + }; + + class EdgeSizeBasedRefinement : public ComputeRefinement + { + public: + void compute() const override; + + private: + }; + + class SamplesBasedRefinement : public ComputeRefinement + { + public: + // Will loop over all elememts in the mesh calling ComputeForFace + void compute() const override; + + private: + virtual void ComputeForFace() const = 0; + }; + + class WaveCourantRefinement : public SamplesBasedRefinement + { + public: + private: + void ComputeForFace(const UInt face) const override; + }; + + class RefinementLevelsRefinement : public SamplesBasedRefinement + { + public: + private: + void ComputeForFace(const UInt face) const override; + }; + + class RidgeDetectionRefinement : public SamplesBasedRefinement + { + public: + private: + void ComputeForFace(const UInt face) const override; + }; + +#endif + /// @brief A class used to refine a Mesh2D instance /// /// Mesh refinement operates on Mesh2D and is based on @@ -89,6 +142,11 @@ namespace meshkernel RidgeDetection = 3 }; + /// @brief Convert the integer value of refinementInt to a valid value of RefinementType. + /// + /// If no such RefinementType value exists then a ConstraintError will be thrown. + static RefinementType GetRefinementTypeValue(const int refinementInt); + public: /// @brief The constructor for refining based on samples /// @param[in] mesh The mesh to be refined @@ -219,13 +277,17 @@ namespace meshkernel /// @brief Compute which edge will be below the minimum size after refinement void ComputeEdgeBelowMinSizeAfterRefinement(); + /// @brief Review the refinement, some elements marked for refinement may not need to be refined. + void ReviewRefinement(const int level); + RTree m_samplesRTree; ///< The sample node RTree std::vector m_faceMask; ///< Compute face without hanging nodes (1), refine face with hanging nodes (2), do not refine cell at all (0) or refine face outside polygon (-2) std::vector m_edgeMask; ///< If 0, edge is not split std::vector m_isEdgeBelowMinSizeAfterRefinement; ///< If an edge is below the minimum size after refinement - std::vector m_nodeMask; ///< The node mask used in the refinement process - std::vector m_brotherEdges; ///< The index of the brother edge for each edge + // TODO add enum for nodeMask values: Or comment describing the meaning of the values -2, -1, 0 and 1. + std::vector m_nodeMask; ///< The node mask used in the refinement process + std::vector m_brotherEdges; ///< The index of the brother edge for each edge /// Local caches std::vector m_isHangingNodeCache; ///< Cache for maintaining if node is hanging diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index f35744696..e9a9fff87 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -35,6 +35,25 @@ using meshkernel::Mesh2D; using meshkernel::MeshRefinement; +MeshRefinement::RefinementType MeshRefinement::GetRefinementTypeValue(const int refinementInt) +{ + static constexpr int WaveCourantInt = static_cast(RefinementType::WaveCourant); + static constexpr int RefinementLevelsInt = static_cast(RefinementType::RefinementLevels); + static constexpr int RidgeDetectionInt = static_cast(RefinementType::RidgeDetection); + + switch (refinementInt) + { + case WaveCourantInt: + return RefinementType::WaveCourant; + case RefinementLevelsInt: + return RefinementType::RefinementLevels; + case RidgeDetectionInt: + return RefinementType::RidgeDetection; + default: + throw ConstraintError("Cannot convert the integer value, {}, for refinement to a valid refinement type enumeration", refinementInt); + } +} + MeshRefinement::MeshRefinement(std::shared_ptr mesh, std::shared_ptr interpolant, const MeshRefinementParameters& meshRefinementParameters) @@ -43,7 +62,7 @@ MeshRefinement::MeshRefinement(std::shared_ptr mesh, { CheckMeshRefinementParameters(meshRefinementParameters); m_meshRefinementParameters = meshRefinementParameters; - m_refinementType = static_cast(m_meshRefinementParameters.refinement_type); + m_refinementType = GetRefinementTypeValue(m_meshRefinementParameters.refinement_type); } MeshRefinement::MeshRefinement(std::shared_ptr mesh, @@ -56,7 +75,7 @@ MeshRefinement::MeshRefinement(std::shared_ptr mesh, { CheckMeshRefinementParameters(meshRefinementParameters); m_meshRefinementParameters = meshRefinementParameters; - m_refinementType = static_cast(m_meshRefinementParameters.refinement_type); + m_refinementType = GetRefinementTypeValue(m_meshRefinementParameters.refinement_type); } MeshRefinement::MeshRefinement(std::shared_ptr mesh, @@ -69,6 +88,65 @@ MeshRefinement::MeshRefinement(std::shared_ptr mesh, m_meshRefinementParameters = meshRefinementParameters; } +void MeshRefinement::ReviewRefinement(const int level) +{ + + // TODO is this comment correct? Looks like its disabling refinement + // if one face node is in polygon enable face refinement + // TODO add enum for nodeMask values + for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) + { + for (UInt n = 0; n < m_mesh->GetNumFaceEdges(f); ++n) + { + const auto nodeIndex = m_mesh->m_facesNodes[f][n]; + + bool disableRefinement = (level == 0 ? m_nodeMask[nodeIndex] != 0 && m_nodeMask[nodeIndex] != -2 : m_nodeMask[nodeIndex] != 1); + + if (disableRefinement) + { + m_faceMask[f] = 0; + break; + } + } + } + +#if 0 + if (level == 0) + { + // TODO is this comment correct? Looks like its disabling refinement + // if one face node is in polygon enable face refinement + for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) + { + for (UInt n = 0; n < m_mesh->GetNumFaceEdges(f); ++n) + { + const auto nodeIndex = m_mesh->m_facesNodes[f][n]; + if (m_nodeMask[nodeIndex] != 0 && m_nodeMask[nodeIndex] != -2) + { + m_faceMask[f] = 0; + break; + } + } + } + } + if (level > 0) + { + // if one face node is not in polygon disable refinement + for (UInt f = 0; f < m_mesh->GetNumFaces(); f++) + { + for (UInt n = 0; n < m_mesh->GetNumFaceEdges(f); n++) + { + const auto nodeIndex = m_mesh->m_facesNodes[f][n]; + if (m_nodeMask[nodeIndex] != 1) + { + m_faceMask[f] = 0; + break; + } + } + } + } +#endif +} + void MeshRefinement::Compute() { // administrate mesh once more @@ -108,7 +186,7 @@ void MeshRefinement::Compute() auto numFacesAfterRefinement = m_mesh->GetNumFaces(); // reserve some extra capacity for the node mask m_nodeMask.reserve(m_nodeMask.size() * 2); - for (auto level = 0; level < m_meshRefinementParameters.max_num_refinement_iterations; level++) + for (int level = 0; level < m_meshRefinementParameters.max_num_refinement_iterations; level++) { // Compute all edge lengths at once ComputeEdgeBelowMinSizeAfterRefinement(); @@ -123,48 +201,11 @@ void MeshRefinement::Compute() ComputeRefinementMaskFromEdgeSize(); } - if (level == 0) - { - // if one face node is in polygon enable face refinement - for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) - { - bool activeNodeFound = false; - for (UInt n = 0; n < m_mesh->GetNumFaceEdges(f); ++n) - { - const auto nodeIndex = m_mesh->m_facesNodes[f][n]; - if (m_nodeMask[nodeIndex] != 0 && m_nodeMask[nodeIndex] != -2) - { - activeNodeFound = true; - break; - } - } - if (!activeNodeFound) - { - m_faceMask[f] = 0; - } - } - } - if (level > 0) - { - // if one face node is not in polygon disable refinement - for (UInt f = 0; f < m_mesh->GetNumFaces(); f++) - { - for (UInt n = 0; n < m_mesh->GetNumFaceEdges(f); n++) - { - const auto nodeIndex = m_mesh->m_facesNodes[f][n]; - if (m_nodeMask[nodeIndex] != 1) - { - m_faceMask[f] = 0; - break; - } - } - } - } - + ReviewRefinement(level); ComputeEdgesRefinementMask(); - ComputeIfFaceShouldBeSplit(); +#if 0 UInt numFacesToRefine = 0; for (UInt f = 0; f < m_mesh->GetNumFaces(); f++) { @@ -173,10 +214,22 @@ void MeshRefinement::Compute() numFacesToRefine++; } } - if (numFacesToRefine == 0) +#endif + + // UInt numFacesToRefine = std::count_if(m_faceMask.begin(), m_faceMask.end(), [](UInt value) + // { return value != 0; }); + + // if (numFacesToRefine == 0) + // { + // break; + // } + + if (std::count_if(m_faceMask.begin(), m_faceMask.end(), [](UInt value) + { return value != 0; }) == 0) { break; } + numFacesAfterRefinement = numFacesAfterRefinement * 4; // spit the edges From 71b938910602378c710cbfac99d21f11efc51cca Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 27 Nov 2023 11:19:50 +0100 Subject: [PATCH 02/15] GRIDEDIT-785 Eliminated the allocation of averaging strategy for each interpolation point --- .../include/MeshKernel/AveragingInterpolation.hpp | 3 ++- .../AveragingStrategies/AveragingStrategy.hpp | 3 +++ .../AveragingStrategyFactory.hpp | 2 -- .../ClosestAveragingStrategy.hpp | 8 +++++--- .../InverseWeightedAveragingStrategy.hpp | 11 ++++++----- .../AveragingStrategies/MaxAveragingStrategy.hpp | 3 +++ .../AveragingStrategies/MinAbsAveragingStrategy.hpp | 3 +++ .../AveragingStrategies/MinAveragingStrategy.hpp | 3 +++ .../AveragingStrategies/SimpleAveragingStrategy.hpp | 3 +++ libs/MeshKernel/src/AveragingInterpolation.cpp | 10 +++++----- .../AveragingStrategyFactory.cpp | 5 ++--- .../ClosestAveragingStrategy.cpp | 11 ++++++++--- .../InverseWeightedAveragingStrategy.cpp | 13 +++++++++---- .../AveragingStrategies/MaxAveragingStrategy.cpp | 6 ++++++ .../AveragingStrategies/MinAbsAveragingStrategy.cpp | 5 +++++ .../AveragingStrategies/MinAveragingStrategy.cpp | 5 +++++ .../AveragingStrategies/SimpleAveragingStrategy.cpp | 6 ++++++ libs/MeshKernel/src/MeshRefinement.cpp | 12 +++--------- .../MeshKernel/tests/src/AveragingStrategyTests.cpp | 8 ++++++-- 19 files changed, 83 insertions(+), 37 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp b/libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp index 02051294e..f8e5d205f 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp @@ -154,6 +154,7 @@ namespace meshkernel std::vector m_visitedSamples; ///< The visited samples - RTree m_samplesRtree; ///< The samples tree + RTree m_samplesRtree; ///< The samples tree + std::unique_ptr m_strategy; ///< Averaging strategy }; } // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp index b1b039be1..fc69d4029 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp @@ -36,6 +36,9 @@ namespace meshkernel::averaging public: virtual ~AveragingStrategy() = default; + /// @brief Reset the state of the averaging strategy, ready for the next calculation. + virtual void Reset(const Point& interpolationPoint) = 0; + /// @brief Adds the specified sample point. /// @param[in] samplePoint The sample point to add to this strategy. /// @param[in] sampleValue The value associated with the sample point. diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategyFactory.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategyFactory.hpp index 0757439bd..2dce6faec 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategyFactory.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategyFactory.hpp @@ -38,12 +38,10 @@ namespace meshkernel::averaging /// @brief The static method returning the strategy /// @param[in] averagingMethod The averaging method enumeration value /// @param[in] minNumSamples The minimum a of samples used for certain interpolation algorithms - /// @param[in] interpolationPoint The interpolation point /// @param[in] projection The projection to use /// @return The interpolation strategy to use [[nodiscard]] std::unique_ptr static GetAveragingStrategy(AveragingInterpolation::Method averagingMethod, size_t minNumSamples, - Point const& interpolationPoint, Projection projection); }; } // namespace meshkernel::averaging diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp index 7a992d71e..39a61291b 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp @@ -39,8 +39,10 @@ namespace meshkernel::averaging /// @brief Construct a new ClosestAveragingStrategy. /// @param[in] interpolationPoint The point for which the average should be calculated. /// @param[in] projection The projection used to calculate distances with. - ClosestAveragingStrategy(Point const& interpolationPoint, - Projection projection); + ClosestAveragingStrategy(Projection projection); + + /// @brief Reset the state of the closest averaging strategy. + void Reset(const Point& interpolationPoint) override; void Add(Point const& samplePoint, double sampleValue) override; [[nodiscard]] double Calculate() const override; @@ -53,7 +55,7 @@ namespace meshkernel::averaging double m_closestSquaredValue = std::numeric_limits::max(); /// @brief The interpolation point from which the closest value is calculated. - Point const& m_interpolationPoint; + Point m_interpolationPoint{constants::missing::doubleValue, constants::missing::doubleValue}; /// @brief The projection used to calculate the squared distance. Projection const m_projection; diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp index b38cbf269..9f3165bd0 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp @@ -37,13 +37,14 @@ namespace meshkernel::averaging { public: /// @brief Construct a new InverseWeightedAveragingStrategy. - /// @param[in] interpolationPoint The point for which the average should be calculated. /// @param[in] minNumSamples The minimum amount of samples for a valid interpolation. /// @param[in] projection The projection used in calculating the distance. - InverseWeightedAveragingStrategy(Point const& interpolationPoint, - size_t minNumSamples, + InverseWeightedAveragingStrategy(size_t minNumSamples, Projection projection); + /// @brief Reset the state of the inverse weighted averaging strategy. + void Reset(const Point& interpolationPoint) override; + void Add(Point const& samplePoint, double sampleValue) override; [[nodiscard]] double Calculate() const override; @@ -55,10 +56,10 @@ namespace meshkernel::averaging double m_wall = 0.0; /// @brief The interpolation point from which the inverse weight is calculated. - Point const& m_interpolationPoint; + Point m_interpolationPoint{constants::missing::doubleValue, constants::missing::doubleValue}; /// @brief The minimum number of samples for a valid interpolation. - size_t m_minNumSamples; + const size_t m_minNumSamples; /// @brief The projection used to calculate the distance. Projection const m_projection; diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MaxAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MaxAveragingStrategy.hpp index 04e2c8fc5..ff3b2a790 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MaxAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MaxAveragingStrategy.hpp @@ -36,6 +36,9 @@ namespace meshkernel::averaging class MaxAveragingStrategy final : public AveragingStrategy { public: + /// @brief Reset the state of the max averaging strategy. + void Reset(const Point& interpolationPoint) override; + void Add(Point const& samplePoint, double sampleValue) override; [[nodiscard]] double Calculate() const override; diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAbsAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAbsAveragingStrategy.hpp index 921a0f461..192eba7b0 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAbsAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAbsAveragingStrategy.hpp @@ -36,6 +36,9 @@ namespace meshkernel::averaging class MinAbsAveragingStrategy final : public AveragingStrategy { public: + /// @brief Reset the state of the absolute-value-of-the-mininimum averaging-strategy. + void Reset(const Point& interpolationPoint) override; + void Add(Point const& samplePoint, double sampleValue) override; [[nodiscard]] double Calculate() const override; diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAveragingStrategy.hpp index e2920d78d..5f09a4739 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAveragingStrategy.hpp @@ -35,6 +35,9 @@ namespace meshkernel::averaging class MinAveragingStrategy final : public AveragingStrategy { public: + /// @brief Reset the state of the minimum value averaging-strategy. + void Reset(const Point& interpolationPoint) override; + void Add(Point const& samplePoint, double sampleValue) override; [[nodiscard]] double Calculate() const override; diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp index 51df35b86..0849056cf 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp @@ -13,6 +13,9 @@ namespace meshkernel::averaging /// @param minNumSamples[in] The minimum amount of samples for a valid interpolation SimpleAveragingStrategy(size_t minNumSamples); + /// @brief Reset the state of the simple averaging-strategy. + void Reset(const Point& interpolationPoint) override; + void Add(Point const& samplePoint, double sampleValue) override; [[nodiscard]] double Calculate() const override; diff --git a/libs/MeshKernel/src/AveragingInterpolation.cpp b/libs/MeshKernel/src/AveragingInterpolation.cpp index 524b5e88a..9d756a775 100644 --- a/libs/MeshKernel/src/AveragingInterpolation.cpp +++ b/libs/MeshKernel/src/AveragingInterpolation.cpp @@ -49,7 +49,8 @@ AveragingInterpolation::AveragingInterpolation(Mesh2D& mesh, m_relativeSearchRadius(relativeSearchRadius), m_useClosestSampleIfNoneAvailable(useClosestSampleIfNoneAvailable), m_transformSamples(transformSamples), - m_minNumSamples(minNumSamples) + m_minNumSamples(minNumSamples), + m_strategy(averaging::AveragingStrategyFactory::GetAveragingStrategy(m_method, m_minNumSamples, m_mesh.m_projection)) { } @@ -199,8 +200,7 @@ double AveragingInterpolation::GetSampleValueFromRTree(UInt const index) double AveragingInterpolation::ComputeInterpolationResultFromNeighbors(const Point& interpolationPoint, std::vector const& searchPolygon) { - - auto strategy = averaging::AveragingStrategyFactory::GetAveragingStrategy(m_method, m_minNumSamples, interpolationPoint, m_mesh.m_projection); + m_strategy->Reset(interpolationPoint); for (UInt i = 0; i < m_samplesRtree.GetQueryResultSize(); ++i) { @@ -215,11 +215,11 @@ double AveragingInterpolation::ComputeInterpolationResultFromNeighbors(const Poi Point samplePoint{m_samples[sampleIndex].x, m_samples[sampleIndex].y}; if (IsPointInPolygonNodes(samplePoint, searchPolygon, m_mesh.m_projection)) { - strategy->Add(samplePoint, sampleValue); + m_strategy->Add(samplePoint, sampleValue); } } - return strategy->Calculate(); + return m_strategy->Calculate(); } double AveragingInterpolation::ComputeOnPolygon(const std::vector& polygon, diff --git a/libs/MeshKernel/src/AveragingStrategies/AveragingStrategyFactory.cpp b/libs/MeshKernel/src/AveragingStrategies/AveragingStrategyFactory.cpp index 359dfae5b..12e5bda20 100644 --- a/libs/MeshKernel/src/AveragingStrategies/AveragingStrategyFactory.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/AveragingStrategyFactory.cpp @@ -36,7 +36,6 @@ std::unique_ptr meshkernel::averaging::AveragingStrategyFactory::GetAveragingStrategy(AveragingInterpolation::Method averagingMethod, size_t minNumSamples, - Point const& interpolationPoint, Projection projection) { switch (averagingMethod) @@ -44,13 +43,13 @@ std::unique_ptr meshkernel::averaging: case AveragingInterpolation::Method::SimpleAveraging: return std::make_unique(minNumSamples); case AveragingInterpolation::Method::Closest: - return std::make_unique(interpolationPoint, projection); + return std::make_unique(projection); case AveragingInterpolation::Method::Max: return std::make_unique(); case AveragingInterpolation::Method::Min: return std::make_unique(); case AveragingInterpolation::Method::InverseWeightedDistance: - return std::make_unique(interpolationPoint, minNumSamples, projection); + return std::make_unique(minNumSamples, projection); case AveragingInterpolation::Method::MinAbsValue: return std::make_unique(); default: diff --git a/libs/MeshKernel/src/AveragingStrategies/ClosestAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/ClosestAveragingStrategy.cpp index 15a4d8166..4a47a0ddf 100644 --- a/libs/MeshKernel/src/AveragingStrategies/ClosestAveragingStrategy.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/ClosestAveragingStrategy.cpp @@ -30,9 +30,14 @@ namespace meshkernel::averaging { - ClosestAveragingStrategy::ClosestAveragingStrategy(Point const& interpolationPoint, - Projection const projection) : m_interpolationPoint(interpolationPoint), - m_projection(projection) {} + ClosestAveragingStrategy::ClosestAveragingStrategy(Projection const projection) : m_projection(projection) {} + + void ClosestAveragingStrategy::Reset(const Point& interpolationPoint) + { + m_interpolationPoint = interpolationPoint; + m_result = constants::missing::doubleValue; + m_closestSquaredValue = std::numeric_limits::max(); + } void ClosestAveragingStrategy::Add(Point const& samplePoint, double const sampleValue) { diff --git a/libs/MeshKernel/src/AveragingStrategies/InverseWeightedAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/InverseWeightedAveragingStrategy.cpp index 81724c98b..1ed4f1d72 100644 --- a/libs/MeshKernel/src/AveragingStrategies/InverseWeightedAveragingStrategy.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/InverseWeightedAveragingStrategy.cpp @@ -3,12 +3,17 @@ namespace meshkernel::averaging { - InverseWeightedAveragingStrategy::InverseWeightedAveragingStrategy(Point const& interpolation_point, - size_t minNumSamples, - Projection const projection) : m_interpolationPoint(interpolation_point), - m_minNumSamples(minNumSamples), + InverseWeightedAveragingStrategy::InverseWeightedAveragingStrategy(size_t minNumSamples, + Projection const projection) : m_minNumSamples(minNumSamples), m_projection(projection) {} + void InverseWeightedAveragingStrategy::Reset(const Point& interpolationPoint) + { + m_interpolationPoint = interpolationPoint; + m_result = 0.0; + m_wall = 0.0; + } + void InverseWeightedAveragingStrategy::Add(Point const& samplePoint, double const sampleValue) { double const distance = std::max(0.01, ComputeDistance(m_interpolationPoint, samplePoint, m_projection)); diff --git a/libs/MeshKernel/src/AveragingStrategies/MaxAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/MaxAveragingStrategy.cpp index 9f058eb1f..4a1d74309 100644 --- a/libs/MeshKernel/src/AveragingStrategies/MaxAveragingStrategy.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/MaxAveragingStrategy.cpp @@ -29,6 +29,12 @@ namespace meshkernel::averaging { + + void MaxAveragingStrategy::Reset(const Point& interpolationPoint [[maybe_unused]]) + { + m_result = std::numeric_limits::lowest(); + } + void MaxAveragingStrategy::Add(Point const& /*samplePoint*/, double const sampleValue) { m_result = std::max(m_result, sampleValue); diff --git a/libs/MeshKernel/src/AveragingStrategies/MinAbsAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/MinAbsAveragingStrategy.cpp index 7b1d62b7a..a7eb3e91d 100644 --- a/libs/MeshKernel/src/AveragingStrategies/MinAbsAveragingStrategy.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/MinAbsAveragingStrategy.cpp @@ -29,6 +29,11 @@ namespace meshkernel::averaging { + void MinAbsAveragingStrategy::Reset(const Point& interpolationPoint [[maybe_unused]]) + { + m_result = std::numeric_limits::max(); + } + void MinAbsAveragingStrategy::Add(Point const& /*samplePoint*/, double const sampleValue) { m_result = std::min(m_result, std::abs(sampleValue)); diff --git a/libs/MeshKernel/src/AveragingStrategies/MinAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/MinAveragingStrategy.cpp index d12458492..9e144dbd0 100644 --- a/libs/MeshKernel/src/AveragingStrategies/MinAveragingStrategy.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/MinAveragingStrategy.cpp @@ -29,6 +29,11 @@ namespace meshkernel::averaging { + + void MinAveragingStrategy::Reset(const Point& interpolationPoint [[maybe_unused]]) + { + m_result = std::numeric_limits::max(); + } void MinAveragingStrategy::Add(Point const& /*samplePoint*/, double const sampleValue) { m_result = std::min(m_result, sampleValue); diff --git a/libs/MeshKernel/src/AveragingStrategies/SimpleAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/SimpleAveragingStrategy.cpp index f66a0cd24..8b6e59d25 100644 --- a/libs/MeshKernel/src/AveragingStrategies/SimpleAveragingStrategy.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/SimpleAveragingStrategy.cpp @@ -32,6 +32,12 @@ namespace meshkernel::averaging SimpleAveragingStrategy::SimpleAveragingStrategy(size_t minNumSamples) : m_minNumPoints(minNumSamples) {} + void SimpleAveragingStrategy::Reset(const Point& interpolationPoint [[maybe_unused]]) + { + m_result = 0.0; + m_nAdds = 0; + } + void SimpleAveragingStrategy::Add(Point const& /*samplePoint*/, double const sampleValue) { m_result += sampleValue; diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index e9a9fff87..280092bf2 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -216,16 +216,10 @@ void MeshRefinement::Compute() } #endif - // UInt numFacesToRefine = std::count_if(m_faceMask.begin(), m_faceMask.end(), [](UInt value) - // { return value != 0; }); + auto notZero = [](UInt value) + { return value != 0; }; - // if (numFacesToRefine == 0) - // { - // break; - // } - - if (std::count_if(m_faceMask.begin(), m_faceMask.end(), [](UInt value) - { return value != 0; }) == 0) + if (std::count_if(m_faceMask.begin(), m_faceMask.end(), notZero) == 0) { break; } diff --git a/libs/MeshKernel/tests/src/AveragingStrategyTests.cpp b/libs/MeshKernel/tests/src/AveragingStrategyTests.cpp index a44d1221b..15244e799 100644 --- a/libs/MeshKernel/tests/src/AveragingStrategyTests.cpp +++ b/libs/MeshKernel/tests/src/AveragingStrategyTests.cpp @@ -31,7 +31,9 @@ namespace meshkernel::averaging Point p = Point(0.0, 0.0); std::unique_ptr pStrategy = AveragingStrategyFactory::GetAveragingStrategy(GetParam().first, 1, - p, Projection::cartesian); + Projection::cartesian); + + pStrategy->Reset(p); // Call double result = pStrategy->Calculate(); @@ -97,7 +99,9 @@ namespace meshkernel::averaging // Setup Point p = Point(0.0, 0.0); std::unique_ptr pStrategy = AveragingStrategyFactory::GetAveragingStrategy(GetParam().method_, 1, - p, Projection::cartesian); + Projection::cartesian); + + pStrategy->Reset(p); for (auto const& p : GetParam().addData_) { From d334d86fdb69887a8b6cedc32524ce8c83e54b38 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 27 Nov 2023 14:17:36 +0100 Subject: [PATCH 03/15] GRIDEDIT-785 Fixed bug when reviewing refinement --- .../SimpleAveragingStrategy.hpp | 2 +- .../include/MeshKernel/MeshRefinement.hpp | 2 +- libs/MeshKernel/src/MeshRefinement.cpp | 78 ++++++------------- 3 files changed, 27 insertions(+), 55 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp index 0849056cf..b2154ed38 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp @@ -21,7 +21,7 @@ namespace meshkernel::averaging private: /// @brief The minimum number of points for a valid interpolation. - size_t m_minNumPoints; + const size_t m_minNumPoints; /// @brief The current result from which Calculate calculates the final value. double m_result = 0.0; diff --git a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp index 916aaee0f..b6ed98f7b 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp @@ -60,7 +60,7 @@ namespace meshkernel class SamplesBasedRefinement : public ComputeRefinement { public: - // Will loop over all elememts in the mesh calling ComputeForFace + // Will loop over all elements in the mesh calling ComputeForFace void compute() const override; private: diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index 280092bf2..ed4b219d0 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -91,60 +91,46 @@ MeshRefinement::MeshRefinement(std::shared_ptr mesh, void MeshRefinement::ReviewRefinement(const int level) { - // TODO is this comment correct? Looks like its disabling refinement - // if one face node is in polygon enable face refinement - // TODO add enum for nodeMask values - for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) + if (level == 0) { - for (UInt n = 0; n < m_mesh->GetNumFaceEdges(f); ++n) + // TODO is this comment correct? Looks like its disabling refinement + // if one face node is in polygon enable face refinement + for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) { - const auto nodeIndex = m_mesh->m_facesNodes[f][n]; - - bool disableRefinement = (level == 0 ? m_nodeMask[nodeIndex] != 0 && m_nodeMask[nodeIndex] != -2 : m_nodeMask[nodeIndex] != 1); + bool activeNodeFound = false; - if (disableRefinement) + for (UInt n = 0; n < m_mesh->GetNumFaceEdges(f); ++n) { - m_faceMask[f] = 0; - break; + const auto nodeIndex = m_mesh->m_facesNodes[f][n]; + if (m_nodeMask[nodeIndex] != 0 && m_nodeMask[nodeIndex] != -2) + { + activeNodeFound = true; + break; + } } - } - } -#if 0 - if (level == 0) - { - // TODO is this comment correct? Looks like its disabling refinement - // if one face node is in polygon enable face refinement - for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) + if (!activeNodeFound) { - for (UInt n = 0; n < m_mesh->GetNumFaceEdges(f); ++n) - { - const auto nodeIndex = m_mesh->m_facesNodes[f][n]; - if (m_nodeMask[nodeIndex] != 0 && m_nodeMask[nodeIndex] != -2) - { - m_faceMask[f] = 0; - break; - } - } + m_faceMask[f] = 0; } } - if (level > 0) + } + if (level > 0) + { + // if one face node is not in polygon disable refinement + for (UInt f = 0; f < m_mesh->GetNumFaces(); f++) { - // if one face node is not in polygon disable refinement - for (UInt f = 0; f < m_mesh->GetNumFaces(); f++) + for (UInt n = 0; n < m_mesh->GetNumFaceEdges(f); n++) { - for (UInt n = 0; n < m_mesh->GetNumFaceEdges(f); n++) + const auto nodeIndex = m_mesh->m_facesNodes[f][n]; + if (m_nodeMask[nodeIndex] != 1) { - const auto nodeIndex = m_mesh->m_facesNodes[f][n]; - if (m_nodeMask[nodeIndex] != 1) - { - m_faceMask[f] = 0; - break; - } + m_faceMask[f] = 0; + break; } } } -#endif + } } void MeshRefinement::Compute() @@ -183,7 +169,6 @@ void MeshRefinement::Compute() // set_initial_mask ComputeNodeMaskAtPolygonPerimeter(); - auto numFacesAfterRefinement = m_mesh->GetNumFaces(); // reserve some extra capacity for the node mask m_nodeMask.reserve(m_nodeMask.size() * 2); for (int level = 0; level < m_meshRefinementParameters.max_num_refinement_iterations; level++) @@ -205,17 +190,6 @@ void MeshRefinement::Compute() ComputeEdgesRefinementMask(); ComputeIfFaceShouldBeSplit(); -#if 0 - UInt numFacesToRefine = 0; - for (UInt f = 0; f < m_mesh->GetNumFaces(); f++) - { - if (m_faceMask[f] != 0) - { - numFacesToRefine++; - } - } -#endif - auto notZero = [](UInt value) { return value != 0; }; @@ -224,8 +198,6 @@ void MeshRefinement::Compute() break; } - numFacesAfterRefinement = numFacesAfterRefinement * 4; - // spit the edges RefineFacesBySplittingEdges(); From 06d4802eaf1f48e467a76a81fb83d38b2ff1ff76 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 27 Nov 2023 16:20:34 +0100 Subject: [PATCH 04/15] GRIDEDIT-785 Removed some member variables, now local variables --- .../MeshKernel/AveragingInterpolation.hpp | 2 -- libs/MeshKernel/src/AveragingInterpolation.cpp | 16 ++++------------ 2 files changed, 4 insertions(+), 14 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp b/libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp index f8e5d205f..5f5860226 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp @@ -152,8 +152,6 @@ namespace meshkernel bool m_transformSamples = false; ///< Wheher to transform samples UInt m_minNumSamples = 1; ///< The minimum amount of samples for a valid interpolation. Used in some interpolation algorithms. - std::vector m_visitedSamples; ///< The visited samples - RTree m_samplesRtree; ///< The samples tree std::unique_ptr m_strategy; ///< Averaging strategy }; diff --git a/libs/MeshKernel/src/AveragingInterpolation.cpp b/libs/MeshKernel/src/AveragingInterpolation.cpp index 9d756a775..296e4b3ca 100644 --- a/libs/MeshKernel/src/AveragingInterpolation.cpp +++ b/libs/MeshKernel/src/AveragingInterpolation.cpp @@ -66,16 +66,10 @@ void AveragingInterpolation::Compute() m_samplesRtree.BuildTree(m_samples); } - if (m_visitedSamples.empty()) - { - m_visitedSamples.resize(m_samples.size()); - } - if (m_interpolationLocation == Mesh::Location::Nodes || m_interpolationLocation == Mesh::Location::Edges) { m_nodeResults.resize(m_mesh.GetNumNodes(), constants::missing::doubleValue); - std::ranges::fill(m_nodeResults, constants::missing::doubleValue); // make sure edge centers are computed m_mesh.ComputeEdgesCenters(); @@ -93,7 +87,6 @@ void AveragingInterpolation::Compute() if (m_interpolationLocation == Mesh::Location::Edges) { m_edgeResults.resize(m_mesh.GetNumEdges(), constants::missing::doubleValue); - std::ranges::fill(m_edgeResults, constants::missing::doubleValue); for (UInt e = 0; e < m_mesh.GetNumEdges(); ++e) { @@ -111,11 +104,10 @@ void AveragingInterpolation::Compute() if (m_interpolationLocation == Mesh::Location::Faces) { + std::vector visitedSamples(m_samples.size(), false); ///< The visited samples + std::vector polygonNodesCache(Mesh::m_maximumNumberOfNodesPerFace + 1); m_faceResults.resize(m_mesh.GetNumFaces(), constants::missing::doubleValue); - std::ranges::fill(m_faceResults, constants::missing::doubleValue); - std::vector polygonNodesCache(Mesh::m_maximumNumberOfNodesPerFace + 1); - std::fill(m_visitedSamples.begin(), m_visitedSamples.end(), false); for (UInt f = 0; f < m_mesh.GetNumFaces(); ++f) { polygonNodesCache.clear(); @@ -134,9 +126,9 @@ void AveragingInterpolation::Compute() // it is difficult to do it otherwise without sharing or caching the query result for (UInt i = 0; i < m_samplesRtree.GetQueryResultSize(); ++i) { - if (const auto sample = m_samplesRtree.GetQueryResult(i); !m_visitedSamples[sample]) + if (const auto sample = m_samplesRtree.GetQueryResult(i); !visitedSamples[sample]) { - m_visitedSamples[sample] = true; + visitedSamples[sample] = true; m_samples[sample].value -= 1; } } From fb263c78152615827e28565555e4f15001693bb0 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 27 Nov 2023 16:31:53 +0100 Subject: [PATCH 05/15] GRIDEDIT-785 Refactoring of mesh refinement class --- .../include/MeshKernel/MeshRefinement.hpp | 169 +++-- libs/MeshKernel/src/MeshRefinement.cpp | 669 +++++++++++------- 2 files changed, 524 insertions(+), 314 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp index b6ed98f7b..4993b8995 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp @@ -37,59 +37,33 @@ namespace meshkernel // Forward declarations class Mesh2D; -#if 0 - class ComputeRefinement { public: - virtual ~ComputeRefinement() = default; - - virtual void compute() const = 0; - - private: - }; - - class EdgeSizeBasedRefinement : public ComputeRefinement - { - public: - void compute() const override; - - private: - }; + ComputeRefinement(const Mesh2D& mesh) : m_mesh(mesh) {} - class SamplesBasedRefinement : public ComputeRefinement - { - public: - // Will loop over all elements in the mesh calling ComputeForFace - void compute() const override; + virtual ~ComputeRefinement() = default; - private: - virtual void ComputeForFace() const = 0; - }; + virtual void Compute(const std::vector& edgeIsBelowminimumSize, + const std::vector& brotherEdges, + std::vector& edgeMask, + std::vector& faceMask) = 0; - class WaveCourantRefinement : public SamplesBasedRefinement - { - public: - private: - void ComputeForFace(const UInt face) const override; - }; + virtual bool IsSampleBased() const + { + return false; + } - class RefinementLevelsRefinement : public SamplesBasedRefinement - { - public: - private: - void ComputeForFace(const UInt face) const override; - }; + protected: + const Mesh2D& GetMesh() const + { + return m_mesh; + } - class RidgeDetectionRefinement : public SamplesBasedRefinement - { - public: private: - void ComputeForFace(const UInt face) const override; + const Mesh2D& m_mesh; }; -#endif - /// @brief A class used to refine a Mesh2D instance /// /// Mesh refinement operates on Mesh2D and is based on @@ -126,6 +100,8 @@ namespace meshkernel /// existing Mesh2D instance. class MeshRefinement { + + public: /// @brief Enumerator describing the face location types enum class FaceLocation { @@ -147,7 +123,6 @@ namespace meshkernel /// If no such RefinementType value exists then a ConstraintError will be thrown. static RefinementType GetRefinementTypeValue(const int refinementInt); - public: /// @brief The constructor for refining based on samples /// @param[in] mesh The mesh to be refined /// @param[in] interpolant The averaging interpolation to use @@ -307,5 +282,113 @@ namespace meshkernel MeshRefinementParameters m_meshRefinementParameters; ///< The mesh refinement parameters bool m_useNodalRefinement = false; ///< Use refinement based on interpolated values at nodes const double m_mergingDistance = 0.001; ///< The distance for merging two edges + + std::unique_ptr m_refinementMasks; + }; + + class EdgeSizeBasedRefinement : public ComputeRefinement + { + public: + EdgeSizeBasedRefinement(Mesh2D& mesh) : ComputeRefinement(mesh) {} + + void Compute(const std::vector& edgeIsBelowMinimumSize, + const std::vector& brotherEdges, + std::vector& edgeMask, + std::vector& faceMask) override; + + private: + }; + + class SamplesBasedRefinement : public ComputeRefinement + { + public: + SamplesBasedRefinement(Mesh2D& mesh, + std::shared_ptr interpolant, + const MeshRefinementParameters& meshRefinementParameters) : ComputeRefinement(mesh), m_interpolant(interpolant), m_meshRefinementParameters(meshRefinementParameters) {} + + // Will loop over all elements in the mesh calling ComputeForFace + void Compute(const std::vector& edgeIsBelowminimumSize, + const std::vector& brotherEdges, + std::vector& edgeMask, + std::vector& faceMask) override; + + bool IsSampleBased() const override + { + return true; + } + + // TODO should all be private + std::shared_ptr m_interpolant; + MeshRefinementParameters m_meshRefinementParameters; ///< The mesh refinement parameters + + std::vector m_isHangingNodeCache; ///< Cache for maintaining if node is hanging + std::vector m_isHangingEdgeCache; ///< Cache for maintaining if edge is hanging + + private: + void SmoothRefinementMasks(const std::vector& brotherEdges, + std::vector& edgeMask, + std::vector& faceMask); + + // What to do here, it is needed by both the samples based refinement and the mesh-refinement + void FindHangingNodes(const std::vector& brotherEdges, UInt face); + + virtual UInt ComputeForFace(const UInt faceId, + const std::vector& edgeIsBelowMinimumSize, + std::vector& refineEdgeCache) = 0; + + std::vector m_refineEdgeCache; ///< Cache for the edges to be refined + }; + + class RefinementLevelsRefinement : public SamplesBasedRefinement + { + public: + RefinementLevelsRefinement(Mesh2D& mesh, + std::shared_ptr interpolant, + const MeshRefinementParameters& meshRefinementParameters) : SamplesBasedRefinement(mesh, interpolant, meshRefinementParameters) {} + + private: + UInt ComputeForFace(const UInt faceId, + const std::vector& edgeIsBelowMinimumSize, + std::vector& refineEdgeCache) override; }; + + class RidgeDetectionRefinement : public SamplesBasedRefinement + { + public: + RidgeDetectionRefinement(Mesh2D& mesh, + std::shared_ptr interpolant, + const MeshRefinementParameters& meshRefinementParameters) : SamplesBasedRefinement(mesh, interpolant, meshRefinementParameters) {} + + private: + UInt ComputeForFace(const UInt faceId, + const std::vector& edgeIsBelowMinimumSize, + std::vector& refineEdgeCache) override; + }; + + class WaveCourantRefinement : public SamplesBasedRefinement + { + public: + WaveCourantRefinement(Mesh2D& mesh, + std::shared_ptr interpolant, + + const MeshRefinementParameters& meshRefinementParameters, + bool useNodalRefinement) : SamplesBasedRefinement(mesh, interpolant, meshRefinementParameters), m_useNodalRefinement(useNodalRefinement) + { + } + + private: + bool IsRefineNeededBasedOnCourantCriteria(UInt edge, double depthValues) const; + + void ComputeFaceLocationTypes(); + + UInt ComputeForFace(const UInt faceId, + const std::vector& edgeIsBelowMinimumSize, + std::vector& refineEdgeCache) override; + + std::vector m_faceLocationType; ///< Cache for the face location types + + bool m_useNodalRefinement = false; ///< Use refinement based on interpolated values at nodes + const double m_mergingDistance = 0.001; ///< The distance for merging two edges + }; + } // namespace meshkernel diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index ed4b219d0..fbae03551 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -32,8 +32,401 @@ #include #include +using meshkernel::EdgeSizeBasedRefinement; using meshkernel::Mesh2D; using meshkernel::MeshRefinement; +using meshkernel::RefinementLevelsRefinement; +using meshkernel::RidgeDetectionRefinement; +using meshkernel::SamplesBasedRefinement; +using meshkernel::WaveCourantRefinement; + +void EdgeSizeBasedRefinement::Compute(const std::vector& edgeIsBelowMinimumSize, + const std::vector& brotherEdges [[maybe_unused]], + std::vector& edgeMask, + std::vector& faceMask) +{ + for (UInt f = 0; f < GetMesh().GetNumFaces(); ++f) + { + for (UInt n = 0; n < GetMesh().GetNumFaceEdges(f); ++n) + { + const auto e = GetMesh().m_facesEdges[f][n]; + if (!edgeIsBelowMinimumSize[e]) + { + edgeMask[e] = -1; + faceMask[f] = 1; + } + } + } +} + +void SamplesBasedRefinement::Compute(const std::vector& edgeIsBelowMinimumSize, + const std::vector& brotherEdges, + std::vector& edgeMask, + std::vector& faceMask) +{ + + m_refineEdgeCache.resize(Mesh::m_maximumNumberOfEdgesPerFace, 0); + + // Compute all interpolated values + m_interpolant->Compute(); + + for (UInt face = 0; face < GetMesh().GetNumFaces(); face++) + { + FindHangingNodes(brotherEdges, face); + UInt numberToBeRefined = ComputeForFace(face, edgeIsBelowMinimumSize, m_refineEdgeCache); + + if (numberToBeRefined > 1) + { + faceMask[face] = 1; + + for (size_t n = 0; n < GetMesh().GetNumFaceEdges(face); n++) + { + if (m_refineEdgeCache[n] == 1) + { + const auto edgeIndex = GetMesh().m_facesEdges[face][n]; + if (edgeIndex != constants::missing::uintValue) + { + edgeMask[edgeIndex] = 1; + } + } + } + } + } + + for (auto& edge : edgeMask) + { + edge = -edge; + } + + SmoothRefinementMasks(brotherEdges, edgeMask, faceMask); +} + +void SamplesBasedRefinement::FindHangingNodes(const std::vector& brotherEdges, UInt face) +{ + const auto numFaceNodes = GetMesh().GetNumFaceEdges(face); + + if (numFaceNodes > Mesh::m_maximumNumberOfEdgesPerNode) + { + throw AlgorithmError("The number of face nodes is greater than the maximum number of edges per node."); + } + + m_isHangingNodeCache.resize(Mesh::m_maximumNumberOfNodesPerFace); + m_isHangingEdgeCache.resize(Mesh::m_maximumNumberOfEdgesPerFace); + std::fill(m_isHangingNodeCache.begin(), m_isHangingNodeCache.end(), false); + std::fill(m_isHangingEdgeCache.begin(), m_isHangingEdgeCache.end(), false); + + auto kknod = numFaceNodes; + for (UInt n = 0; n < numFaceNodes; n++) + { + const auto edgeIndex = GetMesh().m_facesEdges[face][n]; + + // check if the parent edge is in the cell + if (brotherEdges[edgeIndex] != constants::missing::uintValue) + { + const auto e = NextCircularBackwardIndex(n, numFaceNodes); + const auto ee = NextCircularForwardIndex(n, numFaceNodes); + const auto firstEdgeIndex = GetMesh().m_facesEdges[face][e]; + const auto secondEdgeIndex = GetMesh().m_facesEdges[face][ee]; + + UInt commonNode = constants::missing::uintValue; + if (brotherEdges[edgeIndex] == firstEdgeIndex) + { + commonNode = GetMesh().FindCommonNode(edgeIndex, firstEdgeIndex); + if (commonNode == constants::missing::uintValue) + { + throw AlgorithmError("Could not find common node."); + } + } + else if (brotherEdges[edgeIndex] == secondEdgeIndex) + { + commonNode = GetMesh().FindCommonNode(edgeIndex, secondEdgeIndex); + if (commonNode == constants::missing::uintValue) + { + throw AlgorithmError("Could not find common node."); + } + } + + if (commonNode != constants::missing::uintValue) + { + m_isHangingEdgeCache[n] = true; + for (UInt nn = 0; nn < numFaceNodes; nn++) + { + kknod = NextCircularForwardIndex(kknod, numFaceNodes); + + if (GetMesh().m_facesNodes[face][kknod] == commonNode && !m_isHangingNodeCache[kknod]) + { + m_isHangingNodeCache[kknod] = true; + break; + } + } + } + } + } +} + +void SamplesBasedRefinement::SmoothRefinementMasks(const std::vector& brotherEdges, + std::vector& edgeMask, + std::vector& faceMask) +{ + if (m_meshRefinementParameters.directional_refinement == 1) + { + throw AlgorithmError("Directional refinement cannot be used in combination with smoothing. Please set directional refinement to off!"); + } + if (m_meshRefinementParameters.smoothing_iterations == 0) + { + return; + } + + std::vector splitEdge(edgeMask.size(), false); + + for (int iter = 0; iter < m_meshRefinementParameters.smoothing_iterations; ++iter) + { + std::fill(splitEdge.begin(), splitEdge.end(), false); + + for (UInt f = 0; f < GetMesh().GetNumFaces(); ++f) + { + if (faceMask[f] != 1) + { + continue; + } + + const auto numEdges = GetMesh().GetNumFaceEdges(f); + + for (UInt e = 0; e < numEdges; ++e) + { + const auto edgeIndex = GetMesh().m_facesEdges[f][e]; + const auto nextEdgeIndex = NextCircularForwardIndex(e, numEdges); + const auto previousEdgeIndex = NextCircularBackwardIndex(e, numEdges); + const auto split = brotherEdges[edgeIndex] != GetMesh().m_facesEdges[f][nextEdgeIndex] && + brotherEdges[edgeIndex] != GetMesh().m_facesEdges[f][previousEdgeIndex]; + + if (split) + { + splitEdge[edgeIndex] = true; + } + } + } + + // update face refinement mask + for (UInt f = 0; f < GetMesh().GetNumFaces(); ++f) + { + const auto numEdges = GetMesh().GetNumFaceEdges(f); + + for (UInt e = 0; e < numEdges; ++e) + { + const auto edgeIndex = GetMesh().m_facesEdges[f][e]; + if (splitEdge[edgeIndex]) + { + faceMask[f] = 1; + } + } + } + + // update edge refinement mask + for (UInt f = 0; f < GetMesh().GetNumFaces(); ++f) + { + if (faceMask[f] != 1) + { + continue; + } + + const auto numEdges = GetMesh().GetNumFaceEdges(f); + + for (UInt e = 0; e < numEdges; ++e) + { + const auto edgeIndex = GetMesh().m_facesEdges[f][e]; + const auto nextEdgeIndex = NextCircularForwardIndex(e, numEdges); + const auto previousEdgeIndex = NextCircularBackwardIndex(e, numEdges); + const auto split = brotherEdges[edgeIndex] != GetMesh().m_facesEdges[f][nextEdgeIndex] && + brotherEdges[edgeIndex] != GetMesh().m_facesEdges[f][previousEdgeIndex]; + + if (split) + { + edgeMask[edgeIndex] = 1; + } + } + } + } +} + +meshkernel::UInt RefinementLevelsRefinement::ComputeForFace(const UInt face, + const std::vector& edgeIsBelowMinimumSize [[maybe_unused]], + std::vector& edgeToRefine) +{ + if (m_interpolant->GetFaceResult(face) <= 0) + { + return 0; + } + + UInt numberOfEdgesToRefine = 0; + + for (UInt i = 0; i < GetMesh().GetNumFaceEdges(face); i++) + { + numberOfEdgesToRefine++; + edgeToRefine[i] = 1; + } + + return numberOfEdgesToRefine; +} + +bool WaveCourantRefinement::IsRefineNeededBasedOnCourantCriteria(UInt edge, double depthValues) const +{ + const double maxDtCourant = m_meshRefinementParameters.max_courant_time; + const double celerity = constants::physical::sqrt_gravity * std::sqrt(std::abs(depthValues)); + const double waveCourant = celerity * maxDtCourant / GetMesh().m_edgeLengths[edge]; + return waveCourant < 1.0; +} + +void WaveCourantRefinement::ComputeFaceLocationTypes() +{ + m_faceLocationType.resize(GetMesh().GetNumFaces()); + std::ranges::fill(m_faceLocationType, MeshRefinement::FaceLocation::Water); + for (UInt face = 0; face < GetMesh().GetNumFaces(); face++) + { + double maxVal = -std::numeric_limits::max(); + double minVal = std::numeric_limits::max(); + for (UInt e = 0; e < GetMesh().GetNumFaceEdges(face); ++e) + { + const auto node = GetMesh().m_facesNodes[face][e]; + const auto val = m_interpolant->GetNodeResult(node); + maxVal = std::max(maxVal, val); + minVal = std::min(minVal, val); + } + + if (minVal > 0.0) + { + m_faceLocationType[face] = MeshRefinement::FaceLocation::Land; + } + if (maxVal >= 0.0 && (!IsEqual(minVal, constants::missing::doubleValue) && minVal < 0.0)) + { + m_faceLocationType[face] = MeshRefinement::FaceLocation::LandWater; + } + } +} + +meshkernel::UInt WaveCourantRefinement::ComputeForFace(const UInt face, + const std::vector& edgeIsBelowMinimumSize [[maybe_unused]], + std::vector& edgeToRefine) +{ + UInt numberOfEdgesToRefine = 0; + + if (m_useNodalRefinement) + { + ComputeFaceLocationTypes(); + } + for (size_t e = 0; e < GetMesh().GetNumFaceEdges(face); ++e) + { + const auto edge = GetMesh().m_facesEdges[face][e]; + if (GetMesh().m_edgeLengths[edge] < m_mergingDistance) + { + numberOfEdgesToRefine++; + continue; + } + if (edgeIsBelowMinimumSize[edge]) + { + continue; + } + + bool doRefinement; + if (m_useNodalRefinement) + { + if (m_faceLocationType[face] == MeshRefinement::FaceLocation::Land) + { + doRefinement = false; + } + else if (m_faceLocationType[face] == MeshRefinement::FaceLocation::LandWater) + { + doRefinement = true; + } + else + { + const auto edgeDepth = std::abs(m_interpolant->GetEdgeResult(edge)); + doRefinement = IsRefineNeededBasedOnCourantCriteria(edge, edgeDepth); + } + } + else + { + const double faceDepth = m_interpolant->GetFaceResult(face); + doRefinement = IsRefineNeededBasedOnCourantCriteria(edge, faceDepth); + } + + if (doRefinement) + { + numberOfEdgesToRefine++; + edgeToRefine[e] = 1; + } + } + + if (numberOfEdgesToRefine > 0) + { + numberOfEdgesToRefine = 0; + for (UInt i = 0; i < GetMesh().GetNumFaceEdges(face); i++) + { + if (edgeToRefine[i] == 1 || m_isHangingNodeCache[i]) + { + numberOfEdgesToRefine++; + } + } + } + + if (m_meshRefinementParameters.directional_refinement == 0) + { + if (numberOfEdgesToRefine == GetMesh().GetNumFaceEdges(face)) + { + for (UInt i = 0; i < GetMesh().GetNumFaceEdges(face); i++) + { + if (!m_isHangingNodeCache[i]) + { + edgeToRefine[i] = 1; + } + } + } + else + { + numberOfEdgesToRefine = 0; + } + } + + return numberOfEdgesToRefine; +} + +meshkernel::UInt RidgeDetectionRefinement::ComputeForFace(const UInt face, + const std::vector& edgeIsBelowMinimumSize [[maybe_unused]], + std::vector& edgeToRefine) +{ + UInt numberOfEdgesToRefine = 0; + + double maxEdgeLength = 0.0; + UInt numEdges = GetMesh().GetNumFaceEdges(face); + + for (size_t i = 0; i < numEdges; ++i) + { + + const auto& edgeIndex = GetMesh().m_facesEdges[face][i]; + const auto& [firstNode, secondNode] = GetMesh().m_edges[edgeIndex]; + const auto distance = ComputeDistance(GetMesh().m_nodes[firstNode], GetMesh().m_nodes[secondNode], GetMesh().m_projection); + maxEdgeLength = std::max(maxEdgeLength, distance); + } + + double absInterpolatedFaceValue = std::abs(m_interpolant->GetFaceResult(face)); + double threshold = 1.0; // In Fortran code this value is 100 + double thresholdMin = 1.0; + double hmin = m_meshRefinementParameters.min_edge_size; + + double elementSizeWanted = threshold / (absInterpolatedFaceValue + 1.0e-8); + + if (maxEdgeLength > elementSizeWanted && maxEdgeLength > 2.0 * hmin && absInterpolatedFaceValue > thresholdMin) + { + numberOfEdgesToRefine += numEdges; + + for (UInt i = 0; i < numEdges; i++) + { + edgeToRefine[i] = 1; + } + } + + return numberOfEdgesToRefine; +} MeshRefinement::RefinementType MeshRefinement::GetRefinementTypeValue(const int refinementInt) { @@ -153,7 +546,7 @@ void MeshRefinement::Compute() } // select the nodes to refine - auto const isRefinementBasedOnSamples = m_interpolant != nullptr; + auto const isRefinementBasedOnSamples = m_refinementMasks->IsSampleBased(); // m_interpolant != nullptr; if (!isRefinementBasedOnSamples && m_meshRefinementParameters.refine_intersected == 1) { const auto edgeMask = m_mesh->MaskEdgesOfFacesInPolygon(m_polygons, false, true); @@ -176,15 +569,10 @@ void MeshRefinement::Compute() // Compute all edge lengths at once ComputeEdgeBelowMinSizeAfterRefinement(); - // computes the edge and face refinement mask from samples - if (isRefinementBasedOnSamples) - { - ComputeRefinementMasksFromSamples(); - } - else - { - ComputeRefinementMaskFromEdgeSize(); - } + std::ranges::fill(m_edgeMask, 0); + std::ranges::fill(m_faceMask, 0); + + m_refinementMasks->Compute(m_isEdgeBelowMinSizeAfterRefinement, m_brotherEdges, m_edgeMask, m_faceMask); ReviewRefinement(level); ComputeEdgesRefinementMask(); @@ -219,24 +607,6 @@ void MeshRefinement::Compute() } } -void MeshRefinement::ComputeRefinementMaskFromEdgeSize() -{ - std::ranges::fill(m_edgeMask, 0); - std::ranges::fill(m_faceMask, 0); - for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) - { - for (UInt n = 0; n < m_mesh->GetNumFaceEdges(f); ++n) - { - const auto e = m_mesh->m_facesEdges[f][n]; - if (!m_isEdgeBelowMinSizeAfterRefinement[e]) - { - m_edgeMask[e] = -1; - m_faceMask[f] = 1; - } - } - } -} - meshkernel::UInt MeshRefinement::DeleteIsolatedHangingnodes() { @@ -757,33 +1127,6 @@ void MeshRefinement::ComputeNodeMaskAtPolygonPerimeter() } } -void MeshRefinement::ComputeRefinementMasksFromSamples() -{ - std::ranges::fill(m_edgeMask, 0); - std::ranges::fill(m_faceMask, 0); - - m_polygonNodesCache.resize(Mesh::m_maximumNumberOfNodesPerFace + 1); - m_localNodeIndicesCache.resize(Mesh::m_maximumNumberOfNodesPerFace + 1, constants::missing::uintValue); - m_globalEdgeIndicesCache.resize(Mesh::m_maximumNumberOfEdgesPerFace + 1, constants::missing::uintValue); - m_refineEdgeCache.resize(Mesh::m_maximumNumberOfEdgesPerFace, 0); - - // Compute all interpolated values - m_interpolant->Compute(); - - for (UInt f = 0; f < m_mesh->GetNumFaces(); f++) - { - FindHangingNodes(f); - - ComputeRefinementMasksFromSamples(f); - } - - for (auto& edge : m_edgeMask) - { - edge = -edge; - } - SmoothRefinementMasks(); -} - void MeshRefinement::FindHangingNodes(UInt face) { const auto numFaceNodes = m_mesh->GetNumFaceEdges(face); @@ -890,214 +1233,6 @@ meshkernel::UInt MeshRefinement::CountEdgesToRefine(UInt face) const return result; } -void MeshRefinement::ComputeRefinementMasksForRefinementLevels(UInt face, - size_t& numberOfEdgesToRefine, - std::vector& edgeToRefine) const -{ - if (m_interpolant->GetFaceResult(face) <= 0) - { - return; - } - - for (UInt i = 0; i < m_mesh->GetNumFaceEdges(face); i++) - { - numberOfEdgesToRefine++; - edgeToRefine[i] = 1; - } -} - -void MeshRefinement::ComputeRefinementMasksForWaveCourant(UInt face, - size_t& numberOfEdgesToRefine, - std::vector& edgeToRefine) -{ - if (m_useNodalRefinement) - { - ComputeFaceLocationTypes(); - } - for (size_t e = 0; e < m_mesh->GetNumFaceEdges(face); ++e) - { - const auto edge = m_mesh->m_facesEdges[face][e]; - if (m_mesh->m_edgeLengths[edge] < m_mergingDistance) - { - numberOfEdgesToRefine++; - continue; - } - if (m_isEdgeBelowMinSizeAfterRefinement[edge]) - { - continue; - } - - bool doRefinement; - if (m_useNodalRefinement) - { - if (m_faceLocationType[face] == FaceLocation::Land) - { - doRefinement = false; - } - else if (m_faceLocationType[face] == FaceLocation::LandWater) - { - doRefinement = true; - } - else - { - const auto edgeDepth = std::abs(m_interpolant->GetEdgeResult(edge)); - doRefinement = IsRefineNeededBasedOnCourantCriteria(edge, edgeDepth); - } - } - else - { - const double faceDepth = m_interpolant->GetFaceResult(face); - doRefinement = IsRefineNeededBasedOnCourantCriteria(edge, faceDepth); - } - - if (doRefinement) - { - numberOfEdgesToRefine++; - edgeToRefine[e] = 1; - } - } - - if (numberOfEdgesToRefine > 0) - { - numberOfEdgesToRefine = 0; - for (UInt i = 0; i < m_mesh->GetNumFaceEdges(face); i++) - { - if (edgeToRefine[i] == 1 || m_isHangingNodeCache[i]) - { - numberOfEdgesToRefine++; - } - } - } - - if (m_meshRefinementParameters.directional_refinement == 0) - { - if (numberOfEdgesToRefine == m_mesh->GetNumFaceEdges(face)) - { - for (UInt i = 0; i < m_mesh->GetNumFaceEdges(face); i++) - { - if (!m_isHangingNodeCache[i]) - { - edgeToRefine[i] = 1; - } - } - } - else - { - numberOfEdgesToRefine = 0; - } - } -} - -void MeshRefinement::ComputeRefinementMasksForRidgeDetection(UInt face, - size_t& numberOfEdgesToRefine, - std::vector& edgeToRefine) const -{ - - double maxEdgeLength = 0.0; - UInt numEdges = m_mesh->GetNumFaceEdges(face); - - for (size_t i = 0; i < numEdges; ++i) - { - - const auto& edgeIndex = m_mesh->m_facesEdges[face][i]; - const auto& [firstNode, secondNode] = m_mesh->m_edges[edgeIndex]; - const auto distance = ComputeDistance(m_mesh->m_nodes[firstNode], m_mesh->m_nodes[secondNode], m_mesh->m_projection); - maxEdgeLength = std::max(maxEdgeLength, distance); - } - - double absInterpolatedFaceValue = std::abs(m_interpolant->GetFaceResult(face)); - double threshold = 1.0; // In Fortran code this value is 100 - double thresholdMin = 1.0; - double hmin = m_meshRefinementParameters.min_edge_size; - - double elementSizeWanted = threshold / (absInterpolatedFaceValue + 1.0e-8); - - if (maxEdgeLength > elementSizeWanted && maxEdgeLength > 2.0 * hmin && absInterpolatedFaceValue > thresholdMin) - { - numberOfEdgesToRefine += numEdges; - - for (UInt i = 0; i < numEdges; i++) - { - edgeToRefine[i] = 1; - } - } -} - -void MeshRefinement::ComputeRefinementMasksFromSamples(UInt face) -{ - - if (IsEqual(m_interpolant->GetFaceResult(face), constants::missing::doubleValue)) - { - return; - } - - size_t numEdgesToBeRefined = 0; - std::ranges::fill(m_refineEdgeCache, 0); - - switch (m_refinementType) - { - case RefinementType::RefinementLevels: - ComputeRefinementMasksForRefinementLevels(face, numEdgesToBeRefined, m_refineEdgeCache); - break; - - case RefinementType::WaveCourant: - ComputeRefinementMasksForWaveCourant(face, numEdgesToBeRefined, m_refineEdgeCache); - break; - - case RefinementType::RidgeDetection: - ComputeRefinementMasksForRidgeDetection(face, numEdgesToBeRefined, m_refineEdgeCache); - break; - - default: - throw AlgorithmError("Invalid refinement type"); - } - - // Compute face and edge masks - if (numEdgesToBeRefined > 1) - { - m_faceMask[face] = 1; - - for (size_t n = 0; n < m_mesh->GetNumFaceEdges(face); n++) - { - if (m_refineEdgeCache[n] == 1) - { - const auto edgeIndex = m_mesh->m_facesEdges[face][n]; - if (edgeIndex != constants::missing::uintValue) - { - m_edgeMask[edgeIndex] = 1; - } - } - } - } -} - -void MeshRefinement::ComputeFaceLocationTypes() -{ - m_faceLocationType.resize(m_mesh->GetNumFaces()); - std::ranges::fill(m_faceLocationType, FaceLocation::Water); - for (UInt face = 0; face < m_mesh->GetNumFaces(); face++) - { - double maxVal = -std::numeric_limits::max(); - double minVal = std::numeric_limits::max(); - for (UInt e = 0; e < m_mesh->GetNumFaceEdges(face); ++e) - { - const auto node = m_mesh->m_facesNodes[face][e]; - const auto val = m_interpolant->GetNodeResult(node); - maxVal = std::max(maxVal, val); - minVal = std::min(minVal, val); - } - - if (minVal > 0.0) - { - m_faceLocationType[face] = FaceLocation::Land; - } - if (maxVal >= 0.0 && (!IsEqual(minVal, constants::missing::doubleValue) && minVal < 0.0)) - { - m_faceLocationType[face] = FaceLocation::LandWater; - } - } -} - void MeshRefinement::ComputeEdgeBelowMinSizeAfterRefinement() { m_mesh->ComputeEdgesLengths(); @@ -1109,14 +1244,6 @@ void MeshRefinement::ComputeEdgeBelowMinSizeAfterRefinement() } } -bool MeshRefinement::IsRefineNeededBasedOnCourantCriteria(UInt edge, double depthValues) const -{ - const double maxDtCourant = m_meshRefinementParameters.max_courant_time; - const double celerity = constants::physical::sqrt_gravity * std::sqrt(std::abs(depthValues)); - const double waveCourant = celerity * maxDtCourant / m_mesh->m_edgeLengths[edge]; - return waveCourant < 1.0; -} - void MeshRefinement::ComputeEdgesRefinementMask() { bool repeat = true; From 4069e8bdcbb5f960c998eb75436c0c55f5813804 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 27 Nov 2023 17:57:50 +0100 Subject: [PATCH 06/15] GRIDEDIT-785 (Temporarily) Allocated a refinement mask calculator in mesh refinement constructor --- .../include/MeshKernel/MeshRefinement.hpp | 6 +-- libs/MeshKernel/src/MeshRefinement.cpp | 48 ++++++++++++++++++- 2 files changed, 48 insertions(+), 6 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp index 4993b8995..b86b6b7cd 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp @@ -44,7 +44,7 @@ namespace meshkernel virtual ~ComputeRefinement() = default; - virtual void Compute(const std::vector& edgeIsBelowminimumSize, + virtual void Compute(const std::vector& edgeIsBelowMinimumSize, const std::vector& brotherEdges, std::vector& edgeMask, std::vector& faceMask) = 0; @@ -100,7 +100,6 @@ namespace meshkernel /// existing Mesh2D instance. class MeshRefinement { - public: /// @brief Enumerator describing the face location types enum class FaceLocation @@ -307,7 +306,7 @@ namespace meshkernel const MeshRefinementParameters& meshRefinementParameters) : ComputeRefinement(mesh), m_interpolant(interpolant), m_meshRefinementParameters(meshRefinementParameters) {} // Will loop over all elements in the mesh calling ComputeForFace - void Compute(const std::vector& edgeIsBelowminimumSize, + void Compute(const std::vector& edgeIsBelowMinimumSize, const std::vector& brotherEdges, std::vector& edgeMask, std::vector& faceMask) override; @@ -370,7 +369,6 @@ namespace meshkernel public: WaveCourantRefinement(Mesh2D& mesh, std::shared_ptr interpolant, - const MeshRefinementParameters& meshRefinementParameters, bool useNodalRefinement) : SamplesBasedRefinement(mesh, interpolant, meshRefinementParameters), m_useNodalRefinement(useNodalRefinement) { diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index fbae03551..51b1316ad 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -456,8 +456,29 @@ MeshRefinement::MeshRefinement(std::shared_ptr mesh, CheckMeshRefinementParameters(meshRefinementParameters); m_meshRefinementParameters = meshRefinementParameters; m_refinementType = GetRefinementTypeValue(m_meshRefinementParameters.refinement_type); -} + if (interpolant != nullptr) + { + switch (m_refinementType) + { + case RefinementType::WaveCourant: + m_refinementMasks = std::make_unique(*mesh, interpolant, meshRefinementParameters, false); + break; + case RefinementType::RefinementLevels: + m_refinementMasks = std::make_unique(*mesh, interpolant, meshRefinementParameters); + break; + case RefinementType::RidgeDetection: + m_refinementMasks = std::make_unique(*mesh, interpolant, meshRefinementParameters); + break; + default: + throw ConstraintError("Invalid mesh refinement type when creating with interpolant: refinement type {}", m_meshRefinementParameters.refinement_type); + } + } + else + { + m_refinementMasks = std::make_unique(*mesh); + } +} MeshRefinement::MeshRefinement(std::shared_ptr mesh, std::shared_ptr interpolant, const MeshRefinementParameters& meshRefinementParameters, @@ -469,6 +490,28 @@ MeshRefinement::MeshRefinement(std::shared_ptr mesh, CheckMeshRefinementParameters(meshRefinementParameters); m_meshRefinementParameters = meshRefinementParameters; m_refinementType = GetRefinementTypeValue(m_meshRefinementParameters.refinement_type); + + if (interpolant != nullptr) + { + switch (m_refinementType) + { + case RefinementType::WaveCourant: + m_refinementMasks = std::make_unique(*mesh, interpolant, meshRefinementParameters, useNodalRefinement); + break; + case RefinementType::RefinementLevels: + m_refinementMasks = std::make_unique(*mesh, interpolant, meshRefinementParameters); + break; + case RefinementType::RidgeDetection: + m_refinementMasks = std::make_unique(*mesh, interpolant, meshRefinementParameters); + break; + default: + throw ConstraintError("Invalid mesh refinement type when creating with interpolant: refinement type {}", m_meshRefinementParameters.refinement_type); + } + } + else + { + m_refinementMasks = std::make_unique(*mesh); + } } MeshRefinement::MeshRefinement(std::shared_ptr mesh, @@ -479,6 +522,7 @@ MeshRefinement::MeshRefinement(std::shared_ptr mesh, { CheckMeshRefinementParameters(meshRefinementParameters); m_meshRefinementParameters = meshRefinementParameters; + m_refinementMasks = std::make_unique(*mesh); } void MeshRefinement::ReviewRefinement(const int level) @@ -546,7 +590,7 @@ void MeshRefinement::Compute() } // select the nodes to refine - auto const isRefinementBasedOnSamples = m_refinementMasks->IsSampleBased(); // m_interpolant != nullptr; + auto const isRefinementBasedOnSamples = m_refinementMasks != nullptr && m_refinementMasks->IsSampleBased(); // m_interpolant != nullptr; if (!isRefinementBasedOnSamples && m_meshRefinementParameters.refine_intersected == 1) { const auto edgeMask = m_mesh->MaskEdgesOfFacesInPolygon(m_polygons, false, true); From 3aa50de748839409237e049212e0a108a7a12319 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Tue, 28 Nov 2023 13:45:40 +0100 Subject: [PATCH 07/15] GRIDEDIT-785 Fixed bug introduced by refactoring --- .../include/MeshKernel/MeshRefinement.hpp | 23 ++++----- libs/MeshKernel/src/MeshRefinement.cpp | 47 +++++++++++++++++-- 2 files changed, 55 insertions(+), 15 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp index b86b6b7cd..27e0cdd42 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp @@ -37,12 +37,12 @@ namespace meshkernel // Forward declarations class Mesh2D; - class ComputeRefinement + class ComputeRefinementMask { public: - ComputeRefinement(const Mesh2D& mesh) : m_mesh(mesh) {} + ComputeRefinementMask(const Mesh2D& mesh) : m_mesh(mesh) {} - virtual ~ComputeRefinement() = default; + virtual ~ComputeRefinementMask() = default; virtual void Compute(const std::vector& edgeIsBelowMinimumSize, const std::vector& brotherEdges, @@ -282,13 +282,13 @@ namespace meshkernel bool m_useNodalRefinement = false; ///< Use refinement based on interpolated values at nodes const double m_mergingDistance = 0.001; ///< The distance for merging two edges - std::unique_ptr m_refinementMasks; + std::unique_ptr m_refinementMasks; }; - class EdgeSizeBasedRefinement : public ComputeRefinement + class EdgeSizeBasedRefinement : public ComputeRefinementMask { public: - EdgeSizeBasedRefinement(Mesh2D& mesh) : ComputeRefinement(mesh) {} + EdgeSizeBasedRefinement(Mesh2D& mesh) : ComputeRefinementMask(mesh) {} void Compute(const std::vector& edgeIsBelowMinimumSize, const std::vector& brotherEdges, @@ -298,12 +298,14 @@ namespace meshkernel private: }; - class SamplesBasedRefinement : public ComputeRefinement + class SamplesBasedRefinement : public ComputeRefinementMask { public: SamplesBasedRefinement(Mesh2D& mesh, std::shared_ptr interpolant, - const MeshRefinementParameters& meshRefinementParameters) : ComputeRefinement(mesh), m_interpolant(interpolant), m_meshRefinementParameters(meshRefinementParameters) {} + const MeshRefinementParameters& meshRefinementParameters) : ComputeRefinementMask(mesh), + m_interpolant(interpolant), + m_meshRefinementParameters(meshRefinementParameters) {} // Will loop over all elements in the mesh calling ComputeForFace void Compute(const std::vector& edgeIsBelowMinimumSize, @@ -370,9 +372,8 @@ namespace meshkernel WaveCourantRefinement(Mesh2D& mesh, std::shared_ptr interpolant, const MeshRefinementParameters& meshRefinementParameters, - bool useNodalRefinement) : SamplesBasedRefinement(mesh, interpolant, meshRefinementParameters), m_useNodalRefinement(useNodalRefinement) - { - } + bool useNodalRefinement) : SamplesBasedRefinement(mesh, interpolant, meshRefinementParameters), + m_useNodalRefinement(useNodalRefinement) {} private: bool IsRefineNeededBasedOnCourantCriteria(UInt edge, double depthValues) const; diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index 51b1316ad..11a6c75bd 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -73,6 +73,13 @@ void SamplesBasedRefinement::Compute(const std::vector& edgeIsBelowMinimum for (UInt face = 0; face < GetMesh().GetNumFaces(); face++) { FindHangingNodes(brotherEdges, face); + + if (IsEqual(m_interpolant->GetFaceResult(face), constants::missing::doubleValue)) + { + continue; + } + + std::ranges::fill(m_refineEdgeCache, 0); UInt numberToBeRefined = ComputeForFace(face, edgeIsBelowMinimumSize, m_refineEdgeCache); if (numberToBeRefined > 1) @@ -93,6 +100,18 @@ void SamplesBasedRefinement::Compute(const std::vector& edgeIsBelowMinimum } } + // for (UInt i = 0; i < edgeIsBelowMinimumSize.size (); ++i) { + // std::cout << "edgeIsBelowMinimumSize " << i << " " << edgeIsBelowMinimumSize [i] << std::endl; + // } + + // for (UInt i = 0; i < faceMask.size (); ++i) { + // std::cout << " face " << i << " "<< faceMask[i] << std::endl; + // } + + for (UInt i = 0; i < edgeMask.size (); ++i) { + std::cout << " edge " << i << " "<< edgeMask[i] << std::endl; + } + for (auto& edge : edgeMask) { edge = -edge; @@ -310,6 +329,11 @@ meshkernel::UInt WaveCourantRefinement::ComputeForFace(const UInt face, { UInt numberOfEdgesToRefine = 0; + if (face == 288) { + [[maybe_unused]] int dummy; + dummy = 1; + } + if (m_useNodalRefinement) { ComputeFaceLocationTypes(); @@ -387,6 +411,8 @@ meshkernel::UInt WaveCourantRefinement::ComputeForFace(const UInt face, } } + std::cout << " face to refine " << face << " " << numberOfEdgesToRefine << std::endl; + return numberOfEdgesToRefine; } @@ -622,14 +648,27 @@ void MeshRefinement::Compute() ComputeEdgesRefinementMask(); ComputeIfFaceShouldBeSplit(); - auto notZero = [](UInt value) - { return value != 0; }; - - if (std::count_if(m_faceMask.begin(), m_faceMask.end(), notZero) == 0) + UInt numFacesToRefine = 0; + for (UInt f = 0; f < m_mesh->GetNumFaces(); f++) + { + if (m_faceMask[f] != 0) + { + numFacesToRefine++; + } + } + if (numFacesToRefine == 0) { break; } + // auto notZero = [](UInt value) + // { return value != 0; }; + + // if (std::count_if(m_faceMask.begin(), m_faceMask.end(), notZero) == 0) + // { + // break; + // } + // spit the edges RefineFacesBySplittingEdges(); From 9dd8619684b338e9ccb64f6f9129e5751bce38e7 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Tue, 28 Nov 2023 16:30:03 +0100 Subject: [PATCH 08/15] GRIDEDIT-785 FIrst step in removing state in the averaging strategy objects --- .../MeshKernel/AveragingInterpolation.hpp | 5 +++++ .../AveragingStrategies/AveragingStrategy.hpp | 15 +++++++++++++ .../ClosestAveragingStrategy.hpp | 4 ++++ .../InverseWeightedAveragingStrategy.hpp | 4 ++++ .../MeshKernel/src/AveragingInterpolation.cpp | 13 +++++++++--- .../ClosestAveragingStrategy.cpp | 21 +++++++++++++++++++ .../InverseWeightedAveragingStrategy.cpp | 18 ++++++++++++++++ libs/MeshKernel/src/MeshRefinement.cpp | 14 ------------- 8 files changed, 77 insertions(+), 17 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp b/libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp index 5f5860226..a50b33a06 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp @@ -109,6 +109,9 @@ namespace meshkernel void Compute() override; private: + /// @brief Default size for interpolation points and sample caches. + static constexpr UInt DefaultMaximumCacheSize = 100; + /// @brief Compute the averaging results in polygon /// @param[in] polygon The bounding polygon where the samples are included /// @param[in] interpolationPoint The interpolation point @@ -151,6 +154,8 @@ namespace meshkernel bool m_useClosestSampleIfNoneAvailable = false; ///< Whether to use the closest sample if there is none available bool m_transformSamples = false; ///< Wheher to transform samples UInt m_minNumSamples = 1; ///< The minimum amount of samples for a valid interpolation. Used in some interpolation algorithms. + std::vector m_interpolationPointCache; ///< Cache for interpolation points + std::vector m_interpolationSampleCache; ///< Cache for interpolation samples RTree m_samplesRtree; ///< The samples tree std::unique_ptr m_strategy; ///< Averaging strategy diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp index fc69d4029..8df4da2b9 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp @@ -47,5 +47,20 @@ namespace meshkernel::averaging /// @brief Calculates the average value based on the values added. /// @return The calculated average [[nodiscard]] virtual double Calculate() const = 0; + + virtual double Calculate (const Point& interpolationPoint, + const std::vector& samplePoints, + const std::vector& sampleValues) { + + + Reset (interpolationPoint); + + for (UInt i = 0; i < samplePoints.size (); ++i) + { + Add(samplePoints[i], sampleValues[i]); + } + + return Calculate (); + } }; } // namespace meshkernel::averaging diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp index 39a61291b..da305843e 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp @@ -47,6 +47,10 @@ namespace meshkernel::averaging void Add(Point const& samplePoint, double sampleValue) override; [[nodiscard]] double Calculate() const override; + double Calculate (const Point& interpolationPoint, + const std::vector& samplePoints, + const std::vector& sampleValues) override; //const + private: /// @brief The result used to calculate the final value in Calculate. double m_result = constants::missing::doubleValue; diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp index 9f3165bd0..3f81a51a7 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp @@ -48,6 +48,10 @@ namespace meshkernel::averaging void Add(Point const& samplePoint, double sampleValue) override; [[nodiscard]] double Calculate() const override; + double Calculate (const Point& interpolationPoint, + const std::vector& samplePoints, + const std::vector& sampleValues) override; //const + private: /// @brief The current result used in Calculate to calculate the final value. double m_result = 0.0; diff --git a/libs/MeshKernel/src/AveragingInterpolation.cpp b/libs/MeshKernel/src/AveragingInterpolation.cpp index 296e4b3ca..131cac245 100644 --- a/libs/MeshKernel/src/AveragingInterpolation.cpp +++ b/libs/MeshKernel/src/AveragingInterpolation.cpp @@ -52,6 +52,8 @@ AveragingInterpolation::AveragingInterpolation(Mesh2D& mesh, m_minNumSamples(minNumSamples), m_strategy(averaging::AveragingStrategyFactory::GetAveragingStrategy(m_method, m_minNumSamples, m_mesh.m_projection)) { + m_interpolationPointCache.reserve (DefaultMaximumCacheSize); + m_interpolationSampleCache.reserve (DefaultMaximumCacheSize); } void AveragingInterpolation::Compute() @@ -192,7 +194,10 @@ double AveragingInterpolation::GetSampleValueFromRTree(UInt const index) double AveragingInterpolation::ComputeInterpolationResultFromNeighbors(const Point& interpolationPoint, std::vector const& searchPolygon) { - m_strategy->Reset(interpolationPoint); + // m_strategy->Reset(interpolationPoint); + + m_interpolationPointCache.resize (0); + m_interpolationSampleCache.resize (0); for (UInt i = 0; i < m_samplesRtree.GetQueryResultSize(); ++i) { @@ -207,11 +212,13 @@ double AveragingInterpolation::ComputeInterpolationResultFromNeighbors(const Poi Point samplePoint{m_samples[sampleIndex].x, m_samples[sampleIndex].y}; if (IsPointInPolygonNodes(samplePoint, searchPolygon, m_mesh.m_projection)) { - m_strategy->Add(samplePoint, sampleValue); + m_interpolationPointCache.emplace_back (samplePoint); + m_interpolationSampleCache.emplace_back (sampleValue); + // m_strategy->Add(samplePoint, sampleValue); } } - return m_strategy->Calculate(); + return m_strategy->Calculate(interpolationPoint, m_interpolationPointCache, m_interpolationSampleCache); } double AveragingInterpolation::ComputeOnPolygon(const std::vector& polygon, diff --git a/libs/MeshKernel/src/AveragingStrategies/ClosestAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/ClosestAveragingStrategy.cpp index 4a47a0ddf..c293efbef 100644 --- a/libs/MeshKernel/src/AveragingStrategies/ClosestAveragingStrategy.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/ClosestAveragingStrategy.cpp @@ -53,4 +53,25 @@ namespace meshkernel::averaging { return m_result; } + + double ClosestAveragingStrategy::Calculate (const Point& interpolationPoint, + const std::vector& samplePoints, + const std::vector& sampleValues) //const + { + double result = constants::missing::doubleValue; + double closestSquaredValue = std::numeric_limits::max(); + + for (UInt i = 0; i < samplePoints.size (); ++i) + { + if (const auto squaredDistance = ComputeSquaredDistance(interpolationPoint, samplePoints[i], m_projection); + squaredDistance < closestSquaredValue) + { + closestSquaredValue = squaredDistance; + result = sampleValues[i]; + } + } + + return result; + } + } // namespace meshkernel::averaging diff --git a/libs/MeshKernel/src/AveragingStrategies/InverseWeightedAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/InverseWeightedAveragingStrategy.cpp index 1ed4f1d72..17cd71833 100644 --- a/libs/MeshKernel/src/AveragingStrategies/InverseWeightedAveragingStrategy.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/InverseWeightedAveragingStrategy.cpp @@ -26,4 +26,22 @@ namespace meshkernel::averaging { return m_wall >= m_minNumSamples ? m_result / m_wall : constants::missing::doubleValue; } + + double InverseWeightedAveragingStrategy::Calculate (const Point& interpolationPoint, + const std::vector& samplePoints, + const std::vector& sampleValues) //const + { + double result = 0.0; + double wall = 0.0; + + for (UInt i = 0; i < samplePoints.size (); ++i) + { + double weight = 1.0 / std::max(0.01, ComputeDistance(interpolationPoint, samplePoints[i], m_projection)); + wall += weight; + result += weight * sampleValues[i]; + } + + return wall >= static_cast(m_minNumSamples) ? result / wall : constants::missing::doubleValue; + } + } // namespace meshkernel::averaging diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index 11a6c75bd..4cde97f3b 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -100,18 +100,6 @@ void SamplesBasedRefinement::Compute(const std::vector& edgeIsBelowMinimum } } - // for (UInt i = 0; i < edgeIsBelowMinimumSize.size (); ++i) { - // std::cout << "edgeIsBelowMinimumSize " << i << " " << edgeIsBelowMinimumSize [i] << std::endl; - // } - - // for (UInt i = 0; i < faceMask.size (); ++i) { - // std::cout << " face " << i << " "<< faceMask[i] << std::endl; - // } - - for (UInt i = 0; i < edgeMask.size (); ++i) { - std::cout << " edge " << i << " "<< edgeMask[i] << std::endl; - } - for (auto& edge : edgeMask) { edge = -edge; @@ -411,8 +399,6 @@ meshkernel::UInt WaveCourantRefinement::ComputeForFace(const UInt face, } } - std::cout << " face to refine " << face << " " << numberOfEdgesToRefine << std::endl; - return numberOfEdgesToRefine; } From bb6442476a8b860c9194d47d7fc795c51567ffb1 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Tue, 28 Nov 2023 16:48:58 +0100 Subject: [PATCH 09/15] GRIDEDIT-785 Added stateless Calculate to all averaging strategy classes --- .../AveragingStrategies/AveragingStrategy.hpp | 15 +++------------ .../ClosestAveragingStrategy.hpp | 2 +- .../InverseWeightedAveragingStrategy.hpp | 2 +- .../AveragingStrategies/MaxAveragingStrategy.hpp | 4 ++++ .../MinAbsAveragingStrategy.hpp | 4 ++++ .../AveragingStrategies/MinAveragingStrategy.hpp | 4 ++++ .../SimpleAveragingStrategy.hpp | 4 ++++ .../ClosestAveragingStrategy.cpp | 2 +- .../InverseWeightedAveragingStrategy.cpp | 2 +- .../AveragingStrategies/MaxAveragingStrategy.cpp | 15 +++++++++++++++ .../MinAbsAveragingStrategy.cpp | 15 +++++++++++++++ .../AveragingStrategies/MinAveragingStrategy.cpp | 15 +++++++++++++++ .../SimpleAveragingStrategy.cpp | 15 +++++++++++++++ 13 files changed, 83 insertions(+), 16 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp index 8df4da2b9..97b58a6da 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp @@ -48,19 +48,10 @@ namespace meshkernel::averaging /// @return The calculated average [[nodiscard]] virtual double Calculate() const = 0; + /// @brief Calculates the average value based on the values added. + /// @return The calculated average virtual double Calculate (const Point& interpolationPoint, const std::vector& samplePoints, - const std::vector& sampleValues) { - - - Reset (interpolationPoint); - - for (UInt i = 0; i < samplePoints.size (); ++i) - { - Add(samplePoints[i], sampleValues[i]); - } - - return Calculate (); - } + const std::vector& sampleValues) const = 0; }; } // namespace meshkernel::averaging diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp index da305843e..7999bcf4a 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp @@ -49,7 +49,7 @@ namespace meshkernel::averaging double Calculate (const Point& interpolationPoint, const std::vector& samplePoints, - const std::vector& sampleValues) override; //const + const std::vector& sampleValues) const override; private: /// @brief The result used to calculate the final value in Calculate. diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp index 3f81a51a7..bae973bca 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp @@ -50,7 +50,7 @@ namespace meshkernel::averaging double Calculate (const Point& interpolationPoint, const std::vector& samplePoints, - const std::vector& sampleValues) override; //const + const std::vector& sampleValues) const override; private: /// @brief The current result used in Calculate to calculate the final value. diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MaxAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MaxAveragingStrategy.hpp index ff3b2a790..803812127 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MaxAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MaxAveragingStrategy.hpp @@ -43,6 +43,10 @@ namespace meshkernel::averaging [[nodiscard]] double Calculate() const override; + double Calculate (const Point& interpolationPoint, + const std::vector& samplePoints, + const std::vector& sampleValues) const override; + private: /// @brief The current result returned in Calculate. double m_result = std::numeric_limits::lowest(); diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAbsAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAbsAveragingStrategy.hpp index 192eba7b0..ca01c3971 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAbsAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAbsAveragingStrategy.hpp @@ -43,6 +43,10 @@ namespace meshkernel::averaging [[nodiscard]] double Calculate() const override; + double Calculate (const Point& interpolationPoint, + const std::vector& samplePoints, + const std::vector& sampleValues) const override; + private: /// @brief The current result returned in Calculate double m_result = std::numeric_limits::max(); diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAveragingStrategy.hpp index 5f09a4739..1b6d82c81 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAveragingStrategy.hpp @@ -42,6 +42,10 @@ namespace meshkernel::averaging [[nodiscard]] double Calculate() const override; + double Calculate (const Point& interpolationPoint, + const std::vector& samplePoints, + const std::vector& sampleValues) const override; + private: /// @brief The current result returned in Calculate double m_result = std::numeric_limits::max(); diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp index b2154ed38..b0df92f8f 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp @@ -19,6 +19,10 @@ namespace meshkernel::averaging void Add(Point const& samplePoint, double sampleValue) override; [[nodiscard]] double Calculate() const override; + double Calculate (const Point& interpolationPoint, + const std::vector& samplePoints, + const std::vector& sampleValues) const override; + private: /// @brief The minimum number of points for a valid interpolation. const size_t m_minNumPoints; diff --git a/libs/MeshKernel/src/AveragingStrategies/ClosestAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/ClosestAveragingStrategy.cpp index c293efbef..0dce1cb0f 100644 --- a/libs/MeshKernel/src/AveragingStrategies/ClosestAveragingStrategy.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/ClosestAveragingStrategy.cpp @@ -56,7 +56,7 @@ namespace meshkernel::averaging double ClosestAveragingStrategy::Calculate (const Point& interpolationPoint, const std::vector& samplePoints, - const std::vector& sampleValues) //const + const std::vector& sampleValues) const { double result = constants::missing::doubleValue; double closestSquaredValue = std::numeric_limits::max(); diff --git a/libs/MeshKernel/src/AveragingStrategies/InverseWeightedAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/InverseWeightedAveragingStrategy.cpp index 17cd71833..d2c1adc7d 100644 --- a/libs/MeshKernel/src/AveragingStrategies/InverseWeightedAveragingStrategy.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/InverseWeightedAveragingStrategy.cpp @@ -29,7 +29,7 @@ namespace meshkernel::averaging double InverseWeightedAveragingStrategy::Calculate (const Point& interpolationPoint, const std::vector& samplePoints, - const std::vector& sampleValues) //const + const std::vector& sampleValues) const { double result = 0.0; double wall = 0.0; diff --git a/libs/MeshKernel/src/AveragingStrategies/MaxAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/MaxAveragingStrategy.cpp index 4a1d74309..aa82fe891 100644 --- a/libs/MeshKernel/src/AveragingStrategies/MaxAveragingStrategy.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/MaxAveragingStrategy.cpp @@ -44,4 +44,19 @@ namespace meshkernel::averaging { return m_result != std::numeric_limits::lowest() ? m_result : constants::missing::doubleValue; } + + double MaxAveragingStrategy::Calculate(const Point& interpolationPoint [[maybe_unused]], + const std::vector& samplePoints, + const std::vector& sampleValues) const + { + double result = std::numeric_limits::lowest(); + + for (UInt i = 0; i < samplePoints.size (); ++i) + { + result = std::max(result, sampleValues[i]); + } + + return result != std::numeric_limits::lowest() ? result : constants::missing::doubleValue; + } + } // namespace meshkernel::averaging diff --git a/libs/MeshKernel/src/AveragingStrategies/MinAbsAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/MinAbsAveragingStrategy.cpp index a7eb3e91d..96348bcaf 100644 --- a/libs/MeshKernel/src/AveragingStrategies/MinAbsAveragingStrategy.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/MinAbsAveragingStrategy.cpp @@ -43,4 +43,19 @@ namespace meshkernel::averaging { return m_result != std::numeric_limits::max() ? m_result : constants::missing::doubleValue; } + + double MinAbsAveragingStrategy::Calculate(const Point& interpolationPoint [[maybe_unused]], + const std::vector& samplePoints, + const std::vector& sampleValues) const + { + double result = std::numeric_limits::max(); + + for (UInt i = 0; i < samplePoints.size (); ++i) + { + result = std::min(result, std::abs(sampleValues[i])); + } + + return result != std::numeric_limits::max() ? result : constants::missing::doubleValue; + } + } // namespace meshkernel::averaging diff --git a/libs/MeshKernel/src/AveragingStrategies/MinAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/MinAveragingStrategy.cpp index 9e144dbd0..daf93538d 100644 --- a/libs/MeshKernel/src/AveragingStrategies/MinAveragingStrategy.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/MinAveragingStrategy.cpp @@ -43,4 +43,19 @@ namespace meshkernel::averaging { return m_result != std::numeric_limits::max() ? m_result : constants::missing::doubleValue; } + + double MinAveragingStrategy::Calculate(const Point& interpolationPoint [[maybe_unused]], + const std::vector& samplePoints, + const std::vector& sampleValues) const + { + double result = std::numeric_limits::max(); + + for (UInt i = 0; i < samplePoints.size (); ++i) + { + result = std::min(result, sampleValues[i]); + } + + return result != std::numeric_limits::max() ? result : constants::missing::doubleValue; + } + } // namespace meshkernel::averaging diff --git a/libs/MeshKernel/src/AveragingStrategies/SimpleAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/SimpleAveragingStrategy.cpp index 8b6e59d25..91e67505b 100644 --- a/libs/MeshKernel/src/AveragingStrategies/SimpleAveragingStrategy.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/SimpleAveragingStrategy.cpp @@ -48,4 +48,19 @@ namespace meshkernel::averaging { return m_nAdds >= m_minNumPoints ? m_result / static_cast(m_nAdds) : constants::missing::doubleValue; } + + double SimpleAveragingStrategy::Calculate(const Point& interpolationPoint [[maybe_unused]], + const std::vector& samplePoints [[maybe_unused]], + const std::vector& sampleValues) const + { + double result = 0.0; + + for (UInt i = 0; i < sampleValues.size (); ++i) + { + result += sampleValues[i]; + } + + return sampleValues.size () >= m_minNumPoints ? result / static_cast(sampleValues.size ()) : constants::missing::doubleValue; + } + } // namespace meshkernel::averaging From 50e51e6a75dd338de138c19e79541683fe625320 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Tue, 28 Nov 2023 16:55:17 +0100 Subject: [PATCH 10/15] GRIDEDIT-785 Added doxygen comment to calculate function --- .../MeshKernel/AveragingStrategies/AveragingStrategy.hpp | 5 ++++- .../AveragingStrategies/ClosestAveragingStrategy.hpp | 5 +++++ .../AveragingStrategies/InverseWeightedAveragingStrategy.hpp | 5 +++++ .../MeshKernel/AveragingStrategies/MaxAveragingStrategy.hpp | 5 +++++ .../AveragingStrategies/MinAbsAveragingStrategy.hpp | 5 +++++ .../MeshKernel/AveragingStrategies/MinAveragingStrategy.hpp | 5 +++++ .../AveragingStrategies/SimpleAveragingStrategy.hpp | 5 +++++ libs/MeshKernel/src/AveragingInterpolation.cpp | 4 +--- 8 files changed, 35 insertions(+), 4 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp index 97b58a6da..ad9aaad7b 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/AveragingStrategy.hpp @@ -48,7 +48,10 @@ namespace meshkernel::averaging /// @return The calculated average [[nodiscard]] virtual double Calculate() const = 0; - /// @brief Calculates the average value based on the values added. + /// @brief Calculates the average value based on the sample values. + /// @param[in] interpolationPoint The point for which the average should be calculated. + /// @param[in] samplePoints The sample points to used by this strategy. + /// @param[in] sampleValues The sample values associated with each sample point. /// @return The calculated average virtual double Calculate (const Point& interpolationPoint, const std::vector& samplePoints, diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp index 7999bcf4a..f17b936c3 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp @@ -47,6 +47,11 @@ namespace meshkernel::averaging void Add(Point const& samplePoint, double sampleValue) override; [[nodiscard]] double Calculate() const override; + /// @brief Calculates the average value based on the sample values. + /// @param[in] interpolationPoint The point for which the average should be calculated. + /// @param[in] samplePoints The sample points to used by this strategy. + /// @param[in] sampleValues The sample values associated with each sample point. + /// @return The calculated average double Calculate (const Point& interpolationPoint, const std::vector& samplePoints, const std::vector& sampleValues) const override; diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp index bae973bca..2463f9e80 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp @@ -48,6 +48,11 @@ namespace meshkernel::averaging void Add(Point const& samplePoint, double sampleValue) override; [[nodiscard]] double Calculate() const override; + /// @brief Calculates the average value based on the sample values. + /// @param[in] interpolationPoint The point for which the average should be calculated. + /// @param[in] samplePoints The sample points to used by this strategy. + /// @param[in] sampleValues The sample values associated with each sample point. + /// @return The calculated average double Calculate (const Point& interpolationPoint, const std::vector& samplePoints, const std::vector& sampleValues) const override; diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MaxAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MaxAveragingStrategy.hpp index 803812127..d8f12cad8 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MaxAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MaxAveragingStrategy.hpp @@ -43,6 +43,11 @@ namespace meshkernel::averaging [[nodiscard]] double Calculate() const override; + /// @brief Calculates the average value based on the sample values. + /// @param[in] interpolationPoint The point for which the average should be calculated. + /// @param[in] samplePoints The sample points to used by this strategy. + /// @param[in] sampleValues The sample values associated with each sample point. + /// @return The calculated average double Calculate (const Point& interpolationPoint, const std::vector& samplePoints, const std::vector& sampleValues) const override; diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAbsAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAbsAveragingStrategy.hpp index ca01c3971..92feb07e1 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAbsAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAbsAveragingStrategy.hpp @@ -43,6 +43,11 @@ namespace meshkernel::averaging [[nodiscard]] double Calculate() const override; + /// @brief Calculates the average value based on the sample values. + /// @param[in] interpolationPoint The point for which the average should be calculated. + /// @param[in] samplePoints The sample points to used by this strategy. + /// @param[in] sampleValues The sample values associated with each sample point. + /// @return The calculated average double Calculate (const Point& interpolationPoint, const std::vector& samplePoints, const std::vector& sampleValues) const override; diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAveragingStrategy.hpp index 1b6d82c81..d40353a33 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAveragingStrategy.hpp @@ -42,6 +42,11 @@ namespace meshkernel::averaging [[nodiscard]] double Calculate() const override; + /// @brief Calculates the average value based on the sample values. + /// @param[in] interpolationPoint The point for which the average should be calculated. + /// @param[in] samplePoints The sample points to used by this strategy. + /// @param[in] sampleValues The sample values associated with each sample point. + /// @return The calculated average double Calculate (const Point& interpolationPoint, const std::vector& samplePoints, const std::vector& sampleValues) const override; diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp index b0df92f8f..e0c954d4a 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp @@ -19,6 +19,11 @@ namespace meshkernel::averaging void Add(Point const& samplePoint, double sampleValue) override; [[nodiscard]] double Calculate() const override; + /// @brief Calculates the average value based on the sample values. + /// @param[in] interpolationPoint The point for which the average should be calculated. + /// @param[in] samplePoints The sample points to used by this strategy. + /// @param[in] sampleValues The sample values associated with each sample point. + /// @return The calculated average double Calculate (const Point& interpolationPoint, const std::vector& samplePoints, const std::vector& sampleValues) const override; diff --git a/libs/MeshKernel/src/AveragingInterpolation.cpp b/libs/MeshKernel/src/AveragingInterpolation.cpp index 131cac245..a2c109804 100644 --- a/libs/MeshKernel/src/AveragingInterpolation.cpp +++ b/libs/MeshKernel/src/AveragingInterpolation.cpp @@ -194,8 +194,6 @@ double AveragingInterpolation::GetSampleValueFromRTree(UInt const index) double AveragingInterpolation::ComputeInterpolationResultFromNeighbors(const Point& interpolationPoint, std::vector const& searchPolygon) { - // m_strategy->Reset(interpolationPoint); - m_interpolationPointCache.resize (0); m_interpolationSampleCache.resize (0); @@ -210,11 +208,11 @@ double AveragingInterpolation::ComputeInterpolationResultFromNeighbors(const Poi } Point samplePoint{m_samples[sampleIndex].x, m_samples[sampleIndex].y}; + if (IsPointInPolygonNodes(samplePoint, searchPolygon, m_mesh.m_projection)) { m_interpolationPointCache.emplace_back (samplePoint); m_interpolationSampleCache.emplace_back (sampleValue); - // m_strategy->Add(samplePoint, sampleValue); } } From fe20718a96cb0d2514320ff71467407239c5af0b Mon Sep 17 00:00:00 2001 From: BillSenior Date: Tue, 28 Nov 2023 17:10:55 +0100 Subject: [PATCH 11/15] GRIDEDIT-785 Simplified the simple averaging strategy --- .../AveragingStrategies/SimpleAveragingStrategy.cpp | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/libs/MeshKernel/src/AveragingStrategies/SimpleAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/SimpleAveragingStrategy.cpp index 91e67505b..e4cba2b38 100644 --- a/libs/MeshKernel/src/AveragingStrategies/SimpleAveragingStrategy.cpp +++ b/libs/MeshKernel/src/AveragingStrategies/SimpleAveragingStrategy.cpp @@ -25,6 +25,8 @@ // //------------------------------------------------------------------------------ +#include + #include namespace meshkernel::averaging @@ -53,13 +55,7 @@ namespace meshkernel::averaging const std::vector& samplePoints [[maybe_unused]], const std::vector& sampleValues) const { - double result = 0.0; - - for (UInt i = 0; i < sampleValues.size (); ++i) - { - result += sampleValues[i]; - } - + double result = std::accumulate (sampleValues.begin (), sampleValues.end (), 0.0); return sampleValues.size () >= m_minNumPoints ? result / static_cast(sampleValues.size ()) : constants::missing::doubleValue; } From c20a260d19fe03b323b5c4364c4bd0f194b1200f Mon Sep 17 00:00:00 2001 From: BillSenior Date: Wed, 29 Nov 2023 11:32:56 +0100 Subject: [PATCH 12/15] GRIDEDIT-785 Simplified the creation of the mesh refinement mask calculator --- .../include/MeshKernel/MeshRefinement.hpp | 19 +++-- libs/MeshKernel/src/MeshRefinement.cpp | 69 ++++++++----------- 2 files changed, 42 insertions(+), 46 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp index 27e0cdd42..127389f52 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp @@ -164,6 +164,15 @@ namespace meshkernel void Compute(); private: + /// @brief Allocate the refinement mask calculator based on the RefinementType enum + /// + /// \note A ConstraintError exception may be thrown given inconsistent parameters + static std::unique_ptr CreateRefinementMaskCalculator(const RefinementType refinementType, + const Mesh2D& mesh, + std::shared_ptr interpolant, + const MeshRefinementParameters& meshRefinementParameters, + const bool useNodalRefinement); + /// @brief Finds if two edges are brothers, sharing an hanging node. Can be moved to Mesh2D void FindBrotherEdges(); @@ -288,7 +297,7 @@ namespace meshkernel class EdgeSizeBasedRefinement : public ComputeRefinementMask { public: - EdgeSizeBasedRefinement(Mesh2D& mesh) : ComputeRefinementMask(mesh) {} + EdgeSizeBasedRefinement(const Mesh2D& mesh) : ComputeRefinementMask(mesh) {} void Compute(const std::vector& edgeIsBelowMinimumSize, const std::vector& brotherEdges, @@ -301,7 +310,7 @@ namespace meshkernel class SamplesBasedRefinement : public ComputeRefinementMask { public: - SamplesBasedRefinement(Mesh2D& mesh, + SamplesBasedRefinement(const Mesh2D& mesh, std::shared_ptr interpolant, const MeshRefinementParameters& meshRefinementParameters) : ComputeRefinementMask(mesh), m_interpolant(interpolant), @@ -343,7 +352,7 @@ namespace meshkernel class RefinementLevelsRefinement : public SamplesBasedRefinement { public: - RefinementLevelsRefinement(Mesh2D& mesh, + RefinementLevelsRefinement(const Mesh2D& mesh, std::shared_ptr interpolant, const MeshRefinementParameters& meshRefinementParameters) : SamplesBasedRefinement(mesh, interpolant, meshRefinementParameters) {} @@ -356,7 +365,7 @@ namespace meshkernel class RidgeDetectionRefinement : public SamplesBasedRefinement { public: - RidgeDetectionRefinement(Mesh2D& mesh, + RidgeDetectionRefinement(const Mesh2D& mesh, std::shared_ptr interpolant, const MeshRefinementParameters& meshRefinementParameters) : SamplesBasedRefinement(mesh, interpolant, meshRefinementParameters) {} @@ -369,7 +378,7 @@ namespace meshkernel class WaveCourantRefinement : public SamplesBasedRefinement { public: - WaveCourantRefinement(Mesh2D& mesh, + WaveCourantRefinement(const Mesh2D& mesh, std::shared_ptr interpolant, const MeshRefinementParameters& meshRefinementParameters, bool useNodalRefinement) : SamplesBasedRefinement(mesh, interpolant, meshRefinementParameters), diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index 4cde97f3b..f1fb9fae8 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -317,11 +317,6 @@ meshkernel::UInt WaveCourantRefinement::ComputeForFace(const UInt face, { UInt numberOfEdgesToRefine = 0; - if (face == 288) { - [[maybe_unused]] int dummy; - dummy = 1; - } - if (m_useNodalRefinement) { ComputeFaceLocationTypes(); @@ -459,38 +454,51 @@ MeshRefinement::RefinementType MeshRefinement::GetRefinementTypeValue(const int } } -MeshRefinement::MeshRefinement(std::shared_ptr mesh, - std::shared_ptr interpolant, - const MeshRefinementParameters& meshRefinementParameters) - : m_mesh(mesh), - m_interpolant(interpolant) +std::unique_ptr MeshRefinement::CreateRefinementMaskCalculator(const RefinementType refinementType, + const Mesh2D& mesh, + std::shared_ptr interpolant, + const MeshRefinementParameters& meshRefinementParameters, + const bool useNodalRefinement) { - CheckMeshRefinementParameters(meshRefinementParameters); - m_meshRefinementParameters = meshRefinementParameters; - m_refinementType = GetRefinementTypeValue(m_meshRefinementParameters.refinement_type); + std::unique_ptr refinementMaskCalculator; if (interpolant != nullptr) { - switch (m_refinementType) + switch (refinementType) { case RefinementType::WaveCourant: - m_refinementMasks = std::make_unique(*mesh, interpolant, meshRefinementParameters, false); + refinementMaskCalculator = std::make_unique(mesh, interpolant, meshRefinementParameters, useNodalRefinement); break; case RefinementType::RefinementLevels: - m_refinementMasks = std::make_unique(*mesh, interpolant, meshRefinementParameters); + refinementMaskCalculator = std::make_unique(mesh, interpolant, meshRefinementParameters); break; case RefinementType::RidgeDetection: - m_refinementMasks = std::make_unique(*mesh, interpolant, meshRefinementParameters); + refinementMaskCalculator = std::make_unique(mesh, interpolant, meshRefinementParameters); break; default: - throw ConstraintError("Invalid mesh refinement type when creating with interpolant: refinement type {}", m_meshRefinementParameters.refinement_type); + throw ConstraintError("Invalid mesh refinement type when creating with interpolant: refinement type {}", meshRefinementParameters.refinement_type); } } else { - m_refinementMasks = std::make_unique(*mesh); + refinementMaskCalculator = std::make_unique(mesh); } + + return refinementMaskCalculator; } + +MeshRefinement::MeshRefinement(std::shared_ptr mesh, + std::shared_ptr interpolant, + const MeshRefinementParameters& meshRefinementParameters) + : m_mesh(mesh), + m_interpolant(interpolant) +{ + CheckMeshRefinementParameters(meshRefinementParameters); + m_meshRefinementParameters = meshRefinementParameters; + m_refinementType = GetRefinementTypeValue(m_meshRefinementParameters.refinement_type); + m_refinementMasks = CreateRefinementMaskCalculator(m_refinementType, *mesh, interpolant, meshRefinementParameters, false); +} + MeshRefinement::MeshRefinement(std::shared_ptr mesh, std::shared_ptr interpolant, const MeshRefinementParameters& meshRefinementParameters, @@ -502,28 +510,7 @@ MeshRefinement::MeshRefinement(std::shared_ptr mesh, CheckMeshRefinementParameters(meshRefinementParameters); m_meshRefinementParameters = meshRefinementParameters; m_refinementType = GetRefinementTypeValue(m_meshRefinementParameters.refinement_type); - - if (interpolant != nullptr) - { - switch (m_refinementType) - { - case RefinementType::WaveCourant: - m_refinementMasks = std::make_unique(*mesh, interpolant, meshRefinementParameters, useNodalRefinement); - break; - case RefinementType::RefinementLevels: - m_refinementMasks = std::make_unique(*mesh, interpolant, meshRefinementParameters); - break; - case RefinementType::RidgeDetection: - m_refinementMasks = std::make_unique(*mesh, interpolant, meshRefinementParameters); - break; - default: - throw ConstraintError("Invalid mesh refinement type when creating with interpolant: refinement type {}", m_meshRefinementParameters.refinement_type); - } - } - else - { - m_refinementMasks = std::make_unique(*mesh); - } + m_refinementMasks = CreateRefinementMaskCalculator(m_refinementType, *mesh, interpolant, meshRefinementParameters, useNodalRefinement); } MeshRefinement::MeshRefinement(std::shared_ptr mesh, From a7c7a54a0261fe0c4af99149e7df95691891c907 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Wed, 29 Nov 2023 12:10:52 +0100 Subject: [PATCH 13/15] GRIDEDIT-785 Reduced cyclomatic complexity in refinement smoothing --- .../include/MeshKernel/MeshRefinement.hpp | 9 ++ libs/MeshKernel/src/MeshRefinement.cpp | 129 ++++++++++-------- 2 files changed, 80 insertions(+), 58 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp index 127389f52..ab29a250b 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp @@ -233,6 +233,15 @@ namespace meshkernel /// @brief Connect the hanging nodes with triangles (connect_hanging_nodes) void ConnectHangingNodes(); + /// @brief Update edge refinement indicator to ensure smoother refinement transition between elements. + void SmoothEdgeRefinement(std::vector& splitEdge); + + /// @brief Using the latest edge refinement indicator, update face refinement mask. + void UpdateFaceRefinementMask(const std::vector& splitEdge); + + /// @brief Update the edge refinement mask + void UpdateEdgeRefinementMask(); + /// @brief Smooth the face and edge refinement masks (smooth_jarefine) void SmoothRefinementMasks(); diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index f1fb9fae8..33ca9625f 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -1594,85 +1594,98 @@ void MeshRefinement::FindBrotherEdges() } } -void MeshRefinement::SmoothRefinementMasks() +void MeshRefinement::SmoothEdgeRefinement(std::vector& splitEdge) { - if (m_meshRefinementParameters.directional_refinement == 1) - { - throw AlgorithmError("Directional refinement cannot be used in combination with smoothing. Please set directional refinement to off!"); - } - if (m_meshRefinementParameters.smoothing_iterations == 0) - { - return; - } - std::vector splitEdge(m_edgeMask.size(), false); - - for (int iter = 0; iter < m_meshRefinementParameters.smoothing_iterations; ++iter) + for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) { - std::fill(splitEdge.begin(), splitEdge.end(), false); - - for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) + if (m_faceMask[f] != 1) { - if (m_faceMask[f] != 1) - { - continue; - } + continue; + } - const auto numEdges = m_mesh->GetNumFaceEdges(f); + const auto numEdges = m_mesh->GetNumFaceEdges(f); - for (UInt e = 0; e < numEdges; ++e) - { - const auto edgeIndex = m_mesh->m_facesEdges[f][e]; - const auto nextEdgeIndex = NextCircularForwardIndex(e, numEdges); - const auto previousEdgeIndex = NextCircularBackwardIndex(e, numEdges); - const auto split = m_brotherEdges[edgeIndex] != m_mesh->m_facesEdges[f][nextEdgeIndex] && - m_brotherEdges[edgeIndex] != m_mesh->m_facesEdges[f][previousEdgeIndex]; + for (UInt e = 0; e < numEdges; ++e) + { + const auto edgeIndex = m_mesh->m_facesEdges[f][e]; + const auto nextEdgeIndex = NextCircularForwardIndex(e, numEdges); + const auto previousEdgeIndex = NextCircularBackwardIndex(e, numEdges); + const auto split = m_brotherEdges[edgeIndex] != m_mesh->m_facesEdges[f][nextEdgeIndex] && + m_brotherEdges[edgeIndex] != m_mesh->m_facesEdges[f][previousEdgeIndex]; - if (split) - { - splitEdge[edgeIndex] = true; - } + if (split) + { + splitEdge[edgeIndex] = true; } } + } +} - // update face refinement mask - for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) +void MeshRefinement::UpdateFaceRefinementMask(const std::vector& splitEdge) +{ + for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) + { + const auto numEdges = m_mesh->GetNumFaceEdges(f); + + for (UInt e = 0; e < numEdges; ++e) { - const auto numEdges = m_mesh->GetNumFaceEdges(f); + const auto edgeIndex = m_mesh->m_facesEdges[f][e]; - for (UInt e = 0; e < numEdges; ++e) + if (splitEdge[edgeIndex]) { - const auto edgeIndex = m_mesh->m_facesEdges[f][e]; - if (splitEdge[edgeIndex]) - { - m_faceMask[f] = 1; - } + m_faceMask[f] = 1; } } + } +} - // update edge refinement mask - for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) +void MeshRefinement::UpdateEdgeRefinementMask() +{ + for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) + { + if (m_faceMask[f] != 1) { - if (m_faceMask[f] != 1) - { - continue; - } + continue; + } - const auto numEdges = m_mesh->GetNumFaceEdges(f); + const auto numEdges = m_mesh->GetNumFaceEdges(f); - for (UInt e = 0; e < numEdges; ++e) - { - const auto edgeIndex = m_mesh->m_facesEdges[f][e]; - const auto nextEdgeIndex = NextCircularForwardIndex(e, numEdges); - const auto previousEdgeIndex = NextCircularBackwardIndex(e, numEdges); - const auto split = m_brotherEdges[edgeIndex] != m_mesh->m_facesEdges[f][nextEdgeIndex] && - m_brotherEdges[edgeIndex] != m_mesh->m_facesEdges[f][previousEdgeIndex]; + for (UInt e = 0; e < numEdges; ++e) + { + const auto edgeIndex = m_mesh->m_facesEdges[f][e]; + const auto nextEdgeIndex = NextCircularForwardIndex(e, numEdges); + const auto previousEdgeIndex = NextCircularBackwardIndex(e, numEdges); + const auto split = m_brotherEdges[edgeIndex] != m_mesh->m_facesEdges[f][nextEdgeIndex] && + m_brotherEdges[edgeIndex] != m_mesh->m_facesEdges[f][previousEdgeIndex]; - if (split) - { - m_edgeMask[edgeIndex] = 1; - } + if (split) + { + m_edgeMask[edgeIndex] = 1; } } } } + +void MeshRefinement::SmoothRefinementMasks() +{ + if (m_meshRefinementParameters.directional_refinement == 1) + { + throw AlgorithmError("Directional refinement cannot be used in combination with smoothing. Please set directional refinement to off!"); + } + if (m_meshRefinementParameters.smoothing_iterations == 0) + { + return; + } + + std::vector splitEdge(m_edgeMask.size(), false); + + for (int iter = 0; iter < m_meshRefinementParameters.smoothing_iterations; ++iter) + { + std::fill(splitEdge.begin(), splitEdge.end(), false); + + SmoothEdgeRefinement(splitEdge); + UpdateFaceRefinementMask(splitEdge); + UpdateEdgeRefinementMask(); + } +} From 77cea5c3315ce65fb2120cbc92e7fad5110fdd4b Mon Sep 17 00:00:00 2001 From: BillSenior Date: Wed, 29 Nov 2023 12:34:11 +0100 Subject: [PATCH 14/15] GRIDEDIT-785 Reduced cyclomatic complexity in finding element refinement --- .../include/MeshKernel/MeshRefinement.hpp | 3 + libs/MeshKernel/src/MeshRefinement.cpp | 69 ++++++++++--------- 2 files changed, 41 insertions(+), 31 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp index ab29a250b..47d863d3d 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp @@ -245,6 +245,9 @@ namespace meshkernel /// @brief Smooth the face and edge refinement masks (smooth_jarefine) void SmoothRefinementMasks(); + /// @brief Determine if refinement is required for the face. + bool IsRefinementRequired(const UInt face, const UInt numFaceNodes) const; + /// @brief Computes m_faceMask, if a face must be split later on (split_cells) void ComputeIfFaceShouldBeSplit(); diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index 33ca9625f..10c97fc0e 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -1445,6 +1445,42 @@ void MeshRefinement::ComputeEdgesRefinementMask() } } +bool MeshRefinement::IsRefinementRequired(const UInt face, const UInt numFaceNodes) const +{ + bool isRefinementRequired = false; + + for (UInt n = 0; n < numFaceNodes; n++) + { + const auto edgeIndex = m_mesh->m_facesEdges[face][n]; + + if (m_isHangingEdgeCache[n] && m_edgeMask[edgeIndex] > 0) + { + isRefinementRequired = true; + } + } + + const auto numEdgesToRefine = CountEdgesToRefine(face); + const auto numHangingEdges = CountHangingEdges(); + const auto numHangingNodes = CountHangingNodes(); + + // compute the effective face type + const auto numNodesEffective = numFaceNodes - static_cast(static_cast(numHangingEdges) / 2.0); + if (2 * (numFaceNodes - numNodesEffective) != numHangingEdges) + { + // uneven number of brotherlinks + // TODO: ADD DOT + } + + if (numFaceNodes + numEdgesToRefine > Mesh::m_maximumNumberOfEdgesPerFace || // would result in unsupported cells after refinement + numFaceNodes - numHangingNodes - numEdgesToRefine <= 1 || // faces with only one unrefined edge + numNodesEffective == numEdgesToRefine) // refine all edges + { + isRefinementRequired = true; + } + + return isRefinementRequired; +} + void MeshRefinement::ComputeIfFaceShouldBeSplit() { const UInt maxiter = 1000; @@ -1466,47 +1502,18 @@ void MeshRefinement::ComputeIfFaceShouldBeSplit() continue; } - const auto numEdgesToRefine = CountEdgesToRefine(f); - FindHangingNodes(f); - const auto numHangingEdges = CountHangingEdges(); - const auto numHangingNodes = CountHangingNodes(); - - bool isSplittingRequired = false; // check if the edge has a brother edge and needs to be refined const auto numFaceNodes = m_mesh->GetNumFaceEdges(f); if (numFaceNodes > Mesh::m_maximumNumberOfEdgesPerFace) { + // Should this be an error, or at least a warning message logged? return; } - for (UInt n = 0; n < numFaceNodes; n++) - { - const auto edgeIndex = m_mesh->m_facesEdges[f][n]; - if (m_isHangingEdgeCache[n] && m_edgeMask[edgeIndex] > 0) - { - isSplittingRequired = true; - } - } - - // compute the effective face type - const auto numNodesEffective = numFaceNodes - static_cast(static_cast(numHangingEdges) / 2.0); - if (2 * (numFaceNodes - numNodesEffective) != numHangingEdges) - { - // uneven number of brotherlinks - // TODO: ADD DOT - } - - if (numFaceNodes + numEdgesToRefine > Mesh::m_maximumNumberOfEdgesPerFace || // would result in unsupported cells after refinement - numFaceNodes - numHangingNodes - numEdgesToRefine <= 1 || // faces with only one unrefined edge - numNodesEffective == numEdgesToRefine) // refine all edges - { - isSplittingRequired = true; - } - - if (isSplittingRequired) + if (IsRefinementRequired(f, numFaceNodes)) { if (m_faceMask[f] != -1) { From 45494997748018c5b91bd3e52b2938899976ed4d Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 30 Nov 2023 11:16:48 +0100 Subject: [PATCH 15/15] GRIDEDIT-785 Further refactoring of the mesh refinement class --- .../include/MeshKernel/MeshRefinement.hpp | 48 +- libs/MeshKernel/src/MeshRefinement.cpp | 433 +++++++++++------- 2 files changed, 310 insertions(+), 171 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp index 47d863d3d..a2801b38a 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp @@ -49,12 +49,32 @@ namespace meshkernel std::vector& edgeMask, std::vector& faceMask) = 0; + void ComputeIfFaceShouldBeSplit(const std::vector& brotherEdges, + std::vector& edgeMask, + std::vector& faceMask); + virtual bool IsSampleBased() const { return false; } + // TODO SHould be private + std::vector m_isHangingNodeCache; ///< Cache for maintaining if node is hanging + std::vector m_isHangingEdgeCache; ///< Cache for maintaining if edge is hanging + protected: + UInt CountHangingNodes() const; + UInt CountHangingEdges() const; + UInt CountEdgesToRefine(const std::vector& edgeMask, + UInt face) const; + + // What to do here, it is needed by both the samples based refinement and the mesh-refinement + void FindHangingNodes(const std::vector& brotherEdges, UInt face); + + bool IsRefinementRequired(const std::vector& edgeMask, + const UInt face, + const UInt numFaceNodes) const; + const Mesh2D& GetMesh() const { return m_mesh; @@ -233,18 +253,6 @@ namespace meshkernel /// @brief Connect the hanging nodes with triangles (connect_hanging_nodes) void ConnectHangingNodes(); - /// @brief Update edge refinement indicator to ensure smoother refinement transition between elements. - void SmoothEdgeRefinement(std::vector& splitEdge); - - /// @brief Using the latest edge refinement indicator, update face refinement mask. - void UpdateFaceRefinementMask(const std::vector& splitEdge); - - /// @brief Update the edge refinement mask - void UpdateEdgeRefinementMask(); - - /// @brief Smooth the face and edge refinement masks (smooth_jarefine) - void SmoothRefinementMasks(); - /// @brief Determine if refinement is required for the face. bool IsRefinementRequired(const UInt face, const UInt numFaceNodes) const; @@ -343,17 +351,21 @@ namespace meshkernel std::shared_ptr m_interpolant; MeshRefinementParameters m_meshRefinementParameters; ///< The mesh refinement parameters - std::vector m_isHangingNodeCache; ///< Cache for maintaining if node is hanging - std::vector m_isHangingEdgeCache; ///< Cache for maintaining if edge is hanging - private: + void SmoothEdgeRefinement(const std::vector& brotherEdges, + const std::vector& faceMask, + std::vector& splitEdge); + + void UpdateFaceRefinementMask(const std::vector& splitEdge, std::vector& faceMask); + + void UpdateEdgeRefinementMask(const std::vector& brotherEdges, + const std::vector& faceMask, + std::vector& edgeMask); + void SmoothRefinementMasks(const std::vector& brotherEdges, std::vector& edgeMask, std::vector& faceMask); - // What to do here, it is needed by both the samples based refinement and the mesh-refinement - void FindHangingNodes(const std::vector& brotherEdges, UInt face); - virtual UInt ComputeForFace(const UInt faceId, const std::vector& edgeIsBelowMinimumSize, std::vector& refineEdgeCache) = 0; diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index 10c97fc0e..b72406bd1 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -32,6 +32,7 @@ #include #include +using meshkernel::ComputeRefinementMask; using meshkernel::EdgeSizeBasedRefinement; using meshkernel::Mesh2D; using meshkernel::MeshRefinement; @@ -40,75 +41,154 @@ using meshkernel::RidgeDetectionRefinement; using meshkernel::SamplesBasedRefinement; using meshkernel::WaveCourantRefinement; -void EdgeSizeBasedRefinement::Compute(const std::vector& edgeIsBelowMinimumSize, - const std::vector& brotherEdges [[maybe_unused]], - std::vector& edgeMask, - std::vector& faceMask) +meshkernel::UInt ComputeRefinementMask::CountHangingNodes() const { - for (UInt f = 0; f < GetMesh().GetNumFaces(); ++f) + meshkernel::UInt result = 0; + for (const auto& v : m_isHangingNodeCache) { - for (UInt n = 0; n < GetMesh().GetNumFaceEdges(f); ++n) + if (v) { - const auto e = GetMesh().m_facesEdges[f][n]; - if (!edgeIsBelowMinimumSize[e]) - { - edgeMask[e] = -1; - faceMask[f] = 1; - } + result++; } } + return result; } -void SamplesBasedRefinement::Compute(const std::vector& edgeIsBelowMinimumSize, - const std::vector& brotherEdges, - std::vector& edgeMask, - std::vector& faceMask) +meshkernel::UInt ComputeRefinementMask::CountHangingEdges() const { + meshkernel::UInt result = 0; + for (const auto& v : m_isHangingEdgeCache) + { + if (v) + { + result++; + } + } + return result; +} - m_refineEdgeCache.resize(Mesh::m_maximumNumberOfEdgesPerFace, 0); +meshkernel::UInt ComputeRefinementMask::CountEdgesToRefine(const std::vector& edgeMask, + UInt face) const +{ + const auto numFaceNodes = GetMesh().GetNumFaceEdges(face); - // Compute all interpolated values - m_interpolant->Compute(); + UInt result = 0; - for (UInt face = 0; face < GetMesh().GetNumFaces(); face++) + for (UInt n = 0; n < numFaceNodes; n++) { - FindHangingNodes(brotherEdges, face); + const auto edgeIndex = GetMesh().m_facesEdges[face][n]; - if (IsEqual(m_interpolant->GetFaceResult(face), constants::missing::doubleValue)) + if (edgeMask[edgeIndex] != 0) { - continue; + result += 1; } + } + return result; +} - std::ranges::fill(m_refineEdgeCache, 0); - UInt numberToBeRefined = ComputeForFace(face, edgeIsBelowMinimumSize, m_refineEdgeCache); +bool ComputeRefinementMask::IsRefinementRequired(const std::vector& edgeMask, + const UInt face, + const UInt numFaceNodes) const +{ + bool isRefinementRequired = false; - if (numberToBeRefined > 1) + for (UInt n = 0; n < numFaceNodes; n++) + { + const auto edgeIndex = GetMesh().m_facesEdges[face][n]; + + if (m_isHangingEdgeCache[n] && edgeMask[edgeIndex] > 0) { - faceMask[face] = 1; + isRefinementRequired = true; + } + } - for (size_t n = 0; n < GetMesh().GetNumFaceEdges(face); n++) + const auto numEdgesToRefine = CountEdgesToRefine(edgeMask, face); + const auto numHangingEdges = CountHangingEdges(); + const auto numHangingNodes = CountHangingNodes(); + + // compute the effective face type + const auto numNodesEffective = numFaceNodes - static_cast(static_cast(numHangingEdges) / 2.0); + if (2 * (numFaceNodes - numNodesEffective) != numHangingEdges) + { + // uneven number of brotherlinks + // TODO: ADD DOT + } + + if (numFaceNodes + numEdgesToRefine > Mesh::m_maximumNumberOfEdgesPerFace || // would result in unsupported cells after refinement + numFaceNodes - numHangingNodes - numEdgesToRefine <= 1 || // faces with only one unrefined edge + numNodesEffective == numEdgesToRefine) // refine all edges + { + isRefinementRequired = true; + } + + return isRefinementRequired; +} + +void ComputeRefinementMask::ComputeIfFaceShouldBeSplit(const std::vector& brotherEdges, + std::vector& edgeMask, + std::vector& faceMask) +{ + const UInt maxiter = 1000; + UInt num = 1; + UInt iter = 0; + + while (num != 0) + { + iter++; + if (iter > maxiter) + { + break; + } + + num = 0; + for (UInt f = 0; f < GetMesh().GetNumFaces(); f++) + { + if (faceMask[f] != 0 && faceMask[f] != -1) { - if (m_refineEdgeCache[n] == 1) + continue; + } + + FindHangingNodes(brotherEdges, f); + + // check if the edge has a brother edge and needs to be refined + const auto numFaceNodes = GetMesh().GetNumFaceEdges(f); + + if (numFaceNodes > Mesh::m_maximumNumberOfEdgesPerFace) + { + // Should this be an error, or at least a warning message logged? + return; + } + + if (IsRefinementRequired(edgeMask, f, numFaceNodes)) + { + if (faceMask[f] != -1) { - const auto edgeIndex = GetMesh().m_facesEdges[face][n]; - if (edgeIndex != constants::missing::uintValue) + faceMask[f] = 2; + } + else + { + faceMask[f] = -2; + } + + for (UInt n = 0; n < numFaceNodes; n++) + { + const auto edgeIndex = GetMesh().m_facesEdges[f][n]; + if (!m_isHangingEdgeCache[n] && edgeMask[edgeIndex] == 0) { edgeMask[edgeIndex] = 1; + num++; + } + if (iter == maxiter) + { + // TODO: ADD DOT/MESSAGES } } } } } - - for (auto& edge : edgeMask) - { - edge = -edge; - } - - SmoothRefinementMasks(brotherEdges, edgeMask, faceMask); } -void SamplesBasedRefinement::FindHangingNodes(const std::vector& brotherEdges, UInt face) +void ComputeRefinementMask::FindHangingNodes(const std::vector& brotherEdges, UInt face) { const auto numFaceNodes = GetMesh().GetNumFaceEdges(face); @@ -171,6 +251,151 @@ void SamplesBasedRefinement::FindHangingNodes(const std::vector& brotherEd } } +void EdgeSizeBasedRefinement::Compute(const std::vector& edgeIsBelowMinimumSize, + const std::vector& brotherEdges [[maybe_unused]], + std::vector& edgeMask, + std::vector& faceMask) +{ + for (UInt f = 0; f < GetMesh().GetNumFaces(); ++f) + { + for (UInt n = 0; n < GetMesh().GetNumFaceEdges(f); ++n) + { + const auto e = GetMesh().m_facesEdges[f][n]; + if (!edgeIsBelowMinimumSize[e]) + { + edgeMask[e] = -1; + faceMask[f] = 1; + } + } + } +} + +void SamplesBasedRefinement::Compute(const std::vector& edgeIsBelowMinimumSize, + const std::vector& brotherEdges, + std::vector& edgeMask, + std::vector& faceMask) +{ + + m_refineEdgeCache.resize(Mesh::m_maximumNumberOfEdgesPerFace, 0); + + // Compute all interpolated values + m_interpolant->Compute(); + + for (UInt face = 0; face < GetMesh().GetNumFaces(); face++) + { + FindHangingNodes(brotherEdges, face); + + if (IsEqual(m_interpolant->GetFaceResult(face), constants::missing::doubleValue)) + { + continue; + } + + std::ranges::fill(m_refineEdgeCache, 0); + UInt numberToBeRefined = ComputeForFace(face, edgeIsBelowMinimumSize, m_refineEdgeCache); + + if (numberToBeRefined > 1) + { + faceMask[face] = 1; + + for (size_t n = 0; n < GetMesh().GetNumFaceEdges(face); n++) + { + if (m_refineEdgeCache[n] == 1) + { + const auto edgeIndex = GetMesh().m_facesEdges[face][n]; + if (edgeIndex != constants::missing::uintValue) + { + edgeMask[edgeIndex] = 1; + } + } + } + } + } + + for (auto& edge : edgeMask) + { + edge = -edge; + } + + SmoothRefinementMasks(brotherEdges, edgeMask, faceMask); +} + +void SamplesBasedRefinement::SmoothEdgeRefinement(const std::vector& brotherEdges, + const std::vector& faceMask, + std::vector& splitEdge) +{ + + for (UInt f = 0; f < GetMesh().GetNumFaces(); ++f) + { + if (faceMask[f] != 1) + { + continue; + } + + const auto numEdges = GetMesh().GetNumFaceEdges(f); + + for (UInt e = 0; e < numEdges; ++e) + { + const auto edgeIndex = GetMesh().m_facesEdges[f][e]; + const auto nextEdgeIndex = NextCircularForwardIndex(e, numEdges); + const auto previousEdgeIndex = NextCircularBackwardIndex(e, numEdges); + const auto split = brotherEdges[edgeIndex] != GetMesh().m_facesEdges[f][nextEdgeIndex] && + brotherEdges[edgeIndex] != GetMesh().m_facesEdges[f][previousEdgeIndex]; + + if (split) + { + splitEdge[edgeIndex] = true; + } + } + } +} + +void SamplesBasedRefinement::UpdateFaceRefinementMask(const std::vector& splitEdge, std::vector& faceMask) +{ + for (UInt f = 0; f < GetMesh().GetNumFaces(); ++f) + { + const auto numEdges = GetMesh().GetNumFaceEdges(f); + + for (UInt e = 0; e < numEdges; ++e) + { + const auto edgeIndex = GetMesh().m_facesEdges[f][e]; + + if (splitEdge[edgeIndex]) + { + faceMask[f] = 1; + } + } + } +} + +void SamplesBasedRefinement::UpdateEdgeRefinementMask(const std::vector& brotherEdges, + const std::vector& faceMask, + std::vector& edgeMask) +{ + for (UInt f = 0; f < GetMesh().GetNumFaces(); ++f) + { + if (faceMask[f] != 1) + { + continue; + } + + const auto numEdges = GetMesh().GetNumFaceEdges(f); + + for (UInt e = 0; e < numEdges; ++e) + { + const auto edgeIndex = GetMesh().m_facesEdges[f][e]; + const auto nextEdgeIndex = NextCircularForwardIndex(e, numEdges); + const auto previousEdgeIndex = NextCircularBackwardIndex(e, numEdges); + const auto split = brotherEdges[edgeIndex] != GetMesh().m_facesEdges[f][nextEdgeIndex] && + brotherEdges[edgeIndex] != GetMesh().m_facesEdges[f][previousEdgeIndex]; + + if (split) + { + edgeMask[edgeIndex] = 1; + } + } + } +} + void SamplesBasedRefinement::SmoothRefinementMasks(const std::vector& brotherEdges, std::vector& edgeMask, std::vector& faceMask) @@ -185,7 +410,16 @@ void SamplesBasedRefinement::SmoothRefinementMasks(const std::vector& brot } std::vector splitEdge(edgeMask.size(), false); + for (int iter = 0; iter < m_meshRefinementParameters.smoothing_iterations; ++iter) + { + std::fill(splitEdge.begin(), splitEdge.end(), false); + + SmoothEdgeRefinement(brotherEdges, faceMask, splitEdge); + UpdateFaceRefinementMask(splitEdge, faceMask); + UpdateEdgeRefinementMask(brotherEdges, faceMask, edgeMask); + } +#if 0 for (int iter = 0; iter < m_meshRefinementParameters.smoothing_iterations; ++iter) { std::fill(splitEdge.begin(), splitEdge.end(), false); @@ -254,6 +488,7 @@ void SamplesBasedRefinement::SmoothRefinementMasks(const std::vector& brot } } } +#endif } meshkernel::UInt RefinementLevelsRefinement::ComputeForFace(const UInt face, @@ -551,6 +786,7 @@ void MeshRefinement::ReviewRefinement(const int level) } } } + if (level > 0) { // if one face node is not in polygon disable refinement @@ -619,29 +855,16 @@ void MeshRefinement::Compute() ReviewRefinement(level); ComputeEdgesRefinementMask(); - ComputeIfFaceShouldBeSplit(); + m_refinementMasks->ComputeIfFaceShouldBeSplit(m_brotherEdges, m_edgeMask, m_faceMask); - UInt numFacesToRefine = 0; - for (UInt f = 0; f < m_mesh->GetNumFaces(); f++) - { - if (m_faceMask[f] != 0) - { - numFacesToRefine++; - } - } - if (numFacesToRefine == 0) + auto notZero = [](UInt value) + { return value != 0; }; + + if (std::count_if(m_faceMask.begin(), m_faceMask.end(), notZero) == 0) { break; } - // auto notZero = [](UInt value) - // { return value != 0; }; - - // if (std::count_if(m_faceMask.begin(), m_faceMask.end(), notZero) == 0) - // { - // break; - // } - // spit the edges RefineFacesBySplittingEdges(); @@ -1600,99 +1823,3 @@ void MeshRefinement::FindBrotherEdges() } } } - -void MeshRefinement::SmoothEdgeRefinement(std::vector& splitEdge) -{ - - for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) - { - if (m_faceMask[f] != 1) - { - continue; - } - - const auto numEdges = m_mesh->GetNumFaceEdges(f); - - for (UInt e = 0; e < numEdges; ++e) - { - const auto edgeIndex = m_mesh->m_facesEdges[f][e]; - const auto nextEdgeIndex = NextCircularForwardIndex(e, numEdges); - const auto previousEdgeIndex = NextCircularBackwardIndex(e, numEdges); - const auto split = m_brotherEdges[edgeIndex] != m_mesh->m_facesEdges[f][nextEdgeIndex] && - m_brotherEdges[edgeIndex] != m_mesh->m_facesEdges[f][previousEdgeIndex]; - - if (split) - { - splitEdge[edgeIndex] = true; - } - } - } -} - -void MeshRefinement::UpdateFaceRefinementMask(const std::vector& splitEdge) -{ - for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) - { - const auto numEdges = m_mesh->GetNumFaceEdges(f); - - for (UInt e = 0; e < numEdges; ++e) - { - const auto edgeIndex = m_mesh->m_facesEdges[f][e]; - - if (splitEdge[edgeIndex]) - { - m_faceMask[f] = 1; - } - } - } -} - -void MeshRefinement::UpdateEdgeRefinementMask() -{ - for (UInt f = 0; f < m_mesh->GetNumFaces(); ++f) - { - if (m_faceMask[f] != 1) - { - continue; - } - - const auto numEdges = m_mesh->GetNumFaceEdges(f); - - for (UInt e = 0; e < numEdges; ++e) - { - const auto edgeIndex = m_mesh->m_facesEdges[f][e]; - const auto nextEdgeIndex = NextCircularForwardIndex(e, numEdges); - const auto previousEdgeIndex = NextCircularBackwardIndex(e, numEdges); - const auto split = m_brotherEdges[edgeIndex] != m_mesh->m_facesEdges[f][nextEdgeIndex] && - m_brotherEdges[edgeIndex] != m_mesh->m_facesEdges[f][previousEdgeIndex]; - - if (split) - { - m_edgeMask[edgeIndex] = 1; - } - } - } -} - -void MeshRefinement::SmoothRefinementMasks() -{ - if (m_meshRefinementParameters.directional_refinement == 1) - { - throw AlgorithmError("Directional refinement cannot be used in combination with smoothing. Please set directional refinement to off!"); - } - if (m_meshRefinementParameters.smoothing_iterations == 0) - { - return; - } - - std::vector splitEdge(m_edgeMask.size(), false); - - for (int iter = 0; iter < m_meshRefinementParameters.smoothing_iterations; ++iter) - { - std::fill(splitEdge.begin(), splitEdge.end(), false); - - SmoothEdgeRefinement(splitEdge); - UpdateFaceRefinementMask(splitEdge); - UpdateEdgeRefinementMask(); - } -}