refactor: extract core domain dof modules

include/fesa/Boundary/Boundary.hpp

@@ -1,9 +1,23 @@
 #pragma once
 
+#include "fesa/Core/Types.hpp"
 #include "fesa/ModuleInfo.hpp"
 
+#include <string>
+
 namespace fesa::module {
 
 inline constexpr std::string_view kBoundary = "Boundary";
 
 }  // namespace fesa::module
+
+namespace fesa {
+
+struct BoundaryCondition {
+  std::string target;
+  int first_dof = 0;
+  int last_dof = 0;
+  Real magnitude = 0.0;
+};
+
+}  // namespace fesa

include/fesa/Core/AnalysisModel.hpp

@@ -0,0 +1,67 @@
+#pragma once
+
+#include "fesa/Core/Domain.hpp"
+#include "fesa/Util/Diagnostics.hpp"
+
+#include <cstddef>
+#include <string>
+#include <vector>
+
+namespace fesa {
+
+struct AnalysisModel {
+  StepDefinition step;
+  std::vector<GlobalId> active_element_ids;
+  std::vector<std::size_t> active_boundary_condition_indices;
+  std::vector<std::size_t> active_load_indices;
+  std::vector<std::size_t> active_shell_section_indices;
+  std::vector<std::string> active_material_keys;
+  std::vector<Diagnostic> diagnostics;
+
+  bool ok() const {
+    return !hasError(diagnostics);
+  }
+};
+
+inline AnalysisModel buildLinearStaticAnalysisModel(const Domain& domain, LocalIndex step_index = 0) {
+  AnalysisModel model;
+  if (domain.steps.empty()) {
+    model.step = {"Step-1", "linear_static"};
+  } else {
+    if (step_index < 0 || step_index >= static_cast<LocalIndex>(domain.steps.size())) {
+      model.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-ANALYSIS-STEP-INDEX",
+                                                 "Requested analysis step index is out of range", "analysis model"));
+      model.step = domain.steps.front();
+    } else {
+      model.step = domain.steps[static_cast<std::size_t>(step_index)];
+    }
+  }
+  if (domain.steps.size() > 1) {
+    model.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-ANALYSIS-MULTIPLE-STEPS",
+                                               "Phase 1 execution supports one active linear static step", "analysis model"));
+  }
+  if (model.step.analysis_type != "linear_static") {
+    model.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-ANALYSIS-UNSUPPORTED-STEP",
+                                               "Only linear static steps are supported in Phase 1", "analysis model"));
+  }
+  for (const auto& [element_id, element] : domain.elements) {
+    (void)element;
+    model.active_element_ids.push_back(element_id);
+  }
+  for (std::size_t i = 0; i < domain.boundary_conditions.size(); ++i) {
+    model.active_boundary_condition_indices.push_back(i);
+  }
+  for (std::size_t i = 0; i < domain.loads.size(); ++i) {
+    model.active_load_indices.push_back(i);
+  }
+  for (std::size_t i = 0; i < domain.shell_sections.size(); ++i) {
+    model.active_shell_section_indices.push_back(i);
+  }
+  for (const auto& [material_key, material] : domain.materials) {
+    (void)material;
+    model.active_material_keys.push_back(material_key);
+  }
+  return model;
+}
+
+}  // namespace fesa

include/fesa/Core/AnalysisState.hpp

@@ -0,0 +1,17 @@
+#pragma once
+
+#include "fesa/Core/Types.hpp"
+
+#include <vector>
+
+namespace fesa {
+
+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;
+};
+
+}  // namespace fesa

include/fesa/Core/Core.hpp

@@ -1,5 +1,12 @@
 #pragma once
 
