feat(nonlinear-static-newton-raphson): step 5 - load-boundary-path-and-convergence-gates

include/fesa/analysis/nonlinear_static_control.h

@@ -0,0 +1,131 @@
+#ifndef FESA_ANALYSIS_NONLINEAR_STATIC_CONTROL_H_
+#define FESA_ANALYSIS_NONLINEAR_STATIC_CONTROL_H_
+
+#include <cstdint>
+#include <string>
+#include <vector>
+
+#include "fesa/core/diagnostic.h"
+#include "fesa/core/status.h"
+#include "fesa/fem/dof_manager.h"
+#include "fesa/math/sparse_matrix.h"
+#include "fesa/math/vector.h"
+#include "fesa/math/vector3.h"
+#include "fesa/model/domain.h"
+
+namespace fesa {
+
+/// @brief Identifies one deterministic nonlinear load-factor increment.
+struct LoadFactorIncrement {
+  std::uint32_t increment_index;
+  double load_factor;
+  bool is_factor_one;
+};
+
+/// @brief Identifies the physical dimension of a convergence threshold.
+enum class NewtonChannelUnit {
+  kForce,
+  kForceLength,
+  kLength,
+  kDimensionless,
+};
+
+inline constexpr double kNewtonRelativeThreshold = 1.0e-2;
+inline constexpr double kNewtonZeroScaleThreshold = 1.0e-3;
+
+/// @brief Carries one same-trial free residual convergence channel.
+struct NewtonResidualChannelInput {
+  std::string identity;
+  Vector free_residual;
+  Vector free_internal;
+  Vector free_external;
+};
+
+/// @brief Carries one same-iteration free correction convergence channel.
+struct NewtonCorrectionChannelInput {
+  std::string identity;
+  Vector fresh_correction;
+  Vector current_increment;
+};
+
+/// @brief Groups the four approved Newton convergence channels.
+struct NewtonConvergenceInput {
+  NewtonResidualChannelInput force_residual;
+  NewtonResidualChannelInput moment_residual;
+  NewtonCorrectionChannelInput translation_correction;
+  NewtonCorrectionChannelInput rotation_correction;
+};
+
+/// @brief Identifies the only approved Newton convergence comparison branch.
+enum class NewtonConvergenceBranch { kRelative, kZeroScale };
+
+/// @brief Records the deterministic decision for one Newton channel.
+struct NewtonChannelDecision {
+  std::string identity;
+  double numerator_norm;
+  double denominator_scale;
+  NewtonConvergenceBranch branch;
+  double threshold;
+  NewtonChannelUnit threshold_unit;
+  bool passed;
+};
+
+/// @brief Groups all four channel decisions and their logical conjunction.
+struct NewtonConvergenceResult {
+  NewtonChannelDecision force_residual;
+  NewtonChannelDecision moment_residual;
+  NewtonChannelDecision translation_correction;
+  NewtonChannelDecision rotation_correction;
+  bool combined_passed;
+};
+
+/// @brief Stores the full-space fixed-load vector and its exact tangent.
+struct GlobalFixedLoadContribution {
+  Vector generalized_external_force;
+  SparseMatrix generalized_external_tangent;
+  std::vector<Diagnostic> diagnostics;
+};
+
+/// @brief Derives the fixed nonlinear factor path from raw static controls.
+/// @param step Immutable nonlinear static step that owns raw source controls.
+/// @return Ordered factors ending in an exactly assigned factor one, or an
+/// input failure.
+Result<std::vector<LoadFactorIncrement>> BuildFixedLoadFactorSequence(
+    const StepDefinition& step);
+
+/// @brief Scales factor-one prescribed targets from the zero step-start state.
+/// @param step Immutable nonlinear step whose boundary order fixes output
+/// order.
+/// @param load_factor Current finite factor in the closed interval [0, 1].
+/// @return One absolute prescribed value per step boundary component, or an
+/// input failure.
+Result<Vector> ApplyAbsolutePrescribedDisplacementTarget(
+    const StepDefinition& step, double load_factor);
+
+/// @brief Assembles fixed-global CLOAD work once in full generalized space.
+/// @param domain Immutable nonlinear semantic model with one active step.
+/// @param dofs Stable full-DOF mapping for domain nodes.
+/// @param trial_generalized_coordinates Current full global [u, phi] trial.
+/// @param load_factor Current finite factor in the closed interval [0, 1].
+/// @return External vector and tangent for r=f_int-f_ext and K=K_int-K_ext,
+/// or a structured model/solver failure.
+Result<GlobalFixedLoadContribution> AssembleGlobalFixedLoadContribution(
+    const Domain& domain, const DofManager& dofs,
+    const Vector& trial_generalized_coordinates, double load_factor);
+
+/// @brief Rejects a nonzero source moment with director-parallel work.
+/// @param director Current finite unit director in global coordinates.
+/// @param source_moment Immutable fixed-global source moment.
+/// @return Success when rho_M is at most 1e-12 or the source moment is zero.
+Status CheckMomentAdmissibility(const Vector3& director,
+                                const Vector3& source_moment);
+
+/// @brief Evaluates the four immutable, typed Newton convergence gates.
+/// @param input Same-trial free residuals and same-iteration corrections.
+/// @return Per-channel decision records and their logical conjunction.
+NewtonConvergenceResult EvaluateNewtonConvergence(
+    const NewtonConvergenceInput& input);
+
+}  // namespace fesa
+
+#endif  // FESA_ANALYSIS_NONLINEAR_STATIC_CONTROL_H_

src/fesa/CMakeLists.txt

