#include "fesa/elements/mitc4_shell.h"

#include <gtest/gtest.h>

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

#include "fesa/math/vector3.h"

namespace {

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

fesa::Node Node(std::int64_t label, Vector3 coordinates) {
  return {{"Shell-Instance", label, std::to_string(label)},
          coordinates,
          {"mitc4-shell.inp", static_cast<std::size_t>(label + 1)}};
}

std::array<const fesa::Node*, 4> NodePointers(
    const std::array<fesa::Node, 4>& nodes) {
  return {&nodes[0], &nodes[1], &nodes[2], &nodes[3]};
}

fesa::ShellSection Section(double thickness = 2.0) {
  return {"Section", thickness, 0U, {"mitc4-shell.inp", 20U}};
}

fesa::LinearElasticMaterial Material(double youngs_modulus = 120.0,
                                     double poisson_ratio = 0.25) {
  return {"Material", youngs_modulus, poisson_ratio, {"mitc4-shell.inp", 21U}};
}

std::array<Vector3, 4> Directors(Vector3 director = {0.0, 0.0, 1.0}) {
  return {director, director, director, director};
}

double Dot(const Vector3& left, const Vector3& right) {
  return left[0] * right[0] + left[1] * right[1] + left[2] * right[2];
}

Vector3 Cross(const Vector3& left, const Vector3& right) {
  return {left[1] * right[2] - left[2] * right[1],
          left[2] * right[0] - left[0] * right[2],
          left[0] * right[1] - left[1] * right[0]};
}

double Norm(const Vector3& value) { return std::sqrt(Dot(value, value)); }

void ExpectVectorNear(const Vector3& actual, const Vector3& expected,
                      double tolerance = 1.0e-12) {
  for (std::size_t component = 0U; component < actual.size(); ++component) {
    EXPECT_NEAR(actual[component], expected[component], tolerance);
  }
}

void ExpectOrthonormalRightHanded(const fesa::Mitc4LocalFrame& frame) {
  EXPECT_NEAR(Norm(frame.e1), 1.0, 1.0e-12);
  EXPECT_NEAR(Norm(frame.e2), 1.0, 1.0e-12);
  EXPECT_NEAR(Norm(frame.e3), 1.0, 1.0e-12);
  EXPECT_NEAR(Dot(frame.e1, frame.e2), 0.0, 1.0e-12);
  EXPECT_NEAR(Dot(frame.e1, frame.e3), 0.0, 1.0e-12);
  EXPECT_NEAR(Dot(frame.e2, frame.e3), 0.0, 1.0e-12);
  ExpectVectorNear(Cross(frame.e1, frame.e2), frame.e3);
}

void ExpectMatrixNear(const fesa::Matrix& actual, const fesa::Matrix& expected,
                      double tolerance = 1.0e-12) {
  ASSERT_EQ(actual.Rows(), expected.Rows());
  ASSERT_EQ(actual.Columns(), expected.Columns());
  for (std::size_t row = 0U; row < actual.Rows(); ++row) {
    for (std::size_t column = 0U; column < actual.Columns(); ++column) {
      EXPECT_NEAR(actual(row, column), expected(row, column), tolerance)
          << "at (" << row << ", " << column << ")";
    }
  }
}

void ExpectSymmetric(const fesa::Matrix& matrix) {
  ASSERT_EQ(matrix.Rows(), matrix.Columns());
  for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
    for (std::size_t column = 0U; column < matrix.Columns(); ++column) {
      EXPECT_NEAR(matrix(row, column), matrix(column, row), 1.0e-12);
    }
  }
}

bool HasPositiveCholeskyPivots(const fesa::Matrix& matrix) {
  if (matrix.Rows() != matrix.Columns()) {
    return false;
  }
  fesa::Matrix lower{matrix.Rows(), matrix.Columns()};
  for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
    for (std::size_t column = 0U; column <= row; ++column) {
      double value = matrix(row, column);
      for (std::size_t inner = 0U; inner < column; ++inner) {
        value -= lower(row, inner) * lower(column, inner);
      }
      if (row == column) {
        if (!std::isfinite(value) || !(value > 0.0)) {
          return false;
        }
        lower(row, column) = std::sqrt(value);
      } else {
        lower(row, column) = value / lower(column, column);
      }
    }
  }
  return true;
}

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

double ScaledSymmetryError(const fesa::Matrix& matrix,
                           std::size_t dofs_per_node, double element_length) {
  fesa::Matrix difference{matrix.Rows(), matrix.Columns()};
  fesa::Matrix scaled{matrix.Rows(), matrix.Columns()};
  for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
    const double row_scale = row % dofs_per_node < 3U ? element_length : 1.0;
    for (std::size_t column = 0U; column < matrix.Columns(); ++column) {
      const double column_scale =
          column % dofs_per_node < 3U ? element_length : 1.0;
      scaled(row, column) = row_scale * matrix(row, column) * column_scale;
      difference(row, column) = row_scale *
                                (matrix(row, column) - matrix(column, row)) *
                                column_scale;
    }
  }
  return FrobeniusNorm(difference) / FrobeniusNorm(scaled);
}

fesa::Matrix ScaledStiffness(const fesa::Matrix& matrix,
                             std::size_t dofs_per_node, double element_length) {
  fesa::Matrix scaled{matrix.Rows(), matrix.Columns()};
  for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
    const double row_scale = row % dofs_per_node < 3U ? element_length : 1.0;
    for (std::size_t column = 0U; column < matrix.Columns(); ++column) {
      const double column_scale =
          column % dofs_per_node < 3U ? element_length : 1.0;
      scaled(row, column) = row_scale * matrix(row, column) * column_scale;
    }
  }
  return scaled;
}

