feat(fem-and-beam-kernel): step 1 — dof-manager

CMakeLists.txt

@@ -27,6 +27,7 @@ include(cmake/FesaDependencies.cmake)
 
 add_library(fesa_core STATIC
     src/fesa/core/version.cpp
+    src/fesa/fem/dof_manager.cpp
     src/fesa/fem/gauss_rule.cpp
     src/fesa/fem/line2_shape.cpp
     src/fesa/io/abaqus/parser.cpp

include/fesa/fem/dof_manager.hpp

@@ -0,0 +1,52 @@
+#pragma once
+
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <map>
+#include <optional>
+#include <span>
+#include <vector>
+
+#include <fesa/model/domain.hpp>
+
+namespace fesa {
+
+enum class NodeDof : std::uint8_t {
+    ux,
+    uy,
+    uz,
+    rx,
+    ry,
+    rz,
+};
+
+struct DofAddress final {
+    NodeId node;
+    NodeDof dof;
+};
+
+class DofManager final {
+public:
+    [[nodiscard]] static DofManager build(const Domain& domain);
+
+    [[nodiscard]] std::size_t full_dof_count() const noexcept;
+    [[nodiscard]] std::size_t free_equation_count() const noexcept;
+    [[nodiscard]] std::optional<std::size_t> equation(
+        DofAddress address) const;
+    [[nodiscard]] std::array<std::size_t, 12> element_full_dofs(
+        const BeamElement& element) const;
+    [[nodiscard]] std::vector<double> reconstruct_full(
+        std::span<const double> reduced) const;
+
+private:
+    [[nodiscard]] std::size_t full_dof(DofAddress address) const;
+    [[nodiscard]] std::size_t node_full_dof_base(NodeId node) const;
+
+    std::map<std::int64_t, std::size_t> node_full_dof_bases_;
+    std::vector<std::optional<std::size_t>> equations_;
+    std::vector<double> prescribed_values_;
+    std::size_t free_equation_count_ = 0;
+};
+
+}  // namespace fesa

src/fesa/fem/dof_manager.cpp

@@ -0,0 +1,111 @@
+#include <fesa/fem/dof_manager.hpp>
+
+#include <algorithm>
+#include <stdexcept>
+#include <vector>
+
+namespace fesa {
+namespace {
+
+constexpr std::size_t dofs_per_node = 6;
+
+}  // namespace
+
+DofManager DofManager::build(const Domain& domain) {
+    DofManager manager;
+
+    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);
+
+    for (std::size_t index = 0; index < node_ids.size(); ++index) {
+        manager.node_full_dof_bases_.emplace(
+            node_ids[index].value(), index * dofs_per_node);
+    }
+
+    const std::size_t full_dof_count = node_ids.size() * dofs_per_node;
+    manager.equations_.resize(full_dof_count);
+    manager.prescribed_values_.resize(full_dof_count);
+    std::vector<bool> constrained(full_dof_count, false);
+
+    for (const PrescribedDof& prescribed : domain.step().prescribed_dofs) {
+        const NodeDof dof =
+            static_cast<NodeDof>(prescribed.dof - std::uint8_t{1});
+        const std::size_t full =
+            manager.full_dof({prescribed.node, dof});
+        constrained[full] = true;
+        manager.prescribed_values_[full] = prescribed.value;
+    }
+
+    for (std::size_t full = 0; full < full_dof_count; ++full) {
+        if (!constrained[full]) {
+            manager.equations_[full] = manager.free_equation_count_;
+            ++manager.free_equation_count_;
+        }
+    }
+
+    return manager;
+}
+
+std::size_t DofManager::full_dof_count() const noexcept {
+    return equations_.size();
+}
+
+std::size_t DofManager::free_equation_count() const noexcept {
+    return free_equation_count_;
+}
+
+std::optional<std::size_t> DofManager::equation(
+    const DofAddress address) const {
+    return equations_[full_dof(address)];
+}
+
+std::array<std::size_t, 12> DofManager::element_full_dofs(
+    const BeamElement& element) const {
+    std::array<std::size_t, 12> full_dofs{};
+    for (std::size_t node = 0; node < element.nodes.size(); ++node) {
+        const std::size_t base = node_full_dof_base(element.nodes[node]);
+        for (std::size_t component = 0; component < dofs_per_node;
+             ++component) {
+            full_dofs[node * dofs_per_node + component] = base + component;
+        }
+    }
+    return full_dofs;
+}
+
+std::vector<double> DofManager::reconstruct_full(
+    const std::span<const double> reduced) const {
+    if (reduced.size() != free_equation_count_) {
+        throw std::invalid_argument{
+            "Reduced vector size must equal the free equation count."};
+    }
+
+    std::vector<double> full = prescribed_values_;
+    for (std::size_t index = 0; index < equations_.size(); ++index) {
+        if (equations_[index].has_value()) {
+            full[index] = reduced[*equations_[index]];
+        }
+    }
+    return full;
+}
+
+std::size_t DofManager::full_dof(const DofAddress address) const {
+    const auto component = static_cast<std::uint8_t>(address.dof);
+    if (component >= dofs_per_node) {
+        throw std::invalid_argument{"Node DOF must be one of ux through rz."};
+    }
+    return node_full_dof_base(address.node) + component;
+}
+
+std::size_t DofManager::node_full_dof_base(const NodeId node) const {
+    const auto found = node_full_dof_bases_.find(node.value());
+    if (found == node_full_dof_bases_.end()) {
+        throw std::out_of_range{"Node ID is not present in the DofManager."};
+    }
+    return found->second;
+}
+
+}  // namespace fesa

