From 5de2b3ddb2742526016b97805a1751d1057d4237 Mon Sep 17 00:00:00 2001 From: faisal-bhuiyan Date: Tue, 4 Feb 2025 17:34:28 -0700 Subject: [PATCH 1/2] Add derivatives of shape function to ShapeFunctionMatrices() --- src/math/least_squares_fit.hpp | 23 ++++++---- .../math/test_least_squares_fit.cpp | 42 +++++++++++++++++-- 2 files changed, 54 insertions(+), 11 deletions(-) diff --git a/src/math/least_squares_fit.hpp b/src/math/least_squares_fit.hpp index 402fdb70..1f083871 100644 --- a/src/math/least_squares_fit.hpp +++ b/src/math/least_squares_fit.hpp @@ -11,6 +11,8 @@ namespace openturbine { +using Matrix = std::vector>; + /** * @brief Maps input geometric locations -> normalized domain using linear mapping * @@ -38,30 +40,36 @@ inline std::vector MapGeometricLocations(const std::vector& geom } /** - * @brief Computes shape function matrices ϕg at points ξg + * @brief Computes shape function matrices ϕg and their derivatives dϕg at points ξg * * @param n Number of geometric points to fit (>=2) * @param p Number of points representing the polynomial of order p-1 (2 <= p <= n) * @param evaluation_pts Evaluation points in [-1, 1] - * @return Tuple containing shape function matrix and GLL points + * @return Tuple containing shape function matrix, derivative matrix, and GLL points */ -inline std::tuple>, std::vector> ShapeFunctionMatrices( +inline std::tuple> ShapeFunctionMatrices( size_t n, size_t p, const std::vector& evaluation_pts ) { // Compute GLL points which will act as the nodes for the shape functions auto gll_pts = GenerateGLLPoints(p - 1); - // Compute weights for the shape functions at the evaluation points + // Compute weights for the shape functions and its derivatives at the evaluation points std::vector weights(p, 0.); - std::vector> shape_functions(p, std::vector(n, 0.)); + Matrix shape_functions(p, std::vector(n, 0.)); + Matrix derivative_functions(p, std::vector(n, 0.)); for (size_t j = 0; j < n; ++j) { LagrangePolynomialInterpWeights(evaluation_pts[j], gll_pts, weights); for (size_t k = 0; k < p; ++k) { shape_functions[k][j] = weights[k]; } + + LagrangePolynomialDerivWeights(evaluation_pts[j], gll_pts, weights); + for (size_t k = 0; k < p; ++k) { + derivative_functions[k][j] = weights[k]; + } } - return {shape_functions, gll_pts}; + return {shape_functions, derivative_functions, gll_pts}; } /** @@ -77,8 +85,7 @@ inline std::tuple>, std::vector> ShapeFu * @return Interpolation coefficients (p x 3) */ inline std::vector> PerformLeastSquaresFitting( - size_t p, const std::vector>& shape_functions, - const std::vector>& points_to_fit + size_t p, const Matrix& shape_functions, const std::vector>& points_to_fit ) { if (shape_functions.size() != p) { throw std::invalid_argument("shape_functions rows do not match order p."); diff --git a/tests/unit_tests/math/test_least_squares_fit.cpp b/tests/unit_tests/math/test_least_squares_fit.cpp index 456cb6f5..a3426953 100644 --- a/tests/unit_tests/math/test_least_squares_fit.cpp +++ b/tests/unit_tests/math/test_least_squares_fit.cpp @@ -51,7 +51,7 @@ TEST(LeastSquaresFitTest, ShapeFunctionMatrices_FirstOrder) { const size_t n{3}; // Number of pts to fit const size_t p{2}; // Polynomial order + 1 const std::vector xi_g = {-1., 0., 1.}; // Evaluation points - const auto [phi_g, gll_pts] = openturbine::ShapeFunctionMatrices(n, p, xi_g); + const auto [phi_g, dphi_g, gll_pts] = openturbine::ShapeFunctionMatrices(n, p, xi_g); // Check GLL points (2 at -1 and 1) ASSERT_EQ(gll_pts.size(), p); @@ -74,13 +74,30 @@ TEST(LeastSquaresFitTest, ShapeFunctionMatrices_FirstOrder) { EXPECT_NEAR(phi_g[i][j], expected[i][j], 1.e-15); } } + + // Check derivative shape function matrix dimensions (2 x 3) + ASSERT_EQ(dphi_g.size(), p); + ASSERT_EQ(dphi_g[0].size(), n); + ASSERT_EQ(dphi_g[1].size(), n); + + // Check derivative shape function values at evaluation points + const std::vector> expected_dphi_g = { + {-0.5, -0.5, -0.5}, // First derivative shape function + {0.5, 0.5, 0.5} // Second derivative shape function + }; + + for (size_t i = 0; i < dphi_g.size(); ++i) { + for (size_t j = 0; j < dphi_g[i].size(); ++j) { + EXPECT_NEAR(dphi_g[i][j], expected_dphi_g[i][j], 1.e-15); + } + } } TEST(LeastSquaresFitTest, ShapeFunctionMatrices_SecondOrder) { const size_t n{5}; // Number of pts to fit const size_t p{3}; // Polynomial order + 1 const std::vector xi_g = {-1., -0.5, 0., 0.5, 1.}; // Evaluation points - const auto [phi_g, gll_pts] = openturbine::ShapeFunctionMatrices(n, p, xi_g); + const auto [phi_g, dphi_g, gll_pts] = openturbine::ShapeFunctionMatrices(n, p, xi_g); // Check GLL points (3 at -1, 0, and 1) ASSERT_EQ(gll_pts.size(), 3); @@ -106,6 +123,25 @@ TEST(LeastSquaresFitTest, ShapeFunctionMatrices_SecondOrder) { EXPECT_NEAR(phi_g[i][j], expected[i][j], 1.e-15); } } + + // Check derivative shape function matrix dimensions (3 x 5) + ASSERT_EQ(dphi_g.size(), p); + for (const auto& row : dphi_g) { + ASSERT_EQ(row.size(), 5); + } + + // Check derivative shape function values at evaluation points + const std::vector> expected_dphi_g = { + {-1.5, -1.0, -0.5, 0.0, 0.5}, // First derivative shape function + {2.0, 1.0, 0.0, -1.0, -2.0}, // Second derivative shape function + {-0.5, 0.0, 0.5, 1.0, 1.5} // Third derivative shape function + }; + + for (size_t i = 0; i < dphi_g.size(); ++i) { + for (size_t j = 0; j < dphi_g[i].size(); ++j) { + EXPECT_NEAR(dphi_g[i][j], expected_dphi_g[i][j], 1.e-15); + } + } } TEST(LeastSquaresFitTest, FitsParametricCurve) { @@ -125,7 +161,7 @@ TEST(LeastSquaresFitTest, FitsParametricCurve) { // Step 2: Generate shape function matrices (using p = 4 i.e. cubic interpolation) const size_t n = input_points.size(); const size_t p = 4; - const auto [phi_g, gll_points] = ShapeFunctionMatrices(n, p, mapped_locations); + const auto [phi_g, dphi_g, gll_points] = ShapeFunctionMatrices(n, p, mapped_locations); // Step 3: Perform least squares fitting const auto X = PerformLeastSquaresFitting(p, phi_g, input_points); From 0f35398a79c1e5e96ea009fbeb172899f3ee7227 Mon Sep 17 00:00:00 2001 From: faisal-bhuiyan Date: Wed, 5 Feb 2025 11:00:17 -0700 Subject: [PATCH 2/2] Fix some comments in unit tests --- .../math/test_least_squares_fit.cpp | 29 +++++++++---------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/tests/unit_tests/math/test_least_squares_fit.cpp b/tests/unit_tests/math/test_least_squares_fit.cpp index a3426953..29e0feb9 100644 --- a/tests/unit_tests/math/test_least_squares_fit.cpp +++ b/tests/unit_tests/math/test_least_squares_fit.cpp @@ -65,8 +65,8 @@ TEST(LeastSquaresFitTest, ShapeFunctionMatrices_FirstOrder) { // Check shape function values at evaluation points const std::vector> expected = { - {1., 0.5, 0.}, // First shape function - {0., 0.5, 1.} // Second shape function + {1., 0.5, 0.}, // row 1 + {0., 0.5, 1.} // row 2 }; for (size_t i = 0; i < phi_g.size(); ++i) { @@ -75,15 +75,15 @@ TEST(LeastSquaresFitTest, ShapeFunctionMatrices_FirstOrder) { } } - // Check derivative shape function matrix dimensions (2 x 3) + // Check shape function derivative matrix dimensions (2 x 3) ASSERT_EQ(dphi_g.size(), p); ASSERT_EQ(dphi_g[0].size(), n); ASSERT_EQ(dphi_g[1].size(), n); - // Check derivative shape function values at evaluation points + // Check shape function derivative values at evaluation points const std::vector> expected_dphi_g = { - {-0.5, -0.5, -0.5}, // First derivative shape function - {0.5, 0.5, 0.5} // Second derivative shape function + {-0.5, -0.5, -0.5}, // row 1 + {0.5, 0.5, 0.5} // row 2 }; for (size_t i = 0; i < dphi_g.size(); ++i) { @@ -113,9 +113,9 @@ TEST(LeastSquaresFitTest, ShapeFunctionMatrices_SecondOrder) { // Check shape function values at evaluation points const std::vector> expected = { - {1.0, 0.375, 0.0, -0.125, 0.0}, // First shape function - {0.0, 0.75, 1.0, 0.75, 0.0}, // Second shape function - {0.0, -0.125, 0.0, 0.375, 1.0} // Third shape function + {1., 0.375, 0., -0.125, 0.}, // row 1 + {0., 0.75, 1., 0.75, 0.}, // row 2 + {0., -0.125, 0., 0.375, 1.} // row 3 }; for (size_t i = 0; i < phi_g.size(); ++i) { @@ -124,17 +124,17 @@ TEST(LeastSquaresFitTest, ShapeFunctionMatrices_SecondOrder) { } } - // Check derivative shape function matrix dimensions (3 x 5) + // Check shape function derivative matrix dimensions (3 x 5) ASSERT_EQ(dphi_g.size(), p); for (const auto& row : dphi_g) { ASSERT_EQ(row.size(), 5); } - // Check derivative shape function values at evaluation points + // Check shape function derivative values at evaluation points const std::vector> expected_dphi_g = { - {-1.5, -1.0, -0.5, 0.0, 0.5}, // First derivative shape function - {2.0, 1.0, 0.0, -1.0, -2.0}, // Second derivative shape function - {-0.5, 0.0, 0.5, 1.0, 1.5} // Third derivative shape function + {-1.5, -1., -0.5, 0., 0.5}, // row 1 + {2., 1., 0., -1., -2.}, // row 2 + {-0.5, 0., 0.5, 1., 1.5} // row 3 }; for (size_t i = 0; i < dphi_g.size(); ++i) { @@ -174,7 +174,6 @@ TEST(LeastSquaresFitTest, FitsParametricCurve) { {5., 1., -1.} // Last point - same as input }; - // Verify results ASSERT_EQ(X.size(), expected_coefficients.size()); for (size_t i = 0; i < X.size(); ++i) { for (size_t j = 0; j < 3; ++j) {