#include <gtest/gtest.h>

#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <cstring>
#include <stdexcept>
#include <string>

#include "fesa/elements/finite_rotation.h"
#include "fesa/elements/nonlinear_mitc4_inertia.h"
#include "fesa/elements/nonlinear_mitc4_shell.h"
#include "fesa/model/model_types.h"
#include "elements/mitc4_drilling_internal.h"
#include "elements/nonlinear_mitc4_configuration_internal.h"

namespace {

// User-approved physical/numerical comparison tolerance for the local experiment.
constexpr double kNumericalTolerance = 1.0e-3;

using Coordinates = std::array<double, 3>;

fesa::Node Node(const std::int64_t label, const Coordinates coordinates) {
  return {{"Shell-Instance", label, std::to_string(label)},
          coordinates,
          {"nonlinear-mitc4-tangent.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::NonlinearMitc4Shell ProductionShell(const std::array<fesa::Node, 4>& nodes) {
  const fesa::ShellSection section{
      "Section", 0.2, 0U, {"nonlinear-mitc4-tangent.inp", 20U}};
  const fesa::LinearElasticMaterial material{
      "Material", 120.0, 0.25, {"nonlinear-mitc4-tangent.inp", 21U}};
  const auto shell = fesa::NonlinearMitc4Shell::Create(
      NodePointers(nodes), Frames(), section, material);
  if (!shell.HasValue()) {
    throw std::runtime_error{shell.GetStatus().Diagnostics()[0U].code};
  }
  return shell.Value();
}

fesa::NonlinearMitc4Shell LegacyShell(const std::array<fesa::Node, 4>& nodes) {
  return fesa::mitc4_drilling_internal::TestAccess::WithMode(
      ProductionShell(nodes), fesa::Mitc4DrillingMode::kLegacyCorotatedTwist);
}

fesa::NonlinearMitc4TrialDofs Trial() {
  fesa::NonlinearMitc4TrialDofs trial{};
  trial.displacement[0U] = {-0.02, 0.01, 0.03};
  trial.displacement[1U] = {0.01, -0.03, -0.01};
  trial.displacement[2U] = {0.03, 0.02, 0.04};
  trial.displacement[3U] = {-0.01, 0.04, -0.02};
  trial.rotation_vector[0U] = {0.03, -0.02, 0.01};
  trial.rotation_vector[1U] = {-0.01, 0.04, -0.02};
  trial.rotation_vector[2U] = {0.02, 0.01, 0.03};
  trial.rotation_vector[3U] = {-0.03, 0.02, -0.01};
  return trial;
}

fesa::NonlinearMitc4TrialDofs AddScaled(
    const fesa::NonlinearMitc4TrialDofs& trial,
    const fesa::NonlinearMitc4TrialDofs& direction, const double scale) {
  fesa::NonlinearMitc4TrialDofs result = trial;
  for (std::size_t node = 0U; node < 4U; ++node) {
    result.displacement[node] =
        result.displacement[node] + scale * direction.displacement[node];
    result.rotation_vector[node] =
        result.rotation_vector[node] + scale * direction.rotation_vector[node];
  }
  return result;
}

fesa::NonlinearMitc4TrialDofs Direction() {
  fesa::NonlinearMitc4TrialDofs direction{};
  direction.displacement[0U] = {0.5, -0.2, 0.1};
  direction.displacement[1U] = {-0.1, 0.3, -0.4};
  direction.displacement[2U] = {0.2, 0.1, 0.5};
  direction.displacement[3U] = {-0.3, -0.4, -0.2};
  direction.rotation_vector[0U] = {0.2, -0.1, 0.3};
  direction.rotation_vector[1U] = {-0.3, 0.2, -0.1};
  direction.rotation_vector[2U] = {0.1, 0.3, -0.2};
  direction.rotation_vector[3U] = {-0.2, -0.3, 0.1};
  return direction;
}

fesa::Vector ToVector(const fesa::NonlinearMitc4TrialDofs& trial) {
  fesa::Vector result{24U};
  for (std::size_t node = 0U; node < 4U; ++node) {
    for (std::size_t component = 0U; component < 3U; ++component) {
      result[6U * node + component] = trial.displacement[node][component];
      result[6U * node + 3U + component] =
          trial.rotation_vector[node][component];
    }
  }
  return result;
}

double ReferenceLength(const std::array<fesa::Node, 4>& nodes) {
  double length = 0.0;
  for (std::size_t first = 0U; first < nodes.size(); ++first) {
    for (std::size_t second = first + 1U; second < nodes.size(); ++second) {
      length = (std::max)(length, (fesa::Vector3{nodes[first].coordinates} -
                                   fesa::Vector3{nodes[second].coordinates})
                                      .Norm());
    }
  }
  return length;
}

fesa::NonlinearMitc4TrialDofs WorkConjugateDirection(
    const std::array<fesa::Node, 4>& nodes) {
  auto direction = Direction();
  double maximum_component = 0.0;
  for (std::size_t node = 0U; node < 4U; ++node) {
    for (std::size_t component = 0U; component < 3U; ++component) {
      maximum_component =
          (std::max)(maximum_component,
                     std::abs(direction.displacement[node][component]));
      maximum_component =
          (std::max)(maximum_component,
                     std::abs(direction.rotation_vector[node][component]));
    }
  }
  const double length = ReferenceLength(nodes);
  for (std::size_t node = 0U; node < 4U; ++node) {
    direction.displacement[node] =
        length * direction.displacement[node] / maximum_component;
    direction.rotation_vector[node] =
        direction.rotation_vector[node] / maximum_component;
  }
  return direction;
}

fesa::Vector WorkConjugate(const fesa::Vector& vector, const double length) {
  fesa::Vector result = vector;
  for (std::size_t node = 0U; node < 4U; ++node) {
    for (std::size_t component = 0U; component < 3U; ++component) {
      result[6U * node + component] *= length;
    }
  }
  return result;
}

double ChannelNorm(const fesa::Vector& vector, const bool translation) {
  double sum = 0.0;
  const std::size_t offset = translation ? 0U : 3U;
  for (std::size_t node = 0U; node < 4U; ++node) {
    for (std::size_t component = 0U; component < 3U; ++component) {
      const double value = vector[6U * node + offset + component];
      sum += value * value;
    }
  }
  return std::sqrt(sum);
}

double FrobeniusNorm(const fesa::Matrix& matrix) {
  double sum = 0.0;
  for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
    for (std::size_t column = 0U; column < matrix.Columns(); ++column) {
      sum += matrix(row, column) * matrix(row, column);
    }
  }
  return std::sqrt(sum);
}

fesa::Matrix WorkConjugateTangent(const fesa::Matrix& tangent,
                                  const double length) {
  fesa::Matrix result = tangent;
  for (std::size_t row = 0U; row < tangent.Rows(); ++row) {
    const double row_scale = row % 6U < 3U ? length : 1.0;
    for (std::size_t column = 0U; column < tangent.Columns(); ++column) {
      const double column_scale = column % 6U < 3U ? length : 1.0;
      result(row, column) *= row_scale * column_scale;
    }
  }
  return result;
}

void ExpectRelative(const double actual, const double expected) {
  if (expected == 0.0) {
    EXPECT_LE(std::abs(actual), kNumericalTolerance);
    return;
  }
  EXPECT_LE(std::abs(actual - expected) / std::abs(expected), kNumericalTolerance);
}

// NSNR-ELEM-TAN-001
TEST(NonlinearMitc4Residual, EvaluatesResidualTangentAndRecoveryAtOneTrial) {
  const auto nodes = Nodes();
  const auto shell = LegacyShell(nodes);
  const auto trial = Trial();
  const auto recovery = shell.RecoverState(trial);
  const auto result = fesa::EvaluateInternalResidualAndTangent(shell, trial);

  ASSERT_TRUE(recovery.HasValue());
  ASSERT_TRUE(result.HasValue());
  EXPECT_EQ(result.Value().internal_residual.Size(), 24U);
  EXPECT_EQ(result.Value().internal_tangent.Rows(), 24U);
  EXPECT_EQ(result.Value().internal_tangent.Columns(), 24U);
  EXPECT_DOUBLE_EQ(result.Value().physical_strain_energy,
                   result.Value().recovery.physical_strain_energy);
  EXPECT_DOUBLE_EQ(result.Value().recovery.physical_strain_energy,
                   recovery.Value().physical_strain_energy);
  for (std::size_t index = 0U; index < 24U; ++index) {
    EXPECT_DOUBLE_EQ(result.Value().internal_residual[index],
                     result.Value().internal_force[index]);
  }
  for (std::size_t point = 0U;
       point < result.Value().recovery.gauss_points.size(); ++point) {
    EXPECT_EQ(result.Value().recovery.gauss_points[point].generalized_strain,
              recovery.Value().gauss_points[point].generalized_strain);
    EXPECT_EQ(result.Value().recovery.gauss_points[point].section_resultant,
              recovery.Value().gauss_points[point].section_resultant);
  }
}

// NSNR-ELEM-TAN-002 and NSNR-ELEM-TAN-003
TEST(NonlinearMitc4Tangent, MatchesTheApprovedDirectionalSchedule) {
  const auto nodes = Nodes();
  const auto shell = LegacyShell(nodes);
  const auto trial = Trial();
  const auto direction = WorkConjugateDirection(nodes);
  const double length = ReferenceLength(nodes);
  const auto equilibrium =
      fesa::EvaluateInternalResidualAndTangent(shell, trial);
  ASSERT_TRUE(equilibrium.HasValue());

  const fesa::Vector direction_vector = ToVector(direction);
  const fesa::Vector analytical = WorkConjugate(
      equilibrium.Value().internal_tangent.Multiply(direction_vector), length);
  const std::array<double, 8> steps{1.0e-3, 3.0e-4, 1.0e-4, 3.0e-5,
                                    1.0e-5, 3.0e-6, 1.0e-6, 3.0e-7};
  std::array<double, 8> force_errors{};
  std::array<double, 8> moment_errors{};
  ASSERT_GT(ChannelNorm(analytical, true), 0.0);
  ASSERT_GT(ChannelNorm(analytical, false), 0.0);
  for (std::size_t index = 0U; index < steps.size(); ++index) {
    const auto plus = fesa::EvaluateInternalResidualAndTangent(
        shell, AddScaled(trial, direction, steps[index]));
    const auto minus = fesa::EvaluateInternalResidualAndTangent(
        shell, AddScaled(trial, direction, -steps[index]));
    ASSERT_TRUE(plus.HasValue());
    ASSERT_TRUE(minus.HasValue());
    fesa::Vector derivative = plus.Value().internal_residual;
    derivative.Axpy(-1.0, minus.Value().internal_residual);
    derivative.Scale(0.5 / steps[index]);
    derivative = WorkConjugate(derivative, length);
    derivative.Axpy(-1.0, analytical);
    force_errors[index] =
        ChannelNorm(derivative, true) / ChannelNorm(analytical, true);
    moment_errors[index] =
        ChannelNorm(derivative, false) / ChannelNorm(analytical, false);
  }

  std::size_t force_order_count = 0U;
  std::size_t moment_order_count = 0U;
  for (std::size_t index = 1U; index < force_errors.size(); ++index) {
    if (force_errors[index - 1U] > 0.0 && force_errors[index] > 0.0) {
      const double order =
          std::log(force_errors[index - 1U] / force_errors[index]) /
          std::log(steps[index - 1U] / steps[index]);
      if (order >= 1.5) {
        ++force_order_count;
      }
    }
    if (moment_errors[index - 1U] > 0.0 && moment_errors[index] > 0.0) {
      const double order =
          std::log(moment_errors[index - 1U] / moment_errors[index]) /
          std::log(steps[index - 1U] / steps[index]);
      if (order >= 1.5) {
        ++moment_order_count;
      }
    }
  }
  EXPECT_GE(force_order_count, 2U);
  EXPECT_GE(moment_order_count, 2U);
  EXPECT_LE(force_errors[3U], kNumericalTolerance);
  EXPECT_LE(moment_errors[3U], kNumericalTolerance);
}

// NSNR-ELEM-TAN-004 and NSNR-ELEM-TAN-005
TEST(NonlinearMitc4Tangent, PreservesRigidTranslationNullspaceAndSymmetry) {
  const auto nodes = Nodes();
  const auto equilibrium =
      fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), {});
  ASSERT_TRUE(equilibrium.HasValue());
  const double length = ReferenceLength(nodes);

