refinement in spherical coordinates (GRIDEDIT-556 | #170)

libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridCreateUniform.hpp

@@ -43,20 +43,33 @@ namespace meshkernel
     public:
         /// @brief Class constructor
         ///
-        /// @param[in] MakeGridParameters The structure containing the make grid parameters
         /// @param[in] projection The projection to use
-        CurvilinearGridCreateUniform(const MakeGridParameters& MakeMeshParameters, Projection projection);
+        CurvilinearGridCreateUniform(Projection projection);
 
         /// @brief Compute an uniform curvilinear grid using the make mesh parameters
-        CurvilinearGrid Compute() const;
+        /// @param[in] makeGridParameters The parameters to use for the curvilinear grid generation
+        CurvilinearGrid Compute(const MakeGridParameters& makeGridParameters) const;
 
         /// @brief Compute an uniform curvilinear grid in one polygon
+        /// @param[in] makeGridParameters The parameters to use for the curvilinear grid generation
         /// @param[in] polygons The input polygons
         /// @param[in] polygonIndex The polygon index
-        CurvilinearGrid Compute(std::shared_ptr<Polygons> polygons, size_t polygonIndex) const;
+        CurvilinearGrid Compute(const MakeGridParameters& makeGridParameters,
+                                std::shared_ptr<Polygons> polygons,
+                                size_t polygonIndex) const;
 
     private:
-        MakeGridParameters m_makeGridParameters; ///< A copy of the structure containing the parameters used for making the grid
-        Projection m_projection;                 ///< The projection to use
+        /// @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.
+        /// For preventing the creation of small edges, the correction stops when the distance is less than 2000 meters and the grid is generated around the poles.
+        /// This is a customized fix for GTSM models.
+        /// @param[in] makeGridParameters The parameters to use for the curvilinear grid generation
+        static std::vector<std::vector<Point>> ComputeSpherical(const MakeGridParameters& makeGridParameters);
+
+        /// @brief Compute an uniform curvilinear grid on cartesian coordinates.
+        /// @param[in] makeGridParameters The parameters to use for the curvilinear grid generation
+        static std::vector<std::vector<Point>> ComputeCartesian(const MakeGridParameters& makeGridParameters);
+
+        Projection m_projection; ///< The projection to use
     };
 } // namespace meshkernel

libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp

@@ -142,22 +142,31 @@ namespace meshkernel
         /// @brief Computes the edge and face refinement mask from sample values (compute_jarefine_poly)
         void ComputeRefinementMasksFromSamples();
 
-        /// @brief Computes the number of edges that should be refined in a face (compute_jarefine_poly)
+        /// @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 number of face nodes
-        /// @param refineEdgeCache 1 if the edge should be refined, 0 otherwise
-        /// @param numEdgesToBeRefined The computed number of edges to refined
-        void ComputeEdgesRefinementMaskFromSamples(size_t face,
-                                                   std::vector<size_t>& refineEdgeCache,
-                                                   size_t& numEdgesToBeRefined);
+        /// @param face The face index
+        void ComputeRefinementMasksFromSamples(size_t face);
 
         /// Computes the edge refinement mask (comp_jalink)
         void ComputeEdgesRefinementMask();
 
         /// @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
-        std::tuple<size_t, size_t, size_t> FindHangingNodes(size_t face);
+        void FindHangingNodes(size_t face);
+
+        /// @brief Get the number of hanging nodes
+        /// @returns The number of hanging nodes
+        size_t CountHangingNodes() const;
+
+        /// @brief Get the number of hanging nodes
+        /// @returns The number of hanging nodes
+        size_t CountHangingEdges() const;
+
+        /// @brief Get the number of hanging nodes
+        /// @param[in] face The current face index
+        /// @returns The number of hanging nodes
+        size_t CountEdgesToRefine(size_t face) const;
 
         /// Deletes isolated hanging nodes(remove_isolated_hanging_nodes)
         /// @returns Number of deleted isolated hanging nodes
@@ -166,8 +175,8 @@ namespace meshkernel
         /// @brief Connect the hanging nodes with triangles (connect_hanging_nodes)
         void ConnectHangingNodes();
 
-        /// @brief Smooth the face refinement mask (smooth_jarefine)
-        void SmoothEdgeRefinementMask() const;
+        /// @brief Smooth the face and edge refinement masks (smooth_jarefine)
+        void SmoothRefinementMasks();
 
         /// @brief Computes m_faceMask, if a face must be split later on (split_cells)
         void ComputeIfFaceShouldBeSplit();
