#include "fesa/elements/euler_beam_3d.h"

#include <gtest/gtest.h>

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

#include "fesa/math/vector3.h"
#include "fesa/model/model_types.h"

namespace fesa {
namespace {

constexpr std::size_t kElementDofCount = 12U;
constexpr double kMatrixTolerance = 1.0e-12;
constexpr double kRigidTolerance = 1.0e-10;
constexpr double kAnalyticalTolerance = 1.0e-9;

Node MakeNode(std::array<double, 3> coordinates, std::size_t line) {
  return {{"Beam-1", static_cast<std::int64_t>(line), std::to_string(line)},
          coordinates,
          {"beam-test.inp", line}};
}

LinearElasticMaterial MakeMaterial(double youngs_modulus = 210.0e9,
                                   double poisson_ratio = 0.3) {
  return {"Steel", youngs_modulus, poisson_ratio, {"beam-test.inp", 20U}};
}

GeneralBeamSection MakeSection(
    std::array<double, 3> first_axis = {0.0, 1.0, 0.0},
    std::vector<std::array<double, 2>> section_points = {}) {
  return {"Section-1",
          0.012,
          2.5e-5,
          0.0,
          4.0e-5,
          1.5e-5,
          first_axis,
          std::move(section_points),
          {"beam-test.inp", 30U}};
}

EulerBeam3D RequireBeam(const Node& first_node, const Node& second_node,
                        const GeneralBeamSection& section,
                        const LinearElasticMaterial& material) {
  auto result = EulerBeam3D::Create(first_node, second_node, section, material);
  if (!result.HasValue()) {
    throw std::runtime_error{"Expected a valid EulerBeam3D fixture."};
  }
  return std::move(result.Value());
}

EulerBeam3D AlignedBeam(
    double length, const GeneralBeamSection& section = MakeSection(),
    const LinearElasticMaterial& material = MakeMaterial()) {
  return RequireBeam(MakeNode({0.0, 0.0, 0.0}, 1U),
                     MakeNode({length, 0.0, 0.0}, 2U), section, material);
}

double MaximumAbsoluteEntry(const Matrix& matrix) {
  double maximum = 0.0;
  for (std::size_t row = 0; row < matrix.Rows(); ++row) {
    for (std::size_t column = 0; column < matrix.Columns(); ++column) {
      maximum = (std::max)(maximum, std::abs(matrix(row, column)));
    }
  }
  return maximum;
}

bool MatrixIsFinite(const Matrix& matrix) {
  for (std::size_t row = 0; row < matrix.Rows(); ++row) {
    for (std::size_t column = 0; column < matrix.Columns(); ++column) {
      if (!std::isfinite(matrix(row, column))) {
        return false;
      }
    }
  }
  return true;
}

double NormalizedMatrixError(const Matrix& actual, const Matrix& expected) {
  if (actual.Rows() != expected.Rows() ||
      actual.Columns() != expected.Columns()) {
    throw std::invalid_argument{"Matrix comparison requires equal shapes."};
  }

  double maximum_difference = 0.0;
  for (std::size_t row = 0; row < actual.Rows(); ++row) {
    for (std::size_t column = 0; column < actual.Columns(); ++column) {
      maximum_difference =
          (std::max)(maximum_difference,
                     std::abs(actual(row, column) - expected(row, column)));
    }
  }

  const double scale =
      (std::max)(1.0, (std::max)(MaximumAbsoluteEntry(actual),
                                 MaximumAbsoluteEntry(expected)));
  return maximum_difference / scale;
}

double VectorNorm(const Vector& vector) {
  double sum = 0.0;
  for (std::size_t index = 0; index < vector.Size(); ++index) {
    sum += vector[index] * vector[index];
  }
  return std::sqrt(sum);
}

double QuadraticEnergy(const Matrix& matrix, const Vector& vector) {
  const Vector product = matrix.Multiply(vector);
  double value = 0.0;
  for (std::size_t index = 0; index < vector.Size(); ++index) {
    value += vector[index] * product[index];
  }
  return value;
}

void ExpectScaledNear(double actual, double expected,
                      double relative_tolerance) {
  const double scale = (std::max)(1.0, std::abs(expected));
  EXPECT_LE(std::abs(actual - expected), relative_tolerance * scale);
}

void ExpectRelativeNear(double actual, double expected,
                        double relative_tolerance) {
  ASSERT_NE(expected, 0.0);
  EXPECT_LE(std::abs(actual - expected) / std::abs(expected),
            relative_tolerance);
}

Matrix ExpectedClosedStiffness(double length, const GeneralBeamSection& section,
                               const LinearElasticMaterial& material) {
  Matrix expected{kElementDofCount, kElementDofCount};
  const double shear_modulus =
      material.youngs_modulus / (2.0 * (1.0 + material.poisson_ratio));

  const auto add_block = [&expected](const std::vector<std::size_t>& indices,
                                     const std::vector<double>& values) {
    const std::size_t width = indices.size();
    for (std::size_t row = 0; row < width; ++row) {
      for (std::size_t column = 0; column < width; ++column) {
        expected(indices[row], indices[column]) = values[row * width + column];
      }
    }
  };

  const double axial = material.youngs_modulus * section.area / length;
  add_block({0U, 6U}, {axial, -axial, -axial, axial});

  const double torsion = shear_modulus * section.torsional_constant / length;
  add_block({3U, 9U}, {torsion, -torsion, -torsion, torsion});

  const auto bending_block = [length](double flexural_rigidity,
                                      double rotation_sign) {
    const double v = 12.0 * flexural_rigidity / (length * length * length);
    const double c =
        rotation_sign * 6.0 * flexural_rigidity / (length * length);
    const double d = 4.0 * flexural_rigidity / length;
    const double e = 2.0 * flexural_rigidity / length;
    return std::vector<double>{v,  c,  -v, c,  c, d, -c, e,
                               -v, -c, v,  -c, c, e, -c, d};
  };

  add_block({1U, 5U, 7U, 11U},
            bending_block(material.youngs_modulus * section.i22, 1.0));
  add_block({2U, 4U, 8U, 10U},
            bending_block(material.youngs_modulus * section.i11, -1.0));
  return expected;
}

std::array<double, kElementDofCount> SymmetricEigenvalues(Matrix matrix) {
  for (std::size_t iteration = 0; iteration < 100U * kElementDofCount;
       ++iteration) {
    std::size_t p = 0U;
    std::size_t q = 1U;
    double maximum_off_diagonal = 0.0;
    for (std::size_t row = 0; row < kElementDofCount; ++row) {
      for (std::size_t column = row + 1U; column < kElementDofCount; ++column) {
        const double candidate = std::abs(matrix(row, column));
        if (candidate > maximum_off_diagonal) {
          maximum_off_diagonal = candidate;
          p = row;
          q = column;
        }
      }
    }
    if (maximum_off_diagonal <=
        1.0e-14 * (std::max)(1.0, MaximumAbsoluteEntry(matrix))) {
      break;
    }

    const double app = matrix(p, p);
    const double aqq = matrix(q, q);
    const double apq = matrix(p, q);
    const double angle = 0.5 * std::atan2(2.0 * apq, aqq - app);
    const double cosine = std::cos(angle);
    const double sine = std::sin(angle);

    for (std::size_t index = 0; index < kElementDofCount; ++index) {
      if (index == p || index == q) {
        continue;
      }
      const double aip = matrix(index, p);
      const double aiq = matrix(index, q);
      matrix(index, p) = cosine * aip - sine * aiq;
      matrix(p, index) = matrix(index, p);
      matrix(index, q) = sine * aip + cosine * aiq;
      matrix(q, index) = matrix(index, q);
    }

    matrix(p, p) =
        cosine * cosine * app - 2.0 * sine * cosine * apq + sine * sine * aqq;
    matrix(q, q) =
        sine * sine * app + 2.0 * sine * cosine * apq + cosine * cosine * aqq;
    matrix(p, q) = 0.0;
    matrix(q, p) = 0.0;
  }

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

std::size_t SymmetricRank(const Matrix& matrix, double relative_tolerance) {
  const auto eigenvalues = SymmetricEigenvalues(matrix);
  double maximum = 0.0;
  for (const double value : eigenvalues) {
    maximum = (std::max)(maximum, std::abs(value));
  }
  return static_cast<std::size_t>(
      std::count_if(eigenvalues.begin(), eigenvalues.end(),
                    [maximum, relative_tolerance](double value) {
                      return std::abs(value) > relative_tolerance * maximum;
                    }));
}

Matrix TestOnlyOnePointStiffness(double length,
                                 const GeneralBeamSection& section,
                                 const LinearElasticMaterial& material) {
  // At xi=0 the two bending curvature rows retain only the nodal rotations.
  // This deliberately under-integrated negative control is independent of
  // production.
  Matrix b{4U, kElementDofCount};
  b(0U, 0U) = -1.0 / length;
  b(0U, 6U) = 1.0 / length;
  b(1U, 3U) = -1.0 / length;
  b(1U, 9U) = 1.0 / length;
  b(2U, 4U) = -1.0 / length;
  b(2U, 10U) = 1.0 / length;
  b(3U, 5U) = -1.0 / length;
  b(3U, 11U) = 1.0 / length;

  const double shear_modulus =
      material.youngs_modulus / (2.0 * (1.0 + material.poisson_ratio));
  const std::array<double, 4> diagonal = {
      material.youngs_modulus * section.area,
      shear_modulus * section.torsional_constant,
      material.youngs_modulus * section.i11,
      material.youngs_modulus * section.i22};

  Matrix stiffness{kElementDofCount, kElementDofCount};
  for (std::size_t row = 0; row < kElementDofCount; ++row) {
    for (std::size_t column = 0; column < kElementDofCount; ++column) {
      for (std::size_t component = 0; component < diagonal.size();
           ++component) {
        stiffness(row, column) += b(component, row) * diagonal[component] *
                                  b(component, column) * length;
      }
    }
  }
  return stiffness;
}

Vector SolveFixedFirstNode(const Matrix& stiffness,
                           const std::array<double, 6>& free_end_load) {
  std::array<std::array<double, 7>, 6> augmented{};
  for (std::size_t row = 0; row < 6U; ++row) {
    for (std::size_t column = 0; column < 6U; ++column) {
      augmented[row][column] = stiffness(row + 6U, column + 6U);
    }
    augmented[row][6U] = free_end_load[row];
  }

  for (std::size_t pivot = 0; pivot < 6U; ++pivot) {
    std::size_t pivot_row = pivot;
    for (std::size_t row = pivot + 1U; row < 6U; ++row) {
      if (std::abs(augmented[row][pivot]) >
          std::abs(augmented[pivot_row][pivot])) {
        pivot_row = row;
      }
    }
    if (std::abs(augmented[pivot_row][pivot]) <=
        std::numeric_limits<double>::min()) {
      throw std::runtime_error{"Cantilever fixture is singular."};
    }
    std::swap(augmented[pivot], augmented[pivot_row]);

    const double pivot_value = augmented[pivot][pivot];
    for (std::size_t column = pivot; column < 7U; ++column) {
      augmented[pivot][column] /= pivot_value;
    }
    for (std::size_t row = 0; row < 6U; ++row) {
      if (row == pivot) {
        continue;
      }
      const double factor = augmented[row][pivot];
      for (std::size_t column = pivot; column < 7U; ++column) {
        augmented[row][column] -= factor * augmented[pivot][column];
      }
    }
  }

  Vector displacement{kElementDofCount};
  for (std::size_t component = 0; component < 6U; ++component) {
    displacement[component + 6U] = augmented[component][6U];
  }
  return displacement;
}

Vector SolveDenseSystem(Matrix matrix, Vector right_hand_side) {
  if (matrix.Rows() != matrix.Columns() ||
      matrix.Rows() != right_hand_side.Size()) {
    throw std::invalid_argument{"Dense test solve requires a square system."};
  }

  for (std::size_t pivot = 0; pivot < matrix.Rows(); ++pivot) {
    std::size_t pivot_row = pivot;
    for (std::size_t row = pivot + 1U; row < matrix.Rows(); ++row) {
      if (std::abs(matrix(row, pivot)) > std::abs(matrix(pivot_row, pivot))) {
        pivot_row = row;
      }
    }
    if (!(std::abs(matrix(pivot_row, pivot)) > 0.0) ||
        !std::isfinite(matrix(pivot_row, pivot))) {
      throw std::runtime_error{"Uniform-load test fixture is singular."};
    }
    for (std::size_t column = pivot; column < matrix.Columns(); ++column) {
      std::swap(matrix(pivot, column), matrix(pivot_row, column));
    }
    std::swap(right_hand_side[pivot], right_hand_side[pivot_row]);

    const double pivot_value = matrix(pivot, pivot);
    for (std::size_t column = pivot; column < matrix.Columns(); ++column) {
      matrix(pivot, column) /= pivot_value;
    }
    right_hand_side[pivot] /= pivot_value;
    for (std::size_t row = 0; row < matrix.Rows(); ++row) {
      if (row == pivot) {
        continue;
      }
      const double factor = matrix(row, pivot);
      for (std::size_t column = pivot; column < matrix.Columns(); ++column) {
        matrix(row, column) -= factor * matrix(pivot, column);
      }
      right_hand_side[row] -= factor * right_hand_side[pivot];
    }
  }
  return right_hand_side;
}

Vector SolveUniformTransverseCantilever(std::size_t element_count,
                                        double length, double line_load,
                                        const GeneralBeamSection& section,
                                        const LinearElasticMaterial& material) {
  const double element_length = length / static_cast<double>(element_count);
  const std::size_t system_size = 2U * (element_count + 1U);
  Matrix assembled_stiffness{system_size, system_size};
  Vector assembled_load{system_size};
  const std::array<std::size_t, 4> bending_dofs = {1U, 5U, 7U, 11U};

  // Test-only direct assembly keeps this evidence at the formulation boundary:
  // two [v,rz] DOFs per node, with no Domain, parser, or DLOAD path.
  for (std::size_t element = 0; element < element_count; ++element) {
    const EulerBeam3D beam = AlignedBeam(element_length, section, material);
    const Matrix element_stiffness = beam.LocalStiffness();
    const Vector element_load =
        beam.LocalEquivalentLoad({0.0, line_load, 0.0, 0.0});
    const std::array<std::size_t, 4> assembled_dofs = {
        2U * element, 2U * element + 1U, 2U * (element + 1U),
        2U * (element + 1U) + 1U};
    for (std::size_t row = 0; row < bending_dofs.size(); ++row) {
      assembled_load[assembled_dofs[row]] += element_load[bending_dofs[row]];
      for (std::size_t column = 0; column < bending_dofs.size(); ++column) {
        assembled_stiffness(assembled_dofs[row], assembled_dofs[column]) +=
            element_stiffness(bending_dofs[row], bending_dofs[column]);
      }
    }
  }

  const std::size_t free_size = system_size - 2U;
  Matrix free_stiffness{free_size, free_size};
  Vector free_load{free_size};
  for (std::size_t row = 0; row < free_size; ++row) {
    free_load[row] = assembled_load[row + 2U];
    for (std::size_t column = 0; column < free_size; ++column) {
      free_stiffness(row, column) = assembled_stiffness(row + 2U, column + 2U);
    }
  }

  const Vector free_displacement =
      SolveDenseSystem(std::move(free_stiffness), std::move(free_load));
  Vector nodal_displacement{system_size};
  for (std::size_t dof = 0; dof < free_size; ++dof) {
    nodal_displacement[dof + 2U] = free_displacement[dof];
  }
  return nodal_displacement;
}

double UniformLoadInteriorDisplacementError(const Vector& nodal_displacement,
                                            std::size_t element_count,
                                            double length, double line_load,
                                            double flexural_rigidity) {
  const double element_length = length / static_cast<double>(element_count);
  const std::array<double, 5> gauss_points = {
      -0.9061798459386640, -0.5384693101056831, 0.0, 0.5384693101056831,
      0.9061798459386640};
  const std::array<double, 5> gauss_weights = {
      0.2369268850561891, 0.4786286704993665, 0.5688888888888889,
      0.4786286704993665, 0.2369268850561891};
  double squared_error = 0.0;
  double squared_reference = 0.0;

  // Five-point integration is independent of production and exactly integrates
  // the squared error between cubic Hermite interpolation and the quartic beam
  // solution.
  for (std::size_t element = 0; element < element_count; ++element) {
    for (std::size_t point = 0; point < gauss_points.size(); ++point) {
      const double r = 0.5 * (1.0 + gauss_points[point]);
      const double r_squared = r * r;
      const double r_cubed = r_squared * r;
      const double h1 = 1.0 - 3.0 * r_squared + 2.0 * r_cubed;
      const double h2 = element_length * (r - 2.0 * r_squared + r_cubed);
      const double h3 = 3.0 * r_squared - 2.0 * r_cubed;
      const double h4 = element_length * (-r_squared + r_cubed);
      const double interpolated =
          h1 * nodal_displacement[2U * element] +
          h2 * nodal_displacement[2U * element + 1U] +
          h3 * nodal_displacement[2U * (element + 1U)] +
          h4 * nodal_displacement[2U * (element + 1U) + 1U];
      const double x = element_length * (static_cast<double>(element) + r);
      const double analytical =
          line_load * x * x *
          (6.0 * length * length - 4.0 * length * x + x * x) /
          (24.0 * flexural_rigidity);
      const double weight = 0.5 * element_length * gauss_weights[point];
      const double difference = interpolated - analytical;
      squared_error += weight * difference * difference;
      squared_reference += weight * analytical * analytical;
    }
  }
  return std::sqrt(squared_error / squared_reference);
}

Matrix TransformationFromKnownRows(
    const std::array<std::array<double, 3>, 3>& rotation) {
  Matrix transformation{kElementDofCount, kElementDofCount};
  for (std::size_t block = 0; block < 4U; ++block) {
    for (std::size_t row = 0; row < 3U; ++row) {
      for (std::size_t column = 0; column < 3U; ++column) {
        transformation(block * 3U + row, block * 3U + column) =
            rotation[row][column];
      }
    }
  }
  return transformation;
}

Vector TransposeMultiply(const Matrix& matrix, const Vector& vector) {
  if (matrix.Rows() != vector.Size()) {
    throw std::invalid_argument{"Transpose multiply dimension mismatch."};
  }
  Vector result{matrix.Columns()};
  for (std::size_t column = 0; column < matrix.Columns(); ++column) {
    for (std::size_t row = 0; row < matrix.Rows(); ++row) {
      result[column] += matrix(row, column) * vector[row];
    }
  }
  return result;
}

double Determinant(const std::array<std::array<double, 3>, 3>& matrix) {
  return matrix[0][0] *
             (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1]) -
         matrix[0][1] *
             (matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0]) +
         matrix[0][2] *
             (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]);
}