  fesa::NonlinearMitc4TrialDofs rigid_translation{};
  for (auto& displacement : rigid_translation.displacement) {
    displacement = {0.3, -0.2, 0.4};
  }
  const fesa::Vector null_action =
      WorkConjugate(equilibrium.Value().internal_tangent.Multiply(
                        ToVector(rigid_translation)),
                    length);
  EXPECT_LE(null_action.Norm(), kNumericalTolerance);

  const fesa::Matrix scaled_tangent =
      WorkConjugateTangent(equilibrium.Value().internal_tangent, length);
  fesa::Matrix antisymmetric = scaled_tangent;
  for (std::size_t row = 0U; row < antisymmetric.Rows(); ++row) {
    for (std::size_t column = 0U; column < antisymmetric.Columns(); ++column) {
      antisymmetric(row, column) -= scaled_tangent(column, row);
    }
  }
  ASSERT_GT(FrobeniusNorm(scaled_tangent), 0.0);
  EXPECT_LE(FrobeniusNorm(antisymmetric) / FrobeniusNorm(scaled_tangent),
            kNumericalTolerance);
}

// NSNR-ELEM-TAN-006
TEST(NonlinearMitc4Residual, DoesNotAcceptExternalLoadShares) {
  const auto nodes = Nodes();
  const auto first =
      fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), Trial());
  const auto second =
      fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), Trial());
  ASSERT_TRUE(first.HasValue());
  ASSERT_TRUE(second.HasValue());
  for (std::size_t index = 0U; index < 24U; ++index) {
    EXPECT_DOUBLE_EQ(first.Value().internal_residual[index],
                     second.Value().internal_residual[index]);
  }
}

// NSNR-ELEM-TAN-007
TEST(NonlinearMitc4Objectivity, FiniteRigidMotionHasZeroInternalForce) {
  const auto nodes = Nodes();
  const fesa::Vector3 rotation{0.3, -0.2, 0.1};
  const fesa::Matrix tensor = fesa::ExpSkew(rotation);
  fesa::NonlinearMitc4TrialDofs trial{};
  for (std::size_t node = 0U; node < nodes.size(); ++node) {
    const fesa::Vector3 initial{nodes[node].coordinates};
    const fesa::Vector3 rotated{
        tensor(0U, 0U) * initial.X() + tensor(0U, 1U) * initial.Y() +
            tensor(0U, 2U) * initial.Z(),
        tensor(1U, 0U) * initial.X() + tensor(1U, 1U) * initial.Y() +
            tensor(1U, 2U) * initial.Z(),
        tensor(2U, 0U) * initial.X() + tensor(2U, 1U) * initial.Y() +
            tensor(2U, 2U) * initial.Z()};
    trial.displacement[node] =
        rotated + fesa::Vector3{0.7, -0.4, 1.2} - initial;
    trial.rotation_vector[node] = rotation;
  }
  const auto result =
      fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), trial);
  ASSERT_TRUE(result.HasValue());
  EXPECT_LE(result.Value().internal_residual.Norm(), kNumericalTolerance);
  EXPECT_LE(result.Value().physical_strain_energy, kNumericalTolerance);
}

