feat(linear-static-3d-euler-beam): step 18 - sparse-assembly

src/fesa/CMakeLists.txt

@@ -4,6 +4,7 @@ add_library(
     analysis/analysis_model.cpp
     analysis/analysis_state.cpp
     assembly/parallel_for.cpp
+    assembly/sparse_assembler.cpp
     build_info.cpp
     core/diagnostic.cpp
     core/status.cpp
@@ -12,6 +13,7 @@ add_library(
     io/abaqus/domain_mapper.cpp
     io/abaqus/input_reader.cpp
     math/matrix.cpp
+    math/sparse_matrix.cpp
     math/vector.cpp
     model/domain.cpp
 )

src/fesa/assembly/sparse_assembler.cpp

@@ -0,0 +1,183 @@
+#include "fesa/assembly/sparse_assembler.hpp"
+
+#include "fesa/analysis/analysis_model.hpp"
+#include "fesa/assembly/parallel_for.hpp"
+#include "fesa/elements/euler_beam_3d.hpp"
+#include "fesa/fem/dof_manager.hpp"
+
+#include <array>
+#include <limits>
+#include <optional>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace fesa {
+namespace {
+
+constexpr std::size_t kDofsPerNode = 6U;
+constexpr std::size_t kElementDofCount = 12U;
+constexpr std::size_t kContributionCount =
+    kElementDofCount * kElementDofCount;
+
+using ElementBuffer = std::array<CooContribution, kContributionCount>;
+
+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::model,
+        {{Severity::error,
+          code,
+          location,
+          "*ELEMENT",
+          identity,
+          message}}));
+}
+
+} // namespace
+
+Result<SparseMatrix> SparseAssembler::assembleStiffness(
+    const AnalysisModel& model,
+    const DofManager& dofs,
+    const ParallelFor& parallelFor) {
+    const Domain& domain = model.domain();
+    if (domain.nodes().size() >
+        (std::numeric_limits<std::size_t>::max)() / kDofsPerNode ||
+        dofs.fullDofCount() != domain.nodes().size() * kDofsPerNode) {
+        return assemblyFailure(
+            "invalid-assembly-dimensions",
+            {domain.sourcePath(), 0U},
+            std::to_string(dofs.fullDofCount()),
+            "DofManager dimensions do not match the active model nodes.");
+    }
+    if (model.activeElements().size() >
+        (std::numeric_limits<std::size_t>::max)() / kContributionCount) {
+        return assemblyFailure(
+            "invalid-assembly-dimensions",
+            {domain.sourcePath(), 0U},
+            std::to_string(model.activeElements().size()),
+            "Element contribution storage exceeds the addressable range.");
+    }
+
+    std::vector<std::array<std::size_t, kElementDofCount>> scatters;
+    scatters.reserve(model.activeElements().size());
+    for (const EntityIndex elementIndex : model.activeElements()) {
+        if (elementIndex >= domain.elements().size()) {
+            return assemblyFailure(
+                "invalid-assembly-element",
+                {domain.sourcePath(), 0U},
+                std::to_string(elementIndex),
+                "Active element index is outside the Domain.");
+        }
+        const auto& element = domain.elements()[elementIndex];
+        if (element.nodeIndices[0U] >= domain.nodes().size() ||
+            element.nodeIndices[1U] >= domain.nodes().size() ||
+            element.materialIndex >= domain.materials().size() ||
+            element.sectionIndex >= domain.sections().size()) {
+            return assemblyFailure(
+                "invalid-assembly-element",
+                element.location,
+                element.sourceId.sourceLabelText,
+                "Element references an entity outside the Domain.");
+        }
+
+        std::array<std::size_t, kElementDofCount> scatter{};
+        try {
+            scatter = dofs.elementScatter(elementIndex);
+        } catch (const std::out_of_range&) {
+            return assemblyFailure(
+                "invalid-assembly-scatter",
+                element.location,
+                element.sourceId.sourceLabelText,
+                "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.nodeIndices[endpoint]) *
+                        kDofsPerNode +
+                    component;
+                if (scatter[local] != expected ||
+                    scatter[local] >= dofs.fullDofCount()) {
+                    return assemblyFailure(
+                        "invalid-assembly-scatter",
+                        element.location,
+                        element.sourceId.sourceLabelText,
+                        "Element scatter does not match the active model topology.");
+                }
+            }
+        }
+        scatters.push_back(scatter);
+    }
+
+    std::vector<ElementBuffer> localBuffers(model.activeElements().size());
+    std::vector<std::optional<Status>> localFailures(
+        model.activeElements().size());
+    parallelFor.execute(
+        model.activeElements().size(),
+        [&](const std::size_t elementOrder) {
+            const EntityIndex elementIndex = model.activeElements()[elementOrder];
+            const auto& definition = domain.elements()[elementIndex];
+            const auto beam = EulerBeam3D::create(
+                domain.nodes()[definition.nodeIndices[0U]],
+                domain.nodes()[definition.nodeIndices[1U]],
+                domain.sections()[definition.sectionIndex],
+                domain.materials()[definition.materialIndex]);
+            if (!beam.hasValue()) {
+                localFailures[elementOrder] = beam.status();
+                return;
+            }
+
+            const Matrix stiffness = beam.value().globalStiffness();
+            auto& buffer = localBuffers[elementOrder];
+            const auto& scatter = scatters[elementOrder];
+            for (std::size_t localRow = 0U;
+                 localRow < kElementDofCount;
+                 ++localRow) {
+                for (std::size_t localColumn = 0U;
+                     localColumn < kElementDofCount;
+                     ++localColumn) {
+                    const std::size_t localOrder =
+                        localRow * kElementDofCount + localColumn;
+                    buffer[localOrder] = {
+                        scatter[localRow],
+                        scatter[localColumn],
+                        stiffness(localRow, localColumn),
+                        elementOrder,
+                        localOrder};
+                }
+            }
+        });
+
+    for (std::size_t elementOrder = 0U;
+         elementOrder < localFailures.size();
+         ++elementOrder) {
+        if (localFailures[elementOrder]) {
+            return Result<SparseMatrix>::failure(
+                *localFailures[elementOrder]);
+        }
+    }
+
+    std::vector<CooContribution> contributions;
+    contributions.reserve(
+        localBuffers.size() * kContributionCount);
+    // Flatten only after all workers complete; workers never share CSR state.
+    for (const auto& buffer : localBuffers) {
+        contributions.insert(
+            contributions.end(), buffer.begin(), buffer.end());
+    }
+    return SparseMatrix::fromCoo(
+        dofs.fullDofCount(),
+        dofs.fullDofCount(),
+        std::move(contributions),
+        dofs.sparsePattern());
+}
+
+} // namespace fesa