+#include "fesa/Core/AnalysisModel.hpp"
+#include "fesa/Core/AnalysisState.hpp"
+#include "fesa/Core/Dof.hpp"
+#include "fesa/Core/DofManager.hpp"
+#include "fesa/Core/Domain.hpp"
+#include "fesa/Core/Types.hpp"
+#include "fesa/Core/Validation.hpp"
 #include "fesa/ModuleInfo.hpp"
 
 namespace fesa::module {

include/fesa/Core/Dof.hpp

@@ -0,0 +1,59 @@
+#pragma once
+
+#include "fesa/Core/Types.hpp"
+
+#include <array>
+#include <optional>
+#include <string>
+#include <vector>
+
+namespace fesa {
+
+enum class Dof : int { UX = 0, UY = 1, UZ = 2, RX = 3, RY = 4, RZ = 5 };
+
+inline std::array<Dof, 6> allDofs() {
+  return {Dof::UX, Dof::UY, Dof::UZ, Dof::RX, Dof::RY, Dof::RZ};
+}
+
+inline int dofIndex(Dof dof) {
+  return static_cast<int>(dof);
+}
+
+inline int abaqusDofNumber(Dof dof) {
+  return dofIndex(dof) + 1;
+}
+
+inline std::optional<Dof> dofFromAbaqus(int dof) {
+  if (dof < 1 || dof > 6) {
+    return std::nullopt;
+  }
+  return static_cast<Dof>(dof - 1);
+}
+
+inline const char* dofLabel(Dof dof) {
+  switch (dof) {
+    case Dof::UX:
+      return "UX";
+    case Dof::UY:
+      return "UY";
+    case Dof::UZ:
+      return "UZ";
+    case Dof::RX:
+      return "RX";
+    case Dof::RY:
+      return "RY";
+    case Dof::RZ:
+      return "RZ";
+  }
+  return "";
+}
+
+inline std::vector<std::string> displacementComponentLabels() {
+  return {"UX", "UY", "UZ", "RX", "RY", "RZ"};
+}
+
+inline std::vector<std::string> reactionComponentLabels() {
+  return {"RFX", "RFY", "RFZ", "RMX", "RMY", "RMZ"};
+}
+
+}  // namespace fesa

include/fesa/Core/DofManager.hpp

@@ -0,0 +1,145 @@
+#pragma once
+
+#include "fesa/Core/Dof.hpp"
+#include "fesa/Core/Domain.hpp"
+#include "fesa/Core/Validation.hpp"
+
+#include <array>
+#include <map>
+#include <set>
+#include <utility>
+#include <vector>
+
+namespace fesa {
+
+struct DofAddress {
+  GlobalId node_id = 0;
+  Dof dof = Dof::UX;
+};
+
+class DofManager {
+ public:
+  explicit DofManager(const Domain& domain) {
+    for (const auto& [node_id, node] : domain.nodes) {
+      (void)node;
+      node_ids_.push_back(node_id);
+      for (Dof dof : allDofs()) {
+        const LocalIndex full_index = static_cast<LocalIndex>(all_dofs_.size());
+        const auto key = std::make_pair(node_id, dofIndex(dof));
+        all_dofs_.push_back(key);
+        full_index_by_key_[key] = full_index;
+      }
+    }
+    for (const BoundaryCondition& boundary : domain.boundary_conditions) {
+      if (!validAbaqusDofRange(boundary.first_dof, boundary.last_dof)) {
+        continue;
+      }
+      for (GlobalId node_id : resolveNodeTarget(domain, boundary.target)) {
+        for (int dof = boundary.first_dof; dof <= boundary.last_dof; ++dof) {
+          constrained_.insert(std::make_pair(node_id, dof - 1));
+        }
+      }
+    }
+    for (const auto& key : all_dofs_) {
+      const LocalIndex full_index = full_index_by_key_.at(key);
+      if (constrained_.count(key) == 0) {
+        equation_by_key_[key] = static_cast<EquationId>(free_full_indices_.size());
+        free_full_indices_.push_back(full_index);
+      } else {
+        equation_by_key_[key] = -1;
+        constrained_full_indices_.push_back(full_index);
+      }
+    }
+  }
+
+  LocalIndex fullDofCount() const {
+    return static_cast<LocalIndex>(all_dofs_.size());
+  }
+
+  LocalIndex freeDofCount() const {
+    return static_cast<LocalIndex>(free_full_indices_.size());
+  }
+
+  LocalIndex constrainedDofCount() const {
+    return static_cast<LocalIndex>(constrained_full_indices_.size());
+  }
+
+  const std::vector<GlobalId>& nodeIds() const {
+    return node_ids_;
+  }
+
+  const std::vector<LocalIndex>& freeFullIndices() const {
+    return free_full_indices_;
+  }
+
+  const std::vector<LocalIndex>& constrainedFullIndices() const {
+    return constrained_full_indices_;
+  }
+
+  DofAddress fullDof(LocalIndex full_index) const {
+    const auto& key = all_dofs_.at(static_cast<std::size_t>(full_index));
+    return {key.first, static_cast<Dof>(key.second)};
+  }
+
+  LocalIndex fullIndex(GlobalId node_id, Dof dof) const {
+    return full_index_by_key_.at(std::make_pair(node_id, dofIndex(dof)));
+  }
+
+  EquationId equation(GlobalId node_id, Dof dof) const {
+    return equation_by_key_.at(std::make_pair(node_id, dofIndex(dof)));
+  }
+
+  bool isConstrained(GlobalId node_id, Dof dof) const {
+    return constrained_.count(std::make_pair(node_id, dofIndex(dof))) != 0;
+  }
+
+  std::vector<Real> reduceFullVector(const std::vector<Real>& full) const {
+    std::vector<Real> reduced;
+    reduced.reserve(free_full_indices_.size());
+    for (LocalIndex full_index : free_full_indices_) {
+      reduced.push_back(full.at(static_cast<std::size_t>(full_index)));
+    }
+    return reduced;
+  }
+
+  std::vector<Real> reconstructFullVector(const std::vector<Real>& reduced) const {
+    std::vector<Real> full(static_cast<std::size_t>(fullDofCount()), 0.0);
+    for (std::size_t i = 0; i < free_full_indices_.size(); ++i) {
+      full[static_cast<std::size_t>(free_full_indices_[i])] = reduced.at(i);
+    }
+    return full;
+  }
+
+  std::array<LocalIndex, 24> elementFullDofIndices(const Element& element) const {
+    std::array<LocalIndex, 24> indices{};
+    for (LocalIndex node = 0; node < 4; ++node) {
+      for (Dof dof : allDofs()) {
+        const LocalIndex local = 6 * node + dofIndex(dof);
+        indices[static_cast<std::size_t>(local)] = fullIndex(element.node_ids[static_cast<std::size_t>(node)], dof);
+      }
+    }
+    return indices;
+  }
+
+  std::array<EquationId, 24> elementEquationIds(const Element& element) const {
+    std::array<EquationId, 24> equations{};
+    for (LocalIndex node = 0; node < 4; ++node) {
+      for (Dof dof : allDofs()) {
+        const LocalIndex local = 6 * node + dofIndex(dof);
+        equations[static_cast<std::size_t>(local)] = equation(element.node_ids[static_cast<std::size_t>(node)], dof);
+      }
+    }
+    return equations;
+  }
+
+ private:
+  std::vector<GlobalId> node_ids_;
+  std::vector<std::pair<GlobalId, int>> all_dofs_;
+  std::set<std::pair<GlobalId, int>> constrained_;
+  std::map<std::pair<GlobalId, int>, LocalIndex> full_index_by_key_;
+  std::map<std::pair<GlobalId, int>, EquationId> equation_by_key_;
+  std::vector<LocalIndex> free_full_indices_;
+  std::vector<LocalIndex> constrained_full_indices_;
+};
+
+}  // namespace fesa

include/fesa/Core/Domain.hpp

@@ -0,0 +1,81 @@
+#pragma once
+
+#include "fesa/Boundary/Boundary.hpp"
+#include "fesa/Core/Types.hpp"
+#include "fesa/Load/Load.hpp"
+#include "fesa/Property/Property.hpp"
+#include "fesa/Util/String.hpp"
+
+#include <array>
+#include <map>
+#include <string>
+#include <vector>
+
+namespace fesa {
+
+struct Vec3 {
+  Real x = 0.0;
+  Real y = 0.0;
+  Real z = 0.0;
+};
+
+struct Node {
+  GlobalId id = 0;
+  Vec3 coordinates;
+};
+
+enum class ElementType { MITC4 };
+
+inline std::string elementTypeLabel(ElementType type) {
+  switch (type) {
+    case ElementType::MITC4:
+      return "MITC4";
+  }
+  return "UNKNOWN";
+}
+
+struct Element {
+  GlobalId id = 0;
+  ElementType type = ElementType::MITC4;
+  std::array<GlobalId, 4> node_ids{};
+  std::string source_elset;
+};
+
+struct NodeSet {
+  std::string name;
+  std::vector<GlobalId> node_ids;
+};
+
+struct ElementSet {
+  std::string name;
+  std::vector<GlobalId> element_ids;
+};
+
+struct Material {
+  std::string name;
+  Real elastic_modulus = 0.0;
+  Real poisson_ratio = 0.0;
+};
+
+struct StepDefinition {
+  std::string name = "Step-1";
+  std::string analysis_type = "linear_static";
+};
+
+struct Domain {
+  std::map<GlobalId, Node> nodes;
+  std::map<GlobalId, Element> elements;
+  std::map<std::string, NodeSet> node_sets;
+  std::map<std::string, ElementSet> element_sets;
+  std::map<std::string, Material> materials;
+  std::vector<ShellSection> shell_sections;
+  std::vector<BoundaryCondition> boundary_conditions;
+  std::vector<NodalLoad> loads;
+  std::vector<StepDefinition> steps;
+
+  static std::string key(const std::string& label) {
+    return lower(trim(label));
+  }
+};
+
+}  // namespace fesa

include/fesa/Core/Types.hpp

@@ -0,0 +1,13 @@
+#pragma once
+
+#include <cstdint>
+
+namespace fesa {
+
+using Real = double;
+using GlobalId = std::int64_t;
+using LocalIndex = std::int64_t;
+using EquationId = std::int64_t;
+using SparseIndex = std::int64_t;
+
+}  // namespace fesa

include/fesa/Core/Validation.hpp

@@ -0,0 +1,219 @@
+#pragma once
+
+#include "fesa/Core/Dof.hpp"
+#include "fesa/Core/Domain.hpp"
+#include "fesa/Util/Diagnostics.hpp"
+#include "fesa/Util/String.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <optional>
+#include <set>
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace fesa {
+
+inline std::optional<GlobalId> numericTarget(const std::string& target) {
+  return parseInt64(target);
+}
+
+inline std::vector<GlobalId> resolveNodeTarget(const Domain& domain, const std::string& target, std::vector<Diagnostic>* diagnostics = nullptr,
+                                               const std::string& diagnostic_keyword = "node target") {
+  if (auto node_id = numericTarget(target)) {
+    if (domain.nodes.count(*node_id) == 0) {
+      if (diagnostics != nullptr) {
+        diagnostics->push_back(
+            makeDiagnostic(Severity::Error, "FESA-VALIDATION-MISSING-NODE", "Missing node target: " + target, diagnostic_keyword));
+      }
+      return {};
+    }
+    return {*node_id};
+  }
+  auto set_it = domain.node_sets.find(Domain::key(target));
+  if (set_it == domain.node_sets.end()) {
+    if (diagnostics != nullptr) {
+      diagnostics->push_back(
+          makeDiagnostic(Severity::Error, "FESA-VALIDATION-MISSING-NSET", "Missing node set: " + target, diagnostic_keyword));
+    }
+    return {};
+  }
+  return set_it->second.node_ids;
+}
+
+inline const ShellSection* shellSectionForElement(const Domain& domain, GlobalId element_id) {
+  for (const ShellSection& section : domain.shell_sections) {
+    auto set_it = domain.element_sets.find(Domain::key(section.element_set));
+    if (set_it == domain.element_sets.end()) {
+      continue;
+    }
+    if (std::find(set_it->second.element_ids.begin(), set_it->second.element_ids.end(), element_id) != set_it->second.element_ids.end()) {
+      return &section;
+    }
+  }
+  return nullptr;
+}
+
+inline std::string dofNameOrNumber(int abaqus_dof) {
+  auto dof = dofFromAbaqus(abaqus_dof);
+  if (dof) {
+    return dofLabel(*dof);
+  }
+  return "DOF " + std::to_string(abaqus_dof);
+}
+
+inline bool validAbaqusDofRange(int first, int last) {
+  return dofFromAbaqus(first).has_value() && dofFromAbaqus(last).has_value() && first <= last;
+}
+
+inline std::vector<Diagnostic> validateDomain(const Domain& domain) {
+  std::vector<Diagnostic> diagnostics;
+  if (domain.elements.empty()) {
+    diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SINGULAR-NO-ACTIVE-ELEMENTS",
+                                         "No active elements exist in the current model", "analysis model"));
+  }
+  if (domain.boundary_conditions.empty()) {
+    diagnostics.push_back(makeDiagnostic(Severity::Warning, "FESA-SINGULAR-NO-BOUNDARY", "No boundary constraints are defined", "boundary"));
+  }
+  for (const auto& [set_key, set] : domain.node_sets) {
+    (void)set_key;
+    for (GlobalId node_id : set.node_ids) {
+      if (domain.nodes.count(node_id) == 0) {
+        diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-NSET-MISSING-NODE",
+                                             "Node set " + set.name + " references missing node " + std::to_string(node_id),
+                                             "nset"));
+      }
+    }
+  }
+  for (const auto& [set_key, set] : domain.element_sets) {
+    (void)set_key;
+    for (GlobalId element_id : set.element_ids) {
+      if (domain.elements.count(element_id) == 0) {
+        diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-ELSET-MISSING-ELEMENT",
+                                             "Element set " + set.name + " references missing element " + std::to_string(element_id),
+                                             "elset"));
+      }
+    }
+  }
+  for (const auto& [id, element] : domain.elements) {
+    for (GlobalId node_id : element.node_ids) {
+      if (domain.nodes.count(node_id) == 0) {
+        diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-ELEMENT-MISSING-NODE",
+                                             "Element " + std::to_string(id) + " references missing node " + std::to_string(node_id),
+                                             "element"));
+      }
+    }
+    const ShellSection* section = shellSectionForElement(domain, id);
+    if (section == nullptr) {
+      diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-MISSING-PROPERTY",
+                                           "Element " + std::to_string(id) + " has no assigned shell section", "element"));
+    }
+  }
+  for (const ShellSection& section : domain.shell_sections) {
+    if (section.thickness <= 0.0 || !std::isfinite(section.thickness)) {
+      diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-NONPOSITIVE-THICKNESS",
+                                           "Shell section for element set " + section.element_set + " has non-positive thickness",
+                                           "shell section"));
+    }
+    if (domain.element_sets.count(Domain::key(section.element_set)) == 0) {
+      diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-MISSING-ELSET",
+                                           "Shell section references missing element set: " + section.element_set, "shell section"));
+    }
+    auto material_it = domain.materials.find(Domain::key(section.material));
+    if (material_it == domain.materials.end()) {
+      diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-MISSING-MATERIAL",
+                                           "Shell section references missing material: " + section.material, "shell section"));
+    } else if (material_it->second.elastic_modulus <= 0.0) {
+      diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-INCOMPLETE-MATERIAL",
+                                           "Material has no valid elastic constants: " + section.material, "material"));
+    }
+  }
+  for (const BoundaryCondition& boundary : domain.boundary_conditions) {
+    if (!validAbaqusDofRange(boundary.first_dof, boundary.last_dof)) {
+      diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-BOUNDARY-DOF",
+                                           "Boundary target " + boundary.target + " has invalid DOF range " +
+                                               dofNameOrNumber(boundary.first_dof) + " to " + dofNameOrNumber(boundary.last_dof),
+                                           "boundary"));
+    }
+    (void)resolveNodeTarget(domain, boundary.target, &diagnostics, "boundary");
+  }
+  for (const NodalLoad& load : domain.loads) {
+    if (!dofFromAbaqus(load.dof)) {
+      diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-CLOAD-DOF",
+                                           "Load target " + load.target + " has invalid " + dofNameOrNumber(load.dof), "cload"));
+    }
+    (void)resolveNodeTarget(domain, load.target, &diagnostics, "cload");
+  }
+  const bool any_nonzero_load = std::any_of(domain.loads.begin(), domain.loads.end(), [](const NodalLoad& load) {
+    return std::fabs(load.magnitude) > 0.0;
+  });
+  if (!any_nonzero_load) {
+    diagnostics.push_back(makeDiagnostic(Severity::Warning, "FESA-SINGULAR-NO-NONZERO-LOAD", "No nonzero load is defined", "cload"));
+  }
+
+  std::set<std::pair<GlobalId, int>> constrained_dofs;
+  for (const BoundaryCondition& boundary : domain.boundary_conditions) {
+    if (!validAbaqusDofRange(boundary.first_dof, boundary.last_dof)) {
+      continue;
+    }
+    for (GlobalId node_id : resolveNodeTarget(domain, boundary.target)) {
+      if (domain.nodes.count(node_id) == 0) {
+        continue;
+      }
+      for (int dof = boundary.first_dof; dof <= boundary.last_dof; ++dof) {
+        constrained_dofs.insert(std::make_pair(node_id, dof - 1));
+      }
+    }
+  }
+
+  std::set<GlobalId> active_connectivity_nodes;
+  for (const auto& [element_id, element] : domain.elements) {
+    (void)element_id;
+    for (GlobalId node_id : element.node_ids) {
+      if (domain.nodes.count(node_id) != 0) {
+        active_connectivity_nodes.insert(node_id);
+      }
+    }
+  }
+
+  LocalIndex free_dof_count = 0;
+  LocalIndex weak_drilling_count = 0;
+  GlobalId weak_drilling_example = 0;
+  for (const auto& [node_id, node] : domain.nodes) {
+    (void)node;
+    for (Dof dof : allDofs()) {
+      const auto key = std::make_pair(node_id, dofIndex(dof));
+      if (constrained_dofs.count(key) != 0) {
+        continue;
+      }
+      ++free_dof_count;
+      if (!domain.elements.empty() && active_connectivity_nodes.count(node_id) == 0) {
+        diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SINGULAR-DOF-UNTOUCHED",
+                                             "Node " + std::to_string(node_id) + " DOF " + dofLabel(dof) +
+                                                 " is free but is not touched by active element connectivity",
+                                             "dof"));
+      }
+      if (!domain.elements.empty() && active_connectivity_nodes.count(node_id) != 0 && dof == Dof::RZ) {
+        if (weak_drilling_count == 0) {
+          weak_drilling_example = node_id;
+        }
+        ++weak_drilling_count;
+      }
+    }
+  }
+  if (!domain.nodes.empty() && free_dof_count == 0) {
+    diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SINGULAR-NO-FREE-DOFS",
+                                         "No free DOFs exist after applying boundary constraints", "dof"));
+  }
+  if (weak_drilling_count > 0) {
+    diagnostics.push_back(makeDiagnostic(Severity::Warning, "FESA-SINGULAR-WEAK-DRILLING-DOF",
+                                         "Node " + std::to_string(weak_drilling_example) +
+                                             " DOF RZ is free; drilling rotation is weakly stabilized in Phase 1 (" +
+                                             std::to_string(weak_drilling_count) + " free drilling DOF(s))",
+                                         "dof"));
+  }
+  return diagnostics;
+}
+
+}  // namespace fesa