@@ -206,10 +215,10 @@ namespace meshkernel
         std::vector<Point> m_polygonNodesCache;       ///< Cache for maintaining polygon nodes
         std::vector<size_t> m_localNodeIndicesCache;  ///< Cache for maintaining local node indices
         std::vector<size_t> m_globalEdgeIndicesCache; ///< Cache for maintaining edge indices
+        std::vector<size_t> m_refineEdgeCache;        ///< Cache for the edges to be refined
 
         RefinementType m_refinementType = RefinementType::WaveCourant; ///< The type of refinement to use
         bool m_directionalRefinement = false;                          ///< Whether there is directional refinement
-        bool m_useMassCenters = false;                                 ///< Split cells on the mass centers
 
         std::shared_ptr<Mesh2D> m_mesh;                             ///< Pointer to the mesh
         std::shared_ptr<MeshInterpolation> m_interpolant = nullptr; ///< Pointer to the AveragingInterpolation instance

libs/MeshKernel/include/MeshKernel/Parameters.hpp

@@ -124,7 +124,7 @@ namespace meshkernel
         /// Whether to use the mass center when splitting a face in the refinement process (yes=1/no=0)
         int use_mass_center_when_refining = 1;
 
-        /// @brief Minimum edge size
+        /// @brief Minimum edge size in meters
         double min_edge_size = 0.5;
 
         /// @brief Refinement criterion type
@@ -135,6 +135,15 @@ namespace meshkernel
 
         /// @brief Take samples outside face into account , 1 yes 0 no
         int account_for_samples_outside = 0;
+
+        /// @brief The number of smoothing iterations
+        int smoothing_iterations = 5;
+
+        /// @brief Maximum courant time in seconds
+        double max_courant_time = 120.0;
+
+        /// @brief Directional refinement, cannot be used when the number of smoothing iterations is larger than 0
+        int directional_refinement = 0;
     };
 
     /// @brief A struct used to describe the orthogonalization parameters in a C-compatible manner
@@ -159,4 +168,4 @@ namespace meshkernel
         double areal_to_angle_smoothing_factor = 1.0;
     };
 
-} // namespace meshkernel
+} // namespace meshkernel

libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridCreateUniform.cpp

@@ -25,6 +25,8 @@
 //
 //------------------------------------------------------------------------------
 
+#include "MeshKernel/Exceptions.hpp"
+
 #include <MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp>
 #include <MeshKernel/CurvilinearGrid/CurvilinearGridCreateUniform.hpp>
 #include <MeshKernel/Operations.hpp>
@@ -33,43 +35,105 @@
 using meshkernel::CurvilinearGrid;
 using meshkernel::CurvilinearGridCreateUniform;
 
-CurvilinearGridCreateUniform::CurvilinearGridCreateUniform(const MakeGridParameters& makeGridParameters, Projection projection)
-    : m_makeGridParameters(makeGridParameters), m_projection(projection)
+CurvilinearGridCreateUniform::CurvilinearGridCreateUniform(Projection projection) : m_projection(projection)
 {
 }
 
