feat(fem-and-beam-kernel): step 0 — quadrature-and-shape-functions

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

include/fesa/fem/gauss_rule.hpp

@@ -0,0 +1,14 @@
+#pragma once
+
+#include <span>
+
+namespace fesa {
+
+struct GaussPoint1D final {
+    double xi;
+    double weight;
+};
+
+[[nodiscard]] std::span<const GaussPoint1D> gauss_rule_1d(int order);
+
+}  // namespace fesa

include/fesa/fem/line2_shape.hpp

@@ -0,0 +1,13 @@
+#pragma once
+
+#include <array>
+
+namespace fesa {
+
+[[nodiscard]] std::array<double, 2> line2_shape(double xi);
+
+[[nodiscard]] std::array<double, 2> line2_shape_derivative(double xi);
+
+[[nodiscard]] double line2_jacobian(double length);
+
+}  // namespace fesa

src/fesa/fem/gauss_rule.cpp

@@ -0,0 +1,33 @@
+#include <fesa/fem/gauss_rule.hpp>
+
+#include <array>
+#include <stdexcept>
+
+namespace fesa {
+
+namespace {
+
+constexpr std::array<GaussPoint1D, 1> order_one{{
+    {0.0, 2.0},
+}};
+
+constexpr double inverse_sqrt_three = 0.57735026918962576451;
+constexpr std::array<GaussPoint1D, 2> order_two{{
+    {-inverse_sqrt_three, 1.0},
+    {inverse_sqrt_three, 1.0},
+}};
+
+}  // namespace
+
+std::span<const GaussPoint1D> 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

src/fesa/fem/line2_shape.cpp

@@ -0,0 +1,31 @@
+#include <fesa/fem/line2_shape.hpp>
+
+#include <cmath>
+#include <stdexcept>
+
+namespace fesa {
+
+std::array<double, 2> line2_shape(const double xi) {
+    return {0.5 * (1.0 - xi), 0.5 * (1.0 + xi)};
+}
+
+std::array<double, 2> 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

tests/CMakeLists.txt

@@ -170,3 +170,41 @@ add_test(
     COMMAND "$<TARGET_FILE:fesa_deck_to_domain_tests>"
         --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 "$<TARGET_FILE:fesa_fem_primitives_tests>"
+        --gtest_filter=Quadrature.*
+)
+
+add_test(
+    NAME ShapeFunction
+    COMMAND "$<TARGET_FILE:fesa_fem_primitives_tests>"
+        --gtest_filter=ShapeFunction.*
+)
+
+add_test(
+    NAME Jacobian
+    COMMAND "$<TARGET_FILE:fesa_fem_primitives_tests>"
+        --gtest_filter=Jacobian.*
+)

tests/unit/fem/gauss_rule_test.cpp

@@ -0,0 +1,53 @@
+#include <fesa/fem/gauss_rule.hpp>
+
+#include <array>
+#include <limits>
+#include <stdexcept>
+
+#include <gtest/gtest.h>
+
+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<double>::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<void>(fesa::gauss_rule_1d(order)),
+            std::invalid_argument);
+    }
+}
+
+}  // namespace

tests/unit/fem/line2_shape_test.cpp

@@ -0,0 +1,71 @@
+#include <fesa/fem/line2_shape.hpp>
+
+#include <array>
+#include <limits>
+#include <stdexcept>
+
+#include <gtest/gtest.h>
+
+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<double>::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<double>::infinity();
+    const double nan = std::numeric_limits<double>::quiet_NaN();
+
+    for (const double length : std::array{0.0, -1.0, infinity, nan}) {
+        EXPECT_THROW(
+            static_cast<void>(fesa::line2_jacobian(length)),
+            std::invalid_argument);
+    }
+}
+
+TEST(Jacobian, RejectsLengthWhoseHalfUnderflowsToZero) {
+    EXPECT_THROW(
+        static_cast<void>(fesa::line2_jacobian(
+            std::numeric_limits<double>::denorm_min())),
+        std::invalid_argument);
+}
+
+}  // namespace