#include <gtest/gtest.h>

#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <memory>
#include <optional>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>

#include "analysis/hht_dynamic_internal.h"
#include "fesa/analysis/analysis_model.h"
#include "fesa/assembly/parallel_for.h"
#include "fesa/elements/finite_rotation.h"
#include "fesa/solvers/linear/mkl_pardiso_solver.h"
#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h"

namespace {

using Frame = fesa::DynamicAcceptedFrameState;
using Assembly = fesa::hht_internal::DynamicAssembly;
using L3 = std::array<long double, 3>;
constexpr double kEndTime = 0.25;

std::string Message(const fesa::Status& status) {
  return status.Diagnostics().empty()
             ? "Unspecified solver failure"
             : status.Diagnostics().front().code + ": " +
                   status.Diagnostics().front().message;
}

void Require(const fesa::Status& status) {
  if (!status.IsOk()) throw std::runtime_error{Message(status)};
}

template <typename T>
T Require(fesa::Result<T> result) {
  if (!result.HasValue()) throw std::runtime_error{Message(result.GetStatus())};
  return std::move(result.Value());
}

L3 Add(const L3& a, const L3& b) {
  return {a[0] + b[0], a[1] + b[1], a[2] + b[2]};
}
L3 Scale(const L3& a, const long double factor) {
  return {factor * a[0], factor * a[1], factor * a[2]};
}
L3 Cross(const L3& a, const L3& b) {
  return {a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2],
          a[0] * b[1] - a[1] * b[0]};
}
long double Dot(const L3& a, const L3& b) {
  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
long double Norm(const L3& a) {
  return std::hypot(a[0], std::hypot(a[1], a[2]));
}
L3 Triple(const fesa::Vector& vector, const std::size_t offset) {
  return {vector[offset], vector[offset + 1U], vector[offset + 2U]};
}
L3 Coordinates(const fesa::Node& node) {
  return {node.coordinates[0], node.coordinates[1], node.coordinates[2]};
}
void Store(const L3& value, fesa::Vector* vector, const std::size_t offset) {
  for (std::size_t c = 0; c < 3; ++c)
    (*vector)[offset + c] = static_cast<double>(value[c]);
}

// Independent subset-block exponential, with no Rodrigues/Jl scalar formulas.
// The final block gives the mixed derivative in the supplied chart directions.
L3 ExponentialDerivative(const L3& phi, const L3& initial,
                         const std::vector<L3>& directions = {}) {
  std::array<L3, 16> term{};
  std::array<L3, 16> sum{};
  term[0] = sum[0] = initial;
  const std::size_t blocks = std::size_t{1} << directions.size();
  for (std::size_t order = 1U; order <= 64U; ++order) {
    std::array<L3, 16> next{};
    for (std::size_t subset = 0; subset < blocks; ++subset) {
      next[subset] = Cross(phi, term[subset]);
      for (std::size_t direction = 0; direction < directions.size();
           ++direction) {
        const std::size_t bit = std::size_t{1} << direction;
        if ((subset & bit) != 0U)
          next[subset] = Add(next[subset],
                             Cross(directions[direction], term[subset ^ bit]));
      }
      next[subset] =
          Scale(next[subset], 1.0L / static_cast<long double>(order));
      sum[subset] = Add(sum[subset], next[subset]);
    }
    term = next;
  }
  return sum[blocks - 1U];
}

struct Fixture {
  std::unique_ptr<fesa::Domain> domain;
  fesa::DofManager dofs;
  std::vector<fesa::hht_internal::RuntimeShell> shells;
  double length;
};

Fixture MakeFixture(const std::size_t divisions) {
  fesa::ModelDefinition definition{};
  definition.source_path = "hht-shell-physics.inp";
  definition.source_content_identity = "internal-C08-fixture";
  for (std::size_t y = 0; y <= divisions; ++y) {
    for (std::size_t x = 0; x <= divisions; ++x) {
      const auto index = definition.nodes.size();
      const auto label = static_cast<std::int64_t>(index + 1U);
      definition.nodes.push_back({{"Plate", label, std::to_string(label)},
                                  {-1.0 + 2.0 * static_cast<double>(x) /
                                              static_cast<double>(divisions),
                                   -1.0 + 2.0 * static_cast<double>(y) /
                                              static_cast<double>(divisions),
                                   0.0},
                                  {}});
      definition.shell_node_initial_frames.push_back(
          {static_cast<fesa::EntityIndex>(index),
           {0.0, 0.0, 1.0},
           {1.0, 0.0, 0.0},
           {0.0, 1.0, 0.0}});
    }
  }
  definition.materials = {{"Elastic", 1000.0, 0.3, {}}};
  definition.densities = {{0U, 1.0, {}}};
  definition.shell_sections = {{"Centered", 0.2, 0U, {}}};
  for (std::size_t y = 0; y < divisions; ++y) {
    for (std::size_t x = 0; x < divisions; ++x) {
      const auto first =
          static_cast<fesa::EntityIndex>(y * (divisions + 1U) + x);
      const auto stride = static_cast<fesa::EntityIndex>(divisions + 1U);
      const auto label =
          static_cast<std::int64_t>(definition.shell_elements.size() + 1U);
      definition.shell_elements.push_back(
          {{"Plate", label, std::to_string(label)},
           fesa::ShellSourceElementType::kS4,
           {first, first + 1U, first + stride + 1U, first + stride},
           0U,
           0U,
           {}});
    }
  }
  fesa::DynamicStepDefinition step{};
  step.source_step_name_raw = "C08";
  step.source_step_name_canonical = "C08";
  step.initial_increment = 1.0 / 64.0;
  step.time_period = kEndTime;
  step.increment_count = 16U;
  definition.dynamic_step = step;
  auto domain = std::make_unique<fesa::Domain>(
      Require(fesa::Domain::Create(std::move(definition))));
  auto model = Require(fesa::AnalysisModel::Create(*domain));
  auto dofs = Require(fesa::DofManager::Create(model));
  auto shells = Require(fesa::hht_internal::BuildRuntimeShells(*domain, dofs));
  const double length = Require(fesa::hht_internal::ReferenceLength(*domain));
  return {std::move(domain), std::move(dofs), std::move(shells), length};
}

Frame EmptyFrame(const Fixture& fixture) {
  Frame frame{};
  frame.is_initial = true;
  for (auto* vector :
       {&frame.generalized_coordinates, &frame.velocity, &frame.acceleration,
        &frame.algorithmic_acceleration, &frame.external_force}) {
    *vector = fesa::Vector(fixture.dofs.FullDofCount());
  }
  return frame;
}

Assembly Evaluate(const Fixture& fixture, Frame frame, const double h,
                  const fesa::HhtParameters& parameters, const bool third) {
  const fesa::SerialParallelFor parallel;
  fesa::Vector load_scale(frame.external_force.Size());
  for (std::size_t i = 0; i < load_scale.Size(); ++i)
    load_scale[i] = std::abs(frame.external_force[i]);
  return Require(fesa::hht_internal::EvaluateAssembly(
      *fixture.domain, fixture.dofs, fixture.shells, std::move(frame),
      load_scale, h, parameters, parallel, third));
}

struct VolumeOracle {
  fesa::Vector force;
  long double kinetic = 0.0L;
  long double power = 0.0L;
  long double power_scale = 0.0L;
  L3 momentum{};
  L3 angular_momentum{};
  L3 momentum_rate{};
  L3 angular_momentum_rate{};
  long double momentum_scale = 0.0L;
  long double angular_momentum_scale = 0.0L;
  long double momentum_rate_scale = 0.0L;
  long double angular_momentum_rate_scale = 0.0L;
};

// Independent 5^3 rectangle quadrature and derivative-exponential kinematics.
// No production mass, inertia force, momentum, or director derivative is
// reused.
VolumeOracle Integrate(const Fixture& fixture, const Frame& frame) {
  constexpr std::array<long double, 5> points{
      -0.906179845938663992797626878299L, -0.538469310105683091036314420700L,
      0.0L, 0.538469310105683091036314420700L,
      0.906179845938663992797626878299L};
  constexpr std::array<long double, 5> weights{
      0.236926885056189087514264040720L, 0.478628670499366468041291514836L,
      0.568888888888888888888888888889L, 0.478628670499366468041291514836L,
      0.236926885056189087514264040720L};
  constexpr std::array<std::array<int, 2>, 4> signs{
      {{-1, -1}, {1, -1}, {1, 1}, {-1, 1}}};
  VolumeOracle oracle{fesa::Vector(fixture.dofs.FullDofCount())};
  std::vector<long double> force(fixture.dofs.FullDofCount(), 0.0L);
  for (std::size_t element = 0;
       element < fixture.domain->ShellElements().Size(); ++element) {
    const auto& shell = fixture.domain->ShellElements()[element];
    std::array<L3, 4> x{};
    std::array<L3, 4> v{};
    std::array<L3, 4> a{};
    std::array<L3, 4> director{};
    std::array<L3, 4> director_rate{};
    std::array<L3, 4> director_acceleration{};
    std::array<std::array<L3, 3>, 4> derivative{};
    for (std::size_t node = 0; node < 4; ++node) {
      const auto global = shell.node_indices[node];
      const auto offset = 6U * global;
      const L3 phi = Triple(frame.generalized_coordinates, offset + 3U);
      const L3 phidot = Triple(frame.velocity, offset + 3U);
      const L3 phiddot = Triple(frame.acceleration, offset + 3U);
      x[node] = Add(Coordinates(fixture.domain->Nodes()[global]),
                    Triple(frame.generalized_coordinates, offset));
      v[node] = Triple(frame.velocity, offset);
      a[node] = Triple(frame.acceleration, offset);
      director[node] = ExponentialDerivative(phi, {0.0L, 0.0L, 1.0L});
      director_rate[node] =
          ExponentialDerivative(phi, {0.0L, 0.0L, 1.0L}, {phidot});
      director_acceleration[node] =
          Add(ExponentialDerivative(phi, {0.0L, 0.0L, 1.0L}, {phidot, phidot}),
              ExponentialDerivative(phi, {0.0L, 0.0L, 1.0L}, {phiddot}));
      for (std::size_t c = 0; c < 3; ++c) {
        L3 axis{};
        axis[c] = 1.0L;
        derivative[node][c] =
            ExponentialDerivative(phi, {0.0L, 0.0L, 1.0L}, {axis});
      }
    }
    const L3 edge_x =
        Add(Coordinates(fixture.domain->Nodes()[shell.node_indices[1]]),
            Scale(Coordinates(fixture.domain->Nodes()[shell.node_indices[0]]),
                  -1.0L));
    const L3 edge_y =
        Add(Coordinates(fixture.domain->Nodes()[shell.node_indices[3]]),
            Scale(Coordinates(fixture.domain->Nodes()[shell.node_indices[0]]),
                  -1.0L));
    const long double jacobian = Norm(Cross(edge_x, edge_y)) * 0.025L;
    for (std::size_t i = 0; i < 5; ++i) {
      for (std::size_t j = 0; j < 5; ++j) {
        std::array<long double, 4> shape{};
        for (std::size_t node = 0; node < 4; ++node)
          shape[node] = 0.25L * (1.0L + signs[node][0] * points[i]) *
                        (1.0L + signs[node][1] * points[j]);
        for (std::size_t k = 0; k < 5; ++k) {
          const long double z = 0.1L * points[k];
          const long double weight =
              jacobian * weights[i] * weights[j] * weights[k];
          L3 position{};
          L3 velocity{};
          L3 acceleration{};
          for (std::size_t node = 0; node < 4; ++node) {
            position =
                Add(position,
                    Scale(Add(x[node], Scale(director[node], z)), shape[node]));
            velocity =
                Add(velocity, Scale(Add(v[node], Scale(director_rate[node], z)),
                                    shape[node]));
            acceleration =
                Add(acceleration,
                    Scale(Add(a[node], Scale(director_acceleration[node], z)),
                          shape[node]));
          }
          oracle.kinetic += 0.5L * weight * Dot(velocity, velocity);
          oracle.power += weight * Dot(velocity, acceleration);
          oracle.power_scale += weight * Norm(velocity) * Norm(acceleration);
          oracle.momentum = Add(oracle.momentum, Scale(velocity, weight));
          oracle.angular_momentum =
              Add(oracle.angular_momentum,
                  Scale(Cross(position, velocity), weight));
          oracle.momentum_rate =
              Add(oracle.momentum_rate, Scale(acceleration, weight));
          oracle.angular_momentum_rate =
              Add(oracle.angular_momentum_rate,
                  Scale(Cross(position, acceleration), weight));
          oracle.momentum_scale += weight * Norm(velocity);
          oracle.angular_momentum_scale +=
              weight * Norm(Cross(position, velocity));
          oracle.momentum_rate_scale += weight * Norm(acceleration);
          oracle.angular_momentum_rate_scale +=
              weight * Norm(Cross(position, acceleration));
          for (std::size_t node = 0; node < 4; ++node) {
            const auto offset = 6U * shell.node_indices[node];
            for (std::size_t c = 0; c < 3; ++c) {
              force[offset + c] += weight * shape[node] * acceleration[c];
              force[offset + 3U + c] += weight * shape[node] * z *
                                        Dot(derivative[node][c], acceleration);
            }
          }
        }
      }
    }
  }
  for (std::size_t i = 0; i < force.size(); ++i)
    oracle.force[i] = static_cast<double>(force[i]);
  return oracle;
}

std::array<L3, 5> OrbitChart(const double time) {
  const long double t = time;
  return {{{0.4L * std::sin(t), 0.6L * std::sin(0.7L * t),
            0.2L * std::sin(1.3L * t)},
           {0.4L * std::cos(t), 0.42L * std::cos(0.7L * t),
            0.26L * std::cos(1.3L * t)},
           {-0.4L * std::sin(t), -0.294L * std::sin(0.7L * t),
            -0.338L * std::sin(1.3L * t)},
           {-0.4L * std::cos(t), -0.2058L * std::cos(0.7L * t),
            -0.4394L * std::cos(1.3L * t)},
           {0.4L * std::sin(t), 0.14406L * std::sin(0.7L * t),
            0.57122L * std::sin(1.3L * t)}}};
}

Frame RigidOrbit(const Fixture& fixture, const double time) {
  Frame frame = EmptyFrame(fixture);
  frame.time = time;
  frame.is_initial = time == 0.0;
  const auto chart = OrbitChart(time);
  const auto& phi = chart[0];
  const auto& velocity = chart[1];
  const auto& acceleration = chart[2];
  for (std::size_t node = 0; node < fixture.domain->Nodes().size(); ++node) {
    const auto offset = 6U * node;
    const L3 initial = Coordinates(fixture.domain->Nodes()[node]);
    Store(Add(ExponentialDerivative(phi, initial), Scale(initial, -1.0L)),
          &frame.generalized_coordinates, offset);
    Store(phi, &frame.generalized_coordinates, offset + 3U);
    Store(ExponentialDerivative(phi, initial, {velocity}), &frame.velocity,
          offset);
    Store(velocity, &frame.velocity, offset + 3U);
    Store(Add(ExponentialDerivative(phi, initial, {velocity, velocity}),
              ExponentialDerivative(phi, initial, {acceleration})),
          &frame.acceleration, offset);
    Store(acceleration, &frame.acceleration, offset + 3U);
  }
  frame.algorithmic_acceleration = frame.acceleration;
  return frame;
}

std::array<L3, 3> OrbitPhysicalRates(const std::array<L3, 5>& chart,
                                     const L3& initial) {
  const auto& q = chart[0];
  const auto& v = chart[1];
  const auto& a = chart[2];
  const auto& j = chart[3];
  const auto& s = chart[4];
  const auto d2 = Add(ExponentialDerivative(q, initial, {v, v}),
                      ExponentialDerivative(q, initial, {a}));
  const auto d3 =
      Add(Add(ExponentialDerivative(q, initial, {v, v, v}),
              Scale(ExponentialDerivative(q, initial, {v, a}), 3.0L)),
          ExponentialDerivative(q, initial, {j}));
  const auto d4 =
      Add(Add(ExponentialDerivative(q, initial, {v, v, v, v}),
              Scale(ExponentialDerivative(q, initial, {v, v, a}), 6.0L)),
          Add(Scale(ExponentialDerivative(q, initial, {a, a}), 3.0L),
              Add(Scale(ExponentialDerivative(q, initial, {v, j}), 4.0L),
                  ExponentialDerivative(q, initial, {s}))));
  return {d2, d3, d4};
}

// For the single centered rectangle, integral(N_i X)=X_i/3 and
// integral(N_i z^2)=t^3/12. Differentiate those independent volume moments.
std::array<fesa::Vector, 2> RigidForceRates(const Fixture& fixture,
                                            const double time) {
  std::array<fesa::Vector, 2> rates{fesa::Vector(24U), fesa::Vector(24U)};
  const auto chart = OrbitChart(time);
  const auto director = OrbitPhysicalRates(chart, {0.0L, 0.0L, 1.0L});
  for (std::size_t node = 0; node < 4U; ++node) {
    const auto x =
        OrbitPhysicalRates(chart, Coordinates(fixture.domain->Nodes()[node]));
    Store(Scale(x[1], 0.2L / 3.0L), &rates[0], 6U * node);
    Store(Scale(x[2], 0.2L / 3.0L), &rates[1], 6U * node);
    for (std::size_t c = 0; c < 3U; ++c) {
      L3 axis{};
      axis[c] = 1.0L;
      const auto g =
          ExponentialDerivative(chart[0], {0.0L, 0.0L, 1.0L}, {axis});
      const auto dg =
          ExponentialDerivative(chart[0], {0.0L, 0.0L, 1.0L}, {axis, chart[1]});
      const auto ddg = Add(ExponentialDerivative(chart[0], {0.0L, 0.0L, 1.0L},
                                                 {axis, chart[1], chart[1]}),
                           ExponentialDerivative(chart[0], {0.0L, 0.0L, 1.0L},
                                                 {axis, chart[2]}));
      rates[0][6U * node + 3U + c] = static_cast<double>(
          (0.008L / 12.0L) * (Dot(dg, director[0]) + Dot(g, director[1])));
      rates[1][6U * node + 3U + c] =
          static_cast<double>((0.008L / 12.0L) * (Dot(ddg, director[0]) +
                                                  2.0L * Dot(dg, director[1]) +
                                                  Dot(g, director[2])));
    }
  }
  return rates;
}

double SparseNorm(const fesa::SparseMatrix& matrix) {
  double norm = 0.0;
  for (const double value : matrix.Values()) norm = std::hypot(norm, value);
  return norm;
}

// The manufactured F(t) is prescribed in generalized coordinates. Its null
// work vanishes on the exact orbit, but not at a Newton trial. Include every
// explicit-time term of psi=n(q)^T(f_int(q)-F(t)); the public dead-load helper
// deliberately does not implement this verification-only load convention.
void ManufacturedHidden(const Assembly& assembly,
                        const std::array<fesa::Vector, 2>& force_rates,
                        Frame* frame) {
  const auto nodes = assembly.rotations.size();
  for (auto* vector :
       {&frame->hidden_constraint, &frame->hidden_constraint_scale,
        &frame->hidden_velocity_defect, &frame->hidden_velocity_defect_scale,
        &frame->hidden_acceleration_defect,
        &frame->hidden_acceleration_defect_scale})
    *vector = fesa::Vector(nodes);
  for (std::size_t node = 0; node < nodes; ++node) {
    const auto& rotation = assembly.rotations[node];
    double psi = 0.0;
    double psi_scale = 0.0;
    double first = 0.0;
    double first_scale = 0.0;
    double second = assembly.hidden_second[node];
    double second_scale = assembly.hidden_second_scale[node];
    auto add = [](const double term, double* value, double* scale) {
      *value += term;
      *scale += std::abs(term);
    };
    for (std::size_t c = 0; c < 3U; ++c) {
      const auto row = 6U * node + 3U + c;
      const double n = rotation.drill_null[c];
      const double g = frame->internal_force[row] - frame->external_force[row];
      psi += n * g;
      psi_scale += std::abs(n) * frame->residual_scale[row];
      for (std::size_t d = 0; d < 3U; ++d) {
        const auto column = 6U * node + 3U + d;
        const double dn = rotation.null_gradient(c, d);
        add(dn * g * frame->velocity[column], &first, &first_scale);
        add(dn * g * frame->acceleration[column], &second, &second_scale);
        add(-2.0 * dn * frame->velocity[column] * force_rates[0][row], &second,
            &second_scale);
      }
      const auto& stiffness = assembly.internal_tangent;
      for (std::size_t index = stiffness.RowOffsets()[row];
           index < stiffness.RowOffsets()[row + 1U]; ++index) {
        const auto column = stiffness.ColumnIndices()[index];
        const double gradient = n * stiffness.Values()[index];
        add(gradient * frame->velocity[column], &first, &first_scale);
        add(gradient * frame->acceleration[column], &second, &second_scale);
      }
      add(-n * force_rates[0][row], &first, &first_scale);
      add(-rotation.null_directional_second[c] * frame->external_force[row],
          &second, &second_scale);
      add(-n * force_rates[1][row], &second, &second_scale);
    }
    frame->hidden_constraint[node] = psi;
    frame->hidden_constraint_scale[node] = psi_scale;
    frame->hidden_velocity_defect[node] = first;
    frame->hidden_velocity_defect_scale[node] = first_scale;
    frame->hidden_acceleration_defect[node] = second;
    frame->hidden_acceleration_defect_scale[node] = second_scale;
  }
}

void CompleteManufacturedConstraints(
    const Fixture& fixture, const Assembly& assembly,
    const std::vector<std::size_t>& initial_null_nodes,
    fesa::MklPardisoSolver* positive, fesa::MklPardisoTangentSolver* tangent,
    Frame* frame) {
  // This read-only probe exercises current positive mass, null inventory,
  // psi and psi_q regularity. psi_q includes n'*(f_int-F), including moments.
  // Its autonomous hidden-rate outputs are discarded. Explicit time terms
  // and initial rate consistency are evaluated separately below.
  auto probe = assembly.frame;
  probe.is_initial = false;
  Require(fesa::hht_internal::CompleteConstraintEvidence(
      *fixture.domain, fixture.dofs, assembly, fixture.length,
      initial_null_nodes, *positive, *tangent, &probe));
  ManufacturedHidden(assembly, RigidForceRates(fixture, frame->time), frame);
  if (frame->is_initial) {
    for (std::size_t node = 0; node < assembly.rotations.size(); ++node) {
      for (const auto quantity :
           {std::make_pair(&frame->hidden_velocity_defect,
                           &frame->hidden_velocity_defect_scale),
            std::make_pair(&frame->hidden_acceleration_defect,
                           &frame->hidden_acceleration_defect_scale)}) {
        const double value = std::abs((*quantity.first)[node]);
        const double scale = (*quantity.second)[node];
        if (scale == 0.0)
          EXPECT_DOUBLE_EQ(value, 0.0);
        else
          EXPECT_LE(value / scale, 1.0e-10);
      }
    }
  }
}

TEST(HhtShellPhysics,
     C08aRigidOrbitMatchesIndependentVolumeForcePowerAndMomentum) {
  const auto fixture = MakeFixture(1U);
  const auto parameters = Require(fesa::MakeHhtParameters(-0.05));
  const auto initial =
      Evaluate(fixture, EmptyFrame(fixture), 1.0 / 16.0, parameters, false);
  const auto scaled_stiffness = Require(fesa::hht_internal::ScaledFree(
      initial.internal_tangent, fixture.dofs, fixture.length));
  double maximum_rotation = 0.0;
  for (std::size_t sample = 0; sample <= 16; ++sample) {
    const double time = static_cast<double>(sample) / 16.0;
    SCOPED_TRACE(time);
    auto exact = RigidOrbit(fixture, time);
    const auto oracle = Integrate(fixture, exact);
    exact.external_force = oracle.force;
    const auto assembly =
        Evaluate(fixture, exact, 1.0 / 16.0, parameters, true);
    auto difference = assembly.frame.inertia_force;
    difference.Axpy(-1.0, oracle.force);
    EXPECT_LE(difference.Norm() / oracle.force.Norm(), 1.0e-12);
    EXPECT_LE(std::abs(assembly.frame.kinetic_energy -
                       static_cast<double>(oracle.kinetic)) /
                  static_cast<double>(oracle.kinetic),
              1.0e-12);
    EXPECT_LE(std::abs(oracle.force.Dot(exact.velocity) -
                       static_cast<double>(oracle.power)) /
                  static_cast<double>(oracle.power_scale),
              1.0e-10);
    EXPECT_LE(assembly.frame.physical_strain_energy / 800.0, 1.0e-10);
    EXPECT_LE(assembly.frame.drilling_energy / 800.0, 1.0e-10);
    fesa::Vector scaled_q = exact.generalized_coordinates;
    fesa::Vector scaled_force = assembly.frame.internal_force;
    for (std::size_t i = 0; i < scaled_q.Size(); ++i) {
      if (i % 6U < 3U) {
        scaled_q[i] /= fixture.length;
        scaled_force[i] *= fixture.length;
      }
    }
    const double force_scale = SparseNorm(scaled_stiffness) * scaled_q.Norm();
    if (force_scale == 0.0)
      EXPECT_DOUBLE_EQ(scaled_force.Norm(), 0.0);
    else
      EXPECT_LE(scaled_force.Norm() / force_scale, 1.0e-10);
    double momentum_error = 0.0;
    double angular_error = 0.0;
    for (std::size_t c = 0; c < 3; ++c) {
      momentum_error = std::hypot(momentum_error,
                                  assembly.frame.linear_momentum[c] -
                                      static_cast<double>(oracle.momentum[c]));
      angular_error = std::hypot(
          angular_error, assembly.frame.angular_momentum[c] -
                             static_cast<double>(oracle.angular_momentum[c]));
    }
    EXPECT_LE(momentum_error / static_cast<double>(oracle.momentum_scale),
              1.0e-12);
    EXPECT_LE(
        angular_error / static_cast<double>(oracle.angular_momentum_scale),
        1.0e-12);
    auto balance = assembly.frame;
    ManufacturedHidden(assembly, RigidForceRates(fixture, time), &balance);
    for (std::size_t node = 0; node < 4U; ++node) {
      for (const auto quantity :
           {std::make_pair(&balance.hidden_constraint,
                           &balance.hidden_constraint_scale),
            std::make_pair(&balance.hidden_velocity_defect,
                           &balance.hidden_velocity_defect_scale),
            std::make_pair(&balance.hidden_acceleration_defect,
                           &balance.hidden_acceleration_defect_scale)}) {
        const double value = std::abs((*quantity.first)[node]);
        const double scale = (*quantity.second)[node];
        if (scale == 0.0)
          EXPECT_DOUBLE_EQ(value, 0.0);
        else
          EXPECT_LE(value / scale, 1.0e-10);
      }
    }
    EXPECT_TRUE(fesa::hht_internal::CompleteEquilibrium(
                    *fixture.domain, fixture.dofs, assembly, &balance)
                    .IsOk());
    const double mass_norm = SparseNorm(assembly.mass);
    for (std::size_t node = 0; node < 4U; ++node) {
      fesa::Vector null(24U);
      for (std::size_t c = 0; c < 3U; ++c)
        null[6U * node + 3U + c] = assembly.rotations[node].drill_null[c];
      EXPECT_LE(assembly.mass.Multiply(null).Norm() / mass_norm, 1.0e-12);
    }
    maximum_rotation = std::max(
        maximum_rotation,
        static_cast<double>(Norm(Triple(exact.generalized_coordinates, 3U))));
  }
  std::cout << "C08a pointwise independent-volume max_rotation="
            << std::setprecision(17) << maximum_rotation << '\n';
  EXPECT_GT(maximum_rotation, 0.5);
}

// At q0=0 the normalized null is e_phi3 and f_int=0, hence psi_q is the
// corresponding stiffness row. Determine all null velocities, rather than
// assigning the unobservable director-axis chart rate arbitrarily.
void InitializeNullVelocity(const Fixture& fixture, const Assembly& initial,
                            Frame* frame,
                            fesa::MklPardisoTangentSolver* solver) {
  const auto nodes = fixture.domain->Nodes().size();
  std::vector<fesa::CooContribution> entries;
  fesa::Vector rhs(nodes);
  for (std::size_t node = 0; node < nodes; ++node) {
    const auto row = 6U * node + 5U;
    for (auto offset = initial.internal_tangent.RowOffsets()[row];
         offset < initial.internal_tangent.RowOffsets()[row + 1U]; ++offset) {
      const auto column = initial.internal_tangent.ColumnIndices()[offset];
      const double value = initial.internal_tangent.Values()[offset];
      if (column % 6U == 5U)
        entries.push_back({node, column / 6U, value, 0U, entries.size()});
      else
        rhs[node] -= value * frame->velocity[column];
    }
  }
  const auto regularity =
      Require(fesa::hht_internal::MakeSparse(nodes, nodes, std::move(entries)));
  fesa::Vector null_velocity(nodes);
  Require(solver->FactorizeFresh(regularity));
  Require(solver->SolveCorrection(rhs, null_velocity));
  Require(fesa::hht_internal::CheckLinearSolve(regularity, rhs, null_velocity));
  for (std::size_t node = 0; node < nodes; ++node)
    frame->velocity[6U * node + 5U] = null_velocity[node];
}

struct GridEvidence {
  double h;
  double end_time;
  std::vector<Frame> history;
  std::array<double, 6> solve_noise{};
  std::array<double, 6> arithmetic_noise{};
  double energy_error = 0.0;
  double momentum_error = 0.0;
  double angular_momentum_error = 0.0;
  double maximum_rotation = 0.0;
  double algorithmic_acceleration_gap = 0.0;
  double hidden_velocity = 0.0;
  double hidden_acceleration = 0.0;
  struct RecoveredAcceleration {
    fesa::Vector value;
    double pre_hidden2 = 0.0;
    double pre_hidden2_relative = 0.0;
    double post_hidden2 = 0.0;
    double post_hidden2_relative = 0.0;
    double mass_delta_a_relative = 0.0;
  };
  // These are diagnostic-only: RunPlate never substitutes them into Frame.
  std::vector<RecoveredAcceleration> recovered_acceleration;
};

fesa::Vector ScaledResidual(const Fixture& fixture, const Frame& frame) {
  fesa::Vector rhs(fixture.dofs.FreeDofCount());
  for (std::size_t index = 0; index < rhs.Size(); ++index) {
    const auto full = fixture.dofs.FreeDofs()[index];
    rhs[index] =
        -frame.residual[full] * (full % 6U < 3U ? fixture.length : 1.0);
  }
  return rhs;
}

struct LastSolve {
  fesa::SparseMatrix tangent;
  fesa::Vector rhs;
  fesa::Vector correction;
};

void RecordNoise(const Fixture& fixture, const Assembly& accepted,
                 const fesa::HhtPredictor& predictor, const LastSolve& last,
                 const fesa::HhtParameters& p,
                 fesa::MklPardisoTangentSolver* solver,
                 GridEvidence* evidence) {
  auto residual = last.tangent.Multiply(last.correction);
  residual.Axpy(-1.0, last.rhs);
  residual.Scale(-1.0);
  fesa::Vector linear_defect(residual.Size());
  Require(solver->FactorizeFresh(last.tangent));
  Require(solver->SolveCorrection(residual, linear_defect));
  Require(fesa::hht_internal::CheckLinearSolve(last.tangent, residual,
                                               linear_defect));
  const auto endpoint = Require(fesa::hht_internal::ScaledFree(
      accepted.effective_tangent, fixture.dofs, fixture.length));
  const auto rhs = ScaledResidual(fixture, accepted.frame);
  fesa::Vector nonlinear_defect(rhs.Size());
  Require(solver->FactorizeFresh(endpoint));
  Require(solver->SolveCorrection(rhs, nonlinear_defect));
  Require(
      fesa::hht_internal::CheckLinearSolve(endpoint, rhs, nonlinear_defect));
  const double u = std::ldexp(1.0, -53);
  const double gamma32 = 32.0 * u / (1.0 - 32.0 * u);
  const std::array<double, 3> propagation{
      1.0, evidence->end_time * p.gamma / (p.beta * evidence->h),
      evidence->end_time * evidence->end_time /
          ((1.0 + p.alpha) * p.beta * evidence->h * evidence->h)};
  for (std::size_t node = 0; node < fixture.domain->Nodes().size(); ++node) {
    for (std::size_t channel = 0; channel < 2; ++channel) {
      const auto offset = 6U * node + 3U * channel;
      const double scale = channel == 0U ? fixture.length : 1.0;
      const double eta =
          static_cast<double>(Norm(Triple(last.correction, offset)) +
                              Norm(Triple(linear_defect, offset)) +
                              Norm(Triple(nonlinear_defect, offset)));
      const double arithmetic =
          gamma32 *
          static_cast<double>(
              Norm(Triple(accepted.frame.generalized_coordinates, offset)) +
              Norm(Triple(predictor.coordinates, offset))) /
          scale;
      for (std::size_t rate = 0; rate < 3; ++rate) {
        evidence->solve_noise[2U * rate + channel] =
            std::max(evidence->solve_noise[2U * rate + channel],
                     propagation[rate] * eta);
        evidence->arithmetic_noise[2U * rate + channel] =
            std::max(evidence->arithmetic_noise[2U * rate + channel],
                     propagation[rate] * arithmetic);
      }
    }
  }
}

GridEvidence::RecoveredAcceleration RecoverAcceleration(
    const Fixture& fixture, const Assembly& assembly,
    fesa::MklPardisoTangentSolver* solver, const bool manufactured) {
  const auto nodes = assembly.rotations.size();
  GridEvidence::RecoveredAcceleration result{assembly.frame.acceleration};
  std::vector<double> weights(nodes);
  std::vector<fesa::CooContribution> entries;
  fesa::Vector rhs(nodes);
  for (std::size_t row = 0; row < nodes; ++row) {
    const auto& rotation = assembly.rotations[row];
    fesa::Vector null(assembly.frame.acceleration.Size());
    for (std::size_t c = 0; c < 3U; ++c)
      null[6U * row + 3U + c] = rotation.drill_null[c];
    const auto metric_null = assembly.drilling_metric.Multiply(null);
    const double stiffness = null.Dot(metric_null);
    if (!std::isfinite(stiffness) || !(stiffness > 0.0))
      throw std::runtime_error{"C08 recovery requires positive null metric."};
    weights[row] = std::sqrt(stiffness);
  }
  for (std::size_t row = 0; row < nodes; ++row) {
    const auto& rotation = assembly.rotations[row];
    const double raw = assembly.frame.hidden_acceleration_defect[row];
    const double raw_scale =
        assembly.frame.hidden_acceleration_defect_scale[row];
    result.pre_hidden2 = std::max(result.pre_hidden2, std::abs(raw));
    result.pre_hidden2_relative =
        std::max(result.pre_hidden2_relative,
                 raw_scale == 0.0 ? std::abs(raw) : std::abs(raw) / raw_scale);
    rhs[row] = -raw / weights[row];
    for (std::size_t c = 0; c < 3U; ++c) {
      const auto full = 6U * row + 3U + c;
      const double force = assembly.frame.internal_force[full] -
                           assembly.frame.external_force[full];
      for (std::size_t d = 0; d < 3U; ++d)
        entries.push_back({row, row,
                           rotation.null_gradient(c, d) * force *
                               rotation.drill_null[d] / weights[row] /
                               weights[row],
                           0U, entries.size()});
      for (std::size_t index = assembly.internal_tangent.RowOffsets()[full];
           index < assembly.internal_tangent.RowOffsets()[full + 1U]; ++index) {
        const auto column = assembly.internal_tangent.ColumnIndices()[index];
        if (column % 6U < 3U) continue;
        const auto column_node = column / 6U;
        entries.push_back(
            {row, column_node,
             rotation.drill_null[c] *
                 assembly.internal_tangent.Values()[index] *
                 assembly.rotations[column_node].drill_null[column % 6U - 3U] /
                 weights[row] / weights[column_node],
             0U, entries.size()});
      }
    }
  }
  const auto regularity =
      Require(fesa::hht_internal::MakeSparse(nodes, nodes, std::move(entries)));
  fesa::Vector scaled_eta(nodes);
  Require(solver->FactorizeFresh(regularity));
  Require(solver->SolveCorrection(rhs, scaled_eta));
  Require(fesa::hht_internal::CheckLinearSolve(regularity, rhs, scaled_eta));
  fesa::Vector delta(result.value.Size());
  for (std::size_t node = 0; node < nodes; ++node)
    for (std::size_t c = 0; c < 3U; ++c)
      delta[6U * node + 3U + c] = assembly.rotations[node].drill_null[c] *
                                  scaled_eta[node] / weights[node];
  result.value.Axpy(1.0, delta);
  const auto mass_delta = assembly.mass.Multiply(delta);
  const double mass_scale = SparseNorm(assembly.mass) * delta.Norm();
  result.mass_delta_a_relative =
      mass_scale == 0.0 ? mass_delta.Norm() : mass_delta.Norm() / mass_scale;
  if (manufactured) {
    auto post = assembly.frame;
    post.acceleration = result.value;
    ManufacturedHidden(assembly, RigidForceRates(fixture, post.time), &post);
    for (std::size_t row = 0; row < nodes; ++row) {
      const double value = std::abs(post.hidden_acceleration_defect[row]);
      const double scale = post.hidden_acceleration_defect_scale[row];
      result.post_hidden2 = std::max(result.post_hidden2, value);
      result.post_hidden2_relative = std::max(
          result.post_hidden2_relative, scale == 0.0 ? value : value / scale);
    }
    return result;
  }
  for (std::size_t row = 0; row < nodes; ++row) {
    fesa::Vector gradient(result.value.Size());
    const auto& rotation = assembly.rotations[row];
    for (std::size_t c = 0; c < 3U; ++c) {
      const auto full = 6U * row + 3U + c;
      const double force = assembly.frame.internal_force[full] -
                           assembly.frame.external_force[full];
      for (std::size_t d = 0; d < 3U; ++d) {
        const auto column = 6U * row + 3U + d;
        gradient[column] += rotation.null_gradient(c, d) * force;
      }
      for (std::size_t index = assembly.internal_tangent.RowOffsets()[full];
           index < assembly.internal_tangent.RowOffsets()[full + 1U]; ++index) {
        const auto column = assembly.internal_tangent.ColumnIndices()[index];
        gradient[column] +=
            rotation.drill_null[c] * assembly.internal_tangent.Values()[index];
      }
    }
    const double post =
        assembly.hidden_second[row] + gradient.Dot(result.value);
    double post_scale = assembly.hidden_second_scale[row];
    for (std::size_t column = 0; column < gradient.Size(); ++column)
      post_scale += std::abs(gradient[column] * result.value[column]);
    result.post_hidden2 = std::max(result.post_hidden2, std::abs(post));
    result.post_hidden2_relative = std::max(
        result.post_hidden2_relative,
        post_scale == 0.0 ? std::abs(post) : std::abs(post) / post_scale);
  }
  return result;
}

GridEvidence RunPlate(const Fixture& fixture, const double h,
                      const bool manufactured) {
  GridEvidence evidence{h, manufactured ? 1.0 : kEndTime};
  const auto p = Require(fesa::MakeHhtParameters(-0.05));
  fesa::MklPardisoSolver positive;
  fesa::MklPardisoTangentSolver tangent;
  Frame frame = manufactured ? RigidOrbit(fixture, 0.0) : EmptyFrame(fixture);
  const L3 translation{0.125L, -0.25L, 0.0625L};
  const L3 omega{0.8L, 0.6L, 0.4L};
  if (manufactured) {
    frame.external_force = Integrate(fixture, frame).force;
  } else {
    for (std::size_t node = 0; node < fixture.domain->Nodes().size(); ++node) {
      const auto x = Coordinates(fixture.domain->Nodes()[node]);
      L3 velocity = Add(translation, Cross(omega, x));
      velocity[2] += 0.02L * x[0] * x[1];
      Store(velocity, &frame.velocity, 6U * node);
      Store({0.8L, 0.6L, 0.0L}, &frame.velocity, 6U * node + 3U);
    }
  }
  auto assembly = Evaluate(fixture, frame, h, p, true);
  if (!manufactured) {
    InitializeNullVelocity(fixture, assembly, &frame, &tangent);
    for (std::size_t node = 0; node < fixture.domain->Nodes().size(); ++node)
      EXPECT_NEAR(frame.velocity[6U * node + 5U], 0.4, 1.0e-12);
    assembly = Evaluate(fixture, frame, h, p, true);
  }
  const auto null_nodes =
      fesa::hht_internal::NullNodes(fixture.dofs, assembly.rotations);
  if (!manufactured)
    frame.acceleration = Require(fesa::hht_internal::InitializeAcceleration(
        *fixture.domain, fixture.dofs, assembly, fixture.length, positive,
        tangent));
  frame.algorithmic_acceleration = frame.acceleration;
  assembly = Evaluate(fixture, frame, h, p, true);
  frame = assembly.frame;
  if (manufactured)
    CompleteManufacturedConstraints(fixture, assembly, null_nodes, &positive,
                                    &tangent, &frame);
  else
    Require(fesa::hht_internal::CompleteConstraintEvidence(
        *fixture.domain, fixture.dofs, assembly, fixture.length, null_nodes,
        positive, tangent, &frame));
  Require(fesa::hht_internal::CompleteEquilibrium(*fixture.domain, fixture.dofs,
                                                  assembly, &frame));
  const auto initial_oracle = Integrate(fixture, frame);
  const double energy0 = frame.kinetic_energy + frame.physical_strain_energy +
                         frame.drilling_energy;
  EXPECT_GT(energy0, 0.0);
  EXPECT_LE(std::abs(frame.kinetic_energy -
                     static_cast<double>(initial_oracle.kinetic)) /
                energy0,
            1.0e-12);
  const auto initial_momentum = frame.linear_momentum;
  const auto initial_angular_momentum = frame.angular_momentum;
  evidence.history.push_back(frame);
  evidence.recovered_acceleration.push_back({frame.acceleration});
  const auto increments = static_cast<std::uint32_t>(evidence.end_time / h);
  for (std::uint32_t increment = 1U; increment <= increments; ++increment) {
    const Frame old = frame;
    const auto predictor =
        Require(fesa::PredictHht(old.generalized_coordinates, old.velocity,
                                 old.algorithmic_acceleration, h, p));
    frame.generalized_coordinates = old.generalized_coordinates;
    frame.is_initial = false;
    frame.frame_index = increment;
    frame.increment_index = increment;
    frame.time = static_cast<double>(increment) * h;
    if (manufactured)
      frame.external_force =
          Integrate(fixture, RigidOrbit(fixture, frame.time)).force;
    bool correction_pass = true;
    std::optional<LastSolve> last;
    bool converged = false;
    for (std::uint32_t iteration = 0U; iteration <= 100U; ++iteration) {
      auto rates = Require(fesa::EvaluateHhtRates(
          frame.generalized_coordinates, predictor, old.acceleration, h, p));
      frame.velocity = std::move(rates.velocity);
      frame.acceleration = std::move(rates.acceleration);
      frame.algorithmic_acceleration =
          std::move(rates.algorithmic_acceleration);
      assembly = Evaluate(fixture, frame, h, p, false);
      frame = assembly.frame;
      if (correction_pass && frame.free_residual_metrics[0] <= 1.0e-10 &&
          frame.free_residual_metrics[1] <= 1.0e-10) {
        assembly = Evaluate(fixture, frame, h, p, true);
        frame = assembly.frame;
        if (manufactured)
          CompleteManufacturedConstraints(fixture, assembly, null_nodes,
                                          &positive, &tangent, &frame);
        else
          Require(fesa::hht_internal::CompleteConstraintEvidence(
              *fixture.domain, fixture.dofs, assembly, fixture.length,
              null_nodes, positive, tangent, &frame));
        Require(fesa::hht_internal::CompleteEquilibrium(
            *fixture.domain, fixture.dofs, assembly, &frame));
        if (!last)
          throw std::runtime_error{
              "C08 nontrivial step lacks N7 last-solve evidence."};
        RecordNoise(fixture, assembly, predictor, *last, p, &tangent,
                    &evidence);
        auto recovery_assembly = assembly;
        recovery_assembly.frame = frame;
        evidence.recovered_acceleration.push_back(RecoverAcceleration(
            fixture, recovery_assembly, &tangent, manufactured));
        frame.iteration_count = iteration;
        converged = true;
        break;
      }
      if (iteration == 100U) break;
      auto matrix = Require(fesa::hht_internal::ScaledFree(
          assembly.effective_tangent, fixture.dofs, fixture.length));
      auto rhs = ScaledResidual(fixture, frame);
      fesa::Vector correction(rhs.Size());
      Require(tangent.FactorizeFresh(matrix));
      Require(tangent.SolveCorrection(rhs, correction));
      Require(fesa::hht_internal::CheckLinearSolve(matrix, rhs, correction));
      std::array<double, 2> correction_norm{};
      for (std::size_t i = 0; i < correction.Size(); ++i) {
        const auto full = fixture.dofs.FreeDofs()[i];
        frame.generalized_coordinates[full] +=
            correction[i] * (full % 6U < 3U ? fixture.length : 1.0);
        const auto channel = full % 6U < 3U ? 0U : 1U;
        correction_norm[channel] =
            std::max(correction_norm[channel], std::abs(correction[i]));
      }
      correction_pass = true;
      for (const double norm : correction_norm)
        for (const double factor :
             {1.0, p.gamma / p.beta, 1.0 / ((1.0 + p.alpha) * p.beta),
              1.0 / p.beta})
          correction_pass = correction_pass && factor * norm <= 1.0e-10;
      last.emplace(
          LastSolve{std::move(matrix), std::move(rhs), std::move(correction)});
    }
    if (!converged)
      throw std::runtime_error{
          "C08 exhausted 100 Newton corrections at h=" + std::to_string(h) +
          ", increment=" + std::to_string(increment)};
    const double energy = frame.kinetic_energy + frame.physical_strain_energy +
                          frame.drilling_energy;
    if (!manufactured)
      evidence.energy_error =
          std::max(evidence.energy_error, std::abs(energy - energy0) / energy0);
    double dp = 0.0;
    double dl = 0.0;
    for (std::size_t c = 0; c < 3; ++c) {
      dp = std::hypot(dp, frame.linear_momentum[c] - initial_momentum[c]);
      dl = std::hypot(dl,
                      frame.angular_momentum[c] - initial_angular_momentum[c]);
    }
    if (!manufactured) {
      evidence.momentum_error =
          std::max(evidence.momentum_error,
                   dp / static_cast<double>(initial_oracle.momentum_scale));
      evidence.angular_momentum_error = std::max(
          evidence.angular_momentum_error,
          dl / static_cast<double>(initial_oracle.angular_momentum_scale));
    }
    for (std::size_t node = 0; node < fixture.domain->Nodes().size(); ++node)
      evidence.maximum_rotation = std::max(
          evidence.maximum_rotation,
          static_cast<double>(
              Norm(Triple(frame.generalized_coordinates, 6U * node + 3U))));
    for (std::size_t node = 0; node < fixture.domain->Nodes().size(); ++node)
      for (std::size_t channel = 0; channel < 2U; ++channel) {
        const auto offset = 6U * node + 3U * channel;
        const double scale = channel == 0U ? fixture.length : 1.0;
        evidence.algorithmic_acceleration_gap = std::max(
            evidence.algorithmic_acceleration_gap,
            evidence.end_time * evidence.end_time / scale *
                static_cast<double>(Norm(
                    Add(Triple(frame.algorithmic_acceleration, offset),
                        Scale(Triple(frame.acceleration, offset), -1.0L)))));
      }
    for (std::size_t i = 0; i < frame.hidden_constraint.Size(); ++i) {
      if (frame.hidden_constraint_scale[i] == 0.0)
        EXPECT_DOUBLE_EQ(frame.hidden_constraint[i], 0.0);
      else
        EXPECT_LE(std::abs(frame.hidden_constraint[i]) /
                      frame.hidden_constraint_scale[i],
                  1.0e-10);
      const double vscale = frame.hidden_velocity_defect_scale[i];
      const double ascale = frame.hidden_acceleration_defect_scale[i];
      if (vscale > 0.0)
        evidence.hidden_velocity =
            std::max(evidence.hidden_velocity,
                     std::abs(frame.hidden_velocity_defect[i]) / vscale);
      else
        EXPECT_DOUBLE_EQ(frame.hidden_velocity_defect[i], 0.0);
      if (ascale > 0.0)
        evidence.hidden_acceleration =
            std::max(evidence.hidden_acceleration,
                     std::abs(frame.hidden_acceleration_defect[i]) / ascale);
      else
        EXPECT_DOUBLE_EQ(frame.hidden_acceleration_defect[i], 0.0);
    }
    evidence.history.push_back(frame);
  }
  std::cout << std::setprecision(17) << (manufactured ? "C08a" : "C08b")
            << " h=" << h;
  if (!manufactured)
    std::cout << " energy=" << evidence.energy_error
              << " P=" << evidence.momentum_error
              << " L=" << evidence.angular_momentum_error;
  std::cout << " phi_max=" << evidence.maximum_rotation
            << " a_A_gap=" << evidence.algorithmic_acceleration_gap
            << " hidden_v=" << evidence.hidden_velocity
            << " hidden_a=" << evidence.hidden_acceleration << " solve_noise=";
  for (const double value : evidence.solve_noise) std::cout << value << ',';
  std::cout << " round_noise=";
  for (const double value : evidence.arithmetic_noise)
    std::cout << value << ',';
  std::cout << '\n';
  return evidence;
}

std::array<double, 6> Differences(const GridEvidence& coarse,
                                  const GridEvidence& fine,
                                  const double length) {
  std::array<double, 6> error{};
  for (std::size_t frame = 0; frame < coarse.history.size(); ++frame) {
    const auto& a = coarse.history[frame];
    const auto& b = fine.history[2U * frame];
    const std::array<const fesa::Vector*, 3> first{
        &a.generalized_coordinates, &a.velocity, &a.acceleration};
    const std::array<const fesa::Vector*, 3> second{
        &b.generalized_coordinates, &b.velocity, &b.acceleration};
    for (std::size_t node = 0; node < a.nodal_internal_order.size(); ++node) {
      for (std::size_t channel = 0; channel < 2; ++channel) {
        const auto offset = 6U * node + 3U * channel;
        for (std::size_t rate = 0; rate < 3; ++rate) {
          const double normalization =
              (rate == 0   ? 1.0
               : rate == 1 ? coarse.end_time
                           : coarse.end_time * coarse.end_time) /
              (channel == 0U ? length : 1.0);
          const double difference = static_cast<double>(
              Norm(Add(Triple(*first[rate], offset),
                       Scale(Triple(*second[rate], offset), -1.0L))));
          error[2U * rate + channel] =
              std::max(error[2U * rate + channel], normalization * difference);
        }
      }
    }
  }
  return error;
}

std::array<double, 6> FrameError(const Frame& actual, const Frame& exact,
                                 const double length, const double end_time) {
  std::array<double, 6> error{};
  const std::array<const fesa::Vector*, 3> first{
      &actual.generalized_coordinates, &actual.velocity, &actual.acceleration};
  const std::array<const fesa::Vector*, 3> second{
      &exact.generalized_coordinates, &exact.velocity, &exact.acceleration};
  for (std::size_t node = 0; node < actual.generalized_coordinates.Size() / 6U;
       ++node) {
    for (std::size_t channel = 0; channel < 2U; ++channel) {
      const auto offset = 6U * node + 3U * channel;
      for (std::size_t rate = 0; rate < 3U; ++rate) {
        const double factor = (rate == 0U   ? 1.0
                               : rate == 1U ? end_time
                                            : end_time * end_time) /
                              (channel == 0U ? length : 1.0);
        error[2U * rate + channel] = std::max(
            error[2U * rate + channel],
            factor * static_cast<double>(Norm(
                         Add(Triple(*first[rate], offset),
                             Scale(Triple(*second[rate], offset), -1.0L)))));
      }
    }
  }
  return error;
}

std::array<double, 6> RecoveredFrameError(
    const Frame& actual, const fesa::Vector& recovered_acceleration,
    const Frame& exact, const double length, const double end_time) {
  auto recovered = actual;
  recovered.acceleration = recovered_acceleration;
  return FrameError(recovered, exact, length, end_time);
}

std::array<double, 6> RecoveredDifferences(const GridEvidence& coarse,
                                           const GridEvidence& fine,
                                           const double length) {
  std::array<double, 6> error{};
  for (std::size_t frame = 0; frame < coarse.history.size(); ++frame) {
    const auto& a = coarse.history[frame];
    const auto& b = fine.history[2U * frame];
    const std::array<const fesa::Vector*, 3> first{
        &a.generalized_coordinates, &a.velocity,
        &coarse.recovered_acceleration[frame].value};
    const std::array<const fesa::Vector*, 3> second{
        &b.generalized_coordinates, &b.velocity,
        &fine.recovered_acceleration[2U * frame].value};
    for (std::size_t node = 0; node < a.nodal_internal_order.size(); ++node) {
      for (std::size_t channel = 0; channel < 2; ++channel) {
        const auto offset = 6U * node + 3U * channel;
        const double scale = channel == 0U ? length : 1.0;
        for (std::size_t rate = 0; rate < 3; ++rate) {
          const double normalization =
              (rate == 0U   ? 1.0
               : rate == 1U ? coarse.end_time
                            : coarse.end_time * coarse.end_time) /
              scale;
          error[2U * rate + channel] = std::max(
              error[2U * rate + channel],
              normalization *
                  static_cast<double>(
                      Norm(Add(Triple(*first[rate], offset),
                               Scale(Triple(*second[rate], offset), -1.0L)))));
        }
      }
    }
  }
  return error;
}

std::array<double, 2> AccelerationSplit(const Frame& basis,
                                        const fesa::Vector& first,
                                        const fesa::Vector& second,
                                        const double end_time) {
  const auto dofs = basis.generalized_coordinates.Size();
  if (dofs % 6U != 0U || basis.velocity.Size() != dofs ||
      basis.acceleration.Size() != dofs || first.Size() != dofs ||
      second.Size() != dofs) {
    throw std::runtime_error{
        "C08 acceleration split requires matching nodal vectors."};
  }
  std::array<double, 2> split{};
  for (std::size_t node = 0; node < dofs / 6U; ++node) {
    const auto offset = 6U * node + 3U;
    auto vector3 = [](const L3& value) {
      return fesa::Vector3{static_cast<double>(value[0]),
                           static_cast<double>(value[1]),
                           static_cast<double>(value[2])};
    };
    const auto rotation = Require(fesa::EvaluateDynamicRotation(
        vector3(Triple(basis.generalized_coordinates, offset)), {0.0, 0.0, 1.0},
        vector3(Triple(basis.velocity, offset)),
        vector3(Triple(basis.acceleration, offset))));
    const L3 null{rotation.drill_null[0], rotation.drill_null[1],
                  rotation.drill_null[2]};
    const L3 difference =
        Add(Triple(first, offset), Scale(Triple(second, offset), -1.0L));
    split[0] = std::max(
        split[0], end_time * end_time *
                      static_cast<double>(std::abs(Dot(null, difference))));
    split[1] = std::max(split[1],
                        end_time * end_time *
                            static_cast<double>(Norm(Cross(null, difference))));
  }
  return split;
}

TEST(HhtShellPhysics, C08aRigidManufacturedOrbitFullCoordinateTimeConvergence) {
  const auto fixture = MakeFixture(1U);
  std::vector<GridEvidence> grids;
  std::vector<std::array<double, 6>> errors;
  for (const unsigned denominator : {16U, 32U, 64U, 128U, 256U}) {
    SCOPED_TRACE(denominator);
    grids.push_back(
        RunPlate(fixture, 1.0 / static_cast<double>(denominator), true));
    const auto& grid = grids.back();
    EXPECT_GT(grid.maximum_rotation, 0.5);
    std::array<double, 6> maximum{};
    std::array<double, 6> first{};
    std::array<double, 6> late{};
    std::array<double, 6> final{};
    std::array<double, 6> maximum_time{};
    std::array<double, 6> recovered_maximum{};
    std::array<double, 6> recovered_maximum_time{};
    std::array<double, 2> recovered_acceleration_split{};
    double recovery_pre_hidden2 = 0.0;
    double recovery_post_hidden2 = 0.0;
    double recovery_mass_delta_a = 0.0;
    double angular_acceleration_null_error = 0.0;
    double angular_acceleration_physical_error = 0.0;
    for (std::size_t index = 0; index < grid.history.size(); ++index) {
      const auto& frame = grid.history[index];
      const auto exact = RigidOrbit(fixture, frame.time);
      const auto error =
          FrameError(frame, exact, fixture.length, grid.end_time);
      const auto recovered_error =
          RecoveredFrameError(frame, grid.recovered_acceleration[index].value,
                              exact, fixture.length, grid.end_time);
      const auto split =
          AccelerationSplit(exact, grid.recovered_acceleration[index].value,
                            exact.acceleration, grid.end_time);
      for (std::size_t channel = 0; channel < 6U; ++channel) {
        if (recovered_error[channel] > recovered_maximum[channel])
          recovered_maximum_time[channel] = frame.time;
        recovered_maximum[channel] =
            std::max(recovered_maximum[channel], recovered_error[channel]);
      }
      for (std::size_t part = 0; part < 2U; ++part)
        recovered_acceleration_split[part] =
            std::max(recovered_acceleration_split[part], split[part]);
      recovery_pre_hidden2 =
          std::max(recovery_pre_hidden2,
                   grid.recovered_acceleration[index].pre_hidden2_relative);
      recovery_post_hidden2 =
          std::max(recovery_post_hidden2,
                   grid.recovered_acceleration[index].post_hidden2_relative);
      recovery_mass_delta_a =
          std::max(recovery_mass_delta_a,
                   grid.recovered_acceleration[index].mass_delta_a_relative);
      if (error[5] > maximum[5]) {
        double largest = 0.0;
        for (std::size_t node = 0; node < fixture.domain->Nodes().size();
             ++node) {
          const auto offset = 6U * node + 3U;
          const auto difference =
              Add(Triple(frame.acceleration, offset),
                  Scale(Triple(exact.acceleration, offset), -1.0L));
          const double magnitude = static_cast<double>(Norm(difference));
          if (magnitude <= largest) continue;
          largest = magnitude;
          const auto phi = Triple(frame.generalized_coordinates, offset);
          const auto velocity = Triple(frame.velocity, offset);
          const auto acceleration = Triple(frame.acceleration, offset);
          auto vector3 = [](const L3& value) {
            return fesa::Vector3{static_cast<double>(value[0]),
                                 static_cast<double>(value[1]),
                                 static_cast<double>(value[2])};
          };
          const auto rotation = Require(fesa::EvaluateDynamicRotation(
              vector3(phi), {0.0, 0.0, 1.0}, vector3(velocity),
              vector3(acceleration)));
          const L3 null{rotation.drill_null[0], rotation.drill_null[1],
                        rotation.drill_null[2]};
          angular_acceleration_null_error =
              static_cast<double>(std::abs(Dot(null, difference)));
          angular_acceleration_physical_error =
              static_cast<double>(Norm(Cross(null, difference)));
        }
      }
      for (std::size_t channel = 0; channel < 6U; ++channel) {
        if (error[channel] > maximum[channel])
          maximum_time[channel] = frame.time;
        maximum[channel] = std::max(maximum[channel], error[channel]);
        if (index == 1U) first[channel] = error[channel];
        if (frame.time >= grid.end_time / 4.0)
          late[channel] = std::max(late[channel], error[channel]);
      }
      final = error;
    }
    errors.push_back(maximum);
    std::cout << "C08a temporal_error h=" << grid.h;
    for (const auto& row :
         {std::make_pair(" first=", first), std::make_pair(" late=", late),
          std::make_pair(" maximum=", maximum),
          std::make_pair(" final=", final),
          std::make_pair(" maximum_time=", maximum_time)}) {
      std::cout << row.first;
      for (const double value : row.second) std::cout << value << ',';
    }
    std::cout << " max_angular_a_null=" << angular_acceleration_null_error
              << " max_angular_a_physical="
              << angular_acceleration_physical_error << '\n';
    std::cout << "C08a recovered_diagnostic h=" << grid.h << " max_time_error=";
    for (const double value : recovered_maximum) std::cout << value << ',';
    std::cout << " max_time=";
    for (const double value : recovered_maximum_time) std::cout << value << ',';
    std::cout << " angular_a_null=" << recovered_acceleration_split[0]
              << " angular_a_physical=" << recovered_acceleration_split[1]
              << " pre_hidden2=" << recovery_pre_hidden2
              << " post_hidden2=" << recovery_post_hidden2
              << " Mdeltaa=" << recovery_mass_delta_a << '\n';
  }
  for (std::size_t channel = 0; channel < 6U; ++channel) {
    SCOPED_TRACE(channel);
    EXPECT_LE(errors.back()[channel], 1.0e-4);
    for (std::size_t grid = 2U; grid < 4U; ++grid) {
      if (errors[grid][channel] == 0.0 && errors[grid + 1U][channel] == 0.0)
        continue;
      ASSERT_GT(errors[grid][channel], 0.0);
      ASSERT_GT(errors[grid + 1U][channel], 0.0);
      EXPECT_GE(std::log2(errors[grid][channel] / errors[grid + 1U][channel]),
                1.8);
    }
    for (std::size_t grid = 2U; grid < 5U; ++grid)
      EXPECT_GT(errors[grid][channel],
                100.0 * std::max(grids[grid].solve_noise[channel],
                                 grids[grid].arithmetic_noise[channel]))
          << "N7 inconclusive: temporal error does not resolve "
             "solve/arithmetic noise.";
  }
  for (const auto member :
       {&GridEvidence::hidden_velocity, &GridEvidence::hidden_acceleration}) {
    EXPECT_LT(grids[3].*member, grids[2].*member);
    EXPECT_LT(grids[4].*member, grids[3].*member);
  }
}

TEST(HhtShellPhysics, C08bInitialPhysicalSpectrumProbe) {
  const auto fixture = MakeFixture(2U);
  const auto parameters = Require(fesa::MakeHhtParameters(-0.05));
  const auto assembly =
      Evaluate(fixture, EmptyFrame(fixture), 1.0 / 64.0, parameters, false);
  const auto print_matrix = [](const char* prefix,
                               const fesa::SparseMatrix& matrix) {
    ASSERT_EQ(matrix.Rows(), 54U);
    ASSERT_EQ(matrix.Columns(), 54U);
    ASSERT_EQ(matrix.RowOffsets().size(), 55U);
    for (std::size_t row = 0; row < matrix.Rows(); ++row) {
      std::array<double, 54U> values{};
      for (std::size_t index = matrix.RowOffsets()[row];
           index < matrix.RowOffsets()[row + 1U]; ++index) {
        const auto column = matrix.ColumnIndices()[index];
        ASSERT_LT(column, values.size());
        values[column] = matrix.Values()[index];
      }
      std::cout << prefix << " row=" << row
                << " values=" << std::setprecision(17);
      for (const double value : values) std::cout << value << ',';
      std::cout << '\n';
    }
  };
  // Diagnostic matrix dump for offline static condensation/eigensolution only.
  print_matrix("C08b_initial_M", assembly.mass);
  print_matrix("C08b_initial_K", assembly.internal_tangent);
}

TEST(HhtShellPhysics,
     C08bFreeFiniteRotationPlateSelfConvergenceAndPhysicalInvariants) {
  const auto fixture = MakeFixture(2U);
  std::vector<GridEvidence> grids;
  for (const unsigned denominator : {64U, 128U, 256U, 512U, 1024U}) {
    SCOPED_TRACE(denominator);
    grids.push_back(
        RunPlate(fixture, 1.0 / static_cast<double>(denominator), false));
    EXPECT_GE(grids.back().maximum_rotation, 0.2);
  }
  std::vector<std::array<double, 6>> differences;
  for (std::size_t pair = 0; pair + 1U < grids.size(); ++pair) {
    differences.push_back(
        Differences(grids[pair], grids[pair + 1U], fixture.length));
    std::cout << "C08b self_difference h=" << grids[pair].h << ' ';
    for (const double value : differences.back()) std::cout << value << ',';
    std::cout << '\n';
    const auto recovered =
        RecoveredDifferences(grids[pair], grids[pair + 1U], fixture.length);
    std::array<double, 2> acceleration_split{};
    double recovery_pre_hidden2 = 0.0;
    double recovery_post_hidden2 = 0.0;
    double recovery_mass_delta_a = 0.0;
    for (std::size_t frame = 0; frame < grids[pair].history.size(); ++frame) {
      const auto fine_frame = 2U * frame;
      const auto split = AccelerationSplit(
          grids[pair + 1U].history[fine_frame],
          grids[pair].recovered_acceleration[frame].value,
          grids[pair + 1U].recovered_acceleration[fine_frame].value,
          grids[pair].end_time);
      for (std::size_t part = 0; part < 2U; ++part)
        acceleration_split[part] =
            std::max(acceleration_split[part], split[part]);
    }
    for (const auto* grid : {&grids[pair], &grids[pair + 1U]}) {
      for (const auto& recovery : grid->recovered_acceleration) {
        recovery_pre_hidden2 =
            std::max(recovery_pre_hidden2, recovery.pre_hidden2_relative);
        recovery_post_hidden2 =
            std::max(recovery_post_hidden2, recovery.post_hidden2_relative);
        recovery_mass_delta_a =
            std::max(recovery_mass_delta_a, recovery.mass_delta_a_relative);
      }
    }
    std::cout << "C08b recovered_diagnostic h=" << grids[pair].h
              << " common_time_difference=";
    for (const double value : recovered) std::cout << value << ',';
    std::cout << " angular_a_null=" << acceleration_split[0]
              << " angular_a_physical=" << acceleration_split[1]
              << " pre_hidden2=" << recovery_pre_hidden2
              << " post_hidden2=" << recovery_post_hidden2
              << " Mdeltaa=" << recovery_mass_delta_a << '\n';
  }
  for (std::size_t channel = 0; channel < 6; ++channel) {
    SCOPED_TRACE(channel);
    EXPECT_LE(differences.back()[channel], 1.0e-4);
    for (std::size_t pair = 1U; pair < 3U; ++pair) {
      const double first = differences[pair][channel];
      const double second = differences[pair + 1U][channel];
      if (first == 0.0 && second == 0.0) continue;
      ASSERT_GT(first, 0.0);
      ASSERT_GT(second, 0.0);
      EXPECT_GE(std::log2(first / second), 1.8);
    }
    for (std::size_t pair = 1U; pair < differences.size(); ++pair) {
      const double noise = std::max(grids[pair].solve_noise[channel],
                                    grids[pair].arithmetic_noise[channel]) +
                           std::max(grids[pair + 1U].solve_noise[channel],
                                    grids[pair + 1U].arithmetic_noise[channel]);
      EXPECT_GT(differences[pair][channel], 100.0 * noise)
          << "N7 inconclusive: temporal difference does not resolve "
             "solve/arithmetic noise.";
    }
  }
  for (const auto member :
       {&GridEvidence::energy_error, &GridEvidence::momentum_error,
        &GridEvidence::angular_momentum_error}) {
    EXPECT_LT(grids[3].*member, grids[2].*member);
    EXPECT_LT(grids[4].*member, grids[3].*member);
    EXPECT_LE(grids[4].*member, 0.6 * (grids[3].*member));
    EXPECT_LE(grids[4].*member, 1.0e-3);
  }
}

}  // namespace