feat(results-and-pipeline): step 2 — linear-static-analysis

CMakeLists.txt

@@ -26,6 +26,7 @@ include(CTest)
 include(cmake/FesaDependencies.cmake)
 
 add_library(fesa_core STATIC
+    src/fesa/analysis/linear_static_analysis.cpp
     src/fesa/assembly/serial_assembler.cpp
     src/fesa/constraints/essential_bc.cpp
     src/fesa/core/version.cpp

include/fesa/analysis/linear_static_analysis.hpp

@@ -0,0 +1,23 @@
+#pragma once
+
+#include <optional>
+#include <vector>
+
+#include <fesa/core/diagnostic.hpp>
+#include <fesa/model/domain.hpp>
+#include <fesa/results/result_database.hpp>
+
+namespace fesa {
+
+struct AnalysisRunResult final {
+    bool succeeded;
+    std::optional<ResultDatabase> results;
+    std::vector<Diagnostic> diagnostics;
+};
+
+class LinearStaticAnalysis final {
+public:
+    [[nodiscard]] AnalysisRunResult run(const Domain& domain) const;
+};
+
+}  // namespace fesa

src/fesa/analysis/linear_static_analysis.cpp

@@ -0,0 +1,137 @@
+#include <fesa/analysis/linear_static_analysis.hpp>
+
+#include <algorithm>
+#include <array>
+#include <cstddef>
+#include <exception>
+#include <optional>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <fesa/assembly/equation_system.hpp>
+#include <fesa/assembly/serial_assembler.hpp>
+#include <fesa/constraints/essential_bc.hpp>
+#include <fesa/fem/dof_manager.hpp>
+#include <fesa/solvers/linear/pardiso_linear_solver.hpp>
+
+namespace fesa {
+namespace {
+
+AnalysisRunResult failure(std::vector<Diagnostic> diagnostics) {
+    return {false, std::nullopt, std::move(diagnostics)};
+}
+
+AnalysisRunResult equation_failure(
+    std::string code,
+    std::string message) {
+    return failure({{
+        DiagnosticStage::equation,
+        Severity::error,
+        std::move(code),
+        std::move(message),
+        std::nullopt,
+    }});
+}
+
+NodalFrame build_nodal_frame(
+    const Domain& domain,
+    const DofManager& dofs,
+    const std::vector<double>& displacement,
+    const std::vector<double>& reaction) {
+    std::vector<NodeId> node_ids;
+    node_ids.reserve(domain.nodes().size());
+    for (const Node& node : domain.nodes()) {
+        node_ids.push_back(node.id);
+    }
+    std::ranges::sort(node_ids);
+
+    NodalFrame nodal;
+    nodal.node_ids = std::move(node_ids);
+    nodal.displacement.reserve(nodal.node_ids.size());
+    nodal.reaction.reserve(nodal.node_ids.size());
+    for (const NodeId node_id : nodal.node_ids) {
+        std::array<double, 6> node_displacement{};
+        std::array<double, 6> node_reaction{};
+        for (std::size_t component = 0; component < 6; ++component) {
+            const std::size_t full_dof = dofs.full_dof({
+                node_id,
+                static_cast<NodeDof>(component),
+            });
+            node_displacement[component] = displacement[full_dof];
+            node_reaction[component] = reaction[full_dof];
+        }
+        nodal.displacement.push_back(node_displacement);
+        nodal.reaction.push_back(node_reaction);
+    }
+    return nodal;
+}
+
+}  // namespace
+
+AnalysisRunResult LinearStaticAnalysis::run(const Domain& domain) const {
+    const DofManager dofs = DofManager::build(domain);
+
+    std::optional<EquationSystem> original;
+    try {
+        original = assemble_serial(domain, dofs);
+    } catch (const std::exception& error) {
+        return equation_failure(
+            "equation.assembly_failed", error.what());
+    }
+
+    ConstraintResult constrained =
+        eliminate_essential_bcs(*original, dofs);
+    if (!constrained.reduced_system.has_value()) {
+        return failure(std::move(constrained.diagnostics));
+    }
+
+    std::vector<double> reduced_solution;
+    if (dofs.free_equation_count() != 0) {
+        PardisoLinearSolver solver;
+        LinearSolveResult solved = solver.solve(
+            constrained.reduced_system->stiffness,
+            constrained.reduced_system->force);
+        if (!solved.diagnostics.empty()) {
+            return failure(std::move(solved.diagnostics));
+        }
+        reduced_solution = std::move(solved.solution);
+    }
+
+    std::vector<double> displacement;
+    try {
+        displacement = dofs.reconstruct_full(reduced_solution);
+    } catch (const std::exception& error) {
+        return equation_failure(
+            "equation.solution_reconstruction_failed", error.what());
+    }
+
+    std::vector<double> reaction;
+    try {
+        reaction = recover_reaction(*original, displacement);
+    } catch (const std::exception& error) {
+        return equation_failure(
+            "equation.reaction_recovery_failed", error.what());
+    }
+
+    ResultDatabase database{
+        "1.0.0",
+        {{
+            domain.step().name,
+            {{
+                1.0,
+                build_nodal_frame(
+                    domain, dofs, displacement, reaction),
+                {},
+            }},
+        }},
+    };
+    Status status = validate_result_database(database);
+    if (!status.succeeded) {
+        return failure(std::move(status.diagnostics));
+    }
+
+    return {true, std::move(database), {}};
+}
+
+}  // namespace fesa

tests/CMakeLists.txt

@@ -449,3 +449,42 @@ set_property(
     PROPERTY ENVIRONMENT_MODIFICATION
         ${FESA_DEPENDENCY_RUNTIME_MODIFICATIONS}
 )
+
+add_executable(fesa_linear_static_analysis_tests
+    unit/analysis/linear_static_analysis_test.cpp
+)
+
+target_compile_features(
+    fesa_linear_static_analysis_tests PRIVATE cxx_std_20
+)
+target_compile_options(
+    fesa_linear_static_analysis_tests
+    PRIVATE
+        /W4
+        /permissive-
+        /EHsc
+)
+
+target_link_libraries(fesa_linear_static_analysis_tests
+    PRIVATE
+        fesa_core
+        GTest::gtest_main
+)
+
+add_test(
+    NAME LinearStaticAnalysis
+    COMMAND "$<TARGET_FILE:fesa_linear_static_analysis_tests>"
+        --gtest_filter=LinearStaticAnalysis.*
+)
+
+add_test(
+    NAME StaticEquilibrium
+    COMMAND "$<TARGET_FILE:fesa_linear_static_analysis_tests>"
+        --gtest_filter=StaticEquilibrium.*
+)
+
+set_property(
+    TEST LinearStaticAnalysis StaticEquilibrium
+    PROPERTY ENVIRONMENT_MODIFICATION
+        ${FESA_DEPENDENCY_RUNTIME_MODIFICATIONS}
+)

tests/unit/analysis/linear_static_analysis_test.cpp

@@ -0,0 +1,217 @@
+#include <fesa/analysis/linear_static_analysis.hpp>
+
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+#include <fesa/assembly/serial_assembler.hpp>
+#include <fesa/fem/dof_manager.hpp>
+#include <fesa/model/domain_builder.hpp>
+#include <fesa/results/result_database.hpp>
+
+#include <gtest/gtest.h>
+
+namespace {
+
+constexpr std::size_t dofs_per_node = 6;
+
+std::vector<fesa::PrescribedDof> fixed_dofs(
+    const fesa::NodeId node,
+    const double ux = 0.0) {
+    std::vector<fesa::PrescribedDof> prescribed;
+    prescribed.reserve(dofs_per_node);
+    for (std::uint8_t dof = 1; dof <= dofs_per_node; ++dof) {
+        prescribed.push_back({node, dof, dof == 1 ? ux : 0.0});
+    }
+    return prescribed;
+}
+
+fesa::Domain build_axial_domain(const bool all_constrained) {
+    fesa::DomainBuilder builder;
+    builder.add_node({
+        fesa::NodeId{20},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 1},
+        fesa::Vec3{0.0, 0.0, 0.0},
+    });
+    builder.add_node({
+        fesa::NodeId{4},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 2},
+        fesa::Vec3{4.0, 0.0, 0.0},
+    });
+    builder.add_material({
+        fesa::MaterialId{0},
+        "HandCalculated",
+        100.0,
+        0.25,
+    });
+    builder.add_section({
+        fesa::SectionId{0},
+        "General",
+        2.0,
+        1.0,
+        1.0,
+        1.0,
+        2.0,
+        2.0,
+        fesa::ShearPropertySource::input,
+        fesa::Vec3{0.0, 1.0, 0.0},
+        {},
+    });
+    builder.add_beam_element({
+        fesa::ElementId{0},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 1},
+        {fesa::NodeId{20}, fesa::NodeId{4}},
+        fesa::MaterialId{0},
+        fesa::SectionId{0},
+    });
+
+    std::vector<fesa::PrescribedDof> prescribed =
+        fixed_dofs(fesa::NodeId{20});
+    std::vector<fesa::PrescribedDof> end_constraints = fixed_dofs(
+        fesa::NodeId{4}, all_constrained ? 0.1 : 0.0);
+    if (!all_constrained) {
+        end_constraints.erase(end_constraints.begin());
+    }
+    prescribed.insert(
+        prescribed.end(),
+        end_constraints.begin(),
+        end_constraints.end());
+
+    std::vector<fesa::NodalLoad> loads;
+    if (!all_constrained) {
+        loads.push_back({
+            fesa::NodeId{4},
+            {10.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+        });
+    }
+    builder.set_step({"Load", std::move(prescribed), std::move(loads)});
+
+    auto result = std::move(builder).build();
+    if (!result.domain.has_value()) {
+        throw std::runtime_error{"Test Domain failed validation."};
+    }
+    return std::move(*result.domain);
+}
+
+std::vector<double> full_field(
+    const fesa::NodalFrame& nodal,
+    const fesa::DofManager& dofs,
+    const bool reaction) {
+    std::vector<double> field(dofs.full_dof_count(), 0.0);
+    for (std::size_t node = 0; node < nodal.node_ids.size(); ++node) {
+        const auto& values = reaction
+            ? nodal.reaction[node]
+            : nodal.displacement[node];
+        for (std::size_t component = 0; component < dofs_per_node;
+             ++component) {
+            field[dofs.full_dof({
+                nodal.node_ids[node],
+                static_cast<fesa::NodeDof>(component),
+            })] = values[component];
+        }
+    }
+    return field;
+}
+
+std::vector<double> multiply(
+    const fesa::SymmetricCsr& matrix,
+    const std::vector<double>& vector) {
+    std::vector<double> product(matrix.order, 0.0);
+    for (std::size_t row = 0; row < matrix.order; ++row) {
+        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::size_t column = static_cast<std::size_t>(
+                matrix.column_indices[entry]);
+            const double value = matrix.values[entry];
+            product[row] += value * vector[column];
+            if (column != row) {
+                product[column] += value * vector[row];
+            }
+        }
+    }
+    return product;
+}
+
+TEST(LinearStaticAnalysis, SolvesHandCalculatedAxialBeamInNodeIdOrder) {
+    const fesa::Domain domain = build_axial_domain(false);
+
+    const fesa::AnalysisRunResult run =
+        fesa::LinearStaticAnalysis{}.run(domain);
+
+    ASSERT_TRUE(run.succeeded);
+    ASSERT_TRUE(run.results.has_value());
+    EXPECT_TRUE(run.diagnostics.empty());
+    EXPECT_TRUE(fesa::validate_result_database(*run.results).succeeded);
+    EXPECT_EQ(run.results->schema_version, "1.0.0");
+    ASSERT_EQ(run.results->steps.size(), 1);
+    EXPECT_EQ(run.results->steps[0].name, "Load");
+    ASSERT_EQ(run.results->steps[0].frames.size(), 1);
+    const fesa::ResultFrame& frame = run.results->steps[0].frames[0];
+    EXPECT_DOUBLE_EQ(frame.step_time, 1.0);
+    EXPECT_TRUE(frame.diagnostics.empty());
+    ASSERT_EQ(frame.nodal.node_ids.size(), 2);
+    EXPECT_EQ(
+        frame.nodal.node_ids,
+        (std::vector<fesa::NodeId>{
+            fesa::NodeId{4}, fesa::NodeId{20}}));
+    EXPECT_NEAR(frame.nodal.displacement[0][0], 0.2, 1.0e-12);
+    EXPECT_NEAR(frame.nodal.displacement[1][0], 0.0, 1.0e-12);
+    EXPECT_NEAR(frame.nodal.reaction[0][0], 0.0, 1.0e-12);
+    EXPECT_NEAR(frame.nodal.reaction[1][0], -10.0, 1.0e-12);
+}
+
+TEST(StaticEquilibrium, ReturnedFieldsSatisfyOriginalFullEquation) {
+    const fesa::Domain domain = build_axial_domain(false);
+    const fesa::DofManager dofs = fesa::DofManager::build(domain);
+    const fesa::EquationSystem original =
+        fesa::assemble_serial(domain, dofs);
+
+    const fesa::AnalysisRunResult run =
+        fesa::LinearStaticAnalysis{}.run(domain);
+
+    ASSERT_TRUE(run.succeeded);
+    ASSERT_TRUE(run.results.has_value());
+    const fesa::NodalFrame& nodal =
+        run.results->steps[0].frames[0].nodal;
+    const std::vector<double> displacement =
+        full_field(nodal, dofs, false);
+    const std::vector<double> reaction =
+        full_field(nodal, dofs, true);
+    const std::vector<double> internal =
+        multiply(original.stiffness, displacement);
+    for (std::size_t dof = 0; dof < internal.size(); ++dof) {
+        EXPECT_NEAR(
+            internal[dof] - original.force[dof] - reaction[dof],
+            0.0,
+            1.0e-12);
+        EXPECT_TRUE(std::isfinite(displacement[dof]));
+        EXPECT_TRUE(std::isfinite(reaction[dof]));
+    }
+}
+
+TEST(LinearStaticAnalysis, SolvesAllConstrainedSystemWithoutPardiso) {
+    const fesa::Domain domain = build_axial_domain(true);
+
+    const fesa::AnalysisRunResult run =
+        fesa::LinearStaticAnalysis{}.run(domain);
+
+    ASSERT_TRUE(run.succeeded);
+    ASSERT_TRUE(run.results.has_value());
+    EXPECT_TRUE(run.diagnostics.empty());
+    const fesa::NodalFrame& nodal =
+        run.results->steps[0].frames[0].nodal;
+    ASSERT_EQ(nodal.node_ids.size(), 2);
+    ASSERT_EQ(nodal.node_ids[0], fesa::NodeId{4});
+    EXPECT_NEAR(nodal.displacement[0][0], 0.1, 1.0e-12);
+    EXPECT_NEAR(nodal.reaction[0][0], 5.0, 1.0e-12);
+    EXPECT_NEAR(nodal.reaction[1][0], -5.0, 1.0e-12);
+}
+
+}  // namespace