include/fesa/Load/Load.hpp

@@ -1,9 +1,22 @@
 #pragma once
 
+#include "fesa/Core/Types.hpp"
 #include "fesa/ModuleInfo.hpp"
 
+#include <string>
+
 namespace fesa::module {
 
 inline constexpr std::string_view kLoad = "Load";
 
 }  // namespace fesa::module
+
+namespace fesa {
+
+struct NodalLoad {
+  std::string target;
+  int dof = 0;
+  Real magnitude = 0.0;
+};
+
+}  // namespace fesa

include/fesa/Property/Property.hpp

@@ -1,9 +1,22 @@
 #pragma once
 
+#include "fesa/Core/Types.hpp"
 #include "fesa/ModuleInfo.hpp"
 
+#include <string>
+
 namespace fesa::module {
 
 inline constexpr std::string_view kProperty = "Property";
 
 }  // namespace fesa::module
+
+namespace fesa {
+
+struct ShellSection {
+  std::string element_set;
+  std::string material;
+  Real thickness = 0.0;
+};
+
+}  // namespace fesa

include/fesa/Util/Diagnostics.hpp

@@ -0,0 +1,44 @@
+#pragma once
+
+#include "fesa/Core/Types.hpp"
+
+#include <algorithm>
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace fesa {
+
+enum class Severity { Info, Warning, Error };
+
+struct SourceLocation {
+  std::string file;
+  LocalIndex line = 0;
+  std::string keyword;
+};
+
+struct Diagnostic {
+  Severity severity = Severity::Error;
+  std::string code;
+  std::string message;
+  SourceLocation source;
+};
+
+inline bool hasError(const std::vector<Diagnostic>& diagnostics) {
+  return std::any_of(diagnostics.begin(), diagnostics.end(), [](const Diagnostic& diagnostic) {
+    return diagnostic.severity == Severity::Error;
+  });
+}
+
+inline bool containsDiagnostic(const std::vector<Diagnostic>& diagnostics, const std::string& code) {
+  return std::any_of(diagnostics.begin(), diagnostics.end(), [&](const Diagnostic& diagnostic) {
+    return diagnostic.code == code;
+  });
+}
+
+inline Diagnostic makeDiagnostic(Severity severity, std::string code, std::string message, std::string keyword,
+                                 std::string file = "<domain>", LocalIndex line = 0) {
+  return {severity, std::move(code), std::move(message), {std::move(file), line, std::move(keyword)}};
+}
+
+}  // namespace fesa