// NSNR-ELEM-TAN-008
TEST(NonlinearMitc4Patch, RetainsPhysicalRecoveryIdentities) {
  const auto nodes = Nodes();
  const auto shell = LegacyShell(nodes);
  constexpr double kMembraneStrain = 1.0e-2;
  const double green_strain =
      kMembraneStrain + 0.5 * kMembraneStrain * kMembraneStrain;
  constexpr double kYoungsModulus = 120.0;
  constexpr double kPoissonRatio = 0.25;
  constexpr double kThickness = 0.2;
  const double plane_stress =
      kYoungsModulus / (1.0 - kPoissonRatio * kPoissonRatio);
  const double s11 = plane_stress * green_strain;
  const double s22 = plane_stress * kPoissonRatio * green_strain;
  fesa::NonlinearMitc4TrialDofs trial{};
  for (std::size_t node = 0U; node < nodes.size(); ++node) {
    trial.displacement[node] = {kMembraneStrain * nodes[node].coordinates[0U],
                                0.0, 0.0};
  }
  const auto result = fesa::EvaluateInternalResidualAndTangent(shell, trial);
  ASSERT_TRUE(result.HasValue());
  for (const auto& point : result.Value().recovery.gauss_points) {
    ExpectRelative(point.generalized_strain[0U], green_strain);
    for (std::size_t component = 1U; component < 8U; ++component) {
      ExpectRelative(point.generalized_strain[component], 0.0);
    }
    ExpectRelative(point.section_resultant[0U], kThickness * s11);
    ExpectRelative(point.section_resultant[1U], kThickness * s22);
    for (std::size_t component = 2U; component < 8U; ++component) {
      ExpectRelative(point.section_resultant[component], 0.0);
    }
    for (const auto& stress : point.stress) {
      ExpectRelative(stress.components[0U], s11);
      ExpectRelative(stress.components[1U], s22);
      ExpectRelative(stress.components[2U], 0.0);
    }
  }
}

// NSNR-ELEM-TAN-009
TEST(NonlinearMitc4Residual, RejectsNonfiniteTrialData) {
  const auto nodes = Nodes();
  auto trial = Trial();
  trial.rotation_vector[2U] = {std::numeric_limits<double>::quiet_NaN(), 0.0,
                               0.0};
  const auto result =
      fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), trial);
  ASSERT_FALSE(result.HasValue());
  EXPECT_EQ(result.GetStatus().Diagnostics()[0U].code,
            "invalid-nonlinear-shell-trial");

  auto overflow = Trial();
  overflow.displacement[1U] = {1.0e300, 0.0, 0.0};
  const auto overflow_result =
      fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), overflow);
  EXPECT_FALSE(overflow_result.HasValue());
}

void CheckThirdSweep(const fesa::NonlinearMitc4TrialDofs& trial) {
  const auto nodes = Nodes();
  const auto shell = LegacyShell(nodes);
  auto direction = Direction();
  double norm = 0.0;
  const auto unscaled = ToVector(direction);
  for (std::size_t i = 0U; i < 24U; ++i) {
    norm = std::hypot(norm, unscaled[i]);
  }
  for (std::size_t node = 0U; node < 4U; ++node) {
    direction.displacement[node] = 2.0 * direction.displacement[node] / norm;
    direction.rotation_vector[node] = direction.rotation_vector[node] / norm;
  }
  const auto velocity = ToVector(direction);
  const auto exact =
      fesa::EvaluateInternalResidualAndTangent(shell, trial, &velocity);
  ASSERT_TRUE(exact.HasValue());
  unsigned int consecutive = 0U;
  bool passed = false;
  for (int power = 8; power <= 20; power += 2) {
    const double h = std::ldexp(1.0, -power);
    const auto plus = fesa::EvaluateInternalResidualAndTangent(
        shell, AddScaled(trial, direction, h));
    const auto minus = fesa::EvaluateInternalResidualAndTangent(
        shell, AddScaled(trial, direction, -h));
    ASSERT_TRUE(plus.HasValue());
    ASSERT_TRUE(minus.HasValue());
    double error = 0.0;
    double analytic_norm = 0.0;
    double fd_norm = 0.0;
    for (std::size_t i = 0U; i < 24U; ++i) {
      double fd = 0.0;
      for (std::size_t j = 0U; j < 24U; ++j) {
        fd += (plus.Value().internal_tangent(i, j) -
               minus.Value().internal_tangent(i, j)) *
              velocity[j] / (2.0 * h);
      }
      const double d = i % 6U < 3U ? 2.0 : 1.0;
      const double analytic = exact.Value().force_directional_second[i];
      error = std::hypot(error, d * (analytic - fd));
      analytic_norm = std::hypot(analytic_norm, d * analytic);
      fd_norm = std::hypot(fd_norm, d * fd);
    }
    if (power >= 10 && power <= 18 &&
        error / (analytic_norm + fd_norm) <= kNumericalTolerance) {
      ++consecutive;
      passed = passed || consecutive >= 2U;
    } else {
      consecutive = 0U;
    }
  }
  EXPECT_TRUE(passed);
  EXPECT_GE(exact.Value().drilling_energy, 0.0);
  for (std::size_t i = 0U; i < 24U; ++i) {
    EXPECT_DOUBLE_EQ(
        exact.Value().internal_force[i],
        exact.Value().physical_force[i] + exact.Value().drilling_force[i]);
  }
}

TEST(NonlinearMitc4Residual,
     ZeroTrialOnTranslatedDecimalMeshHasNoArtificialForce) {
  auto nodes = Nodes();
  nodes[0].coordinates = {3.84615385, 2.30769231, 0.0};
  nodes[1].coordinates = {4.23076923, 2.30769231, 0.0};
  nodes[2].coordinates = {4.23076923, 2.69230769, 0.0};
  nodes[3].coordinates = {3.84615385, 2.69230769, 0.0};
  const auto shell = LegacyShell(nodes);
  const auto evaluated = fesa::EvaluateInternalResidualAndTangent(shell, {});
  ASSERT_TRUE(evaluated.HasValue());
  for (std::size_t i = 0; i < 24U; ++i)
    EXPECT_DOUBLE_EQ(evaluated.Value().internal_force[i], 0.0);
  EXPECT_DOUBLE_EQ(evaluated.Value().drilling_energy, 0.0);
}

TEST(NonlinearMitc4Residual, ExactDirectionalThirdGradientMatchesTangentSweep) {
  CheckThirdSweep(Trial());
}

TEST(NonlinearMitc4Residual, SmallRotationPhysicalForceRetainsDrillingNullIdentity) {
  // Exact local state of C08b h=1/32768, increment4: the rotation has just
  // crossed the original 1e-4 coefficient branch and its physical null leaked.
  auto nodes = Nodes();
  nodes[1U].coordinates = {0.0, -1.0, 0.0};
  nodes[2U].coordinates = {0.0, 0.0, 0.0};
  nodes[3U].coordinates = {-1.0, 0.0, 0.0};
  const fesa::ShellSection section{"Centered", 0.2, 0U, {}};
  const fesa::LinearElasticMaterial material{"Elastic", 1000.0, 0.3, {}};
  const auto shell = fesa::NonlinearMitc4Shell::Create(NodePointers(nodes), Frames(),
                                                      section, material);
  ASSERT_TRUE(shell.HasValue());
  const std::array<std::array<double, 6>, 4> coordinates{{
      {6.4086914072366947e-05, -7.9345703110834861e-05, -1.4343267109683367e-05, 9.7657633876829583e-05, 7.324020770488069e-05, 4.8828020264083388e-05},
      {6.4086914054701282e-05, -3.0517578102879528e-05, -9.0026855471989683e-05, 9.7658037630745924e-05, 7.3240207702551634e-05, 4.8828064949797786e-05},
      {1.5258789062740791e-05, -3.0517578125180581e-05, 7.6293945312499653e-06, 9.7658037632949931e-05, 7.3239803947312715e-05, 4.8828124536813954e-05},
      {1.5258789080115012e-05, -7.9345703132764625e-05, 8.0871582026927324e-05, 9.7657633879033603e-05, 7.3239803949641825e-05, 4.8828079851099556e-05}
  }};
  fesa::NonlinearMitc4TrialDofs trial{};
  for (std::size_t node = 0U; node < 4U; ++node) {
    const auto& q = coordinates[node];
    trial.displacement[node] = {q[0U], q[1U], q[2U]};
    trial.rotation_vector[node] = {q[3U], q[4U], q[5U]};
  }
  const fesa::Vector direction(24U);
  const auto response = fesa::EvaluateInternalResidualAndTangent(shell.Value(), trial, &direction);
  ASSERT_TRUE(response.HasValue());
  for (std::size_t node = 0U; node < 4U; ++node) {
    SCOPED_TRACE(node);
    const auto rotation = fesa::EvaluateDynamicRotation(trial.rotation_vector[node],
                                                       {0.0, 0.0, 1.0}, {}, {});
    ASSERT_TRUE(rotation.HasValue());
    double physical = 0.0, physical_scale = 0.0, projected = 0.0, total_scale = 0.0;
    for (std::size_t c = 0U; c < 3U; ++c) {
      const auto full = 6U * node + 3U + c;
      const double n = rotation.Value().drill_null[c];
      physical += n * response.Value().physical_force[full];
      physical_scale += std::abs(n * response.Value().physical_force[full]);
      projected += n * response.Value().internal_force[full];
      total_scale += std::abs(n * response.Value().internal_force[full]);
    }
    ASSERT_GT(physical_scale, 0.0);
    ASSERT_GT(total_scale, 0.0);
    EXPECT_LE(std::abs(physical) / physical_scale, kNumericalTolerance);
    EXPECT_LE(std::abs(projected - response.Value().drilling_constraint[node]) / total_scale,
              kNumericalTolerance);
  }
}

