diff --git a/include/fesa/elements/nonlinear_mitc4_shell.h b/include/fesa/elements/nonlinear_mitc4_shell.h
new file mode 100644
index 0000000..e9fcb73
--- /dev/null
+++ b/include/fesa/elements/nonlinear_mitc4_shell.h
@@ -0,0 +1,86 @@
+#ifndef FESA_ELEMENTS_NONLINEAR_MITC4_SHELL_H_
+#define FESA_ELEMENTS_NONLINEAR_MITC4_SHELL_H_
+
+#include <array>
+#include <string>
+
+#include "fesa/core/status.h"
+#include "fesa/elements/mitc4_shell.h"
+#include "fesa/math/matrix.h"
+#include "fesa/math/vector3.h"
+#include "fesa/model/model_types.h"
+#include "fesa/results/result_records.h"
+
+namespace fesa {
+
+/// @brief Stores one four-node nonlinear shell trial in global coordinates.
+struct NonlinearMitc4TrialDofs {
+  std::array<Vector3, 4> displacement;
+  std::array<Vector3, 4> rotation_vector;
+};
+
+/// @brief Stores physical shell recovery at one midsurface Gauss location.
+struct NonlinearMitc4RecoveryPoint {
+  ShellMidsurfaceLocation location;
+  std::array<double, 2> natural_coordinates;
+  Mitc4LocalFrame local_frame;
+  std::array<double, 8> generalized_strain;
+  std::array<double, 8> section_resultant;
+  std::array<ShellSectionStressRow, 3> stress;
+};
+
+/// @brief Stores the fixed four-row nonlinear physical recovery inventory.
+struct NonlinearMitc4Recovery {
+  std::array<NonlinearMitc4RecoveryPoint, 4> gauss_points;
+  double physical_strain_energy;
+};
+
+/// @brief Evaluates fixed-reference Total Lagrangian MITC4 kinematics.
+/// @note This class owns no mutable trial or accepted analysis state.
+class NonlinearMitc4Shell {
+ public:
+  /// @brief Creates one checked nonlinear MITC4 element kinematics object.
+  /// @param nodes Four non-owning source nodes in natural source order.
+  /// @param initial_frames Matching fixed initial node frames in source order.
+  /// @param section Centered positive-thickness shell section.
+  /// @param material Homogeneous isotropic material with shell-valid constants.
+  /// @return A checked evaluator or a structured model failure.
+  static Result<NonlinearMitc4Shell> Create(
+      std::array<const Node*, 4> nodes,
+      std::array<ShellNodeInitialFrame, 4> initial_frames,
+      const ShellSection& section, const LinearElasticMaterial& material);
+
+  /// @brief Evaluates trial physical kinematics and fixed recovery rows.
+  /// @param trial Global nodal translation and fixed-global rotation vectors.
+  /// @return Physical recovery or a fail-closed trial diagnostic.
+  [[nodiscard]] Result<NonlinearMitc4Recovery> EvaluateKinematics(
+      const NonlinearMitc4TrialDofs& trial) const;
+
+  /// @brief Recovers the same physical quantities for one complete trial.
+  /// @param trial Global nodal translation and fixed-global rotation vectors.
+  /// @return Four ordered Gauss rows without nodal stress projection.
+  [[nodiscard]] Result<NonlinearMitc4Recovery> RecoverState(
+      const NonlinearMitc4TrialDofs& trial) const;
+
+ private:
+  /// @brief Retains validated fixed-reference data without analysis ownership.
+  NonlinearMitc4Shell(std::array<Vector3, 4> coordinates,
+                      std::array<Matrix, 4> initial_frames,
+                      std::array<Vector3, 4> initial_directors,
+                      double thickness, double youngs_modulus,
+                      double poisson_ratio, SourceLocation source_location,
+                      std::string identity);
+
+  std::array<Vector3, 4> coordinates_;
+  std::array<Matrix, 4> initial_frames_;
+  std::array<Vector3, 4> initial_directors_;
+  double thickness_;
+  double youngs_modulus_;
+  double poisson_ratio_;
+  SourceLocation source_location_;
+  std::string identity_;
+};
+
+}  // namespace fesa
+
+#endif  // FESA_ELEMENTS_NONLINEAR_MITC4_SHELL_H_
diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt
index 725a8b1..de0913d 100644
--- a/src/fesa/CMakeLists.txt
+++ b/src/fesa/CMakeLists.txt
@@ -19,6 +19,7 @@ add_library(
     elements/euler_beam_3d.cpp
     elements/finite_rotation.cpp
     elements/mitc4_shell.cpp
+    elements/nonlinear_mitc4_shell.cpp
     fem/dof_manager.cpp
     io/abaqus/domain_builder.cpp
     io/abaqus/domain_mapping_context.cpp
diff --git a/src/fesa/elements/nonlinear_mitc4_shell.cpp b/src/fesa/elements/nonlinear_mitc4_shell.cpp
new file mode 100644
index 0000000..0fa7548
--- /dev/null
+++ b/src/fesa/elements/nonlinear_mitc4_shell.cpp
@@ -0,0 +1,614 @@
+#include "fesa/elements/nonlinear_mitc4_shell.h"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <string>
+#include <utility>
+
+#include "fesa/elements/finite_rotation.h"
+
+namespace fesa {
+namespace {
+
+constexpr std::size_t kNodeCount = 4U;
+constexpr double kFrameTolerance = 1.0e-12;
+constexpr double kDirectorOppositionTolerance = 1.0e-8;
+constexpr double kShearCorrection = 5.0 / 6.0;
+
+using Tensor = std::array<std::array<double, 3>, 3>;
+using EngineeringStrain = std::array<double, 5>;
+using EngineeringStress = std::array<double, 5>;
+
+struct Geometry {
+  std::array<Vector3, 3> reference_covariant;
+  std::array<Vector3, 3> current_covariant;
+  std::array<Vector3, 3> reciprocal;
+  Mitc4LocalFrame local_frame;
+  double reference_jacobian;
+};
+
+struct PointEvaluation {
+  Geometry geometry;
+  EngineeringStrain strain;
+  EngineeringStress stress;
+};
+
+/// @brief Builds one deterministic structured nonlinear-shell failure.
+template <class T>
+Result<T> Failure(const std::string& code, const SourceLocation& location,
+                  const std::string& identity, std::string message) {
+  return Result<T>::Failure(Status::Failure(
+      FailureCategory::kModel,
+      {{Severity::kError, code, location, "*ELEMENT", identity,
+        std::move(message)}}));
+}
+
+/// @brief Builds the stable source label concatenation used by diagnostics.
+std::string ElementIdentity(const std::array<const Node*, kNodeCount>& nodes) {
+  std::string identity;
+  for (const Node* node : nodes) {
+    if (node == nullptr) {
+      continue;
+    }
+    if (!identity.empty()) {
+      identity += "-";
+    }
+    identity += node->source_id.source_label_text;
+  }
+  return identity;
+}
+
+/// @brief Forms a weighted vector sum in source-node order.
+Vector3 WeightedSum(const std::array<double, kNodeCount>& weights,
+                    const std::array<Vector3, kNodeCount>& values) {
+  Vector3 result{};
+  for (std::size_t node = 0U; node < kNodeCount; ++node) {
+    result = result + weights[node] * values[node];
+  }
+  return result;
+}
+
+/// @brief Forms a stable derivative sum after cancelling the constant term.
+Vector3 DerivativeSum(const std::array<double, kNodeCount>& derivatives,
+                      const std::array<Vector3, kNodeCount>& values) {
+  std::array<Vector3, kNodeCount> relative{};
+  for (std::size_t node = 0U; node < kNodeCount; ++node) {
+    relative[node] = values[node] - values[0U];
+  }
+  return WeightedSum(derivatives, relative);
+}
+
+/// @brief Returns a frame matrix whose columns are a, b and the director.
+Matrix FrameMatrix(const ShellNodeInitialFrame& frame) {
+  Matrix result{3U, 3U};
+  for (std::size_t row = 0U; row < 3U; ++row) {
+    result(row, 0U) = frame.tangent_a[row];
+    result(row, 1U) = frame.tangent_b[row];
+    result(row, 2U) = frame.director[row];
+  }
+  return result;
+}
+
+/// @brief Checks every scalar in a tensor before recovery commits it.
+bool IsFinite(const Tensor& tensor) {
+  for (const auto& row : tensor) {
+    for (const double value : row) {
+      if (!std::isfinite(value)) {
+        return false;
+      }
+    }
+  }
+  return true;
+}
+
+/// @brief Checks every scalar in one fixed-size engineering quantity.
+template <std::size_t N>
+bool IsFinite(const std::array<double, N>& values) {
+  return std::all_of(values.begin(), values.end(),
+                     [](const double value) { return std::isfinite(value); });
+}
+
+/// @brief Evaluates a determinant from covariant basis columns.
+double Determinant(const std::array<Vector3, 3>& covariant) {
+  return covariant[0U].Dot(covariant[1U].Cross(covariant[2U]));
+}
+
+/// @brief Builds one proper fixed reference local frame at a midsurface point.
+bool MakeLocalFrame(const Vector3& midsurface_xi, const Vector3& midsurface_eta,
+                    const Vector3& director, Mitc4LocalFrame& frame) {
+  const Vector3 area = midsurface_xi.Cross(midsurface_eta);
+  const double director_norm = director.Norm();
+  if (!area.IsFinite() || !director.IsFinite() ||
+      !(director_norm > kDirectorOppositionTolerance) ||
+      !std::isfinite(director_norm)) {
+    return false;
+  }
+
+  const Vector3 e3 = director / director_norm;
+  if (!(area.Dot(e3) > 0.0)) {
+    return false;
+  }
+  const Vector3 e1_candidate = midsurface_xi - midsurface_xi.Dot(e3) * e3;
+  const double e1_norm = e1_candidate.Norm();
+  if (!e1_candidate.IsFinite() || !std::isfinite(e1_norm) ||
+      !(e1_norm > 0.0)) {
+    return false;
+  }
+  const Vector3 e1 = e1_candidate / e1_norm;
+  const Vector3 e2 = e3.Cross(e1);
+  if (!e2.IsFinite()) {
+    return false;
+  }
+  frame = {e1.Components(), e2.Components(), e3.Components()};
+  return true;
+}
+
+/// @brief Evaluates fixed-reference and trial covariant geometry together.
+bool EvaluateGeometry(const std::array<Vector3, kNodeCount>& coordinates,
+                      const std::array<Vector3, kNodeCount>& initial_directors,
+                      const std::array<Vector3, kNodeCount>& directors,
+                      const NonlinearMitc4TrialDofs& trial,
+                      const double thickness, const double xi, const double eta,
+                      const double zeta, Geometry& geometry) {
+  const Mitc4ShapeFunctions shape = Mitc4Shell::ShapeFunctions(xi, eta);
+  const Vector3 reference_xi =
+      DerivativeSum(shape.xi_derivatives, coordinates);
+  const Vector3 reference_eta =
+      DerivativeSum(shape.eta_derivatives, coordinates);
+  const Vector3 initial_director_xi =
+      DerivativeSum(shape.xi_derivatives, initial_directors);
+  const Vector3 initial_director_eta =
+      DerivativeSum(shape.eta_derivatives, initial_directors);
+  const Vector3 initial_director = WeightedSum(shape.values, initial_directors);
+  const Vector3 displacement_xi =
+      DerivativeSum(shape.xi_derivatives, trial.displacement);
+  const Vector3 displacement_eta =
+      DerivativeSum(shape.eta_derivatives, trial.displacement);
+  const Vector3 director_xi = DerivativeSum(shape.xi_derivatives, directors);
+  const Vector3 director_eta = DerivativeSum(shape.eta_derivatives, directors);
+  const Vector3 director = WeightedSum(shape.values, directors);
+  const double half_thickness = 0.5 * thickness;
+
+  geometry.reference_covariant[0U] =
+      reference_xi + half_thickness * zeta * initial_director_xi;
+  geometry.reference_covariant[1U] =
+      reference_eta + half_thickness * zeta * initial_director_eta;
+  geometry.reference_covariant[2U] = half_thickness * initial_director;
+  geometry.current_covariant[0U] =
+      reference_xi + displacement_xi + half_thickness * zeta * director_xi;
+  geometry.current_covariant[1U] =
+      reference_eta + displacement_eta + half_thickness * zeta * director_eta;
+  geometry.current_covariant[2U] = half_thickness * director;
+  geometry.reference_jacobian = Determinant(geometry.reference_covariant);
+  const double current_jacobian = Determinant(geometry.current_covariant);
+
+  if (!reference_xi.IsFinite() || !reference_eta.IsFinite() ||
+      !initial_director.IsFinite() || !director.IsFinite() ||
+      !(director.Norm() > kDirectorOppositionTolerance) ||
+      !std::isfinite(geometry.reference_jacobian) ||
+      !(geometry.reference_jacobian > 0.0) ||
+      !std::isfinite(current_jacobian) || !(current_jacobian > 0.0)) {
+    return false;
+  }
+
+  geometry.reciprocal[0U] =
+      geometry.reference_covariant[1U].Cross(geometry.reference_covariant[2U]) /
+      geometry.reference_jacobian;
+  geometry.reciprocal[1U] =
+      geometry.reference_covariant[2U].Cross(geometry.reference_covariant[0U]) /
+      geometry.reference_jacobian;
+  geometry.reciprocal[2U] =
+      geometry.reference_covariant[0U].Cross(geometry.reference_covariant[1U]) /
+      geometry.reference_jacobian;
+  if (!geometry.reciprocal[0U].IsFinite() ||
+      !geometry.reciprocal[1U].IsFinite() ||
+      !geometry.reciprocal[2U].IsFinite()) {
+    return false;
+  }
+  return MakeLocalFrame(reference_xi, reference_eta, initial_director,
+                        geometry.local_frame);
+}
+
+/// @brief Forms the direct Green-Lagrange covariant strain with E33 excluded.
+Tensor DirectStrain(const Geometry& geometry) {
+  Tensor strain{};
+  for (std::size_t row = 0U; row < 3U; ++row) {
+    for (std::size_t column = 0U; column < 3U; ++column) {
+      strain[row][column] =
+          0.5 * (geometry.current_covariant[row].Dot(
+                     geometry.current_covariant[column]) -
+                 geometry.reference_covariant[row].Dot(
+                     geometry.reference_covariant[column]));
+    }
+  }
+  strain[2U][2U] = 0.0;
+  return strain;
+}
+
+/// @brief Reconstructs a fixed-reference Cartesian strain from covariant data.
+Tensor CartesianStrain(const Tensor& covariant,
+                       const std::array<Vector3, 3>& reciprocal) {
+  Tensor cartesian{};
+  for (std::size_t row = 0U; row < 3U; ++row) {
+    for (std::size_t column = 0U; column < 3U; ++column) {
+      for (std::size_t first = 0U; first < 3U; ++first) {
+        for (std::size_t second = 0U; second < 3U; ++second) {
+          cartesian[row][column] += covariant[first][second] *
+                                    reciprocal[first][row] *
+                                    reciprocal[second][column];
+        }
+      }
+    }
+  }
+  return cartesian;
+}
+
+/// @brief Evaluates a tensor component in a fixed local reference frame.
+double FrameComponent(const Vector3& left, const Tensor& tensor,
+                      const Vector3& right) {
+  double value = 0.0;
+  for (std::size_t row = 0U; row < 3U; ++row) {
+    for (std::size_t column = 0U; column < 3U; ++column) {
+      value += left[row] * tensor[row][column] * right[column];
+    }
+  }
+  return value;
+}
+
+/// @brief Converts fixed-reference Cartesian strain to engineering components.
+EngineeringStrain LocalEngineering(const Tensor& tensor,
+                                   const Mitc4LocalFrame& frame) {
+  const Vector3 e1{frame.e1};
+  const Vector3 e2{frame.e2};
+  const Vector3 e3{frame.e3};
+  return {FrameComponent(e1, tensor, e1), FrameComponent(e2, tensor, e2),
+          2.0 * FrameComponent(e1, tensor, e2),
+          2.0 * FrameComponent(e1, tensor, e3),
+          2.0 * FrameComponent(e2, tensor, e3)};
+}
+
+/// @brief Applies the approved homogeneous plane-stress plus shear law.
+EngineeringStress Stress(const EngineeringStrain& strain, const double youngs,
+                         const double poisson) {
+  const double factor = youngs / (1.0 - poisson * poisson);
+  const double shear = kShearCorrection * youngs / (2.0 * (1.0 + poisson));
+  return {factor * (strain[0U] + poisson * strain[1U]),
+          factor * (poisson * strain[0U] + strain[1U]),
+          factor * (1.0 - poisson) * 0.5 * strain[2U],
+          shear * strain[3U], shear * strain[4U]};
+}
+
+/// @brief Evaluates direct strain at one mandatory edge-midpoint tying point.
+bool EvaluateTyingSample(const std::array<Vector3, kNodeCount>& coordinates,
+                         const std::array<Vector3, kNodeCount>& initial_directors,
+                         const std::array<Vector3, kNodeCount>& directors,
+                         const NonlinearMitc4TrialDofs& trial,
+                         const double thickness, const double xi,
+                         const double eta, Tensor& strain) {
+  Geometry geometry{};
+  if (!EvaluateGeometry(coordinates, initial_directors, directors, trial,
+                        thickness, xi, eta, 0.0, geometry)) {
+    return false;
+  }
+  strain = DirectStrain(geometry);
+  return IsFinite(strain);
+}
+
+/// @brief Computes one MITC-projected local strain and material stress point.
+bool EvaluatePoint(const std::array<Vector3, kNodeCount>& coordinates,
+                   const std::array<Vector3, kNodeCount>& initial_directors,
+                   const std::array<Vector3, kNodeCount>& directors,
+                   const NonlinearMitc4TrialDofs& trial,
+                   const std::array<Tensor, 4>& tying_samples,
+                   const double thickness, const double youngs,
+                   const double poisson, const double xi, const double eta,
+                   const double zeta, PointEvaluation& point) {
+  if (!EvaluateGeometry(coordinates, initial_directors, directors, trial,
+                        thickness, xi, eta, zeta, point.geometry)) {
+    return false;
+  }
+  Tensor covariant = DirectStrain(point.geometry);
+  const Mitc4TyingWeights weights = Mitc4Shell::TyingWeights(xi, eta);
+  covariant[0U][2U] = weights.xi_zeta[0U] * tying_samples[0U][0U][2U] +
+                       weights.xi_zeta[1U] * tying_samples[1U][0U][2U];
+  covariant[2U][0U] = covariant[0U][2U];
+  covariant[1U][2U] = weights.eta_zeta[0U] * tying_samples[2U][1U][2U] +
+                       weights.eta_zeta[1U] * tying_samples[3U][1U][2U];
+  covariant[2U][1U] = covariant[1U][2U];
+  const Tensor cartesian = CartesianStrain(covariant, point.geometry.reciprocal);
+  point.strain = LocalEngineering(cartesian, point.geometry.local_frame);
+  point.stress = Stress(point.strain, youngs, poisson);
+  return IsFinite(covariant) && IsFinite(cartesian) && IsFinite(point.strain) &&
+         IsFinite(point.stress);
+}
+
+/// @brief Tests the fixed-order physical recovery inventory for finite values.
+bool IsFinite(const NonlinearMitc4Recovery& recovery) {
+  if (!std::isfinite(recovery.physical_strain_energy)) {
+    return false;
+  }
+  for (const auto& point : recovery.gauss_points) {
+    if (!IsFinite(point.generalized_strain) ||
+        !IsFinite(point.section_resultant)) {
+      return false;
+    }
+    for (const auto& stress : point.stress) {
+      if (!IsFinite(stress.components)) {
+        return false;
+      }
+    }
+  }
+  return true;
+}
+
+}  // namespace
+
+NonlinearMitc4Shell::NonlinearMitc4Shell(
+    std::array<Vector3, 4> coordinates, std::array<Matrix, 4> initial_frames,
+    std::array<Vector3, 4> initial_directors, const double thickness,
+    const double youngs_modulus, const double poisson_ratio,
+    SourceLocation source_location, std::string identity)
+    : coordinates_{std::move(coordinates)},
+      initial_frames_{std::move(initial_frames)},
+      initial_directors_{std::move(initial_directors)},
+      thickness_{thickness},
+      youngs_modulus_{youngs_modulus},
+      poisson_ratio_{poisson_ratio},
+      source_location_{std::move(source_location)},
+      identity_{std::move(identity)} {}
+
+Result<NonlinearMitc4Shell> NonlinearMitc4Shell::Create(
+    std::array<const Node*, 4> nodes,
+    std::array<ShellNodeInitialFrame, 4> initial_frames,
+    const ShellSection& section, const LinearElasticMaterial& material) {
+  const std::string identity = ElementIdentity(nodes);
+  if (std::any_of(nodes.begin(), nodes.end(),
+                  [](const Node* node) { return node == nullptr; })) {
+    return Failure<NonlinearMitc4Shell>(
+        "invalid-nonlinear-shell-geometry", section.location, identity,
+        "Nonlinear MITC4 creation requires four valid node references.");
+  }
+
+  std::array<Vector3, kNodeCount> coordinates{};
+  for (std::size_t node = 0U; node < kNodeCount; ++node) {
+    coordinates[node] = Vector3{nodes[node]->coordinates};
+    if (!coordinates[node].IsFinite()) {
+      return Failure<NonlinearMitc4Shell>(
+          "invalid-nonlinear-shell-geometry", nodes[node]->location, identity,
+          "Nonlinear MITC4 node coordinates must be finite.");
+    }
+    for (std::size_t previous = 0U; previous < node; ++previous) {
+      if (coordinates[node] == coordinates[previous]) {
+        return Failure<NonlinearMitc4Shell>(
+            "invalid-nonlinear-shell-geometry", nodes[node]->location,
+            identity, "Nonlinear MITC4 nodes must have distinct coordinates.");
+      }
+      if (initial_frames[node].node_index == initial_frames[previous].node_index) {
+        return Failure<NonlinearMitc4Shell>(
+            "invalid-nonlinear-shell-frame", nodes[node]->location, identity,
+            "Nonlinear MITC4 initial frames must identify four distinct nodes.");
+      }
+    }
+  }
+
+  std::array<Matrix, kNodeCount> frames{Matrix{3U, 3U}, Matrix{3U, 3U},
+                                        Matrix{3U, 3U}, Matrix{3U, 3U}};
+  std::array<Vector3, kNodeCount> directors{};
+  for (std::size_t node = 0U; node < kNodeCount; ++node) {
+    frames[node] = FrameMatrix(initial_frames[node]);
+    const auto state = MakeFiniteRotationState(Vector3{}, frames[node]);
+    if (!state.HasValue()) {
+      return Failure<NonlinearMitc4Shell>(
+          "invalid-nonlinear-shell-frame", nodes[node]->location, identity,
+          "Nonlinear MITC4 initial frames must be finite proper triads.");
+    }
+    directors[node] = state.Value().director;
+    if ((directors[node] - Vector3{initial_frames[node].director}).Norm() >
+        kFrameTolerance) {
+      return Failure<NonlinearMitc4Shell>(
+          "invalid-nonlinear-shell-frame", nodes[node]->location, identity,
+          "Nonlinear MITC4 initial frame director must be its third axis.");
+    }
+  }
+
+  const Mitc4ShapeFunctions center = Mitc4Shell::ShapeFunctions(0.0, 0.0);
+  const Vector3 center_xi = DerivativeSum(center.xi_derivatives, coordinates);
+  const Vector3 center_eta = DerivativeSum(center.eta_derivatives, coordinates);
+  const Vector3 center_area = center_xi.Cross(center_eta);
+  const double center_measure = center_area.Norm();
+  if (!center_area.IsFinite() || !std::isfinite(center_measure) ||
+      !(center_measure > 0.0)) {
+    return Failure<NonlinearMitc4Shell>(
+        "invalid-nonlinear-shell-geometry", nodes[0U]->location, identity,
+        "Nonlinear MITC4 center surface basis must be finite and nonzero.");
+  }
+  const Vector3 positive_face = center_area / center_measure;
+  if (std::any_of(directors.begin(), directors.end(),
+                  [&positive_face](const Vector3& director) {
+                    return !(positive_face.Dot(director) > 0.0);
+                  })) {
+    return Failure<NonlinearMitc4Shell>(
+        "invalid-nonlinear-shell-frame", nodes[0U]->location, identity,
+        "Nonlinear MITC4 initial directors must follow the source positive face.");
+  }
+
+  if (!std::isfinite(section.thickness) || !(section.thickness > 0.0)) {
+    return Failure<NonlinearMitc4Shell>(
+        "invalid-nonlinear-shell-section", section.location, identity,
+        "Nonlinear MITC4 shell thickness must be finite and positive.");
+  }
+  if (!std::isfinite(material.youngs_modulus) ||
+      !(material.youngs_modulus > 0.0) ||
+      !std::isfinite(material.poisson_ratio) ||
+      !(material.poisson_ratio > -1.0) || !(material.poisson_ratio < 0.5)) {
+    return Failure<NonlinearMitc4Shell>(
+        "invalid-nonlinear-shell-material", material.location, identity,
+        "Nonlinear MITC4 requires finite E>0 and -1<nu<0.5.");
+  }
+
+  NonlinearMitc4Shell shell{coordinates,
+                             std::move(frames),
+                             directors,
+                             section.thickness,
+                             material.youngs_modulus,
+                             material.poisson_ratio,
+                             nodes[0U]->location,
+                             identity};
+  const auto reference = shell.EvaluateKinematics(NonlinearMitc4TrialDofs{});
+  if (!reference.HasValue()) {
+    return Result<NonlinearMitc4Shell>::Failure(reference.GetStatus());
+  }
+  return Result<NonlinearMitc4Shell>::Success(std::move(shell));
+}
+
+Result<NonlinearMitc4Recovery> NonlinearMitc4Shell::EvaluateKinematics(
+    const NonlinearMitc4TrialDofs& trial) const {
+  std::array<Vector3, kNodeCount> directors{};
+  for (std::size_t node = 0U; node < kNodeCount; ++node) {
+    if (!trial.displacement[node].IsFinite() ||
+        !trial.rotation_vector[node].IsFinite()) {
+      return Failure<NonlinearMitc4Recovery>(
+          "invalid-nonlinear-shell-trial", source_location_, identity_,
+          "Nonlinear MITC4 trial displacements and rotations must be finite.");
+    }
+    const auto state =
+        MakeFiniteRotationState(trial.rotation_vector[node], initial_frames_[node]);
+    if (!state.HasValue()) {
+      return Failure<NonlinearMitc4Recovery>(
+          "invalid-nonlinear-shell-rotation", source_location_, identity_,
+          "Nonlinear MITC4 trial rotations must remain in the principal chart.");
+    }
+    directors[node] = state.Value().director;
+  }
+  for (std::size_t node = 0U; node < kNodeCount; ++node) {
+    for (std::size_t other = 0U; other < node; ++other) {
+      if (!(directors[node].Dot(directors[other]) >
+            -1.0 + kDirectorOppositionTolerance)) {
+        return Failure<NonlinearMitc4Recovery>(
+            "invalid-nonlinear-shell-director", source_location_, identity_,
+            "Nonlinear MITC4 trial directors must not be opposed.");
+      }
+    }
+  }
+
+  Geometry center{};
+  if (!EvaluateGeometry(coordinates_, initial_directors_, directors, trial,
+                        thickness_, 0.0, 0.0, 0.0, center)) {
+    return Failure<NonlinearMitc4Recovery>(
+        "invalid-nonlinear-shell-geometry", source_location_, identity_,
+        "Nonlinear MITC4 center reference or trial geometry is invalid.");
+  }
+
+  constexpr std::array<std::array<double, 2>, 4> tying_locations{
+      std::array<double, 2>{0.0, -1.0}, std::array<double, 2>{0.0, 1.0},
+      std::array<double, 2>{-1.0, 0.0}, std::array<double, 2>{1.0, 0.0}};
+  std::array<Tensor, 4> tying_samples{};
+  for (std::size_t point = 0U; point < tying_locations.size(); ++point) {
+    if (!EvaluateTyingSample(coordinates_, initial_directors_, directors, trial,
+                             thickness_, tying_locations[point][0U],
+                             tying_locations[point][1U], tying_samples[point])) {
+      return Failure<NonlinearMitc4Recovery>(
+          "invalid-nonlinear-shell-geometry", source_location_, identity_,
+          "Nonlinear MITC4 tying-point reference or trial geometry is invalid.");
+    }
+  }
+
+  const double gauss = 1.0 / std::sqrt(3.0);
+  constexpr std::array<ShellMidsurfaceLocation, 4> locations{
+      ShellMidsurfaceLocation::kGp1, ShellMidsurfaceLocation::kGp2,
+      ShellMidsurfaceLocation::kGp3, ShellMidsurfaceLocation::kGp4};
+  const std::array<std::array<double, 2>, 4> surface_points{
+      std::array<double, 2>{-gauss, -gauss},
+      std::array<double, 2>{gauss, -gauss},
+      std::array<double, 2>{gauss, gauss},
+      std::array<double, 2>{-gauss, gauss}};
+  constexpr std::array<double, 2> thickness_gauss{-1.0, 1.0};
+  constexpr std::array<double, 3> stress_zeta{-1.0, 0.0, 1.0};
+  constexpr std::array<ShellSectionPosition, 3> stress_position{
+      ShellSectionPosition::kBottom, ShellSectionPosition::kMiddle,
+      ShellSectionPosition::kTop};
+
+  NonlinearMitc4Recovery recovery{};
+  for (std::size_t surface = 0U; surface < surface_points.size(); ++surface) {
+    auto& row = recovery.gauss_points[surface];
+    row.location = locations[surface];
+    row.natural_coordinates = surface_points[surface];
+    Geometry midpoint{};
+    if (!EvaluateGeometry(coordinates_, initial_directors_, directors, trial,
+                          thickness_, row.natural_coordinates[0U],
+                          row.natural_coordinates[1U], 0.0, midpoint)) {
+      return Failure<NonlinearMitc4Recovery>(
+          "invalid-nonlinear-shell-geometry", source_location_, identity_,
+          "Nonlinear MITC4 midsurface recovery geometry is invalid.");
+    }
+    row.local_frame = midpoint.local_frame;
+
+    for (const double thickness_sign : thickness_gauss) {
+      PointEvaluation point{};
+      if (!EvaluatePoint(coordinates_, initial_directors_, directors, trial,
+                         tying_samples, thickness_, youngs_modulus_,
+                         poisson_ratio_, row.natural_coordinates[0U],
+                         row.natural_coordinates[1U], thickness_sign * gauss,
+                         point)) {
+        return Failure<NonlinearMitc4Recovery>(
+            "invalid-nonlinear-shell-geometry", source_location_, identity_,
+            "Nonlinear MITC4 volume recovery geometry is invalid.");
+      }
+      for (std::size_t component = 0U; component < 3U; ++component) {
+        row.generalized_strain[component] += 0.5 * point.strain[component];
+        row.generalized_strain[3U + component] +=
+            3.0 * thickness_sign * gauss * point.strain[component] /
+            thickness_;
+        row.section_resultant[component] +=
+            0.5 * thickness_ * point.stress[component];
+        row.section_resultant[3U + component] +=
+            0.25 * thickness_ * thickness_sign * gauss * point.stress[component];
+      }
+      for (std::size_t component = 0U; component < 2U; ++component) {
+        row.generalized_strain[6U + component] +=
+            0.5 * point.strain[3U + component];
+        row.section_resultant[6U + component] +=
+            0.5 * thickness_ * point.stress[3U + component];
+      }
+      double density = 0.0;
+      for (std::size_t component = 0U; component < point.strain.size();
+           ++component) {
+        density += point.strain[component] * point.stress[component];
+      }
+      recovery.physical_strain_energy +=
+          0.5 * density * point.geometry.reference_jacobian;
+    }
+
+    for (std::size_t position = 0U; position < stress_zeta.size(); ++position) {
+      PointEvaluation point{};
+      if (!EvaluatePoint(coordinates_, initial_directors_, directors, trial,
+                         tying_samples, thickness_, youngs_modulus_,
+                         poisson_ratio_, row.natural_coordinates[0U],
+                         row.natural_coordinates[1U], stress_zeta[position],
+                         point)) {
+        return Failure<NonlinearMitc4Recovery>(
+            "invalid-nonlinear-shell-geometry", source_location_, identity_,
+            "Nonlinear MITC4 section-position recovery geometry is invalid.");
+      }
+      row.stress[position] = {stress_position[position], stress_zeta[position],
+                              {point.stress[0U], point.stress[1U],
+                               point.stress[2U]}};
+    }
+  }
+
+  if (!IsFinite(recovery)) {
+    return Failure<NonlinearMitc4Recovery>(
+        "invalid-nonlinear-shell-recovery", source_location_, identity_,
+        "Nonlinear MITC4 physical recovery must remain finite.");
+  }
+  return Result<NonlinearMitc4Recovery>::Success(std::move(recovery));
+}
+
+Result<NonlinearMitc4Recovery> NonlinearMitc4Shell::RecoverState(
+    const NonlinearMitc4TrialDofs& trial) const {
+  return EvaluateKinematics(trial);
+}
+
+}  // namespace fesa
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index e7444f5..69d92e9 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -20,6 +20,7 @@ add_executable(
     unit/elements/euler_beam_3d_test.cpp
     unit/elements/finite_rotation_test.cpp
     unit/elements/mitc4_shell_test.cpp
+    unit/elements/nonlinear_mitc4_shell_test.cpp
     unit/fem/dof_manager_test.cpp
     unit/math/dense_blas_internal_test.cpp
     unit/math/matrix_test.cpp
diff --git a/tests/unit/elements/nonlinear_mitc4_shell_test.cpp b/tests/unit/elements/nonlinear_mitc4_shell_test.cpp
new file mode 100644
index 0000000..f293d18
--- /dev/null
+++ b/tests/unit/elements/nonlinear_mitc4_shell_test.cpp
@@ -0,0 +1,263 @@
+#include "fesa/elements/nonlinear_mitc4_shell.h"
+
+#include <gtest/gtest.h>
+
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <string>
+
+#include "fesa/elements/finite_rotation.h"
+#include "fesa/elements/mitc4_shell.h"
+#include "fesa/model/model_types.h"
+
+namespace {
+
+using Vector3 = std::array<double, 3>;
+
+fesa::Node Node(const std::int64_t label, const Vector3 coordinates) {
+  return {{"Shell-Instance", label, std::to_string(label)}, coordinates,
+          {"nonlinear-mitc4.inp", static_cast<std::size_t>(label)}};
+}
+
+std::array<fesa::Node, 4> Nodes() {
+  return {Node(1, {-1.0, -1.0, 0.0}), Node(2, {1.0, -1.0, 0.0}),
+          Node(3, {1.0, 1.0, 0.0}), Node(4, {-1.0, 1.0, 0.0})};
+}
+
+std::array<const fesa::Node*, 4> NodePointers(
+    const std::array<fesa::Node, 4>& nodes) {
+  return {&nodes[0], &nodes[1], &nodes[2], &nodes[3]};
+}
+
+std::array<fesa::ShellNodeInitialFrame, 4> Frames() {
+  return {{{0U, {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}},
+           {1U, {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}},
+           {2U, {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}},
+           {3U, {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}}}};
+}
+
+fesa::ShellSection Section() {
+  return {"Section", 0.2, 0U, {"nonlinear-mitc4.inp", 20U}};
+}
+
+fesa::LinearElasticMaterial Material() {
+  return {"Material", 120.0, 0.25, {"nonlinear-mitc4.inp", 21U}};
+}
+
+fesa::NonlinearMitc4TrialDofs ZeroTrial() { return {}; }
+
+fesa::Vector3 Multiply(const fesa::Matrix& matrix, const fesa::Vector3& value) {
+  return {matrix(0U, 0U) * value.X() + matrix(0U, 1U) * value.Y() +
+              matrix(0U, 2U) * value.Z(),
+          matrix(1U, 0U) * value.X() + matrix(1U, 1U) * value.Y() +
+              matrix(1U, 2U) * value.Z(),
+          matrix(2U, 0U) * value.X() + matrix(2U, 1U) * value.Y() +
+              matrix(2U, 2U) * value.Z()};
+}
+
+void ExpectZeroRecovery(const fesa::NonlinearMitc4Recovery& recovery) {
+  EXPECT_DOUBLE_EQ(recovery.physical_strain_energy, 0.0);
+  for (const auto& point : recovery.gauss_points) {
+    for (const double value : point.generalized_strain) {
+      EXPECT_DOUBLE_EQ(value, 0.0);
+    }
+    for (const double value : point.section_resultant) {
+      EXPECT_DOUBLE_EQ(value, 0.0);
+    }
+    for (const auto& stress : point.stress) {
+      for (const double value : stress.components) {
+        EXPECT_DOUBLE_EQ(value, 0.0);
+      }
+    }
+  }
+}
+
+fesa::NonlinearMitc4Shell Shell(const std::array<fesa::Node, 4>& nodes) {
+  const auto result = fesa::NonlinearMitc4Shell::Create(
+      NodePointers(nodes), Frames(), Section(), Material());
+  EXPECT_TRUE(result.HasValue());
+  return result.Value();
+}
+
+// NSNR-ELEM-KIN-001
+TEST(NonlinearMitc4Shell, UndeformedStateHasZeroPhysicalRecovery) {
+  const auto nodes = Nodes();
+  const auto recovery = Shell(nodes).RecoverState(ZeroTrial());
+
+  ASSERT_TRUE(recovery.HasValue());
+  ExpectZeroRecovery(recovery.Value());
+}
+
+// NSNR-ELEM-KIN-001A
+TEST(NonlinearMitc4Shell, RejectsInvalidNodeAndInitialFrameContracts) {
+  const auto nodes = Nodes();
+  auto node_pointers = NodePointers(nodes);
+  node_pointers[3U] = nullptr;
+  const auto null_node = fesa::NonlinearMitc4Shell::Create(
+      node_pointers, Frames(), Section(), Material());
+  ASSERT_FALSE(null_node.HasValue());
+  EXPECT_EQ(null_node.GetStatus().Diagnostics()[0U].code,
+            "invalid-nonlinear-shell-geometry");
+
+  auto frames = Frames();
+  frames[2U].node_index = frames[1U].node_index;
+  const auto duplicate_frame = fesa::NonlinearMitc4Shell::Create(
+      NodePointers(nodes), frames, Section(), Material());
+  ASSERT_FALSE(duplicate_frame.HasValue());
+  EXPECT_EQ(duplicate_frame.GetStatus().Diagnostics()[0U].code,
+            "invalid-nonlinear-shell-frame");
+
+  frames = Frames();
+  frames[0U].tangent_b = {1.0, 0.0, 0.0};
+  const auto invalid_frame = fesa::NonlinearMitc4Shell::Create(
+      NodePointers(nodes), frames, Section(), Material());
+  ASSERT_FALSE(invalid_frame.HasValue());
+  EXPECT_EQ(invalid_frame.GetStatus().Diagnostics()[0U].code,
+            "invalid-nonlinear-shell-frame");
+}
+
+// NSNR-ELEM-KIN-002
+TEST(NonlinearMitc4Shell, UsesFixedInitialReferenceGeometry) {
+  const auto nodes = Nodes();
+  auto trial = ZeroTrial();
+  for (auto& displacement : trial.displacement) {
+    displacement = {0.2, -0.3, 0.4};
+  }
+
+  const auto recovery = Shell(nodes).EvaluateKinematics(trial);
+  ASSERT_TRUE(recovery.HasValue());
+  ExpectZeroRecovery(recovery.Value());
+}
+
+// NSNR-ELEM-KIN-003
+TEST(NonlinearMitc4Shell, PreservesUnnormalizedInterpolatedDirector) {
+  const auto nodes = Nodes();
+  auto trial = ZeroTrial();
+  trial.rotation_vector[0U] = {0.8, 0.0, 0.0};
+  trial.rotation_vector[1U] = {-0.8, 0.0, 0.0};
+  trial.rotation_vector[2U] = {-0.8, 0.0, 0.0};
+  trial.rotation_vector[3U] = {0.8, 0.0, 0.0};
+
+  const auto recovery = Shell(nodes).EvaluateKinematics(trial);
+  ASSERT_TRUE(recovery.HasValue());
+  EXPECT_GT(recovery.Value().physical_strain_energy, 0.0);
+  const double expected_membrane =
+      0.01 * std::sin(0.8) * std::sin(0.8) / 6.0;
+  EXPECT_NEAR(recovery.Value().gauss_points[0U].generalized_strain[0U],
+              expected_membrane, 1.0e-12);
+}
+
+// NSNR-ELEM-KIN-004
+TEST(NonlinearMitc4Shell, RejectsOpposedOrDegenerateCurrentDirectors) {
+  const auto nodes = Nodes();
+  auto trial = ZeroTrial();
+  trial.rotation_vector[0U] = {std::acos(-1.0) - 1.0e-5, 0.0, 0.0};
+  trial.rotation_vector[1U] = {-(std::acos(-1.0) - 1.0e-5), 0.0, 0.0};
+
+  const auto recovery = Shell(nodes).EvaluateKinematics(trial);
+  ASSERT_FALSE(recovery.HasValue());
+  EXPECT_EQ(recovery.GetStatus().Diagnostics()[0U].code,
+            "invalid-nonlinear-shell-director");
+}
+
+// NSNR-ELEM-KIN-005
+TEST(NonlinearMitc4Shell, MatchesLinearMitc4RecoveryInSmallRotationLimit) {
+  const auto nodes = Nodes();
+  const auto linear = fesa::Mitc4Shell::Create(
+      NodePointers(nodes),
+      std::array<Vector3, 4>{{{0.0, 0.0, 1.0}, {0.0, 0.0, 1.0},
+                               {0.0, 0.0, 1.0}, {0.0, 0.0, 1.0}}},
+      Section(), Material());
+  ASSERT_TRUE(linear.HasValue());
+
+  auto trial = ZeroTrial();
+  fesa::Vector linear_dofs{24U};
+  for (std::size_t node = 0U; node < 4U; ++node) {
+    const double x = nodes[node].coordinates[0U];
+    const double y = nodes[node].coordinates[1U];
+    trial.displacement[node] = {1.0e-7 * x, -2.0e-7 * y, 3.0e-7 * x * y};
+    trial.rotation_vector[node] = {-2.0e-7 * y, 1.0e-7 * x, 0.0};
+    for (std::size_t component = 0U; component < 3U; ++component) {
+      linear_dofs[6U * node + component] = trial.displacement[node][component];
+      linear_dofs[6U * node + 3U + component] =
+          trial.rotation_vector[node][component];
+    }
+  }
+
+  const auto nonlinear_recovery = Shell(nodes).EvaluateKinematics(trial);
+  const auto linear_recovery = linear.Value().RecoverPhysical(linear_dofs);
+  ASSERT_TRUE(nonlinear_recovery.HasValue());
+  ASSERT_TRUE(linear_recovery.HasValue());
+  for (std::size_t point = 0U; point < 4U; ++point) {
+    for (std::size_t component = 0U; component < 8U; ++component) {
+      EXPECT_NEAR(nonlinear_recovery.Value().gauss_points[point]
+                      .generalized_strain[component],
+                  linear_recovery.Value().points[point]
+                      .generalized_strain[component],
+                  1.0e-11);
+    }
+  }
+}
+
+// NSNR-ELEM-KIN-006
+TEST(NonlinearMitc4Shell, FiniteRigidMotionIsObjective) {
+  const auto nodes = Nodes();
+  const fesa::Vector3 rotation{0.3, -0.2, 0.1};
+  const fesa::Matrix rigid_rotation = fesa::ExpSkew(rotation);
+  const fesa::Vector3 translation{0.7, -0.4, 1.2};
+  auto trial = ZeroTrial();
+  for (std::size_t node = 0U; node < 4U; ++node) {
+    const fesa::Vector3 initial{nodes[node].coordinates};
+    trial.displacement[node] =
+        Multiply(rigid_rotation, initial) + translation - initial;
+    trial.rotation_vector[node] = rotation;
+  }
+
+  const auto recovery = Shell(nodes).EvaluateKinematics(trial);
+  ASSERT_TRUE(recovery.HasValue());
+  for (const auto& point : recovery.Value().gauss_points) {
+    for (const double value : point.generalized_strain) {
+      EXPECT_NEAR(value, 0.0, 1.0e-12);
+    }
+    for (const auto& stress : point.stress) {
+      for (const double value : stress.components) {
+        EXPECT_NEAR(value, 0.0, 1.0e-12);
+      }
+    }
+  }
+  EXPECT_NEAR(recovery.Value().physical_strain_energy, 0.0, 1.0e-12);
+}
+
+// NSNR-ELEM-KIN-007 and NSNR-ELEM-KIN-008
+TEST(NonlinearMitc4Recovery, EmitsOnlyFixedPhysicalGaussAndStressRows) {
+  const auto nodes = Nodes();
+  const auto recovery = Shell(nodes).RecoverState(ZeroTrial());
+  ASSERT_TRUE(recovery.HasValue());
+  const double gauss = 1.0 / std::sqrt(3.0);
+  constexpr std::array<fesa::ShellMidsurfaceLocation, 4> locations{
+      fesa::ShellMidsurfaceLocation::kGp1, fesa::ShellMidsurfaceLocation::kGp2,
+      fesa::ShellMidsurfaceLocation::kGp3, fesa::ShellMidsurfaceLocation::kGp4};
+  constexpr std::array<fesa::ShellSectionPosition, 3> positions{
+      fesa::ShellSectionPosition::kBottom, fesa::ShellSectionPosition::kMiddle,
+      fesa::ShellSectionPosition::kTop};
+  constexpr std::array<double, 3> zeta{-1.0, 0.0, 1.0};
+
+  for (std::size_t point = 0U; point < recovery.Value().gauss_points.size();
+       ++point) {
+    const auto& row = recovery.Value().gauss_points[point];
+    EXPECT_EQ(row.location, locations[point]);
+    EXPECT_EQ(row.natural_coordinates[0U], point == 0U || point == 3U ? -gauss
+                                                                       : gauss);
+    EXPECT_EQ(row.natural_coordinates[1U], point < 2U ? -gauss : gauss);
+    for (std::size_t stress = 0U; stress < row.stress.size(); ++stress) {
+      EXPECT_EQ(row.stress[stress].position, positions[stress]);
+      EXPECT_EQ(row.stress[stress].zeta, zeta[stress]);
+      EXPECT_EQ(row.stress[stress].components.size(), 3U);
+    }
+  }
+}
+
+}  // namespace