TEST(EulerBeam3D, HermiteAndBMatrixMatchReviewedSigns) {
  const double length = 2.5;
  const auto section = MakeSection();
  const auto material = MakeMaterial();
  const auto beam = AlignedBeam(length, section, material);

  const double axial_strain = 0.012;
  const double twist = -0.021;
  const std::array<double, 4> v = {1.2, -0.4, 0.3, -0.07};
  const std::array<double, 4> w = {-0.8, 0.6, -0.2, 0.05};
  const auto value = [](const std::array<double, 4>& coefficients, double x) {
    return coefficients[0] + coefficients[1] * x + coefficients[2] * x * x +
           coefficients[3] * x * x * x;
  };
  const auto slope = [](const std::array<double, 4>& coefficients, double x) {
    return coefficients[1] + 2.0 * coefficients[2] * x +
           3.0 * coefficients[3] * x * x;
  };
  const auto curvature = [](const std::array<double, 4>& coefficients,
                            double x) {
    return 2.0 * coefficients[2] + 6.0 * coefficients[3] * x;
  };

  Vector displacement{kElementDofCount};
  displacement[0U] = 0.2;
  displacement[1U] = value(v, 0.0);
  displacement[2U] = value(w, 0.0);
  displacement[3U] = -0.1;
  displacement[4U] = -slope(w, 0.0);
  displacement[5U] = slope(v, 0.0);
  displacement[6U] = displacement[0U] + axial_strain * length;
  displacement[7U] = value(v, length);
  displacement[8U] = value(w, length);
  displacement[9U] = displacement[3U] + twist * length;
  displacement[10U] = -slope(w, length);
  displacement[11U] = slope(v, length);

  const BeamRecovery recovery = beam.RecoverBeam(displacement);
  const double inverse_sqrt_three = 1.0 / std::sqrt(3.0);
  const std::array<double, 2> gauss_xi = {-inverse_sqrt_three,
                                          inverse_sqrt_three};
  for (std::size_t point = 0; point < gauss_xi.size(); ++point) {
    const double x = 0.5 * length * (1.0 + gauss_xi[point]);
    EXPECT_NEAR(recovery.gauss_generalized_strains[point][0U], axial_strain,
                1.0e-14);
    EXPECT_NEAR(recovery.gauss_generalized_strains[point][1U], twist, 1.0e-14);
    EXPECT_NEAR(recovery.gauss_generalized_strains[point][2U], -curvature(w, x),
                1.0e-13);
    EXPECT_NEAR(recovery.gauss_generalized_strains[point][3U], curvature(v, x),
                1.0e-13);
  }

  EXPECT_NEAR(recovery.endpoint_section_resultants[0U][2U],
              material.youngs_modulus * section.i11 * -curvature(w, 0.0),
              1.0e-5);
  EXPECT_NEAR(recovery.endpoint_section_resultants[1U][2U],
              material.youngs_modulus * section.i11 * -curvature(w, length),
              1.0e-5);
  EXPECT_NEAR(recovery.endpoint_section_resultants[0U][3U],
              material.youngs_modulus * section.i22 * curvature(v, 0.0),
              1.0e-5);
  EXPECT_NEAR(recovery.endpoint_section_resultants[1U][3U],
              material.youngs_modulus * section.i22 * curvature(v, length),
              1.0e-5);
}