TEST(NonlinearMitc4Residual, ThirdGradientIsStableAcrossSmallAngleBranches) {
  for (const double angle : {1e-4 * (1 - std::ldexp(1., -10)),
                             1e-4 * (1 + std::ldexp(1., -10)), 2e-4, 1e-3,
                             .5 * (1 - std::ldexp(1., -12)),
                             .5 * (1 + std::ldexp(1., -12))}) {
    SCOPED_TRACE(angle);
    auto trial = Trial();
    for (auto& phi : trial.rotation_vector) phi = {angle, 0., 0.};
    CheckThirdSweep(trial);
  }
}

TEST(NonlinearMitc4Residual,
     ExactDrillingConstraintHasIndependentForceGradientAndHessian) {
  for (bool warped : {false, true}) {
    auto nodes = Nodes();
    if (warped) nodes[2U].coordinates[2U] = 0.1;
    const auto shell = LegacyShell(nodes);
    auto trial = Trial();
    trial.rotation_vector[0U] = {.6, -.2, .1};
    const auto direction = Direction();
    const auto velocity = ToVector(direction);
    const auto exact =
        fesa::EvaluateInternalResidualAndTangent(shell, trial, &velocity);
    ASSERT_TRUE(exact.HasValue());
    for (std::size_t node = 0U; node < 4U; ++node) {
      const auto rotation = fesa::EvaluateDynamicRotation(
          trial.rotation_vector[node], {0., 0., 1.},
          direction.rotation_vector[node], {});
      ASSERT_TRUE(rotation.HasValue());
      double psi = 0.;
      for (std::size_t c = 0U; c < 3U; ++c)
        psi += rotation.Value().drill_null[c] *
               exact.Value().internal_force[6U * node + 3U + c];
      EXPECT_NEAR(exact.Value().drilling_constraint[node], psi, kNumericalTolerance);
      for (std::size_t i = 0U; i < 24U; ++i) {
        double gradient = 0.;
        for (std::size_t c = 0U; c < 3U; ++c) {
          gradient += rotation.Value().drill_null[c] *
                      exact.Value().internal_tangent(6U * node + 3U + c, i);
          if (i / 6U == node && i % 6U >= 3U)
            gradient += rotation.Value().null_gradient(c, i % 6U - 3U) *
                        exact.Value().internal_force[6U * node + 3U + c];
        }
        EXPECT_NEAR(exact.Value().drilling_constraint_gradient(node, i),
                    gradient, kNumericalTolerance);
      }
    }
    unsigned consecutive = 0U;
    bool passed = false;
    for (int power = 8; power <= 20; power += 2) {
      const double h = std::ldexp(1., -power);
      const auto plus = fesa::EvaluateInternalResidualAndTangent(
          shell, AddScaled(trial, direction, h), &velocity);
      const auto minus = fesa::EvaluateInternalResidualAndTangent(
          shell, AddScaled(trial, direction, -h), &velocity);
      ASSERT_TRUE(plus.HasValue());
      ASSERT_TRUE(minus.HasValue());
      double error = 0., anorm = 0., fnorm = 0.;
      for (std::size_t node = 0U; node < 4U; ++node)
        for (std::size_t i = 0U; i < 24U; ++i) {
          double analytic = 0.;
          for (std::size_t j = 0U; j < 24U; ++j)
            analytic += exact.Value().drilling_constraint_hessian[node](i, j) *
                        velocity[j];
          const double fd =
              (plus.Value().drilling_constraint_gradient(node, i) -
               minus.Value().drilling_constraint_gradient(node, i)) /
              (2 * h);
          const double scale = i % 6U < 3U ? 2. : 1.;
          error = std::hypot(error, scale * (analytic - fd));
          anorm = std::hypot(anorm, scale * analytic);
          fnorm = std::hypot(fnorm, scale * fd);
        }
      if (power >= 10 && power <= 18 && error / (anorm + fnorm) <= kNumericalTolerance) {
        ++consecutive;
        passed = passed || consecutive >= 2U;
      } else
        consecutive = 0U;
    }
    EXPECT_TRUE(passed);
  }
}

TEST(NonlinearMitc4Residual,
     ExposesExactDrillingConstraintThirdDirectionalGradient) {
  const auto nodes = Nodes();
  const auto shell = LegacyShell(nodes);
  const auto trial = Trial();
  const auto direction = Direction();
  const auto velocity = ToVector(direction);
  const auto exact =
      fesa::EvaluateInternalResidualAndTangent(shell, trial, &velocity);
  ASSERT_TRUE(exact.HasValue());
  constexpr double h = 1.0e-4;
  const auto plus = fesa::EvaluateInternalResidualAndTangent(
      shell, AddScaled(trial, direction, h), &velocity);
  const auto minus = fesa::EvaluateInternalResidualAndTangent(
      shell, AddScaled(trial, direction, -h), &velocity);
  ASSERT_TRUE(plus.HasValue());
  ASSERT_TRUE(minus.HasValue());

  for (std::size_t node = 0U; node < 4U; ++node) {
    ASSERT_EQ(exact.Value().drilling_constraint_third_directional_gradient[node]
                  .Size(),
              24U);
    for (std::size_t row = 0U; row < 24U; ++row) {
      double plus_hessian_direction = 0.0;
      double minus_hessian_direction = 0.0;
      for (std::size_t column = 0U; column < 24U; ++column) {
        plus_hessian_direction +=
            plus.Value().drilling_constraint_hessian[node](row, column) *
            velocity[column];
        minus_hessian_direction +=
            minus.Value().drilling_constraint_hessian[node](row, column) *
            velocity[column];
      }
      const double finite_difference =
          (plus_hessian_direction - minus_hessian_direction) / (2.0 * h);
      EXPECT_NEAR(
          exact.Value().drilling_constraint_third_directional_gradient[node]
              [row],
          finite_difference, kNumericalTolerance);
    }
  }
}