@@ -4,6 +4,7 @@ add_library(
     analysis/analysis_model.cpp
     analysis/analysis_state.cpp
     analysis/linear_static_analysis.cpp
+    analysis/nonlinear_static_control.cpp
     app/fesa_application.cpp
     app/output_paths.cpp
     assembly/load_assembler.cpp

src/fesa/analysis/nonlinear_static_control.cpp

@@ -0,0 +1,548 @@
+#include "fesa/analysis/nonlinear_static_control.h"
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "fesa/analysis/analysis_model.h"
+#include "fesa/assembly/load_assembler.h"
+#include "fesa/elements/finite_rotation.h"
+
+namespace fesa {
+namespace {
+
+constexpr std::size_t kDofsPerNode = 6U;
+constexpr std::size_t kTranslationDofCount = 3U;
+constexpr double kMomentAdmissibilityThreshold = 1.0e-12;
+
+/// @brief Creates a deterministic control-layer diagnostic status.
+Status ControlFailure(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}});
+}
+
+/// @brief Verifies raw nonlinear static controls without deriving state.
+Status ValidateNonlinearStaticControls(const StepDefinition& step) {
+  if (!step.IsNonlinearStatic() || !std::isfinite(step.InitialIncrement()) ||
+      !(step.InitialIncrement() > 0.0) || !std::isfinite(step.TimePeriod()) ||
+      step.TimePeriod() != 1.0 || !std::isfinite(step.MinimumIncrement()) ||
+      !(step.MinimumIncrement() > 0.0) ||
+      !std::isfinite(step.MaximumIncrement()) ||
+      !(step.MaximumIncrement() > 0.0)) {
+    return ControlFailure(
+        FailureCategory::kInput, "invalid-nonlinear-static-data",
+        step.Location(), "STATIC", step.Name(),
+        "A nonlinear static step requires four positive finite controls and "
+        "an exactly factor-one time period.");
+  }
+  return Status::Ok();
+}
+
+/// @brief Returns a finite three-component vector from full trial storage.
+Result<Vector3> RotationVectorForNode(const Vector& trial,
+                                      const std::size_t node) {
+  const std::size_t first = node * kDofsPerNode + kTranslationDofCount;
+  if (first + kTranslationDofCount > trial.Size()) {
+    return Result<Vector3>::Failure(ControlFailure(
+        FailureCategory::kSolver, "nonfinite-newton-data", {{}, 0U},
+        "NONLINEAR_STATIC_CONTROL", std::to_string(node),
+        "The full trial coordinate vector does not contain every nodal "
+        "rotation component."));
+  }
+  const Vector3 rotation{trial[first], trial[first + 1U], trial[first + 2U]};
+  if (!rotation.IsFinite()) {
+    return Result<Vector3>::Failure(ControlFailure(
+        FailureCategory::kSolver, "nonfinite-newton-data", {{}, 0U},
+        "NONLINEAR_STATIC_CONTROL", std::to_string(node),
+        "The current trial rotation coordinate must be finite."));
+  }
+  return Result<Vector3>::Success(rotation);
+}
+
+/// @brief Builds a proper initial frame from its semantic nodal record.
+Matrix InitialFrame(const ShellNodeInitialFrame& frame) {
+  Matrix result{3U, 3U};
+  for (std::size_t row = 0U; row < kTranslationDofCount; ++row) {
+    result(row, 0U) = frame.tangent_a[row];
+    result(row, 1U) = frame.tangent_b[row];
+    result(row, 2U) = frame.director[row];
+  }
+  return result;
+}
+
+/// @brief Multiplies a 3-by-3 matrix by a Vector3.
+Vector3 Multiply(const Matrix& matrix, const Vector3& vector) {
+  return {matrix(0U, 0U) * vector.X() + matrix(0U, 1U) * vector.Y() +
+              matrix(0U, 2U) * vector.Z(),
+          matrix(1U, 0U) * vector.X() + matrix(1U, 1U) * vector.Y() +
+              matrix(1U, 2U) * vector.Z(),
+          matrix(2U, 0U) * vector.X() + matrix(2U, 1U) * vector.Y() +
+              matrix(2U, 2U) * vector.Z()};
+}
+
+/// @brief Multiplies a 3-by-3 matrix transpose by a Vector3.
+Vector3 MultiplyTranspose(const Matrix& matrix, const Vector3& vector) {
+  return {matrix(0U, 0U) * vector.X() + matrix(1U, 0U) * vector.Y() +
+              matrix(2U, 0U) * vector.Z(),
+          matrix(0U, 1U) * vector.X() + matrix(1U, 1U) * vector.Y() +
+              matrix(2U, 1U) * vector.Z(),
+          matrix(0U, 2U) * vector.X() + matrix(1U, 2U) * vector.Y() +
+              matrix(2U, 2U) * vector.Z()};
+}
+
+/// @brief Forms the three-dimensional skew matrix of a global vector.
+Matrix Skew(const Vector3& vector) {
+  Matrix result{3U, 3U};
+  result(0U, 1U) = -vector.Z();
+  result(0U, 2U) = vector.Y();
+  result(1U, 0U) = vector.Z();
+  result(1U, 2U) = -vector.X();
+  result(2U, 0U) = -vector.Y();
+  result(2U, 1U) = vector.X();
+  return result;
+}
+
+/// @brief Scales a small dense matrix without changing its shape.
+Matrix Scale(const Matrix& matrix, const double scale) {
+  Matrix result{matrix.Rows(), matrix.Columns()};
+  for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
+    for (std::size_t column = 0U; column < matrix.Columns(); ++column) {
+      result(row, column) = scale * matrix(row, column);
+    }
+  }
+  return result;
+}
+
+/// @brief Adds matching small dense matrices.
+Matrix Add(const Matrix& left, const Matrix& right) {
+  Matrix result{left.Rows(), left.Columns()};
+  for (std::size_t row = 0U; row < left.Rows(); ++row) {
+    for (std::size_t column = 0U; column < left.Columns(); ++column) {
+      result(row, column) = left(row, column) + right(row, column);
+    }
+  }
+  return result;
+}
+
+/// @brief Evaluates left-Jacobian coefficients and derivatives by theta
+/// squared.
+void LeftJacobianCoefficients(const double theta_squared, double& a, double& b,
+                              double& da_dtheta_squared,
+                              double& db_dtheta_squared) {
+  if (theta_squared < 1.0e-8) {
+    const double theta_fourth = theta_squared * theta_squared;
+    const double theta_sixth = theta_fourth * theta_squared;
+    a = 0.5 - theta_squared / 24.0 + theta_fourth / 720.0 -
+        theta_sixth / 40320.0;
+    b = 1.0 / 6.0 - theta_squared / 120.0 + theta_fourth / 5040.0 -
+        theta_sixth / 362880.0;
+    da_dtheta_squared =
+        -1.0 / 24.0 + theta_squared / 360.0 - theta_fourth / 13440.0;
+    db_dtheta_squared =
+        -1.0 / 120.0 + theta_squared / 2520.0 - theta_fourth / 120960.0;
+    return;
+  }
+
+  const double theta = std::sqrt(theta_squared);
+  const double sine = std::sin(theta);
+  const double cosine = std::cos(theta);
+  a = (1.0 - cosine) / theta_squared;
+  b = (theta - sine) / (theta_squared * theta);
+  da_dtheta_squared = (theta * sine - 2.0 * (1.0 - cosine)) /
+                      (2.0 * theta_squared * theta_squared);
+  db_dtheta_squared =
+      (1.0 - cosine) / (2.0 * theta_squared * theta_squared) -
+      3.0 * (theta - sine) / (2.0 * theta_squared * theta_squared * theta);
+}
+
+/// @brief Differentiates the SO(3) left Jacobian along one coordinate vector.
+Matrix LeftJacobianDerivative(const Vector3& phi, const Vector3& direction) {
+  const Matrix skew = Skew(phi);
+  const Matrix direction_skew = Skew(direction);
+  const Matrix skew_squared = skew.Multiply(skew);
+  double a = 0.0;
+  double b = 0.0;
+  double da_dtheta_squared = 0.0;
+  double db_dtheta_squared = 0.0;
+  LeftJacobianCoefficients(phi.Dot(phi), a, b, da_dtheta_squared,
+                           db_dtheta_squared);
+  const double derivative_scale = 2.0 * phi.Dot(direction);
+  return Add(Add(Scale(skew, derivative_scale * da_dtheta_squared),
+                 Scale(direction_skew, a)),
+             Add(Scale(skew_squared, derivative_scale * db_dtheta_squared),
+                 Scale(Add(direction_skew.Multiply(skew),
+                           skew.Multiply(direction_skew)),
+                       b)));
+}
+
+/// @brief Projects the immutable source moment into current director tangent
+/// space.
+Vector3 ProjectMoment(const Vector3& director, const Vector3& source_moment) {
+  return source_moment - director.Dot(source_moment) * director;
+}
+
+/// @brief Differentiates the projected moment along a coordinate direction.
+Vector3 ProjectMomentDerivative(const Vector3& director,
+                                const Vector3& director_derivative,
+                                const Vector3& source_moment) {
+  return -director.Dot(source_moment) * director_derivative -
+         director_derivative.Dot(source_moment) * director;
+}
+
+/// @brief Computes the infinity norm and reports finite input validity.
+double InfinityNorm(const Vector& values, bool& is_finite) {
+  double result = 0.0;
+  for (std::size_t index = 0U; index < values.Size(); ++index) {
+    if (!std::isfinite(values[index])) {
+      is_finite = false;
+      return std::numeric_limits<double>::quiet_NaN();
+    }
+    result = (std::max)(result, std::abs(values[index]));
+  }
+  return result;
+}
+
+/// @brief Evaluates one residual channel with its fixed physical zero-scale
+/// unit.
+NewtonChannelDecision EvaluateResidualChannel(
+    const NewtonResidualChannelInput& input, const NewtonChannelUnit unit) {
+  bool is_finite = input.free_residual.Size() == input.free_internal.Size() &&
+                   input.free_residual.Size() == input.free_external.Size();
+  const double numerator = InfinityNorm(input.free_residual, is_finite);
+  const double internal = InfinityNorm(input.free_internal, is_finite);
+  const double external = InfinityNorm(input.free_external, is_finite);
+  const double denominator = (std::max)(internal, external);
+  if (!is_finite) {
+    return {input.identity,
+            numerator,
+            denominator,
+            NewtonConvergenceBranch::kZeroScale,
+            kNewtonZeroScaleThreshold,
+            unit,
+            false};
+  }
+  if (denominator == 0.0) {
+    return {input.identity,
+            numerator,
+            denominator,
+            NewtonConvergenceBranch::kZeroScale,
+            kNewtonZeroScaleThreshold,
+            unit,
+            numerator <= kNewtonZeroScaleThreshold};
+  }
+  return {input.identity,
+          numerator,
+          denominator,
+          NewtonConvergenceBranch::kRelative,
+          kNewtonRelativeThreshold,
+          NewtonChannelUnit::kDimensionless,
+          numerator / denominator <= kNewtonRelativeThreshold};
+}
+
+/// @brief Evaluates one correction channel with its fixed physical zero-scale
+/// unit.
+NewtonChannelDecision EvaluateCorrectionChannel(
+    const NewtonCorrectionChannelInput& input, const NewtonChannelUnit unit) {
+  bool is_finite =
+      input.fresh_correction.Size() == input.current_increment.Size();
+  const double numerator = InfinityNorm(input.fresh_correction, is_finite);
+  const double denominator = InfinityNorm(input.current_increment, is_finite);
+  if (!is_finite) {
+    return {input.identity,
+            numerator,
+            denominator,
+            NewtonConvergenceBranch::kZeroScale,
+            kNewtonZeroScaleThreshold,
+            unit,
+            false};
+  }
+  if (denominator == 0.0) {
+    return {input.identity,
+            numerator,
+            denominator,
+            NewtonConvergenceBranch::kZeroScale,
+            kNewtonZeroScaleThreshold,
+            unit,
+            numerator <= kNewtonZeroScaleThreshold};
+  }
+  return {input.identity,
+          numerator,
+          denominator,
+          NewtonConvergenceBranch::kRelative,
+          kNewtonRelativeThreshold,
+          NewtonChannelUnit::kDimensionless,
+          numerator / denominator <= kNewtonRelativeThreshold};
+}
+
+}  // namespace
+
+Result<std::vector<LoadFactorIncrement>> BuildFixedLoadFactorSequence(
+    const StepDefinition& step) {
+  const Status control_status = ValidateNonlinearStaticControls(step);
+  if (!control_status.IsOk()) {
+    return Result<std::vector<LoadFactorIncrement>>::Failure(control_status);
+  }
+
+  std::vector<LoadFactorIncrement> increments;
+  std::uint32_t increment_index = 1U;
+  while (true) {
+    const double ordinal_factor =
+        static_cast<double>(increment_index) * step.InitialIncrement();
+    if (!std::isfinite(ordinal_factor)) {
+      return Result<std::vector<LoadFactorIncrement>>::Failure(ControlFailure(
+          FailureCategory::kInput, "invalid-nonlinear-static-data",
+          step.Location(), "STATIC", step.Name(),
+          "The fixed load-factor sequence must remain finite."));
+    }
+    if (!(ordinal_factor < 1.0)) {
+      increments.push_back({increment_index, 1.0, true});
+      break;
+    }
+    increments.push_back({increment_index, ordinal_factor, false});
+    if (increment_index == (std::numeric_limits<std::uint32_t>::max)()) {
+      return Result<std::vector<LoadFactorIncrement>>::Failure(ControlFailure(
+          FailureCategory::kInput, "invalid-nonlinear-static-data",
+          step.Location(), "STATIC", step.Name(),
+          "The fixed load-factor sequence exceeds the supported increment "
+          "identity range."));
+    }
+    ++increment_index;
+  }
+  return Result<std::vector<LoadFactorIncrement>>::Success(
+      std::move(increments));
+}
+
+Result<Vector> ApplyAbsolutePrescribedDisplacementTarget(
+    const StepDefinition& step, const double load_factor) {
+  const Status control_status = ValidateNonlinearStaticControls(step);
+  if (!control_status.IsOk()) {
+    return Result<Vector>::Failure(control_status);
+  }
+  if (!std::isfinite(load_factor) || load_factor < 0.0 || load_factor > 1.0) {
+    return Result<Vector>::Failure(ControlFailure(
+        FailureCategory::kInput, "invalid-nonlinear-static-data",
+        step.Location(), "BOUNDARY", step.Name(),
+        "The prescribed-displacement path requires a finite load factor in "
+        "the closed interval from zero to one."));
+  }
+
+  const auto& boundaries = step.PrescribedDisplacements();
+  Vector targets{boundaries.Size()};
+  for (std::size_t index = 0U; index < boundaries.Size(); ++index) {
+    const double target = boundaries[index].PrescribedValue();
+    const double value = load_factor * target;
+    if (!std::isfinite(target) || !std::isfinite(value)) {
+      return Result<Vector>::Failure(ControlFailure(
+          FailureCategory::kInput, "invalid-nonlinear-static-data",
+          boundaries[index].Location(), "BOUNDARY",
+          boundaries[index].Target().target_name_or_label,
+          "The factor-one prescribed displacement target must remain finite."));
+    }
+    targets[index] = value;
+  }
+  return Result<Vector>::Success(std::move(targets));
+}
+
+Status CheckMomentAdmissibility(const Vector3& director,
+                                const Vector3& source_moment) {
+  if (!source_moment.IsFinite()) {
+    return ControlFailure(FailureCategory::kModel, "nonfinite-newton-data",
+                          {{}, 0U}, "CLOAD", "",
+                          "The fixed-global source moment must be finite.");
+  }
+  if (source_moment == Vector3{}) {
+    return Status::Ok();
+  }
+  if (!director.IsFinite() || !(director.Norm() > 0.0)) {
+    return ControlFailure(
+        FailureCategory::kSolver, "nonfinite-newton-data", {{}, 0U}, "CLOAD",
+        "", "The current shell director must be finite and nonzero.");
+  }
+  const double source_norm = source_moment.Norm();
+  const double rho = std::abs(director.Dot(source_moment)) / source_norm;
+  if (!std::isfinite(source_norm) || !(source_norm > 0.0) ||
+      !std::isfinite(rho)) {
+    return ControlFailure(
+        FailureCategory::kModel, "nonfinite-newton-data", {{}, 0U}, "CLOAD", "",
+        "The fixed-global source moment norm must be finite.");
+  }
+  if (!(rho <= kMomentAdmissibilityThreshold)) {
+    return ControlFailure(
+        FailureCategory::kModel, "unsupported-drilling-load", {{}, 0U}, "CLOAD",
+        "",
+        "The aggregate nodal moment has an unsupported director-parallel "
+        "component.");
+  }
+  return Status::Ok();
+}
+
+Result<GlobalFixedLoadContribution> AssembleGlobalFixedLoadContribution(
+    const Domain& domain, const DofManager& dofs,
+    const Vector& trial_generalized_coordinates, const double load_factor) {
+  if (!std::isfinite(load_factor) || load_factor < 0.0 || load_factor > 1.0 ||
+      domain.Steps().Size() != 1U || !domain.Steps()[0U].IsNonlinearStatic()) {
+    return Result<GlobalFixedLoadContribution>::Failure(ControlFailure(
+        FailureCategory::kModel, "invalid-nonlinear-static-data",
+        {domain.SourcePath(), 0U}, "NONLINEAR_STATIC_CONTROL",
+        domain.SourceContentIdentity(),
+        "Fixed-load assembly requires one nonlinear step and a finite load "
+        "factor in the closed interval from zero to one."));
+  }
+  const Status dof_status = dofs.ValidateInvariants();
+  if (!dof_status.IsOk()) {
+    return Result<GlobalFixedLoadContribution>::Failure(dof_status);
+  }
+  if (dofs.FullDofCount() != trial_generalized_coordinates.Size() ||
+      domain.Nodes().size() >
+          (std::numeric_limits<std::size_t>::max)() / kDofsPerNode ||
+      dofs.FullDofCount() != domain.Nodes().size() * kDofsPerNode) {
+    return Result<GlobalFixedLoadContribution>::Failure(ControlFailure(
+        FailureCategory::kSolver, "nonfinite-newton-data",
+        {domain.SourcePath(), 0U}, "NONLINEAR_STATIC_CONTROL",
+        domain.SourceContentIdentity(),
+        "The full trial coordinates must match the stable Domain DOF "
+        "dimension."));
+  }
+  for (std::size_t index = 0U; index < trial_generalized_coordinates.Size();
+       ++index) {
+    if (!std::isfinite(trial_generalized_coordinates[index])) {
+      return Result<GlobalFixedLoadContribution>::Failure(ControlFailure(
+          FailureCategory::kSolver, "nonfinite-newton-data",
+          {domain.SourcePath(), 0U}, "NONLINEAR_STATIC_CONTROL",
+          std::to_string(index),
+          "The full trial coordinate vector must contain finite values."));
+    }
+  }
+
+  auto model = AnalysisModel::Create(domain);
+  if (!model.HasValue()) {
+    return Result<GlobalFixedLoadContribution>::Failure(model.GetStatus());
+  }
+  auto source_load = LoadAssembler::AssembleFullNodalLoad(model.Value(), dofs);
+  if (!source_load.HasValue()) {
+    return Result<GlobalFixedLoadContribution>::Failure(
+        source_load.GetStatus());
+  }
+
+  Vector external_force{dofs.FullDofCount()};
+  for (std::size_t index = 0U; index < external_force.Size(); ++index) {
+    external_force[index] = load_factor * source_load.Value()[index];
+  }
+
+  std::vector<const ShellNodeInitialFrame*> frames(domain.Nodes().size(),
+                                                   nullptr);
+  for (const auto& frame : domain.ShellNodeInitialFrames()) {
+    if (frame.node_index >= frames.size() ||
+        frames[frame.node_index] != nullptr) {
+      return Result<GlobalFixedLoadContribution>::Failure(ControlFailure(
+          FailureCategory::kModel, "invalid-shell-director",
+          {domain.SourcePath(), 0U}, "NODE", std::to_string(frame.node_index),
+          "Shell nodal frames must provide one unique record per node."));
+    }
+    frames[frame.node_index] = &frame;
+  }
+
+  std::vector<CooContribution> tangent_contributions;
+  for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) {
+    const std::size_t first_rotation =
+        node * kDofsPerNode + kTranslationDofCount;
+    const Vector3 source_moment{source_load.Value()[first_rotation],
+                                source_load.Value()[first_rotation + 1U],
+                                source_load.Value()[first_rotation + 2U]};
+    if (source_moment == Vector3{}) {
+      continue;
+    }
+    if (frames[node] == nullptr) {
+      return Result<GlobalFixedLoadContribution>::Failure(ControlFailure(
+          FailureCategory::kModel, "invalid-shell-director",
+          domain.Nodes()[node].location, "NODE",
+          domain.Nodes()[node].source_id.source_label_text,
+          "A node carrying a moment requires an initial shell frame."));
+    }
+    auto phi = RotationVectorForNode(trial_generalized_coordinates, node);
+    if (!phi.HasValue()) {
+      return Result<GlobalFixedLoadContribution>::Failure(phi.GetStatus());
+    }
+    auto rotation =
+        MakeFiniteRotationState(phi.Value(), InitialFrame(*frames[node]));
+    if (!rotation.HasValue()) {
+      return Result<GlobalFixedLoadContribution>::Failure(rotation.GetStatus());
+    }
+    const Status admissibility =
+        CheckMomentAdmissibility(rotation.Value().director, source_moment);
+    if (!admissibility.IsOk()) {
+      return Result<GlobalFixedLoadContribution>::Failure(admissibility);
+    }
+
+    const Matrix left_jacobian = LeftJacobianSO3(phi.Value());
+    const Vector3 projected =
+        ProjectMoment(rotation.Value().director, source_moment);
+    const Vector3 rotational_force =
+        load_factor * MultiplyTranspose(left_jacobian, projected);
+    for (std::size_t row = 0U; row < kTranslationDofCount; ++row) {
+      external_force[first_rotation + row] = rotational_force[row];
+    }
+
+    for (std::size_t column = 0U; column < kTranslationDofCount; ++column) {
+      Vector3 direction{};
+      if (column == 0U) {
+        direction = {1.0, 0.0, 0.0};
+      } else if (column == 1U) {
+        direction = {0.0, 1.0, 0.0};
+      } else {
+        direction = {0.0, 0.0, 1.0};
+      }
+      const Matrix left_jacobian_derivative =
+          LeftJacobianDerivative(phi.Value(), direction);
+      const Vector3 director_derivative =
+          Multiply(left_jacobian, direction).Cross(rotation.Value().director);
+      const Vector3 projected_derivative = ProjectMomentDerivative(
+          rotation.Value().director, director_derivative, source_moment);
+      const Vector3 tangent_column =
+          load_factor *
+          (MultiplyTranspose(left_jacobian_derivative, projected) +
+           MultiplyTranspose(left_jacobian, projected_derivative));
+      for (std::size_t row = 0U; row < kTranslationDofCount; ++row) {
+        tangent_contributions.push_back(
+            {first_rotation + row, first_rotation + column, tangent_column[row],
+             node, row * kTranslationDofCount + column});
+      }
+    }
+  }
+
+  auto tangent = SparseMatrix::FromCoo(dofs.FullDofCount(), dofs.FullDofCount(),
+                                       std::move(tangent_contributions),
+                                       dofs.GetSparsePattern());
+  if (!tangent.HasValue()) {
+    return Result<GlobalFixedLoadContribution>::Failure(tangent.GetStatus());
+  }
+  return Result<GlobalFixedLoadContribution>::Success(
+      {std::move(external_force), std::move(tangent.Value()), {}});
+}
+
+NewtonConvergenceResult EvaluateNewtonConvergence(
+    const NewtonConvergenceInput& input) {
+  NewtonConvergenceResult result{
+      EvaluateResidualChannel(input.force_residual, NewtonChannelUnit::kForce),
+      EvaluateResidualChannel(input.moment_residual,
+                              NewtonChannelUnit::kForceLength),
+      EvaluateCorrectionChannel(input.translation_correction,
+                                NewtonChannelUnit::kLength),
+      EvaluateCorrectionChannel(input.rotation_correction,
+                                NewtonChannelUnit::kDimensionless),
+      false};
+  result.combined_passed =
+      result.force_residual.passed && result.moment_residual.passed &&
+      result.translation_correction.passed && result.rotation_correction.passed;
+  return result;
+}
+
+}  // namespace fesa

