Skip to content

Commit

Permalink
Implement ridge mesh refinement (#235 | GRIDEDIT-502)
Browse files Browse the repository at this point in the history
  • Loading branch information
BillSenior authored Nov 13, 2023
1 parent f003e73 commit 5bf8a0f
Show file tree
Hide file tree
Showing 19 changed files with 1,129 additions and 111 deletions.
4 changes: 4 additions & 0 deletions libs/MeshKernel/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ set(
${SRC_DIR}/Definitions.cpp
${SRC_DIR}/Entities.cpp
${SRC_DIR}/FlipEdges.cpp
${SRC_DIR}/Hessian.cpp
${SRC_DIR}/LandBoundaries.cpp
${SRC_DIR}/LandBoundary.cpp
${SRC_DIR}/Mesh.cpp
Expand All @@ -50,6 +51,7 @@ set(
${SRC_DIR}/PolygonalEnclosure.cpp
${SRC_DIR}/Polygons.cpp
${SRC_DIR}/RemoveDisconnectedRegions.cpp
${SRC_DIR}/SamplesHessianCalculator.cpp
${SRC_DIR}/Smoother.cpp
${SRC_DIR}/SplineAlgorithms.cpp
${SRC_DIR}/Splines.cpp
Expand Down Expand Up @@ -108,6 +110,7 @@ set(
${DOMAIN_INC_DIR}/Exceptions.hpp
${DOMAIN_INC_DIR}/FlipEdges.hpp
${DOMAIN_INC_DIR}/Formatting.hpp
${DOMAIN_INC_DIR}/Hessian.hpp
${DOMAIN_INC_DIR}/LandBoundaries.hpp
${DOMAIN_INC_DIR}/LandBoundary.hpp
${DOMAIN_INC_DIR}/Mesh.hpp
Expand All @@ -128,6 +131,7 @@ set(
${DOMAIN_INC_DIR}/Polygons.hpp
${DOMAIN_INC_DIR}/RangeCheck.hpp
${DOMAIN_INC_DIR}/RemoveDisconnectedRegions.hpp
${DOMAIN_INC_DIR}/SamplesHessianCalculator.hpp
${DOMAIN_INC_DIR}/Smoother.hpp
${DOMAIN_INC_DIR}/SplineAlgorithms.hpp
${DOMAIN_INC_DIR}/Splines.hpp
Expand Down
14 changes: 7 additions & 7 deletions libs/MeshKernel/include/MeshKernel/AveragingInterpolation.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -27,12 +27,12 @@

#pragma once

#include "MeshInterpolation.hpp"

#include <MeshKernel/AveragingStrategies/AveragingStrategy.hpp>
#include <MeshKernel/Constants.hpp>
#include <MeshKernel/Mesh2D.hpp>
#include <MeshKernel/Utilities/RTree.hpp>
#include "MeshKernel/AveragingStrategies/AveragingStrategy.hpp"
#include "MeshKernel/Constants.hpp"
#include "MeshKernel/Mesh2D.hpp"
#include "MeshKernel/MeshInterpolation.hpp"
#include "MeshKernel/Utilities/LinearAlgebra.hpp"
#include "MeshKernel/Utilities/RTree.hpp"

namespace meshkernel
{
Expand Down Expand Up @@ -129,7 +129,7 @@ namespace meshkernel
/// param[in] strategy The input strategy
/// param[in] searchPolygon The bounding polygon
/// @return The interpolated result
[[nodiscard]] double ComputeInterpolationResultFromNeighbors(std::unique_ptr<averaging::AveragingStrategy> strategy, std::vector<Point> const& searchPolygon);
[[nodiscard]] double ComputeInterpolationResultFromNeighbors(const Point& interpolationPoint, std::vector<Point> const& searchPolygon);

/// @brief Gets the sample value from an r-tree query
/// param[in] index The query index
Expand Down
9 changes: 9 additions & 0 deletions libs/MeshKernel/include/MeshKernel/Constants.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -33,8 +33,17 @@
#include <limits>
#include <math.h>

#if defined(__linux__) || defined(__APPLE__)
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
#endif

#include <Eigen/Core>

#if defined(__linux__) || defined(__APPLE__)
#pragma GCC diagnostic pop
#endif

namespace meshkernel
{

Expand Down
139 changes: 139 additions & 0 deletions libs/MeshKernel/include/MeshKernel/Hessian.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,139 @@
//---- GPL ---------------------------------------------------------------------
//
// Copyright (C) Stichting Deltares, 2011-2023.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation version 3.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// contact: [email protected]
// Stichting Deltares
// P.O. Box 177
// 2600 MH Delft, The Netherlands
//
// All indications and logos of, and references to, "Delft3D" and "Deltares"
// are registered trademarks of Stichting Deltares, and remain the property of
// Stichting Deltares. All rights reserved.
//
//------------------------------------------------------------------------------

#pragma once

#include <array>

#include "MeshKernel/Definitions.hpp"
#include "MeshKernel/Utilities/LinearAlgebra.hpp"

namespace meshkernel
{

/// @brief Array containing dimensions of the hessian
using HessianDimension = std::array<UInt, 3>;

/// @brief Define column major orientation
using MatrixColMajor = lin_alg::Matrix<double, Eigen::ColMajor>;

/// @brief The hessian values
///
/// Implemented as an array of matrices.
/// Not sure what is the best implementation yet for performance.
class Hessian
{
public:
/// @brief Default constructor
Hessian() = default;

/// @brief Constructor taking 3 parameters
Hessian(const UInt dim1, const UInt dim2, const UInt dim3);

/// @brief Resize taking 3 parameters
void resize(const UInt dim1, const UInt dim2, const UInt dim3);

/// @brief Get the dimension for each dimension
///
/// @param [in] dim For which dimension is the size required, dim in range [0,2]
UInt size(const UInt dim) const;

/// @brief Get all the Hessian dimensions
const HessianDimension& size() const;

/// @brief Get the 1-dimension index of
UInt get1DIndex(const UInt dim2, const UInt dim3) const;

/// @brief Get the value of the hessian
double operator()(const UInt dim1, const UInt dim2, const UInt dim3) const;

/// @brief Get the value of the hessian
double& operator()(const UInt dim1, const UInt dim2, const UInt dim3);

/// @brief Access the matrix in 'dim1' as though it were a 1 dimensional array
///
/// dim2 = i + size(1) * j?
double operator()(const UInt dim1, const UInt dim2) const;

/// @brief Get the matrix for a dimension
const MatrixColMajor& getMatrix(const UInt dim) const;

/// @brief Get the matrix for a dimension
MatrixColMajor& getMatrix(const UInt dim);

/// @brief Set all entries to zero.
void zero();

private:
/// @brief a vector of matrices for storing the hessian calculations
std::vector<MatrixColMajor> m_hessian;

/// @brief the initial hessian matrix dimension
HessianDimension m_dimensions{0, 0, 0};
};

} // namespace meshkernel

inline meshkernel::UInt meshkernel::Hessian::size(const UInt dim) const
{
return m_dimensions[dim];
}

inline const meshkernel::HessianDimension& meshkernel::Hessian::size() const
{
return m_dimensions;
}

inline meshkernel::UInt meshkernel::Hessian::get1DIndex(const UInt dim2, const UInt dim3) const
{
return dim2 + m_dimensions[1] * dim3;
}

inline double meshkernel::Hessian::operator()(const UInt dim1, const UInt dim2, const UInt dim3) const
{
return m_hessian[dim1](dim2, dim3);
}

inline double& meshkernel::Hessian::operator()(const UInt dim1, const UInt dim2, const UInt dim3)
{
return m_hessian[dim1](dim2, dim3);
}

inline double meshkernel::Hessian::operator()(const UInt dim1, const UInt dim2) const
{
return m_hessian[dim1](dim2);
}

inline const meshkernel::MatrixColMajor& meshkernel::Hessian::getMatrix(const UInt dim) const
{
return m_hessian[dim];
}

inline meshkernel::MatrixColMajor& meshkernel::Hessian::getMatrix(const UInt dim)
{
return m_hessian[dim];
}
19 changes: 17 additions & 2 deletions libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,8 @@ namespace meshkernel
enum class RefinementType
{
WaveCourant = 1,
RefinementLevels = 2
RefinementLevels = 2,
RidgeDetection = 3
};

public:
Expand Down Expand Up @@ -144,6 +145,21 @@ namespace meshkernel
/// @brief Computes the edge and face refinement mask from the minimum edge size
void ComputeRefinementMaskFromEdgeSize();

/// @brief Computes the refinement mask for refinement based on levels
void ComputeRefinementMasksForRefinementLevels(UInt face,
size_t& numberOfEdgesToRefine,
std::vector<UInt>& edgeToRefine) const;

/// @brief Computes the refinement mask for refinement based on wave courant criteria
void ComputeRefinementMasksForWaveCourant(UInt face,
size_t& numberOfEdgesToRefine,
std::vector<UInt>& edgeToRefine);

/// @brief Computes the refinement mask for refinement based on ridge detection
void ComputeRefinementMasksForRidgeDetection(UInt face,
size_t& numberOfEdgesToRefine,
std::vector<UInt>& edgeToRefine) const;

/// @brief Computes refinement masks (compute_jarefine_poly)
/// Face nodes, edge and edge lengths are stored in local caches. See Mesh2D.FaceClosedPolygon method
/// @param face The face index
Expand All @@ -154,7 +170,6 @@ namespace meshkernel

/// @brief Finds the hanging nodes in a face (find_hangingnodes)
/// @param[in] face The current face index
/// @returns The number of hanging edges on the face, the number of hanging nodes and the number of edges to refine
void FindHangingNodes(UInt face);

/// @brief Get the number of hanging nodes
Expand Down
2 changes: 1 addition & 1 deletion libs/MeshKernel/include/MeshKernel/Parameters.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -191,7 +191,7 @@ namespace meshkernel
range_check::CheckGreater(parameters.min_edge_size, 0.0, "Min edge size");
// Move enum meshkernel::MeshRefinement::RefinementType out of meshkernel::MeshRefinement
// then use the underlying type of the enums to define ValidMeshRefinementTypes
static std::vector<int> const ValidMeshRefinementTypes{1, 2};
static std::vector<int> const ValidMeshRefinementTypes{1, 2, 3};
range_check::CheckOneOf(parameters.refinement_type, ValidMeshRefinementTypes, "Refinement type");
range_check::CheckOneOf(parameters.connect_hanging_nodes, {0, 1}, "Connect hanging nodes");
range_check::CheckOneOf(parameters.account_for_samples_outside, {0, 1}, "Account for samples outside");
Expand Down
40 changes: 31 additions & 9 deletions libs/MeshKernel/include/MeshKernel/Point.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,11 @@ namespace meshkernel
{
}

/// @brief Set the point to be invalid
///
/// Both x and y values are set to the null, missing value
void SetInvalid();

/// @brief Inplace add point to point.
Point& operator+=(const Point& p);

Expand Down Expand Up @@ -117,11 +122,18 @@ namespace meshkernel

return !isInvalid;
}

/// @brief Set the point to be invalid.
void SetInvalid();
};

/// @brief Compute the dot product of a point and a vector.
///
/// This is mainly a convenience function
double dot(const Point& p, const Vector& v);

/// @brief Compute the dot product of a point and a vector.
///
/// This is mainly a convenience function
double dot(const Vector& v, const Point& p);

/// @brief Unary minus
///
/// @returns \f$ (-p1.x, -p1.y)\f$
Expand Down Expand Up @@ -242,6 +254,12 @@ namespace meshkernel

} // namespace meshkernel

inline void meshkernel::Point::SetInvalid()
{
x = constants::missing::doubleValue;
y = constants::missing::doubleValue;
}

inline meshkernel::Point& meshkernel::Point::operator+=(const Point& p)
{
x += p.x;
Expand Down Expand Up @@ -312,6 +330,16 @@ inline meshkernel::Point& meshkernel::Point::operator*=(const double p)
return *this;
}

inline double meshkernel::dot(const Point& p, const Vector& v)
{
return p.x * v.x() + p.y * v.y();
}

inline double meshkernel::dot(const Vector& v, const Point& p)
{
return dot(p, v);
}

inline meshkernel::Point meshkernel::operator-(const Point& pnt)
{
return Point(-pnt.x, -pnt.y);
Expand Down Expand Up @@ -416,9 +444,3 @@ inline meshkernel::Point meshkernel::PointAlongLine(const Point& startPoint, con
{
return (1.0 - lambda) * startPoint + lambda * endPoint;
}

void inline meshkernel::Point::SetInvalid()
{
x = constants::missing::doubleValue;
y = constants::missing::doubleValue;
}
Loading

0 comments on commit 5bf8a0f

Please sign in to comment.