diff --git a/CMakeLists.txt b/CMakeLists.txt index 500731a..98f165b 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/include/fesa/fem/dof_manager.hpp b/include/fesa/fem/dof_manager.hpp new file mode 100644 index 0000000..77736df --- /dev/null +++ b/include/fesa/fem/dof_manager.hpp @@ -0,0 +1,52 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include + +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 equation( + DofAddress address) const; + [[nodiscard]] std::array element_full_dofs( + const BeamElement& element) const; + [[nodiscard]] std::vector reconstruct_full( + std::span 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 node_full_dof_bases_; + std::vector> equations_; + std::vector prescribed_values_; + std::size_t free_equation_count_ = 0; +}; + +} // namespace fesa diff --git a/src/fesa/fem/dof_manager.cpp b/src/fesa/fem/dof_manager.cpp new file mode 100644 index 0000000..072fd1a --- /dev/null +++ b/src/fesa/fem/dof_manager.cpp @@ -0,0 +1,111 @@ +#include + +#include +#include +#include + +namespace fesa { +namespace { + +constexpr std::size_t dofs_per_node = 6; + +} // namespace + +DofManager DofManager::build(const Domain& domain) { + DofManager manager; + + std::vector 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 constrained(full_dof_count, false); + + for (const PrescribedDof& prescribed : domain.step().prescribed_dofs) { + const NodeDof dof = + static_cast(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 DofManager::equation( + const DofAddress address) const { + return equations_[full_dof(address)]; +} + +std::array DofManager::element_full_dofs( + const BeamElement& element) const { + std::array 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 DofManager::reconstruct_full( + const std::span reduced) const { + if (reduced.size() != free_equation_count_) { + throw std::invalid_argument{ + "Reduced vector size must equal the free equation count."}; + } + + std::vector 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(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 diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 1edf7ce..be416c7 100644 --- a/tests/CMakeLists.txt +++ b/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 "$" --gtest_filter=Jacobian.* ) + +add_test( + NAME DofManager + COMMAND "$" + --gtest_filter=DofManager.* +) + +add_test( + NAME EquationNumbering + COMMAND "$" + --gtest_filter=EquationNumbering.* +) diff --git a/tests/unit/fem/dof_manager_test.cpp b/tests/unit/fem/dof_manager_test.cpp new file mode 100644 index 0000000..6dd9fc0 --- /dev/null +++ b/tests/unit/fem/dof_manager_test.cpp @@ -0,0 +1,199 @@ +#include + +#include +#include +#include +#include +#include + +#include + +#include + +namespace { + +fesa::Domain build_domain( + std::vector nodes, + fesa::StepDefinition step, + const std::array 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 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{0}); + EXPECT_EQ( + dofs.equation({fesa::NodeId{4}, fesa::NodeDof::rz}), + std::optional{5}); + EXPECT_EQ( + dofs.equation({fesa::NodeId{10}, fesa::NodeDof::ux}), + std::optional{6}); + EXPECT_EQ( + dofs.equation({fesa::NodeId{20}, fesa::NodeDof::rz}), + std::optional{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{0}); + EXPECT_EQ( + dofs.equation({fesa::NodeId{10}, fesa::NodeDof::rx}), + std::optional{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{10}); + + std::vector reduced(dofs.free_equation_count()); + for (std::size_t equation = 0; equation < reduced.size(); ++equation) { + reduced[equation] = 100.0 + static_cast(equation); + } + + EXPECT_EQ( + dofs.reconstruct_full(reduced), + (std::vector{ + 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{ + 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(dofs.equation({ + fesa::NodeId{4}, + static_cast(6), + })), + std::invalid_argument); + EXPECT_THROW( + static_cast( + dofs.equation({fesa::NodeId{99}, fesa::NodeDof::ux})), + std::out_of_range); + EXPECT_THROW( + static_cast( + dofs.reconstruct_full(std::vector(17, 0.0))), + std::invalid_argument); +} + +} // namespace