std::vector<double> SymmetricEigenvalues(fesa::Matrix matrix) {
  if (matrix.Rows() != matrix.Columns()) {
    throw std::invalid_argument{
        "Symmetric eigensolve requires a square matrix."};
  }
  const std::size_t size = matrix.Rows();
  double matrix_scale = 0.0;
  for (std::size_t row = 0U; row < size; ++row) {
    for (std::size_t column = 0U; column < size; ++column) {
      matrix_scale = (std::max)(matrix_scale, std::abs(matrix(row, column)));
    }
  }
  if (matrix_scale != 0.0) {
    const double convergence_tolerance = 1.0e-14 * matrix_scale;
    const std::size_t iteration_limit = 100U * size * size;
    for (std::size_t iteration = 0U; iteration < iteration_limit; ++iteration) {
      std::size_t pivot_row = 0U;
      std::size_t pivot_column = 0U;
      double largest_off_diagonal = 0.0;
      for (std::size_t row = 0U; row < size; ++row) {
        for (std::size_t column = row + 1U; column < size; ++column) {
          const double candidate = std::abs(matrix(row, column));
          if (candidate > largest_off_diagonal) {
            largest_off_diagonal = candidate;
            pivot_row = row;
            pivot_column = column;
          }
        }
      }
      if (largest_off_diagonal <= convergence_tolerance) {
        break;
      }

      const double pivot = matrix(pivot_row, pivot_column);
      const double tau =
          (matrix(pivot_column, pivot_column) - matrix(pivot_row, pivot_row)) /
          (2.0 * pivot);
      const double tangent = tau >= 0.0
                                 ? 1.0 / (tau + std::sqrt(1.0 + tau * tau))
                                 : -1.0 / (-tau + std::sqrt(1.0 + tau * tau));
      const double cosine = 1.0 / std::sqrt(1.0 + tangent * tangent);
      const double sine = tangent * cosine;
      const double row_diagonal = matrix(pivot_row, pivot_row);
      const double column_diagonal = matrix(pivot_column, pivot_column);
      matrix(pivot_row, pivot_row) = row_diagonal - tangent * pivot;
      matrix(pivot_column, pivot_column) = column_diagonal + tangent * pivot;
      matrix(pivot_row, pivot_column) = 0.0;
      matrix(pivot_column, pivot_row) = 0.0;
      for (std::size_t index = 0U; index < size; ++index) {
        if (index == pivot_row || index == pivot_column) {
          continue;
        }
        const double row_value = matrix(index, pivot_row);
        const double column_value = matrix(index, pivot_column);
        const double rotated_row = cosine * row_value - sine * column_value;
        const double rotated_column = sine * row_value + cosine * column_value;
        matrix(index, pivot_row) = rotated_row;
        matrix(pivot_row, index) = rotated_row;
        matrix(index, pivot_column) = rotated_column;
        matrix(pivot_column, index) = rotated_column;
      }
    }
  }

  std::vector<double> eigenvalues(size);
  for (std::size_t index = 0U; index < size; ++index) {
    eigenvalues[index] = matrix(index, index);
  }
  return eigenvalues;
}

std::size_t NumericalRank(const fesa::Matrix& scaled_matrix) {
  const auto eigenvalues = SymmetricEigenvalues(scaled_matrix);
  double spectral_scale = 0.0;
  for (double eigenvalue : eigenvalues) {
    spectral_scale = (std::max)(spectral_scale, std::abs(eigenvalue));
  }
  return static_cast<std::size_t>(
      std::count_if(eigenvalues.begin(), eigenvalues.end(),
                    [spectral_scale](double eigenvalue) {
                      return std::abs(eigenvalue) > 1.0e-9 * spectral_scale;
                    }));
}

double SymmetricOperatorNorm(const fesa::Matrix& matrix) {
  const auto eigenvalues = SymmetricEigenvalues(matrix);
  double result = 0.0;
  for (double eigenvalue : eigenvalues) {
    result = (std::max)(result, std::abs(eigenvalue));
  }
  return result;
}

double QuadraticEnergy(const fesa::Matrix& stiffness,
                       const fesa::Vector& vector) {
  return 0.5 * vector.Dot(stiffness.Multiply(vector));
}

fesa::Vector PhysicalField(const std::array<std::array<double, 5>, 4>& values) {
  fesa::Vector result{20U};
  for (std::size_t node_index = 0U; node_index < 4U; ++node_index) {
    for (std::size_t component = 0U; component < 5U; ++component) {
      result[5U * node_index + component] = values[node_index][component];
    }
  }
  return result;
}

void ExpectStrain(const fesa::Mitc4Shell& shell, const fesa::Vector& field,
                  double xi, double eta, double zeta,
                  const std::array<double, 5>& expected) {
  const auto actual = shell.StrainDisplacement20(xi, eta, zeta).Multiply(field);
  for (std::size_t component = 0U; component < expected.size(); ++component) {
    EXPECT_NEAR(actual[component], expected[component], 1.0e-12)
        << "component " << component;
  }
}

fesa::Vector PhysicalRigidMode(const std::array<fesa::Node, 4>& nodes,
                               const Vector3& translation,
                               const Vector3& rotation) {
  fesa::Vector mode{24U};
  for (std::size_t node_index = 0U; node_index < nodes.size(); ++node_index) {
    const Vector3 rotational_translation =
        Cross(rotation, nodes[node_index].coordinates);
    const std::size_t offset = 6U * node_index;
    for (std::size_t component = 0U; component < 3U; ++component) {
      mode[offset + component] =
          translation[component] + rotational_translation[component];
    }
    // Remove the director-parallel component: it is numerical drilling,
    // not part of the five-DOF physical rigid motion.
    mode[offset + 3U] = rotation[0];
    mode[offset + 4U] = rotation[1];
    mode[offset + 5U] = 0.0;
  }
  return mode;
}

std::array<fesa::Node, 4> PlanarNodes() {
  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})};
}

}  // namespace

// MITC4-KIN-001
TEST(Mitc4ShellKinematics, ShapeFunctionsSatisfyNodalAndDerivativeIdentities) {
  constexpr std::array<Vector3, 4> kNaturalNodes{
      Vector3{-1.0, -1.0, 0.0}, Vector3{1.0, -1.0, 0.0}, Vector3{1.0, 1.0, 0.0},
      Vector3{-1.0, 1.0, 0.0}};

  for (std::size_t point = 0U; point < kNaturalNodes.size(); ++point) {
    const auto shape = fesa::Mitc4Shell::ShapeFunctions(
        kNaturalNodes[point][0], kNaturalNodes[point][1]);
    for (std::size_t node_index = 0U; node_index < kNaturalNodes.size();
         ++node_index) {
      EXPECT_DOUBLE_EQ(shape.values[node_index],
                       point == node_index ? 1.0 : 0.0);
    }
  }

  const auto shape = fesa::Mitc4Shell::ShapeFunctions(0.25, -0.5);
  double value_sum = 0.0;
  double xi_derivative_sum = 0.0;
  double eta_derivative_sum = 0.0;
  for (std::size_t node_index = 0U; node_index < 4U; ++node_index) {
    value_sum += shape.values[node_index];
    xi_derivative_sum += shape.xi_derivatives[node_index];
    eta_derivative_sum += shape.eta_derivatives[node_index];
  }
  EXPECT_DOUBLE_EQ(value_sum, 1.0);
  EXPECT_DOUBLE_EQ(xi_derivative_sum, 0.0);
  EXPECT_DOUBLE_EQ(eta_derivative_sum, 0.0);
  EXPECT_EQ(shape.values,
            (std::array<double, 4>{0.28125, 0.46875, 0.15625, 0.09375}));
}

