Skip to content

Commit

Permalink
Merge pull request #349 from Exawind/refactor-LSFE
Browse files Browse the repository at this point in the history 
Refactor LSFE logic to add derivative of shape functions
  • Loading branch information
@deslaughter
deslaughter authored Feb 6, 2025
2 parents 4dd0e26 + 0f35398 commit 16fd9ea
Show file tree
Hide file tree
Showing 2 changed files with 59 additions and 17 deletions.
23 changes: 15 additions & 8 deletions src/math/least_squares_fit.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,8 @@

namespace openturbine {

using Matrix = std::vector<std::vector<double>>;

/**
* @brief Maps input geometric locations -> normalized domain using linear mapping
*
Expand Down Expand Up @@ -38,30 +40,36 @@ inline std::vector<double> MapGeometricLocations(const std::vector<double>& 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<std::vector<double>>, std::vector<double>> ShapeFunctionMatrices(
inline std::tuple<Matrix, Matrix, std::vector<double>> ShapeFunctionMatrices(
size_t n, size_t p, const std::vector<double>& 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<double> weights(p, 0.);
std::vector<std::vector<double>> shape_functions(p, std::vector<double>(n, 0.));
Matrix shape_functions(p, std::vector<double>(n, 0.));
Matrix derivative_functions(p, std::vector<double>(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};
}

/**
Expand All @@ -77,8 +85,7 @@ inline std::tuple<std::vector<std::vector<double>>, std::vector<double>> ShapeFu
* @return Interpolation coefficients (p x 3)
*/
inline std::vector<std::array<double, 3>> PerformLeastSquaresFitting(
size_t p, const std::vector<std::vector<double>>& shape_functions,
const std::vector<std::array<double, 3>>& points_to_fit
size_t p, const Matrix& shape_functions, const std::vector<std::array<double, 3>>& points_to_fit
) {
if (shape_functions.size() != p) {
throw std::invalid_argument("shape_functions rows do not match order p.");
Expand Down
53 changes: 44 additions & 9 deletions tests/unit_tests/math/test_least_squares_fit.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double> 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);
Expand All @@ -65,22 +65,39 @@ TEST(LeastSquaresFitTest, ShapeFunctionMatrices_FirstOrder) {

// Check shape function values at evaluation points
const std::vector<std::vector<double>> 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) {
for (size_t j = 0; j < phi_g[i].size(); ++j) {
EXPECT_NEAR(phi_g[i][j], expected[i][j], 1.e-15);
}
}

// 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 shape function derivative values at evaluation points
const std::vector<std::vector<double>> expected_dphi_g = {
{-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) {
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<double> 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);
Expand All @@ -96,16 +113,35 @@ TEST(LeastSquaresFitTest, ShapeFunctionMatrices_SecondOrder) {

// Check shape function values at evaluation points
const std::vector<std::vector<double>> 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) {
for (size_t j = 0; j < phi_g[i].size(); ++j) {
EXPECT_NEAR(phi_g[i][j], expected[i][j], 1.e-15);
}
}

// 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 shape function derivative values at evaluation points
const std::vector<std::vector<double>> expected_dphi_g = {
{-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) {
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) {
Expand All @@ -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);
Expand All @@ -138,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) {
Expand Down

0 comments on commit 16fd9ea

Please sign in to comment.