TEST(EulerBeam3D, TwoPointGaussMatchesClosedStiffness) {
  const double length = 3.7;
  auto section = MakeSection();
  section.area = 0.019;
  section.i11 = 3.1e-5;
  section.i22 = 7.4e-5;
  section.torsional_constant = 2.2e-5;
  const auto material = MakeMaterial(73.0e9, 0.27);
  const auto beam = AlignedBeam(length, section, material);

  const Matrix actual = beam.LocalStiffness();
  const Matrix closed = ExpectedClosedStiffness(length, section, material);
  EXPECT_LE(NormalizedMatrixError(actual, closed), kMatrixTolerance);

  Matrix transpose{actual.Rows(), actual.Columns()};
  for (std::size_t row = 0; row < actual.Rows(); ++row) {
    for (std::size_t column = 0; column < actual.Columns(); ++column) {
      transpose(row, column) = actual(column, row);
    }
  }
  EXPECT_LE(NormalizedMatrixError(actual, transpose), kMatrixTolerance);
}

TEST(EulerBeam3D, HasSixRigidModesRankSixAndPositiveDeformationEnergy) {
  const double length = 2.0;
  auto section = MakeSection();
  section.area = 1.4;
  section.i11 = 0.8;
  section.i22 = 1.1;
  section.torsional_constant = 0.6;
  const auto material = MakeMaterial(5.0, 0.25);
  const Matrix stiffness =
      AlignedBeam(length, section, material).LocalStiffness();

  std::array<Vector, 6> rigid_modes = {
      Vector{kElementDofCount}, Vector{kElementDofCount},
      Vector{kElementDofCount}, Vector{kElementDofCount},
      Vector{kElementDofCount}, Vector{kElementDofCount}};
  rigid_modes[0U][0U] = rigid_modes[0U][6U] = 1.0;
  rigid_modes[1U][1U] = rigid_modes[1U][7U] = 1.0;
  rigid_modes[2U][2U] = rigid_modes[2U][8U] = 1.0;
  rigid_modes[3U][3U] = rigid_modes[3U][9U] = 1.0;
  rigid_modes[4U][4U] = rigid_modes[4U][10U] = 1.0;
  rigid_modes[4U][8U] = -length;
  rigid_modes[5U][5U] = rigid_modes[5U][11U] = 1.0;
  rigid_modes[5U][7U] = length;

  const double stiffness_scale =
      (std::max)(1.0, MaximumAbsoluteEntry(stiffness));
  for (const Vector& mode : rigid_modes) {
    const double normalizedResidual =
        VectorNorm(stiffness.Multiply(mode)) /
        (stiffness_scale * (std::max)(1.0, VectorNorm(mode)));
    EXPECT_LE(normalizedResidual, kRigidTolerance);
  }

  const std::array<double, kElementDofCount> q = {
      1.0, 1.0, 1.0, length, length, length,
      1.0, 1.0, 1.0, length, length, length};
  Matrix scaled{kElementDofCount, kElementDofCount};
  for (std::size_t row = 0; row < kElementDofCount; ++row) {
    for (std::size_t column = 0; column < kElementDofCount; ++column) {
      scaled(row, column) = stiffness(row, column) / (q[row] * q[column]);
    }
  }
  const auto eigenvalues = SymmetricEigenvalues(scaled);
  double maximum_singular_value = 0.0;
  for (const double value : eigenvalues) {
    maximum_singular_value =
        (std::max)(maximum_singular_value, std::abs(value));
  }
  const auto positive_count = std::count_if(
      eigenvalues.begin(), eigenvalues.end(),
      [maximum_singular_value](double value) {
        return std::abs(value) > kRigidTolerance * maximum_singular_value;
      });
  EXPECT_EQ(positive_count, 6);
  for (const double value : eigenvalues) {
    EXPECT_GE(value, -kRigidTolerance * maximum_singular_value);
  }

  for (std::size_t component = 0; component < 6U; ++component) {
    Vector deformation{kElementDofCount};
    deformation[6U + component] = 1.0;
    EXPECT_GT(QuadraticEnergy(stiffness, deformation), 0.0);
  }
}