// MITC4-KIN-002
TEST(Mitc4ShellKinematics,
     BuildsRightHandedFramesAndSeparatePhysicalDrillingMaps) {
  const std::array<fesa::Node, 4> nodes{
      Node(1, {0.0, -1.0, -1.0}), Node(2, {0.0, 1.0, -1.0}),
      Node(3, {0.0, 1.0, 1.0}), Node(4, {0.0, -1.0, 1.0})};
  const auto candidate = fesa::Mitc4Shell::Create(
      NodePointers(nodes), Directors({1.0, 0.0, 0.0}), Section(), Material());
  ASSERT_TRUE(candidate.HasValue());
  const auto& shell = candidate.Value();

  const auto frame = shell.LocalFrame(0.0, 0.0);
  ExpectVectorNear(frame.e1, {0.0, 1.0, 0.0});
  ExpectVectorNear(frame.e2, {0.0, 0.0, 1.0});
  ExpectVectorNear(frame.e3, {1.0, 0.0, 0.0});
  ExpectOrthonormalRightHanded(frame);

  const auto physical = shell.PhysicalTransformation20();
  const auto drilling = shell.DrillingTransformation4();
  ASSERT_EQ(physical.Rows(), 20U);
  ASSERT_EQ(physical.Columns(), 24U);
  ASSERT_EQ(drilling.Rows(), 4U);
  ASSERT_EQ(drilling.Columns(), 24U);
  for (std::size_t node_index = 0U; node_index < 4U; ++node_index) {
    const std::size_t physical_offset = 5U * node_index;
    const std::size_t global_offset = 6U * node_index;
    for (std::size_t component = 0U; component < 3U; ++component) {
      EXPECT_DOUBLE_EQ(
          physical(physical_offset + component, global_offset + component),
          1.0);
    }
    EXPECT_DOUBLE_EQ(physical(physical_offset + 3U, global_offset + 4U), 1.0);
    EXPECT_DOUBLE_EQ(physical(physical_offset + 4U, global_offset + 5U), 1.0);
    EXPECT_DOUBLE_EQ(drilling(node_index, global_offset + 3U), 1.0);

    for (std::size_t global_dof = 0U; global_dof < 24U; ++global_dof) {
      if (global_dof != global_offset + 4U) {
        EXPECT_DOUBLE_EQ(physical(physical_offset + 3U, global_dof), 0.0);
      }
      if (global_dof != global_offset + 5U) {
        EXPECT_DOUBLE_EQ(physical(physical_offset + 4U, global_dof), 0.0);
      }
      if (global_dof != global_offset + 3U) {
        EXPECT_DOUBLE_EQ(drilling(node_index, global_dof), 0.0);
      }
    }
  }

  auto invalid_directors = Directors({1.0, 0.0, 0.0});
  invalid_directors[2] = {0.0, 0.0, 0.0};
  EXPECT_FALSE(fesa::Mitc4Shell::Create(NodePointers(nodes), invalid_directors,
                                        Section(), Material())
                   .HasValue());
}

// MITC4-KIN-003
TEST(Mitc4ShellKinematics, FormsDirectColumnsAndAllCovariantTyingSamples) {
  const auto nodes = PlanarNodes();
  const auto candidate = fesa::Mitc4Shell::Create(
      NodePointers(nodes), Directors(), Section(), Material());
  ASSERT_TRUE(candidate.HasValue());
  const auto& shell = candidate.Value();

  const auto direct = shell.DirectStrainDisplacement20(0.0, 0.0, 0.5);
  ASSERT_EQ(direct.Rows(), 5U);
  ASSERT_EQ(direct.Columns(), 20U);
  EXPECT_DOUBLE_EQ(direct(0U, 0U), -0.25);
  EXPECT_DOUBLE_EQ(direct(0U, 4U), -0.125);
  EXPECT_DOUBLE_EQ(direct(1U, 1U), -0.25);
  EXPECT_DOUBLE_EQ(direct(1U, 3U), 0.125);
  EXPECT_DOUBLE_EQ(direct(2U, 0U), -0.25);
  EXPECT_DOUBLE_EQ(direct(2U, 1U), -0.25);
  EXPECT_DOUBLE_EQ(direct(2U, 3U), 0.125);
  EXPECT_DOUBLE_EQ(direct(2U, 4U), -0.125);
  EXPECT_DOUBLE_EQ(direct(3U, 2U), -0.25);
  EXPECT_DOUBLE_EQ(direct(3U, 4U), 0.25);
  EXPECT_DOUBLE_EQ(direct(4U, 2U), -0.25);
  EXPECT_DOUBLE_EQ(direct(4U, 3U), -0.25);

  const auto samples = shell.CovariantTyingShearSamples20();
  ASSERT_EQ(samples.Rows(), 4U);
  ASSERT_EQ(samples.Columns(), 20U);
  EXPECT_DOUBLE_EQ(samples(0U, 2U), -0.25);
  EXPECT_DOUBLE_EQ(samples(0U, 4U), 0.25);
  EXPECT_DOUBLE_EQ(samples(0U, 7U), 0.25);
  EXPECT_DOUBLE_EQ(samples(0U, 9U), 0.25);
  EXPECT_DOUBLE_EQ(samples(1U, 12U), 0.25);
  EXPECT_DOUBLE_EQ(samples(1U, 14U), 0.25);
  EXPECT_DOUBLE_EQ(samples(1U, 17U), -0.25);
  EXPECT_DOUBLE_EQ(samples(1U, 19U), 0.25);
  EXPECT_DOUBLE_EQ(samples(2U, 2U), -0.25);
  EXPECT_DOUBLE_EQ(samples(2U, 3U), -0.25);
  EXPECT_DOUBLE_EQ(samples(2U, 17U), 0.25);
  EXPECT_DOUBLE_EQ(samples(2U, 18U), -0.25);
  EXPECT_DOUBLE_EQ(samples(3U, 7U), -0.25);
  EXPECT_DOUBLE_EQ(samples(3U, 8U), -0.25);
  EXPECT_DOUBLE_EQ(samples(3U, 12U), 0.25);
  EXPECT_DOUBLE_EQ(samples(3U, 13U), -0.25);

  const auto weights = fesa::Mitc4Shell::TyingWeights(0.25, -0.5);
  EXPECT_EQ(weights.xi_zeta, (std::array<double, 2>{0.75, 0.25}));
  EXPECT_EQ(weights.eta_zeta, (std::array<double, 2>{0.375, 0.625}));

  const auto tied = shell.StrainDisplacement20(0.0, 0.0, 0.0);
  EXPECT_DOUBLE_EQ(tied(3U, 4U),
                   2.0 * (0.5 * samples(0U, 4U) + 0.5 * samples(1U, 4U)));
  EXPECT_DOUBLE_EQ(tied(4U, 3U),
                   2.0 * (0.5 * samples(2U, 3U) + 0.5 * samples(3U, 3U)));
}

