From 1bf277cbf186716bff88d5afd1074711e12f671e Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Fri, 31 Jul 2026 16:00:45 +0900 Subject: [PATCH] =?UTF-8?q?feat(equation-and-linear-solve):=20step=201=20?= =?UTF-8?q?=E2=80=94=20essential-bc-elimination?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CMakeLists.txt | 1 + include/fesa/constraints/essential_bc.hpp | 36 +++ include/fesa/fem/dof_manager.hpp | 2 + src/fesa/constraints/essential_bc.cpp | 283 +++++++++++++++++++ src/fesa/fem/dof_manager.cpp | 9 + tests/CMakeLists.txt | 37 +++ tests/unit/constraints/essential_bc_test.cpp | 178 ++++++++++++ 7 files changed, 546 insertions(+) create mode 100644 include/fesa/constraints/essential_bc.hpp create mode 100644 src/fesa/constraints/essential_bc.cpp create mode 100644 tests/unit/constraints/essential_bc_test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index bffc356..a33195d 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/include/fesa/constraints/essential_bc.hpp b/include/fesa/constraints/essential_bc.hpp new file mode 100644 index 0000000..a95c8d8 --- /dev/null +++ b/include/fesa/constraints/essential_bc.hpp @@ -0,0 +1,36 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace fesa { + +struct ReducedSystem final { + SymmetricCsr stiffness; + std::vector force; + std::vector free_to_full; + std::vector prescribed_full; +}; + +struct ConstraintResult final { + std::optional reduced_system; + std::vector diagnostics; +}; + +[[nodiscard]] ConstraintResult eliminate_essential_bcs( + const EquationSystem& original, + const DofManager& dofs, + std::span prescribed); + +[[nodiscard]] std::vector recover_reaction( + const EquationSystem& original, + std::span full_displacement); + +} // namespace fesa diff --git a/include/fesa/fem/dof_manager.hpp b/include/fesa/fem/dof_manager.hpp index a844426..341ae55 100644 --- a/include/fesa/fem/dof_manager.hpp +++ b/include/fesa/fem/dof_manager.hpp @@ -34,6 +34,8 @@ public: [[nodiscard]] std::size_t free_equation_count() const noexcept; [[nodiscard]] std::optional equation( DofAddress address) const; + [[nodiscard]] std::optional prescribed_value( + DofAddress address) const; [[nodiscard]] std::size_t full_dof(DofAddress address) const; [[nodiscard]] std::array element_full_dofs( const BeamElement& element) const; diff --git a/src/fesa/constraints/essential_bc.cpp b/src/fesa/constraints/essential_bc.cpp new file mode 100644 index 0000000..1c327bb --- /dev/null +++ b/src/fesa/constraints/essential_bc.cpp @@ -0,0 +1,283 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace fesa { +namespace { + +std::optional 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::numeric_limits::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(offset) > + matrix.column_indices.size()) { + return "CSR row offsets must be nondecreasing and in range."; + } + previous_offset = offset; + } + if (static_cast(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(matrix.row_offsets[row]); + const std::size_t end = + static_cast(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(row) || + column >= static_cast(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::numeric_limits::max())) { + throw std::overflow_error{ + "Reduced system exceeds the 32-bit CSR index range."}; + } + return static_cast(value); +} + +} // namespace + +ConstraintResult eliminate_essential_bcs( + const EquationSystem& original, + const DofManager& dofs, + const std::span 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 constrained(full_count, false); + std::vector 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(value.dof - std::uint8_t{1}), + }; + std::size_t full = 0; + std::optional 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 full_to_free(full_count, full_count); + std::vector 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 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(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( + original.stiffness.row_offsets[full_row]); + const std::size_t end = static_cast( + original.stiffness.row_offsets[full_row + 1]); + for (std::size_t entry = begin; entry < end; ++entry) { + const std::size_t full_column = + static_cast( + 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 recover_reaction( + const EquationSystem& original, + const std::span 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 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( + original.stiffness.row_offsets[row]); + const std::size_t end = static_cast( + original.stiffness.row_offsets[row + 1]); + for (std::size_t entry = begin; entry < end; ++entry) { + const std::size_t column = + static_cast( + 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 diff --git a/src/fesa/fem/dof_manager.cpp b/src/fesa/fem/dof_manager.cpp index 072fd1a..edebd02 100644 --- a/src/fesa/fem/dof_manager.cpp +++ b/src/fesa/fem/dof_manager.cpp @@ -63,6 +63,15 @@ std::optional DofManager::equation( return equations_[full_dof(address)]; } +std::optional 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 DofManager::element_full_dofs( const BeamElement& element) const { std::array full_dofs{}; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index cc0796a..113a566 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -308,3 +308,40 @@ add_test( COMMAND "$" --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 "$" + --gtest_filter=EssentialBc.* +) + +add_test( + NAME ConstraintElimination + COMMAND "$" + --gtest_filter=ConstraintElimination.* +) + +add_test( + NAME Reaction + COMMAND "$" + --gtest_filter=Reaction.* +) diff --git a/tests/unit/constraints/essential_bc_test.cpp b/tests/unit/constraints/essential_bc_test.cpp new file mode 100644 index 0000000..49674a8 --- /dev/null +++ b/tests/unit/constraints/essential_bc_test.cpp @@ -0,0 +1,178 @@ +#include + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +namespace { + +fesa::DofManager build_dofs( + std::vector 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 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{7.0, 9.0, 0.0, 0.0, 0.0})); +} + +TEST( + ConstraintElimination, + ShiftsNonzeroPrescribedValueAndPreservesOriginalSystem) { + const std::vector 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{0, 2, 3, 4, 5, 6})); + EXPECT_EQ( + reduced.stiffness.column_indices, + (std::vector{0, 1, 1, 2, 3, 4})); + EXPECT_EQ( + reduced.stiffness.values, + (std::vector{4.0, 2.0, 5.0, 1.0, 1.0, 1.0})); + EXPECT_EQ( + reduced.force, + (std::vector{5.0, -3.0, 0.0, 0.0, 0.0})); + EXPECT_EQ( + reduced.free_to_full, + (std::vector{0, 2, 3, 4, 5})); + EXPECT_EQ( + reduced.prescribed_full, + (std::vector{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{1.0, 3.0, 0.0, 0.0, 0.0}), + (std::vector{1.0, 2.0, 3.0, 0.0, 0.0, 0.0})); +} + +TEST(ConstraintElimination, AcceptsAllDofsConstrained) { + const std::vector 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{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{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 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 reaction = fesa::recover_reaction( + hand_system(), + std::vector{1.0, 2.0, 3.0, 0.0, 0.0, 0.0}); + + EXPECT_EQ( + reaction, + (std::vector{5.0, 17.0, 20.0, 0.0, 0.0, 0.0})); +} + +} // namespace