feat: add assembly reduced solver boundary

2026-05-04 23:35:41 +09:00
parent ebdb2519dc
commit d373969732
5 changed files with 328 additions and 22 deletions
@@ -1083,6 +1083,52 @@ class DofManager {
  std::vector<LocalIndex> constrained_full_indices_;
 };

+struct SparsePatternEntry {
+  EquationId row = 0;
+  EquationId col = 0;
+};
+
+struct SparsePattern {
+  EquationId equation_count = 0;
+  std::vector<SparsePatternEntry> entries;
+
+  SparseIndex nonzeroCount() const {
+    return static_cast<SparseIndex>(entries.size());
+  }
+
+  bool contains(EquationId row, EquationId col) const {
+    return std::any_of(entries.begin(), entries.end(), [&](const SparsePatternEntry& entry) {
+      return entry.row == row && entry.col == col;
+    });
+  }
+};
+
+inline SparsePattern buildReducedSparsePattern(const Domain& domain, const DofManager& dofs) {
+  SparsePattern pattern;
+  pattern.equation_count = dofs.freeDofCount();
+  std::set<std::pair<EquationId, EquationId>> ordered_entries;
+  for (const auto& [element_id, element] : domain.elements) {
+    (void)element_id;
+    const auto equations = dofs.elementEquationIds(element);
+    for (EquationId row : equations) {
+      if (row < 0) {
+        continue;
+      }
+      for (EquationId col : equations) {
+        if (col < 0) {
+          continue;
+        }
+        ordered_entries.insert({row, col});
+      }
+    }
+  }
+  pattern.entries.reserve(ordered_entries.size());
+  for (const auto& entry : ordered_entries) {
+    pattern.entries.push_back({entry.first, entry.second});
+  }
+  return pattern;
+}
+
 class DenseMatrix {
 public:
  DenseMatrix() = default;
@@ -2102,12 +2148,34 @@ class MITC4ElementKernel {
 struct AssemblyResult {
  DenseMatrix k_full;
  std::vector<Real> f_full;
+  SparsePattern reduced_pattern;
  std::vector<Diagnostic> diagnostics;
+
+  bool ok() const {
+    return !hasError(diagnostics);
+  }
+};
+
+struct ReducedSystem {
+  DenseMatrix k;
+  std::vector<Real> f;
+  std::vector<LocalIndex> free_full_indices;
+  std::vector<Diagnostic> diagnostics;
+
+  bool ok() const {
+    return !hasError(diagnostics);
+  }
 };

 inline AssemblyResult assembleSystem(const Domain& domain, const DofManager& dofs, ElementStiffnessOptions options = {}) {
-  AssemblyResult result{DenseMatrix(dofs.fullDofCount(), dofs.fullDofCount()), std::vector<Real>(static_cast<std::size_t>(dofs.fullDofCount()), 0.0), {}};
-  MITC4ElementKernel kernel;
+  AssemblyResult result;
+  result.k_full = DenseMatrix(dofs.fullDofCount(), dofs.fullDofCount());
+  result.f_full = std::vector<Real>(static_cast<std::size_t>(dofs.fullDofCount()), 0.0);
+  result.reduced_pattern = buildReducedSparsePattern(domain, dofs);
+  if (dofs.freeDofCount() > 0 && result.reduced_pattern.nonzeroCount() == 0) {
+    result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SINGULAR-SPARSE-PATTERN",
+                                                "Reduced sparse pattern has no stiffness entries", "assembly"));
+  }
  for (const auto& [element_id, element] : domain.elements) {
    const ShellSection* section = shellSectionForElement(domain, element_id);
    if (section == nullptr) {
@@ -2123,19 +2191,63 @@ inline AssemblyResult assembleSystem(const Domain& domain, const DofManager& dof
    for (std::size_t i = 0; i < 4; ++i) {
      coordinates[i] = domain.nodes.at(element.node_ids[i]).coordinates;
    }
-    DenseMatrix ke = kernel.stiffness(coordinates, material_it->second.elastic_modulus, material_it->second.poisson_ratio, section->thickness, options);
+    const auto stiffness = mitc4ElementStiffness(coordinates, material_it->second.elastic_modulus,
+                                                material_it->second.poisson_ratio, section->thickness, options);
+    result.diagnostics.insert(result.diagnostics.end(), stiffness.diagnostics.begin(), stiffness.diagnostics.end());
+    if (!stiffness.ok()) {
+      continue;
+    }
    const auto element_full_indices = dofs.elementFullDofIndices(element);
    for (LocalIndex a = 0; a < 24; ++a) {
      const LocalIndex ia = element_full_indices[static_cast<std::size_t>(a)];
      for (LocalIndex b = 0; b < 24; ++b) {
        const LocalIndex ib = element_full_indices[static_cast<std::size_t>(b)];
-        result.k_full.add(ia, ib, ke(a, b));
+        result.k_full.add(ia, ib, stiffness.global(a, b));
      }
    }
  }
  for (const NodalLoad& load : domain.loads) {
+    const auto dof = dofFromAbaqus(load.dof);
+    if (!dof) {
+      result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-ASSEMBLY-LOAD-DOF",
+                                                  "Nodal load references an invalid DOF", "cload"));
+      continue;
+    }
    for (GlobalId node_id : resolveNodeTarget(domain, load.target, &result.diagnostics)) {
-      result.f_full[static_cast<std::size_t>(dofs.fullIndex(node_id, *dofFromAbaqus(load.dof)))] += load.magnitude;
+      result.f_full[static_cast<std::size_t>(dofs.fullIndex(node_id, *dof))] += load.magnitude;
+    }
+  }
+  return result;
+}
+
+inline ReducedSystem projectToReducedSystem(const AssemblyResult& assembly, const DofManager& dofs) {
+  ReducedSystem result;
+  result.k = DenseMatrix(dofs.freeDofCount(), dofs.freeDofCount());
+  result.f = std::vector<Real>(static_cast<std::size_t>(dofs.freeDofCount()), 0.0);
+  result.free_full_indices = dofs.freeFullIndices();
+  if (assembly.k_full.rows() != assembly.k_full.cols() ||
+      assembly.k_full.rows() != dofs.fullDofCount() ||
+      static_cast<LocalIndex>(assembly.f_full.size()) != dofs.fullDofCount()) {
+    result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-ASSEMBLY-SIZE",
+                                                "Full-space stiffness/load sizes do not match DofManager", "assembly"));
+    return result;
+  }
+  if (dofs.freeDofCount() == 0) {
+    result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SINGULAR-NO-FREE-DOFS",
+                                                "No free DOFs exist after applying constraints", "dof"));
+    return result;
+  }
+  if (assembly.reduced_pattern.equation_count != dofs.freeDofCount() || assembly.reduced_pattern.nonzeroCount() == 0) {
+    result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SINGULAR-SPARSE-PATTERN",
+                                                "Reduced sparse pattern is empty or inconsistent with free DOFs", "assembly"));
+    return result;
+  }
+  for (LocalIndex i = 0; i < dofs.freeDofCount(); ++i) {
+    const LocalIndex full_i = dofs.freeFullIndices()[static_cast<std::size_t>(i)];
+    result.f[static_cast<std::size_t>(i)] = assembly.f_full[static_cast<std::size_t>(full_i)];
+    for (LocalIndex j = 0; j < dofs.freeDofCount(); ++j) {
+      const LocalIndex full_j = dofs.freeFullIndices()[static_cast<std::size_t>(j)];
+      result.k(i, j) = assembly.k_full(full_i, full_j);
    }
  }
  return result;