tests/CMakeLists.txt

@@ -6,6 +6,7 @@ add_executable(
     unit/analysis/analysis_test.cpp
     unit/analysis/analysis_model_test.cpp
     unit/analysis/analysis_state_test.cpp
+    unit/analysis/nonlinear_static_control_test.cpp
     unit/assembly/parallel_for_test.cpp
     unit/assembly/load_assembler_test.cpp
     unit/assembly/sparse_assembler_test.cpp

tests/unit/analysis/nonlinear_static_control_test.cpp

@@ -0,0 +1,365 @@
+#include "fesa/analysis/nonlinear_static_control.h"
+
+#include <gtest/gtest.h>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <filesystem>
+#include <memory>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "fesa/analysis/analysis_model.h"
+#include "fesa/fem/dof_manager.h"
+#include "fesa/model/domain.h"
+#include "fesa/model/model_types.h"
+
+namespace {
+
+fesa::Vector VectorFrom(const std::vector<double>& values) {
+  fesa::Vector result{values.size()};
+  for (std::size_t index = 0U; index < values.size(); ++index) {
+    result[index] = values[index];
+  }
+  return result;
+}
+
+struct Fixture {
+  std::unique_ptr<fesa::Domain> domain;
+  std::unique_ptr<fesa::AnalysisModel> model;
+  std::unique_ptr<fesa::DofManager> dofs;
+};
+
+Fixture MakeFixture(
+    std::vector<fesa::PrescribedDisplacementDefinition> boundaries = {},
+    std::vector<fesa::NodalLoad> loads = {},
+    const double initial_increment = 0.1, const double time_period = 1.0,
+    const double minimum_increment = 1.0e-9,
+    const double maximum_increment = 5.0) {
+  const std::filesystem::path source{"nonlinear-static-control.inp"};
+  fesa::ModelDefinition definition{};
+  definition.source_path = source;
+  definition.source_content_identity = "control-fixture";
+  definition.nodes = {
+      {{"Shell-1", 10, "10"}, {0.0, 0.0, 0.0}, {source, 1U}},
+      {{"Shell-1", 20, "20"}, {1.0, 0.0, 0.0}, {source, 2U}},
+      {{"Shell-1", 30, "30"}, {1.0, 1.0, 0.0}, {source, 3U}},
+      {{"Shell-1", 40, "40"}, {0.0, 1.0, 0.0}, {source, 4U}},
+      {{"Shell-1", 50, "50"}, {-1.0, 1.0, 0.0}, {source, 5U}},
+      {{"Shell-1", 60, "60"}, {-1.0, 0.0, 0.0}, {source, 6U}},
+  };
+  definition.materials = {{"Material", 1000.0, 0.25, {source, 7U}}};
+  definition.shell_sections = {{"Shell", 0.1, 0U, {source, 8U}}};
+  definition.shell_elements = {
+      {{"Shell-1", 1, "1"},
+       fesa::ShellSourceElementType::kS4,
+       {0U, 1U, 2U, 3U},
+       0U,
+       0U,
+       {source, 9U}},
+      {{"Shell-1", 2, "2"},
+       fesa::ShellSourceElementType::kS4,
+       {0U, 3U, 4U, 5U},
+       0U,
+       0U,
+       {source, 10U}},
+  };
+  for (std::size_t node = 0U; node < definition.nodes.size(); ++node) {
+    definition.shell_node_initial_frames.push_back(
+        {static_cast<fesa::EntityIndex>(node),
+         {0.0, 0.0, 1.0},
+         {1.0, 0.0, 0.0},
+         {0.0, 1.0, 0.0}});
+  }
+  definition.steps = {{"Step-1",
+                       std::move(boundaries),
+                       std::move(loads),
+                       fesa::StaticProcedureKind::kNonlinearFullNewton,
+                       true,
+                       initial_increment,
+                       time_period,
+                       minimum_increment,
+                       maximum_increment,
+                       {source, 11U}}};
+
+  auto domain_result = fesa::Domain::Create(std::move(definition));
+  if (!domain_result.HasValue()) {
+    throw std::runtime_error{"Control fixture Domain construction failed."};
+  }
+  auto domain =
+      std::make_unique<fesa::Domain>(std::move(domain_result.Value()));
+  auto model_result = fesa::AnalysisModel::Create(*domain);
+  if (!model_result.HasValue()) {
+    throw std::runtime_error{
+        "Control fixture AnalysisModel construction failed."};
+  }
+  auto model =
+      std::make_unique<fesa::AnalysisModel>(std::move(model_result.Value()));
+  auto dofs_result = fesa::DofManager::Create(*model);
+  if (!dofs_result.HasValue()) {
+    throw std::runtime_error{"Control fixture DofManager construction failed."};
+  }
+  auto dofs =
+      std::make_unique<fesa::DofManager>(std::move(dofs_result.Value()));
+  return {std::move(domain), std::move(model), std::move(dofs)};
+}
+
+fesa::Vector TrialCoordinates(const Fixture& fixture) {
+  return fesa::Vector{fixture.dofs->FullDofCount()};
+}
+
+fesa::Vector AddScaled(const fesa::Vector& values,
+                       const fesa::Vector& direction, const double scale) {
+  fesa::Vector result = values;
+  result.Axpy(scale, direction);
+  return result;
+}
+
+double MaximumAbsoluteDifference(const fesa::Vector& left,
+                                 const fesa::Vector& right) {
+  if (left.Size() != right.Size()) {
+    throw std::invalid_argument{"Vector dimensions must match."};
+  }
+  double maximum = 0.0;
+  for (std::size_t index = 0U; index < left.Size(); ++index) {
+    maximum = (std::max)(maximum, std::abs(left[index] - right[index]));
+  }
+  return maximum;
+}
+
+fesa::NewtonConvergenceInput RelativeConvergenceInput() {
+  return {
+      {"residual_force", VectorFrom({1.0e-2}), VectorFrom({1.0}),
+       VectorFrom({0.5})},
+      {"residual_moment", VectorFrom({1.0e-2}), VectorFrom({1.0}),
+       VectorFrom({0.25})},
+      {"correction_translation", VectorFrom({1.0e-2}), VectorFrom({1.0})},
+      {"correction_rotation", VectorFrom({1.0e-2}), VectorFrom({1.0})},
+  };
+}
+
+// NSNR-CTRL-001, NSNR-CTRL-002, and NSNR-CTRL-017
+TEST(NonlinearStaticControl, BuildsFixedFactorSequenceWithExactEndpoint) {
+  const auto fixture = MakeFixture();
+  const auto result =
+      fesa::BuildFixedLoadFactorSequence(fixture.domain->Steps()[0U]);
+
+  ASSERT_TRUE(result.HasValue());
+  ASSERT_EQ(result.Value().size(), 10U);
+  for (std::size_t index = 0U; index < result.Value().size(); ++index) {
+    const auto& increment = result.Value()[index];
+    EXPECT_EQ(increment.increment_index,
+              static_cast<std::uint32_t>(index + 1U));
+    EXPECT_DOUBLE_EQ(increment.load_factor,
+                     static_cast<double>(index + 1U) * 0.1);
+    EXPECT_EQ(increment.is_factor_one, index + 1U == result.Value().size());
+  }
+  EXPECT_DOUBLE_EQ(result.Value().back().load_factor, 1.0);
+
+  const auto one_hundred = MakeFixture({}, {}, 0.01);
+  const auto one_hundred_result =
+      fesa::BuildFixedLoadFactorSequence(one_hundred.domain->Steps()[0U]);
+  ASSERT_TRUE(one_hundred_result.HasValue());
+  ASSERT_EQ(one_hundred_result.Value().size(), 100U);
+  EXPECT_EQ(one_hundred_result.Value().back().increment_index, 100U);
+  EXPECT_TRUE(one_hundred_result.Value().back().is_factor_one);
+}
+
+// NSNR-CTRL-003
+TEST(NonlinearStaticControl, IgnoresPositiveMinimumAndMaximumDuringStepping) {
+  const auto first = MakeFixture({}, {}, 0.3, 1.0, 1.0e-12, 100.0);
+  const auto second = MakeFixture({}, {}, 0.3, 1.0, 0.29, 0.31);
+  const auto first_result =
+      fesa::BuildFixedLoadFactorSequence(first.domain->Steps()[0U]);
+  const auto second_result =
+      fesa::BuildFixedLoadFactorSequence(second.domain->Steps()[0U]);
+
+  ASSERT_TRUE(first_result.HasValue());
+  ASSERT_TRUE(second_result.HasValue());
+  ASSERT_EQ(first_result.Value().size(), second_result.Value().size());
+  for (std::size_t index = 0U; index < first_result.Value().size(); ++index) {
+    EXPECT_DOUBLE_EQ(first_result.Value()[index].load_factor,
+                     second_result.Value()[index].load_factor);
+    EXPECT_EQ(first_result.Value()[index].increment_index,
+              second_result.Value()[index].increment_index);
+  }
+}
+
+// NSNR-CTRL-004 and NSNR-CTRL-005
+TEST(NonlinearLoadBoundary, AppliesAbsoluteBoundaryTargetsWithoutCorrection) {
+  const std::filesystem::path source{"nonlinear-static-control.inp"};
+  const auto fixture = MakeFixture(
+      {{"10", 1, 1, 4.0, {source, 12U}}, {"20", 4, 4, -0.5, {source, 13U}}});
+  const auto result = fesa::ApplyAbsolutePrescribedDisplacementTarget(
+      fixture.domain->Steps()[0U], 0.25);
+
+  ASSERT_TRUE(result.HasValue());
+  ASSERT_EQ(result.Value().Size(), 2U);
+  EXPECT_DOUBLE_EQ(result.Value()[0U], 1.0);
+  EXPECT_DOUBLE_EQ(result.Value()[1U], -0.125);
+  for (const std::size_t constrained : fixture.dofs->ConstrainedDofs()) {
+    EXPECT_EQ(TrialCoordinates(fixture)[constrained], 0.0);
+  }
+}
+
+// NSNR-CTRL-006 and NSNR-CTRL-008
+TEST(NonlinearLoadBoundary, AssemblesEverySharedNodeLoadOnceInFixedGlobalAxes) {
+  const std::filesystem::path source{"nonlinear-static-control.inp"};
+  const std::vector<fesa::NodalLoad> loads{
+      {"10", 1, 8.0, {source, 20U}},
+      {"10", 2, -4.0, {source, 21U}},
+      {"10", 4, 2.0, {source, 22U}},
+  };
+  const auto fixture = MakeFixture({}, loads);
+  const fesa::Vector trial = TrialCoordinates(fixture);
+  const auto result = fesa::AssembleGlobalFixedLoadContribution(
+      *fixture.domain, *fixture.dofs, trial, 0.25);
+
+  ASSERT_TRUE(result.HasValue());
+  ASSERT_EQ(result.Value().generalized_external_force.Size(),
+            fixture.dofs->FullDofCount());
+  EXPECT_DOUBLE_EQ(result.Value().generalized_external_force[0U], 2.0);
+  EXPECT_DOUBLE_EQ(result.Value().generalized_external_force[1U], -1.0);
+  EXPECT_DOUBLE_EQ(result.Value().generalized_external_force[3U], 0.5);
+  EXPECT_DOUBLE_EQ(fixture.domain->Steps()[0U]
+                       .ConcentratedLoads()[0U]
+                       .GlobalComponents()[0U],
+                   8.0);
+  EXPECT_DOUBLE_EQ(fixture.domain->Steps()[0U]
+                       .ConcentratedLoads()[2U]
+                       .GlobalComponents()[3U],
+                   2.0);
+}
+
+// NSNR-CTRL-007
+TEST(NonlinearLoadBoundary, RejectsDirectorAlignedMomentAboveTheExactBand) {
+  const fesa::Vector3 director{0.0, 0.0, 1.0};
+  EXPECT_TRUE(
+      fesa::CheckMomentAdmissibility(director, {1.0, 0.0, 1.0e-12}).IsOk());
+  const auto rejected =
+      fesa::CheckMomentAdmissibility(director, {1.0, 0.0, 2.0e-12});
+  EXPECT_FALSE(rejected.IsOk());
+  ASSERT_EQ(rejected.Diagnostics().size(), 1U);
+  EXPECT_EQ(rejected.Diagnostics()[0U].code, "unsupported-drilling-load");
+}
+
+// NSNR-CTRL-009
+TEST(NonlinearLoadBoundary, PreservesSourceOrderDuringGlobalAccumulation) {
+  const std::filesystem::path source{"nonlinear-static-control.inp"};
+  const auto fixture = MakeFixture({}, {{"10", 1, 1.0e16, {source, 30U}},
+                                        {"10", 1, -1.0e16, {source, 31U}},
+                                        {"10", 1, 1.0, {source, 32U}}});
+  const fesa::Vector trial = TrialCoordinates(fixture);
+  const auto first = fesa::AssembleGlobalFixedLoadContribution(
+      *fixture.domain, *fixture.dofs, trial, 1.0);
+  const auto second = fesa::AssembleGlobalFixedLoadContribution(
+      *fixture.domain, *fixture.dofs, trial, 1.0);
+
+  ASSERT_TRUE(first.HasValue());
+  ASSERT_TRUE(second.HasValue());
+  EXPECT_DOUBLE_EQ(first.Value().generalized_external_force[0U], 1.0);
+  EXPECT_EQ(first.Value().generalized_external_tangent.RowOffsets(),
+            second.Value().generalized_external_tangent.RowOffsets());
+  EXPECT_EQ(first.Value().generalized_external_tangent.ColumnIndices(),
+            second.Value().generalized_external_tangent.ColumnIndices());
+  EXPECT_EQ(first.Value().generalized_external_tangent.Values(),
+            second.Value().generalized_external_tangent.Values());
+}
+
+// NSNR-CTRL-010 and NSNR-CTRL-011
+TEST(NonlinearLoadBoundary,
+     DifferentiatesOneFixedSpatialMomentAndUsesSubtraction) {
+  const std::filesystem::path source{"nonlinear-static-control.inp"};
+  const auto fixture = MakeFixture({}, {{"10", 4, 3.0, {source, 40U}}});
+  fesa::Vector trial = TrialCoordinates(fixture);
+  trial[5U] = 0.4;
+  const auto result = fesa::AssembleGlobalFixedLoadContribution(
+      *fixture.domain, *fixture.dofs, trial, 0.7);
+  ASSERT_TRUE(result.HasValue());
+
+  fesa::Vector direction{fixture.dofs->FullDofCount()};
+  direction[5U] = 1.0;
+  constexpr double kStep = 1.0e-6;
+  const auto plus = fesa::AssembleGlobalFixedLoadContribution(
+      *fixture.domain, *fixture.dofs, AddScaled(trial, direction, kStep), 0.7);
+  const auto minus = fesa::AssembleGlobalFixedLoadContribution(
+      *fixture.domain, *fixture.dofs, AddScaled(trial, direction, -kStep), 0.7);
+  ASSERT_TRUE(plus.HasValue());
+  ASSERT_TRUE(minus.HasValue());
+
+  fesa::Vector derivative = plus.Value().generalized_external_force;
+  derivative.Axpy(-1.0, minus.Value().generalized_external_force);
+  derivative.Scale(0.5 / kStep);
+  const fesa::Vector tangent_action =
+      result.Value().generalized_external_tangent.Multiply(direction);
+  EXPECT_LE(MaximumAbsoluteDifference(derivative, tangent_action), 1.0e-8);
+
+  fesa::Vector internal{fixture.dofs->FullDofCount()};
+  internal[3U] = 5.0;
+  fesa::Vector residual = internal;
+  residual.Axpy(-1.0, result.Value().generalized_external_force);
+  EXPECT_DOUBLE_EQ(residual[3U],
+                   5.0 - result.Value().generalized_external_force[3U]);
+  EXPECT_DOUBLE_EQ(
+      -tangent_action[3U],
+      -result.Value().generalized_external_tangent.Multiply(direction)[3U]);
+}
+
+// NSNR-CTRL-012, NSNR-CTRL-013, and NSNR-CTRL-014
+TEST(NonlinearConvergence, UsesOwnedRelativeThresholdsForEveryChannel) {
+  static_assert(fesa::kNewtonRelativeThreshold == 1.0e-2);
+  static_assert(fesa::kNewtonZeroScaleThreshold == 1.0e-3);
+  const auto result =
+      fesa::EvaluateNewtonConvergence(RelativeConvergenceInput());
+
+  EXPECT_TRUE(result.combined_passed);
+  const std::array<fesa::NewtonChannelDecision, 4> decisions{
+      result.force_residual, result.moment_residual,
+      result.translation_correction, result.rotation_correction};
+  for (const auto& decision : decisions) {
+    EXPECT_FALSE(decision.identity.empty());
+    EXPECT_DOUBLE_EQ(decision.numerator_norm, 1.0e-2);
+    EXPECT_DOUBLE_EQ(decision.denominator_scale, 1.0);
+    EXPECT_EQ(decision.branch, fesa::NewtonConvergenceBranch::kRelative);
+    EXPECT_DOUBLE_EQ(decision.threshold, fesa::kNewtonRelativeThreshold);
+    EXPECT_EQ(decision.threshold_unit, fesa::NewtonChannelUnit::kDimensionless);
+    EXPECT_TRUE(decision.passed);
+  }
+}
+
+// NSNR-CTRL-015 and NSNR-CTRL-016
+TEST(NonlinearConvergence, AppliesTypedZeroScaleGatesToEmptyChannels) {
+  const fesa::NewtonConvergenceInput input{
+      {"residual_force", fesa::Vector{0U}, fesa::Vector{0U}, fesa::Vector{0U}},
+      {"residual_moment", fesa::Vector{0U}, fesa::Vector{0U}, fesa::Vector{0U}},
+      {"correction_translation", fesa::Vector{0U}, fesa::Vector{0U}},
+      {"correction_rotation", fesa::Vector{0U}, fesa::Vector{0U}},
+  };
+  const auto result = fesa::EvaluateNewtonConvergence(input);
+  const std::array<fesa::NewtonChannelDecision, 4> decisions{
+      result.force_residual, result.moment_residual,
+      result.translation_correction, result.rotation_correction};
+  const std::array<fesa::NewtonChannelUnit, 4> units{
+      fesa::NewtonChannelUnit::kForce, fesa::NewtonChannelUnit::kForceLength,
+      fesa::NewtonChannelUnit::kLength,
+      fesa::NewtonChannelUnit::kDimensionless};
+
+  EXPECT_TRUE(result.combined_passed);
+  for (std::size_t index = 0U; index < decisions.size(); ++index) {
+    EXPECT_DOUBLE_EQ(decisions[index].numerator_norm, 0.0);
+    EXPECT_DOUBLE_EQ(decisions[index].denominator_scale, 0.0);
+    EXPECT_EQ(decisions[index].branch,
+              fesa::NewtonConvergenceBranch::kZeroScale);
+    EXPECT_DOUBLE_EQ(decisions[index].threshold,
+                     fesa::kNewtonZeroScaleThreshold);
+    EXPECT_EQ(decisions[index].threshold_unit, units[index]);
+    EXPECT_TRUE(decisions[index].passed);
+  }
+}
+
+}  // namespace