feat(linear-static-3d-euler-beam): step 14 - dof-manager

2026-08-09 17:10:34 +09:00
parent a362d9938a
commit ac0e6b566e
6 changed files with 542 additions and 0 deletions
@@ -0,0 +1,230 @@
+#include "fesa/fem/dof_manager.hpp"
+
+#include <algorithm>
+#include <charconv>
+#include <stdexcept>
+#include <string>
+#include <system_error>
+#include <utility>
+
+namespace fesa {
+namespace {
+
+constexpr std::size_t dofsPerNode = 6U;
+
+char asciiLower(char value) {
+    if (value >= 'A' && value <= 'Z') {
+        return static_cast<char>(value + ('a' - 'A'));
+    }
+    return value;
+}
+
+bool equalName(const std::string& left, const std::string& right) {
+    return left.size() == right.size() &&
+           std::equal(
+               left.begin(), left.end(), right.begin(),
+               [](char leftValue, char rightValue) {
+                   return asciiLower(leftValue) == asciiLower(rightValue);
+               });
+}
+
+bool tryPositiveInteger(const std::string& text, std::int64_t& value) {
+    const char* const first = text.data();
+    const char* const last = first + text.size();
+    const auto parsed = std::from_chars(first, last, value);
+    return parsed.ec == std::errc{} && parsed.ptr == last && value > 0;
+}
+
+std::vector<EntityIndex> expandBoundaryTarget(
+    const Domain& domain, const BoundaryCondition& boundary) {
+    for (const auto& set : domain.nodeSets()) {
+        if (equalName(set.name, boundary.target)) {
+            return set.nodeIndices;
+        }
+    }
+
+    std::int64_t sourceLabel = 0;
+    if (tryPositiveInteger(boundary.target, sourceLabel)) {
+        for (std::size_t node = 0U; node < domain.nodes().size(); ++node) {
+            if (domain.nodes()[node].sourceId.sourceLabel == sourceLabel) {
+                return {static_cast<EntityIndex>(node)};
+            }
+        }
+    }
+    return {};
+}
+
+SparsePattern buildSparsePattern(
+    std::size_t fullDofCount,
+    const std::vector<EntityIndex>& activeElements,
+    const std::vector<std::array<std::size_t, 12>>& elementScatters) {
+    std::vector<std::vector<std::size_t>> columnsByRow(fullDofCount);
+    for (const EntityIndex element : activeElements) {
+        const auto& scatter = elementScatters.at(element);
+        for (const std::size_t row : scatter) {
+            auto& columns = columnsByRow[row];
+            columns.insert(columns.end(), scatter.begin(), scatter.end());
+        }
+    }
+
+    SparsePattern pattern;
+    pattern.rowOffsets.reserve(fullDofCount + 1U);
+    pattern.rowOffsets.push_back(0U);
+    for (auto& columns : columnsByRow) {
+        // Stable CSR structure is independent of element traversal duplicates.
+        std::sort(columns.begin(), columns.end());
+        columns.erase(std::unique(columns.begin(), columns.end()), columns.end());
+        pattern.columnIndices.insert(
+            pattern.columnIndices.end(), columns.begin(), columns.end());
+        pattern.rowOffsets.push_back(pattern.columnIndices.size());
+    }
+    return pattern;
+}
+
+} // namespace
+
+Result<DofManager> DofManager::create(const AnalysisModel& model) {
+    const Domain& domain = model.domain();
+    const std::size_t fullCount = domain.nodes().size() * dofsPerNode;
+
+    std::vector<std::optional<double>> prescribedByFullDof(fullCount);
+    for (const EntityIndex boundaryIndex : model.activeBoundaryConditions()) {
+        const auto& boundary = model.step().boundaries.at(boundaryIndex);
+        const auto target = expandBoundaryTarget(domain, boundary);
+        for (const EntityIndex node : target) {
+            for (int component = boundary.firstDof;
+                 component <= boundary.lastDof;
+                 ++component) {
+                const std::size_t fullDof =
+                    static_cast<std::size_t>(node) * dofsPerNode +
+                    static_cast<std::size_t>(component - 1);
+                auto& prescribed = prescribedByFullDof[fullDof];
+                if (prescribed && *prescribed != boundary.value) {
+                    return Result<DofManager>::failure(Status::failure(
+                        FailureCategory::input,
+                        {{Severity::error,
+                          "conflicting-boundary-condition",
+                          boundary.location,
+                          "BOUNDARY",
+                          boundary.target,
+                          "Expanded boundary rows prescribe different values to one node/DOF."}}));
+                }
+                prescribed = boundary.value;
+            }
+        }
+    }
+
+    std::vector<std::size_t> freeDofs;
+    std::vector<std::size_t> constrainedDofs;
+    std::vector<double> constrainedValues;
+    std::vector<std::optional<std::size_t>> freeEquations(fullCount);
+    freeDofs.reserve(fullCount);
+    constrainedDofs.reserve(fullCount);
+    constrainedValues.reserve(fullCount);
+    // A full-DOF scan fixes free equations, constrained DOFs, and dc in the
+    // same stable order regardless of boundary declaration overlap.
+    for (std::size_t fullDof = 0U; fullDof < fullCount; ++fullDof) {
+        if (prescribedByFullDof[fullDof]) {
+            constrainedDofs.push_back(fullDof);
+            constrainedValues.push_back(*prescribedByFullDof[fullDof]);
+        } else {
+            freeEquations[fullDof] = freeDofs.size();
+            freeDofs.push_back(fullDof);
+        }
+    }
+
+    Vector prescribedValues{constrainedValues.size()};
+    for (std::size_t index = 0U; index < constrainedValues.size(); ++index) {
+        prescribedValues[index] = constrainedValues[index];
+    }
+
+    std::vector<std::array<std::size_t, 12>> elementScatters(
+        domain.elements().size());
+    for (const EntityIndex elementIndex : model.activeElements()) {
+        const auto& element = domain.elements().at(elementIndex);
+        auto& scatter = elementScatters.at(elementIndex);
+        for (std::size_t endpoint = 0U; endpoint < element.nodeIndices.size(); ++endpoint) {
+            const std::size_t node = element.nodeIndices[endpoint];
+            for (std::size_t component = 0U; component < dofsPerNode; ++component) {
+                scatter[endpoint * dofsPerNode + component] =
+                    node * dofsPerNode + component;
+            }
+        }
+    }
+
+    auto pattern = buildSparsePattern(
+        fullCount, model.activeElements(), elementScatters);
+    return Result<DofManager>::success(DofManager{
+        fullCount,
+        std::move(freeEquations),
+        std::move(elementScatters),
+        std::move(freeDofs),
+        std::move(constrainedDofs),
+        std::move(prescribedValues),
+        std::move(pattern)});
+}
+
+std::size_t DofManager::fullDofCount() const noexcept {
+    return fullDofCount_;
+}
+
+std::size_t DofManager::freeDofCount() const noexcept {
+    return freeDofs_.size();
+}
+
+std::size_t DofManager::constrainedDofCount() const noexcept {
+    return constrainedDofs_.size();
+}
+
+std::size_t DofManager::fullDof(
+    EntityIndex node, DofComponent component) const {
+    const std::size_t componentIndex = static_cast<std::size_t>(component);
+    if (node >= fullDofCount_ / dofsPerNode || componentIndex >= dofsPerNode) {
+        throw std::out_of_range{"Node or DOF component is out of range."};
+    }
+    return static_cast<std::size_t>(node) * dofsPerNode + componentIndex;
+}
+
+std::optional<std::size_t> DofManager::freeEquation(
+    std::size_t fullDof) const {
+    return freeEquations_.at(fullDof);
+}
+
+const std::array<std::size_t, 12>& DofManager::elementScatter(
+    EntityIndex element) const {
+    return elementScatters_.at(element);
+}
+
+const std::vector<std::size_t>& DofManager::freeDofs() const noexcept {
+    return freeDofs_;
+}
+
+const std::vector<std::size_t>& DofManager::constrainedDofs() const noexcept {
+    return constrainedDofs_;
+}
+
+const Vector& DofManager::prescribedValues() const noexcept {
+    return prescribedValues_;
+}
+
+const SparsePattern& DofManager::sparsePattern() const noexcept {
+    return sparsePattern_;
+}
+
+DofManager::DofManager(
+    std::size_t fullDofCount,
+    std::vector<std::optional<std::size_t>> freeEquations,
+    std::vector<std::array<std::size_t, 12>> elementScatters,
+    std::vector<std::size_t> freeDofs,
+    std::vector<std::size_t> constrainedDofs,
+    Vector prescribedValues,
+    SparsePattern sparsePattern)
+    : fullDofCount_{fullDofCount},
+      freeEquations_{std::move(freeEquations)},
+      elementScatters_{std::move(elementScatters)},
+      freeDofs_{std::move(freeDofs)},
+      constrainedDofs_{std::move(constrainedDofs)},
+      prescribedValues_{std::move(prescribedValues)},
+      sparsePattern_{std::move(sparsePattern)} {}
+
+} // namespace fesa