src/fesa/math/sparse_matrix.cpp

@@ -0,0 +1,239 @@
+#include "fesa/math/sparse_matrix.hpp"
+
+#include "fesa/fem/dof_manager.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+#include <stdexcept>
+#include <string>
+#include <tuple>
+#include <utility>
+
+namespace fesa {
+namespace {
+
+Status sparseFailure(
+    const std::string& code,
+    const std::string& identity,
+    const std::string& message) {
+    return Status::failure(
+        FailureCategory::model,
+        {{Severity::error,
+          code,
+          {{}, 0U},
+          "SPARSE_MATRIX",
+          identity,
+          message}});
+}
+
+Status validateCsr(
+    const std::size_t rows,
+    const std::size_t columns,
+    const std::vector<std::size_t>& rowOffsets,
+    const std::vector<std::size_t>& columnIndices,
+    const std::vector<double>* const values) {
+    if (rows == (std::numeric_limits<std::size_t>::max)() ||
+        rowOffsets.size() != rows + 1U) {
+        return sparseFailure(
+            "invalid-sparse-shape",
+            "row-offset-count",
+            "CSR row offsets must contain exactly rows plus one entries.");
+    }
+    if (rowOffsets.empty() || rowOffsets.front() != 0U ||
+        rowOffsets.back() != columnIndices.size()) {
+        return sparseFailure(
+            "invalid-sparse-pattern",
+            "row-offset-range",
+            "CSR row offsets must start at zero and end at the column count.");
+    }
+    if (values != nullptr && values->size() != columnIndices.size()) {
+        return sparseFailure(
+            "invalid-sparse-shape",
+            "value-count",
+            "CSR column and value arrays must have equal sizes.");
+    }
+
+    for (std::size_t row = 0U; row < rows; ++row) {
+        const std::size_t begin = rowOffsets[row];
+        const std::size_t end = rowOffsets[row + 1U];
+        if (begin > end || end > columnIndices.size()) {
+            return sparseFailure(
+                "invalid-sparse-pattern",
+                std::to_string(row),
+                "CSR row offsets must be nondecreasing and remain in range.");
+        }
+        for (std::size_t position = begin; position < end; ++position) {
+            if (columnIndices[position] >= columns) {
+                return sparseFailure(
+                    "invalid-sparse-index",
+                    std::to_string(position),
+                    "CSR column index is outside the matrix dimensions.");
+            }
+            if (position > begin &&
+                columnIndices[position - 1U] >= columnIndices[position]) {
+                return sparseFailure(
+                    "invalid-sparse-pattern",
+                    std::to_string(row),
+                    "CSR columns must be sorted and unique within each row.");
+            }
+            if (values != nullptr && !std::isfinite((*values)[position])) {
+                return sparseFailure(
+                    "nonfinite-sparse-value",
+                    std::to_string(position),
+                    "CSR values must be finite.");
+            }
+        }
+    }
+    return Status::ok();
+}
+
+} // namespace
+
+Result<SparseMatrix> SparseMatrix::fromCoo(
+    const std::size_t rows,
+    const std::size_t columns,
+    std::vector<CooContribution> contributions,
+    const SparsePattern& expectedPattern) {
+    const Status patternStatus = validateCsr(
+        rows,
+        columns,
+        expectedPattern.rowOffsets,
+        expectedPattern.columnIndices,
+        nullptr);
+    if (!patternStatus.isOk()) {
+        return Result<SparseMatrix>::failure(patternStatus);
+    }
+
+    for (const auto& contribution : contributions) {
+        if (contribution.row >= rows || contribution.column >= columns) {
+            return Result<SparseMatrix>::failure(sparseFailure(
+                "invalid-sparse-index",
+                std::to_string(contribution.row) + ":" +
+                    std::to_string(contribution.column),
+                "COO contribution index is outside the matrix dimensions."));
+        }
+        if (!std::isfinite(contribution.value)) {
+            return Result<SparseMatrix>::failure(sparseFailure(
+                "nonfinite-sparse-value",
+                std::to_string(contribution.elementOrder) + ":" +
+                    std::to_string(contribution.localOrder),
+                "COO contribution values must be finite."));
+        }
+    }
+
+    // The complete tuple fixes duplicate summation order independently of
+    // worker completion order. stable_sort also preserves exact tuple ties.
+    std::stable_sort(
+        contributions.begin(),
+        contributions.end(),
+        [](const CooContribution& left, const CooContribution& right) {
+            return std::tie(
+                       left.row,
+                       left.column,
+                       left.elementOrder,
+                       left.localOrder) <
+                   std::tie(
+                       right.row,
+                       right.column,
+                       right.elementOrder,
+                       right.localOrder);
+        });
+
+    std::vector<double> values(expectedPattern.columnIndices.size(), 0.0);
+    for (const auto& contribution : contributions) {
+        const std::size_t begin = expectedPattern.rowOffsets[contribution.row];
+        const std::size_t end = expectedPattern.rowOffsets[contribution.row + 1U];
+        const auto first = expectedPattern.columnIndices.begin() + begin;
+        const auto last = expectedPattern.columnIndices.begin() + end;
+        const auto found = std::lower_bound(first, last, contribution.column);
+        if (found == last || *found != contribution.column) {
+            return Result<SparseMatrix>::failure(sparseFailure(
+                "sparse-pattern-mismatch",
+                std::to_string(contribution.row) + ":" +
+                    std::to_string(contribution.column),
+                "COO contribution is absent from the expected sparse pattern."));
+        }
+
+        const std::size_t position = static_cast<std::size_t>(
+            std::distance(expectedPattern.columnIndices.begin(), found));
+        values[position] += contribution.value;
+        if (!std::isfinite(values[position])) {
+            return Result<SparseMatrix>::failure(sparseFailure(
+                "nonfinite-sparse-value",
+                std::to_string(contribution.row) + ":" +
+                    std::to_string(contribution.column),
+                "Ordered COO duplicate summation produced a nonfinite value."));
+        }
+    }
+
+    SparseMatrix matrix{
+        rows,
+        columns,
+        expectedPattern.rowOffsets,
+        expectedPattern.columnIndices,
+        std::move(values)};
+    const Status status = matrix.validate();
+    if (!status.isOk()) {
+        return Result<SparseMatrix>::failure(status);
+    }
+    return Result<SparseMatrix>::success(std::move(matrix));
+}
+
+std::size_t SparseMatrix::rows() const noexcept {
+    return rows_;
+}
+
+std::size_t SparseMatrix::columns() const noexcept {
+    return columns_;
+}
+
+const std::vector<std::size_t>& SparseMatrix::rowOffsets() const noexcept {
+    return rowOffsets_;
+}
+
+const std::vector<std::size_t>& SparseMatrix::columnIndices() const noexcept {
+    return columnIndices_;
+}
+
+const std::vector<double>& SparseMatrix::values() const noexcept {
+    return values_;
+}
+
+Vector SparseMatrix::multiply(const Vector& rhs) const {
+    if (columns_ != rhs.size()) {
+        throw std::invalid_argument{
+            "Sparse matrix-vector multiplication has incompatible dimensions."};
+    }
+
+    Vector result{rows_};
+    for (std::size_t row = 0U; row < rows_; ++row) {
+        double value = 0.0;
+        for (std::size_t position = rowOffsets_[row];
+             position < rowOffsets_[row + 1U];
+             ++position) {
+            value += values_[position] * rhs[columnIndices_[position]];
+        }
+        result[row] = value;
+    }
+    return result;
+}
+
+Status SparseMatrix::validate() const {
+    return validateCsr(
+        rows_, columns_, rowOffsets_, columnIndices_, &values_);
+}
+
+SparseMatrix::SparseMatrix(
+    const std::size_t rows,
+    const std::size_t columns,
+    std::vector<std::size_t> rowOffsets,
+    std::vector<std::size_t> columnIndices,
+    std::vector<double> values)
+    : rows_{rows},
+      columns_{columns},
+      rowOffsets_{std::move(rowOffsets)},
+      columnIndices_{std::move(columnIndices)},
+      values_{std::move(values)} {}
+
+} // namespace fesa