tests/CMakeLists.txt

@@ -172,6 +172,7 @@ add_test(
 )
 
 add_executable(fesa_fem_primitives_tests
+    unit/fem/dof_manager_test.cpp
     unit/fem/gauss_rule_test.cpp
     unit/fem/line2_shape_test.cpp
 )
@@ -208,3 +209,15 @@ add_test(
     COMMAND "$<TARGET_FILE:fesa_fem_primitives_tests>"
         --gtest_filter=Jacobian.*
 )
+
+add_test(
+    NAME DofManager
+    COMMAND "$<TARGET_FILE:fesa_fem_primitives_tests>"
+        --gtest_filter=DofManager.*
+)
+
+add_test(
+    NAME EquationNumbering
+    COMMAND "$<TARGET_FILE:fesa_fem_primitives_tests>"
+        --gtest_filter=EquationNumbering.*
+)

tests/unit/fem/dof_manager_test.cpp

@@ -0,0 +1,199 @@
+#include <fesa/fem/dof_manager.hpp>
+
+#include <array>
+#include <cstdint>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+#include <fesa/model/domain_builder.hpp>
+
+#include <gtest/gtest.h>
+
+namespace {
+
+fesa::Domain build_domain(
+    std::vector<fesa::Node> nodes,
+    fesa::StepDefinition step,
+    const std::array<fesa::NodeId, 2> element_nodes = {
+        fesa::NodeId{4},
+        fesa::NodeId{10},
+    }) {
+    fesa::DomainBuilder builder;
+    for (fesa::Node& node : nodes) {
+        builder.add_node(std::move(node));
+    }
+    builder.add_material({
+        fesa::MaterialId{0},
+        "Steel",
+        210.0e9,
+        0.3,
+    });
+    builder.add_section({
+        fesa::SectionId{0},
+        "General",
+        0.04,
+        1.2e-4,
+        1.4e-4,
+        2.0e-4,
+        0.03,
+        0.031,
+        fesa::ShearPropertySource::input,
+        fesa::Vec3{0.0, 1.0, 0.0},
+        {},
+    });
+    builder.add_beam_element({
+        fesa::ElementId{0},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 1},
+        element_nodes,
+        fesa::MaterialId{0},
+        fesa::SectionId{0},
+    });
+    builder.set_step(std::move(step));
+
+    auto result = std::move(builder).build();
+    EXPECT_TRUE(result.diagnostics.empty());
+    EXPECT_TRUE(result.domain.has_value());
+    return std::move(*result.domain);
+}
+
+std::vector<fesa::Node> nodes_in_external_label_order() {
+    return {
+        {
+            fesa::NodeId{20},
+            fesa::EntityOrigin{"BeamPart", "Beam-1", 2},
+            fesa::Vec3{2.0, 0.0, 0.0},
+        },
+        {
+            fesa::NodeId{4},
+            fesa::EntityOrigin{"BeamPart", "Beam-1", 100},
+            fesa::Vec3{0.0, 0.0, 0.0},
+        },
+        {
+            fesa::NodeId{10},
+            fesa::EntityOrigin{"BeamPart", "Beam-1", 50},
+            fesa::Vec3{1.0, 0.0, 0.0},
+        },
+    };
+}
+
+TEST(EquationNumbering, UsesInternalNodeIdOrderNotExternalLabelOrStorageOrder) {
+    const fesa::Domain domain =
+        build_domain(nodes_in_external_label_order(), {"Load", {}, {}});
+
+    const fesa::DofManager dofs = fesa::DofManager::build(domain);
+
+    EXPECT_EQ(dofs.full_dof_count(), 18);
+    EXPECT_EQ(dofs.free_equation_count(), 18);
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{4}, fesa::NodeDof::ux}),
+        std::optional<std::size_t>{0});
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{4}, fesa::NodeDof::rz}),
+        std::optional<std::size_t>{5});
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{10}, fesa::NodeDof::ux}),
+        std::optional<std::size_t>{6});
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{20}, fesa::NodeDof::rz}),
+        std::optional<std::size_t>{17});
+}
+
+TEST(DofManager, NumbersOnlyFreeDofsAndReconstructsPrescribedValues) {
+    const fesa::Domain domain = build_domain(
+        nodes_in_external_label_order(),
+        {
+            "Load",
+            {
+                {fesa::NodeId{4}, 1, 1.25},
+                {fesa::NodeId{10}, 6, -0.5},
+            },
+            {},
+        });
+
+    const fesa::DofManager dofs = fesa::DofManager::build(domain);
+
+    EXPECT_EQ(dofs.full_dof_count(), 18);
+    EXPECT_EQ(dofs.free_equation_count(), 16);
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{4}, fesa::NodeDof::ux}),
+        std::nullopt);
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{4}, fesa::NodeDof::uy}),
+        std::optional<std::size_t>{0});
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{10}, fesa::NodeDof::rx}),
+        std::optional<std::size_t>{8});
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{10}, fesa::NodeDof::rz}),
+        std::nullopt);
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{20}, fesa::NodeDof::ux}),
+        std::optional<std::size_t>{10});
+
+    std::vector<double> reduced(dofs.free_equation_count());
+    for (std::size_t equation = 0; equation < reduced.size(); ++equation) {
+        reduced[equation] = 100.0 + static_cast<double>(equation);
+    }
+
+    EXPECT_EQ(
+        dofs.reconstruct_full(reduced),
+        (std::vector<double>{
+            1.25,
+            100.0,
+            101.0,
+            102.0,
+            103.0,
+            104.0,
+            105.0,
+            106.0,
+            107.0,
+            108.0,
+            109.0,
+            -0.5,
+            110.0,
+            111.0,
+            112.0,
+            113.0,
+            114.0,
+            115.0,
+        }));
+}
+
+TEST(DofManager, MapsElementNodesToTwelveFullDofsInComponentOrder) {
+    const fesa::Domain domain = build_domain(
+        nodes_in_external_label_order(),
+        {"Load", {}, {}},
+        {fesa::NodeId{10}, fesa::NodeId{4}});
+    const fesa::DofManager dofs = fesa::DofManager::build(domain);
+
+    EXPECT_EQ(
+        dofs.element_full_dofs(domain.beam_elements().front()),
+        (std::array<std::size_t, 12>{
+            6, 7, 8, 9, 10, 11,
+            0, 1, 2, 3, 4, 5,
+        }));
+}
+
+TEST(DofManager, RejectsInvalidAddressesAndReducedVectorSize) {
+    const fesa::Domain domain =
+        build_domain(nodes_in_external_label_order(), {"Load", {}, {}});
+    const fesa::DofManager dofs = fesa::DofManager::build(domain);
+
+    EXPECT_THROW(
+        static_cast<void>(dofs.equation({
+            fesa::NodeId{4},
+            static_cast<fesa::NodeDof>(6),
+        })),
+        std::invalid_argument);
+    EXPECT_THROW(
+        static_cast<void>(
+            dofs.equation({fesa::NodeId{99}, fesa::NodeDof::ux})),
+        std::out_of_range);
+    EXPECT_THROW(
+        static_cast<void>(
+            dofs.reconstruct_full(std::vector<double>(17, 0.0))),
+        std::invalid_argument);
+}
+
+}  // namespace