@@ -2237,6 +2349,7 @@ class InMemoryResultsWriter {
 struct AnalysisState {
  std::vector<Real> u_full;
  std::vector<Real> f_external_full;
+  std::vector<Real> f_internal_full;
  std::vector<Real> reaction_full;
  bool converged = false;
 };
@@ -2273,6 +2386,9 @@ class Analysis {
 };

 class LinearStaticAnalysis final : public Analysis {
+ public:
+  explicit LinearStaticAnalysis(const LinearSolver* solver = nullptr) : solver_(solver) {}
+
 protected:
  void solve(const Domain& domain, AnalysisResult& result) const override {
    DofManager dofs(domain);
@@ -2286,30 +2402,39 @@ class LinearStaticAnalysis final : public Analysis {
    if (hasError(result.diagnostics)) {
      return;
    }
-    DenseMatrix k_reduced(dofs.freeDofCount(), dofs.freeDofCount());
-    std::vector<Real> f_reduced(static_cast<std::size_t>(dofs.freeDofCount()), 0.0);
-    for (LocalIndex i = 0; i < dofs.freeDofCount(); ++i) {
-      const LocalIndex full_i = dofs.freeFullIndices()[static_cast<std::size_t>(i)];
-      f_reduced[static_cast<std::size_t>(i)] = assembly.f_full[static_cast<std::size_t>(full_i)];
-      for (LocalIndex j = 0; j < dofs.freeDofCount(); ++j) {
-        const LocalIndex full_j = dofs.freeFullIndices()[static_cast<std::size_t>(j)];
-        k_reduced(i, j) = assembly.k_full(full_i, full_j);
-      }
+    const auto reduced = projectToReducedSystem(assembly, dofs);
+    result.diagnostics.insert(result.diagnostics.end(), reduced.diagnostics.begin(), reduced.diagnostics.end());
+    if (hasError(result.diagnostics)) {
+      return;
    }
-    GaussianEliminationSolver solver;
-    SolveResult solved = solver.solve(k_reduced, f_reduced);
+    const LinearSolver& active_solver = solver_ == nullptr ? defaultSolver() : *solver_;
+    SolveResult solved = active_solver.solve(reduced.k, reduced.f);
    result.diagnostics.insert(result.diagnostics.end(), solved.diagnostics.begin(), solved.diagnostics.end());
    if (!solved.ok()) {
      return;
    }
+    if (static_cast<LocalIndex>(solved.x.size()) != dofs.freeDofCount()) {
+      result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SOLVER-SIZE",
+                                                  "Linear solver returned a vector with the wrong size", "solver"));
+      return;
+    }
    result.state.u_full = dofs.reconstructFullVector(solved.x);
    result.state.f_external_full = assembly.f_full;
+    result.state.f_internal_full = assembly.k_full.multiply(result.state.u_full);
    result.state.reaction_full = recoverFullReaction(assembly.k_full, result.state.u_full, result.state.f_external_full);
    result.state.converged = true;
    InMemoryResultsWriter writer;
    writer.writeLinearStatic(domain, dofs, result.state.u_full, result.state.reaction_full);
    result.result_file = writer.result();
  }
+
+ private:
+  static const LinearSolver& defaultSolver() {
+    static const GaussianEliminationSolver solver;
+    return solver;
+  }
+
+  const LinearSolver* solver_ = nullptr;
 };

 struct CsvDisplacementRow {