TEST(NonlinearMitc4Residual,
     ExternalDrillingConstraintIncludesForceTimeAndCurvatureTerms) {
  const auto nodes = Nodes();
  const auto shell = LegacyShell(nodes);
  const fesa::Vector3 phi{.6, -.2, .1};
  const fesa::Vector3 phi_dot{.2, -.1, .3};
  const fesa::Vector3 force{1.7, -.8, .4};
  const fesa::Vector3 force_rate{-.3, .6, .2};
  const fesa::Vector3 force_acceleration{.5, -.4, .1};
  const auto exact = shell.EvaluateExternalDrillingConstraint(
      0U, phi, phi_dot, force, force_rate, force_acceleration);
  ASSERT_TRUE(exact.HasValue());
  constexpr double h = 1.0e-4;
  const auto at_time = [&](const double time) {
    return shell.EvaluateExternalDrillingConstraint(
        0U, phi + time * phi_dot, phi_dot,
        force + time * force_rate + 0.5 * time * time * force_acceleration,
        force_rate + time * force_acceleration, force_acceleration);
  };
  const auto plus = at_time(h);
  const auto minus = at_time(-h);
  ASSERT_TRUE(plus.HasValue());
  ASSERT_TRUE(minus.HasValue());
  EXPECT_NEAR(exact.Value().directional_first,
              (plus.Value().value - minus.Value().value) / (2.0 * h),
              kNumericalTolerance);
  EXPECT_NEAR(exact.Value().directional_second,
              (plus.Value().value - 2.0 * exact.Value().value +
               minus.Value().value) /
                  (h * h),
              kNumericalTolerance);
  for (std::size_t row = 0U; row < 3U; ++row) {
    EXPECT_NEAR(exact.Value().gradient_directional_first[row],
                (plus.Value().gradient[row] - minus.Value().gradient[row]) /
                    (2.0 * h),
                kNumericalTolerance);
    EXPECT_NEAR(
        exact.Value().gradient_directional_second[row],
        (plus.Value().gradient[row] - 2.0 * exact.Value().gradient[row] +
         minus.Value().gradient[row]) /
            (h * h),
        kNumericalTolerance);
  }
  for (std::size_t column = 0U; column < 3U; ++column) {
    fesa::Vector3 basis{};
    if (column == 0U) basis = {1.0, 0.0, 0.0};
    if (column == 1U) basis = {0.0, 1.0, 0.0};
    if (column == 2U) basis = {0.0, 0.0, 1.0};
    const auto positive = shell.EvaluateExternalDrillingConstraint(
        0U, phi + h * basis, phi_dot, force, force_rate, force_acceleration);
    const auto negative = shell.EvaluateExternalDrillingConstraint(
        0U, phi - h * basis, phi_dot, force, force_rate, force_acceleration);
    ASSERT_TRUE(positive.HasValue());
    ASSERT_TRUE(negative.HasValue());
    for (std::size_t row = 0U; row < 3U; ++row) {
      EXPECT_NEAR(exact.Value().hessian(row, column),
                  (positive.Value().gradient[row] -
                   negative.Value().gradient[row]) /
                      (2.0 * h),
                  kNumericalTolerance);
    }
  }
}

TEST(NonlinearMitc4Residual,
     ExternalDrillingConstraintRejectsInvalidRotationAndHasZeroForceValue) {
  const auto shell = LegacyShell(Nodes());
  const auto zero = shell.EvaluateExternalDrillingConstraint(
      0U, {.2, -.1, .3}, {.1, .2, -.3}, {}, {}, {});
  ASSERT_TRUE(zero.HasValue());
  EXPECT_DOUBLE_EQ(zero.Value().value, 0.0);
  EXPECT_DOUBLE_EQ(zero.Value().directional_first, 0.0);
  EXPECT_DOUBLE_EQ(zero.Value().directional_second, 0.0);
  for (std::size_t row = 0U; row < 3U; ++row) {
    EXPECT_DOUBLE_EQ(zero.Value().gradient[row], 0.0);
    EXPECT_DOUBLE_EQ(zero.Value().gradient_directional_first[row], 0.0);
    EXPECT_DOUBLE_EQ(zero.Value().gradient_directional_second[row], 0.0);
  }
  const auto invalid = shell.EvaluateExternalDrillingConstraint(
      0U, {std::acos(-1.0), 0.0, 0.0}, {}, {}, {}, {});
  EXPECT_FALSE(invalid.HasValue());
  EXPECT_EQ(invalid.GetStatus().Diagnostics()[0U].code,
            "invalid-nonlinear-shell-rotation");
}

TEST(NonlinearMitc4Residual, TinyRotationPreservesPhysicalForceAndRecoveryEnergy) {
  for (const bool rotated : {false, true}) {
    for (const bool warped : {false, true}) {
      SCOPED_TRACE(rotated);
      SCOPED_TRACE(warped);
      auto nodes = Nodes();
      if (warped) nodes[2U].coordinates[2U] = 0.1;
      auto frames = Frames();
      const auto rotation = fesa::ExpSkew(
          rotated ? fesa::Vector3{0.4, -0.2, 0.7} : fesa::Vector3{});
      const auto transform = [&](const Coordinates& values) {
        Coordinates result{};
        for (std::size_t i = 0U; i < 3U; ++i)
          for (std::size_t j = 0U; j < 3U; ++j)
            result[i] += rotation(i, j) * values[j];
        return result;
      };
      for (std::size_t node = 0U; node < 4U; ++node) {
        nodes[node].coordinates = transform(nodes[node].coordinates);
        frames[node].director = transform(frames[node].director);
        frames[node].tangent_a = transform(frames[node].tangent_a);
        frames[node].tangent_b = transform(frames[node].tangent_b);
      }
      const fesa::ShellSection section{"Section", 0.2, 0U, {"tiny.inp", 1U}};
      const fesa::LinearElasticMaterial material{
          "Material", 120.0, 0.25, {"tiny.inp", 2U}};
      const auto shell = fesa::NonlinearMitc4Shell::Create(
          NodePointers(nodes), frames, section, material);
      ASSERT_TRUE(shell.HasValue());
      fesa::NonlinearMitc4TrialDofs tiny{};
      fesa::Vector increment{24U};
      for (std::size_t node = 0U; node < 4U; ++node) {
        tiny.rotation_vector[node] =
            1e-17 * fesa::Vector3{frames[node].tangent_a};
        for (std::size_t c = 0U; c < 3U; ++c)
          increment[6U * node + 3U + c] = tiny.rotation_vector[node][c];
      }
      const auto zero =
          fesa::EvaluateInternalResidualAndTangent(shell.Value(), {});
      const auto actual =
          fesa::EvaluateInternalResidualAndTangent(shell.Value(), tiny);
      ASSERT_TRUE(zero.HasValue());
      ASSERT_TRUE(actual.HasValue());
      EXPECT_DOUBLE_EQ(zero.Value().physical_strain_energy, 0.0);
      double error = 0.0;
      double scale = 0.0;
      double quadratic_energy = 0.0;
      for (std::size_t row = 0U; row < 24U; ++row) {
        EXPECT_DOUBLE_EQ(zero.Value().physical_force[row], 0.0);
        double expected = 0.0;
        for (std::size_t column = 0U; column < 24U; ++column)
          expected += (zero.Value().internal_tangent(row, column) -
                       zero.Value().drilling_tangent(row, column)) *
                      increment[column];
        error = std::max(error,
                         std::abs(actual.Value().physical_force[row] - expected));
        scale = std::max(scale, std::abs(expected));
        quadratic_energy += 0.5 * increment[row] * expected;
      }
      ASSERT_GT(scale, 0.0);
      ASSERT_GT(quadratic_energy, 0.0);
      EXPECT_LE(error / scale, kNumericalTolerance);
      EXPECT_GT(actual.Value().physical_strain_energy, 0.0);
      EXPECT_LE(std::abs(actual.Value().physical_strain_energy - quadratic_energy) /
                    quadratic_energy,
                kNumericalTolerance);
      if (!warped) {
        // Uniform director tilt has .5 * (5/6) * G * t * area * theta^2.
        constexpr double shear_energy = 1.6e-33;
        EXPECT_LE(std::abs(actual.Value().physical_strain_energy - shear_energy) /
                      shear_energy,
                  kNumericalTolerance);
      }
    }
  }
}

