From 80ea9759da90a29e042c3d3e9bbe508f67674940 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Fri, 31 Jul 2026 01:48:33 +0900 Subject: [PATCH] =?UTF-8?q?feat(fem-and-beam-kernel):=20step=200=20?= =?UTF-8?q?=E2=80=94=20quadrature-and-shape-functions?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CMakeLists.txt | 2 + include/fesa/fem/gauss_rule.hpp | 14 ++++++ include/fesa/fem/line2_shape.hpp | 13 ++++++ src/fesa/fem/gauss_rule.cpp | 33 ++++++++++++++ src/fesa/fem/line2_shape.cpp | 31 +++++++++++++ tests/CMakeLists.txt | 38 +++++++++++++++ tests/unit/fem/gauss_rule_test.cpp | 53 +++++++++++++++++++++ tests/unit/fem/line2_shape_test.cpp | 71 +++++++++++++++++++++++++++++ 8 files changed, 255 insertions(+) create mode 100644 include/fesa/fem/gauss_rule.hpp create mode 100644 include/fesa/fem/line2_shape.hpp create mode 100644 src/fesa/fem/gauss_rule.cpp create mode 100644 src/fesa/fem/line2_shape.cpp create mode 100644 tests/unit/fem/gauss_rule_test.cpp create mode 100644 tests/unit/fem/line2_shape_test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 29e37f0..500731a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -27,6 +27,8 @@ include(cmake/FesaDependencies.cmake) add_library(fesa_core STATIC src/fesa/core/version.cpp + src/fesa/fem/gauss_rule.cpp + src/fesa/fem/line2_shape.cpp src/fesa/io/abaqus/parser.cpp src/fesa/io/abaqus/semantic_mapper.cpp src/fesa/model/domain.cpp diff --git a/include/fesa/fem/gauss_rule.hpp b/include/fesa/fem/gauss_rule.hpp new file mode 100644 index 0000000..f33d45b --- /dev/null +++ b/include/fesa/fem/gauss_rule.hpp @@ -0,0 +1,14 @@ +#pragma once + +#include + +namespace fesa { + +struct GaussPoint1D final { + double xi; + double weight; +}; + +[[nodiscard]] std::span gauss_rule_1d(int order); + +} // namespace fesa diff --git a/include/fesa/fem/line2_shape.hpp b/include/fesa/fem/line2_shape.hpp new file mode 100644 index 0000000..a95dd1f --- /dev/null +++ b/include/fesa/fem/line2_shape.hpp @@ -0,0 +1,13 @@ +#pragma once + +#include + +namespace fesa { + +[[nodiscard]] std::array line2_shape(double xi); + +[[nodiscard]] std::array line2_shape_derivative(double xi); + +[[nodiscard]] double line2_jacobian(double length); + +} // namespace fesa diff --git a/src/fesa/fem/gauss_rule.cpp b/src/fesa/fem/gauss_rule.cpp new file mode 100644 index 0000000..cdf4876 --- /dev/null +++ b/src/fesa/fem/gauss_rule.cpp @@ -0,0 +1,33 @@ +#include + +#include +#include + +namespace fesa { + +namespace { + +constexpr std::array order_one{{ + {0.0, 2.0}, +}}; + +constexpr double inverse_sqrt_three = 0.57735026918962576451; +constexpr std::array order_two{{ + {-inverse_sqrt_three, 1.0}, + {inverse_sqrt_three, 1.0}, +}}; + +} // namespace + +std::span gauss_rule_1d(const int order) { + switch (order) { + case 1: + return order_one; + case 2: + return order_two; + default: + throw std::invalid_argument{"Gauss rule order must be 1 or 2."}; + } +} + +} // namespace fesa diff --git a/src/fesa/fem/line2_shape.cpp b/src/fesa/fem/line2_shape.cpp new file mode 100644 index 0000000..960a12a --- /dev/null +++ b/src/fesa/fem/line2_shape.cpp @@ -0,0 +1,31 @@ +#include + +#include +#include + +namespace fesa { + +std::array line2_shape(const double xi) { + return {0.5 * (1.0 - xi), 0.5 * (1.0 + xi)}; +} + +std::array line2_shape_derivative(const double) { + return {-0.5, 0.5}; +} + +double line2_jacobian(const double length) { + if (!std::isfinite(length) || length <= 0.0) { + throw std::invalid_argument{ + "Line2 Jacobian requires a positive finite length."}; + } + + const double jacobian = length / 2.0; + if (jacobian == 0.0) { + throw std::invalid_argument{ + "Line2 Jacobian must be representable as a positive double."}; + } + + return jacobian; +} + +} // namespace fesa diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 519a78a..1edf7ce 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -170,3 +170,41 @@ add_test( COMMAND "$" --gtest_filter=ActiveInstance.* ) + +add_executable(fesa_fem_primitives_tests + unit/fem/gauss_rule_test.cpp + unit/fem/line2_shape_test.cpp +) + +target_compile_features(fesa_fem_primitives_tests PRIVATE cxx_std_20) +target_compile_options( + fesa_fem_primitives_tests + PRIVATE + /W4 + /permissive- + /EHsc +) + +target_link_libraries(fesa_fem_primitives_tests + PRIVATE + fesa_core + GTest::gtest_main +) + +add_test( + NAME Quadrature + COMMAND "$" + --gtest_filter=Quadrature.* +) + +add_test( + NAME ShapeFunction + COMMAND "$" + --gtest_filter=ShapeFunction.* +) + +add_test( + NAME Jacobian + COMMAND "$" + --gtest_filter=Jacobian.* +) diff --git a/tests/unit/fem/gauss_rule_test.cpp b/tests/unit/fem/gauss_rule_test.cpp new file mode 100644 index 0000000..3370112 --- /dev/null +++ b/tests/unit/fem/gauss_rule_test.cpp @@ -0,0 +1,53 @@ +#include + +#include +#include +#include + +#include + +namespace { + +TEST(Quadrature, OnePointRuleIntegratesDegreeZeroAndOneExactly) { + const auto rule = fesa::gauss_rule_1d(1); + + ASSERT_EQ(rule.size(), 1); + EXPECT_DOUBLE_EQ(rule[0].weight, 2.0); + EXPECT_DOUBLE_EQ(rule[0].weight * rule[0].xi, 0.0); +} + +TEST(Quadrature, TwoPointRuleIntegratesThroughDegreeThreeWithinRoundoff) { + const auto rule = fesa::gauss_rule_1d(2); + + ASSERT_EQ(rule.size(), 2); + + double degree_zero = 0.0; + double degree_one = 0.0; + double degree_two = 0.0; + double degree_three = 0.0; + for (const auto& point : rule) { + const double xi_squared = point.xi * point.xi; + degree_zero += point.weight; + degree_one += point.weight * point.xi; + degree_two += point.weight * xi_squared; + degree_three += point.weight * xi_squared * point.xi; + } + + // Sixteen ulps at unit scale cover only the rounding in these short sums. + constexpr double tolerance = + 16.0 * std::numeric_limits::epsilon(); + EXPECT_NEAR(degree_zero, 2.0, tolerance); + EXPECT_NEAR(degree_one, 0.0, tolerance); + EXPECT_NEAR(degree_two, 2.0 / 3.0, tolerance); + EXPECT_NEAR(degree_three, 0.0, tolerance); +} + +TEST(Quadrature, RejectsUnsupportedOrders) { + for (const int order : std::array{0, 3, -1}) { + EXPECT_THROW( + static_cast(fesa::gauss_rule_1d(order)), + std::invalid_argument); + } +} + +} // namespace diff --git a/tests/unit/fem/line2_shape_test.cpp b/tests/unit/fem/line2_shape_test.cpp new file mode 100644 index 0000000..30df276 --- /dev/null +++ b/tests/unit/fem/line2_shape_test.cpp @@ -0,0 +1,71 @@ +#include + +#include +#include +#include + +#include + +namespace { + +TEST(ShapeFunction, InterpolatesBothEndpoints) { + const auto first_endpoint = fesa::line2_shape(-1.0); + const auto second_endpoint = fesa::line2_shape(1.0); + + EXPECT_DOUBLE_EQ(first_endpoint[0], 1.0); + EXPECT_DOUBLE_EQ(first_endpoint[1], 0.0); + EXPECT_DOUBLE_EQ(second_endpoint[0], 0.0); + EXPECT_DOUBLE_EQ(second_endpoint[1], 1.0); +} + +TEST(ShapeFunction, FormsPartitionOfUnityWithinRoundoff) { + for (const double xi : std::array{-1.0, -0.25, 0.0, 0.4, 1.0}) { + const auto shape = fesa::line2_shape(xi); + + // Four ulps at unit scale cover the two affine evaluations and sum. + constexpr double tolerance = + 4.0 * std::numeric_limits::epsilon(); + EXPECT_NEAR(shape[0] + shape[1], 1.0, tolerance) << "xi=" << xi; + } +} + +TEST(ShapeFunction, DerivativesSumToZeroAtEveryNaturalCoordinate) { + for (const double xi : std::array{-1.0, -0.25, 0.0, 0.4, 1.0}) { + const auto derivative = fesa::line2_shape_derivative(xi); + + EXPECT_DOUBLE_EQ(derivative[0], -0.5); + EXPECT_DOUBLE_EQ(derivative[1], 0.5); + EXPECT_DOUBLE_EQ(derivative[0] + derivative[1], 0.0); + } +} + +TEST(Jacobian, ReturnsHalfThePhysicalLength) { + EXPECT_DOUBLE_EQ(fesa::line2_jacobian(4.0), 2.0); + EXPECT_DOUBLE_EQ(fesa::line2_jacobian(0.25), 0.125); +} + +TEST(Jacobian, DoesNotAdjustSmallPositiveFiniteLengths) { + constexpr double length = 1.0e-300; + + EXPECT_DOUBLE_EQ(fesa::line2_jacobian(length), length / 2.0); +} + +TEST(Jacobian, RejectsNonpositiveAndNonfiniteLengths) { + const double infinity = std::numeric_limits::infinity(); + const double nan = std::numeric_limits::quiet_NaN(); + + for (const double length : std::array{0.0, -1.0, infinity, nan}) { + EXPECT_THROW( + static_cast(fesa::line2_jacobian(length)), + std::invalid_argument); + } +} + +TEST(Jacobian, RejectsLengthWhoseHalfUnderflowsToZero) { + EXPECT_THROW( + static_cast(fesa::line2_jacobian( + std::numeric_limits::denorm_min())), + std::invalid_argument); +} + +} // namespace