-CurvilinearGrid CurvilinearGridCreateUniform::Compute() const
+CurvilinearGrid CurvilinearGridCreateUniform::Compute(const MakeGridParameters& makeGridParameters) const
+{
+    if (m_projection == Projection::spherical)
+    {
+        return CurvilinearGrid{ComputeSpherical(makeGridParameters), m_projection};
+    }
+    if (m_projection == Projection::cartesian)
+    {
+        return CurvilinearGrid{ComputeCartesian(makeGridParameters), m_projection};
+    }
+
+    const std::string message = "Projection value: " + std::to_string(static_cast<int>(m_projection)) + " not supported";
+    throw NotImplemented(message);
+}
+
+std::vector<std::vector<meshkernel::Point>> CurvilinearGridCreateUniform::ComputeCartesian(const MakeGridParameters& makeGridParameters)
+
 {
-    // regular grid
-    auto numM = static_cast<size_t>(m_makeGridParameters.num_columns + 1);
-    auto numN = static_cast<size_t>(m_makeGridParameters.num_rows + 1);
-    const double XGridBlockSize = m_makeGridParameters.block_size_x;
-    const double YGridBlockSize = m_makeGridParameters.block_size_y;
-    const double cosineAngle = std::cos(m_makeGridParameters.angle * constants::conversion::degToRad);
-    const double sinAngle = std::sin(m_makeGridParameters.angle * constants::conversion::degToRad);
-    double OriginXCoordinate = m_makeGridParameters.origin_x;
-    double OriginYCoordinate = m_makeGridParameters.origin_y;
-
-    // No polygon, use MakeGridParameters as is
-
-    std::vector<std::vector<Point>> gridNodes(numN, std::vector<Point>(numM));
-    for (size_t n = 0; n < numN; ++n)
+    if (makeGridParameters.num_columns <= 0)
     {
+        throw AlgorithmError("Number of columns cannot be <= 0");
+    }
+    if (makeGridParameters.num_rows <= 0)
+    {
+        throw AlgorithmError("Number of rows cannot be <= 0");
+    }
+
+    const auto angleInRad = makeGridParameters.angle * constants::conversion::degToRad;
+    const auto cosineAngle = std::cos(angleInRad);
+    const auto sinAngle = std::sin(angleInRad);
+    const auto numM = makeGridParameters.num_columns + 1;
+    const auto numN = makeGridParameters.num_rows + 1;
+
+    std::vector<std::vector<Point>> result(numN, std::vector<Point>(numM));
+    const auto blockSizeXByCos = makeGridParameters.block_size_x * cosineAngle;
+    const auto blockSizeYbySin = makeGridParameters.block_size_y * sinAngle;
+    const auto blockSizeXBySin = makeGridParameters.block_size_x * sinAngle;
+    const auto blockSizeYByCos = makeGridParameters.block_size_y * cosineAngle;
+    for (int n = 0; n < numN; ++n)
+    {
+        for (int m = 0; m < numM; ++m)
+        {
+            const double newPointXCoordinate = makeGridParameters.origin_x + m * blockSizeXByCos - n * blockSizeYbySin;
+            const double newPointYCoordinate = makeGridParameters.origin_y + m * blockSizeXBySin + n * blockSizeYByCos;
+            result[n][m] = {newPointXCoordinate, newPointYCoordinate};
+        }
+    }
+    return result;
+}
+
+std::vector<std::vector<meshkernel::Point>> CurvilinearGridCreateUniform::ComputeSpherical(const MakeGridParameters& makeGridParameters)
+{
+    std::vector result = ComputeCartesian(makeGridParameters);
+
+    const auto numM = result[0].size();
+    const auto numN = result.size();
+    constexpr double latitudeCloseToPoles = 88.0; // The latitude defining close to poles
+    constexpr double minimumDistance = 2000;      // When the real distance along the latitude becomes smaller than minimumDistance and the location is close to the poles, snap the next point to the poles.
+    bool onPoles = false;
+
+    for (size_t n = 1; n < numN; ++n)
+    {
+        size_t lastRowOnPole = numM;
         for (size_t m = 0; m < numM; ++m)
         {
-            const double newPointXCoordinate = OriginXCoordinate + m * XGridBlockSize * cosineAngle - n * YGridBlockSize * sinAngle;
-            double newPointYCoordinate = OriginYCoordinate + m * XGridBlockSize * sinAngle + n * YGridBlockSize * cosineAngle;
-            if (m_projection == Projection::spherical && n > 0)
+            // adjust the latitude to preserve an aspect ratio of 1
+            const auto latitudeInRadiants = cos(constants::conversion::degToRad * result[n - 1][m].y);
+            const auto asp = latitudeInRadiants + (1.0 - latitudeInRadiants) * 0.3;
+            const auto dy = makeGridParameters.block_size_x * latitudeInRadiants * asp;
+
+            const auto newPointYCoordinate = result[n - 1][m].y + dy;
+            result[n][m].y = newPointYCoordinate;
+
+            // prevent too small dy increments in case we are on the poles
+            const double distance = dy * constants::conversion::degToRad * constants::geometric::earth_radius;
+            if (abs(result[n][m].y) > latitudeCloseToPoles && distance < minimumDistance)
             {
-                newPointYCoordinate = XGridBlockSize * cos(constants::conversion::degToRad * gridNodes[n - 1][m].y);
+                const double sign = result[n][m].y < 0 ? -1.0 : 1.0;
+                result[n][m].y = 90.0 * sign;
+                onPoles = true;
+                lastRowOnPole = n;
             }
-            gridNodes[n][m] = {newPointXCoordinate, newPointYCoordinate};
+        }
+        if (onPoles)
+        {
+            result.erase(result.begin() + lastRowOnPole + 1, result.end());
+            break;
         }
     }
-    return CurvilinearGrid{gridNodes, m_projection};
+    return result;
 }
 