// MITC4-KIN-004
TEST(Mitc4ShellConstitutive,
     BuildsExactPositiveDefiniteSectionMatricesAndRescalesUnits) {
  const auto nodes = PlanarNodes();
  const auto candidate = fesa::Mitc4Shell::Create(
      NodePointers(nodes), Directors(), Section(), Material());
  ASSERT_TRUE(candidate.HasValue());
  const auto& shell = candidate.Value();

  const auto cps = shell.PlaneStressConstitutive();
  const auto c5 = shell.MaterialConstitutive5();
  const auto a = shell.MembraneSectionMatrix();
  const auto d = shell.BendingSectionMatrix();
  const auto as = shell.TransverseShearSectionMatrix();
  EXPECT_EQ(cps.Rows(), 3U);
  EXPECT_EQ(cps.Columns(), 3U);
  EXPECT_EQ(c5.Rows(), 5U);
  EXPECT_EQ(c5.Columns(), 5U);
  EXPECT_EQ(a.Rows(), 3U);
  EXPECT_EQ(d.Rows(), 3U);
  EXPECT_EQ(as.Rows(), 2U);
  EXPECT_DOUBLE_EQ(cps(0U, 0U), 128.0);
  EXPECT_DOUBLE_EQ(cps(0U, 1U), 32.0);
  EXPECT_DOUBLE_EQ(cps(2U, 2U), 48.0);
  EXPECT_DOUBLE_EQ(c5(3U, 3U), 40.0);
  EXPECT_DOUBLE_EQ(c5(4U, 4U), 40.0);
  EXPECT_DOUBLE_EQ(a(0U, 0U), 256.0);
  EXPECT_NEAR(d(0U, 0U), 256.0 / 3.0, 1.0e-12);
  EXPECT_DOUBLE_EQ(as(0U, 0U), 80.0);
  ExpectSymmetric(cps);
  ExpectSymmetric(c5);
  EXPECT_TRUE(HasPositiveCholeskyPivots(cps));
  EXPECT_TRUE(HasPositiveCholeskyPivots(c5));
  EXPECT_TRUE(HasPositiveCholeskyPivots(a));
  EXPECT_TRUE(HasPositiveCholeskyPivots(d));
  EXPECT_TRUE(HasPositiveCholeskyPivots(as));

  constexpr double kForceScale = 7.0;
  constexpr double kLengthScale = 3.0;
  const auto scaled_candidate = fesa::Mitc4Shell::Create(
      NodePointers(nodes), Directors(), Section(2.0 * kLengthScale),
      Material(120.0 * kForceScale / (kLengthScale * kLengthScale), 0.25));
  ASSERT_TRUE(scaled_candidate.HasValue());
  const auto& scaled = scaled_candidate.Value();
  fesa::Matrix expected_cps{3U, 3U};
  fesa::Matrix expected_c5{5U, 5U};
  fesa::Matrix expected_a{3U, 3U};
  fesa::Matrix expected_d{3U, 3U};
  fesa::Matrix expected_as{2U, 2U};
  for (std::size_t row = 0U; row < 3U; ++row) {
    for (std::size_t column = 0U; column < 3U; ++column) {
      expected_cps(row, column) =
          cps(row, column) * kForceScale / (kLengthScale * kLengthScale);
      expected_a(row, column) = a(row, column) * kForceScale / kLengthScale;
      expected_d(row, column) = d(row, column) * kForceScale * kLengthScale;
    }
  }
  for (std::size_t row = 0U; row < 2U; ++row) {
    for (std::size_t column = 0U; column < 2U; ++column) {
      expected_as(row, column) = as(row, column) * kForceScale / kLengthScale;
    }
  }
  for (std::size_t row = 0U; row < 5U; ++row) {
    for (std::size_t column = 0U; column < 5U; ++column) {
      expected_c5(row, column) =
          c5(row, column) * kForceScale / (kLengthScale * kLengthScale);
    }
  }
  ExpectMatrixNear(scaled.PlaneStressConstitutive(), expected_cps);
  ExpectMatrixNear(scaled.MaterialConstitutive5(), expected_c5);
  ExpectMatrixNear(scaled.MembraneSectionMatrix(), expected_a);
  ExpectMatrixNear(scaled.BendingSectionMatrix(), expected_d);
  ExpectMatrixNear(scaled.TransverseShearSectionMatrix(), expected_as);

  EXPECT_FALSE(fesa::Mitc4Shell::Create(NodePointers(nodes), Directors(),
                                        Section(0.0), Material())
                   .HasValue());
  EXPECT_FALSE(fesa::Mitc4Shell::Create(NodePointers(nodes), Directors(),
                                        Section(), Material(0.0, 0.25))
                   .HasValue());
  EXPECT_FALSE(fesa::Mitc4Shell::Create(NodePointers(nodes), Directors(),
                                        Section(), Material(120.0, 0.5))
                   .HasValue());
}

