fix(equation-and-linear-solve): address review findings

.gitignore

@@ -2,6 +2,5 @@ __pycache__/
 .pytest_cache/
 *.py[cod]
 .harness/build/
-*pytest-basetemp*/
 out/
 .worktrees/

include/fesa/constraints/essential_bc.hpp

@@ -1,6 +1,5 @@
 #pragma once
 
-#include <cstddef>
 #include <optional>
 #include <span>
 #include <vector>
@@ -8,15 +7,12 @@
 #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 {
@@ -26,8 +22,7 @@ struct ConstraintResult final {
 
 [[nodiscard]] ConstraintResult eliminate_essential_bcs(
     const EquationSystem& original,
-    const DofManager& dofs,
+    const DofManager& dofs);
-    std::span<const PrescribedDof> prescribed);
 
 [[nodiscard]] std::vector<double> recover_reaction(
     const EquationSystem& original,

include/fesa/fem/dof_manager.hpp

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

include/fesa/solvers/linear/pardiso_linear_solver.hpp

@@ -10,6 +10,9 @@ class PardisoLinearSolver final : public LinearSolver {
 public:
     PardisoLinearSolver();
     ~PardisoLinearSolver() override;
+    PardisoLinearSolver(const PardisoLinearSolver&) = delete;
+    PardisoLinearSolver& operator=(
+        const PardisoLinearSolver&) = delete;
 
     [[nodiscard]] LinearSolveResult solve(
         const SymmetricCsr& matrix,

src/fesa/assembly/serial_assembler.cpp

@@ -47,7 +47,13 @@ SymmetricCsr build_sparsity_pattern(
     const DofManager& dofs) {
     using Coordinate = std::pair<std::size_t, std::size_t>;
     std::vector<Coordinate> coordinates;
-    coordinates.reserve(domain.beam_elements().size() * 78);
+    coordinates.reserve(
+        dofs.full_dof_count() +
+        domain.beam_elements().size() * 78);
+
+    for (std::size_t full = 0; full < dofs.full_dof_count(); ++full) {
+        coordinates.emplace_back(full, full);
+    }
 
     for (const BeamElement& element : domain.beam_elements()) {
         const std::array<std::size_t, 12> full_dofs =

src/fesa/constraints/essential_bc.cpp

@@ -1,5 +1,6 @@
 #include <fesa/constraints/essential_bc.hpp>
 
+#include <algorithm>
 #include <cmath>
 #include <cstddef>
 #include <cstdint>
@@ -33,6 +34,16 @@ std::optional<std::string> validate_system(
     if (matrix.column_indices.size() != matrix.values.size()) {
         return "CSR column and value counts must match.";
     }
+    if (!std::ranges::all_of(
+            matrix.values, [](const double value) {
+                return std::isfinite(value);
+            }) ||
+        !std::ranges::all_of(
+            system.force, [](const double value) {
+                return std::isfinite(value);
+            })) {
+        return "Matrix and force values must be finite.";
+    }
 
     std::int32_t previous_offset = 0;
     for (const std::int32_t offset : matrix.row_offsets) {
@@ -92,8 +103,7 @@ std::int32_t csr_index(const std::size_t value) {
 
 ConstraintResult eliminate_essential_bcs(
     const EquationSystem& original,
-    const DofManager& dofs,
+    const DofManager& dofs) {
-    const std::span<const PrescribedDof> prescribed) {
     ConstraintResult result;
     if (const auto error = validate_system(original);
         error.has_value()) {
@@ -109,91 +119,31 @@ ConstraintResult eliminate_essential_bcs(
     }
 
     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;
+    std::vector<double> prescribed_full(full_count, 0.0);
-    free_to_full.reserve(dofs.free_equation_count());
     for (std::size_t full = 0; full < full_count; ++full) {
-        if (!constrained[full]) {
+        const std::optional<std::size_t> equation =
-            full_to_free[full] = free_to_full.size();
+            dofs.equation(full);
-            free_to_full.push_back(full);
+        if (equation.has_value()) {
+            full_to_free[full] = *equation;
+        } else {
+            prescribed_full[full] =
+                dofs.prescribed_value(full).value();
         }
     }
-    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);
+        dofs.free_equation_count(), 0.0);
-    for (std::size_t free = 0; free < free_to_full.size(); ++free) {
+    for (std::size_t full = 0; full < full_count; ++full) {
-        reduced_force[free] = original.force[free_to_full[free]];
+        if (full_to_free[full] != full_count) {
+            reduced_force[full_to_free[full]] = original.force[full];
+        }
     }
 
     SymmetricCsr reduced{
-        free_to_full.size(),
+        dofs.free_equation_count(),
-        std::vector<std::int32_t>(free_to_full.size() + 1, 0),
+        std::vector<std::int32_t>(
+            dofs.free_equation_count() + 1, 0),
         {},
         {},
     };
@@ -212,12 +162,12 @@ ConstraintResult eliminate_essential_bcs(
                     original.stiffness.column_indices[entry]);
             const double stiffness = original.stiffness.values[entry];
 
-            if (constrained[full_row]) {
+            if (full_to_free[full_row] == full_count) {
-                if (!constrained[full_column]) {
+                if (full_to_free[full_column] != full_count) {
                     reduced_force[full_to_free[full_column]] -=
                         stiffness * prescribed_full[full_row];
                 }
-            } else if (constrained[full_column]) {
+            } else if (full_to_free[full_column] == full_count) {
                 reduced_force[full_to_free[full_row]] -=
                     stiffness * prescribed_full[full_column];
             } else {
@@ -226,17 +176,24 @@ ConstraintResult eliminate_essential_bcs(
                 reduced.values.push_back(stiffness);
             }
         }
-        if (!constrained[full_row]) {
+        if (full_to_free[full_row] != full_count) {
             reduced.row_offsets[full_to_free[full_row] + 1] =
                 csr_index(reduced.column_indices.size());
         }
     }
+    if (!std::ranges::all_of(
+            reduced_force, [](const double value) {
+                return std::isfinite(value);
+            })) {
+        result.diagnostics.push_back(equation_error(
+            "equation.nonfinite_result",
+            "Essential BC elimination produced a nonfinite force."));
+        return result;
+    }
 
     result.reduced_system = ReducedSystem{
         std::move(reduced),
         std::move(reduced_force),
-        std::move(free_to_full),
-        std::move(prescribed_full),
     };
     return result;
 }
@@ -252,6 +209,13 @@ std::vector<double> recover_reaction(
         throw std::invalid_argument{
             "Full displacement size must equal the matrix order."};
     }
+    if (!std::ranges::all_of(
+            full_displacement, [](const double value) {
+                return std::isfinite(value);
+            })) {
+        throw std::invalid_argument{
+            "Full displacement values must be finite."};
+    }
 
     std::vector<double> reaction = original.force;
     for (double& value : reaction) {
@@ -277,6 +241,13 @@ std::vector<double> recover_reaction(
             }
         }
     }
+    if (!std::ranges::all_of(
+            reaction, [](const double value) {
+                return std::isfinite(value);
+            })) {
+        throw std::invalid_argument{
+            "Reaction recovery produced a nonfinite value."};
+    }
     return reaction;
 }
 

src/fesa/fem/dof_manager.cpp

@@ -60,16 +60,25 @@ std::size_t DofManager::free_equation_count() const noexcept {
 
 std::optional<std::size_t> DofManager::equation(
     const DofAddress address) const {
-    return equations_[full_dof(address)];
+    return equation(full_dof(address));
+}
+
+std::optional<std::size_t> DofManager::equation(
+    const std::size_t full_dof) const {
+    return equations_.at(full_dof);
 }
 
 std::optional<double> DofManager::prescribed_value(
     const DofAddress address) const {
-    const std::size_t full = full_dof(address);
+    return prescribed_value(full_dof(address));
-    if (equations_[full].has_value()) {
+}
+
+std::optional<double> DofManager::prescribed_value(
+    const std::size_t full_dof) const {
+    if (equations_.at(full_dof).has_value()) {
         return std::nullopt;
     }
-    return prescribed_values_[full];
+    return prescribed_values_.at(full_dof);
 }
 
 std::array<std::size_t, 12> DofManager::element_full_dofs(

src/fesa/solvers/linear/pardiso_linear_solver.cpp

@@ -106,9 +106,14 @@ public:
     PardisoSession& operator=(const PardisoSession&) = delete;
 
     ~PardisoSession() noexcept {
+        static_cast<void>(release());
+    }
+
+    MKL_INT release() noexcept {
         if (!active_) {
-            return;
+            return 0;
         }
+        active_ = false;
 
         constexpr MKL_INT release_all = -1;
         MKL_INT error = 0;
@@ -129,6 +134,7 @@ public:
             right_hand_side_->data(),
             solution_->data(),
             &error);
+        return error;
     }
 
     MKL_INT execute(
@@ -260,37 +266,45 @@ LinearSolveResult PardisoLinearSolver::solve(
         PardisoSession session;
 
         constexpr MKL_INT analyze = 11;
-        if (const MKL_INT error =
+        MKL_INT error = session.execute(
-                session.execute(
+            analyze, matrix, right_hand_side, solution);
-                    analyze, matrix, right_hand_side, solution);
+        if (error != 0) {
-            error != 0) {
             result.diagnostics.push_back(solver_error(
                 "solver.analysis_failed",
                 pardiso_failure("analysis", error)));
-            return result;
         }
 
-        constexpr MKL_INT factorize = 22;
+        if (error == 0) {
-        if (const MKL_INT error =
+            constexpr MKL_INT factorize = 22;
-                session.execute(
+            error = session.execute(
-                    factorize, matrix, right_hand_side, solution);
+                factorize, matrix, right_hand_side, solution);
-            error != 0) {
+            if (error != 0) {
-            result.diagnostics.push_back(solver_error(
+                result.diagnostics.push_back(solver_error(
-                "solver.factorization_failed",
+                    "solver.factorization_failed",
-                pardiso_failure("factorization", error)));
+                    pardiso_failure("factorization", error)));
-            return result;
+            }
         }
 
-        constexpr MKL_INT solve_system = 33;
+        if (error == 0) {
-        if (const MKL_INT error =
+            constexpr MKL_INT solve_system = 33;
-                session.execute(
+            error = session.execute(
-                    solve_system, matrix, right_hand_side, solution);
+                solve_system, matrix, right_hand_side, solution);
-            error != 0) {
+            if (error != 0) {
-            result.diagnostics.push_back(solver_error(
+                result.diagnostics.push_back(solver_error(
-                "solver.solve_failed",
+                    "solver.solve_failed",
-                pardiso_failure("solve", error)));
+                    pardiso_failure("solve", error)));
-            return result;
+            }
         }
+
+        if (const MKL_INT release_error = session.release();
+            release_error != 0) {
+            result.diagnostics.push_back(solver_error(
+                "solver.release_failed",
+                pardiso_failure("release", release_error)));
+        }
+    }
+    if (!result.diagnostics.empty()) {
+        return result;
     }
     if (!std::ranges::all_of(solution, [](const double value) {
             return std::isfinite(value);

tests/unit/assembly/serial_assembler_test.cpp

@@ -2,6 +2,7 @@
 
 #include <algorithm>
 #include <array>
+#include <bit>
 #include <cstddef>
 #include <cstdint>
 #include <limits>
@@ -158,6 +159,78 @@ fesa::Domain build_parallel_domain(
     return finish_domain(std::move(builder), {"Load", {}, {}});
 }
 
+fesa::Domain build_rounding_domain(const bool large_element_first) {
+    fesa::DomainBuilder builder;
+    builder.add_node({
+        fesa::NodeId{4},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 1},
+        fesa::Vec3{0.0, 0.0, 0.0},
+    });
+    builder.add_node({
+        fesa::NodeId{10},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 2},
+        fesa::Vec3{1.0, 0.0, 0.0},
+    });
+
+    for (std::size_t index = 0; index < 3; ++index) {
+        fesa::IsotropicElastic material = unit_material();
+        material.id =
+            fesa::MaterialId{static_cast<std::int64_t>(index)};
+        material.name = "Material-" + std::to_string(index);
+        material.young = index == 2 ? 1.0e16 : 1.0;
+        builder.add_material(std::move(material));
+    }
+    builder.add_section(unit_section());
+
+    std::array<std::size_t, 3> storage_order{0, 1, 2};
+    if (large_element_first) {
+        storage_order = {2, 0, 1};
+    }
+    for (const std::size_t index : storage_order) {
+        builder.add_beam_element({
+            fesa::ElementId{static_cast<std::int64_t>(index)},
+            fesa::EntityOrigin{
+                "BeamPart",
+                "Beam-1",
+                static_cast<std::int64_t>((index + 1) * 10),
+            },
+            {fesa::NodeId{4}, fesa::NodeId{10}},
+            fesa::MaterialId{static_cast<std::int64_t>(index)},
+            fesa::SectionId{0},
+        });
+    }
+    return finish_domain(std::move(builder), {"Load", {}, {}});
+}
+
+fesa::Domain build_orphan_node_domain() {
+    fesa::DomainBuilder builder;
+    builder.add_node({
+        fesa::NodeId{0},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 1},
+        fesa::Vec3{0.0, 0.0, 0.0},
+    });
+    builder.add_node({
+        fesa::NodeId{1},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 2},
+        fesa::Vec3{1.0, 0.0, 0.0},
+    });
+    builder.add_node({
+        fesa::NodeId{2},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 3},
+        fesa::Vec3{2.0, 0.0, 0.0},
+    });
+    builder.add_material(unit_material());
+    builder.add_section(unit_section());
+    builder.add_beam_element({
+        fesa::ElementId{0},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 1},
+        {fesa::NodeId{0}, fesa::NodeId{1}},
+        fesa::MaterialId{0},
+        fesa::SectionId{0},
+    });
+    return finish_domain(std::move(builder), {"Load", {}, {}});
+}
+
 double csr_value(
     const fesa::SymmetricCsr& matrix,
     std::size_t row,
@@ -246,6 +319,29 @@ TEST(SymmetricCsr, StoresSortedUpperTriangleWithValidOffsets) {
     }
 }
 
+TEST(SparsePattern, IncludesZeroDiagonalForEveryOrphanNodeDof) {
+    const fesa::Domain domain = build_orphan_node_domain();
+    const fesa::SymmetricCsr matrix =
+        fesa::assemble_serial(
+            domain, fesa::DofManager::build(domain))
+            .stiffness;
+
+    ASSERT_EQ(matrix.order, 18);
+    for (std::size_t row = 0; row < matrix.order; ++row) {
+        const auto begin =
+            matrix.column_indices.begin() + matrix.row_offsets[row];
+        const auto end =
+            matrix.column_indices.begin() + matrix.row_offsets[row + 1];
+        const auto diagonal = std::lower_bound(
+            begin, end, static_cast<std::int32_t>(row));
+        ASSERT_NE(diagonal, end);
+        EXPECT_EQ(*diagonal, static_cast<std::int32_t>(row));
+    }
+    for (std::size_t row = 12; row < 18; ++row) {
+        EXPECT_DOUBLE_EQ(csr_value(matrix, row, row), 0.0);
+    }
+}
+
 TEST(SerialAssembly, AssemblesHandCalculatedAxialChainAndFullLoad) {
     const fesa::Domain domain = build_chain_domain(false);
     const fesa::DofManager dofs = fesa::DofManager::build(domain);
@@ -307,6 +403,36 @@ TEST(SerialAssembly, MergesDuplicateElementContributions) {
     }
 }
 
+TEST(SerialAssembly, ReducesByElementOriginNotDomainStorageOrder) {
+    const fesa::Domain origin_order_domain =
+        build_rounding_domain(false);
+    const fesa::Domain storage_order_domain =
+        build_rounding_domain(true);
+
+    const double origin_order_value = csr_value(
+        fesa::assemble_serial(
+            origin_order_domain,
+            fesa::DofManager::build(origin_order_domain))
+            .stiffness,
+        0,
+        0);
+    const double storage_order_value = csr_value(
+        fesa::assemble_serial(
+            storage_order_domain,
+            fesa::DofManager::build(storage_order_domain))
+            .stiffness,
+        0,
+        0);
+    const double expected = (1.0 + 1.0) + 1.0e16;
+
+    EXPECT_EQ(
+        std::bit_cast<std::uint64_t>(origin_order_value),
+        std::bit_cast<std::uint64_t>(expected));
+    EXPECT_EQ(
+        std::bit_cast<std::uint64_t>(storage_order_value),
+        std::bit_cast<std::uint64_t>(expected));
+}
+
 TEST(SerialAssembly, IsIndependentOfDomainStorageAndExternalLabelOrder) {
     const fesa::Domain first_domain = build_chain_domain(false);
     const fesa::Domain second_domain = build_chain_domain(true);

tests/unit/constraints/essential_bc_test.cpp

@@ -4,6 +4,7 @@
 #include <array>
 #include <cstddef>
 #include <cstdint>
+#include <limits>
 #include <stdexcept>
 #include <utility>
 #include <vector>
@@ -49,7 +50,7 @@ TEST(EssentialBc, ZeroPrescribedValueDoesNotShiftReducedForce) {
     const fesa::DofManager dofs = build_dofs(prescribed);
 
     const fesa::ConstraintResult result =
-        fesa::eliminate_essential_bcs(hand_system(), dofs, prescribed);
+        fesa::eliminate_essential_bcs(hand_system(), dofs);
 
     ASSERT_TRUE(result.reduced_system.has_value());
     EXPECT_TRUE(result.diagnostics.empty());
@@ -69,7 +70,7 @@ TEST(
     const fesa::EquationSystem before = original;
 
     const fesa::ConstraintResult result =
-        fesa::eliminate_essential_bcs(original, dofs, prescribed);
+        fesa::eliminate_essential_bcs(original, dofs);
 
     ASSERT_TRUE(result.reduced_system.has_value());
     EXPECT_TRUE(result.diagnostics.empty());
@@ -87,12 +88,6 @@ TEST(
     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(
@@ -122,7 +117,7 @@ TEST(ConstraintElimination, AcceptsAllDofsConstrained) {
     const fesa::DofManager dofs = build_dofs(prescribed);
 
     const fesa::ConstraintResult result =
-        fesa::eliminate_essential_bcs(hand_system(), dofs, prescribed);
+        fesa::eliminate_essential_bcs(hand_system(), dofs);
 
     ASSERT_TRUE(result.reduced_system.has_value());
     const fesa::ReducedSystem& reduced = *result.reduced_system;
@@ -133,36 +128,51 @@ TEST(ConstraintElimination, AcceptsAllDofsConstrained) {
     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);
+        (std::vector<double>{1.0, 2.0, 3.0, 4.0, 5.0, 6.0}));
 }
 
-TEST(ConstraintElimination, RejectsConflictingPrescribedValues) {
+TEST(ConstraintElimination, RejectsNonfiniteSystemData) {
+    const fesa::DofManager dofs = build_dofs({});
+    fesa::EquationSystem nonfinite_matrix = hand_system();
+    nonfinite_matrix.stiffness.values[0] =
+        std::numeric_limits<double>::quiet_NaN();
+    fesa::EquationSystem nonfinite_force = hand_system();
+    nonfinite_force.force[0] =
+        std::numeric_limits<double>::infinity();
+
+    const fesa::ConstraintResult matrix_result =
+        fesa::eliminate_essential_bcs(nonfinite_matrix, dofs);
+    const fesa::ConstraintResult force_result =
+        fesa::eliminate_essential_bcs(nonfinite_force, dofs);
+
+    ASSERT_FALSE(matrix_result.diagnostics.empty());
+    EXPECT_EQ(
+        matrix_result.diagnostics.front().code,
+        "equation.invalid_system");
+    ASSERT_FALSE(force_result.diagnostics.empty());
+    EXPECT_EQ(
+        force_result.diagnostics.front().code,
+        "equation.invalid_system");
+}
+
+TEST(ConstraintElimination, RejectsNonfiniteReducedForce) {
     const fesa::DofManager dofs = build_dofs({
         {fesa::NodeId{0}, 2, 2.0},
     });
-    const std::vector<fesa::PrescribedDof> conflicting{
+    fesa::EquationSystem overflowing = hand_system();
-        {fesa::NodeId{0}, 2, 2.0},
+    overflowing.stiffness.values[1] =
-        {fesa::NodeId{0}, 2, 3.0},
+        std::numeric_limits<double>::max();
-    };
 
     const fesa::ConstraintResult result =
-        fesa::eliminate_essential_bcs(
+        fesa::eliminate_essential_bcs(overflowing, dofs);
-            hand_system(), dofs, conflicting);
 
     EXPECT_FALSE(result.reduced_system.has_value());
     ASSERT_FALSE(result.diagnostics.empty());
-    EXPECT_NE(
+    EXPECT_EQ(
-        std::ranges::find(
+        result.diagnostics.front().code,
-            result.diagnostics,
+        "equation.nonfinite_result");
-            "equation.conflicting_prescribed_dof",
-            &fesa::Diagnostic::code),
-        result.diagnostics.end());
 }
 
 TEST(Reaction, UsesOriginalFullEquilibriumEquation) {
@@ -175,4 +185,26 @@ TEST(Reaction, UsesOriginalFullEquilibriumEquation) {
         (std::vector<double>{5.0, 17.0, 20.0, 0.0, 0.0, 0.0}));
 }
 
+TEST(Reaction, RejectsNonfiniteFullDisplacement) {
+    std::vector<double> displacement(6, 0.0);
+    displacement[0] = std::numeric_limits<double>::quiet_NaN();
+
+    EXPECT_THROW(
+        static_cast<void>(
+            fesa::recover_reaction(hand_system(), displacement)),
+        std::invalid_argument);
+}
+
+TEST(Reaction, RejectsNonfiniteRecoveredValue) {
+    fesa::EquationSystem overflowing = hand_system();
+    overflowing.stiffness.values[0] =
+        std::numeric_limits<double>::max();
+
+    EXPECT_THROW(
+        static_cast<void>(fesa::recover_reaction(
+            overflowing,
+            std::vector<double>{2.0, 0.0, 0.0, 0.0, 0.0, 0.0})),
+        std::invalid_argument);
+}
+
 }  // namespace

tests/unit/fem/dof_manager_test.cpp

@@ -130,6 +130,16 @@ TEST(DofManager, NumbersOnlyFreeDofsAndReconstructsPrescribedValues) {
     EXPECT_EQ(
         dofs.equation({fesa::NodeId{20}, fesa::NodeDof::ux}),
         std::optional<std::size_t>{10});
+    EXPECT_EQ(dofs.equation(std::size_t{0}), std::nullopt);
+    EXPECT_EQ(
+        dofs.prescribed_value(std::size_t{0}),
+        std::optional<double>{1.25});
+    EXPECT_EQ(
+        dofs.equation(std::size_t{1}),
+        std::optional<std::size_t>{0});
+    EXPECT_EQ(
+        dofs.prescribed_value(std::size_t{1}),
+        std::nullopt);
 
     std::vector<double> reduced(dofs.free_equation_count());
     for (std::size_t equation = 0; equation < reduced.size(); ++equation) {
@@ -190,6 +200,13 @@ TEST(DofManager, RejectsInvalidAddressesAndReducedVectorSize) {
         static_cast<void>(
             dofs.equation({fesa::NodeId{99}, fesa::NodeDof::ux})),
         std::out_of_range);
+    EXPECT_THROW(
+        static_cast<void>(dofs.equation(dofs.full_dof_count())),
+        std::out_of_range);
+    EXPECT_THROW(
+        static_cast<void>(
+            dofs.prescribed_value(dofs.full_dof_count())),
+        std::out_of_range);
     EXPECT_THROW(
         static_cast<void>(
             dofs.reconstruct_full(std::vector<double>(17, 0.0))),

tests/unit/solvers/linear/pardiso_linear_solver_test.cpp

@@ -6,12 +6,18 @@
 #include <cstdint>
 #include <span>
 #include <string_view>
+#include <type_traits>
 #include <vector>
 
 #include <gtest/gtest.h>
 
 namespace {
 
+static_assert(
+    !std::is_copy_constructible_v<fesa::PardisoLinearSolver>);
+static_assert(
+    !std::is_copy_assignable_v<fesa::PardisoLinearSolver>);
+
 fesa::SymmetricCsr spd_matrix() {
     return {
         3,