-CurvilinearGrid CurvilinearGridCreateUniform::Compute(std::shared_ptr<Polygons> polygons, size_t polygonIndex) const
+CurvilinearGrid CurvilinearGridCreateUniform::Compute(const MakeGridParameters& makeGridParameters,
+                                                      std::shared_ptr<Polygons> polygons,
+                                                      size_t polygonIndex) const
 {
     if (polygons->GetProjection() != m_projection)
     {
@@ -99,11 +163,9 @@ CurvilinearGrid CurvilinearGridCreateUniform::Compute(std::shared_ptr<Polygons>
     double etamin = std::numeric_limits<double>::max();
     double etamax = -etamin;
 
-    const double XGridBlockSize = m_makeGridParameters.block_size_x;
-    const double YGridBlockSize = m_makeGridParameters.block_size_y;
-    const double cosineAngle = std::cos(m_makeGridParameters.angle * constants::conversion::degToRad);
-    const double sinAngle = std::sin(m_makeGridParameters.angle * constants::conversion::degToRad);
-
+    const auto angleInRad = makeGridParameters.angle * constants::conversion::degToRad;
+    const auto cosineAngle = std::cos(angleInRad);
+    const auto sinAngle = std::sin(angleInRad);
     Projection const polygonProjection = polygons->GetProjection();
 
     for (auto i = startPolygonIndex; i <= endPolygonIndex; ++i)
@@ -130,28 +192,32 @@ CurvilinearGrid CurvilinearGridCreateUniform::Compute(std::shared_ptr<Polygons>
         yShift = yShift / (constants::geometric::earth_radius * std::cos(referencePoint.y * constants::conversion::degToRad)) * constants::conversion::radToDeg;
     }
 
-    auto const OriginXCoordinate = referencePoint.x + xShift;
-    auto const OriginYCoordinate = referencePoint.y + yShift;
-    auto const numN = static_cast<size_t>(std::ceil((etamax - etamin) / XGridBlockSize) + 1);
-    auto const numM = static_cast<size_t>(std::ceil((xmax - xmin) / YGridBlockSize) + 1);
+    MakeGridParameters makeGridParametersInPolygon(makeGridParameters);
+
+    makeGridParametersInPolygon.origin_x = referencePoint.x + xShift;
+    makeGridParametersInPolygon.origin_y = referencePoint.y + yShift;
+    const int numM = static_cast<int>(std::ceil((etamax - etamin) / makeGridParameters.block_size_x) + 1);
+    makeGridParametersInPolygon.num_columns = numM > 0 ? numM - 1 : 0;
+    const int numN = static_cast<int>(std::ceil((xmax - xmin) / makeGridParameters.block_size_y) + 1);
+    makeGridParametersInPolygon.num_rows = numN > 0 ? numN - 1 : 0;
 
-    std::vector<std::vector<Point>> gridNodes(numN, std::vector<Point>(numM));
-    for (size_t n = 0; n < numN; ++n)
+    CurvilinearGrid curvilinearGrid;
+    switch (m_projection)
     {
-        for (size_t m = 0; m < numM; ++m)
-        {
-            const double newPointXCoordinate = OriginXCoordinate + m * XGridBlockSize * cosineAngle - n * YGridBlockSize * sinAngle;
-            double newPointYCoordinate = OriginYCoordinate + m * XGridBlockSize * sinAngle + n * YGridBlockSize * cosineAngle;
-            if (polygonProjection == Projection::spherical && n > 0)
-            {
-                newPointYCoordinate = XGridBlockSize * cos(constants::conversion::degToRad * gridNodes[n - 1][m].y);
-            }
-            gridNodes[n][m] = {newPointXCoordinate, newPointYCoordinate};
-        }
+    case Projection::spherical:
+        curvilinearGrid = CurvilinearGrid{ComputeSpherical(makeGridParametersInPolygon),
+                                          m_projection};
+        break;
+    case Projection::cartesian:
+        curvilinearGrid = CurvilinearGrid{ComputeCartesian(makeGridParametersInPolygon),
+                                          m_projection};
+        break;
+    case Projection::sphericalAccurate:
+    default:
+        const std::string message = "Projection value: " + std::to_string(static_cast<int>(m_projection)) + " not supported";
+        throw NotImplemented(message);
     }
 
-    CurvilinearGrid curvilinearGrid(gridNodes, m_projection);
-
     // remove nodes outside the polygon
     curvilinearGrid.Delete(polygons, polygonIndex);