// Independent fixed-initial-director spring oracle: a total-chart rigid spin
// carries numerical energy even though the physical shell is undeformed.
TEST(Mitc4DrillingSelection, ProductionUsesDirectInitialDirectorSpring) {
  const auto nodes = Nodes();
  const auto shell = ProductionShell(nodes);
  const fesa::ShellSection section{"Section", 0.2, 0U, {"spring.inp", 1U}};
  const fesa::LinearElasticMaterial material{
      "Material", 120.0, 0.25, {"spring.inp", 2U}};
  const auto linear = fesa::Mitc4Shell::Create(
      NodePointers(nodes),
      {{{0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}}}, section,
      material);
  ASSERT_TRUE(linear.HasValue());
  const auto stiffness = linear.Value().Stiffness();
  ASSERT_TRUE(stiffness.HasValue());
  double rotary_trace = 0.0;
  for (std::size_t node = 0; node < 4; ++node)
    for (std::size_t c = 3; c < 5; ++c)
      rotary_trace +=
          stiffness.Value().physical_local20(5 * node + c, 5 * node + c);
  const double kd = 1e-3 * rotary_trace / 8.;
  fesa::NonlinearMitc4TrialDofs trial{};
  const fesa::Vector3 phi{0., 0., 0.4};
  const auto rotation = fesa::ExpSkew(phi);
  for (std::size_t node = 0; node < 4; ++node) {
    const fesa::Vector3 x{nodes[node].coordinates};
    trial.rotation_vector[node] = phi;
    fesa::Vector3 rotated{
        rotation(0, 0) * x[0] + rotation(0, 1) * x[1] + rotation(0, 2) * x[2],
        rotation(1, 0) * x[0] + rotation(1, 1) * x[1] + rotation(1, 2) * x[2],
        rotation(2, 0) * x[0] + rotation(2, 1) * x[1] + rotation(2, 2) * x[2]};
    trial.displacement[node] = rotated - x;
  }
  const auto actual = fesa::EvaluateInternalResidualAndTangent(shell, trial);
  ASSERT_TRUE(actual.HasValue());
  EXPECT_NEAR(actual.Value().physical_strain_energy, 0., 1e-25);
  EXPECT_NEAR(actual.Value().drilling_energy, .5 * kd * 4 * .4 * .4, 1e-14);
  for (std::size_t row = 0; row < 24; ++row) {
    const double force = row % 6 == 5 ? kd * .4 : 0.;
    EXPECT_NEAR(actual.Value().drilling_force[row], force, 1e-14);
    for (std::size_t col = 0; col < 24; ++col) {
      const double tangent = row == col && row % 6 == 5 ? kd : 0.;
      EXPECT_NEAR(actual.Value().drilling_tangent(row, col), tangent, 1e-14);
    }
  }
}

TEST(Mitc4DrillingSelection,
     DirectSpringUsesEachInitialDirectorAndConstantBlocks) {
  for (const bool warped : {false, true}) {
    auto nodes = Nodes();
    if (warped) nodes[2].coordinates[2] = .1;
    auto frames = Frames();
    const auto transform = [](const fesa::Matrix& r, const Coordinates& x) {
      Coordinates y{};
      for (std::size_t i = 0; i < 3; ++i)
        for (std::size_t j = 0; j < 3; ++j) y[i] += r(i, j) * x[j];
      return y;
    };
    const auto global = fesa::ExpSkew({.4, -.3, .2});
    std::array<Coordinates, 4> directors{};
    for (std::size_t n = 0; n < 4; ++n) {
      const auto r =
          fesa::ExpSkew({0., warped ? .03 * static_cast<double>(n) : 0., 0.});
      frames[n].director = transform(global, transform(r, frames[n].director));
      frames[n].tangent_a =
          transform(global, transform(r, frames[n].tangent_a));
      frames[n].tangent_b =
          transform(global, transform(r, frames[n].tangent_b));
      nodes[n].coordinates = transform(global, nodes[n].coordinates);
      directors[n] = frames[n].director;
    }
    const fesa::ShellSection section{"Section", .2, 0, {"spring.inp", 1}};
    const fesa::LinearElasticMaterial material{
        "Material", 120., .25, {"spring.inp", 2}};
    const auto shell = fesa::NonlinearMitc4Shell::Create(
        NodePointers(nodes), frames, section, material);
    const auto linear = fesa::Mitc4Shell::Create(NodePointers(nodes), directors,
                                                 section, material);
    ASSERT_TRUE(shell.HasValue());
    ASSERT_TRUE(linear.HasValue());
    const auto k = linear.Value().Stiffness();
    ASSERT_TRUE(k.HasValue());
    double trace = 0.;
    for (std::size_t n = 0; n < 4; ++n)
      for (std::size_t c = 3; c < 5; ++c)
        trace += k.Value().physical_local20(5 * n + c, 5 * n + c);
    const double kd = 1e-3 * trace / 8.;
    auto trial = Trial();
    trial.rotation_vector[0] = {.65, -.4, .3};
    const auto direction = ToVector(Direction());
    const auto response = fesa::EvaluateInternalResidualAndTangent(
        shell.Value(), trial, &direction);
    const auto shifted = fesa::EvaluateInternalResidualAndTangent(
        shell.Value(), AddScaled(trial, Direction(), .2), &direction);
    ASSERT_TRUE(response.HasValue());
    ASSERT_TRUE(shifted.HasValue());
    double energy = 0.;
    for (std::size_t n = 0; n < 4; ++n) {
      const double gamma =
          fesa::Vector3{directors[n]}.Dot(trial.rotation_vector[n]);
      energy += .5 * kd * gamma * gamma;
      const fesa::Vector3 d{directors[n]};
      const auto rotation = fesa::EvaluateDynamicRotation(
          trial.rotation_vector[n], d, Direction().rotation_vector[n], {});
      ASSERT_TRUE(rotation.HasValue());
      EXPECT_NEAR(response.Value().drilling_constraint[n],
                  kd * gamma * rotation.Value().drill_null.Dot(d), 1e-14);
      for (std::size_t c = 0; c < 3; ++c) {
        double ds = 0.;
        for (std::size_t j = 0; j < 3; ++j)
          ds += d[j] * rotation.Value().null_gradient(j, c);
        EXPECT_NEAR(
            response.Value().drilling_constraint_gradient(n, 6 * n + 3 + c),
            kd * (d[c] * rotation.Value().drill_null.Dot(d) + gamma * ds),
            1e-14);
      }
      for (std::size_t c = 0; c < 3; ++c)
        EXPECT_NEAR(response.Value().drilling_force[6 * n + 3 + c],
                    kd * gamma * directors[n][c], 1e-14);
    }
    EXPECT_NEAR(response.Value().drilling_energy, energy, 1e-14);
    for (std::size_t i = 0; i < 24; ++i) {
      EXPECT_DOUBLE_EQ(response.Value().drilling_force_directional_second[i],
                       0.);
      if (i % 6 < 3) EXPECT_DOUBLE_EQ(response.Value().drilling_force[i], 0.);
      for (std::size_t j = 0; j < 24; ++j) {
        const double expected =
            i / 6 == j / 6 && i % 6 >= 3 && j % 6 >= 3
                ? kd * directors[i / 6][i % 6 - 3] * directors[i / 6][j % 6 - 3]
                : 0.;
        EXPECT_NEAR(response.Value().drilling_tangent(i, j), expected, 1e-14);
        EXPECT_NEAR(response.Value().drilling_metric(i, j), expected, 1e-14);
        EXPECT_DOUBLE_EQ(response.Value().drilling_tangent(i, j),
                         shifted.Value().drilling_tangent(i, j));
        if (expected == 0.)
          EXPECT_DOUBLE_EQ(response.Value().drilling_tangent(i, j), 0.);
      }
    }
  }
}

