diff --git a/libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp b/libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp index 02051294e..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,9 +154,10 @@ 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 - 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..ad9aaad7b 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. @@ -44,5 +47,14 @@ namespace meshkernel::averaging /// @brief Calculates the average value based on the values added. /// @return The calculated average [[nodiscard]] virtual double Calculate() const = 0; + + /// @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, + const std::vector& sampleValues) const = 0; }; } // namespace meshkernel::averaging 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..f17b936c3 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp @@ -39,12 +39,23 @@ 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; + /// @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; + private: /// @brief The result used to calculate the final value in Calculate. double m_result = constants::missing::doubleValue; @@ -53,7 +64,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..2463f9e80 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/InverseWeightedAveragingStrategy.hpp @@ -37,16 +37,26 @@ 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; + /// @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; + private: /// @brief The current result used in Calculate to calculate the final value. double m_result = 0.0; @@ -55,10 +65,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..d8f12cad8 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MaxAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MaxAveragingStrategy.hpp @@ -36,10 +36,22 @@ 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; + /// @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; + 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 921a0f461..92feb07e1 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAbsAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAbsAveragingStrategy.hpp @@ -36,10 +36,22 @@ 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; + /// @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; + 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 e2920d78d..d40353a33 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/MinAveragingStrategy.hpp @@ -35,10 +35,22 @@ 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; + /// @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; + 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 51df35b86..e0c954d4a 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp @@ -13,12 +13,24 @@ 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; + /// @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; + 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 38604049a..a2801b38a 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp @@ -37,6 +37,53 @@ namespace meshkernel // Forward declarations class Mesh2D; + class ComputeRefinementMask + { + public: + ComputeRefinementMask(const Mesh2D& mesh) : m_mesh(mesh) {} + + virtual ~ComputeRefinementMask() = default; + + virtual void Compute(const std::vector& edgeIsBelowMinimumSize, + const std::vector& brotherEdges, + 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; + } + + private: + const Mesh2D& m_mesh; + }; + /// @brief A class used to refine a Mesh2D instance /// /// Mesh refinement operates on Mesh2D and is based on @@ -73,6 +120,7 @@ namespace meshkernel /// existing Mesh2D instance. class MeshRefinement { + public: /// @brief Enumerator describing the face location types enum class FaceLocation { @@ -89,7 +137,11 @@ namespace meshkernel RidgeDetection = 3 }; - public: + /// @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); + /// @brief The constructor for refining based on samples /// @param[in] mesh The mesh to be refined /// @param[in] interpolant The averaging interpolation to use @@ -132,6 +184,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(); @@ -192,8 +253,8 @@ namespace meshkernel /// @brief Connect the hanging nodes with triangles (connect_hanging_nodes) void ConnectHangingNodes(); - /// @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(); @@ -219,13 +280,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 @@ -245,5 +310,117 @@ 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 ComputeRefinementMask + { + public: + EdgeSizeBasedRefinement(const Mesh2D& mesh) : ComputeRefinementMask(mesh) {} + + void Compute(const std::vector& edgeIsBelowMinimumSize, + const std::vector& brotherEdges, + std::vector& edgeMask, + std::vector& faceMask) override; + + private: + }; + + class SamplesBasedRefinement : public ComputeRefinementMask + { + public: + SamplesBasedRefinement(const Mesh2D& mesh, + std::shared_ptr interpolant, + 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, + 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 + + 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); + + 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(const 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(const 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(const 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/AveragingInterpolation.cpp b/libs/MeshKernel/src/AveragingInterpolation.cpp index 524b5e88a..a2c109804 100644 --- a/libs/MeshKernel/src/AveragingInterpolation.cpp +++ b/libs/MeshKernel/src/AveragingInterpolation.cpp @@ -49,8 +49,11 @@ 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)) { + m_interpolationPointCache.reserve (DefaultMaximumCacheSize); + m_interpolationSampleCache.reserve (DefaultMaximumCacheSize); } void AveragingInterpolation::Compute() @@ -65,16 +68,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(); @@ -92,7 +89,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) { @@ -110,11 +106,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(); @@ -133,9 +128,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; } } @@ -199,8 +194,8 @@ 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_interpolationPointCache.resize (0); + m_interpolationSampleCache.resize (0); for (UInt i = 0; i < m_samplesRtree.GetQueryResultSize(); ++i) { @@ -213,13 +208,15 @@ 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_interpolationPointCache.emplace_back (samplePoint); + m_interpolationSampleCache.emplace_back (sampleValue); } } - return 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/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..0dce1cb0f 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) { @@ -48,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 81724c98b..d2c1adc7d 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)); @@ -21,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/AveragingStrategies/MaxAveragingStrategy.cpp b/libs/MeshKernel/src/AveragingStrategies/MaxAveragingStrategy.cpp index 9f058eb1f..aa82fe891 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); @@ -38,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 7b1d62b7a..96348bcaf 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)); @@ -38,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 d12458492..daf93538d 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); @@ -38,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 f66a0cd24..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 @@ -32,6 +34,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; @@ -42,4 +50,13 @@ 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 = std::accumulate (sampleValues.begin (), sampleValues.end (), 0.0); + return sampleValues.size () >= m_minNumPoints ? result / static_cast(sampleValues.size ()) : constants::missing::doubleValue; + } + } // namespace meshkernel::averaging diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index f35744696..b72406bd1 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -32,8 +32,695 @@ #include #include +using meshkernel::ComputeRefinementMask; +using meshkernel::EdgeSizeBasedRefinement; using meshkernel::Mesh2D; using meshkernel::MeshRefinement; +using meshkernel::RefinementLevelsRefinement; +using meshkernel::RidgeDetectionRefinement; +using meshkernel::SamplesBasedRefinement; +using meshkernel::WaveCourantRefinement; + +meshkernel::UInt ComputeRefinementMask::CountHangingNodes() const +{ + meshkernel::UInt result = 0; + for (const auto& v : m_isHangingNodeCache) + { + if (v) + { + result++; + } + } + return result; +} + +meshkernel::UInt ComputeRefinementMask::CountHangingEdges() const +{ + meshkernel::UInt result = 0; + for (const auto& v : m_isHangingEdgeCache) + { + if (v) + { + result++; + } + } + return result; +} + +meshkernel::UInt ComputeRefinementMask::CountEdgesToRefine(const std::vector& edgeMask, + UInt face) const +{ + const auto numFaceNodes = GetMesh().GetNumFaceEdges(face); + + UInt result = 0; + + for (UInt n = 0; n < numFaceNodes; n++) + { + const auto edgeIndex = GetMesh().m_facesEdges[face][n]; + + if (edgeMask[edgeIndex] != 0) + { + result += 1; + } + } + return result; +} + +bool ComputeRefinementMask::IsRefinementRequired(const std::vector& edgeMask, + const UInt face, + const UInt numFaceNodes) const +{ + bool isRefinementRequired = false; + + for (UInt n = 0; n < numFaceNodes; n++) + { + const auto edgeIndex = GetMesh().m_facesEdges[face][n]; + + if (m_isHangingEdgeCache[n] && edgeMask[edgeIndex] > 0) + { + isRefinementRequired = true; + } + } + + 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) + { + 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) + { + 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 + } + } + } + } + } +} + +void ComputeRefinementMask::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 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) +{ + 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); + + 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); + + 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; + } + } + } + } +#endif +} + +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) +{ + 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); + } +} + +std::unique_ptr MeshRefinement::CreateRefinementMaskCalculator(const RefinementType refinementType, + const Mesh2D& mesh, + std::shared_ptr interpolant, + const MeshRefinementParameters& meshRefinementParameters, + const bool useNodalRefinement) +{ + std::unique_ptr refinementMaskCalculator; + + if (interpolant != nullptr) + { + switch (refinementType) + { + case RefinementType::WaveCourant: + refinementMaskCalculator = std::make_unique(mesh, interpolant, meshRefinementParameters, useNodalRefinement); + break; + case RefinementType::RefinementLevels: + refinementMaskCalculator = std::make_unique(mesh, interpolant, meshRefinementParameters); + break; + case RefinementType::RidgeDetection: + refinementMaskCalculator = std::make_unique(mesh, interpolant, meshRefinementParameters); + break; + default: + throw ConstraintError("Invalid mesh refinement type when creating with interpolant: refinement type {}", meshRefinementParameters.refinement_type); + } + } + else + { + refinementMaskCalculator = std::make_unique(mesh); + } + + return refinementMaskCalculator; +} MeshRefinement::MeshRefinement(std::shared_ptr mesh, std::shared_ptr interpolant, @@ -43,7 +730,8 @@ 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); + m_refinementMasks = CreateRefinementMaskCalculator(m_refinementType, *mesh, interpolant, meshRefinementParameters, false); } MeshRefinement::MeshRefinement(std::shared_ptr mesh, @@ -56,7 +744,8 @@ 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); + m_refinementMasks = CreateRefinementMaskCalculator(m_refinementType, *mesh, interpolant, meshRefinementParameters, useNodalRefinement); } MeshRefinement::MeshRefinement(std::shared_ptr mesh, @@ -67,6 +756,53 @@ MeshRefinement::MeshRefinement(std::shared_ptr mesh, { CheckMeshRefinementParameters(meshRefinementParameters); m_meshRefinementParameters = meshRefinementParameters; + m_refinementMasks = std::make_unique(*mesh); +} + +void MeshRefinement::ReviewRefinement(const int level) +{ + + 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) + { + 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; + } + } + } + } } void MeshRefinement::Compute() @@ -89,7 +825,7 @@ void MeshRefinement::Compute() } // select the nodes to refine - auto const isRefinementBasedOnSamples = 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); @@ -105,79 +841,29 @@ 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 (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(); - // 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); - 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; - } - } - } - } + m_refinementMasks->Compute(m_isEdgeBelowMinSizeAfterRefinement, m_brotherEdges, m_edgeMask, m_faceMask); + ReviewRefinement(level); ComputeEdgesRefinementMask(); + m_refinementMasks->ComputeIfFaceShouldBeSplit(m_brotherEdges, m_edgeMask, m_faceMask); - ComputeIfFaceShouldBeSplit(); + auto notZero = [](UInt value) + { return value != 0; }; - UInt numFacesToRefine = 0; - for (UInt f = 0; f < m_mesh->GetNumFaces(); f++) - { - if (m_faceMask[f] != 0) - { - numFacesToRefine++; - } - } - if (numFacesToRefine == 0) + if (std::count_if(m_faceMask.begin(), m_faceMask.end(), notZero) == 0) { break; } - numFacesAfterRefinement = numFacesAfterRefinement * 4; // spit the edges RefineFacesBySplittingEdges(); @@ -200,24 +886,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() { @@ -698,385 +1366,150 @@ void MeshRefinement::RefineFacesBySplittingEdges() { if (m_edgeMask[e] > 0) { - const auto newEdgeIndex = m_mesh->ConnectNodes(m_edgeMask[e], m_mesh->m_edges[e].second); - m_mesh->m_edges[e].second = m_edgeMask[e]; - m_brotherEdges.resize(m_mesh->GetNumEdges()); - m_brotherEdges[newEdgeIndex] = e; - m_brotherEdges[e] = newEdgeIndex; - } - } -} - -void MeshRefinement::ComputeNodeMaskAtPolygonPerimeter() -{ - for (UInt f = 0; f < m_mesh->GetNumFaces(); f++) - { - bool crossing = false; - const auto numnodes = m_mesh->GetNumFaceEdges(f); - for (UInt n = 0; n < numnodes; n++) - { - const auto nodeIndex = m_mesh->m_facesNodes[f][n]; - if (m_nodeMask[nodeIndex] == 0) - { - crossing = true; - break; - } - } - - if (crossing) - { - m_faceMask[f] = 0; - for (UInt n = 0; n < numnodes; n++) - { - const auto nodeIndex = m_mesh->m_facesNodes[f][n]; - if (m_nodeMask[nodeIndex] == 1) - { - m_nodeMask[nodeIndex] = -2; - } - } - } - } -} - -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); - - 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 = m_mesh->m_facesEdges[face][n]; - - // check if the parent edge is in the cell - if (m_brotherEdges[edgeIndex] != constants::missing::uintValue) - { - const auto e = NextCircularBackwardIndex(n, numFaceNodes); - const auto ee = NextCircularForwardIndex(n, numFaceNodes); - const auto firstEdgeIndex = m_mesh->m_facesEdges[face][e]; - const auto secondEdgeIndex = m_mesh->m_facesEdges[face][ee]; - - UInt commonNode = constants::missing::uintValue; - if (m_brotherEdges[edgeIndex] == firstEdgeIndex) - { - commonNode = m_mesh->FindCommonNode(edgeIndex, firstEdgeIndex); - if (commonNode == constants::missing::uintValue) - { - throw AlgorithmError("Could not find common node."); - } - } - else if (m_brotherEdges[edgeIndex] == secondEdgeIndex) - { - commonNode = m_mesh->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 (m_mesh->m_facesNodes[face][kknod] == commonNode && !m_isHangingNodeCache[kknod]) - { - m_isHangingNodeCache[kknod] = true; - break; - } - } - } - } - } -} - -meshkernel::UInt MeshRefinement::CountHangingNodes() const -{ - meshkernel::UInt result = 0; - for (const auto& v : m_isHangingNodeCache) - { - if (v) - { - result++; - } - } - return result; -} - -meshkernel::UInt MeshRefinement::CountHangingEdges() const -{ - meshkernel::UInt result = 0; - for (const auto& v : m_isHangingEdgeCache) - { - if (v) - { - result++; - } - } - return result; -} - -meshkernel::UInt MeshRefinement::CountEdgesToRefine(UInt face) const -{ - const auto numFaceNodes = m_mesh->GetNumFaceEdges(face); - - UInt result = 0; - - for (UInt n = 0; n < numFaceNodes; n++) - { - const auto edgeIndex = m_mesh->m_facesEdges[face][n]; - if (m_edgeMask[edgeIndex] != 0) - { - result += 1; - } - } - 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; + const auto newEdgeIndex = m_mesh->ConnectNodes(m_edgeMask[e], m_mesh->m_edges[e].second); + m_mesh->m_edges[e].second = m_edgeMask[e]; + m_brotherEdges.resize(m_mesh->GetNumEdges()); + m_brotherEdges[newEdgeIndex] = e; + m_brotherEdges[e] = newEdgeIndex; } } +} - if (numberOfEdgesToRefine > 0) +void MeshRefinement::ComputeNodeMaskAtPolygonPerimeter() +{ + for (UInt f = 0; f < m_mesh->GetNumFaces(); f++) { - numberOfEdgesToRefine = 0; - for (UInt i = 0; i < m_mesh->GetNumFaceEdges(face); i++) + bool crossing = false; + const auto numnodes = m_mesh->GetNumFaceEdges(f); + for (UInt n = 0; n < numnodes; n++) { - if (edgeToRefine[i] == 1 || m_isHangingNodeCache[i]) + const auto nodeIndex = m_mesh->m_facesNodes[f][n]; + if (m_nodeMask[nodeIndex] == 0) { - numberOfEdgesToRefine++; + crossing = true; + break; } } - } - if (m_meshRefinementParameters.directional_refinement == 0) - { - if (numberOfEdgesToRefine == m_mesh->GetNumFaceEdges(face)) + if (crossing) { - for (UInt i = 0; i < m_mesh->GetNumFaceEdges(face); i++) + m_faceMask[f] = 0; + for (UInt n = 0; n < numnodes; n++) { - if (!m_isHangingNodeCache[i]) + const auto nodeIndex = m_mesh->m_facesNodes[f][n]; + if (m_nodeMask[nodeIndex] == 1) { - edgeToRefine[i] = 1; + m_nodeMask[nodeIndex] = -2; } } } - else - { - numberOfEdgesToRefine = 0; - } } } -void MeshRefinement::ComputeRefinementMasksForRidgeDetection(UInt face, - size_t& numberOfEdgesToRefine, - std::vector& edgeToRefine) const +void MeshRefinement::FindHangingNodes(UInt face) { + const auto numFaceNodes = m_mesh->GetNumFaceEdges(face); - double maxEdgeLength = 0.0; - UInt numEdges = m_mesh->GetNumFaceEdges(face); - - for (size_t i = 0; i < numEdges; ++i) + if (numFaceNodes > Mesh::m_maximumNumberOfEdgesPerNode) { - - 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); + throw AlgorithmError("The number of face nodes is greater than the maximum number of edges per node."); } - 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); + 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); - if (maxEdgeLength > elementSizeWanted && maxEdgeLength > 2.0 * hmin && absInterpolatedFaceValue > thresholdMin) + auto kknod = numFaceNodes; + for (UInt n = 0; n < numFaceNodes; n++) { - numberOfEdgesToRefine += numEdges; + const auto edgeIndex = m_mesh->m_facesEdges[face][n]; - for (UInt i = 0; i < numEdges; i++) + // check if the parent edge is in the cell + if (m_brotherEdges[edgeIndex] != constants::missing::uintValue) { - 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"); - } + const auto e = NextCircularBackwardIndex(n, numFaceNodes); + const auto ee = NextCircularForwardIndex(n, numFaceNodes); + const auto firstEdgeIndex = m_mesh->m_facesEdges[face][e]; + const auto secondEdgeIndex = m_mesh->m_facesEdges[face][ee]; - // Compute face and edge masks - if (numEdgesToBeRefined > 1) - { - m_faceMask[face] = 1; + UInt commonNode = constants::missing::uintValue; + if (m_brotherEdges[edgeIndex] == firstEdgeIndex) + { + commonNode = m_mesh->FindCommonNode(edgeIndex, firstEdgeIndex); + if (commonNode == constants::missing::uintValue) + { + throw AlgorithmError("Could not find common node."); + } + } + else if (m_brotherEdges[edgeIndex] == secondEdgeIndex) + { + commonNode = m_mesh->FindCommonNode(edgeIndex, secondEdgeIndex); + if (commonNode == constants::missing::uintValue) + { + throw AlgorithmError("Could not find common node."); + } + } - for (size_t n = 0; n < m_mesh->GetNumFaceEdges(face); n++) - { - if (m_refineEdgeCache[n] == 1) + if (commonNode != constants::missing::uintValue) { - const auto edgeIndex = m_mesh->m_facesEdges[face][n]; - if (edgeIndex != constants::missing::uintValue) + m_isHangingEdgeCache[n] = true; + for (UInt nn = 0; nn < numFaceNodes; nn++) { - m_edgeMask[edgeIndex] = 1; + kknod = NextCircularForwardIndex(kknod, numFaceNodes); + + if (m_mesh->m_facesNodes[face][kknod] == commonNode && !m_isHangingNodeCache[kknod]) + { + m_isHangingNodeCache[kknod] = true; + break; + } } } } } } -void MeshRefinement::ComputeFaceLocationTypes() +meshkernel::UInt MeshRefinement::CountHangingNodes() const { - m_faceLocationType.resize(m_mesh->GetNumFaces()); - std::ranges::fill(m_faceLocationType, FaceLocation::Water); - for (UInt face = 0; face < m_mesh->GetNumFaces(); face++) + meshkernel::UInt result = 0; + for (const auto& v : m_isHangingNodeCache) { - double maxVal = -std::numeric_limits::max(); - double minVal = std::numeric_limits::max(); - for (UInt e = 0; e < m_mesh->GetNumFaceEdges(face); ++e) + if (v) { - 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); + result++; } + } + return result; +} - if (minVal > 0.0) +meshkernel::UInt MeshRefinement::CountHangingEdges() const +{ + meshkernel::UInt result = 0; + for (const auto& v : m_isHangingEdgeCache) + { + if (v) { - m_faceLocationType[face] = FaceLocation::Land; + result++; } - if (maxVal >= 0.0 && (!IsEqual(minVal, constants::missing::doubleValue) && minVal < 0.0)) + } + return result; +} + +meshkernel::UInt MeshRefinement::CountEdgesToRefine(UInt face) const +{ + const auto numFaceNodes = m_mesh->GetNumFaceEdges(face); + + UInt result = 0; + + for (UInt n = 0; n < numFaceNodes; n++) + { + const auto edgeIndex = m_mesh->m_facesEdges[face][n]; + if (m_edgeMask[edgeIndex] != 0) { - m_faceLocationType[face] = FaceLocation::LandWater; + result += 1; } } + return result; } void MeshRefinement::ComputeEdgeBelowMinSizeAfterRefinement() @@ -1090,14 +1523,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; @@ -1243,6 +1668,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; @@ -1264,47 +1725,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) { @@ -1391,86 +1823,3 @@ void MeshRefinement::FindBrotherEdges() } } } - -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); - - 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; - } - } - } - - // update face refinement mask - 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; - } - } - } - - // update edge refinement mask - 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; - } - } - } - } -} 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_) {