include/fesa/Util/String.hpp

@@ -0,0 +1,94 @@
+#pragma once
+
+#include "fesa/Core/Types.hpp"
+
+#include <algorithm>
+#include <cctype>
+#include <optional>
+#include <sstream>
+#include <string>
+#include <vector>
+
+namespace fesa {
+
+inline std::string trim(std::string text) {
+  auto is_space = [](unsigned char c) { return std::isspace(c) != 0; };
+  text.erase(text.begin(), std::find_if(text.begin(), text.end(), [&](unsigned char c) { return !is_space(c); }));
+  text.erase(std::find_if(text.rbegin(), text.rend(), [&](unsigned char c) { return !is_space(c); }).base(), text.end());
+  return text;
+}
+
+inline std::string lower(std::string text) {
+  std::transform(text.begin(), text.end(), text.begin(), [](unsigned char c) {
+    return static_cast<char>(std::tolower(c));
+  });
+  return text;
+}
+
+inline std::vector<std::string> splitCsv(const std::string& line) {
+  std::vector<std::string> fields;
+  std::string field;
+  std::istringstream stream(line);
+  while (std::getline(stream, field, ',')) {
+    fields.push_back(trim(field));
+  }
+  if (!line.empty() && line.back() == ',') {
+    fields.emplace_back();
+  }
+  return fields;
+}
+
+inline std::optional<Real> parseReal(std::string token) {
+  token = trim(token);
+  if (token.empty()) {
+    return std::nullopt;
+  }
+  std::replace(token.begin(), token.end(), 'D', 'E');
+  std::replace(token.begin(), token.end(), 'd', 'e');
+  try {
+    std::size_t used = 0;
+    Real value = std::stod(token, &used);
+    if (used != token.size()) {
+      return std::nullopt;
+    }
+    return value;
+  } catch (...) {
+    return std::nullopt;
+  }
+}
+
+inline std::optional<GlobalId> parseInt64(const std::string& token) {
+  std::string value_text = trim(token);
+  if (value_text.empty()) {
+    return std::nullopt;
+  }
+  try {
+    std::size_t used = 0;
+    long long value = std::stoll(value_text, &used);
+    if (used != value_text.size()) {
+      return std::nullopt;
+    }
+    return static_cast<GlobalId>(value);
+  } catch (...) {
+    return std::nullopt;
+  }
+}
+
+inline void addUnique(std::vector<GlobalId>& values, GlobalId value) {
+  if (std::find(values.begin(), values.end(), value) == values.end()) {
+    values.push_back(value);
+  }
+}
+
+inline std::vector<GlobalId> generatedRange(GlobalId first, GlobalId last, GlobalId increment) {
+  std::vector<GlobalId> values;
+  if (increment <= 0) {
+    return values;
+  }
+  for (GlobalId value = first; value <= last; value += increment) {
+    values.push_back(value);
+  }
+  return values;
+}
+
+}  // namespace fesa