TEST(EulerBeam3D, RotatedTransformPreservesWorkAndEnergy) {
  const double inverse_sqrt_two = 1.0 / std::sqrt(2.0);
  const std::array<std::array<double, 3>, 3> rotation = {
      {{{2.0 / 3.0, 2.0 / 3.0, 1.0 / 3.0}},
       {{-inverse_sqrt_two, inverse_sqrt_two, 0.0}},
       {{-inverse_sqrt_two / 3.0, -inverse_sqrt_two / 3.0,
         4.0 * inverse_sqrt_two / 3.0}}}};
  const Matrix transformation = TransformationFromKnownRows(rotation);

  for (std::size_t row = 0; row < 3U; ++row) {
    for (std::size_t column = 0; column < 3U; ++column) {
      double dot = 0.0;
      for (std::size_t component = 0; component < 3U; ++component) {
        dot += rotation[row][component] * rotation[column][component];
      }
      EXPECT_NEAR(dot, row == column ? 1.0 : 0.0, 1.0e-14);
    }
  }
  EXPECT_NEAR(Determinant(rotation), 1.0, 1.0e-14);

  const auto section = MakeSection({-2.0, 2.0, 0.0});
  const auto material = MakeMaterial();
  const auto beam =
      RequireBeam(MakeNode({1.0, -2.0, 0.5}, 1U), MakeNode({3.0, 0.0, 1.5}, 2U),
                  section, material);
  const Matrix local = beam.LocalStiffness();
  const Matrix global = beam.GlobalStiffness();

  Matrix expected_global{kElementDofCount, kElementDofCount};
  const Matrix local_times_transform = local.Multiply(transformation);
  for (std::size_t row = 0; row < kElementDofCount; ++row) {
    for (std::size_t column = 0; column < kElementDofCount; ++column) {
      for (std::size_t inner = 0; inner < kElementDofCount; ++inner) {
        expected_global(row, column) +=
            transformation(inner, row) * local_times_transform(inner, column);
      }
    }
  }
  EXPECT_LE(NormalizedMatrixError(global, expected_global), kMatrixTolerance);

  Vector local_displacement{kElementDofCount};
  for (std::size_t index = 0; index < local_displacement.Size(); ++index) {
    local_displacement[index] = 0.01 * static_cast<double>(index + 1U) - 0.04;
  }
  const Vector global_displacement =
      TransposeMultiply(transformation, local_displacement);
  const Vector local_force = local.Multiply(local_displacement);
  const Vector global_force = global.Multiply(global_displacement);
  const Vector expected_global_force =
      TransposeMultiply(transformation, local_force);
  for (std::size_t index = 0; index < kElementDofCount; ++index) {
    ExpectScaledNear(global_force[index], expected_global_force[index],
                     kMatrixTolerance);
  }
  ExpectScaledNear(QuadraticEnergy(global, global_displacement),
                   QuadraticEnergy(local, local_displacement),
                   kMatrixTolerance);

  Vector global_variation{kElementDofCount};
  for (std::size_t index = 0; index < global_variation.Size(); ++index) {
    global_variation[index] = 0.03 - 0.002 * static_cast<double>(index);
  }
  const Vector local_variation = transformation.Multiply(global_variation);
  ExpectScaledNear(global_variation.Dot(global_force),
                   local_variation.Dot(local_force), kMatrixTolerance);

  const BeamRecovery recovery = beam.RecoverBeam(global_displacement);
  EXPECT_NEAR(recovery.gauss_generalized_strains[0U][0U],
              (local_displacement[6U] - local_displacement[0U]) / 3.0, 1.0e-14);
}

