feat(nonlinear-static-newton-raphson): step 10 - prepared-analysis-flow-and-full-procedure

include/fesa/analysis/linear_static_analysis.h

@@ -1,18 +1,7 @@
 #ifndef FESA_ANALYSIS_LINEAR_STATIC_ANALYSIS_H_
 #define FESA_ANALYSIS_LINEAR_STATIC_ANALYSIS_H_
 
-#include <memory>
-#include <vector>
-
 #include "fesa/analysis/analysis.h"
-#include "fesa/analysis/analysis_model.h"
-#include "fesa/analysis/analysis_state.h"
-#include "fesa/constraints/essential_constraint_policy.h"
-#include "fesa/elements/element.h"
-#include "fesa/fem/dof_manager.h"
-#include "fesa/math/sparse_matrix.h"
-#include "fesa/math/vector.h"
-#include "fesa/model/domain.h"
 
 namespace fesa {
 
@@ -20,9 +9,8 @@ class LinearSolver;
 class ParallelFor;
 class ResultsWriter;
 
-/// @brief Orchestrates the single-step linear-static procedure.
-/// @note Factorization, substitution, recovery, and writing remain separately
-/// observable through injected backend boundaries.
+/// @brief Retains the legacy Analysis facade for linear-static callers.
+/// @note Numerical stages are owned by LinearStaticProcedureExecutor.
 class LinearStaticAnalysis final : public Analysis {
  public:
   /// @brief Creates a procedure using non-owned backend adapters.
@@ -31,42 +19,14 @@ class LinearStaticAnalysis final : public Analysis {
                        LinearSolver& linear_solver,
                        ResultsWriter& results_writer);
 
-  /// @brief Runs the approved eight-stage linear-static lifecycle.
+  /// @brief Prepares, executes, and writes a linear-static analysis.
   /// @return The first stage failure or successful result finalization.
   Status Run(const AnalysisRequest& request) override;
 
  private:
-  /// @brief Initializes owned input and Domain state for a run candidate.
-  Status InitializeCandidate(const AnalysisRequest& request);
-  /// @brief Builds the non-owning active-model view.
-  Status BuildAnalysisModel();
-  /// @brief Creates runtime elements, stable DOFs, and the sparse pattern.
-  Status BuildDofMapAndSparsePattern();
-  /// @brief Assembles full stiffness and stable constraint partitions.
-  Status AssembleAndPartitionStiffness();
-  /// @brief Factorizes Kff before any load assembly.
-  Status FactorizeFreeSystem();
-  /// @brief Assembles loads and forms Ff-Kfc*dc without solving.
-  Status AssembleLoadsAndEffectiveRhs();
-  /// @brief Substitutes the retained factorization and reconstructs full d.
-  Status SubstituteAndReconstruct();
-  /// @brief Recovers a complete candidate before writing authoritative output.
-  Status RecoverAndWrite();
-
   const ParallelFor& parallel_for_;
   LinearSolver& linear_solver_;
   ResultsWriter& results_writer_;
-  AnalysisRequest request_;
-  std::unique_ptr<Domain> domain_;
-  std::unique_ptr<AnalysisModel> model_;
-  std::vector<std::unique_ptr<Element>> elements_;
-  ElementView element_view_;
-  std::unique_ptr<DofManager> dofs_;
-  std::unique_ptr<AnalysisState> state_;
-  std::unique_ptr<SparseMatrix> full_stiffness_;
-  std::unique_ptr<PartitionedStiffness> partitioned_stiffness_;
-  std::unique_ptr<Vector> effective_rhs_;
-  std::vector<Diagnostic> diagnostics_;
 };
 
 }  // namespace fesa

include/fesa/analysis/linear_static_procedure.h

@@ -0,0 +1,31 @@
+#ifndef FESA_ANALYSIS_LINEAR_STATIC_PROCEDURE_H_
+#define FESA_ANALYSIS_LINEAR_STATIC_PROCEDURE_H_
+
+#include "fesa/analysis/prepared_analysis.h"
+
+namespace fesa {
+
+class LinearSolver;
+class ParallelFor;
+
+/// @brief Executes the approved linear-static procedure on a prepared Domain.
+class LinearStaticProcedureExecutor {
+ public:
+  /// @brief Creates a procedure executor with caller-owned backend adapters.
+  /// @param parallel_for Deterministic element-local execution adapter.
+  /// @param linear_solver Factorization and substitution adapter.
+  LinearStaticProcedureExecutor(const ParallelFor& parallel_for,
+                                LinearSolver& linear_solver);
+
+  /// @brief Computes linear-static state without writing HDF5 output.
+  /// @return A validated state and sorted diagnostics, or a procedure failure.
+  Result<LinearStaticExecutionResult> Execute(const Domain& domain) const;
+
+ private:
+  const ParallelFor& parallel_for_;
+  LinearSolver& linear_solver_;
+};
+
+}  // namespace fesa
+
+#endif  // FESA_ANALYSIS_LINEAR_STATIC_PROCEDURE_H_

include/fesa/analysis/nonlinear_static_procedure.h

@@ -0,0 +1,33 @@
+#ifndef FESA_ANALYSIS_NONLINEAR_STATIC_PROCEDURE_H_
+#define FESA_ANALYSIS_NONLINEAR_STATIC_PROCEDURE_H_
+
+#include "fesa/analysis/prepared_analysis.h"
+
+namespace fesa {
+
+class NewtonIterationLogWriter;
+class ParallelFor;
+class TangentSolver;
+
+/// @brief Executes fixed-load-factor full-Newton nonlinear static analysis.
+class NonlinearStaticProcedureExecutor {
+ public:
+  /// @brief Creates a procedure executor with caller-owned backend adapters.
+  /// @param parallel_for Deterministic element-local execution adapter.
+  /// @param tangent_solver Fresh unsymmetric tangent solve adapter.
+  NonlinearStaticProcedureExecutor(const ParallelFor& parallel_for,
+                                   TangentSolver& tangent_solver);
+
+  /// @brief Computes accepted nonlinear frames and appends iteration records.
+  /// @return A validated nonlinear state without writing HDF5 output.
+  Result<NonlinearStaticExecutionResult> Execute(
+      const Domain& domain, NewtonIterationLogWriter& log) const;
+
+ private:
+  const ParallelFor& parallel_for_;
+  TangentSolver& tangent_solver_;
+};
+
+}  // namespace fesa
+
+#endif  // FESA_ANALYSIS_NONLINEAR_STATIC_PROCEDURE_H_

include/fesa/analysis/prepared_analysis.h

@@ -0,0 +1,71 @@
+#ifndef FESA_ANALYSIS_PREPARED_ANALYSIS_H_
+#define FESA_ANALYSIS_PREPARED_ANALYSIS_H_
+
+#include <filesystem>
+#include <variant>
+#include <vector>
+
+#include "fesa/analysis/analysis_state.h"
+#include "fesa/analysis/nonlinear_analysis_state.h"
+#include "fesa/core/diagnostic.h"
+#include "fesa/core/status.h"
+#include "fesa/model/domain.h"
+
+namespace fesa {
+
+class LinearStaticProcedureExecutor;
+class NewtonIterationLogWriter;
+class NonlinearResultsWriter;
+class NonlinearStaticProcedureExecutor;
+class ResultsWriter;
+struct DerivedOutputPaths;
+
+/// @brief Owns one parsed and semantically mapped analysis input.
+struct PreparedAnalysisInput {
+  Domain domain;
+  std::vector<Diagnostic> diagnostics;
+};
+
+/// @brief Parses and maps an input file exactly once for procedure execution.
+/// @param input_path Source Abaqus input path.
+/// @return A movable Domain and sorted mapping diagnostics, or an input
+/// failure.
+Result<PreparedAnalysisInput> PrepareAnalysisInput(
+    const std::filesystem::path& input_path);
+
+/// @brief Stores a completed linear-static calculation before result writing.
+struct LinearStaticExecutionResult {
+  AnalysisState state;
+  std::vector<Diagnostic> diagnostics;
+};
+
+/// @brief Stores a completed nonlinear-static calculation before result
+/// writing.
+struct NonlinearStaticExecutionResult {
+  NonlinearAnalysisState state;
+  std::vector<Diagnostic> diagnostics;
+};
+
+using PreparedAnalysisOutcome =
+    std::variant<LinearStaticExecutionResult, NonlinearStaticExecutionResult>;
+
+/// @brief Selects the approved procedure from a prepared Domain.
+/// @pre The executors and log outlive this call.
+/// @return A validated procedure outcome without opening a result file.
+Result<PreparedAnalysisOutcome> ExecutePreparedAnalysis(
+    const Domain& domain, LinearStaticProcedureExecutor& linear_executor,
+    NonlinearStaticProcedureExecutor& nonlinear_executor,
+    NewtonIterationLogWriter& log);
+
+/// @brief Writes an already computed outcome to its authoritative result file.
+/// @pre The Newton log has closed successfully before this function is called.
+/// @return Success after the selected writer atomically finalizes its result.
+Status FinalizePreparedAnalysisResult(
+    const DerivedOutputPaths& paths, const Domain& domain,
+    const PreparedAnalysisOutcome& outcome,
+    const std::vector<Diagnostic>& diagnostics, ResultsWriter& linear_writer,
+    NonlinearResultsWriter& nonlinear_writer);
+
+}  // namespace fesa
+
+#endif  // FESA_ANALYSIS_PREPARED_ANALYSIS_H_