// MITC4-KIN-005
TEST(Mitc4ShellKinematics, UsesOneFixedTwoByTwoByTwoQuadratureOrder) {
  const auto& points = fesa::Mitc4Shell::VolumeQuadrature();
  ASSERT_EQ(points.size(), 8U);
  const double g = 1.0 / std::sqrt(3.0);
  const std::array<Vector3, 8> expected{Vector3{-g, -g, -g}, Vector3{-g, -g, g},
                                        Vector3{g, -g, -g},  Vector3{g, -g, g},
                                        Vector3{g, g, -g},   Vector3{g, g, g},
                                        Vector3{-g, g, -g},  Vector3{-g, g, g}};
  for (std::size_t point = 0U; point < points.size(); ++point) {
    EXPECT_EQ(points[point].natural_coordinates, expected[point]);
    EXPECT_DOUBLE_EQ(points[point].weight, 1.0);
  }
}

// MITC4-KERNEL-001
TEST(Mitc4ShellKernel,
     FormsFiniteScaledSymmetricPhysicalAndStabilizedStiffness) {
  const auto nodes = PlanarNodes();
  const auto shell_candidate = fesa::Mitc4Shell::Create(
      NodePointers(nodes), Directors(), Section(), Material());
  ASSERT_TRUE(shell_candidate.HasValue());

  const auto stiffness_candidate = shell_candidate.Value().Stiffness();
  ASSERT_TRUE(stiffness_candidate.HasValue());
  const auto& stiffness = stiffness_candidate.Value();
  EXPECT_EQ(stiffness.physical_local20.Rows(), 20U);
  EXPECT_EQ(stiffness.physical_local20.Columns(), 20U);
  EXPECT_EQ(stiffness.physical_global24.Rows(), 24U);
  EXPECT_EQ(stiffness.drilling_global24.Rows(), 24U);
  EXPECT_EQ(stiffness.stabilized_global24.Rows(), 24U);
  for (const fesa::Matrix* matrix :
       {&stiffness.physical_local20, &stiffness.physical_global24,
        &stiffness.drilling_global24, &stiffness.stabilized_global24}) {
    for (std::size_t row = 0U; row < matrix->Rows(); ++row) {
      for (std::size_t column = 0U; column < matrix->Columns(); ++column) {
        EXPECT_TRUE(std::isfinite((*matrix)(row, column)));
      }
    }
  }
  EXPECT_LE(ScaledSymmetryError(stiffness.physical_local20, 5U, 2.0), 1.0e-12);
  EXPECT_LE(ScaledSymmetryError(stiffness.physical_global24, 6U, 2.0), 1.0e-12);
  EXPECT_LE(ScaledSymmetryError(stiffness.drilling_global24, 6U, 2.0), 1.0e-12);
  EXPECT_LE(ScaledSymmetryError(stiffness.stabilized_global24, 6U, 2.0),
            1.0e-12);

  const auto repeated_candidate = shell_candidate.Value().Stiffness();
  ASSERT_TRUE(repeated_candidate.HasValue());
  const auto& repeated = repeated_candidate.Value();
  ExpectMatrixNear(repeated.physical_local20, stiffness.physical_local20, 0.0);
  ExpectMatrixNear(repeated.physical_global24, stiffness.physical_global24,
                   0.0);
  ExpectMatrixNear(repeated.drilling_global24, stiffness.drilling_global24,
                   0.0);
  ExpectMatrixNear(repeated.stabilized_global24, stiffness.stabilized_global24,
                   0.0);
  EXPECT_DOUBLE_EQ(repeated.drilling_stiffness, stiffness.drilling_stiffness);
}

TEST(Mitc4ShellKernel,
     PreservesExactStiffnessRecoveryAndPatchAcrossVector3Migration) {
  const auto nodes = PlanarNodes();
  const auto shell_candidate = fesa::Mitc4Shell::Create(
      NodePointers(nodes), Directors(), Section(), Material());
  ASSERT_TRUE(shell_candidate.HasValue());
  const auto& shell = shell_candidate.Value();
  const auto stiffness_candidate = shell.Stiffness();
  ASSERT_TRUE(stiffness_candidate.HasValue());
  const auto& stiffness = stiffness_candidate.Value();
  const fesa::Vector3 frame_e3{shell.LocalFrame(0.0, 0.0).e3};
  EXPECT_DOUBLE_EQ(frame_e3.X(), 0.0);
  EXPECT_DOUBLE_EQ(frame_e3.Y(), 0.0);
  EXPECT_DOUBLE_EQ(frame_e3.Z(), 1.0);

  constexpr std::array<double, 8> kGeneralized{0.1,   -0.05, 0.2, 0.3,
                                               -0.15, 0.25,  0.4, -0.3};
  fesa::Vector global_field{24U};
  for (std::size_t node_index = 0U; node_index < nodes.size(); ++node_index) {
    const double x = nodes[node_index].coordinates[0];
    const double y = nodes[node_index].coordinates[1];
    const std::size_t offset = 6U * node_index;
    global_field[offset] = kGeneralized[0] * x + 0.5 * kGeneralized[2] * y;
    global_field[offset + 1U] = kGeneralized[1] * y + 0.5 * kGeneralized[2] * x;
    global_field[offset + 2U] = kGeneralized[6] * x + kGeneralized[7] * y -
                                0.5 * kGeneralized[5] * x * y;
    global_field[offset + 3U] =
        -kGeneralized[4] * y - 0.5 * kGeneralized[5] * x;
    global_field[offset + 4U] = kGeneralized[3] * x + 0.5 * kGeneralized[5] * y;
  }
  const auto recovery_candidate = shell.RecoverPhysical(global_field);
  ASSERT_TRUE(recovery_candidate.HasValue());
  const auto& recovery = recovery_candidate.Value();

  std::array<std::array<double, 5>, 4> e11_values{};
  for (std::size_t node_index = 0U; node_index < nodes.size(); ++node_index) {
    e11_values[node_index][0U] = 0.2 * nodes[node_index].coordinates[0U];
  }
  const double gauss = 1.0 / std::sqrt(3.0);
  const fesa::Vector patch = shell.StrainDisplacement20(gauss, -gauss, gauss)
                                 .Multiply(PhysicalField(e11_values));

  EXPECT_DOUBLE_EQ(stiffness.physical_local20(0U, 0U), 0x1.d555555555554p+6);
  EXPECT_DOUBLE_EQ(stiffness.physical_local20(0U, 4U), 0.0);
  EXPECT_DOUBLE_EQ(stiffness.physical_global24(2U, 2U), 0x1.aaaaaaaaaaaadp+5);
  EXPECT_DOUBLE_EQ(stiffness.drilling_stiffness, 0x1.0d6cffc5beeb5p-4);
  EXPECT_DOUBLE_EQ(recovery.strain_energy, 0x1.20a3d70a3d70bp+6);
  constexpr std::array<double, 8> kExpectedStrain{
      0x1.999999999999ap-4, -0x1.999999999999bp-5, 0x1.9999999999998p-3,
      0x1.3333333333335p-2, -0x1.3333333333335p-3, 0x1.0000000000000p-2,
      0x1.999999999999cp-2, -0x1.3333333333334p-2};
  EXPECT_EQ(recovery.points[0U].generalized_strain, kExpectedStrain);
  EXPECT_DOUBLE_EQ(patch[0U], 0x1.9999999999999p-3);
  for (std::size_t component = 1U; component < patch.Size(); ++component) {
    EXPECT_DOUBLE_EQ(patch[component], 0.0);
  }
}

