feat(equation-and-linear-solve): step 1 — essential-bc-elimination

CMakeLists.txt

@@ -27,6 +27,7 @@ include(cmake/FesaDependencies.cmake)
 
 add_library(fesa_core STATIC
     src/fesa/assembly/serial_assembler.cpp
+    src/fesa/constraints/essential_bc.cpp
     src/fesa/core/version.cpp
     src/fesa/elements/beam/beam3d2.cpp
     src/fesa/fem/beam_frame.cpp

include/fesa/constraints/essential_bc.hpp

@@ -0,0 +1,36 @@
+#pragma once
+
+#include <cstddef>
+#include <optional>
+#include <span>
+#include <vector>
+
+#include <fesa/assembly/equation_system.hpp>
+#include <fesa/core/diagnostic.hpp>
+#include <fesa/fem/dof_manager.hpp>
+#include <fesa/model/step_definition.hpp>
+
+namespace fesa {
+
+struct ReducedSystem final {
+    SymmetricCsr stiffness;
+    std::vector<double> force;
+    std::vector<std::size_t> free_to_full;
+    std::vector<double> prescribed_full;
+};
+
+struct ConstraintResult final {
+    std::optional<ReducedSystem> reduced_system;
+    std::vector<Diagnostic> diagnostics;
+};
+
+[[nodiscard]] ConstraintResult eliminate_essential_bcs(
+    const EquationSystem& original,
+    const DofManager& dofs,
+    std::span<const PrescribedDof> prescribed);
+
+[[nodiscard]] std::vector<double> recover_reaction(
+    const EquationSystem& original,
+    std::span<const double> full_displacement);
+
+}  // namespace fesa

include/fesa/fem/dof_manager.hpp

@@ -34,6 +34,8 @@ public:
     [[nodiscard]] std::size_t free_equation_count() const noexcept;
     [[nodiscard]] std::optional<std::size_t> equation(
         DofAddress address) const;
+    [[nodiscard]] std::optional<double> prescribed_value(
+        DofAddress address) const;
     [[nodiscard]] std::size_t full_dof(DofAddress address) const;
     [[nodiscard]] std::array<std::size_t, 12> element_full_dofs(
         const BeamElement& element) const;

src/fesa/constraints/essential_bc.cpp

@@ -0,0 +1,283 @@
+#include <fesa/constraints/essential_bc.hpp>
+
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <optional>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace fesa {
+namespace {
+
+std::optional<std::string> validate_system(
+    const EquationSystem& system) {
+    const SymmetricCsr& matrix = system.stiffness;
+    if (system.force.size() != matrix.order) {
+        return "Force size must equal the matrix order.";
+    }
+    if (matrix.order >
+        static_cast<std::size_t>(
+            std::numeric_limits<std::int32_t>::max())) {
+        return "Matrix order exceeds the 32-bit CSR index range.";
+    }
+    if (matrix.row_offsets.size() != matrix.order + 1 ||
+        matrix.row_offsets.empty() ||
+        matrix.row_offsets.front() != 0) {
+        return "CSR row offsets must contain order + 1 entries "
+               "starting at zero.";
+    }
+    if (matrix.column_indices.size() != matrix.values.size()) {
+        return "CSR column and value counts must match.";
+    }
+
+    std::int32_t previous_offset = 0;
+    for (const std::int32_t offset : matrix.row_offsets) {
+        if (offset < previous_offset || offset < 0 ||
+            static_cast<std::size_t>(offset) >
+                matrix.column_indices.size()) {
+            return "CSR row offsets must be nondecreasing and in range.";
+        }
+        previous_offset = offset;
+    }
+    if (static_cast<std::size_t>(matrix.row_offsets.back()) !=
+        matrix.column_indices.size()) {
+        return "The final CSR row offset must equal the entry count.";
+    }
+
+    for (std::size_t row = 0; row < matrix.order; ++row) {
+        std::int32_t previous_column = -1;
+        const std::size_t begin =
+            static_cast<std::size_t>(matrix.row_offsets[row]);
+        const std::size_t end =
+            static_cast<std::size_t>(matrix.row_offsets[row + 1]);
+        for (std::size_t entry = begin; entry < end; ++entry) {
+            const std::int32_t column = matrix.column_indices[entry];
+            if (column < static_cast<std::int32_t>(row) ||
+                column >= static_cast<std::int32_t>(matrix.order) ||
+                column <= previous_column) {
+                return "CSR rows must contain sorted unique "
+                       "upper-triangle columns.";
+            }
+            previous_column = column;
+        }
+    }
+    return std::nullopt;
+}
+
+Diagnostic equation_error(std::string code, std::string message) {
+    return {
+        DiagnosticStage::equation,
+        Severity::error,
+        std::move(code),
+        std::move(message),
+        std::nullopt,
+    };
+}
+
+std::int32_t csr_index(const std::size_t value) {
+    if (value >
+        static_cast<std::size_t>(
+            std::numeric_limits<std::int32_t>::max())) {
+        throw std::overflow_error{
+            "Reduced system exceeds the 32-bit CSR index range."};
+    }
+    return static_cast<std::int32_t>(value);
+}
+
+}  // namespace
+
+ConstraintResult eliminate_essential_bcs(
+    const EquationSystem& original,
+    const DofManager& dofs,
+    const std::span<const PrescribedDof> prescribed) {
+    ConstraintResult result;
+    if (const auto error = validate_system(original);
+        error.has_value()) {
+        result.diagnostics.push_back(equation_error(
+            "equation.invalid_system", *error));
+        return result;
+    }
+    if (original.stiffness.order != dofs.full_dof_count()) {
+        result.diagnostics.push_back(equation_error(
+            "equation.dof_count_mismatch",
+            "Matrix order must equal the DofManager full DOF count."));
+        return result;
+    }
+
+    const std::size_t full_count = original.stiffness.order;
+    std::vector<bool> constrained(full_count, false);
+    std::vector<double> prescribed_full(full_count, 0.0);
+    for (const PrescribedDof& value : prescribed) {
+        if (value.dof < 1 || value.dof > 6 ||
+            !std::isfinite(value.value)) {
+            result.diagnostics.push_back(equation_error(
+                "equation.invalid_prescribed_dof",
+                "Prescribed DOF data must contain a DOF in [1, 6] "
+                "and a finite value."));
+            continue;
+        }
+
+        const DofAddress address{
+            value.node,
+            static_cast<NodeDof>(value.dof - std::uint8_t{1}),
+        };
+        std::size_t full = 0;
+        std::optional<double> owned_value;
+        try {
+            full = dofs.full_dof(address);
+            owned_value = dofs.prescribed_value(address);
+        } catch (const std::out_of_range&) {
+            result.diagnostics.push_back(equation_error(
+                "equation.invalid_prescribed_dof",
+                "Prescribed DOF references a node absent from the "
+                "DofManager."));
+            continue;
+        }
+
+        if (constrained[full]) {
+            result.diagnostics.push_back(equation_error(
+                "equation.conflicting_prescribed_dof",
+                "A full DOF has more than one prescribed value."));
+            continue;
+        }
+        constrained[full] = true;
+
+        if (!owned_value.has_value() ||
+            *owned_value != value.value) {
+            result.diagnostics.push_back(equation_error(
+                "equation.prescribed_dof_mismatch",
+                "Prescribed DOFs must match the state owned by the "
+                "DofManager."));
+            continue;
+        }
+        prescribed_full[full] = *owned_value;
+    }
+
+    const std::size_t constrained_count =
+        full_count - dofs.free_equation_count();
+    if (prescribed.size() != constrained_count) {
+        result.diagnostics.push_back(equation_error(
+            "equation.prescribed_dof_mismatch",
+            "Prescribed DOF count must match the DofManager "
+            "constraint state."));
+    }
+    if (!result.diagnostics.empty()) {
+        return result;
+    }
+
+    std::vector<std::size_t> full_to_free(full_count, full_count);
+    std::vector<std::size_t> free_to_full;
+    free_to_full.reserve(dofs.free_equation_count());
+    for (std::size_t full = 0; full < full_count; ++full) {
+        if (!constrained[full]) {
+            full_to_free[full] = free_to_full.size();
+            free_to_full.push_back(full);
+        }
+    }
+    if (free_to_full.size() != dofs.free_equation_count()) {
+        result.diagnostics.push_back(equation_error(
+            "equation.prescribed_dof_mismatch",
+            "Prescribed DOFs do not identify the DofManager constraints."));
+        return result;
+    }
+
+    std::vector<double> reduced_force(
+        free_to_full.size(), 0.0);
+    for (std::size_t free = 0; free < free_to_full.size(); ++free) {
+        reduced_force[free] = original.force[free_to_full[free]];
+    }
+
+    SymmetricCsr reduced{
+        free_to_full.size(),
+        std::vector<std::int32_t>(free_to_full.size() + 1, 0),
+        {},
+        {},
+    };
+    reduced.column_indices.reserve(
+        original.stiffness.column_indices.size());
+    reduced.values.reserve(original.stiffness.values.size());
+
+    for (std::size_t full_row = 0; full_row < full_count; ++full_row) {
+        const std::size_t begin = static_cast<std::size_t>(
+            original.stiffness.row_offsets[full_row]);
+        const std::size_t end = static_cast<std::size_t>(
+            original.stiffness.row_offsets[full_row + 1]);
+        for (std::size_t entry = begin; entry < end; ++entry) {
+            const std::size_t full_column =
+                static_cast<std::size_t>(
+                    original.stiffness.column_indices[entry]);
+            const double stiffness = original.stiffness.values[entry];
+
+            if (constrained[full_row]) {
+                if (!constrained[full_column]) {
+                    reduced_force[full_to_free[full_column]] -=
+                        stiffness * prescribed_full[full_row];
+                }
+            } else if (constrained[full_column]) {
+                reduced_force[full_to_free[full_row]] -=
+                    stiffness * prescribed_full[full_column];
+            } else {
+                reduced.column_indices.push_back(
+                    csr_index(full_to_free[full_column]));
+                reduced.values.push_back(stiffness);
+            }
+        }
+        if (!constrained[full_row]) {
+            reduced.row_offsets[full_to_free[full_row] + 1] =
+                csr_index(reduced.column_indices.size());
+        }
+    }
+
+    result.reduced_system = ReducedSystem{
+        std::move(reduced),
+        std::move(reduced_force),
+        std::move(free_to_full),
+        std::move(prescribed_full),
+    };
+    return result;
+}
+
+std::vector<double> recover_reaction(
+    const EquationSystem& original,
+    const std::span<const double> full_displacement) {
+    if (const auto error = validate_system(original);
+        error.has_value()) {
+        throw std::invalid_argument{*error};
+    }
+    if (full_displacement.size() != original.stiffness.order) {
+        throw std::invalid_argument{
+            "Full displacement size must equal the matrix order."};
+    }
+
+    std::vector<double> reaction = original.force;
+    for (double& value : reaction) {
+        value = -value;
+    }
+    for (std::size_t row = 0;
+         row < original.stiffness.order;
+         ++row) {
+        const std::size_t begin = static_cast<std::size_t>(
+            original.stiffness.row_offsets[row]);
+        const std::size_t end = static_cast<std::size_t>(
+            original.stiffness.row_offsets[row + 1]);
+        for (std::size_t entry = begin; entry < end; ++entry) {
+            const std::size_t column =
+                static_cast<std::size_t>(
+                    original.stiffness.column_indices[entry]);
+            const double stiffness = original.stiffness.values[entry];
+            reaction[row] +=
+                stiffness * full_displacement[column];
+            if (column != row) {
+                reaction[column] +=
+                    stiffness * full_displacement[row];
+            }
+        }
+    }
+    return reaction;
+}
+
+}  // namespace fesa

src/fesa/fem/dof_manager.cpp

@@ -63,6 +63,15 @@ std::optional<std::size_t> DofManager::equation(
     return equations_[full_dof(address)];
 }
 
+std::optional<double> DofManager::prescribed_value(
+    const DofAddress address) const {
+    const std::size_t full = full_dof(address);
+    if (equations_[full].has_value()) {
+        return std::nullopt;
+    }
+    return prescribed_values_[full];
+}
+
 std::array<std::size_t, 12> DofManager::element_full_dofs(
     const BeamElement& element) const {
     std::array<std::size_t, 12> full_dofs{};

tests/CMakeLists.txt

@@ -308,3 +308,40 @@ add_test(
     COMMAND "$<TARGET_FILE:fesa_serial_assembly_tests>"
         --gtest_filter=SymmetricCsr.*
 )
+
+add_executable(fesa_constraint_tests
+    unit/constraints/essential_bc_test.cpp
+)
+
+target_compile_features(fesa_constraint_tests PRIVATE cxx_std_20)
+target_compile_options(
+    fesa_constraint_tests
+    PRIVATE
+        /W4
+        /permissive-
+        /EHsc
+)
+
+target_link_libraries(fesa_constraint_tests
+    PRIVATE
+        fesa_core
+        GTest::gtest_main
+)
+
+add_test(
+    NAME EssentialBc
+    COMMAND "$<TARGET_FILE:fesa_constraint_tests>"
+        --gtest_filter=EssentialBc.*
+)
+
+add_test(
+    NAME ConstraintElimination
+    COMMAND "$<TARGET_FILE:fesa_constraint_tests>"
+        --gtest_filter=ConstraintElimination.*
+)
+
+add_test(
+    NAME Reaction
+    COMMAND "$<TARGET_FILE:fesa_constraint_tests>"
+        --gtest_filter=Reaction.*
+)

tests/unit/constraints/essential_bc_test.cpp

@@ -0,0 +1,178 @@
+#include <fesa/constraints/essential_bc.hpp>
+
+#include <algorithm>
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+#include <fesa/model/domain_builder.hpp>
+
+#include <gtest/gtest.h>
+
+namespace {
+
+fesa::DofManager build_dofs(
+    std::vector<fesa::PrescribedDof> prescribed) {
+    fesa::DomainBuilder builder;
+    builder.add_node({
+        fesa::NodeId{0},
+        fesa::EntityOrigin{"", "", 1},
+        fesa::Vec3{0.0, 0.0, 0.0},
+    });
+    builder.set_step({"Load", std::move(prescribed), {}});
+    auto result = std::move(builder).build();
+    if (!result.domain.has_value()) {
+        throw std::runtime_error{"Test Domain failed validation."};
+    }
+    return fesa::DofManager::build(*result.domain);
+}
+
+fesa::EquationSystem hand_system() {
+    return {
+        {
+            6,
+            {0, 3, 5, 6, 7, 8, 9},
+            {0, 1, 2, 1, 2, 2, 3, 4, 5},
+            {4.0, 1.0, 2.0, 3.0, 6.0, 5.0, 1.0, 1.0, 1.0},
+        },
+        {7.0, 8.0, 9.0, 0.0, 0.0, 0.0},
+    };
+}
+
+TEST(EssentialBc, ZeroPrescribedValueDoesNotShiftReducedForce) {
+    const std::vector<fesa::PrescribedDof> prescribed{
+        {fesa::NodeId{0}, 2, 0.0},
+    };
+    const fesa::DofManager dofs = build_dofs(prescribed);
+
+    const fesa::ConstraintResult result =
+        fesa::eliminate_essential_bcs(hand_system(), dofs, prescribed);
+
+    ASSERT_TRUE(result.reduced_system.has_value());
+    EXPECT_TRUE(result.diagnostics.empty());
+    EXPECT_EQ(
+        result.reduced_system->force,
+        (std::vector<double>{7.0, 9.0, 0.0, 0.0, 0.0}));
+}
+
+TEST(
+    ConstraintElimination,
+    ShiftsNonzeroPrescribedValueAndPreservesOriginalSystem) {
+    const std::vector<fesa::PrescribedDof> prescribed{
+        {fesa::NodeId{0}, 2, 2.0},
+    };
+    const fesa::DofManager dofs = build_dofs(prescribed);
+    const fesa::EquationSystem original = hand_system();
+    const fesa::EquationSystem before = original;
+
+    const fesa::ConstraintResult result =
+        fesa::eliminate_essential_bcs(original, dofs, prescribed);
+
+    ASSERT_TRUE(result.reduced_system.has_value());
+    EXPECT_TRUE(result.diagnostics.empty());
+    const fesa::ReducedSystem& reduced = *result.reduced_system;
+    EXPECT_EQ(reduced.stiffness.order, 5);
+    EXPECT_EQ(
+        reduced.stiffness.row_offsets,
+        (std::vector<std::int32_t>{0, 2, 3, 4, 5, 6}));
+    EXPECT_EQ(
+        reduced.stiffness.column_indices,
+        (std::vector<std::int32_t>{0, 1, 1, 2, 3, 4}));
+    EXPECT_EQ(
+        reduced.stiffness.values,
+        (std::vector<double>{4.0, 2.0, 5.0, 1.0, 1.0, 1.0}));
+    EXPECT_EQ(
+        reduced.force,
+        (std::vector<double>{5.0, -3.0, 0.0, 0.0, 0.0}));
+    EXPECT_EQ(
+        reduced.free_to_full,
+        (std::vector<std::size_t>{0, 2, 3, 4, 5}));
+    EXPECT_EQ(
+        reduced.prescribed_full,
+        (std::vector<double>{0.0, 2.0, 0.0, 0.0, 0.0, 0.0}));
+
+    EXPECT_EQ(original.stiffness.order, before.stiffness.order);
+    EXPECT_EQ(
+        original.stiffness.row_offsets,
+        before.stiffness.row_offsets);
+    EXPECT_EQ(
+        original.stiffness.column_indices,
+        before.stiffness.column_indices);
+    EXPECT_EQ(original.stiffness.values, before.stiffness.values);
+    EXPECT_EQ(original.force, before.force);
+
+    EXPECT_EQ(
+        dofs.reconstruct_full(
+            std::vector<double>{1.0, 3.0, 0.0, 0.0, 0.0}),
+        (std::vector<double>{1.0, 2.0, 3.0, 0.0, 0.0, 0.0}));
+}
+
+TEST(ConstraintElimination, AcceptsAllDofsConstrained) {
+    const std::vector<fesa::PrescribedDof> prescribed{
+        {fesa::NodeId{0}, 1, 1.0},
+        {fesa::NodeId{0}, 2, 2.0},
+        {fesa::NodeId{0}, 3, 3.0},
+        {fesa::NodeId{0}, 4, 4.0},
+        {fesa::NodeId{0}, 5, 5.0},
+        {fesa::NodeId{0}, 6, 6.0},
+    };
+    const fesa::DofManager dofs = build_dofs(prescribed);
+
+    const fesa::ConstraintResult result =
+        fesa::eliminate_essential_bcs(hand_system(), dofs, prescribed);
+
+    ASSERT_TRUE(result.reduced_system.has_value());
+    const fesa::ReducedSystem& reduced = *result.reduced_system;
+    EXPECT_EQ(reduced.stiffness.order, 0);
+    EXPECT_EQ(
+        reduced.stiffness.row_offsets,
+        (std::vector<std::int32_t>{0}));
+    EXPECT_TRUE(reduced.stiffness.column_indices.empty());
+    EXPECT_TRUE(reduced.stiffness.values.empty());
+    EXPECT_TRUE(reduced.force.empty());
+    EXPECT_TRUE(reduced.free_to_full.empty());
+    EXPECT_EQ(
+        reduced.prescribed_full,
+        (std::vector<double>{1.0, 2.0, 3.0, 4.0, 5.0, 6.0}));
+    EXPECT_EQ(
+        dofs.reconstruct_full({}),
+        reduced.prescribed_full);
+}
+
+TEST(ConstraintElimination, RejectsConflictingPrescribedValues) {
+    const fesa::DofManager dofs = build_dofs({
+        {fesa::NodeId{0}, 2, 2.0},
+    });
+    const std::vector<fesa::PrescribedDof> conflicting{
+        {fesa::NodeId{0}, 2, 2.0},
+        {fesa::NodeId{0}, 2, 3.0},
+    };
+
+    const fesa::ConstraintResult result =
+        fesa::eliminate_essential_bcs(
+            hand_system(), dofs, conflicting);
+
+    EXPECT_FALSE(result.reduced_system.has_value());
+    ASSERT_FALSE(result.diagnostics.empty());
+    EXPECT_NE(
+        std::ranges::find(
+            result.diagnostics,
+            "equation.conflicting_prescribed_dof",
+            &fesa::Diagnostic::code),
+        result.diagnostics.end());
+}
+
+TEST(Reaction, UsesOriginalFullEquilibriumEquation) {
+    const std::vector<double> reaction = fesa::recover_reaction(
+        hand_system(),
+        std::vector<double>{1.0, 2.0, 3.0, 0.0, 0.0, 0.0});
+
+    EXPECT_EQ(
+        reaction,
+        (std::vector<double>{5.0, 17.0, 20.0, 0.0, 0.0, 0.0}));
+}
+
+}  // namespace