Feature/gridedit 1360 improve clg generation from spline

extern/CMakeLists.txt

@@ -1 +1 @@
-add_subdirectory(triangle)
+add_subdirectory(triangle)

libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridFromSplines.hpp

@@ -132,30 +132,61 @@ namespace meshkernel
         std::shared_ptr<Splines> m_splines; ///< A pointer to spline class instance
 
     private:
+        /// @brief Tolerance used throughout.
+        const double tolerance = 1e-8;
+
+        static constexpr int MaxDegree = 2;
+
+        static constexpr int MaxDegreeP1 = MaxDegree + 1;
+
         /// @brief From the layerIndex index gets the next grid layerIndex and the transversal sublayer index (get_isub)
         /// @param[in] layerIndex The current grid layerIndex index
         /// @returns The next grid layerIndex and the sub layerIndex index
         std::tuple<UInt, UInt> ComputeGridLayerAndSubLayer(UInt layerIndex);
 
+        /// @brief Compute the time node x1 will cross thee segment with end points x3 and x4
+        double ComputeNodeSegmentCrossingTime(const Point& x1, const Point& x3, const Point& x4,
+                                              const Point& v1, const Point& v3, const Point& v4) const;
+
+        /// @brief Should the direct neighbour be included, circular grids also checked
+        bool IncludeDirectNeighbours(const UInt layerIndex,
+                                     const UInt loopIndex,
+                                     const lin_alg::Matrix<UInt>& gridPointsIndices,
+                                     const UInt indexLeftOfLeft,
+                                     const UInt indexRightOfRight) const;
+
+        /// @brief Compute maximum allowable grid layer growth time
+        void ComputeMaximumTimeStep(const UInt layerIndex,
+                                    const lin_alg::RowVector<Point>& activeLayerPoints,
+                                    const std::vector<Point>& velocityVectorAtGridPoints,
+                                    const std::vector<Point>& frontGridPoints,
+                                    const std::vector<Point>& frontVelocities,
+                                    const lin_alg::Matrix<UInt>& gridPointsIndices,
+                                    const double timeStep,
+                                    double& otherTimeStep,
+                                    std::vector<double>& otherTimeStepMax) const;
+
         /// @brief Grow a layer starting from a given layer index
         void GrowLayer(UInt layerIndex);
 
         /// @brief Compute the maximum allowable time step when growing a grid (comp_tmax_self)
         /// @param[in] coordinates The starting point coordinates
         /// @param[in] velocities  The velocities at the coordinates
         /// @returns The maximum allowable time step for each edge
-        std::vector<double> ComputeMaximumEdgeGrowTime(const std::vector<Point>& coordinates,
+        std::vector<double> ComputeMaximumEdgeGrowTime(const lin_alg::RowVector<Point>& coordinates,
                                                        const std::vector<Point>& velocities) const;
 
         /// @brief Copy growth velocities to the advancing front, add points at front corners corners (copy_vel_to_front)
         /// @param[in] layerIndex The current grid layerIndex index
         /// @param[in] previousFrontVelocities The previous front velocities
         /// @returns The front velocities for the next front
-        std::vector<Point> CopyVelocitiesToFront(UInt layerIndex, const std::vector<Point>& previousFrontVelocities);
+        // TODO update return comment and perhaps rename, it now does more than the name suggests
+        std::tuple<std::vector<Point>, std::vector<Point>, lin_alg::Matrix<UInt>> CopyVelocitiesToFront(UInt layerIndex,
+                                                                                                        const std::vector<Point>& previousFrontVelocities);
 
         /// @brief Computes the points at front, which have to be moved.
         /// @returns The indices of the grid points, the front grid points and the number of front points
-        std::tuple<lin_alg::Matrix<UInt>, lin_alg::RowVector<Point>, UInt> FindFront();
+        std::tuple<lin_alg::Matrix<UInt>, std::vector<Point>, UInt> FindFront() const;
 
         /// @brief Compute growth velocity vectors at grid points (comp_vel)
         /// @param[in] layerIndex The current grid layerIndex index
@@ -262,6 +293,21 @@ namespace meshkernel
         /// @brief Delete skewed cells and cells whose aspect ratio exceeds a prescribed value (postgrid)
         void DeleteSkinnyTriangles();
 
+        /// @brief Get the index of the centre spline
+        UInt GetCentralSplineIndex() const;
+
+        /// @brief Determine if segments crosses the centre spline
+        bool SegmentCrossesCentreSpline(const Point& x1, const Point& x2) const;
+
+        void RootFinder(const std::array<double, MaxDegreeP1>& coefficients,
+                        int& degree,
+                        std::array<double, MaxDegree>& realRoots,
+                        std::array<double, MaxDegree>& imaginaryRoots) const;
+
+        /// @brief Find the roots of up to a quadratic polynomial
+        void SolveQuadratic(const std::array<double, MaxDegreeP1>& coeffs,
+                            std::array<double, MaxDegree>& roots) const;
+
         /// @brief Allocate spline properties arrays
         void AllocateSplinesProperties();
 

libs/MeshKernel/include/MeshKernel/Point.hpp

@@ -122,6 +122,16 @@ namespace meshkernel
         }
     };
 
+    /// @brief Compute the dot product of a point with itself.
+    ///
+    /// This is mainly a convenience function
+    double length(const Point& p1);
+
+    /// @brief Compute the dot product of two points.
+    ///
+    /// This is mainly a convenience function
+    double dot(const Point& p1, const Point& p2);
+
     /// @brief Compute the dot product of a point and a vector.
     ///
     /// This is mainly a convenience function
@@ -132,6 +142,11 @@ namespace meshkernel
     /// This is mainly a convenience function
     double dot(const Vector& v, const Point& p);
 
+    /// @brief Compute the cross product of two points (as vectors).
+    ///
+    /// This is mainly a convenience function
+    double cross(const Point& p1, const Point& p2);
+
     /// @brief Unary minus
     ///
     /// @returns \f$ (-p1.x, -p1.y)\f$
@@ -328,6 +343,16 @@ inline meshkernel::Point& meshkernel::Point::operator*=(const double p)
     return *this;
 }
 
+inline double meshkernel::length(const Point& p)
+{
+    return p.x * p.x + p.y * p.y;
+}
+
+inline double meshkernel::dot(const Point& p1, const Point& p2)
+{
+    return p1.x * p2.x + p1.y * p2.y;
+}
+
 inline double meshkernel::dot(const Point& p, const Vector& v)
 {
     return p.x * v.x() + p.y * v.y();
@@ -338,6 +363,11 @@ inline double meshkernel::dot(const Vector& v, const Point& p)
     return dot(p, v);
 }
 
+inline double meshkernel::cross(const Point& p1, const Point& p2)
+{
+    return p1.x * p2.y - p1.y * p2.x;
+}
+
 inline meshkernel::Point meshkernel::operator-(const Point& pnt)
 {
     return Point(-pnt.x, -pnt.y);