// MITC4-KERNEL-002
TEST(Mitc4ShellKernel,
     PreservesPhysicalEnergyUnderTwentyToTwentyFourCongruence) {
  const auto nodes = PlanarNodes();
  const auto shell_candidate = fesa::Mitc4Shell::Create(
      NodePointers(nodes), Directors(), Section(), Material());
  ASSERT_TRUE(shell_candidate.HasValue());
  const auto stiffness_candidate = shell_candidate.Value().Stiffness();
  ASSERT_TRUE(stiffness_candidate.HasValue());
  const auto& stiffness = stiffness_candidate.Value();

  fesa::Vector global_field{24U};
  for (std::size_t index = 0U; index < global_field.Size(); ++index) {
    global_field[index] =
        0.125 * static_cast<double>(static_cast<int>(index % 7U) - 3);
  }
  const auto physical_field20 =
      shell_candidate.Value().PhysicalTransformation20().Multiply(global_field);
  const double local_energy =
      QuadraticEnergy(stiffness.physical_local20, physical_field20);
  const double global_energy =
      QuadraticEnergy(stiffness.physical_global24, global_field);
  ASSERT_NE(local_energy, 0.0);
  ASSERT_NE(global_energy, 0.0);
  EXPECT_LE(std::abs(global_energy - local_energy) /
                (std::abs(global_energy) + std::abs(local_energy)),
            1.0e-12);
}

// MITC4-KERNEL-003
TEST(Mitc4ShellKernel,
     RetainsSixRigidModesAndHasExpectedPhysicalAndStabilizedRank) {
  const auto nodes = PlanarNodes();
  const auto shell_candidate = fesa::Mitc4Shell::Create(
      NodePointers(nodes), Directors(), Section(), Material());
  ASSERT_TRUE(shell_candidate.HasValue());
  const auto stiffness_candidate = shell_candidate.Value().Stiffness();
  ASSERT_TRUE(stiffness_candidate.HasValue());
  const auto& stiffness = stiffness_candidate.Value();

  const auto scaled_physical20 =
      ScaledStiffness(stiffness.physical_local20, 5U, 2.0);
  const auto scaled_physical24 =
      ScaledStiffness(stiffness.physical_global24, 6U, 2.0);
  const auto scaled_stabilized24 =
      ScaledStiffness(stiffness.stabilized_global24, 6U, 2.0);
  EXPECT_EQ(NumericalRank(scaled_physical20), 14U);
  EXPECT_EQ(NumericalRank(scaled_stabilized24), 18U);

  const std::array<fesa::Vector, 6> rigid_modes{
      PhysicalRigidMode(nodes, {1.0, 0.0, 0.0}, {}),
      PhysicalRigidMode(nodes, {0.0, 1.0, 0.0}, {}),
      PhysicalRigidMode(nodes, {0.0, 0.0, 1.0}, {}),
      PhysicalRigidMode(nodes, {}, {1.0, 0.0, 0.0}),
      PhysicalRigidMode(nodes, {}, {0.0, 1.0, 0.0}),
      PhysicalRigidMode(nodes, {}, {0.0, 0.0, 1.0})};
  const double physical_norm = SymmetricOperatorNorm(scaled_physical24);
  const double stabilized_norm = SymmetricOperatorNorm(scaled_stabilized24);
  ASSERT_GT(physical_norm, 0.0);
  ASSERT_GT(stabilized_norm, 0.0);
  for (const auto& rigidMode : rigid_modes) {
    fesa::Vector scaledMode = rigidMode;
    for (std::size_t node_index = 0U; node_index < 4U; ++node_index) {
      for (std::size_t component = 0U; component < 3U; ++component) {
        scaledMode[6U * node_index + component] /= 2.0;
      }
    }
    const double modeNorm = scaledMode.Norm();
    ASSERT_GT(modeNorm, 0.0);
    EXPECT_LE(scaled_physical24.Multiply(scaledMode).Norm() /
                  (physical_norm * modeNorm),
              1.0e-10);
    EXPECT_LE(scaled_stabilized24.Multiply(scaledMode).Norm() /
                  (stabilized_norm * modeNorm),
              1.0e-10);
  }
}