TEST(EulerBeam3D, PreservesExactRotatedResultsAcrossVector3Migration) {
  const auto section = MakeSection({-2.0, 2.0, 0.0});
  const auto material = MakeMaterial();
  const Node first_node = MakeNode({1.0, -2.0, 0.5}, 1U);
  const Node second_node = MakeNode({3.0, 0.0, 1.5}, 2U);
  const Vector3 delta =
      Vector3(second_node.coordinates) - Vector3(first_node.coordinates);
  EXPECT_DOUBLE_EQ(delta.Norm(), 3.0);

  const auto beam = RequireBeam(first_node, second_node, section, material);
  const Matrix global = beam.GlobalStiffness();
  Vector displacement{kElementDofCount};
  for (std::size_t index = 0U; index < displacement.Size(); ++index) {
    displacement[index] = 0.01 * static_cast<double>(index + 1U) - 0.04;
  }
  const BeamRecovery recovery = beam.RecoverBeam(displacement);

  EXPECT_DOUBLE_EQ(global(0U, 0U), 0x1.65f135da12f68p+28);
  EXPECT_DOUBLE_EQ(global(0U, 1U), 0x1.6261c084bda12p+28);
  EXPECT_DOUBLE_EQ(global(0U, 2U), 0x1.630ca684bda13p+27);
  EXPECT_DOUBLE_EQ(global(4U, 4U), 0x1.068e359b59b58p+22);
  EXPECT_DOUBLE_EQ(global(5U, 11U), 0x1.2d14a7ee7ee7bp+22);
  EXPECT_DOUBLE_EQ(recovery.gauss_generalized_strains[0U][0U],
                   0x1.1111111111110p-5);
  EXPECT_DOUBLE_EQ(recovery.gauss_generalized_strains[0U][1U],
                   0x1.1111111111111p-5);
  EXPECT_DOUBLE_EQ(recovery.gauss_generalized_strains[0U][2U],
                   -0x1.382425a7d2473p-6);
  EXPECT_DOUBLE_EQ(recovery.gauss_generalized_strains[0U][3U],
                   -0x1.a9389137b051dp-6);
  EXPECT_DOUBLE_EQ(recovery.endpoint_section_resultants[1U][2U],
                   0x1.525c94a87359fp+17);
}