TEST(Mitc4DrillingSelection,
     DirectConstraintUsesPhysicalNullAndItsHighDerivatives) {
  const auto nodes = Nodes();
  const auto shell = ProductionShell(nodes);
  auto trial = Trial();
  trial.rotation_vector[0] = {.65, -.4, .3};
  const auto direction = Direction();
  const auto velocity = ToVector(direction);
  const auto response =
      fesa::EvaluateInternalResidualAndTangent(shell, trial, &velocity);
  ASSERT_TRUE(response.HasValue());
  const auto zero =
      fesa::EvaluateInternalResidualAndTangent(shell, {}, &velocity);
  ASSERT_TRUE(zero.HasValue());
  const double kd = zero.Value().drilling_tangent(5, 5);
  // Independent inertia API provides n, Dn and D2n, never a twist angle.
  const fesa::Vector3 d{0., 0., 1.};
  const auto scalar_gradient = [&](const fesa::Vector3& phi, std::size_t c) {
    const auto r = fesa::EvaluateDynamicRotation(phi, d, {}, {});
    EXPECT_TRUE(r.HasValue());
    double projection_derivative = 0.;
    for (std::size_t k = 0; k < 3; ++k)
      projection_derivative += r.Value().null_gradient(k, c) * d[k];
    return kd * (d[c] * r.Value().drill_null.Dot(d) +
                 d.Dot(phi) * projection_derivative);
  };
  constexpr double h = .002;
  for (std::size_t n = 0; n < 4; ++n) {
    const auto phi = trial.rotation_vector[n];
    const auto v = direction.rotation_vector[n];
    const auto r = fesa::EvaluateDynamicRotation(phi, d, v, {});
    ASSERT_TRUE(r.HasValue());
    const double gamma = d.Dot(phi);
    const double gamma_dot = d.Dot(v);
    double projection_dot = 0.;
    for (std::size_t i = 0; i < 3; ++i)
      for (std::size_t j = 0; j < 3; ++j)
        projection_dot += d[i] * r.Value().null_gradient(i, j) * v[j];
    const double expected_second =
        kd * (2. * gamma_dot * projection_dot +
              gamma * r.Value().null_directional_second.Dot(d));
    double actual_second = 0.;
    EXPECT_NEAR(response.Value().drilling_constraint[n],
                kd * gamma * r.Value().drill_null.Dot(d), 1e-14);
    for (std::size_t i = 0; i < 24; ++i) {
      const double gradient =
          i / 6 == n && i % 6 >= 3 ? scalar_gradient(phi, i % 6 - 3) : 0.;
      EXPECT_NEAR(response.Value().drilling_constraint_gradient(n, i), gradient,
                  1e-14);
      for (std::size_t j = 0; j < 24; ++j)
        actual_second += velocity[i] *
                         response.Value().drilling_constraint_hessian[n](i, j) *
                         velocity[j];
      if (i / 6 == n && i % 6 >= 3) {
        const auto c = i % 6 - 3;
        const double fd =
            (scalar_gradient(phi + h * v, c) - 2. * scalar_gradient(phi, c) +
             scalar_gradient(phi - h * v, c)) /
            (h * h);
        const double actual =
            response.Value()
                .drilling_constraint_third_directional_gradient[n][i];
        EXPECT_NEAR(actual, fd,
                    2e-5 * (std::abs(actual) + std::abs(fd) + kd * 1e-4));
      } else {
        EXPECT_DOUBLE_EQ(
            response.Value()
                .drilling_constraint_third_directional_gradient[n][i],
            0.);
      }
    }
    EXPECT_NEAR(actual_second, expected_second, 1e-14);
  }
}

TEST(Mitc4DrillingSelection, DirectExternalConstraintUsesActualInertiaNull) {
  const auto nodes = Nodes();
  const auto shell = ProductionShell(nodes);
  const fesa::Vector3 phi{.65, -.4, .3}, velocity{.2, .1, -.15};
  const fesa::Vector3 force{.7, -.3, .2}, rate{-.1, .2, .4},
      acceleration{.5, -.2, .1};
  const auto actual = shell.EvaluateExternalDrillingConstraint(
      0, phi, velocity, force, rate, acceleration);
  const auto rotation =
      fesa::EvaluateDynamicRotation(phi, {0., 0., 1.}, velocity, {});
  ASSERT_TRUE(actual.HasValue());
  ASSERT_TRUE(rotation.HasValue());
  const double h = .002;
  const auto plus = shell.EvaluateExternalDrillingConstraint(
      0, phi + h * velocity, velocity,
      force + h * rate + .5 * h * h * acceleration, rate + h * acceleration,
      acceleration);
  const auto minus = shell.EvaluateExternalDrillingConstraint(
      0, phi - h * velocity, velocity,
      force - h * rate + .5 * h * h * acceleration, rate - h * acceleration,
      acceleration);
  ASSERT_TRUE(plus.HasValue());
  ASSERT_TRUE(minus.HasValue());
  for (std::size_t c = 0; c < 3; ++c) {
    EXPECT_NEAR(
        actual.Value().gradient_directional_first[c],
        (plus.Value().gradient[c] - minus.Value().gradient[c]) / (2. * h),
        1e-7);
    EXPECT_NEAR(actual.Value().gradient_directional_second[c],
                (plus.Value().gradient[c] - 2. * actual.Value().gradient[c] +
                 minus.Value().gradient[c]) /
                    (h * h),
                2e-7);
    std::array<double, 3> basis{};
    basis[c] = 1.;
    constexpr double coordinate_h = 1e-4;
    const auto positive = shell.EvaluateExternalDrillingConstraint(
        0, phi + coordinate_h * fesa::Vector3{basis}, velocity, force, rate,
        acceleration);
    const auto negative = shell.EvaluateExternalDrillingConstraint(
        0, phi - coordinate_h * fesa::Vector3{basis}, velocity, force, rate,
        acceleration);
    ASSERT_TRUE(positive.HasValue());
    ASSERT_TRUE(negative.HasValue());
    for (std::size_t j = 0; j < 3; ++j)
      EXPECT_NEAR(
          actual.Value().hessian(j, c),
          (positive.Value().gradient[j] - negative.Value().gradient[j]) /
              (2. * coordinate_h),
          1e-7);
  }
  const auto& n = rotation.Value();
  EXPECT_NEAR(actual.Value().value, -n.drill_null.Dot(force), 1e-14);
  fesa::Vector3 n_dot{};
  std::array<double, 3> values{};
  for (std::size_t i = 0; i < 3; ++i)
    for (std::size_t j = 0; j < 3; ++j)
      values[i] += n.null_gradient(i, j) * velocity[j];
  n_dot = fesa::Vector3{values};
  EXPECT_NEAR(actual.Value().directional_first,
              -n_dot.Dot(force) - n.drill_null.Dot(rate), 1e-14);
  EXPECT_NEAR(actual.Value().directional_second,
              -n.null_directional_second.Dot(force) - 2. * n_dot.Dot(rate) -
                  n.drill_null.Dot(acceleration),
              1e-14);
  for (std::size_t c = 0; c < 3; ++c) {
    double expected = 0.;
    for (std::size_t i = 0; i < 3; ++i)
      expected -= n.null_gradient(i, c) * force[i];
    EXPECT_NEAR(actual.Value().gradient[c], expected, 1e-14);
  }
}