// MITC4-KERNEL-004
TEST(Mitc4ShellPatch, ReproducesIndependentMembraneBendingShearAndTwistFields) {
  const auto nodes = PlanarNodes();
  const auto shell_candidate = fesa::Mitc4Shell::Create(
      NodePointers(nodes), Directors(), Section(), Material());
  ASSERT_TRUE(shell_candidate.HasValue());
  const auto& shell = shell_candidate.Value();
  const auto stiffness_candidate = shell.Stiffness();
  ASSERT_TRUE(stiffness_candidate.HasValue());
  const auto& stiffness = stiffness_candidate.Value().physical_local20;
  constexpr double kMagnitude = 0.2;
  const double gauss = 1.0 / std::sqrt(3.0);

  std::array<std::array<double, 5>, 4> e11_values{};
  std::array<std::array<double, 5>, 4> e22_values{};
  std::array<std::array<double, 5>, 4> g12_values{};
  std::array<std::array<double, 5>, 4> k11_values{};
  std::array<std::array<double, 5>, 4> k22_values{};
  std::array<std::array<double, 5>, 4> g13_values{};
  std::array<std::array<double, 5>, 4> g23_values{};
  std::array<std::array<double, 5>, 4> k12_values{};
  for (std::size_t node_index = 0U; node_index < nodes.size(); ++node_index) {
    const double x = nodes[node_index].coordinates[0];
    const double y = nodes[node_index].coordinates[1];
    e11_values[node_index][0] = kMagnitude * x;
    e22_values[node_index][1] = kMagnitude * y;
    g12_values[node_index][0] = 0.5 * kMagnitude * y;
    g12_values[node_index][1] = 0.5 * kMagnitude * x;
    k11_values[node_index][4] = kMagnitude * x;
    k22_values[node_index][3] = -kMagnitude * y;
    g13_values[node_index][2] = kMagnitude * x;
    g23_values[node_index][2] = kMagnitude * y;
    k12_values[node_index][2] = -0.5 * kMagnitude * x * y;
    k12_values[node_index][3] = -0.5 * kMagnitude * x;
    k12_values[node_index][4] = 0.5 * kMagnitude * y;
  }

  const std::array<fesa::Vector, 8> fields{
      PhysicalField(e11_values), PhysicalField(e22_values),
      PhysicalField(g12_values), PhysicalField(k11_values),
      PhysicalField(k22_values), PhysicalField(g13_values),
      PhysicalField(g23_values), PhysicalField(k12_values)};
  ExpectStrain(shell, fields[0], gauss, -gauss, gauss,
               {kMagnitude, 0.0, 0.0, 0.0, 0.0});
  ExpectStrain(shell, fields[1], gauss, -gauss, gauss,
               {0.0, kMagnitude, 0.0, 0.0, 0.0});
  ExpectStrain(shell, fields[2], gauss, -gauss, gauss,
               {0.0, 0.0, kMagnitude, 0.0, 0.0});
  ExpectStrain(shell, fields[3], gauss, -gauss, gauss,
               {gauss * kMagnitude, 0.0, 0.0, 0.0, 0.0});
  ExpectStrain(shell, fields[4], gauss, -gauss, gauss,
               {0.0, gauss * kMagnitude, 0.0, 0.0, 0.0});
  ExpectStrain(shell, fields[5], gauss, -gauss, gauss,
               {0.0, 0.0, 0.0, kMagnitude, 0.0});
  ExpectStrain(shell, fields[6], gauss, -gauss, gauss,
               {0.0, 0.0, 0.0, 0.0, kMagnitude});
  ExpectStrain(shell, fields[7], gauss, -gauss, gauss,
               {0.0, 0.0, gauss * kMagnitude, 0.0, 0.0});
  for (const auto& field : fields) {
    EXPECT_GT(QuadraticEnergy(stiffness, field), 0.0);
  }
}

// MITC4-KERNEL-005
TEST(Mitc4ShellDrilling,
     UsesOnlyEightPositivePhysicalRotationDiagonalsAndFixedFactor) {
  const std::array<fesa::Node, 4> nodes{
      Node(1, {-50.0, -50.0, 0.0}), Node(2, {50.0, -50.0, 0.0}),
      Node(3, {50.0, 50.0, 0.0}), Node(4, {-50.0, 50.0, 0.0})};
  const auto shell_candidate = fesa::Mitc4Shell::Create(
      NodePointers(nodes), Directors(), Section(0.1), Material());
  ASSERT_TRUE(shell_candidate.HasValue());
  const auto stiffness_candidate = shell_candidate.Value().Stiffness();
  ASSERT_TRUE(stiffness_candidate.HasValue());
  const auto& stiffness = stiffness_candidate.Value();

  double expected_reference = (std::numeric_limits<double>::max)();
  double all_diagonal_minimum = (std::numeric_limits<double>::max)();
  for (std::size_t node_index = 0U; node_index < 4U; ++node_index) {
    for (std::size_t rotation = 3U; rotation < 5U; ++rotation) {
      const std::size_t index = 5U * node_index + rotation;
      const double diagonal = stiffness.physical_local20(index, index);
      ASSERT_TRUE(std::isfinite(diagonal));
      ASSERT_GT(diagonal, 0.0);
      expected_reference = (std::min)(expected_reference, diagonal);
    }
    for (std::size_t component = 0U; component < 5U; ++component) {
      const std::size_t index = 5U * node_index + component;
      const double diagonal = stiffness.physical_local20(index, index);
      if (std::isfinite(diagonal) && diagonal > 0.0) {
        all_diagonal_minimum = (std::min)(all_diagonal_minimum, diagonal);
      }
    }
  }
  EXPECT_DOUBLE_EQ(stiffness.drilling_stiffness, 1.0e-3 * expected_reference);
  EXPECT_LT(all_diagonal_minimum, expected_reference);
  EXPECT_NE(stiffness.drilling_stiffness, 1.0e-3 * all_diagonal_minimum);
}

// MITC4-KERNEL-006
TEST(Mitc4ShellDrilling,
     FailsNonfiniteReferenceAndStabilizesEachPureDrillCoordinate) {
  const auto nodes = PlanarNodes();
  const auto shell_candidate = fesa::Mitc4Shell::Create(
      NodePointers(nodes), Directors(), Section(), Material());
  ASSERT_TRUE(shell_candidate.HasValue());
  const auto stiffness_candidate = shell_candidate.Value().Stiffness();
  ASSERT_TRUE(stiffness_candidate.HasValue());
  const auto& stiffness = stiffness_candidate.Value();

  for (std::size_t node_index = 0U; node_index < 4U; ++node_index) {
    fesa::Vector pure_drill{24U};
    pure_drill[6U * node_index + 5U] = 1.0;
    EXPECT_DOUBLE_EQ(stiffness.physical_global24.Multiply(pure_drill).Norm(),
                     0.0);
    const auto drill_action = stiffness.drilling_global24.Multiply(pure_drill);
    EXPECT_DOUBLE_EQ(drill_action[6U * node_index + 5U],
                     stiffness.drilling_stiffness);
    EXPECT_GT(QuadraticEnergy(stiffness.drilling_global24, pure_drill), 0.0);
  }

  const std::array<fesa::Node, 4> extreme_nodes{
      Node(1, {-5.0e9, -5.0e9, 0.0}), Node(2, {5.0e9, -5.0e9, 0.0}),
      Node(3, {5.0e9, 5.0e9, 0.0}), Node(4, {-5.0e9, 5.0e9, 0.0})};
  const auto extreme_shell =
      fesa::Mitc4Shell::Create(NodePointers(extreme_nodes), Directors(),
                               Section(1.0), Material(1.0e300));
  ASSERT_TRUE(extreme_shell.HasValue());
  const auto failure = extreme_shell.Value().Stiffness();
  ASSERT_FALSE(failure.HasValue());
  ASSERT_EQ(failure.GetStatus().Diagnostics().size(), 1U);
  EXPECT_EQ(failure.GetStatus().Diagnostics()[0].code,
            "invalid-shell-stiffness");
  const auto repeated_failure = extreme_shell.Value().Stiffness();
  ASSERT_FALSE(repeated_failure.HasValue());
  ASSERT_EQ(repeated_failure.GetStatus().Diagnostics().size(), 1U);
  EXPECT_EQ(repeated_failure.GetStatus().Diagnostics()[0].code,
            failure.GetStatus().Diagnostics()[0].code);
  EXPECT_EQ(repeated_failure.GetStatus().Diagnostics()[0].message,
            failure.GetStatus().Diagnostics()[0].message);
}