TEST(EulerBeam3D, ConstantLineLoadMatchesAllSignedComponents) {
  const double length = 4.0;
  const ConstantLocalLineLoad load{2.5, -3.0, 5.5, -7.0};
  const Vector equivalent = AlignedBeam(length).LocalEquivalentLoad(load);
  const std::array<double, kElementDofCount> expected = {
      5.0, -6.0, 11.0, -14.0, -22.0 / 3.0, -4.0,
      5.0, -6.0, 11.0, -14.0, 22.0 / 3.0,  4.0};
  ASSERT_EQ(equivalent.Size(), expected.size());
  for (std::size_t index = 0; index < expected.size(); ++index) {
    ExpectScaledNear(equivalent[index], expected[index], kMatrixTolerance);
  }

  auto convergence_section = MakeSection();
  convergence_section.area = 1.0;
  convergence_section.i11 = 1.0;
  convergence_section.i22 = 1.0;
  convergence_section.torsional_constant = 1.0;
  const auto convergence_material = MakeMaterial(5.0, 0.25);
  const double transverse_load = -3.0;
  const std::array<std::size_t, 3> element_counts = {1U, 2U, 4U};
  std::array<double, 3> relative_errors{};
  for (std::size_t mesh = 0; mesh < element_counts.size(); ++mesh) {
    const double element_length =
        length / static_cast<double>(element_counts[mesh]);
    const Vector element_load =
        AlignedBeam(element_length, convergence_section, convergence_material)
            .LocalEquivalentLoad({0.0, transverse_load, 0.0, 0.0});
    ExpectRelativeNear(element_load[1U], transverse_load * element_length / 2.0,
                       kMatrixTolerance);
    ExpectRelativeNear(element_load[5U],
                       transverse_load * element_length * element_length / 12.0,
                       kMatrixTolerance);
    ExpectRelativeNear(element_load[7U], transverse_load * element_length / 2.0,
                       kMatrixTolerance);
    ExpectRelativeNear(
        element_load[11U],
        -transverse_load * element_length * element_length / 12.0,
        kMatrixTolerance);

    const Vector nodal_displacement = SolveUniformTransverseCantilever(
        element_counts[mesh], length, transverse_load, convergence_section,
        convergence_material);
    relative_errors[mesh] = UniformLoadInteriorDisplacementError(
        nodal_displacement, element_counts[mesh], length, transverse_load,
        convergence_material.youngs_modulus * convergence_section.i22);
  }
  EXPECT_GT(relative_errors[0U], relative_errors[1U]);
  EXPECT_GT(relative_errors[1U], relative_errors[2U]);
  EXPECT_NEAR(
      std::log(relative_errors[0U] / relative_errors[1U]) / std::log(2.0), 4.0,
      1.0e-8);
  EXPECT_NEAR(
      std::log(relative_errors[1U] / relative_errors[2U]) / std::log(2.0), 4.0,
      1.0e-8);
}

TEST(EulerBeam3D, AnalyticalAxialTorsionAndTwoPlaneBendingRecover) {
  const double length = 3.0;
  const auto section = MakeSection();
  const auto material = MakeMaterial();
  const auto beam = AlignedBeam(length, section, material);
  const Matrix stiffness = beam.LocalStiffness();
  const double shear_modulus =
      material.youngs_modulus / (2.0 * (1.0 + material.poisson_ratio));

  const double axial_force = 1250.0;
  const Vector axial =
      SolveFixedFirstNode(stiffness, {axial_force, 0.0, 0.0, 0.0, 0.0, 0.0});
  ExpectRelativeNear(
      axial[6U],
      axial_force * length / (material.youngs_modulus * section.area),
      kAnalyticalTolerance);
  const BeamRecovery axial_recovery = beam.RecoverBeam(axial);
  ExpectRelativeNear(axial_recovery.equilibrium_end_actions[0U][0U],
                     -axial_force, kAnalyticalTolerance);
  ExpectRelativeNear(axial_recovery.equilibrium_end_actions[1U][0U],
                     axial_force, kAnalyticalTolerance);
  ExpectRelativeNear(axial_recovery.endpoint_section_resultants[0U][0U],
                     axial_force, kAnalyticalTolerance);
  ExpectRelativeNear(axial_recovery.endpoint_section_resultants[1U][0U],
                     axial_force, kAnalyticalTolerance);

  const double torque = -870.0;
  const Vector torsion =
      SolveFixedFirstNode(stiffness, {0.0, 0.0, 0.0, torque, 0.0, 0.0});
  ExpectRelativeNear(
      torsion[9U],
      torque * length / (shear_modulus * section.torsional_constant),
      kAnalyticalTolerance);
  const BeamRecovery torsion_recovery = beam.RecoverBeam(torsion);
  ExpectRelativeNear(torsion_recovery.equilibrium_end_actions[0U][3U], -torque,
                     kAnalyticalTolerance);
  ExpectRelativeNear(torsion_recovery.equilibrium_end_actions[1U][3U], torque,
                     kAnalyticalTolerance);
  ExpectRelativeNear(torsion_recovery.endpoint_section_resultants[0U][1U],
                     torque, kAnalyticalTolerance);

  const double local_yforce = 640.0;
  const Vector local_y =
      SolveFixedFirstNode(stiffness, {0.0, local_yforce, 0.0, 0.0, 0.0, 0.0});
  ExpectRelativeNear(local_y[7U],
                     local_yforce * length * length * length /
                         (3.0 * material.youngs_modulus * section.i22),
                     kAnalyticalTolerance);
  ExpectRelativeNear(local_y[11U],
                     local_yforce * length * length /
                         (2.0 * material.youngs_modulus * section.i22),
                     kAnalyticalTolerance);
  const BeamRecovery local_yrecovery = beam.RecoverBeam(local_y);
  ExpectRelativeNear(local_yrecovery.equilibrium_end_actions[0U][1U],
                     -local_yforce, kAnalyticalTolerance);
  ExpectRelativeNear(local_yrecovery.equilibrium_end_actions[1U][1U],
                     local_yforce, kAnalyticalTolerance);
  ExpectRelativeNear(local_yrecovery.equilibrium_end_actions[0U][5U],
                     -local_yforce * length, kAnalyticalTolerance);
  ExpectRelativeNear(local_yrecovery.endpoint_section_resultants[0U][3U],
                     local_yforce * length, kAnalyticalTolerance);
  EXPECT_NEAR(local_yrecovery.endpoint_section_resultants[1U][3U], 0.0, 1.0e-8);

  const double local_zforce = -510.0;
  const Vector local_z =
      SolveFixedFirstNode(stiffness, {0.0, 0.0, local_zforce, 0.0, 0.0, 0.0});
  ExpectRelativeNear(local_z[8U],
                     local_zforce * length * length * length /
                         (3.0 * material.youngs_modulus * section.i11),
                     kAnalyticalTolerance);
  ExpectRelativeNear(local_z[10U],
                     -local_zforce * length * length /
                         (2.0 * material.youngs_modulus * section.i11),
                     kAnalyticalTolerance);
  const BeamRecovery local_zrecovery = beam.RecoverBeam(local_z);
  ExpectRelativeNear(local_zrecovery.equilibrium_end_actions[0U][2U],
                     -local_zforce, kAnalyticalTolerance);
  ExpectRelativeNear(local_zrecovery.equilibrium_end_actions[1U][2U],
                     local_zforce, kAnalyticalTolerance);
  ExpectRelativeNear(local_zrecovery.equilibrium_end_actions[0U][4U],
                     local_zforce * length, kAnalyticalTolerance);
  ExpectRelativeNear(local_zrecovery.endpoint_section_resultants[0U][2U],
                     -local_zforce * length, kAnalyticalTolerance);
  EXPECT_NEAR(local_zrecovery.endpoint_section_resultants[1U][2U], 0.0, 1.0e-8);
}

