Skip to content

Commit

Permalink
Use Eigen library in curvilinear grid generation (#214 | GRIDEDIT-533…
Browse files Browse the repository at this point in the history 
…, GRIDEDIT-535, GRIDEDIT-594)
  • Loading branch information
@ahmad-el-sayed
ahmad-el-sayed authored Sep 13, 2023
1 parent 5a2f27c commit 8b98cdd
Show file tree
Hide file tree
Showing 44 changed files with 5,381 additions and 4,132 deletions.
13 changes: 8 additions & 5 deletions libs/MeshKernel/include/MeshKernel/Constants.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,8 @@
#include <limits>
#include <math.h>

#include <Eigen/Core>

namespace meshkernel
{

Expand All @@ -41,11 +43,12 @@ namespace meshkernel
// missing values
namespace missing
{
constexpr double innerOuterSeparator = -998.0; ///< Double value used to separate the inner part of a polygon from its outer part
constexpr double doubleValue = -999.0; ///< Double value used as missing value
constexpr int intValue = -999; ///< Integer value used as missing value
constexpr UInt uintValue = std::numeric_limits<UInt>::max(); ///< missing value used for invalid indices
} // namespace missing
constexpr double innerOuterSeparator = -998.0; ///< Double value used to separate the inner part of a polygon from its outer part
constexpr double doubleValue = -999.0; ///< Double value used as missing value
constexpr int intValue = -999; ///< Integer value used as missing value
constexpr UInt uintValue = std::numeric_limits<UInt>::max(); ///< missing value used for invalid indices
constexpr Eigen::Index EigenIndexValue = std::numeric_limits<Eigen::Index>::max(); ///< missing value used for invalid eigen indices
} // namespace missing

// often used values
namespace numeric
Expand Down
88 changes: 32 additions & 56 deletions libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@
#include <MeshKernel/Entities.hpp>
#include <MeshKernel/Exceptions.hpp>
#include <MeshKernel/Mesh.hpp>
#include <MeshKernel/Utilities/LinearAlgebra.hpp>

namespace meshkernel
{
Expand Down Expand Up @@ -74,16 +75,16 @@ namespace meshkernel
/// @brief Constructor taking only a projection
CurvilinearGrid(Projection projection);

/// @brief Lvalue constructor. Creates a new curvilinear grid from a given set of points
/// @param[in] grid The input grid points
/// @param[in] projection The projection to use
CurvilinearGrid(lin_alg::Matrix<Point> const& grid, Projection projection);

/// @brief Deletes a curvilinear grid inside a polygon
/// @param[in] polygons The polygons
/// @param[in] polygonIndex The index of the polygon to use for deletion
void Delete(std::shared_ptr<Polygons> polygons, UInt polygonIndex);

/// @brief Lvalue constructor. Creates a new curvilinear grid from a given set of points
/// @param[in] grid The input grid points
/// @param[in] projection The projection to use
CurvilinearGrid(std::vector<std::vector<Point>> const& grid, Projection projection);

/// @brief Check if current curvilinear grid instance is valid
/// @return True if valid, false otherwise
[[nodiscard]] bool IsValid() const;
Expand All @@ -99,23 +100,35 @@ namespace meshkernel
/// @param[in] point The input grid points
[[nodiscard]] CurvilinearGridNodeIndices GetNodeIndices(Point point);

/// @brief Get the grid node at the (i,j) location
Point& GetNode(const UInt i, const UInt j);
/// @brief Gets a constant reference to the grid node at the (i,j) location
[[nodiscard]] meshkernel::Point& GetNode(const UInt i, const UInt j) { return m_gridNodes(i, j); }

/// @brief Get the grid node at the (i,j) location
Point GetNode(const UInt i, const UInt j) const;
/// @brief Gets a reference to the grid node at the (i,j) location
[[nodiscard]] meshkernel::Point const& GetNode(const UInt i, const UInt j) const { return m_gridNodes(i, j); }

/// @brief Get the grid node at the location specified by the index.
///
/// @brief Gets a reference to the grid node at the location specified by the index.
/// @note Exception will be raised for a non-valid index
/// This is just a helper function, it calls GetNode with (index.m_m, index.m_n)
Point& GetNode(const CurvilinearGridNodeIndices& index);

/// @brief Get the grid node at the location specified by the index.
///
[[nodiscard]] meshkernel::Point& GetNode(const CurvilinearGridNodeIndices& index)
{
if (!index.IsValid())
{
throw ConstraintError("Invalid node index");
}
return m_gridNodes(index.m_m, index.m_n);
}

/// @brief Get a constant reference to the grid node at the location specified by the index.
/// @note Exception will be raised for a non-valid index
/// This is just a helper function, it calls GetNode with (index.m_m, index.m_n)
Point GetNode(const CurvilinearGridNodeIndices& index) const;
[[nodiscard]] meshkernel::Point const& GetNode(const CurvilinearGridNodeIndices& index) const
{
if (!index.IsValid())
{
throw ConstraintError("Invalid node index");
}
return m_gridNodes(index.m_m, index.m_n);
}

/// @brief From a point gets the node indices of the closest edges
/// @param[in] point The input point
Expand Down Expand Up @@ -203,9 +216,9 @@ namespace meshkernel

UInt m_numM = 0; ///< The number of m coordinates (vertical lines)
UInt m_numN = 0; ///< The number of n coordinates (horizontal lines)
std::vector<std::vector<Point>> m_gridNodes; ///< Member variable storing the grid
std::vector<std::vector<bool>> m_gridFacesMask; ///< The mask of the grid faces (true/false)
std::vector<std::vector<NodeType>> m_gridNodesTypes; ///< The grid node types
lin_alg::Matrix<Point> m_gridNodes; ///< Member variable storing the grid
lin_alg::Matrix<bool> m_gridFacesMask; ///< The mask of the grid faces (true/false)
lin_alg::Matrix<NodeType> m_gridNodesTypes; ///< The grid node types
std::vector<CurvilinearGridNodeIndices> m_gridIndices; ///< The original mapping of the flatten nodes in the curvilinear grid

private:
Expand All @@ -217,11 +230,6 @@ namespace meshkernel
/// @brief Computes the valid grid faces
void ComputeGridFacesMask();

/// @brief Compute spline derivatives along a gridline, also accounting for missing values
/// @param[in] gridLine The input gridline
/// @returns The spline derivatives
[[nodiscard]] std::vector<Point> ComputeSplineDerivatesAlongGridLine(const std::vector<Point>& gridLine) const;

/// @brief Adds an edge at the boundary forming a new face. Increase the grid if required (MODGR1)
/// The position of the new edge depends on the type of \p firstNode or \p secondNode.
/// For example, if one of the node types is 'Left' the new edge will be inserted on the left.
Expand All @@ -232,35 +240,3 @@ namespace meshkernel
CurvilinearGridNodeIndices const& secondNode);
};
} // namespace meshkernel