TEST(Mitc4DrillingSelection,
     ModesAreImmutableReferenceInputsAndPreservePhysicalResponse) {
  const auto nodes = Nodes();
  const auto direct = ProductionShell(nodes);
  const auto legacy = LegacyShell(nodes);
  EXPECT_EQ(direct.DrillingMode(), fesa::Mitc4DrillingMode::kDirectSpring);
  EXPECT_EQ(legacy.DrillingMode(),
            fesa::Mitc4DrillingMode::kLegacyCorotatedTwist);
  EXPECT_STREQ(fesa::Mitc4DrillingModeName(direct.DrillingMode()),
               "direct-initial-director-spring");
  EXPECT_STREQ(fesa::Mitc4DrillingModeName(legacy.DrillingMode()),
               "legacy-corotated-twist");
  const auto a =
      fesa::nonlinear_mitc4_internal::ConfigurationSnapshot::Read(direct);
  const auto b =
      fesa::nonlinear_mitc4_internal::ConfigurationSnapshot::Read(legacy);
  EXPECT_EQ(a.values, b.values);
  EXPECT_NE(a.drilling_mode, b.drilling_mode);
  const auto direction = ToVector(Direction());
  const auto r =
      fesa::EvaluateInternalResidualAndTangent(direct, Trial(), &direction);
  const auto l =
      fesa::EvaluateInternalResidualAndTangent(legacy, Trial(), &direction);
  ASSERT_TRUE(r.HasValue());
  ASSERT_TRUE(l.HasValue());
  EXPECT_DOUBLE_EQ(r.Value().physical_strain_energy,
                   l.Value().physical_strain_energy);
  EXPECT_NE(r.Value().drilling_energy, l.Value().drilling_energy);
  for (std::size_t i = 0; i < 24; ++i) {
    EXPECT_DOUBLE_EQ(r.Value().physical_force[i], l.Value().physical_force[i]);
    for (std::size_t j = 0; j < 24; ++j)
      EXPECT_NEAR(
          r.Value().internal_tangent(i, j) - r.Value().drilling_tangent(i, j),
          l.Value().internal_tangent(i, j) - l.Value().drilling_tangent(i, j),
          1e-14);
  }
  for (std::size_t gp = 0; gp < 4; ++gp)
    EXPECT_EQ(r.Value().recovery.gauss_points[gp].section_resultant,
              l.Value().recovery.gauss_points[gp].section_resultant);
  const fesa::ShellSection section{"Section", .2, 0, {"spring.inp", 1}};
  const fesa::LinearElasticMaterial material{
      "Material", 120., .25, {"spring.inp", 2}};
  const auto linear = fesa::Mitc4Shell::Create(
      NodePointers(nodes),
      {{{0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}}}, section,
      material);
  ASSERT_TRUE(linear.HasValue());
  EXPECT_EQ(linear.Value().DrillingMode(), direct.DrillingMode());
  const auto alternate = fesa::mitc4_drilling_internal::TestAccess::WithMode(
      linear.Value(), fesa::Mitc4DrillingMode::kLegacyCorotatedTwist);
  const auto k1 = linear.Value().Stiffness();
  const auto k2 = alternate.Stiffness();
  ASSERT_TRUE(k1.HasValue());
  ASSERT_TRUE(k2.HasValue());
  for (std::size_t i = 0; i < 24; ++i)
    for (std::size_t j = 0; j < 24; ++j)
      EXPECT_DOUBLE_EQ(k1.Value().stabilized_global24(i, j),
                       k2.Value().stabilized_global24(i, j));
}

TEST(Mitc4DrillingSelection, DirectRejectsInvalidScaleAndNonfiniteTrial) {
  const auto nodes = Nodes();
  const fesa::ShellSection section{"Section", .2, 0, {"spring.inp", 1}};
  for (double youngs_modulus : {1e-320, (std::numeric_limits<double>::max)()}) {
    const fesa::LinearElasticMaterial material{
        "Material", youngs_modulus, .25, {"spring.inp", 2}};
    const auto invalid = fesa::NonlinearMitc4Shell::Create(
        NodePointers(nodes), Frames(), section, material);
    ASSERT_FALSE(invalid.HasValue());
    EXPECT_EQ(invalid.GetStatus().Diagnostics()[0].code,
              youngs_modulus < 1.0 ? "invalid-nonlinear-shell-drilling"
                                  : "invalid-shell-material");
  }
  const auto shell = ProductionShell(nodes);
  for (double value : {std::numeric_limits<double>::quiet_NaN(),
                       std::numeric_limits<double>::infinity()}) {
    auto trial = Trial();
    trial.rotation_vector[0] = {value, 0., 0.};
    EXPECT_FALSE(
        fesa::EvaluateInternalResidualAndTangent(shell, trial).HasValue());
  }
  const auto valid = fesa::EvaluateInternalResidualAndTangent(shell, {});
  ASSERT_TRUE(valid.HasValue());
  EXPECT_DOUBLE_EQ(valid.Value().drilling_energy, 0.);
}

TEST(Mitc4DrillingSelection, DirectRetainsPositiveSubnormalSpringScale) {
  const auto nodes = Nodes();
  const fesa::ShellSection section{"Section", .2, 0, {"spring.inp", 1}};
  unsigned accepted = 0;
  for (double youngs_modulus :
       {1e-315, 1e-316, 1e-317, 1e-318, 2e-318, 3e-318, 1e-319}) {
    const fesa::LinearElasticMaterial material{
        "Material", youngs_modulus, .25, {"spring.inp", 2}};
    const auto linear = fesa::Mitc4Shell::Create(
        NodePointers(nodes),
        {{{0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}}}, section,
        material);
    ASSERT_TRUE(linear.HasValue());
    const auto stiffness = linear.Value().Stiffness();
    const auto shell = fesa::NonlinearMitc4Shell::Create(
        NodePointers(nodes), Frames(), section, material);
    if (!stiffness.HasValue()) {
      EXPECT_FALSE(shell.HasValue());
      continue;
    }
    ASSERT_TRUE(shell.HasValue());
    ++accepted;
    fesa::NonlinearMitc4TrialDofs trial{};
    trial.rotation_vector[0] = {0., 0., 1.};
    const auto actual =
        fesa::EvaluateInternalResidualAndTangent(shell.Value(), trial);
    ASSERT_TRUE(actual.HasValue());
    // Scale is already validated. No later half-scale rounding may change K.
    const double kd = stiffness.Value().drilling_stiffness;
    EXPECT_EQ(actual.Value().drilling_tangent(5, 5), kd);
    EXPECT_EQ(actual.Value().drilling_force[5], kd);
    trial.rotation_vector[0] = {0., 0., 2.};
    const auto double_angle =
        fesa::EvaluateInternalResidualAndTangent(shell.Value(), trial);
    ASSERT_TRUE(double_angle.HasValue());
    // U(2) = 2*kd is representable even for the least positive spring scale.
    EXPECT_EQ(double_angle.Value().drilling_energy, 2. * kd);
    EXPECT_EQ(double_angle.Value().drilling_force[5], 2. * kd);
  }
  EXPECT_GT(accepted, 0U);
}

TEST(Mitc4DrillingSelection, LegacyHeadFrozenEvidence) {
  const auto nodes = Nodes();
  const auto shell = LegacyShell(nodes);
  const auto direction = ToVector(Direction());
  const auto response =
      fesa::EvaluateInternalResidualAndTangent(shell, Trial(), &direction);
  ASSERT_TRUE(response.HasValue());
  const auto& r = response.Value();
  std::array<std::uint64_t, 8> hashes{};
  hashes.fill(14695981039346656037ULL);
  const auto add = [&](std::size_t channel, double value) {
    std::uint64_t bits{};
    std::memcpy(&bits, &value, sizeof(bits));
    for (unsigned b = 0; b < 8; ++b) {
      hashes[channel] ^= (bits >> (8 * b)) & 255ULL;
      hashes[channel] *= 1099511628211ULL;
    }
  };
  add(0, r.drilling_energy);
  for (std::size_t i = 0; i < 24; ++i) {
    add(1, r.drilling_force[i]);
    add(4, r.force_directional_second[i]);
    for (std::size_t j = 0; j < 24; ++j) {
      add(2, r.drilling_tangent(i, j));
      add(3, r.drilling_metric(i, j));
    }
  }
  for (std::size_t n = 0; n < 4; ++n) {
    add(5, r.drilling_constraint[n]);
    for (std::size_t i = 0; i < 24; ++i) {
      add(5, r.drilling_constraint_gradient(n, i));
      add(7, r.drilling_constraint_third_directional_gradient[n][i]);
      for (std::size_t j = 0; j < 24; ++j)
        add(6, r.drilling_constraint_hessian[n](i, j));
    }
  }
  // Captured before implementation from HEAD
  // 52e321e15c70377d4318aace3f3bc64b3ea1a65c, native MSVC x64 Debug. Entire
  // energy/force/tangent/metric/directional/psi inventories are pinned
  // independently, including the historical polar1 constraint inconsistency. Do
  // not regenerate these from the selected mode.
  constexpr std::array<std::uint64_t, 8> frozen{
      16395275265220000061ULL, 12405324468590059764ULL, 79034762627498281ULL,
      7305417493946306304ULL,  7191018337857616990ULL,  4845391414006668707ULL,
      11613163243225266145ULL, 3872161451216566265ULL};
  EXPECT_EQ(hashes, frozen);
}

}  // namespace