include/fesa/Util/Util.hpp

@@ -1,5 +1,7 @@
 #pragma once
 
+#include "fesa/Util/Diagnostics.hpp"
+#include "fesa/Util/String.hpp"
 #include "fesa/ModuleInfo.hpp"
 
 namespace fesa::module {

include/fesa/fesa.hpp

@@ -1,6 +1,11 @@
 #pragma once
 
+#include "fesa/Boundary/Boundary.hpp"
+#include "fesa/Core/Core.hpp"
+#include "fesa/Load/Load.hpp"
 #include "fesa/ModuleInfo.hpp"
+#include "fesa/Property/Property.hpp"
+#include "fesa/Util/Util.hpp"
 
 #include <algorithm>
 #include <array>
@@ -22,160 +27,6 @@
 
 namespace fesa {
 
-using Real = double;
-using GlobalId = std::int64_t;
-using LocalIndex = std::int64_t;
-using EquationId = std::int64_t;
-using SparseIndex = std::int64_t;
-
-enum class Severity { Info, Warning, Error };
-
-struct SourceLocation {
-  std::string file;
-  LocalIndex line = 0;
-  std::string keyword;
-};
-
-struct Diagnostic {
-  Severity severity = Severity::Error;
-  std::string code;
-  std::string message;
-  SourceLocation source;
-};
-
-inline bool hasError(const std::vector<Diagnostic>& diagnostics) {
-  return std::any_of(diagnostics.begin(), diagnostics.end(), [](const Diagnostic& diagnostic) {
-    return diagnostic.severity == Severity::Error;
-  });
-}
-
-inline bool containsDiagnostic(const std::vector<Diagnostic>& diagnostics, const std::string& code) {
-  return std::any_of(diagnostics.begin(), diagnostics.end(), [&](const Diagnostic& diagnostic) {
-    return diagnostic.code == code;
-  });
-}
-
-inline Diagnostic makeDiagnostic(Severity severity, std::string code, std::string message, std::string keyword,
-                                 std::string file = "<domain>", LocalIndex line = 0) {
-  return {severity, std::move(code), std::move(message), {std::move(file), line, std::move(keyword)}};
-}
-
-inline std::string trim(std::string text) {
-  auto is_space = [](unsigned char c) { return std::isspace(c) != 0; };
-  text.erase(text.begin(), std::find_if(text.begin(), text.end(), [&](unsigned char c) { return !is_space(c); }));
-  text.erase(std::find_if(text.rbegin(), text.rend(), [&](unsigned char c) { return !is_space(c); }).base(), text.end());
-  return text;
-}
-
-inline std::string lower(std::string text) {
-  std::transform(text.begin(), text.end(), text.begin(), [](unsigned char c) {
-    return static_cast<char>(std::tolower(c));
-  });
-  return text;
-}
-
-inline std::vector<std::string> splitCsv(const std::string& line) {
-  std::vector<std::string> fields;
-  std::string field;
-  std::istringstream stream(line);
-  while (std::getline(stream, field, ',')) {
-    fields.push_back(trim(field));
-  }
-  if (!line.empty() && line.back() == ',') {
-    fields.emplace_back();
-  }
-  return fields;
-}
-
-inline std::optional<Real> parseReal(std::string token) {
-  token = trim(token);
-  if (token.empty()) {
-    return std::nullopt;
-  }
-  std::replace(token.begin(), token.end(), 'D', 'E');
-  std::replace(token.begin(), token.end(), 'd', 'e');
-  try {
-    std::size_t used = 0;
-    Real value = std::stod(token, &used);
-    if (used != token.size()) {
-      return std::nullopt;
-    }
-    return value;
-  } catch (...) {
-    return std::nullopt;
-  }
-}
-
-inline std::optional<GlobalId> parseInt64(const std::string& token) {
-  std::string value_text = trim(token);
-  if (value_text.empty()) {
-    return std::nullopt;
-  }
-  try {
-    std::size_t used = 0;
-    long long value = std::stoll(value_text, &used);
-    if (used != value_text.size()) {
-      return std::nullopt;
-    }
-    return static_cast<GlobalId>(value);
-  } catch (...) {
-    return std::nullopt;
-  }
-}
-
-enum class Dof : int { UX = 0, UY = 1, UZ = 2, RX = 3, RY = 4, RZ = 5 };
-
-inline std::array<Dof, 6> allDofs() {
-  return {Dof::UX, Dof::UY, Dof::UZ, Dof::RX, Dof::RY, Dof::RZ};
-}
-
-inline int dofIndex(Dof dof) {
-  return static_cast<int>(dof);
-}
-
-inline int abaqusDofNumber(Dof dof) {
-  return dofIndex(dof) + 1;
-}
-
-inline std::optional<Dof> dofFromAbaqus(int dof) {
-  if (dof < 1 || dof > 6) {
-    return std::nullopt;
-  }
-  return static_cast<Dof>(dof - 1);
-}
-
-inline const char* dofLabel(Dof dof) {
-  switch (dof) {
-    case Dof::UX:
-      return "UX";
-    case Dof::UY:
-      return "UY";
-    case Dof::UZ:
-      return "UZ";
-    case Dof::RX:
-      return "RX";
-    case Dof::RY:
-      return "RY";
-    case Dof::RZ:
-      return "RZ";
-  }
-  return "";
-}
-
-inline std::vector<std::string> displacementComponentLabels() {
-  return {"UX", "UY", "UZ", "RX", "RY", "RZ"};
-}
-
-inline std::vector<std::string> reactionComponentLabels() {
-  return {"RFX", "RFY", "RFZ", "RMX", "RMY", "RMZ"};
-}
-
-struct Vec3 {
-  Real x = 0.0;
-  Real y = 0.0;
-  Real z = 0.0;
-};
-
 inline Vec3 operator+(const Vec3& a, const Vec3& b) {
   return {a.x + b.x, a.y + b.y, a.z + b.z};
 }
@@ -224,101 +75,6 @@ inline Vec3 normalized(const Vec3& value) {
   return (1.0 / length) * value;
 }
 
-struct Node {
-  GlobalId id = 0;
-  Vec3 coordinates;
-};
-
-enum class ElementType { MITC4 };
-
-inline std::string elementTypeLabel(ElementType type) {
-  switch (type) {
-    case ElementType::MITC4:
-      return "MITC4";
-  }
-  return "UNKNOWN";
-}
-
-struct Element {
-  GlobalId id = 0;
-  ElementType type = ElementType::MITC4;
-  std::array<GlobalId, 4> node_ids{};
-  std::string source_elset;
-};
-
-struct NodeSet {
-  std::string name;
-  std::vector<GlobalId> node_ids;
-};
-
-struct ElementSet {
-  std::string name;
-  std::vector<GlobalId> element_ids;
-};
-
-struct Material {
-  std::string name;
-  Real elastic_modulus = 0.0;
-  Real poisson_ratio = 0.0;
-};
-
-struct ShellSection {
-  std::string element_set;
-  std::string material;
-  Real thickness = 0.0;
-};
-
-struct BoundaryCondition {
-  std::string target;
-  int first_dof = 0;
-  int last_dof = 0;
-  Real magnitude = 0.0;
-};
-
-struct NodalLoad {
-  std::string target;
-  int dof = 0;
-  Real magnitude = 0.0;
-};
-
-struct StepDefinition {
-  std::string name = "Step-1";
-  std::string analysis_type = "linear_static";
-};
-
-struct Domain {
-  std::map<GlobalId, Node> nodes;
-  std::map<GlobalId, Element> elements;
-  std::map<std::string, NodeSet> node_sets;
-  std::map<std::string, ElementSet> element_sets;
-  std::map<std::string, Material> materials;
-  std::vector<ShellSection> shell_sections;
-  std::vector<BoundaryCondition> boundary_conditions;
-  std::vector<NodalLoad> loads;
-  std::vector<StepDefinition> steps;
-
-  static std::string key(const std::string& label) {
-    return lower(trim(label));
-  }
-};
-
-inline void addUnique(std::vector<GlobalId>& values, GlobalId value) {
-  if (std::find(values.begin(), values.end(), value) == values.end()) {
-    values.push_back(value);
-  }
-}
-
-inline std::vector<GlobalId> generatedRange(GlobalId first, GlobalId last, GlobalId increment) {
-  std::vector<GlobalId> values;
-  if (increment <= 0) {
-    return values;
-  }
-  for (GlobalId value = first; value <= last; value += increment) {
-    values.push_back(value);
-  }
-  return values;
-}
-
 struct KeywordLine {
   std::string name;
   std::map<std::string, std::string> parameters;
@@ -762,392 +518,6 @@ class AbaqusInputParser {
   }
 };
 
-inline std::optional<GlobalId> numericTarget(const std::string& target) {
-  return parseInt64(target);
-}
-
-inline std::vector<GlobalId> resolveNodeTarget(const Domain& domain, const std::string& target, std::vector<Diagnostic>* diagnostics = nullptr,
-                                               const std::string& diagnostic_keyword = "node target") {
-  if (auto node_id = numericTarget(target)) {
-    if (domain.nodes.count(*node_id) == 0) {
-      if (diagnostics != nullptr) {
-        diagnostics->push_back(
-            makeDiagnostic(Severity::Error, "FESA-VALIDATION-MISSING-NODE", "Missing node target: " + target, diagnostic_keyword));
-      }
-      return {};
-    }
-    return {*node_id};
-  }
-  auto set_it = domain.node_sets.find(Domain::key(target));
-  if (set_it == domain.node_sets.end()) {
-    if (diagnostics != nullptr) {
-      diagnostics->push_back(
-          makeDiagnostic(Severity::Error, "FESA-VALIDATION-MISSING-NSET", "Missing node set: " + target, diagnostic_keyword));
-    }
-    return {};
-  }
-  return set_it->second.node_ids;
-}
-
-inline const ShellSection* shellSectionForElement(const Domain& domain, GlobalId element_id) {
-  for (const ShellSection& section : domain.shell_sections) {
-    auto set_it = domain.element_sets.find(Domain::key(section.element_set));
-    if (set_it == domain.element_sets.end()) {
-      continue;
-    }
-    if (std::find(set_it->second.element_ids.begin(), set_it->second.element_ids.end(), element_id) != set_it->second.element_ids.end()) {
-      return &section;
-    }
-  }
-  return nullptr;
-}
-
-inline std::string dofNameOrNumber(int abaqus_dof) {
-  auto dof = dofFromAbaqus(abaqus_dof);
-  if (dof) {
-    return dofLabel(*dof);
-  }
-  return "DOF " + std::to_string(abaqus_dof);
-}
-
-inline bool validAbaqusDofRange(int first, int last) {
-  return dofFromAbaqus(first).has_value() && dofFromAbaqus(last).has_value() && first <= last;
-}
-
-inline std::vector<Diagnostic> validateDomain(const Domain& domain) {
-  std::vector<Diagnostic> diagnostics;
-  if (domain.elements.empty()) {
-    diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SINGULAR-NO-ACTIVE-ELEMENTS",
-                                         "No active elements exist in the current model", "analysis model"));
-  }
-  if (domain.boundary_conditions.empty()) {
-    diagnostics.push_back(makeDiagnostic(Severity::Warning, "FESA-SINGULAR-NO-BOUNDARY", "No boundary constraints are defined", "boundary"));
-  }
-  for (const auto& [set_key, set] : domain.node_sets) {
-    (void)set_key;
-    for (GlobalId node_id : set.node_ids) {
-      if (domain.nodes.count(node_id) == 0) {
-        diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-NSET-MISSING-NODE",
-                                             "Node set " + set.name + " references missing node " + std::to_string(node_id),
-                                             "nset"));
-      }
-    }
-  }
-  for (const auto& [set_key, set] : domain.element_sets) {
-    (void)set_key;
-    for (GlobalId element_id : set.element_ids) {
-      if (domain.elements.count(element_id) == 0) {
-        diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-ELSET-MISSING-ELEMENT",
-                                             "Element set " + set.name + " references missing element " + std::to_string(element_id),
-                                             "elset"));
-      }
-    }
-  }
-  for (const auto& [id, element] : domain.elements) {
-    for (GlobalId node_id : element.node_ids) {
-      if (domain.nodes.count(node_id) == 0) {
-        diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-ELEMENT-MISSING-NODE",
-                                             "Element " + std::to_string(id) + " references missing node " + std::to_string(node_id),
-                                             "element"));
-      }
-    }
-    const ShellSection* section = shellSectionForElement(domain, id);
-    if (section == nullptr) {
-      diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-MISSING-PROPERTY",
-                                           "Element " + std::to_string(id) + " has no assigned shell section", "element"));
-    }
-  }
-  for (const ShellSection& section : domain.shell_sections) {
-    if (section.thickness <= 0.0 || !std::isfinite(section.thickness)) {
-      diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-NONPOSITIVE-THICKNESS",
-                                           "Shell section for element set " + section.element_set + " has non-positive thickness",
-                                           "shell section"));
-    }
-    if (domain.element_sets.count(Domain::key(section.element_set)) == 0) {
-      diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-MISSING-ELSET",
-                                           "Shell section references missing element set: " + section.element_set, "shell section"));
-    }
-    auto material_it = domain.materials.find(Domain::key(section.material));
-    if (material_it == domain.materials.end()) {
-      diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-MISSING-MATERIAL",
-                                           "Shell section references missing material: " + section.material, "shell section"));
-    } else if (material_it->second.elastic_modulus <= 0.0) {
-      diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-INCOMPLETE-MATERIAL",
-                                           "Material has no valid elastic constants: " + section.material, "material"));
-    }
-  }
-  for (const BoundaryCondition& boundary : domain.boundary_conditions) {
-    if (!validAbaqusDofRange(boundary.first_dof, boundary.last_dof)) {
-      diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-BOUNDARY-DOF",
-                                           "Boundary target " + boundary.target + " has invalid DOF range " +
-                                               dofNameOrNumber(boundary.first_dof) + " to " + dofNameOrNumber(boundary.last_dof),
-                                           "boundary"));
-    }
-    (void)resolveNodeTarget(domain, boundary.target, &diagnostics, "boundary");
-  }
-  for (const NodalLoad& load : domain.loads) {
-    if (!dofFromAbaqus(load.dof)) {
-      diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-VALIDATION-CLOAD-DOF",
-                                           "Load target " + load.target + " has invalid " + dofNameOrNumber(load.dof), "cload"));
-    }
-    (void)resolveNodeTarget(domain, load.target, &diagnostics, "cload");
-  }
-  const bool any_nonzero_load = std::any_of(domain.loads.begin(), domain.loads.end(), [](const NodalLoad& load) {
-    return std::fabs(load.magnitude) > 0.0;
-  });
-  if (!any_nonzero_load) {
-    diagnostics.push_back(makeDiagnostic(Severity::Warning, "FESA-SINGULAR-NO-NONZERO-LOAD", "No nonzero load is defined", "cload"));
-  }
-
-  std::set<std::pair<GlobalId, int>> constrained_dofs;
-  for (const BoundaryCondition& boundary : domain.boundary_conditions) {
-    if (!validAbaqusDofRange(boundary.first_dof, boundary.last_dof)) {
-      continue;
-    }
-    for (GlobalId node_id : resolveNodeTarget(domain, boundary.target)) {
-      if (domain.nodes.count(node_id) == 0) {
-        continue;
-      }
-      for (int dof = boundary.first_dof; dof <= boundary.last_dof; ++dof) {
-        constrained_dofs.insert(std::make_pair(node_id, dof - 1));
-      }
-    }
-  }
-
-  std::set<GlobalId> active_connectivity_nodes;
-  for (const auto& [element_id, element] : domain.elements) {
-    (void)element_id;
-    for (GlobalId node_id : element.node_ids) {
-      if (domain.nodes.count(node_id) != 0) {
-        active_connectivity_nodes.insert(node_id);
-      }
-    }
-  }
-
-  LocalIndex free_dof_count = 0;
-  LocalIndex weak_drilling_count = 0;
-  GlobalId weak_drilling_example = 0;
-  for (const auto& [node_id, node] : domain.nodes) {
-    (void)node;
-    for (Dof dof : allDofs()) {
-      const auto key = std::make_pair(node_id, dofIndex(dof));
-      if (constrained_dofs.count(key) != 0) {
-        continue;
-      }
-      ++free_dof_count;
-      if (!domain.elements.empty() && active_connectivity_nodes.count(node_id) == 0) {
-        diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SINGULAR-DOF-UNTOUCHED",
-                                             "Node " + std::to_string(node_id) + " DOF " + dofLabel(dof) +
-                                                 " is free but is not touched by active element connectivity",
-                                             "dof"));
-      }
-      if (!domain.elements.empty() && active_connectivity_nodes.count(node_id) != 0 && dof == Dof::RZ) {
-        if (weak_drilling_count == 0) {
-          weak_drilling_example = node_id;
-        }
-        ++weak_drilling_count;
-      }
-    }
-  }
-  if (!domain.nodes.empty() && free_dof_count == 0) {
-    diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SINGULAR-NO-FREE-DOFS",
-                                         "No free DOFs exist after applying boundary constraints", "dof"));
-  }
-  if (weak_drilling_count > 0) {
-    diagnostics.push_back(makeDiagnostic(Severity::Warning, "FESA-SINGULAR-WEAK-DRILLING-DOF",
-                                         "Node " + std::to_string(weak_drilling_example) +
-                                             " DOF RZ is free; drilling rotation is weakly stabilized in Phase 1 (" +
-                                             std::to_string(weak_drilling_count) + " free drilling DOF(s))",
-                                         "dof"));
-  }
-  return diagnostics;
-}
-
-struct AnalysisModel {
-  StepDefinition step;
-  std::vector<GlobalId> active_element_ids;
-  std::vector<std::size_t> active_boundary_condition_indices;
-  std::vector<std::size_t> active_load_indices;
-  std::vector<std::size_t> active_shell_section_indices;
-  std::vector<std::string> active_material_keys;
-  std::vector<Diagnostic> diagnostics;
-
-  bool ok() const {
-    return !hasError(diagnostics);
-  }
-};
-
-inline AnalysisModel buildLinearStaticAnalysisModel(const Domain& domain, LocalIndex step_index = 0) {
-  AnalysisModel model;
-  if (domain.steps.empty()) {
-    model.step = {"Step-1", "linear_static"};
-  } else {
-    if (step_index < 0 || step_index >= static_cast<LocalIndex>(domain.steps.size())) {
-      model.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-ANALYSIS-STEP-INDEX",
-                                                 "Requested analysis step index is out of range", "analysis model"));
-      model.step = domain.steps.front();
-    } else {
-      model.step = domain.steps[static_cast<std::size_t>(step_index)];
-    }
-  }
-  if (domain.steps.size() > 1) {
-    model.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-ANALYSIS-MULTIPLE-STEPS",
-                                               "Phase 1 execution supports one active linear static step", "analysis model"));
-  }
-  if (model.step.analysis_type != "linear_static") {
-    model.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-ANALYSIS-UNSUPPORTED-STEP",
-                                               "Only linear static steps are supported in Phase 1", "analysis model"));
-  }
-  for (const auto& [element_id, element] : domain.elements) {
-    (void)element;
-    model.active_element_ids.push_back(element_id);
-  }
-  for (std::size_t i = 0; i < domain.boundary_conditions.size(); ++i) {
-    model.active_boundary_condition_indices.push_back(i);
-  }
-  for (std::size_t i = 0; i < domain.loads.size(); ++i) {
-    model.active_load_indices.push_back(i);
-  }
-  for (std::size_t i = 0; i < domain.shell_sections.size(); ++i) {
-    model.active_shell_section_indices.push_back(i);
-  }
-  for (const auto& [material_key, material] : domain.materials) {
-    (void)material;
-    model.active_material_keys.push_back(material_key);
-  }
-  return model;
-}
-
-struct DofAddress {
-  GlobalId node_id = 0;
-  Dof dof = Dof::UX;
-};
-
-class DofManager {
- public:
-  explicit DofManager(const Domain& domain) {
-    for (const auto& [node_id, node] : domain.nodes) {
-      (void)node;
-      node_ids_.push_back(node_id);
-      for (Dof dof : allDofs()) {
-        const LocalIndex full_index = static_cast<LocalIndex>(all_dofs_.size());
-        const auto key = std::make_pair(node_id, dofIndex(dof));
-        all_dofs_.push_back(key);
-        full_index_by_key_[key] = full_index;
-      }
-    }
-    for (const BoundaryCondition& boundary : domain.boundary_conditions) {
-      if (!validAbaqusDofRange(boundary.first_dof, boundary.last_dof)) {
-        continue;
-      }
-      for (GlobalId node_id : resolveNodeTarget(domain, boundary.target)) {
-        for (int dof = boundary.first_dof; dof <= boundary.last_dof; ++dof) {
-          constrained_.insert(std::make_pair(node_id, dof - 1));
-        }
-      }
-    }
-    for (const auto& key : all_dofs_) {
-      const LocalIndex full_index = full_index_by_key_.at(key);
-      if (constrained_.count(key) == 0) {
-        equation_by_key_[key] = static_cast<EquationId>(free_full_indices_.size());
-        free_full_indices_.push_back(full_index);
-      } else {
-        equation_by_key_[key] = -1;
-        constrained_full_indices_.push_back(full_index);
-      }
-    }
-  }
-
-  LocalIndex fullDofCount() const {
-    return static_cast<LocalIndex>(all_dofs_.size());
-  }
-
-  LocalIndex freeDofCount() const {
-    return static_cast<LocalIndex>(free_full_indices_.size());
-  }
-
-  LocalIndex constrainedDofCount() const {
-    return static_cast<LocalIndex>(constrained_full_indices_.size());
-  }
-
-  const std::vector<GlobalId>& nodeIds() const {
-    return node_ids_;
-  }
-
-  const std::vector<LocalIndex>& freeFullIndices() const {
-    return free_full_indices_;
-  }
-
-  const std::vector<LocalIndex>& constrainedFullIndices() const {
-    return constrained_full_indices_;
-  }
-
-  DofAddress fullDof(LocalIndex full_index) const {
-    const auto& key = all_dofs_.at(static_cast<std::size_t>(full_index));
-    return {key.first, static_cast<Dof>(key.second)};
-  }
-
-  LocalIndex fullIndex(GlobalId node_id, Dof dof) const {
-    return full_index_by_key_.at(std::make_pair(node_id, dofIndex(dof)));
-  }
-
-  EquationId equation(GlobalId node_id, Dof dof) const {
-    return equation_by_key_.at(std::make_pair(node_id, dofIndex(dof)));
-  }
-
-  bool isConstrained(GlobalId node_id, Dof dof) const {
-    return constrained_.count(std::make_pair(node_id, dofIndex(dof))) != 0;
-  }
-
-  std::vector<Real> reduceFullVector(const std::vector<Real>& full) const {
-    std::vector<Real> reduced;
-    reduced.reserve(free_full_indices_.size());
-    for (LocalIndex full_index : free_full_indices_) {
-      reduced.push_back(full.at(static_cast<std::size_t>(full_index)));
-    }
-    return reduced;
-  }
-
-  std::vector<Real> reconstructFullVector(const std::vector<Real>& reduced) const {
-    std::vector<Real> full(static_cast<std::size_t>(fullDofCount()), 0.0);
-    for (std::size_t i = 0; i < free_full_indices_.size(); ++i) {
-      full[static_cast<std::size_t>(free_full_indices_[i])] = reduced.at(i);
-    }
-    return full;
-  }
-
-  std::array<LocalIndex, 24> elementFullDofIndices(const Element& element) const {
-    std::array<LocalIndex, 24> indices{};
-    for (LocalIndex node = 0; node < 4; ++node) {
-      for (Dof dof : allDofs()) {
-        const LocalIndex local = 6 * node + dofIndex(dof);
-        indices[static_cast<std::size_t>(local)] = fullIndex(element.node_ids[static_cast<std::size_t>(node)], dof);
-      }
-    }
-    return indices;
-  }
-
-  std::array<EquationId, 24> elementEquationIds(const Element& element) const {
-    std::array<EquationId, 24> equations{};
-    for (LocalIndex node = 0; node < 4; ++node) {
-      for (Dof dof : allDofs()) {
-        const LocalIndex local = 6 * node + dofIndex(dof);
-        equations[static_cast<std::size_t>(local)] = equation(element.node_ids[static_cast<std::size_t>(node)], dof);
-      }
-    }
-    return equations;
-  }
-
- private:
-  std::vector<GlobalId> node_ids_;
-  std::vector<std::pair<GlobalId, int>> all_dofs_;
-  std::set<std::pair<GlobalId, int>> constrained_;
-  std::map<std::pair<GlobalId, int>, LocalIndex> full_index_by_key_;
-  std::map<std::pair<GlobalId, int>, EquationId> equation_by_key_;
-  std::vector<LocalIndex> free_full_indices_;
-  std::vector<LocalIndex> constrained_full_indices_;
-};
-
 struct SparsePatternEntry {
   EquationId row = 0;
   EquationId col = 0;
@@ -2419,14 +1789,6 @@ class InMemoryResultsWriter {
   ResultFile result_;
 };
 
-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;
-};
-
 struct AnalysisResult {
   AnalysisModel model;
   AnalysisState state;