inline meshkernel::Point& meshkernel::CurvilinearGrid::GetNode(const UInt i, const UInt j)
{
return m_gridNodes[i][j];
}

inline meshkernel::Point meshkernel::CurvilinearGrid::GetNode(const UInt i, const UInt j) const
{
return m_gridNodes[i][j];
}

inline meshkernel::Point meshkernel::CurvilinearGrid::GetNode(const CurvilinearGridNodeIndices& index) const
{

if (!index.IsValid())
{
throw ConstraintError("Invalid node index");
}

return m_gridNodes[index.m_m][index.m_n];
}

inline meshkernel::Point& meshkernel::CurvilinearGrid::GetNode(const CurvilinearGridNodeIndices& index)
{

if (!index.IsValid())
{
throw ConstraintError("Invalid node index");
}

return m_gridNodes[index.m_m][index.m_n];
}
29 changes: 15 additions & 14 deletions libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridCreateUniform.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@

#include <MeshKernel/Entities.hpp>
#include <MeshKernel/Parameters.hpp>
#include <MeshKernel/Utilities/LinearAlgebra.hpp>

#include <memory>

Expand Down Expand Up @@ -99,13 +100,13 @@ namespace meshkernel
/// @param[in] blockSizeX The grid block size in x dimension
/// @param[in] blockSizeY The grid block size in y dimension
/// @returns[in] The coordinates of the grid point
static std::vector<std::vector<Point>> ComputeCartesian(const int numColumns,
const int numRows,
const double originX,
const double originY,
const double angle,
const double blockSizeX,
const double blockSizeY);
static lin_alg::Matrix<Point> ComputeCartesian(const int numColumns,
const int numRows,
const double originX,
const double originY,
const double angle,
const double blockSizeX,
const double blockSizeY);

/// @brief Compute an uniform curvilinear grid on spherical coordinates.
/// A correction to the longitudinal discretization is applied to preserve an aspect ratio ds/dy = 1 on real distances.
Expand All @@ -119,13 +120,13 @@ namespace meshkernel
/// @param[in] blockSizeX The grid block size in x dimension
/// @param[in] blockSizeY The grid block size in y dimension
/// @returns[in] The coordinates of the grid point
static std::vector<std::vector<Point>> ComputeSpherical(const int numColumns,
const int numRows,
const double originX,
const double originY,
const double angle,
const double blockSizeX,
const double blockSizeY);
static lin_alg::Matrix<Point> ComputeSpherical(const int numColumns,
const int numRows,
const double originX,
const double originY,
const double angle,
const double blockSizeX,
const double blockSizeY);

/// @brief Compute the adjusted latitude for keeping an aspect ratio of 1, considering the spherical coordinates
/// @param[in] blockSize The grid block size in y dimension
Expand Down
Loading

0 comments on commit 8b98cdd

Please sign in to comment.