TEST(EulerBeam3D, RejectsInvalidGeometryAndProperties) {
  const Node origin = MakeNode({0.0, 0.0, 0.0}, 1U);
  const Node unit_x = MakeNode({1.0, 0.0, 0.0}, 2U);
  const auto valid_section = MakeSection();
  const auto valid_material = MakeMaterial();

  const auto expect_failure = [](const Result<EulerBeam3D>& result,
                                 const std::string& code) {
    if (result.HasValue()) {
      const Matrix stiffness = result.Value().LocalStiffness();
      ADD_FAILURE() << "Invalid fixture was accepted; local stiffness finite="
                    << MatrixIsFinite(stiffness) << ", maximum absolute entry="
                    << MaximumAbsoluteEntry(stiffness);
      return;
    }
    EXPECT_EQ(result.GetStatus().Category(), FailureCategory::kModel);
    ASSERT_EQ(result.GetStatus().Diagnostics().size(), 1U);
    EXPECT_EQ(result.GetStatus().Diagnostics()[0U].code, code);
  };

  expect_failure(
      EulerBeam3D::Create(origin, origin, valid_section, valid_material),
      "invalid-beam-length");
  expect_failure(EulerBeam3D::Create(origin, MakeNode({1.0e-12, 0.0, 0.0}, 2U),
                                     valid_section, valid_material),
                 "invalid-beam-length");

  const double coordinate = 1048576.0;
  const Node scaled_first = MakeNode({coordinate, 0.0, 0.0}, 1U);
  const Node below_threshold =
      MakeNode({coordinate + coordinate * 0.5e-12, 0.0, 0.0}, 2U);
  const Node above_threshold =
      MakeNode({coordinate + coordinate * 2.0e-12, 0.0, 0.0}, 2U);
  expect_failure(EulerBeam3D::Create(scaled_first, below_threshold,
                                     valid_section, valid_material),
                 "invalid-beam-length");
  EXPECT_TRUE(EulerBeam3D::Create(scaled_first, above_threshold, valid_section,
                                  valid_material)
                  .HasValue());

  auto parallel_guide = valid_section;
  parallel_guide.first_axis = {1.0, 0.0, 0.0};
  expect_failure(
      EulerBeam3D::Create(origin, unit_x, parallel_guide, valid_material),
      "invalid-beam-guide-vector");
  auto guide_at_threshold = valid_section;
  guide_at_threshold.first_axis = {1.0, 1.0e-12, 0.0};
  expect_failure(
      EulerBeam3D::Create(origin, unit_x, guide_at_threshold, valid_material),
      "invalid-beam-guide-vector");
  auto guide_above_threshold = valid_section;
  guide_above_threshold.first_axis = {1.0, 2.0e-12, 0.0};
  EXPECT_TRUE(
      EulerBeam3D::Create(origin, unit_x, guide_above_threshold, valid_material)
          .HasValue());

  auto invalid_material = valid_material;
  invalid_material.youngs_modulus = 0.0;
  expect_failure(
      EulerBeam3D::Create(origin, unit_x, valid_section, invalid_material),
      "invalid-beam-property");
  invalid_material = valid_material;
  invalid_material.poisson_ratio = -2.0;
  expect_failure(
      EulerBeam3D::Create(origin, unit_x, valid_section, invalid_material),
      "invalid-beam-property");

  for (std::size_t property = 0; property < 4U; ++property) {
    auto invalid_section = valid_section;
    double* properties[] = {&invalid_section.area, &invalid_section.i11,
                            &invalid_section.i22,
                            &invalid_section.torsional_constant};
    *properties[property] = 0.0;
    expect_failure(
        EulerBeam3D::Create(origin, unit_x, invalid_section, valid_material),
        "invalid-beam-property");
  }

  auto overflow_material =
      MakeMaterial(std::numeric_limits<double>::max() / 4.0, 0.25);
  auto overflow_section = valid_section;
  overflow_section.area = 8.0;
  expect_failure(
      EulerBeam3D::Create(origin, unit_x, overflow_section, overflow_material),
      "invalid-beam-property");

  auto underflow_section = valid_section;
  underflow_section.area = std::numeric_limits<double>::denorm_min();
  underflow_section.i11 = std::numeric_limits<double>::denorm_min();
  underflow_section.i22 = std::numeric_limits<double>::denorm_min();
  underflow_section.torsional_constant =
      std::numeric_limits<double>::denorm_min();
  expect_failure(EulerBeam3D::Create(origin, unit_x, underflow_section,
                                     MakeMaterial(0.5, 0.25)),
                 "invalid-beam-property");

  auto length_scaled_section = valid_section;
  length_scaled_section.area = 1.0;
  length_scaled_section.i11 = 1.0;
  length_scaled_section.i22 = 1.0;
  length_scaled_section.torsional_constant = 1.0;
  expect_failure(
      EulerBeam3D::Create(origin, MakeNode({1.0e103, 0.0, 0.0}, 2U),
                          length_scaled_section, MakeMaterial(1.0, 0.25)),
      "invalid-beam-property");

  auto coupled_section = valid_section;
  coupled_section.i12 = 1.0e-9;
  expect_failure(
      EulerBeam3D::Create(origin, unit_x, coupled_section, valid_material),
      "unsupported-coupled-section");
}