// MITC4-KERNEL-007
TEST(Mitc4ShellDrilling, ExcludesPureDrillFromPhysicalRecoveryAndEnergy) {
  const auto nodes = PlanarNodes();
  const auto shell_candidate = fesa::Mitc4Shell::Create(
      NodePointers(nodes), Directors(), Section(), Material());
  ASSERT_TRUE(shell_candidate.HasValue());
  const auto& shell = shell_candidate.Value();
  const auto stiffness_candidate = shell.Stiffness();
  ASSERT_TRUE(stiffness_candidate.HasValue());

  for (std::size_t node_index = 0U; node_index < nodes.size(); ++node_index) {
    fesa::Vector pure_drill{24U};
    pure_drill[6U * node_index + 5U] = 1.0;
    EXPECT_GT(stiffness_candidate.Value()
                  .stabilized_global24.Multiply(pure_drill)
                  .Norm(),
              0.0);

    const auto recovery_candidate = shell.RecoverPhysical(pure_drill);
    ASSERT_TRUE(recovery_candidate.HasValue());
    const auto& recovery = recovery_candidate.Value();
    EXPECT_DOUBLE_EQ(recovery.strain_energy, 0.0);
    for (const auto& point : recovery.points) {
      for (double value : point.generalized_strain) {
        EXPECT_DOUBLE_EQ(value, 0.0);
      }
      for (double value : point.section_resultant) {
        EXPECT_DOUBLE_EQ(value, 0.0);
      }
      for (const auto& stress : point.in_plane_stress) {
        for (double value : stress) {
          EXPECT_DOUBLE_EQ(value, 0.0);
        }
      }
    }
  }
}

// MITC4-PHYSREC-001
TEST(Mitc4ShellPhysicalRecovery,
     RecoversHandFieldAtFixedLocationsAndSectionPositions) {
  const auto nodes = PlanarNodes();
  const auto shell_candidate = fesa::Mitc4Shell::Create(
      NodePointers(nodes), Directors(), Section(), Material());
  ASSERT_TRUE(shell_candidate.HasValue());
  const auto& shell = shell_candidate.Value();

  constexpr std::array<double, 8> kGeneralized{0.1,   -0.05, 0.2, 0.3,
                                               -0.15, 0.25,  0.4, -0.3};
  fesa::Vector global_field{24U};
  for (std::size_t node_index = 0U; node_index < nodes.size(); ++node_index) {
    const double x = nodes[node_index].coordinates[0];
    const double y = nodes[node_index].coordinates[1];
    const std::size_t offset = 6U * node_index;
    global_field[offset] = kGeneralized[0] * x + 0.5 * kGeneralized[2] * y;
    global_field[offset + 1U] = kGeneralized[1] * y + 0.5 * kGeneralized[2] * x;
    global_field[offset + 2U] = kGeneralized[6] * x + kGeneralized[7] * y -
                                0.5 * kGeneralized[5] * x * y;
    global_field[offset + 3U] =
        -kGeneralized[4] * y - 0.5 * kGeneralized[5] * x;
    global_field[offset + 4U] = kGeneralized[3] * x + 0.5 * kGeneralized[5] * y;
  }

  const auto recovery_candidate = shell.RecoverPhysical(global_field);
  ASSERT_TRUE(recovery_candidate.HasValue());
  const auto& recovery = recovery_candidate.Value();
  const double gauss = 1.0 / std::sqrt(3.0);
  const std::array<std::array<double, 2>, 4> expected_coordinates{
      std::array<double, 2>{-gauss, -gauss},
      std::array<double, 2>{gauss, -gauss}, std::array<double, 2>{gauss, gauss},
      std::array<double, 2>{-gauss, gauss}};
  constexpr std::array<double, 8> kExpectedResultant{22.4, -6.4, 19.2, 22.4,
                                                     -6.4, 8.0,  32.0, -24.0};
  constexpr std::array<std::array<double, 3>, 3> kExpectedStress{
      std::array<double, 3>{-22.4, 6.4, -2.4},
      std::array<double, 3>{11.2, -3.2, 9.6},
      std::array<double, 3>{44.8, -12.8, 21.6}};

  ASSERT_EQ(recovery.points.size(), expected_coordinates.size());
  for (std::size_t point_index = 0U; point_index < recovery.points.size();
       ++point_index) {
    const auto& point = recovery.points[point_index];
    EXPECT_EQ(point.natural_coordinates, expected_coordinates[point_index]);
    ExpectOrthonormalRightHanded(point.local_frame);
    ExpectVectorNear(point.local_frame.e1, {1.0, 0.0, 0.0});
    ExpectVectorNear(point.local_frame.e2, {0.0, 1.0, 0.0});
    ExpectVectorNear(point.local_frame.e3, {0.0, 0.0, 1.0});
    for (std::size_t component = 0U; component < kGeneralized.size();
         ++component) {
      EXPECT_NEAR(point.generalized_strain[component], kGeneralized[component],
                  1.0e-12);
      EXPECT_NEAR(point.section_resultant[component],
                  kExpectedResultant[component], 1.0e-12);
    }
    for (std::size_t position = 0U; position < kExpectedStress.size();
         ++position) {
      for (std::size_t component = 0U;
           component < kExpectedStress[position].size(); ++component) {
        EXPECT_NEAR(point.in_plane_stress[position][component],
                    kExpectedStress[position][component], 1.0e-12);
      }
    }
  }
  EXPECT_NEAR(recovery.strain_energy, 72.16, 1.0e-12);
}