From 6c2e1f3ab1e40bb28727b7dfc5874c9924cef57f Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Sat, 1 Aug 2026 23:21:27 +0900 Subject: [PATCH] =?UTF-8?q?feat(deterministic-parallel-assembly):=20step?= =?UTF-8?q?=201=20=E2=80=94=20tbb-element-evaluation?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CMakeLists.txt | 3 +- include/fesa/assembly/assembler.hpp | 21 ++ src/fesa/assembly/parallel_assembler.cpp | 208 ++++++++++++++++++ tests/CMakeLists.txt | 37 ++++ .../unit/assembly/parallel_assembler_test.cpp | 189 ++++++++++++++++ 5 files changed, 457 insertions(+), 1 deletion(-) create mode 100644 include/fesa/assembly/assembler.hpp create mode 100644 src/fesa/assembly/parallel_assembler.cpp create mode 100644 tests/unit/assembly/parallel_assembler_test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 4418637..388f957 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -29,6 +29,7 @@ add_library(fesa_core STATIC src/fesa/analysis/linear_static_analysis.cpp src/fesa/analysis/run_solver.cpp src/fesa/assembly/contribution.cpp + src/fesa/assembly/parallel_assembler.cpp src/fesa/assembly/serial_assembler.cpp src/fesa/constraints/essential_bc.cpp src/fesa/core/version.cpp @@ -55,7 +56,7 @@ target_include_directories(fesa_core target_compile_features(fesa_core PUBLIC cxx_std_20) target_compile_options(fesa_core PRIVATE /W4 /permissive- /EHsc) -target_link_libraries(fesa_core PRIVATE MKL::MKL HDF5::HDF5) +target_link_libraries(fesa_core PRIVATE MKL::MKL TBB::tbb HDF5::HDF5) add_executable(fesa src/fesa/cli/main.cpp diff --git a/include/fesa/assembly/assembler.hpp b/include/fesa/assembly/assembler.hpp new file mode 100644 index 0000000..fea708b --- /dev/null +++ b/include/fesa/assembly/assembler.hpp @@ -0,0 +1,21 @@ +#pragma once + +#include + +#include +#include +#include + +namespace fesa { + +struct AssemblyOptions final { + std::size_t max_threads; + std::size_t grain_size; +}; + +[[nodiscard]] EquationSystem assemble_parallel( + const Domain& domain, + const DofManager& dofs, + AssemblyOptions options); + +} // namespace fesa diff --git a/src/fesa/assembly/parallel_assembler.cpp b/src/fesa/assembly/parallel_assembler.cpp new file mode 100644 index 0000000..c7285e8 --- /dev/null +++ b/src/fesa/assembly/parallel_assembler.cpp @@ -0,0 +1,208 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +namespace fesa { +namespace { + +struct ElementEvaluation final { + std::vector contributions; + std::optional error; +}; + +auto origin_key(const EntityOrigin& origin) { + return std::tie( + origin.instance_name, + origin.local_label, + origin.part_name); +} + +std::string kernel_error_message( + const BeamElement& element, + const BeamKernelResult& result) { + std::string message = + "Beam element " + std::to_string(element.origin.local_label) + + " kernel failed"; + for (const Diagnostic& diagnostic : result.diagnostics) { + message += ": " + diagnostic.code + " - " + diagnostic.message; + } + return message; +} + +std::vector canonical_element_order(const Domain& domain) { + std::vector order(domain.beam_elements().size()); + std::iota(order.begin(), order.end(), std::size_t{0}); + std::ranges::sort( + order, + [&domain](const std::size_t left, const std::size_t right) { + return origin_key(domain.beam_elements()[left].origin) < + origin_key(domain.beam_elements()[right].origin); + }); + return order; +} + +std::vector canonical_contribution_element_ids( + const std::vector& order) { + // Domain ElementIds are not ordered by input identity. These tie-break + // IDs encode the existing serial assembler's element-origin order. + std::vector ids(order.size(), ElementId{0}); + for (std::size_t rank = 0; rank < order.size(); ++rank) { + ids[order[rank]] = ElementId{static_cast(rank)}; + } + return ids; +} + +ElementEvaluation evaluate_element( + const Domain& domain, + const DofManager& dofs, + const BeamElement& element, + const ElementId canonical_id) { + const BeamKernelResult result = compute_beam3d2({ + { + domain.node(element.nodes[0]).position, + domain.node(element.nodes[1]).position, + }, + domain.material(element.material), + domain.section(element.section), + }); + if (!result.contribution.has_value()) { + return {{}, kernel_error_message(element, result)}; + } + + ElementEvaluation evaluation; + evaluation.contributions.reserve(78); + const std::array full_dofs = + dofs.element_full_dofs(element); + for (std::size_t local_row = 0; local_row < full_dofs.size(); + ++local_row) { + for (std::size_t local_column = local_row; + local_column < full_dofs.size(); + ++local_column) { + evaluation.contributions.push_back({ + std::min( + full_dofs[local_row], + full_dofs[local_column]), + std::max( + full_dofs[local_row], + full_dofs[local_column]), + canonical_id, + static_cast( + local_row * full_dofs.size() + local_column), + result.contribution + ->global_stiffness[local_row][local_column], + }); + } + } + return evaluation; +} + +std::vector assemble_force( + const Domain& domain, + const DofManager& dofs) { + std::vector force(dofs.full_dof_count(), 0.0); + for (const NodalLoad& load : domain.step().nodal_loads) { + for (std::size_t component = 0; component < load.values.size(); + ++component) { + const auto dof = static_cast(component); + force[dofs.full_dof({load.node, dof})] += + load.values[component]; + } + } + return force; +} + +void validate_options(const AssemblyOptions options) { + if (options.max_threads == 0) { + throw std::invalid_argument{ + "Assembly max_threads must be greater than zero."}; + } + if (options.max_threads > + static_cast(std::numeric_limits::max())) { + throw std::invalid_argument{ + "Assembly max_threads exceeds the TBB task arena range."}; + } + if (options.grain_size == 0) { + throw std::invalid_argument{ + "Assembly grain_size must be greater than zero."}; + } +} + +} // namespace + +EquationSystem assemble_parallel( + const Domain& domain, + const DofManager& dofs, + const AssemblyOptions options) { + validate_options(options); + + const std::vector element_order = + canonical_element_order(domain); + const std::vector element_ids = + canonical_contribution_element_ids(element_order); + std::vector evaluations( + domain.beam_elements().size()); + + oneapi::tbb::task_arena arena{ + static_cast(options.max_threads)}; + arena.execute([&] { + oneapi::tbb::parallel_for( + oneapi::tbb::blocked_range{ + 0, + evaluations.size(), + options.grain_size, + }, + [&](const oneapi::tbb::blocked_range& range) { + for (std::size_t index = range.begin(); + index != range.end(); + ++index) { + ElementEvaluation local = evaluate_element( + domain, + dofs, + domain.beam_elements()[index], + element_ids[index]); + evaluations[index] = std::move(local); + } + }); + }); + + std::vector contributions; + contributions.reserve(domain.beam_elements().size() * 78); + for (const std::size_t element_index : element_order) { + ElementEvaluation& evaluation = evaluations[element_index]; + if (evaluation.error.has_value()) { + throw std::runtime_error{std::move(*evaluation.error)}; + } + contributions.insert( + contributions.end(), + std::make_move_iterator(evaluation.contributions.begin()), + std::make_move_iterator(evaluation.contributions.end())); + } + + const std::vector canonical = + canonicalize_contributions(contributions); + return { + merge_contributions(dofs.full_dof_count(), canonical), + assemble_force(domain, dofs), + }; +} + +} // namespace fesa diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 7fc1dda..1aad74c 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -409,6 +409,43 @@ add_test( --gtest_filter=DeterministicMerge.* ) +add_executable(fesa_parallel_assembly_tests + unit/assembly/parallel_assembler_test.cpp +) + +target_compile_features(fesa_parallel_assembly_tests PRIVATE cxx_std_20) +target_compile_options( + fesa_parallel_assembly_tests + PRIVATE + /W4 + /permissive- + /EHsc +) + +target_link_libraries(fesa_parallel_assembly_tests + PRIVATE + fesa_core + GTest::gtest_main +) + +add_test( + NAME ParallelAssembly + COMMAND "$" + --gtest_filter=ParallelAssembly.* +) + +add_test( + NAME TbbElementEvaluation + COMMAND "$" + --gtest_filter=TbbElementEvaluation.* +) + +set_property( + TEST ParallelAssembly TbbElementEvaluation + PROPERTY ENVIRONMENT_MODIFICATION + ${FESA_DEPENDENCY_RUNTIME_MODIFICATIONS} +) + add_executable(fesa_constraint_tests unit/constraints/essential_bc_test.cpp ) diff --git a/tests/unit/assembly/parallel_assembler_test.cpp b/tests/unit/assembly/parallel_assembler_test.cpp new file mode 100644 index 0000000..e10ff15 --- /dev/null +++ b/tests/unit/assembly/parallel_assembler_test.cpp @@ -0,0 +1,189 @@ +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +namespace { + +fesa::Domain finish_domain( + fesa::DomainBuilder builder, + fesa::StepDefinition step) { + builder.set_step(std::move(step)); + auto result = std::move(builder).build(); + if (!result.domain.has_value()) { + throw std::runtime_error{"Test Domain failed validation."}; + } + return std::move(*result.domain); +} + +fesa::Domain build_beam_domain(const bool branched) { + fesa::DomainBuilder builder; + const std::array positions{{ + {0.0, 0.0, 0.0}, + {1.0, 0.0, 0.0}, + {2.0, 0.0, 0.0}, + {3.0, 0.0, 0.0}, + {1.0, 0.0, 1.0}, + {1.0, 1.0, 1.0}, + }}; + const std::size_t node_count = branched ? positions.size() : 4; + for (std::size_t index = 0; index < node_count; ++index) { + builder.add_node({ + fesa::NodeId{static_cast(index)}, + fesa::EntityOrigin{ + "BeamPart", + "Beam-1", + static_cast(index + 1), + }, + positions[index], + }); + } + builder.add_material({ + fesa::MaterialId{0}, + "Steel", + 210.0e9, + 0.3, + }); + builder.add_section({ + fesa::SectionId{0}, + "General", + 0.02, + 3.0e-5, + 4.0e-5, + 2.0e-5, + 0.015, + 0.016, + fesa::ShearPropertySource::input, + fesa::Vec3{0.0, 1.0, 0.0}, + {}, + }); + + const std::array elements{{ + { + fesa::ElementId{40}, + fesa::EntityOrigin{"BeamPart", "Beam-1", 300}, + {fesa::NodeId{2}, fesa::NodeId{3}}, + fesa::MaterialId{0}, + fesa::SectionId{0}, + }, + { + fesa::ElementId{10}, + fesa::EntityOrigin{"BeamPart", "Beam-1", 100}, + {fesa::NodeId{0}, fesa::NodeId{1}}, + fesa::MaterialId{0}, + fesa::SectionId{0}, + }, + { + fesa::ElementId{30}, + fesa::EntityOrigin{"BeamPart", "Beam-1", 200}, + {fesa::NodeId{1}, fesa::NodeId{2}}, + fesa::MaterialId{0}, + fesa::SectionId{0}, + }, + { + fesa::ElementId{0}, + fesa::EntityOrigin{"BeamPart", "Beam-1", 500}, + {fesa::NodeId{1}, fesa::NodeId{5}}, + fesa::MaterialId{0}, + fesa::SectionId{0}, + }, + { + fesa::ElementId{20}, + fesa::EntityOrigin{"BeamPart", "Beam-1", 400}, + {fesa::NodeId{1}, fesa::NodeId{4}}, + fesa::MaterialId{0}, + fesa::SectionId{0}, + }, + }}; + const std::size_t element_count = branched ? elements.size() : 3; + for (std::size_t index = 0; index < element_count; ++index) { + builder.add_beam_element(elements[index]); + } + + return finish_domain( + std::move(builder), + { + "Load", + {}, + { + {fesa::NodeId{3}, {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}}, + {fesa::NodeId{1}, {1.0e16, 1.0, 0.0, 0.0, 0.0, 0.0}}, + {fesa::NodeId{1}, {-1.0e16, 2.0, 0.0, 0.0, 0.0, 0.0}}, + }, + }); +} + +std::vector bits(const std::vector& values) { + std::vector result; + result.reserve(values.size()); + for (const double value : values) { + result.push_back(std::bit_cast(value)); + } + return result; +} + +void expect_bitwise_equal( + const fesa::EquationSystem& expected, + const fesa::EquationSystem& actual) { + EXPECT_EQ(actual.stiffness.order, expected.stiffness.order); + EXPECT_EQ( + actual.stiffness.row_offsets, + expected.stiffness.row_offsets); + EXPECT_EQ( + actual.stiffness.column_indices, + expected.stiffness.column_indices); + EXPECT_EQ(bits(actual.stiffness.values), bits(expected.stiffness.values)); + EXPECT_EQ(bits(actual.force), bits(expected.force)); +} + +TEST(ParallelAssembly, MatchesSerialForFixedBeamChain) { + const fesa::Domain domain = build_beam_domain(false); + const fesa::DofManager dofs = fesa::DofManager::build(domain); + const fesa::EquationSystem serial = fesa::assemble_serial(domain, dofs); + + expect_bitwise_equal( + serial, + fesa::assemble_parallel(domain, dofs, {1, 1})); + expect_bitwise_equal( + serial, + fesa::assemble_parallel(domain, dofs, {4, 2})); +} + +TEST(ParallelAssembly, MatchesSerialForFixedBranchedDomain) { + const fesa::Domain domain = build_beam_domain(true); + const fesa::DofManager dofs = fesa::DofManager::build(domain); + const fesa::EquationSystem serial = fesa::assemble_serial(domain, dofs); + + expect_bitwise_equal( + serial, + fesa::assemble_parallel(domain, dofs, {1, 2})); + expect_bitwise_equal( + serial, + fesa::assemble_parallel(domain, dofs, {3, 1})); +} + +TEST(TbbElementEvaluation, RejectsZeroExecutionLimits) { + const fesa::Domain domain = build_beam_domain(false); + const fesa::DofManager dofs = fesa::DofManager::build(domain); + + EXPECT_THROW( + static_cast( + fesa::assemble_parallel(domain, dofs, {0, 1})), + std::invalid_argument); + EXPECT_THROW( + static_cast( + fesa::assemble_parallel(domain, dofs, {1, 0})), + std::invalid_argument); +} + +} // namespace