TEST(EulerBeam3D, RecoversSectionPointAndDefaultCentroidS11) {
  const double length = 2.0;
  const double epsilon = 0.01;
  const double kappa_y = 0.02;
  const double kappa_z = -0.03;
  const auto material = MakeMaterial();
  auto section = MakeSection({0.0, 1.0, 0.0}, {{0.25, -0.5}, {-0.4, 0.3}});
  const auto beam = AlignedBeam(length, section, material);

  Vector displacement{kElementDofCount};
  displacement[6U] = epsilon * length;
  displacement[7U] = 0.5 * kappa_z * length * length;
  displacement[8U] = -0.5 * kappa_y * length * length;
  displacement[10U] = kappa_y * length;
  displacement[11U] = kappa_z * length;

  const BeamRecovery recovery = beam.RecoverBeam(displacement);
  ASSERT_EQ(recovery.stress_points.size(), 4U);
  for (std::size_t gauss_point = 0; gauss_point < 2U; ++gauss_point) {
    for (std::size_t point = 0; point < section.section_points.size();
         ++point) {
      const BeamStressPoint& stress =
          recovery.stress_points[gauss_point * section.section_points.size() +
                                 point];
      const double x1 = section.section_points[point][0U];
      const double x2 = section.section_points[point][1U];
      EXPECT_EQ(stress.gauss_point, static_cast<int>(gauss_point + 1U));
      EXPECT_EQ(stress.section_point, point + 1U);
      EXPECT_DOUBLE_EQ(stress.x1, x1);
      EXPECT_DOUBLE_EQ(stress.x2, x2);
      ExpectScaledNear(
          stress.s11,
          material.youngs_modulus * (epsilon + x2 * kappa_y - x1 * kappa_z),
          kMatrixTolerance);
    }
  }

  const auto default_beam = AlignedBeam(length, MakeSection(), material);
  const BeamRecovery default_recovery = default_beam.RecoverBeam(displacement);
  ASSERT_EQ(default_recovery.stress_points.size(), 2U);
  for (std::size_t gauss_point = 0; gauss_point < 2U; ++gauss_point) {
    const BeamStressPoint& stress = default_recovery.stress_points[gauss_point];
    EXPECT_EQ(stress.gauss_point, static_cast<int>(gauss_point + 1U));
    EXPECT_EQ(stress.section_point, 0U);
    EXPECT_DOUBLE_EQ(stress.x1, 0.0);
    EXPECT_DOUBLE_EQ(stress.x2, 0.0);
    EXPECT_EQ(stress.source, "fesa-default");
    ExpectScaledNear(stress.s11, material.youngs_modulus * epsilon,
                     kMatrixTolerance);
  }
}

TEST(EulerBeam3D, ReproducesConstantStrainTwistAndCurvaturePatches) {
  const double length = 2.8;
  const double epsilon = -0.014;
  const double twist = 0.023;
  const double kappa_y = -0.031;
  const double kappa_z = 0.047;
  const auto section = MakeSection();
  const auto material = MakeMaterial();
  const double shear_modulus =
      material.youngs_modulus / (2.0 * (1.0 + material.poisson_ratio));
  const auto beam = AlignedBeam(length, section, material);

  Vector displacement{kElementDofCount};
  displacement[6U] = epsilon * length;
  displacement[7U] = 0.5 * kappa_z * length * length;
  displacement[8U] = -0.5 * kappa_y * length * length;
  displacement[9U] = twist * length;
  displacement[10U] = kappa_y * length;
  displacement[11U] = kappa_z * length;

  const std::array<double, 4> expected_strain = {epsilon, twist, kappa_y,
                                                 kappa_z};
  const std::array<double, 4> expected_resultant = {
      material.youngs_modulus * section.area * epsilon,
      shear_modulus * section.torsional_constant * twist,
      material.youngs_modulus * section.i11 * kappa_y,
      material.youngs_modulus * section.i22 * kappa_z};
  const BeamRecovery recovery = beam.RecoverBeam(displacement);
  for (std::size_t point = 0; point < 2U; ++point) {
    for (std::size_t component = 0; component < 4U; ++component) {
      ExpectScaledNear(recovery.gauss_generalized_strains[point][component],
                       expected_strain[component], kMatrixTolerance);
      ExpectScaledNear(recovery.gauss_generalized_resultants[point][component],
                       expected_resultant[component], kMatrixTolerance);
      ExpectScaledNear(recovery.endpoint_section_resultants[point][component],
                       expected_resultant[component], kMatrixTolerance);
    }
  }

  const std::array<std::size_t, 4> end_action_components = {0U, 3U, 4U, 5U};
  for (std::size_t component = 0; component < expected_resultant.size();
       ++component) {
    ExpectScaledNear(
        recovery.equilibrium_end_actions[0U][end_action_components[component]],
        -expected_resultant[component], kMatrixTolerance);
    ExpectScaledNear(
        recovery.equilibrium_end_actions[1U][end_action_components[component]],
        expected_resultant[component], kMatrixTolerance);
  }
  EXPECT_NEAR(recovery.equilibrium_end_actions[0U][1U], 0.0, 1.0e-8);
  EXPECT_NEAR(recovery.equilibrium_end_actions[0U][2U], 0.0, 1.0e-8);
  EXPECT_NEAR(recovery.equilibrium_end_actions[1U][1U], 0.0, 1.0e-8);
  EXPECT_NEAR(recovery.equilibrium_end_actions[1U][2U], 0.0, 1.0e-8);
}

TEST(EulerBeam3D, OnePointNegativeControlHasRankFour) {
  const double length = 3.7;
  auto section = MakeSection();
  section.area = 1.0;
  section.i11 = 0.7;
  section.i22 = 1.2;
  section.torsional_constant = 0.9;
  const auto material = MakeMaterial(4.0, 0.25);

  const Matrix one_point = TestOnlyOnePointStiffness(length, section, material);
  const Matrix production =
      AlignedBeam(length, section, material).LocalStiffness();
  EXPECT_EQ(SymmetricRank(one_point, kRigidTolerance), 4U);
  EXPECT_EQ(SymmetricRank(production, kRigidTolerance), 6U);
}

}  // namespace
}  // namespace fesa