feat(cpp-object-oriented-modular-refactoring): step 16 - generic-sparse-assembler

This commit is contained in:
KOKO\Mimi
2026-08-16 10:24:28 +09:00
parent a21b991ef9
commit 6b10f8a7d2
3 changed files with 219 additions and 259 deletions
+14 -1
View File
@@ -2,6 +2,7 @@
#define FESA_ASSEMBLY_SPARSE_ASSEMBLER_H_ #define FESA_ASSEMBLY_SPARSE_ASSEMBLER_H_
#include "fesa/core/status.h" #include "fesa/core/status.h"
#include "fesa/elements/element.h"
#include "fesa/math/sparse_matrix.h" #include "fesa/math/sparse_matrix.h"
namespace fesa { namespace fesa {
@@ -13,10 +14,22 @@ class ParallelFor;
/// @brief Owns deterministic element-contribution reduction into global CSR. /// @brief Owns deterministic element-contribution reduction into global CSR.
class SparseAssembler { class SparseAssembler {
public: public:
/// @brief Assembles runtime element stiffness into validated full-DOF CSR.
/// @param elements Non-owning elements in stable active source order.
/// @param dofs Owner of the matching scatter and structural pattern.
/// @param parallel_for Backend for index-owned element-local computation.
/// @return A validated matrix or a structured model failure.
/// @note Workers write only their element-owned COO buffers; flattening and
/// duplicate reduction retain fixed element and local-entry order.
static Result<SparseMatrix> Assemble(const ElementView& elements,
const DofManager& dofs,
const ParallelFor& parallel_for);
/// @brief Assembles stiffness in stable source-element and local-entry order. /// @brief Assembles stiffness in stable source-element and local-entry order.
/// @return A validated full-DOF CSR matrix or structured model failure. /// @return A validated full-DOF CSR matrix or structured model failure.
/// @note Parallel workers produce index-owned local buffers; serial reduction /// @note Parallel workers produce index-owned local buffers; serial reduction
/// remains the sole global CSR writer. /// remains the sole global CSR writer. This compatibility facade creates
/// runtime candidates until LinearStaticAnalysis owns them directly.
static Result<SparseMatrix> AssembleStiffness( static Result<SparseMatrix> AssembleStiffness(
const AnalysisModel& model, const DofManager& dofs, const AnalysisModel& model, const DofManager& dofs,
const ParallelFor& parallel_for); const ParallelFor& parallel_for);
+107 -258
View File
@@ -1,286 +1,97 @@
#include "fesa/assembly/sparse_assembler.h" #include "fesa/assembly/sparse_assembler.h"
#include <array> #include <cstddef>
#include <functional>
#include <limits> #include <limits>
#include <memory>
#include <optional> #include <optional>
#include <stdexcept>
#include <string> #include <string>
#include <utility> #include <utility>
#include <vector> #include <vector>
#include "fesa/analysis/analysis_model.h" #include "fesa/analysis/analysis_model.h"
#include "fesa/assembly/parallel_for.h" #include "fesa/assembly/parallel_for.h"
#include "fesa/elements/euler_beam_3d.h" #include "fesa/elements/element_factory.h"
#include "fesa/elements/mitc4_shell.h"
#include "fesa/fem/dof_manager.h" #include "fesa/fem/dof_manager.h"
#include "fesa/model/domain.h"
namespace fesa { namespace fesa {
namespace { namespace {
constexpr std::size_t kDofsPerNode = 6U; Status AssemblyFailure(const std::string& code, const SourceLocation& location,
constexpr std::size_t kBeamElementDofCount = 12U; const std::string& identity,
constexpr std::size_t kBeamContributionCount = const std::string& message) {
kBeamElementDofCount * kBeamElementDofCount; return Status::Failure(
constexpr std::size_t kShellElementDofCount = 24U;
constexpr std::size_t kShellContributionCount =
kShellElementDofCount * kShellElementDofCount;
using BeamElementBuffer = std::array<CooContribution, kBeamContributionCount>;
using ShellElementBuffer = std::array<CooContribution, kShellContributionCount>;
Result<SparseMatrix> AssemblyFailure(const std::string& code,
const SourceLocation& location,
const std::string& identity,
const std::string& message) {
return Result<SparseMatrix>::Failure(Status::Failure(
FailureCategory::kModel, FailureCategory::kModel,
{{Severity::kError, code, location, "*ELEMENT", identity, message}})); {{Severity::kError, code, location, "*ELEMENT", identity, message}});
}
bool SameSourceIdentity(const SourceEntityId& left,
const SourceEntityId& right) {
return left.instance_name == right.instance_name &&
left.source_label == right.source_label &&
left.source_label_text == right.source_label_text;
}
bool SameLayout(const ElementDofLayout& left, const ElementDofLayout& right) {
return SameSourceIdentity(left.source_id, right.source_id) &&
left.node_indices == right.node_indices &&
left.components_per_node == right.components_per_node;
} }
} // namespace } // namespace
Result<SparseMatrix> SparseAssembler::AssembleStiffness( Result<SparseMatrix> SparseAssembler::Assemble(
const AnalysisModel& model, const DofManager& dofs, const ElementView& elements, const DofManager& dofs,
const ParallelFor& parallel_for) { const ParallelFor& parallel_for) {
const Domain& domain = model.GetDomain(); std::vector<std::vector<CooContribution>> local_buffers(elements.size());
if (domain.Nodes().size() > std::vector<std::optional<Status>> local_failures(elements.size());
(std::numeric_limits<std::size_t>::max)() / kDofsPerNode || parallel_for.Execute(elements.size(), [&](const std::size_t element_order) {
dofs.FullDofCount() != domain.Nodes().size() * kDofsPerNode) { const Element& element = elements[element_order].get();
return AssemblyFailure( auto stiffness = element.ComputeStiffness();
"invalid-assembly-dimensions", {domain.SourcePath(), 0U}, if (!stiffness.HasValue()) {
std::to_string(dofs.FullDofCount()), local_failures[element_order] = stiffness.GetStatus();
"DofManager dimensions do not match the active model nodes."); return;
}
if (!model.ActiveBeamElements().empty() && !domain.ShellElements().Empty()) {
return AssemblyFailure(
"unsupported-mixed-element-model", {domain.SourcePath(), 0U},
"B33:FESA-MITC4",
"Sparse assembly does not support mixed beam and shell models.");
}
if (!domain.ShellElements().Empty()) {
if (domain.ShellElements().Size() >
(std::numeric_limits<std::size_t>::max)() / kShellContributionCount) {
return AssemblyFailure(
"invalid-assembly-dimensions", {domain.SourcePath(), 0U},
std::to_string(domain.ShellElements().Size()),
"Shell contribution storage exceeds the addressable range.");
} }
std::vector<std::optional<std::array<double, 3>>> directors_by_node( const auto& contribution = stiffness.Value();
domain.Nodes().size()); const auto& layout = element.DofLayout();
for (const auto& frame : domain.ShellNodeInitialFrames()) { if (!SameLayout(contribution.layout, layout)) {
if (frame.node_index >= directors_by_node.size() || local_failures[element_order] = AssemblyFailure(
directors_by_node[frame.node_index]) { "invalid-assembly-layout", {}, layout.source_id.source_label_text,
return AssemblyFailure( "Element stiffness layout must match its declared runtime layout.");
"invalid-assembly-element", {domain.SourcePath(), 0U}, return;
std::to_string(frame.node_index), }
"Shell initial frames must map uniquely to model nodes."); auto scatter = dofs.ElementScatter(layout);
} if (!scatter.HasValue()) {
directors_by_node[frame.node_index] = frame.director; local_failures[element_order] = scatter.GetStatus();
return;
}
const std::size_t local_dof_count = scatter.Value().size();
if (contribution.values.Rows() != local_dof_count ||
contribution.values.Columns() != local_dof_count ||
local_dof_count >
(std::numeric_limits<std::size_t>::max)() / local_dof_count) {
local_failures[element_order] = AssemblyFailure(
"invalid-assembly-dimensions", {}, layout.source_id.source_label_text,
"Element stiffness dimensions must match its declared DOF layout.");
return;
} }
struct ShellInput { auto& buffer = local_buffers[element_order];
std::array<const Node*, 4> nodes; buffer.reserve(local_dof_count * local_dof_count);
std::array<std::array<double, 3>, 4> directors; for (std::size_t local_row = 0U; local_row < local_dof_count; ++local_row) {
const ShellSection* section; for (std::size_t local_column = 0U; local_column < local_dof_count;
const LinearElasticMaterial* material; ++local_column) {
std::array<std::size_t, kShellElementDofCount> scatter; const std::size_t local_order =
}; local_row * local_dof_count + local_column;
std::vector<ShellInput> inputs; buffer.push_back({scatter.Value()[local_row],
inputs.reserve(domain.ShellElements().Size()); scatter.Value()[local_column],
for (std::size_t element_order = 0U; contribution.values(local_row, local_column),
element_order < domain.ShellElements().Size(); ++element_order) { element_order, local_order});
const auto& element = domain.ShellElements()[element_order];
if (element.material_index >= domain.LinearElasticMaterials().Size() ||
element.section_index >= domain.ShellSections().Size()) {
return AssemblyFailure(
"invalid-assembly-element", element.location,
element.source_id.source_label_text,
"Shell element references an entity outside the Domain.");
}
ShellInput input{};
input.section = &domain.ShellSections()[element.section_index];
input.material = &domain.LinearElasticMaterials()[element.material_index];
try {
input.scatter =
dofs.ShellElementScatter(static_cast<EntityIndex>(element_order));
} catch (const std::out_of_range&) {
return AssemblyFailure(
"invalid-assembly-scatter", element.location,
element.source_id.source_label_text,
"DofManager does not contain the active shell scatter.");
}
for (std::size_t node_position = 0U;
node_position < element.node_indices.size(); ++node_position) {
const EntityIndex node_index = element.node_indices[node_position];
if (node_index >= domain.Nodes().size() ||
!directors_by_node[node_index]) {
return AssemblyFailure(
"invalid-assembly-element", element.location,
element.source_id.source_label_text,
"Shell element requires a valid node and initial director.");
}
input.nodes[node_position] = &domain.Nodes()[node_index];
input.directors[node_position] = *directors_by_node[node_index];
for (std::size_t component = 0U; component < kDofsPerNode;
++component) {
const std::size_t local = node_position * kDofsPerNode + component;
const std::size_t expected =
static_cast<std::size_t>(node_index) * kDofsPerNode + component;
if (input.scatter[local] != expected ||
input.scatter[local] >= dofs.FullDofCount()) {
return AssemblyFailure(
"invalid-assembly-scatter", element.location,
element.source_id.source_label_text,
"Shell scatter does not match the active model topology.");
}
}
}
inputs.push_back(input);
}
std::vector<ShellElementBuffer> local_buffers(inputs.size());
std::vector<std::optional<Status>> local_failures(inputs.size());
parallel_for.Execute(inputs.size(), [&](const std::size_t element_order) {
const auto& input = inputs[element_order];
const auto shell = Mitc4Shell::Create(input.nodes, input.directors,
*input.section, *input.material);
if (!shell.HasValue()) {
local_failures[element_order] = shell.GetStatus();
return;
}
const auto stiffness = shell.Value().Stiffness();
if (!stiffness.HasValue()) {
local_failures[element_order] = stiffness.GetStatus();
return;
}
auto& buffer = local_buffers[element_order];
for (std::size_t local_row = 0U; local_row < kShellElementDofCount;
++local_row) {
for (std::size_t local_column = 0U;
local_column < kShellElementDofCount; ++local_column) {
const std::size_t local_order =
local_row * kShellElementDofCount + local_column;
buffer[local_order] = {
input.scatter[local_row], input.scatter[local_column],
stiffness.Value().stabilized_global24(local_row, local_column),
element_order, local_order};
}
}
});
for (std::size_t element_order = 0U; element_order < local_failures.size();
++element_order) {
if (local_failures[element_order]) {
return Result<SparseMatrix>::Failure(*local_failures[element_order]);
} }
} }
});
std::vector<CooContribution> contributions;
contributions.reserve(local_buffers.size() * kShellContributionCount);
// Flatten in source-element order after workers complete. The canonical
// COO reduction remains the sole writer of global CSR values.
for (const auto& buffer : local_buffers) {
contributions.insert(contributions.end(), buffer.begin(), buffer.end());
}
return SparseMatrix::FromCoo(dofs.FullDofCount(), dofs.FullDofCount(),
std::move(contributions),
dofs.GetSparsePattern());
}
if (model.ActiveBeamElements().size() >
(std::numeric_limits<std::size_t>::max)() / kBeamContributionCount) {
return AssemblyFailure(
"invalid-assembly-dimensions", {domain.SourcePath(), 0U},
std::to_string(model.ActiveBeamElements().size()),
"Element contribution storage exceeds the addressable range.");
}
std::vector<std::array<std::size_t, kBeamElementDofCount>> scatters;
scatters.reserve(model.ActiveBeamElements().size());
for (const EntityIndex element_index : model.ActiveBeamElements()) {
if (element_index >= domain.BeamElements().Size()) {
return AssemblyFailure("invalid-assembly-element",
{domain.SourcePath(), 0U},
std::to_string(element_index),
"Active element index is outside the Domain.");
}
const auto& element = domain.BeamElements()[element_index];
if (element.node_indices[0U] >= domain.Nodes().size() ||
element.node_indices[1U] >= domain.Nodes().size() ||
element.material_index >= domain.LinearElasticMaterials().Size() ||
element.section_index >= domain.Sections().Size()) {
return AssemblyFailure(
"invalid-assembly-element", element.location,
element.source_id.source_label_text,
"Element references an entity outside the Domain.");
}
std::array<std::size_t, kBeamElementDofCount> scatter{};
try {
scatter = dofs.ElementScatter(element_index);
} catch (const std::out_of_range&) {
return AssemblyFailure(
"invalid-assembly-scatter", element.location,
element.source_id.source_label_text,
"DofManager does not contain the active element scatter.");
}
for (std::size_t endpoint = 0U; endpoint < 2U; ++endpoint) {
for (std::size_t component = 0U; component < kDofsPerNode; ++component) {
const std::size_t local = endpoint * kDofsPerNode + component;
const std::size_t expected =
static_cast<std::size_t>(element.node_indices[endpoint]) *
kDofsPerNode +
component;
if (scatter[local] != expected ||
scatter[local] >= dofs.FullDofCount()) {
return AssemblyFailure(
"invalid-assembly-scatter", element.location,
element.source_id.source_label_text,
"Element scatter does not match the active model topology.");
}
}
}
scatters.push_back(scatter);
}
std::vector<BeamElementBuffer> local_buffers(
model.ActiveBeamElements().size());
std::vector<std::optional<Status>> local_failures(
model.ActiveBeamElements().size());
parallel_for.Execute(
model.ActiveBeamElements().size(), [&](const std::size_t element_order) {
const EntityIndex element_index =
model.ActiveBeamElements()[element_order];
const auto& definition = domain.BeamElements()[element_index];
const auto beam = EulerBeam3D::Create(
domain.Nodes()[definition.node_indices[0U]],
domain.Nodes()[definition.node_indices[1U]],
domain.Sections()[definition.section_index],
domain.LinearElasticMaterials()[definition.material_index]);
if (!beam.HasValue()) {
local_failures[element_order] = beam.GetStatus();
return;
}
const Matrix stiffness = beam.Value().GlobalStiffness();
auto& buffer = local_buffers[element_order];
const auto& scatter = scatters[element_order];
for (std::size_t local_row = 0U; local_row < kBeamElementDofCount;
++local_row) {
for (std::size_t local_column = 0U;
local_column < kBeamElementDofCount; ++local_column) {
const std::size_t local_order =
local_row * kBeamElementDofCount + local_column;
buffer[local_order] = {scatter[local_row], scatter[local_column],
stiffness(local_row, local_column),
element_order, local_order};
}
}
});
for (std::size_t element_order = 0U; element_order < local_failures.size(); for (std::size_t element_order = 0U; element_order < local_failures.size();
++element_order) { ++element_order) {
@@ -289,9 +100,20 @@ Result<SparseMatrix> SparseAssembler::AssembleStiffness(
} }
} }
std::size_t contribution_count = 0U;
for (const auto& buffer : local_buffers) {
if (buffer.size() >
(std::numeric_limits<std::size_t>::max)() - contribution_count) {
return Result<SparseMatrix>::Failure(AssemblyFailure(
"invalid-assembly-dimensions", {}, std::to_string(elements.size()),
"Element contribution storage exceeds the addressable range."));
}
contribution_count += buffer.size();
}
std::vector<CooContribution> contributions; std::vector<CooContribution> contributions;
contributions.reserve(local_buffers.size() * kBeamContributionCount); contributions.reserve(contribution_count);
// Flatten only after all workers complete; workers never share CSR state. // Flatten only after all workers finish. This fixed element order remains
// independent of serial, reverse, or TBB task completion order.
for (const auto& buffer : local_buffers) { for (const auto& buffer : local_buffers) {
contributions.insert(contributions.end(), buffer.begin(), buffer.end()); contributions.insert(contributions.end(), buffer.begin(), buffer.end());
} }
@@ -300,4 +122,31 @@ Result<SparseMatrix> SparseAssembler::AssembleStiffness(
dofs.GetSparsePattern()); dofs.GetSparsePattern());
} }
Result<SparseMatrix> SparseAssembler::AssembleStiffness(
const AnalysisModel& model, const DofManager& dofs,
const ParallelFor& parallel_for) {
const Domain& domain = model.GetDomain();
std::vector<std::unique_ptr<Element>> owned_elements;
ElementView elements;
owned_elements.reserve(model.ActiveElements().size());
elements.reserve(model.ActiveElements().size());
const ElementFactory factory;
for (const EntityIndex element_index : model.ActiveElements()) {
if (element_index >= domain.Elements().Size()) {
return Result<SparseMatrix>::Failure(
AssemblyFailure("invalid-assembly-element", {domain.SourcePath(), 0U},
std::to_string(element_index),
"Active element index is outside the Domain."));
}
auto candidate = factory.Create(domain.Elements()[element_index], domain);
if (!candidate.HasValue()) {
return Result<SparseMatrix>::Failure(candidate.GetStatus());
}
owned_elements.push_back(std::move(candidate.Value()));
elements.push_back(std::cref(*owned_elements.back()));
}
return Assemble(elements, dofs, parallel_for);
}
} // namespace fesa } // namespace fesa
@@ -6,11 +6,13 @@
#include <array> #include <array>
#include <cstring> #include <cstring>
#include <filesystem> #include <filesystem>
#include <functional>
#include <utility> #include <utility>
#include <vector> #include <vector>
#include "fesa/analysis/analysis_model.h" #include "fesa/analysis/analysis_model.h"
#include "fesa/assembly/parallel_for.h" #include "fesa/assembly/parallel_for.h"
#include "fesa/elements/element.h"
#include "fesa/elements/mitc4_shell.h" #include "fesa/elements/mitc4_shell.h"
#include "fesa/fem/dof_manager.h" #include "fesa/fem/dof_manager.h"
#include "fesa/model/domain.h" #include "fesa/model/domain.h"
@@ -180,6 +182,42 @@ class ReverseParallelFor final : public fesa::ParallelFor {
mutable std::size_t observed_count_{0U}; mutable std::size_t observed_count_{0U};
}; };
class FakeStiffnessElement final : public fesa::Element {
public:
FakeStiffnessElement(fesa::ElementDofLayout layout, fesa::Matrix stiffness)
: layout_{std::move(layout)}, stiffness_{std::move(stiffness)} {}
const fesa::ElementDofLayout& DofLayout() const noexcept override {
return layout_;
}
fesa::Result<fesa::ElementStiffnessContribution> ComputeStiffness()
const override {
return fesa::Result<fesa::ElementStiffnessContribution>::Success(
{layout_, stiffness_});
}
fesa::Result<fesa::ElementResultBundle> Recover(
const fesa::Vector&) const override {
return fesa::Result<fesa::ElementResultBundle>::Success(
{layout_.source_id, fesa::BeamElementResultRows{}});
}
private:
fesa::ElementDofLayout layout_;
fesa::Matrix stiffness_;
};
fesa::Matrix MatrixFromRows(const std::vector<std::vector<double>>& rows) {
fesa::Matrix matrix{rows.size(), rows.empty() ? 0U : rows.front().size()};
for (std::size_t row = 0U; row < rows.size(); ++row) {
for (std::size_t column = 0U; column < rows[row].size(); ++column) {
matrix(row, column) = rows[row][column];
}
}
return matrix;
}
void ExpectByteIdentical(const fesa::SparseMatrix& actual, void ExpectByteIdentical(const fesa::SparseMatrix& actual,
const fesa::SparseMatrix& expected) { const fesa::SparseMatrix& expected) {
EXPECT_TRUE(ByteIdentical(actual.RowOffsets(), expected.RowOffsets())); EXPECT_TRUE(ByteIdentical(actual.RowOffsets(), expected.RowOffsets()));
@@ -187,6 +225,66 @@ void ExpectByteIdentical(const fesa::SparseMatrix& actual,
EXPECT_TRUE(ByteIdentical(actual.Values(), expected.Values())); EXPECT_TRUE(ByteIdentical(actual.Values(), expected.Values()));
} }
// C-ASSEMBLY-001
TEST(SparseAssembly,
AssemblesFakeRuntimeContributionsIntoExactDeterministicCsr) {
auto domain = fesa::Domain::Create(MakeDefinition());
ASSERT_TRUE(domain.HasValue());
auto model = fesa::AnalysisModel::Create(domain.Value());
ASSERT_TRUE(model.HasValue());
fesa::ElementDofLayout first_layout{
{"Beam-1", 10, "10"}, {0U, 2U}, {fesa::DofComponent::kUx}};
FakeStiffnessElement first{std::move(first_layout),
MatrixFromRows({{1.0, 2.0}, {3.0, 4.0}})};
fesa::ElementDofLayout second_layout{
{"Beam-1", 20, "20"},
{2U},
{fesa::DofComponent::kUx, fesa::DofComponent::kUz,
fesa::DofComponent::kUrz}};
FakeStiffnessElement second{
std::move(second_layout),
MatrixFromRows({{5.0, 6.0, 7.0}, {8.0, 9.0, 10.0}, {11.0, 12.0, 13.0}})};
const fesa::ElementView elements{std::cref(first), std::cref(second)};
fesa::DofManager dofs;
ASSERT_TRUE(dofs.Build(model.Value(), elements).IsOk());
fesa::SerialParallelFor serial_executor;
fesa::TbbParallelFor tbb_executor;
ReverseParallelFor reverse_executor;
auto serial =
fesa::SparseAssembler::Assemble(elements, dofs, serial_executor);
auto tbb = fesa::SparseAssembler::Assemble(elements, dofs, tbb_executor);
auto reversed =
fesa::SparseAssembler::Assemble(elements, dofs, reverse_executor);
ASSERT_TRUE(serial.HasValue());
ASSERT_TRUE(tbb.HasValue());
ASSERT_TRUE(reversed.HasValue());
EXPECT_EQ(reverse_executor.Calls(), 1U);
EXPECT_EQ(reverse_executor.ObservedCount(), elements.size());
EXPECT_EQ(serial.Value().Rows(), 18U);
EXPECT_EQ(serial.Value().Columns(), 18U);
EXPECT_EQ(serial.Value().RowOffsets(),
(std::vector<std::size_t>{0U, 2U, 2U, 2U, 2U, 2U, 2U, 2U, 2U, 2U,
2U, 2U, 2U, 6U, 6U, 9U, 9U, 9U, 12U}));
EXPECT_EQ(serial.Value().ColumnIndices(),
(std::vector<std::size_t>{0U, 12U, 0U, 12U, 14U, 17U, 12U, 14U, 17U,
12U, 14U, 17U}));
EXPECT_EQ(serial.Value().Values(),
(std::vector<double>{1.0, 2.0, 3.0, 9.0, 6.0, 7.0, 8.0, 9.0, 10.0,
11.0, 12.0, 13.0}));
ExpectByteIdentical(tbb.Value(), serial.Value());
ExpectByteIdentical(reversed.Value(), serial.Value());
for (std::size_t repetition = 0U; repetition < 8U; ++repetition) {
auto repeated =
fesa::SparseAssembler::Assemble(elements, dofs, tbb_executor);
ASSERT_TRUE(repeated.HasValue());
ExpectByteIdentical(repeated.Value(), serial.Value());
}
}
TEST(SparseAssembly, SerialTbbAndRepeatedRunsAreByteIdentical) { TEST(SparseAssembly, SerialTbbAndRepeatedRunsAreByteIdentical) {
auto domain_result = fesa::Domain::Create(MakeDefinition()); auto domain_result = fesa::Domain::Create(MakeDefinition());
ASSERT_TRUE(domain_result.HasValue()); ASSERT_TRUE(domain_result.HasValue());