src/fesa/CMakeLists.txt

@@ -4,9 +4,12 @@ add_library(
     analysis/analysis_model.cpp
     analysis/analysis_state.cpp
     analysis/linear_static_analysis.cpp
+    analysis/linear_static_procedure.cpp
     analysis/newton_iteration_log.cpp
     analysis/nonlinear_analysis_state.cpp
     analysis/nonlinear_static_control.cpp
+    analysis/nonlinear_static_procedure.cpp
+    analysis/prepared_analysis.cpp
     app/fesa_application.cpp
     app/output_paths.cpp
     assembly/load_assembler.cpp

src/fesa/analysis/linear_static_analysis.cpp

@@ -2,50 +2,12 @@
 
 #include <utility>
 
-#include "fesa/assembly/load_assembler.h"
-#include "fesa/assembly/parallel_for.h"
-#include "fesa/assembly/sparse_assembler.h"
-#include "fesa/elements/element_factory.h"
-#include "fesa/io/abaqus/domain_mapper.h"
-#include "fesa/io/abaqus/input_reader.h"
-#include "fesa/results/result_recovery.h"
+#include "fesa/analysis/linear_static_procedure.h"
+#include "fesa/analysis/prepared_analysis.h"
 #include "fesa/results/results_writer.h"
-#include "fesa/solvers/linear/linear_solver.h"
 
 namespace fesa {
 
-Status LinearStaticAnalysis::Run(const AnalysisRequest& request) {
-  Status status = InitializeCandidate(request);
-  if (!status.IsOk()) {
-    return status;
-  }
-  status = BuildAnalysisModel();
-  if (!status.IsOk()) {
-    return status;
-  }
-  status = BuildDofMapAndSparsePattern();
-  if (!status.IsOk()) {
-    return status;
-  }
-  status = AssembleAndPartitionStiffness();
-  if (!status.IsOk()) {
-    return status;
-  }
-  status = FactorizeFreeSystem();
-  if (!status.IsOk()) {
-    return status;
-  }
-  status = AssembleLoadsAndEffectiveRhs();
-  if (!status.IsOk()) {
-    return status;
-  }
-  status = SubstituteAndReconstruct();
-  if (!status.IsOk()) {
-    return status;
-  }
-  return RecoverAndWrite();
-}
-
 LinearStaticAnalysis::LinearStaticAnalysis(const ParallelFor& parallel_for,
                                            LinearSolver& linear_solver,
                                            ResultsWriter& results_writer)
@@ -53,138 +15,23 @@ LinearStaticAnalysis::LinearStaticAnalysis(const ParallelFor& parallel_for,
       linear_solver_{linear_solver},
       results_writer_{results_writer} {}
 
-Status LinearStaticAnalysis::InitializeCandidate(
-    const AnalysisRequest& request) {
-  // Clear dependent objects in reverse ownership order so a reused analysis
-  // never exposes a view into a Domain from an earlier run.
-  effective_rhs_.reset();
-  partitioned_stiffness_.reset();
-  full_stiffness_.reset();
-  state_.reset();
-  dofs_.reset();
-  element_view_.clear();
-  elements_.clear();
-  model_.reset();
-  domain_.reset();
-  diagnostics_.clear();
-  request_ = request;
-
-  const auto parsed = AbaqusInputReader{}.Read(request_.input_path);
-  if (!parsed.HasValue()) {
-    return parsed.GetStatus();
-  }
-  auto domain = AbaqusDomainMapper{}.Map(parsed.Value());
-  if (!domain.HasValue()) {
-    return domain.GetStatus();
+Status LinearStaticAnalysis::Run(const AnalysisRequest& request) {
+  Result<PreparedAnalysisInput> prepared =
+      PrepareAnalysisInput(request.input_path);
+  if (!prepared.HasValue()) {
+    return prepared.GetStatus();
   }
 
-  domain_ = std::make_unique<Domain>(std::move(domain.Value()));
-  diagnostics_ = domain_->Warnings();
-  SortDiagnostics(diagnostics_);
-  return Status::Ok();
-}
-
-Status LinearStaticAnalysis::BuildAnalysisModel() {
-  auto model = AnalysisModel::Create(*domain_);
-  if (!model.HasValue()) {
-    return model.GetStatus();
-  }
-  model_ = std::make_unique<AnalysisModel>(std::move(model.Value()));
-  return Status::Ok();
-}
-
-Status LinearStaticAnalysis::BuildDofMapAndSparsePattern() {
-  std::vector<std::unique_ptr<Element>> elements;
-  ElementView element_view;
-  elements.reserve(model_->ActiveElements().size());
-  element_view.reserve(model_->ActiveElements().size());
-  const ElementFactory factory;
-  for (const EntityIndex element_index : model_->ActiveElements()) {
-    auto element = factory.Create(domain_->Elements()[element_index], *domain_);
-    if (!element.HasValue()) {
-      return element.GetStatus();
-    }
-    element_view.emplace_back(*element.Value());
-    elements.push_back(std::move(element.Value()));
+  LinearStaticProcedureExecutor executor(parallel_for_, linear_solver_);
+  Result<LinearStaticExecutionResult> execution =
+      executor.Execute(prepared.Value().domain);
+  if (!execution.HasValue()) {
+    return execution.GetStatus();
   }
 
-  DofManager dofs;
-  const Status dof_status = dofs.Build(*model_, element_view);
-  if (!dof_status.IsOk()) {
-    return dof_status;
-  }
-
-  elements_ = std::move(elements);
-  element_view_ = std::move(element_view);
-  dofs_ = std::make_unique<DofManager>(std::move(dofs));
-  state_ = std::make_unique<AnalysisState>(
-      AnalysisState::Create(*dofs_, {"Step-1", 0U}));
-  return Status::Ok();
-}
-
-Status LinearStaticAnalysis::AssembleAndPartitionStiffness() {
-  auto stiffness =
-      SparseAssembler::Assemble(element_view_, *dofs_, parallel_for_);
-  if (!stiffness.HasValue()) {
-    return stiffness.GetStatus();
-  }
-  full_stiffness_ =
-      std::make_unique<SparseMatrix>(std::move(stiffness.Value()));
-
-  auto partitioned =
-      EssentialConstraintPolicy{}.Partition(*full_stiffness_, *dofs_);
-  if (!partitioned.HasValue()) {
-    return partitioned.GetStatus();
-  }
-  partitioned_stiffness_ =
-      std::make_unique<PartitionedStiffness>(std::move(partitioned.Value()));
-  return Status::Ok();
-}
-
-Status LinearStaticAnalysis::FactorizeFreeSystem() {
-  // This call intentionally precedes all load assembly in Run.
-  return linear_solver_.Factorize(partitioned_stiffness_->k_ff);
-}
-
-Status LinearStaticAnalysis::AssembleLoadsAndEffectiveRhs() {
-  auto full_load = LoadAssembler::AssembleFullNodalLoad(*model_, *dofs_);
-  if (!full_load.HasValue()) {
-    return full_load.GetStatus();
-  }
-  state_->ExternalForce() = std::move(full_load.Value());
-
-  auto rhs = LoadAssembler::EffectiveFreeRhs(state_->ExternalForce(),
-                                             partitioned_stiffness_->k_fc,
-                                             dofs_->PrescribedValues(), *dofs_);
-  if (!rhs.HasValue()) {
-    return rhs.GetStatus();
-  }
-  effective_rhs_ = std::make_unique<Vector>(std::move(rhs.Value()));
-  return Status::Ok();
-}
-
-Status LinearStaticAnalysis::SubstituteAndReconstruct() {
-  Vector free_displacement{dofs_->FreeDofCount()};
-  const Status solve_status =
-      linear_solver_.Solve(*effective_rhs_, free_displacement);
-  if (!solve_status.IsOk()) {
-    return solve_status;
-  }
-
-  state_->Displacement() = EssentialConstraintPolicy{}.ReconstructFull(
-      free_displacement, dofs_->PrescribedValues(), *dofs_);
-  return Status::Ok();
-}
-
-Status LinearStaticAnalysis::RecoverAndWrite() {
-  const Status recovery_status = ResultRecovery::Recover(
-      *model_, element_view_, *dofs_, *full_stiffness_, state_->Displacement(),
-      state_->ExternalForce(), *state_);
-  if (!recovery_status.IsOk()) {
-    return recovery_status;
-  }
-  return results_writer_.Write(request_.output_path, *domain_, *state_,
-                               diagnostics_);
+  return results_writer_.Write(request.output_path, prepared.Value().domain,
+                               execution.Value().state,
+                               execution.Value().diagnostics);
 }
 
 }  // namespace fesa

src/fesa/analysis/linear_static_procedure.cpp

@@ -0,0 +1,105 @@
+#include "fesa/analysis/linear_static_procedure.h"
+
+#include <memory>
+#include <utility>
+#include <vector>
+
+#include "fesa/analysis/analysis_model.h"
+#include "fesa/assembly/load_assembler.h"
+#include "fesa/assembly/parallel_for.h"
+#include "fesa/assembly/sparse_assembler.h"
+#include "fesa/constraints/essential_constraint_policy.h"
+#include "fesa/elements/element.h"
+#include "fesa/elements/element_factory.h"
+#include "fesa/fem/dof_manager.h"
+#include "fesa/results/result_recovery.h"
+#include "fesa/solvers/linear/linear_solver.h"
+
+namespace fesa {
+
+LinearStaticProcedureExecutor::LinearStaticProcedureExecutor(
+    const ParallelFor& parallel_for, LinearSolver& linear_solver)
+    : parallel_for_(parallel_for), linear_solver_(linear_solver) {}
+
+Result<LinearStaticExecutionResult> LinearStaticProcedureExecutor::Execute(
+    const Domain& domain) const {
+  Result<AnalysisModel> model = AnalysisModel::Create(domain);
+  if (!model.HasValue()) {
+    return Result<LinearStaticExecutionResult>::Failure(model.GetStatus());
+  }
+
+  std::vector<std::unique_ptr<Element>> elements;
+  ElementView element_view;
+  elements.reserve(model.Value().ActiveElements().size());
+  element_view.reserve(model.Value().ActiveElements().size());
+  const ElementFactory factory;
+  for (const EntityIndex element_index : model.Value().ActiveElements()) {
+    Result<std::unique_ptr<Element>> element =
+        factory.Create(domain.Elements()[element_index], domain);
+    if (!element.HasValue()) {
+      return Result<LinearStaticExecutionResult>::Failure(element.GetStatus());
+    }
+    element_view.emplace_back(*element.Value());
+    elements.push_back(std::move(element.Value()));
+  }
+
+  DofManager dofs;
+  const Status dof_status = dofs.Build(model.Value(), element_view);
+  if (!dof_status.IsOk()) {
+    return Result<LinearStaticExecutionResult>::Failure(dof_status);
+  }
+
+  AnalysisState state = AnalysisState::Create(dofs, {"Step-1", 0U});
+  Result<SparseMatrix> stiffness =
+      SparseAssembler::Assemble(element_view, dofs, parallel_for_);
+  if (!stiffness.HasValue()) {
+    return Result<LinearStaticExecutionResult>::Failure(stiffness.GetStatus());
+  }
+
+  Result<PartitionedStiffness> partitions =
+      EssentialConstraintPolicy{}.Partition(stiffness.Value(), dofs);
+  if (!partitions.HasValue()) {
+    return Result<LinearStaticExecutionResult>::Failure(partitions.GetStatus());
+  }
+
+  const Status factorized = linear_solver_.Factorize(partitions.Value().k_ff);
+  if (!factorized.IsOk()) {
+    return Result<LinearStaticExecutionResult>::Failure(factorized);
+  }
+
+  Result<Vector> external =
+      LoadAssembler::AssembleFullNodalLoad(model.Value(), dofs);
+  if (!external.HasValue()) {
+    return Result<LinearStaticExecutionResult>::Failure(external.GetStatus());
+  }
+  state.ExternalForce() = std::move(external.Value());
+
+  Result<Vector> rhs = LoadAssembler::EffectiveFreeRhs(
+      state.ExternalForce(), partitions.Value().k_fc, dofs.PrescribedValues(),
+      dofs);
+  if (!rhs.HasValue()) {
+    return Result<LinearStaticExecutionResult>::Failure(rhs.GetStatus());
+  }
+
+  Vector free_displacement{dofs.FreeDofCount()};
+  const Status solved = linear_solver_.Solve(rhs.Value(), free_displacement);
+  if (!solved.IsOk()) {
+    return Result<LinearStaticExecutionResult>::Failure(solved);
+  }
+  state.Displacement() = EssentialConstraintPolicy{}.ReconstructFull(
+      free_displacement, dofs.PrescribedValues(), dofs);
+
+  const Status recovered = ResultRecovery::Recover(
+      model.Value(), element_view, dofs, stiffness.Value(),
+      state.Displacement(), state.ExternalForce(), state);
+  if (!recovered.IsOk()) {
+    return Result<LinearStaticExecutionResult>::Failure(recovered);
+  }
+
+  std::vector<Diagnostic> diagnostics = domain.Warnings();
+  SortDiagnostics(diagnostics);
+  return Result<LinearStaticExecutionResult>::Success(
+      {std::move(state), std::move(diagnostics)});
+}
+
+}  // namespace fesa

src/fesa/analysis/nonlinear_static_procedure.cpp

@@ -0,0 +1,807 @@
+#include "fesa/analysis/nonlinear_static_procedure.h"
+
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <optional>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "fesa/analysis/analysis_model.h"
+#include "fesa/analysis/newton_iteration_log.h"
+#include "fesa/analysis/nonlinear_static_control.h"
+#include "fesa/assembly/parallel_for.h"
+#include "fesa/constraints/essential_constraint_policy.h"
+#include "fesa/elements/finite_rotation.h"
+#include "fesa/elements/mitc4_shell.h"
+#include "fesa/elements/nonlinear_mitc4_shell.h"
+#include "fesa/fem/dof_manager.h"
+#include "fesa/math/sparse_matrix.h"
+#include "fesa/properties/shell_section.h"
+#include "fesa/solvers/linear/tangent_solver.h"
+
+namespace fesa {
+namespace {
+
+constexpr std::size_t kDofsPerNode = 6U;
+constexpr std::size_t kTranslationsPerNode = 3U;
+constexpr std::size_t kShellDofCount = 24U;
+constexpr std::uint32_t kMaximumNewtonIterations = 100U;
+
+struct RuntimeShell {
+  EntityIndex element_index;
+  std::array<std::size_t, kShellDofCount> scatter;
+  NonlinearMitc4Shell element;
+};
+
+struct InternalAssembly {
+  Vector force;
+  SparseMatrix tangent;
+  std::vector<ShellResultRow> shell_rows;
+  double physical_strain_energy;
+};
+
+Status ProcedureFailure(const FailureCategory category, const std::string& code,
+                        const SourceLocation& location,
+                        const std::string& keyword,
+                        const std::string& identity,
+                        const std::string& message) {
+  return Status::Failure(category, {{Severity::kError, code, location, keyword,
+                                     identity, message}});
+}
+
+std::string DiagnosticCode(const Status& status) {
+  if (status.Diagnostics().empty()) {
+    return "nonlinear-procedure-failure";
+  }
+  return status.Diagnostics().front().code;
+}
+
+std::string DiagnosticReason(const Status& status) {
+  if (status.Diagnostics().empty()) {
+    return "The nonlinear procedure returned no diagnostic.";
+  }
+  return status.Diagnostics().front().message;
+}
+
+NewtonChannelDecision UnevaluatedDecision(const std::string& identity,
+                                          const NewtonChannelUnit unit) {
+  return {identity,
+          0.0,
+          0.0,
+          NewtonConvergenceBranch::kZeroScale,
+          kNewtonZeroScaleThreshold,
+          unit,
+          false};
+}
+
+NewtonIterationLogRecord NewLogRecord(const StepDefinition& step,
+                                      const LoadFactorIncrement& increment,
+                                      const std::uint32_t iteration) {
+  return {step.Name(),
+          increment.increment_index,
+          increment.load_factor,
+          iteration,
+          NewtonIterationStage::kEntered,
+          {false, NewtonIterationStage::kEntered, {}, {}},
+          "entered",
+          "pending",
+          {},
+          UnevaluatedDecision("residual_force", NewtonChannelUnit::kForce),
+          UnevaluatedDecision("residual_moment",
+                              NewtonChannelUnit::kForceLength),
+          UnevaluatedDecision("correction_translation",
+                              NewtonChannelUnit::kLength),
+          UnevaluatedDecision("correction_rotation",
+                              NewtonChannelUnit::kDimensionless),
+          false};
+}
+
+Status WriteAndFlush(NewtonIterationLogWriter& log,
+                     const NewtonIterationLogRecord& record) {
+  const Status write_status = log.Write(record);
+  if (!write_status.IsOk()) {
+    return write_status;
+  }
+  return log.Flush();
+}
+
+Status RecordFailure(NewtonIterationLogWriter& log,
+                     NewtonIterationLogRecord* record, const Status& failure,
+                     const NewtonIterationStage stage) {
+  record->last_completed_stage = stage;
+  record->failure = {true, stage, DiagnosticCode(failure),
+                     DiagnosticReason(failure)};
+  record->transition = "failed";
+  record->outcome = "failed";
+  record->diagnostic_code = DiagnosticCode(failure);
+  const Status log_status = WriteAndFlush(log, *record);
+  return log_status.IsOk() ? failure : log_status;
+}
+
+Matrix InitialFrame(const ShellNodeInitialFrame& frame) {
+  Matrix result{3U, 3U};
+  for (std::size_t row = 0U; row < kTranslationsPerNode; ++row) {
+    result(row, 0U) = frame.tangent_a[row];
+    result(row, 1U) = frame.tangent_b[row];
+    result(row, 2U) = frame.director[row];
+  }
+  return result;
+}
+
+Result<DofManager> BuildDofManager(const Domain& domain) {
+  Result<AnalysisModel> model = AnalysisModel::Create(domain);
+  if (!model.HasValue()) {
+    return Result<DofManager>::Failure(model.GetStatus());
+  }
+  return DofManager::Create(model.Value());
+}
+
+Result<NonlinearStateBaselineInput> BuildBaseline(const Domain& domain) {
+  if (domain.ShellNodeInitialFrames().size() != domain.Nodes().size()) {
+    return Result<NonlinearStateBaselineInput>::Failure(ProcedureFailure(
+        FailureCategory::kModel, "invalid-nonlinear-state",
+        {domain.SourcePath(), 0U},
+        "NODE", domain.SourceContentIdentity(),
+        "Nonlinear analysis requires one initial shell frame per node."));
+  }
+
+  NonlinearStateBaselineInput baseline;
+  baseline.node_baselines.reserve(domain.Nodes().size());
+  for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) {
+    const ShellNodeInitialFrame& frame = domain.ShellNodeInitialFrames()[node];
+    if (frame.node_index != static_cast<EntityIndex>(node)) {
+      return Result<NonlinearStateBaselineInput>::Failure(ProcedureFailure(
+          FailureCategory::kModel, "invalid-nonlinear-state",
+          domain.Nodes()[node].location, "NODE",
+          domain.Nodes()[node].source_id.source_label_text,
+          "Initial shell frames must follow stable Domain node order."));
+    }
+    baseline.node_baselines.push_back(
+        {static_cast<EntityIndex>(node), domain.Nodes()[node].source_id,
+         {frame.tangent_a, frame.tangent_b, frame.director}, frame.director});
+  }
+  return Result<NonlinearStateBaselineInput>::Success(std::move(baseline));
+}
+
+Result<std::vector<RuntimeShell>> BuildRuntimeShells(const Domain& domain,
+                                                     const DofManager& dofs) {
+  std::vector<RuntimeShell> shells;
+  shells.reserve(domain.Elements().Size());
+  for (std::size_t index = 0U; index < domain.Elements().Size(); ++index) {
+    const auto* definition =
+        dynamic_cast<const Mitc4ShellDefinition*>(&domain.Elements()[index]);
+    if (definition == nullptr) {
+      return Result<std::vector<RuntimeShell>>::Failure(ProcedureFailure(
+          FailureCategory::kModel, "unsupported-nonlinear-element",
+          {}, "ELEMENT",
+          domain.Elements()[index].SourceId().source_label_text,
+           "The nonlinear static procedure supports MITC4 shell elements "
+           "only."));
+    }
+    if (definition->section_index >= domain.Properties().Size() ||
+        definition->material_index >= domain.Materials().Size()) {
+      return Result<std::vector<RuntimeShell>>::Failure(ProcedureFailure(
+          FailureCategory::kModel, "invalid-nonlinear-static-data",
+           definition->location, "ELEMENT",
+           definition->source_id.source_label_text,
+           "The nonlinear shell element property or material index is "
+           "invalid."));
+    }
+
+    const auto* section =
+        dynamic_cast<const ShellSection*>(
+            &domain.Properties()[definition->section_index]);
+    const auto* material = dynamic_cast<const LinearElasticMaterial*>(
+        &domain.Materials()[definition->material_index]);
+    if (section == nullptr || material == nullptr) {
+      return Result<std::vector<RuntimeShell>>::Failure(ProcedureFailure(
+          FailureCategory::kModel, "incompatible-nonlinear-shell-property",
+           definition->location, "ELEMENT",
+           definition->source_id.source_label_text,
+           "A nonlinear MITC4 shell requires a shell section and linear "
+           "elastic material."));
+    }
+
+    std::array<const Node*, 4U> nodes{};
+    std::array<ShellNodeInitialFrame, 4U> frames{};
+    for (std::size_t node = 0U; node < definition->node_indices.size();
+         ++node) {
+      const EntityIndex node_index = definition->node_indices[node];
+      if (node_index >= domain.Nodes().size() ||
+          node_index >= domain.ShellNodeInitialFrames().size()) {
+        return Result<std::vector<RuntimeShell>>::Failure(ProcedureFailure(
+            FailureCategory::kModel, "invalid-shell-director",
+            definition->location, "ELEMENT",
+            definition->source_id.source_label_text,
+            "A nonlinear shell node does not have a stable initial frame."));
+      }
+      nodes[node] = &domain.Nodes()[node_index];
+      frames[node] = domain.ShellNodeInitialFrames()[node_index];
+      if (frames[node].node_index != node_index) {
+        return Result<std::vector<RuntimeShell>>::Failure(ProcedureFailure(
+            FailureCategory::kModel, "invalid-shell-director",
+            definition->location, "ELEMENT",
+            definition->source_id.source_label_text,
+            "A nonlinear shell initial frame does not match its node "
+            "identity."));
+      }
+    }
+
+    Result<NonlinearMitc4Shell> shell =
+        NonlinearMitc4Shell::Create(nodes, frames, *section, *material);
+    if (!shell.HasValue()) {
+      return Result<std::vector<RuntimeShell>>::Failure(shell.GetStatus());
+    }
+    const EntityIndex element_index = static_cast<EntityIndex>(index);
+    shells.push_back(
+        {element_index, dofs.ShellElementScatter(element_index),
+         std::move(shell.Value())});
+  }
+  return Result<std::vector<RuntimeShell>>::Success(std::move(shells));
+}
+
+Status PopulateNodalKinematics(const Domain& domain,
+                               NonlinearTrialState* trial) {
+  if (trial->generalized_coordinates.Size() !=
+      domain.Nodes().size() * kDofsPerNode) {
+    return ProcedureFailure(
+        FailureCategory::kSolver, "nonfinite-newton-data",
+        {domain.SourcePath(), 0U},
+        "NONLINEAR_STATIC_PROCEDURE", domain.SourceContentIdentity(),
+        "Nonlinear trial coordinates do not match the Domain node inventory.");
+  }
+
+  trial->displacement = trial->generalized_coordinates;
+  trial->rotation_vector = Vector{trial->generalized_coordinates.Size()};
+  trial->nodal_triad.clear();
+  trial->nodal_director.clear();
+  trial->nodal_triad.reserve(domain.Nodes().size());
+  trial->nodal_director.reserve(domain.Nodes().size());
+  for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) {
+    const std::size_t first = node * kDofsPerNode;
+    const Vector3 phi{trial->generalized_coordinates[first + 3U],
+                      trial->generalized_coordinates[first + 4U],
+                      trial->generalized_coordinates[first + 5U]};
+    Result<FiniteRotationState> rotation = MakeFiniteRotationState(
+        phi, InitialFrame(domain.ShellNodeInitialFrames()[node]));
+    if (!rotation.HasValue()) {
+      return rotation.GetStatus();
+    }
+    for (std::size_t component = 0U; component < kTranslationsPerNode;
+         ++component) {
+      trial->rotation_vector[first + 3U + component] = phi[component];
+    }
+    const Matrix& q = rotation.Value().q;
+    trial->nodal_triad.push_back(
+        {{{q(0U, 0U), q(1U, 0U), q(2U, 0U)},
+          {q(0U, 1U), q(1U, 1U), q(2U, 1U)},
+          {q(0U, 2U), q(1U, 2U), q(2U, 2U)}}});
+    trial->nodal_director.push_back(rotation.Value().director.Components());
+  }
+  return Status::Ok();
+}
+
+NonlinearMitc4TrialDofs ExtractElementTrial(
+    const RuntimeShell& shell, const Vector& coordinates) {
+  NonlinearMitc4TrialDofs trial{};
+  for (std::size_t node = 0U; node < 4U; ++node) {
+    const std::size_t first = node * kDofsPerNode;
+    trial.displacement[node] =
+        {coordinates[shell.scatter[first]],
+         coordinates[shell.scatter[first + 1U]],
+         coordinates[shell.scatter[first + 2U]]};
+    trial.rotation_vector[node] =
+        {coordinates[shell.scatter[first + 3U]],
+         coordinates[shell.scatter[first + 4U]],
+         coordinates[shell.scatter[first + 5U]]};
+  }
+  return trial;
+}
+
+bool IsFinite(const Vector& values) {
+  for (std::size_t index = 0U; index < values.Size(); ++index) {
+    if (!std::isfinite(values[index])) {
+      return false;
+    }
+  }
+  return true;
+}
+
+bool IsFinite(const Matrix& values) {
+  for (std::size_t row = 0U; row < values.Rows(); ++row) {
+    for (std::size_t column = 0U; column < values.Columns(); ++column) {
+      if (!std::isfinite(values(row, column))) {
+        return false;
+      }
+    }
+  }
+  return true;
+}
+
+Result<InternalAssembly> AssembleInternal(
+    const std::vector<RuntimeShell>& shells, const DofManager& dofs,
+    const Vector& coordinates, const ParallelFor& parallel_for) {
+  std::vector<std::optional<NonlinearMitc4InternalEquilibrium>> local(
+      shells.size());
+  std::vector<std::optional<Status>> failures(shells.size());
+  parallel_for.Execute(shells.size(), [&](const std::size_t index) {
+    Result<NonlinearMitc4InternalEquilibrium> result =
+        EvaluateInternalResidualAndTangent(shells[index].element,
+                                           ExtractElementTrial(shells[index],
+                                                               coordinates));
+    if (!result.HasValue()) {
+      failures[index] = result.GetStatus();
+      return;
+    }
+    local[index] = std::move(result.Value());
+  });
+
+  Vector force{dofs.FullDofCount()};
+  std::vector<CooContribution> contributions;
+  contributions.reserve(shells.size() * kShellDofCount * kShellDofCount);
+  std::vector<ShellResultRow> rows;
+  rows.reserve(shells.size() * 4U);
+  double energy = 0.0;
+  for (std::size_t index = 0U; index < shells.size(); ++index) {
+    if (failures[index].has_value()) {
+      return Result<InternalAssembly>::Failure(*failures[index]);
+    }
+    if (!local[index].has_value() ||
+        !IsFinite(local[index]->internal_residual) ||
+        !IsFinite(local[index]->internal_tangent) ||
+        !std::isfinite(local[index]->physical_strain_energy)) {
+      return Result<InternalAssembly>::Failure(ProcedureFailure(
+          FailureCategory::kSolver, "nonfinite-newton-data", {},
+          "NONLINEAR_STATIC_PROCEDURE", std::to_string(index),
+          "A nonlinear shell internal residual or tangent is nonfinite."));
+    }
+    for (std::size_t row = 0U; row < kShellDofCount; ++row) {
+      force[shells[index].scatter[row]] += local[index]->internal_residual[row];
+      for (std::size_t column = 0U; column < kShellDofCount; ++column) {
+        contributions.push_back(
+            {shells[index].scatter[row], shells[index].scatter[column],
+             local[index]->internal_tangent(row, column), index,
+             row * kShellDofCount + column});
+      }
+    }
+    for (const NonlinearMitc4RecoveryPoint& point :
+         local[index]->recovery.gauss_points) {
+      rows.push_back({shells[index].element_index,
+                      point.location,
+                      point.natural_coordinates,
+                      {point.local_frame.e1, point.local_frame.e2,
+                       point.local_frame.e3},
+                      point.generalized_strain,
+                      point.section_resultant,
+                      point.stress});
+    }
+    energy += local[index]->physical_strain_energy;
+  }
+  if (!IsFinite(force) || !std::isfinite(energy)) {
+    return Result<InternalAssembly>::Failure(ProcedureFailure(
+        FailureCategory::kSolver, "nonfinite-newton-data", {},
+        "NONLINEAR_STATIC_PROCEDURE", {},
+        "Global nonlinear internal force or energy is nonfinite."));
+  }
+
+  Result<SparseMatrix> tangent = SparseMatrix::FromCoo(
+      dofs.FullDofCount(), dofs.FullDofCount(), std::move(contributions),
+      dofs.GetSparsePattern());
+  if (!tangent.HasValue()) {
+    return Result<InternalAssembly>::Failure(tangent.GetStatus());
+  }
+  return Result<InternalAssembly>::Success(
+      {std::move(force), std::move(tangent.Value()), std::move(rows), energy});
+}
+
+Result<Vector> Difference(const Vector& left, const Vector& right) {
+  if (left.Size() != right.Size()) {
+    return Result<Vector>::Failure(ProcedureFailure(
+        FailureCategory::kSolver, "nonfinite-newton-data", {},
+        "NONLINEAR_STATIC_PROCEDURE", {},
+        "Nonlinear internal and external vectors must have equal dimensions."));
+  }
+  Vector result{left.Size()};
+  for (std::size_t index = 0U; index < result.Size(); ++index) {
+    result[index] = left[index] - right[index];
+  }
+  if (!IsFinite(result)) {
+    return Result<Vector>::Failure(ProcedureFailure(
+        FailureCategory::kSolver, "nonfinite-newton-data", {},
+        "NONLINEAR_STATIC_PROCEDURE", {},
+        "Nonlinear residual must remain finite."));
+  }
+  return Result<Vector>::Success(std::move(result));
+}
+
+Result<SparseMatrix> Difference(const SparseMatrix& internal,
+                                const SparseMatrix& external,
+                                const DofManager& dofs) {
+  if (internal.Rows() != external.Rows() ||
+      internal.Columns() != external.Columns() ||
+      internal.RowOffsets() != external.RowOffsets() ||
+      internal.ColumnIndices() != external.ColumnIndices()) {
+    return Result<SparseMatrix>::Failure(ProcedureFailure(
+        FailureCategory::kSolver, "nonfinite-newton-data", {},
+        "NONLINEAR_STATIC_PROCEDURE", {},
+        "Nonlinear internal and external tangent patterns must match."));
+  }
+  std::vector<CooContribution> contributions;
+  contributions.reserve(internal.Values().size());
+  for (std::size_t row = 0U; row < internal.Rows(); ++row) {
+    for (std::size_t offset = internal.RowOffsets()[row];
+         offset < internal.RowOffsets()[row + 1U]; ++offset) {
+      const double value =
+          internal.Values()[offset] - external.Values()[offset];
+      if (!std::isfinite(value)) {
+        return Result<SparseMatrix>::Failure(ProcedureFailure(
+            FailureCategory::kSolver, "nonfinite-newton-data", {},
+            "NONLINEAR_STATIC_PROCEDURE", {},
+            "Nonlinear tangent must remain finite."));
+      }
+      contributions.push_back(
+          {row, internal.ColumnIndices()[offset], value, row, offset});
+    }
+  }
+  return SparseMatrix::FromCoo(internal.Rows(), internal.Columns(),
+                               std::move(contributions),
+                               dofs.GetSparsePattern());
+}
+
+Vector SelectFreeComponents(const Vector& full, const DofManager& dofs,
+                            const bool rotations) {
+  std::size_t count = 0U;
+  for (const std::size_t full_dof : dofs.FreeDofs()) {
+    if ((full_dof % kDofsPerNode >= kTranslationsPerNode) == rotations) {
+      ++count;
+    }
+  }
+  Vector selected{count};
+  std::size_t output = 0U;
+  for (const std::size_t full_dof : dofs.FreeDofs()) {
+    if ((full_dof % kDofsPerNode >= kTranslationsPerNode) == rotations) {
+      selected[output++] = full[full_dof];
+    }
+  }
+  return selected;
+}
+
+Vector SelectCorrectionComponents(const Vector& correction,
+                                  const DofManager& dofs,
+                                  const bool rotations) {
+  std::size_t count = 0U;
+  for (const std::size_t full_dof : dofs.FreeDofs()) {
+    if ((full_dof % kDofsPerNode >= kTranslationsPerNode) == rotations) {
+      ++count;
+    }
+  }
+  Vector selected{count};
+  std::size_t output = 0U;
+  for (std::size_t equation = 0U; equation < dofs.FreeDofs().size();
+       ++equation) {
+    if ((dofs.FreeDofs()[equation] % kDofsPerNode >=
+         kTranslationsPerNode) == rotations) {
+      selected[output++] = correction[equation];
+    }
+  }
+  return selected;
+}
+
+Vector CurrentIncrementComponents(const Vector& current, const Vector& start,
+                                  const DofManager& dofs,
+                                  const bool rotations) {
+  Vector difference{current.Size()};
+  for (std::size_t index = 0U; index < difference.Size(); ++index) {
+    difference[index] = current[index] - start[index];
+  }
+  return SelectFreeComponents(difference, dofs, rotations);
+}
+
+std::array<double, 6> GlobalEquilibrium(const Domain& domain,
+                                        const Vector& residual) {
+  std::array<double, 6> equilibrium{};
+  for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) {
+    const std::size_t first = node * kDofsPerNode;
+    const Vector3 force{residual[first], residual[first + 1U],
+                        residual[first + 2U]};
+    const Vector3 position{domain.Nodes()[node].coordinates};
+    const Vector3 moment = position.Cross(force);
+    equilibrium[0U] += force.X();
+    equilibrium[1U] += force.Y();
+    equilibrium[2U] += force.Z();
+    equilibrium[3U] += moment.X() + residual[first + 3U];
+    equilibrium[4U] += moment.Y() + residual[first + 4U];
+    equilibrium[5U] += moment.Z() + residual[first + 5U];
+  }
+  return equilibrium;
+}
+
+Status ApplyIncrementConstraints(const StepDefinition& step,
+                                 const LoadFactorIncrement& increment,
+                                 const DofManager& dofs,
+                                 NonlinearTrialState* trial) {
+  Result<Vector> source_targets =
+      ApplyAbsolutePrescribedDisplacementTarget(step, increment.load_factor);
+  if (!source_targets.HasValue()) {
+    return source_targets.GetStatus();
+  }
+  for (std::size_t index = 0U; index < dofs.ConstrainedDofs().size();
+       ++index) {
+    const double target =
+        increment.load_factor * dofs.PrescribedValues()[index];
+    if (!std::isfinite(target)) {
+      return ProcedureFailure(
+          FailureCategory::kSolver, "nonfinite-newton-data", step.Location(),
+          "BOUNDARY", step.Name(),
+          "A nonlinear prescribed displacement target is nonfinite.");
+    }
+    trial->generalized_coordinates[dofs.ConstrainedDofs()[index]] = target;
+  }
+  return Status::Ok();
+}
+
+AcceptedFrameState AcceptedFrame(const NonlinearAnalysisState& state,
+                                 const NonlinearTrialState& trial,
+                                 const LoadFactorIncrement& increment,
+                                 const std::uint32_t iteration) {
+  return {static_cast<std::uint32_t>(state.AcceptedFrames().size()),
+          increment.increment_index,
+          increment.load_factor,
+          iteration,
+          increment.is_factor_one,
+          trial.nodal_internal_order,
+          trial.nodal_source_ids,
+          trial.generalized_coordinates,
+          trial.displacement,
+          trial.rotation_vector,
+          trial.nodal_triad,
+          trial.nodal_director,
+          trial.internal_force,
+          trial.external_force,
+          trial.residual,
+          trial.reaction,
+          trial.shell_results,
+          trial.physical_strain_energy,
+          trial.global_equilibrium};
+}
+
+}  // namespace
+
+NonlinearStaticProcedureExecutor::NonlinearStaticProcedureExecutor(
+    const ParallelFor& parallel_for, TangentSolver& tangent_solver)
+    : parallel_for_(parallel_for), tangent_solver_(tangent_solver) {}
+
+Result<NonlinearStaticExecutionResult>
+NonlinearStaticProcedureExecutor::Execute(
+    const Domain& domain, NewtonIterationLogWriter& log) const {
+  if (domain.Steps().Size() != 1U || !domain.Steps()[0U].IsNonlinearStatic()) {
+    return Result<NonlinearStaticExecutionResult>::Failure(ProcedureFailure(
+        FailureCategory::kInput, "invalid-nonlinear-static-data",
+        {domain.SourcePath(), 0U}, "STATIC", domain.SourceContentIdentity(),
+        "Nonlinear procedure execution requires exactly one NLGEOM static "
+        "step."));
+  }
+  const StepDefinition& step = domain.Steps()[0U];
+
+  Result<DofManager> dofs = BuildDofManager(domain);
+  if (!dofs.HasValue()) {
+    return Result<NonlinearStaticExecutionResult>::Failure(dofs.GetStatus());
+  }
+  Result<NonlinearStateBaselineInput> baseline = BuildBaseline(domain);
+  if (!baseline.HasValue()) {
+    return Result<NonlinearStaticExecutionResult>::Failure(
+        baseline.GetStatus());
+  }
+  Result<NonlinearAnalysisState> state =
+      NonlinearAnalysisState::Create(domain, dofs.Value(), baseline.Value());
+  if (!state.HasValue()) {
+    return Result<NonlinearStaticExecutionResult>::Failure(state.GetStatus());
+  }
+  Result<std::vector<RuntimeShell>> shells =
+      BuildRuntimeShells(domain, dofs.Value());
+  if (!shells.HasValue()) {
+    return Result<NonlinearStaticExecutionResult>::Failure(shells.GetStatus());
+  }
+  Result<std::vector<LoadFactorIncrement>> increments =
+      BuildFixedLoadFactorSequence(step);
+  if (!increments.HasValue()) {
+    return Result<NonlinearStaticExecutionResult>::Failure(
+        increments.GetStatus());
+  }
+
+  const EssentialConstraintPolicy constraints;
+  for (const LoadFactorIncrement& increment : increments.Value()) {
+    Result<NonlinearTrialState> trial = state.Value().BeginTrialFromCommitted();
+    if (!trial.HasValue()) {
+      return Result<NonlinearStaticExecutionResult>::Failure(trial.GetStatus());
+    }
+    const Vector increment_start = trial.Value().generalized_coordinates;
+    Status constraint_status =
+        ApplyIncrementConstraints(step, increment, dofs.Value(),
+                                  &trial.Value());
+    if (!constraint_status.IsOk()) {
+      return Result<NonlinearStaticExecutionResult>::Failure(constraint_status);
+    }
+
+    for (std::uint32_t iteration = 1U; iteration <= kMaximumNewtonIterations;
+         ++iteration) {
+      NewtonIterationLogRecord record =
+          NewLogRecord(step, increment, iteration);
+      trial.Value().increment_index = increment.increment_index;
+      trial.Value().load_factor = increment.load_factor;
+      trial.Value().is_factor_zero_baseline = false;
+
+      Status kinematics = PopulateNodalKinematics(domain, &trial.Value());
+      if (!kinematics.IsOk()) {
+        return Result<NonlinearStaticExecutionResult>::Failure(
+            RecordFailure(log, &record, kinematics,
+                          NewtonIterationStage::kKinematicsEvaluated));
+      }
+      record.last_completed_stage = NewtonIterationStage::kKinematicsEvaluated;
+
+      Result<InternalAssembly> internal = AssembleInternal(
+          shells.Value(), dofs.Value(), trial.Value().generalized_coordinates,
+          parallel_for_);
+      if (!internal.HasValue()) {
+        return Result<NonlinearStaticExecutionResult>::Failure(RecordFailure(
+            log, &record, internal.GetStatus(),
+            NewtonIterationStage::kInternalAssembled));
+      }
+      record.last_completed_stage = NewtonIterationStage::kInternalAssembled;
+
+      Result<GlobalFixedLoadContribution> external =
+          AssembleGlobalFixedLoadContribution(
+              domain, dofs.Value(), trial.Value().generalized_coordinates,
+              increment.load_factor);
+      if (!external.HasValue()) {
+        return Result<NonlinearStaticExecutionResult>::Failure(RecordFailure(
+            log, &record, external.GetStatus(),
+            NewtonIterationStage::kExternalAssembled));
+      }
+      record.last_completed_stage = NewtonIterationStage::kExternalAssembled;
+
+      Result<Vector> residual =
+          Difference(internal.Value().force,
+                     external.Value().generalized_external_force);
+      if (!residual.HasValue()) {
+        return Result<NonlinearStaticExecutionResult>::Failure(RecordFailure(
+            log, &record, residual.GetStatus(),
+            NewtonIterationStage::kExternalAssembled));
+      }
+      Result<SparseMatrix> tangent =
+          Difference(internal.Value().tangent,
+                     external.Value().generalized_external_tangent,
+                     dofs.Value());
+      if (!tangent.HasValue()) {
+        return Result<NonlinearStaticExecutionResult>::Failure(RecordFailure(
+            log, &record, tangent.GetStatus(),
+            NewtonIterationStage::kExternalAssembled));
+      }
+      trial.Value().internal_force = std::move(internal.Value().force);
+      trial.Value().external_force =
+          std::move(external.Value().generalized_external_force);
+      trial.Value().residual = std::move(residual.Value());
+      trial.Value().reaction = trial.Value().residual;
+      trial.Value().shell_results = std::move(internal.Value().shell_rows);
+      trial.Value().physical_strain_energy =
+          internal.Value().physical_strain_energy;
+      trial.Value().global_equilibrium =
+          GlobalEquilibrium(domain, trial.Value().residual);
+
+      Result<PartitionedStiffness> partition =
+          constraints.Partition(tangent.Value(), dofs.Value());
+      if (!partition.HasValue()) {
+        return Result<NonlinearStaticExecutionResult>::Failure(RecordFailure(
+            log, &record, partition.GetStatus(),
+            NewtonIterationStage::kExternalAssembled));
+      }
+      const Status factorization = tangent_solver_.FactorizeFresh(
+          partition.Value().k_ff);
+      if (!factorization.IsOk()) {
+        return Result<NonlinearStaticExecutionResult>::Failure(RecordFailure(
+            log, &record, factorization,
+            NewtonIterationStage::kSolvedCorrection));
+      }
+
+      Vector correction_rhs = constraints.GatherFree(trial.Value().residual,
+                                                      dofs.Value());
+      for (std::size_t index = 0U; index < correction_rhs.Size(); ++index) {
+        correction_rhs[index] = -correction_rhs[index];
+      }
+      Vector correction{dofs.Value().FreeDofCount()};
+      const Status solve =
+          tangent_solver_.SolveCorrection(correction_rhs, correction);
+      if (!solve.IsOk()) {
+        return Result<NonlinearStaticExecutionResult>::Failure(RecordFailure(
+            log, &record, solve, NewtonIterationStage::kSolvedCorrection));
+      }
+      record.last_completed_stage = NewtonIterationStage::kSolvedCorrection;
+
+      const NewtonConvergenceResult convergence = EvaluateNewtonConvergence({
+          {"residual_force", SelectFreeComponents(trial.Value().residual,
+                                                    dofs.Value(), false),
+           SelectFreeComponents(trial.Value().internal_force, dofs.Value(),
+                                false),
+           SelectFreeComponents(trial.Value().external_force, dofs.Value(),
+                                false)},
+          {"residual_moment", SelectFreeComponents(trial.Value().residual,
+                                                     dofs.Value(), true),
+           SelectFreeComponents(trial.Value().internal_force, dofs.Value(),
+                                true),
+           SelectFreeComponents(trial.Value().external_force, dofs.Value(),
+                                true)},
+          {"correction_translation",
+           SelectCorrectionComponents(correction, dofs.Value(), false),
+           CurrentIncrementComponents(trial.Value().generalized_coordinates,
+                                      increment_start, dofs.Value(), false)},
+          {"correction_rotation",
+           SelectCorrectionComponents(correction, dofs.Value(), true),
+           CurrentIncrementComponents(trial.Value().generalized_coordinates,
+                                      increment_start, dofs.Value(), true)}});
+      record.force_residual = convergence.force_residual;
+      record.moment_residual = convergence.moment_residual;
+      record.translation_correction = convergence.translation_correction;
+      record.rotation_correction = convergence.rotation_correction;
+      record.combined_passed = convergence.combined_passed;
+      record.last_completed_stage = NewtonIterationStage::kConvergenceEvaluated;
+
+      if (convergence.combined_passed) {
+        const Status committed = state.Value().CommitAcceptedFrame(
+            AcceptedFrame(state.Value(), trial.Value(), increment, iteration));
+        if (!committed.IsOk()) {
+          return Result<NonlinearStaticExecutionResult>::Failure(RecordFailure(
+              log, &record, committed, NewtonIterationStage::kAccepted));
+        }
+        record.last_completed_stage = NewtonIterationStage::kAccepted;
+        record.transition = "accepted";
+        record.outcome = "accepted";
+        const Status log_status = WriteAndFlush(log, record);
+        if (!log_status.IsOk()) {
+          return Result<NonlinearStaticExecutionResult>::Failure(log_status);
+        }
+        break;
+      }
+
+      if (iteration == kMaximumNewtonIterations) {
+        const Status failure = ProcedureFailure(
+            FailureCategory::kSolver, "nonlinear-nonconvergence",
+            step.Location(), "STATIC", step.Name(),
+            "The nonlinear procedure did not converge within 100 iterations.");
+        return Result<NonlinearStaticExecutionResult>::Failure(RecordFailure(
+            log, &record, failure,
+            NewtonIterationStage::kConvergenceEvaluated));
+      }
+
+      for (std::size_t equation = 0U; equation < dofs.Value().FreeDofs().size();
+           ++equation) {
+        const std::size_t full_dof = dofs.Value().FreeDofs()[equation];
+        trial.Value().generalized_coordinates[full_dof] += correction[equation];
+        if (!std::isfinite(trial.Value().generalized_coordinates[full_dof])) {
+          const Status failure = ProcedureFailure(
+              FailureCategory::kSolver, "nonfinite-newton-data", {},
+              "NONLINEAR_STATIC_PROCEDURE", std::to_string(full_dof),
+              "A nonlinear free-coordinate correction is nonfinite.");
+          return Result<NonlinearStaticExecutionResult>::Failure(RecordFailure(
+              log, &record, failure, NewtonIterationStage::kSolvedCorrection));
+        }
+      }
+      record.transition = "continue";
+      record.outcome = "not-converged";
+      const Status log_status = WriteAndFlush(log, record);
+      if (!log_status.IsOk()) {
+        return Result<NonlinearStaticExecutionResult>::Failure(log_status);
+      }
+    }
+  }
+
+  std::vector<Diagnostic> diagnostics = domain.Warnings();
+  SortDiagnostics(diagnostics);
+  return Result<NonlinearStaticExecutionResult>::Success(
+      {std::move(state.Value()), std::move(diagnostics)});
+}
+
+}  // namespace fesa

src/fesa/analysis/prepared_analysis.cpp

@@ -0,0 +1,79 @@
+#include "fesa/analysis/prepared_analysis.h"
+
+#include <utility>
+
+#include "fesa/analysis/linear_static_procedure.h"
+#include "fesa/analysis/newton_iteration_log.h"
+#include "fesa/analysis/nonlinear_static_procedure.h"
+#include "fesa/app/output_paths.h"
+#include "fesa/io/abaqus/domain_mapper.h"
+#include "fesa/io/abaqus/input_reader.h"
+#include "fesa/results/nonlinear_results_writer.h"
+#include "fesa/results/results_writer.h"
+
+namespace fesa {
+
+Result<PreparedAnalysisInput> PrepareAnalysisInput(
+    const std::filesystem::path& input_path) {
+  const Result<ParsedInput> syntax = AbaqusInputReader{}.Read(input_path);
+  if (!syntax.HasValue()) {
+    return Result<PreparedAnalysisInput>::Failure(syntax.GetStatus());
+  }
+
+  Result<Domain> mapped = AbaqusDomainMapper{}.Map(syntax.Value());
+  if (!mapped.HasValue()) {
+    return Result<PreparedAnalysisInput>::Failure(mapped.GetStatus());
+  }
+
+  PreparedAnalysisInput prepared{std::move(mapped.Value()), {}};
+  prepared.diagnostics = prepared.domain.Warnings();
+  SortDiagnostics(prepared.diagnostics);
+  return Result<PreparedAnalysisInput>::Success(std::move(prepared));
+}
+
+Result<PreparedAnalysisOutcome> ExecutePreparedAnalysis(
+    const Domain& domain, LinearStaticProcedureExecutor& linear_executor,
+    NonlinearStaticProcedureExecutor& nonlinear_executor,
+    NewtonIterationLogWriter& log) {
+  if (domain.Steps().Size() != 1U) {
+    return Result<PreparedAnalysisOutcome>::Failure(Status::Failure(
+        FailureCategory::kInput,
+        {{Severity::kError, "invalid-analysis-step-count", {}, "*STEP", {},
+          "Prepared analysis requires exactly one static step."}}));
+  }
+
+  if (domain.Steps()[0U].IsNonlinearStatic()) {
+    Result<NonlinearStaticExecutionResult> nonlinear =
+        nonlinear_executor.Execute(domain, log);
+    if (!nonlinear.HasValue()) {
+      return Result<PreparedAnalysisOutcome>::Failure(nonlinear.GetStatus());
+    }
+    return Result<PreparedAnalysisOutcome>::Success(
+        PreparedAnalysisOutcome{std::move(nonlinear.Value())});
+  }
+
+  Result<LinearStaticExecutionResult> linear = linear_executor.Execute(domain);
+  if (!linear.HasValue()) {
+    return Result<PreparedAnalysisOutcome>::Failure(linear.GetStatus());
+  }
+  return Result<PreparedAnalysisOutcome>::Success(
+      PreparedAnalysisOutcome{std::move(linear.Value())});
+}
+
+Status FinalizePreparedAnalysisResult(
+    const DerivedOutputPaths& paths, const Domain& domain,
+    const PreparedAnalysisOutcome& outcome,
+    const std::vector<Diagnostic>& diagnostics, ResultsWriter& linear_writer,
+    NonlinearResultsWriter& nonlinear_writer) {
+  if (std::holds_alternative<LinearStaticExecutionResult>(outcome)) {
+    return linear_writer.Write(
+        paths.results_path, domain,
+        std::get<LinearStaticExecutionResult>(outcome).state, diagnostics);
+  }
+
+  return nonlinear_writer.Write(
+      paths.results_path, domain,
+      std::get<NonlinearStaticExecutionResult>(outcome).state, diagnostics);
+}
+
+}  // namespace fesa

src/fesa/app/fesa_application.cpp

@@ -2,17 +2,20 @@
 
 #include <filesystem>
 #include <iostream>
-#include <memory>
 #include <string>
 #include <vector>
 
-#include "fesa/analysis/analysis.h"
-#include "fesa/analysis/linear_static_analysis.h"
+#include "fesa/analysis/linear_static_procedure.h"
+#include "fesa/analysis/newton_iteration_log.h"
+#include "fesa/analysis/nonlinear_static_procedure.h"
+#include "fesa/analysis/prepared_analysis.h"
 #include "fesa/app/output_paths.h"
 #include "fesa/assembly/parallel_for.h"
 #include "fesa/core/diagnostic.h"
 #include "fesa/io/hdf5/hdf5_results_writer.h"
+#include "fesa/io/hdf5/nonlinear_hdf5_results_writer.h"
 #include "fesa/solvers/linear/mkl_pardiso_solver.h"
+#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h"
 
 namespace fesa {
 namespace {
@@ -73,17 +76,53 @@ int FesaApplication::Run(const std::vector<std::string>& arguments) {
     return kOutputExitCode;
   }
 
-  AnalysisRequest request;
-  request.input_path = paths.Value().input_path;
-  request.output_path = paths.Value().results_path;
-  request.log_path = paths.Value().log_path;
-
   TbbParallelFor parallel_for;
   MklPardisoSolver linear_solver;
-  Hdf5ResultsWriter results_writer;
-  std::unique_ptr<Analysis> analysis = std::make_unique<LinearStaticAnalysis>(
-      parallel_for, linear_solver, results_writer);
-  const Status status = analysis->Run(request);
+  MklPardisoTangentSolver tangent_solver;
+  Hdf5ResultsWriter linear_writer;
+  NonlinearHdf5ResultsWriter nonlinear_writer;
+  NewtonIterationLogWriter log;
+
+  Status status = log.Open(paths.Value().log_path);
+  if (!status.IsOk()) {
+    WriteDiagnostics(status.Diagnostics());
+    return ExitCodeFor(status);
+  }
+
+  Result<PreparedAnalysisInput> prepared =
+      PrepareAnalysisInput(paths.Value().input_path);
+  if (!prepared.HasValue()) {
+    status = prepared.GetStatus();
+  } else {
+    LinearStaticProcedureExecutor linear_executor(parallel_for, linear_solver);
+    NonlinearStaticProcedureExecutor nonlinear_executor(parallel_for,
+                                                         tangent_solver);
+    Result<PreparedAnalysisOutcome> outcome = ExecutePreparedAnalysis(
+        prepared.Value().domain, linear_executor, nonlinear_executor, log);
+    if (!outcome.HasValue()) {
+      status = outcome.GetStatus();
+    } else {
+      const Status close_status = log.Close();
+      if (!close_status.IsOk()) {
+        WriteDiagnostics(close_status.Diagnostics());
+        return ExitCodeFor(close_status);
+      }
+      status = FinalizePreparedAnalysisResult(
+          paths.Value(), prepared.Value().domain, outcome.Value(),
+          prepared.Value().diagnostics, linear_writer, nonlinear_writer);
+      if (status.IsOk()) {
+        return kSuccessExitCode;
+      }
+      WriteDiagnostics(status.Diagnostics());
+      return ExitCodeFor(status);
+    }
+  }
+
+  const Status close_status = log.Close();
+  if (!close_status.IsOk()) {
+    WriteDiagnostics(close_status.Diagnostics());
+    return ExitCodeFor(close_status);
+  }
   if (status.IsOk()) {
     return kSuccessExitCode;
   }

tests/CMakeLists.txt

@@ -72,6 +72,7 @@ target_link_libraries(
 add_executable(
     fesa_integration_tests
     integration/analysis/linear_static_analysis_test.cpp
+    integration/analysis/nonlinear_static_procedure_test.cpp
     integration/app/fesa_application_test.cpp
 )
 

tests/integration/analysis/linear_static_analysis_test.cpp

@@ -16,6 +16,7 @@
 
 #include "fesa/analysis/analysis.h"
 #include "fesa/assembly/parallel_for.h"
+#include "fesa/math/sparse_matrix.h"
 #include "fesa/results/results_writer.h"
 #include "fesa/solvers/linear/mkl_pardiso_solver.h"
 

tests/integration/analysis/nonlinear_static_procedure_test.cpp

@@ -0,0 +1,157 @@
+#include "fesa/app/fesa_application.h"
+
+#include <gtest/gtest.h>
+#include <hdf5.h>
+
+#include <atomic>
+#include <chrono>
+#include <filesystem>
+#include <fstream>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+namespace {
+
+class TempDirectory {
+ public:
+  explicit TempDirectory(const std::string& label) {
+    static std::atomic<unsigned long long> sequence{0U};
+    const auto tick =
+        std::chrono::steady_clock::now().time_since_epoch().count();
+    path_ = std::filesystem::temp_directory_path() /
+            ("fesa-step10-nonlinear-" + label + "-" +
+             std::to_string(tick) + "-" +
+             std::to_string(sequence.fetch_add(1U)));
+    std::error_code error;
+    if (!std::filesystem::create_directory(path_, error) || error) {
+      throw std::runtime_error{"Unable to create nonlinear procedure fixture."};
+    }
+  }
+
+  TempDirectory(const TempDirectory&) = delete;
+  TempDirectory& operator=(const TempDirectory&) = delete;
+
+  ~TempDirectory() {
+    std::error_code ignored;
+    std::filesystem::remove_all(path_, ignored);
+  }
+
+  const std::filesystem::path& Path() const noexcept { return path_; }
+
+ private:
+  std::filesystem::path path_;
+};
+
+void WriteText(const std::filesystem::path& path, const std::string& text) {
+  std::ofstream stream{path, std::ios::binary | std::ios::trunc};
+  stream.write(text.data(), static_cast<std::streamsize>(text.size()));
+  if (!stream) {
+    throw std::runtime_error{"Unable to write nonlinear procedure input."};
+  }
+}
+
+std::string AllConstrainedNonlinearShellDeck() {
+  return R"inp(*Part, name=ShellPart
+*Node
+1, 0., 0., 0.
+2, 1., 0., 0.
+3, 1., 1., 0.
+4, 0., 1., 0.
+*Element, type=S4
+1, 1, 2, 3, 4
+*Elset, elset=ShellSet
+1
+*Shell Section, elset=ShellSet, material=Steel
+0.1
+*End Part
+*Assembly, name=Assembly
+*Instance, name=Shell-1, part=ShellPart
+*End Instance
+*Nset, nset=All, instance=Shell-1
+1, 2, 3, 4
+*End Assembly
+*Material, name=Steel
+*Elastic
+1000., 0.25
+*Boundary
+All, 1, 6
+*Step, name=Step-1, nlgeom=YES
+*Static
+0.1, 1., 0.001, 0.1
+*End Step
+)inp";
+}
+
+class Hdf5File {
+ public:
+  explicit Hdf5File(const std::filesystem::path& path)
+      : handle_{H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT)} {
+    if (handle_ < 0) {
+      throw std::runtime_error{"Unable to open nonlinear results HDF5."};
+    }
+  }
+
+  Hdf5File(const Hdf5File&) = delete;
+  Hdf5File& operator=(const Hdf5File&) = delete;
+
+  ~Hdf5File() { (void)H5Fclose(handle_); }
+
+  hid_t Get() const noexcept { return handle_; }
+
+ private:
+  hid_t handle_;
+};
+
+std::string ReadStringAttribute(const hid_t object, const char* name) {
+  const hid_t attribute = H5Aopen(object, name, H5P_DEFAULT);
+  const hid_t type = attribute < 0 ? -1 : H5Aget_type(attribute);
+  if (attribute < 0 || type < 0) {
+    if (type >= 0) {
+      (void)H5Tclose(type);
+    }
+    if (attribute >= 0) {
+      (void)H5Aclose(attribute);
+    }
+    throw std::runtime_error{"Unable to read nonlinear metadata attribute."};
+  }
+  char* value = nullptr;
+  const herr_t read = H5Aread(attribute, type, &value);
+  const std::string result = read >= 0 && value != nullptr ? value : "";
+  if (value != nullptr) {
+    (void)H5free_memory(value);
+  }
+  (void)H5Tclose(type);
+  (void)H5Aclose(attribute);
+  if (read < 0) {
+    throw std::runtime_error{"Unable to read nonlinear metadata text."};
+  }
+  return result;
+}
+
+// NSNR-PROC-004, NSNR-PROC-007, and NSNR-PROC-010.
+TEST(NonlinearStaticProcedure,
+     AllConstrainedNlgeomRunWritesSchemaV1AfterAcceptedFrames) {
+  TempDirectory directory{"schema-v1"};
+  const auto input = directory.Path() / "shellnl.inp";
+  const auto output = directory.Path() / "shellnl_results.h5";
+  WriteText(input, AllConstrainedNonlinearShellDeck());
+
+  testing::internal::CaptureStderr();
+  const int exit_code = fesa::FesaApplication{}.Run({input.string()});
+  const std::string standard_error = testing::internal::GetCapturedStderr();
+
+  ASSERT_EQ(exit_code, 0) << standard_error;
+  ASSERT_TRUE(std::filesystem::exists(output));
+  ASSERT_GT(H5Fis_hdf5(output.string().c_str()), 0);
+  const Hdf5File file{output};
+  const hid_t metadata = H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT);
+  ASSERT_GE(metadata, 0);
+  EXPECT_EQ(ReadStringAttribute(metadata, "feature_id"),
+            "nonlinear-static-newton-raphson");
+  (void)H5Gclose(metadata);
+  EXPECT_GT(H5Lexists(file.Get(), "/steps/Step-1/accepted_frames", H5P_DEFAULT),
+            0);
+}
+
+}  // namespace

tests/integration/app/fesa_application_test.cpp

@@ -368,6 +368,20 @@ TEST(LinearStaticCli, DefaultOutputUsesInputSiblingDerivedName) {
   EXPECT_FALSE(std::filesystem::exists(directory.Path() / "results.h5"));
 }
 
+// NSNR-PROC-001: log lifecycle begins after leasing and before parse/map.
+TEST(FesaApplication, OpensDerivedLogBeforeInputMapping) {
+  TempDirectory directory{"log-before-input"};
+  const auto input = directory.Path() / "missing.inp";
+  const auto log = directory.Path() / "missing_log.txt";
+
+  const auto result = RunApplication({input.string()});
+
+  EXPECT_EQ(result.exit_code, 3);
+  ExpectDiagnosticFieldOrder(result.standard_error);
+  ASSERT_TRUE(std::filesystem::exists(log));
+  EXPECT_EQ(std::filesystem::file_size(log), 0U);
+}
+
 TEST(LinearStaticCli, ReturnsEveryExactExitCodeAndOrderedDiagnostic) {
   TempDirectory directory{"exit-codes"};
   const auto valid_input = directory.Path() / "valid.inp";

tests/reference/b33_reference_comparison_test.cpp

@@ -252,6 +252,7 @@ TEST(B33ReferenceComparison, GeneratesAuthoritativeHdf5AndComparisonEvidence) {
   std::sort(generated_names.begin(), generated_names.end());
   EXPECT_EQ(generated_names,
             (std::vector<std::string>{"cantilever beam.inp",
+                                      "cantilever beam_log.txt",
                                       "cantilever beam_results.h5",
                                       "comparison.json"}));
 

tests/reference/mitc4_reference_cases_test.cpp

@@ -111,6 +111,7 @@ CaseEvidence RunCase(const std::string& case_id,
   EXPECT_EQ(generated,
             (std::vector<std::string>{"comparison.json",
                                       input.filename().string(),
+                                      input.stem().string() + "_log.txt",
                                       results.filename().string()}));
   EXPECT_TRUE(std::filesystem::is_directory(reference_directory));
   ExpectUnchanged(input, input_before);