#include "fesa/elements/mitc4_shell.hpp"

#include <gtest/gtest.h>

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

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 youngsModulus = 120.0,
    double poissonRatio = 0.25) {
    return {
        "Material",
        youngsModulus,
        poissonRatio,
        {"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 squaredNorm = 0.0;
    for (std::size_t row = 0U; row < matrix.rows(); ++row) {
        for (std::size_t column = 0U; column < matrix.columns(); ++column) {
            squaredNorm += matrix(row, column) * matrix(row, column);
        }
    }
    return std::sqrt(squaredNorm);
}

double scaledSymmetryError(
    const fesa::Matrix& matrix,
    std::size_t dofsPerNode,
    double elementLength) {
    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 rowScale = row % dofsPerNode < 3U ? elementLength : 1.0;
        for (std::size_t column = 0U; column < matrix.columns(); ++column) {
            const double columnScale =
                column % dofsPerNode < 3U ? elementLength : 1.0;
            scaled(row, column) =
                rowScale * matrix(row, column) * columnScale;
            difference(row, column) = rowScale *
                (matrix(row, column) - matrix(column, row)) * columnScale;
        }
    }
    return frobeniusNorm(difference) / frobeniusNorm(scaled);
}

fesa::Matrix scaledStiffness(
    const fesa::Matrix& matrix,
    std::size_t dofsPerNode,
    double elementLength) {
    fesa::Matrix scaled{matrix.rows(), matrix.columns()};
    for (std::size_t row = 0U; row < matrix.rows(); ++row) {
        const double rowScale = row % dofsPerNode < 3U ? elementLength : 1.0;
        for (std::size_t column = 0U; column < matrix.columns(); ++column) {
            const double columnScale =
                column % dofsPerNode < 3U ? elementLength : 1.0;
            scaled(row, column) =
                rowScale * matrix(row, column) * columnScale;
        }
    }
    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 matrixScale = 0.0;
    for (std::size_t row = 0U; row < size; ++row) {
        for (std::size_t column = 0U; column < size; ++column) {
            matrixScale = (std::max)(matrixScale, std::abs(matrix(row, column)));
        }
    }
    if (matrixScale != 0.0) {
        const double convergenceTolerance = 1.0e-14 * matrixScale;
        const std::size_t iterationLimit = 100U * size * size;
        for (std::size_t iteration = 0U; iteration < iterationLimit; ++iteration) {
            std::size_t pivotRow = 0U;
            std::size_t pivotColumn = 0U;
            double largestOffDiagonal = 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 > largestOffDiagonal) {
                        largestOffDiagonal = candidate;
                        pivotRow = row;
                        pivotColumn = column;
                    }
                }
            }
            if (largestOffDiagonal <= convergenceTolerance) {
                break;
            }

            const double pivot = matrix(pivotRow, pivotColumn);
            const double tau =
                (matrix(pivotColumn, pivotColumn) - matrix(pivotRow, pivotRow)) /
                (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 rowDiagonal = matrix(pivotRow, pivotRow);
            const double columnDiagonal = matrix(pivotColumn, pivotColumn);
            matrix(pivotRow, pivotRow) = rowDiagonal - tangent * pivot;
            matrix(pivotColumn, pivotColumn) = columnDiagonal + tangent * pivot;
            matrix(pivotRow, pivotColumn) = 0.0;
            matrix(pivotColumn, pivotRow) = 0.0;
            for (std::size_t index = 0U; index < size; ++index) {
                if (index == pivotRow || index == pivotColumn) {
                    continue;
                }
                const double rowValue = matrix(index, pivotRow);
                const double columnValue = matrix(index, pivotColumn);
                const double rotatedRow = cosine * rowValue - sine * columnValue;
                const double rotatedColumn = sine * rowValue + cosine * columnValue;
                matrix(index, pivotRow) = rotatedRow;
                matrix(pivotRow, index) = rotatedRow;
                matrix(index, pivotColumn) = rotatedColumn;
                matrix(pivotColumn, index) = rotatedColumn;
            }
        }
    }

    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& scaledMatrix) {
    const auto eigenvalues = symmetricEigenvalues(scaledMatrix);
    double spectralScale = 0.0;
    for (double eigenvalue : eigenvalues) {
        spectralScale = (std::max)(spectralScale, std::abs(eigenvalue));
    }
    return static_cast<std::size_t>(std::count_if(
        eigenvalues.begin(), eigenvalues.end(), [spectralScale](double eigenvalue) {
            return std::abs(eigenvalue) > 1.0e-9 * spectralScale;
        }));
}

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 nodeIndex = 0U; nodeIndex < 4U; ++nodeIndex) {
        for (std::size_t component = 0U; component < 5U; ++component) {
            result[5U * nodeIndex + component] = values[nodeIndex][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 nodeIndex = 0U; nodeIndex < nodes.size(); ++nodeIndex) {
        const Vector3 rotationalTranslation = cross(
            rotation, nodes[nodeIndex].coordinates);
        const std::size_t offset = 6U * nodeIndex;
        for (std::size_t component = 0U; component < 3U; ++component) {
            mode[offset + component] =
                translation[component] + rotationalTranslation[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> naturalNodes{
        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 < naturalNodes.size(); ++point) {
        const auto shape = fesa::Mitc4Shell::shapeFunctions(
            naturalNodes[point][0], naturalNodes[point][1]);
        for (std::size_t nodeIndex = 0U; nodeIndex < naturalNodes.size(); ++nodeIndex) {
            EXPECT_DOUBLE_EQ(shape.values[nodeIndex], point == nodeIndex ? 1.0 : 0.0);
        }
    }

    const auto shape = fesa::Mitc4Shell::shapeFunctions(0.25, -0.5);
    double valueSum = 0.0;
    double xiDerivativeSum = 0.0;
    double etaDerivativeSum = 0.0;
    for (std::size_t nodeIndex = 0U; nodeIndex < 4U; ++nodeIndex) {
        valueSum += shape.values[nodeIndex];
        xiDerivativeSum += shape.xiDerivatives[nodeIndex];
        etaDerivativeSum += shape.etaDerivatives[nodeIndex];
    }
    EXPECT_DOUBLE_EQ(valueSum, 1.0);
    EXPECT_DOUBLE_EQ(xiDerivativeSum, 0.0);
    EXPECT_DOUBLE_EQ(etaDerivativeSum, 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 nodeIndex = 0U; nodeIndex < 4U; ++nodeIndex) {
        const std::size_t physicalOffset = 5U * nodeIndex;
        const std::size_t globalOffset = 6U * nodeIndex;
        for (std::size_t component = 0U; component < 3U; ++component) {
            EXPECT_DOUBLE_EQ(
                physical(physicalOffset + component, globalOffset + component),
                1.0);
        }
        EXPECT_DOUBLE_EQ(physical(physicalOffset + 3U, globalOffset + 4U), 1.0);
        EXPECT_DOUBLE_EQ(physical(physicalOffset + 4U, globalOffset + 5U), 1.0);
        EXPECT_DOUBLE_EQ(drilling(nodeIndex, globalOffset + 3U), 1.0);

        for (std::size_t globalDof = 0U; globalDof < 24U; ++globalDof) {
            if (globalDof != globalOffset + 4U) {
                EXPECT_DOUBLE_EQ(physical(physicalOffset + 3U, globalDof), 0.0);
            }
            if (globalDof != globalOffset + 5U) {
                EXPECT_DOUBLE_EQ(physical(physicalOffset + 4U, globalDof), 0.0);
            }
            if (globalDof != globalOffset + 3U) {
                EXPECT_DOUBLE_EQ(drilling(nodeIndex, globalDof), 0.0);
            }
        }
    }

    auto invalidDirectors = directors({1.0, 0.0, 0.0});
    invalidDirectors[2] = {0.0, 0.0, 0.0};
    EXPECT_FALSE(fesa::Mitc4Shell::create(
                     nodePointers(nodes), invalidDirectors, 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.xiZeta, (std::array<double, 2>{0.75, 0.25}));
    EXPECT_EQ(weights.etaZeta, (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 forceScale = 7.0;
    constexpr double lengthScale = 3.0;
    const auto scaledCandidate = fesa::Mitc4Shell::create(
        nodePointers(nodes),
        directors(),
        section(2.0 * lengthScale),
        material(120.0 * forceScale / (lengthScale * lengthScale), 0.25));
    ASSERT_TRUE(scaledCandidate.hasValue());
    const auto& scaled = scaledCandidate.value();
    fesa::Matrix expectedCps{3U, 3U};
    fesa::Matrix expectedC5{5U, 5U};
    fesa::Matrix expectedA{3U, 3U};
    fesa::Matrix expectedD{3U, 3U};
    fesa::Matrix expectedAs{2U, 2U};
    for (std::size_t row = 0U; row < 3U; ++row) {
        for (std::size_t column = 0U; column < 3U; ++column) {
            expectedCps(row, column) =
                cps(row, column) * forceScale / (lengthScale * lengthScale);
            expectedA(row, column) = a(row, column) * forceScale / lengthScale;
            expectedD(row, column) = d(row, column) * forceScale * lengthScale;
        }
    }
    for (std::size_t row = 0U; row < 2U; ++row) {
        for (std::size_t column = 0U; column < 2U; ++column) {
            expectedAs(row, column) = as(row, column) * forceScale / lengthScale;
        }
    }
    for (std::size_t row = 0U; row < 5U; ++row) {
        for (std::size_t column = 0U; column < 5U; ++column) {
            expectedC5(row, column) =
                c5(row, column) * forceScale / (lengthScale * lengthScale);
        }
    }
    expectMatrixNear(scaled.planeStressConstitutive(), expectedCps);
    expectMatrixNear(scaled.materialConstitutive5(), expectedC5);
    expectMatrixNear(scaled.membraneSectionMatrix(), expectedA);
    expectMatrixNear(scaled.bendingSectionMatrix(), expectedD);
    expectMatrixNear(scaled.transverseShearSectionMatrix(), expectedAs);

    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].naturalCoordinates, expected[point]);
        EXPECT_DOUBLE_EQ(points[point].weight, 1.0);
    }
}

// MITC4-KERNEL-001
TEST(Mitc4ShellKernel, FormsFiniteScaledSymmetricPhysicalAndStabilizedStiffness) {
    const auto nodes = planarNodes();
    const auto shellCandidate = fesa::Mitc4Shell::create(
        nodePointers(nodes), directors(), section(), material());
    ASSERT_TRUE(shellCandidate.hasValue());

    const auto stiffnessCandidate = shellCandidate.value().stiffness();
    ASSERT_TRUE(stiffnessCandidate.hasValue());
    const auto& stiffness = stiffnessCandidate.value();
    EXPECT_EQ(stiffness.physicalLocal20.rows(), 20U);
    EXPECT_EQ(stiffness.physicalLocal20.columns(), 20U);
    EXPECT_EQ(stiffness.physicalGlobal24.rows(), 24U);
    EXPECT_EQ(stiffness.drillingGlobal24.rows(), 24U);
    EXPECT_EQ(stiffness.stabilizedGlobal24.rows(), 24U);
    for (const fesa::Matrix* matrix : {
             &stiffness.physicalLocal20,
             &stiffness.physicalGlobal24,
             &stiffness.drillingGlobal24,
             &stiffness.stabilizedGlobal24}) {
        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.physicalLocal20, 5U, 2.0), 1.0e-12);
    EXPECT_LE(scaledSymmetryError(stiffness.physicalGlobal24, 6U, 2.0), 1.0e-12);
    EXPECT_LE(scaledSymmetryError(stiffness.drillingGlobal24, 6U, 2.0), 1.0e-12);
    EXPECT_LE(scaledSymmetryError(stiffness.stabilizedGlobal24, 6U, 2.0), 1.0e-12);

    const auto repeatedCandidate = shellCandidate.value().stiffness();
    ASSERT_TRUE(repeatedCandidate.hasValue());
    const auto& repeated = repeatedCandidate.value();
    expectMatrixNear(repeated.physicalLocal20, stiffness.physicalLocal20, 0.0);
    expectMatrixNear(repeated.physicalGlobal24, stiffness.physicalGlobal24, 0.0);
    expectMatrixNear(repeated.drillingGlobal24, stiffness.drillingGlobal24, 0.0);
    expectMatrixNear(repeated.stabilizedGlobal24, stiffness.stabilizedGlobal24, 0.0);
    EXPECT_DOUBLE_EQ(repeated.drillingStiffness, stiffness.drillingStiffness);
}

// MITC4-KERNEL-002
TEST(Mitc4ShellKernel, PreservesPhysicalEnergyUnderTwentyToTwentyFourCongruence) {
    const auto nodes = planarNodes();
    const auto shellCandidate = fesa::Mitc4Shell::create(
        nodePointers(nodes), directors(), section(), material());
    ASSERT_TRUE(shellCandidate.hasValue());
    const auto stiffnessCandidate = shellCandidate.value().stiffness();
    ASSERT_TRUE(stiffnessCandidate.hasValue());
    const auto& stiffness = stiffnessCandidate.value();

    fesa::Vector globalField{24U};
    for (std::size_t index = 0U; index < globalField.size(); ++index) {
        globalField[index] = 0.125 * static_cast<double>(
            static_cast<int>(index % 7U) - 3);
    }
    const auto physicalField20 =
        shellCandidate.value().physicalTransformation20().multiply(globalField);
    const double localEnergy = quadraticEnergy(
        stiffness.physicalLocal20, physicalField20);
    const double globalEnergy = quadraticEnergy(
        stiffness.physicalGlobal24, globalField);
    ASSERT_NE(localEnergy, 0.0);
    ASSERT_NE(globalEnergy, 0.0);
    EXPECT_LE(
        std::abs(globalEnergy - localEnergy) /
            (std::abs(globalEnergy) + std::abs(localEnergy)),
        1.0e-12);
}

// MITC4-KERNEL-003
TEST(Mitc4ShellKernel, RetainsSixRigidModesAndHasExpectedPhysicalAndStabilizedRank) {
    const auto nodes = planarNodes();
    const auto shellCandidate = fesa::Mitc4Shell::create(
        nodePointers(nodes), directors(), section(), material());
    ASSERT_TRUE(shellCandidate.hasValue());
    const auto stiffnessCandidate = shellCandidate.value().stiffness();
    ASSERT_TRUE(stiffnessCandidate.hasValue());
    const auto& stiffness = stiffnessCandidate.value();

    const auto scaledPhysical20 =
        scaledStiffness(stiffness.physicalLocal20, 5U, 2.0);
    const auto scaledPhysical24 =
        scaledStiffness(stiffness.physicalGlobal24, 6U, 2.0);
    const auto scaledStabilized24 =
        scaledStiffness(stiffness.stabilizedGlobal24, 6U, 2.0);
    EXPECT_EQ(numericalRank(scaledPhysical20), 14U);
    EXPECT_EQ(numericalRank(scaledStabilized24), 18U);

    const std::array<fesa::Vector, 6> rigidModes{
        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 physicalNorm = symmetricOperatorNorm(scaledPhysical24);
    const double stabilizedNorm = symmetricOperatorNorm(scaledStabilized24);
    ASSERT_GT(physicalNorm, 0.0);
    ASSERT_GT(stabilizedNorm, 0.0);
    for (const auto& rigidMode : rigidModes) {
        fesa::Vector scaledMode = rigidMode;
        for (std::size_t nodeIndex = 0U; nodeIndex < 4U; ++nodeIndex) {
            for (std::size_t component = 0U; component < 3U; ++component) {
                scaledMode[6U * nodeIndex + component] /= 2.0;
            }
        }
        const double modeNorm = scaledMode.norm();
        ASSERT_GT(modeNorm, 0.0);
        EXPECT_LE(
            scaledPhysical24.multiply(scaledMode).norm() /
                (physicalNorm * modeNorm),
            1.0e-10);
        EXPECT_LE(
            scaledStabilized24.multiply(scaledMode).norm() /
                (stabilizedNorm * modeNorm),
            1.0e-10);
    }
}

// MITC4-KERNEL-004
TEST(Mitc4ShellPatch, ReproducesIndependentMembraneBendingShearAndTwistFields) {
    const auto nodes = planarNodes();
    const auto shellCandidate = fesa::Mitc4Shell::create(
        nodePointers(nodes), directors(), section(), material());
    ASSERT_TRUE(shellCandidate.hasValue());
    const auto& shell = shellCandidate.value();
    const auto stiffnessCandidate = shell.stiffness();
    ASSERT_TRUE(stiffnessCandidate.hasValue());
    const auto& stiffness = stiffnessCandidate.value().physicalLocal20;
    constexpr double magnitude = 0.2;
    const double gauss = 1.0 / std::sqrt(3.0);

    std::array<std::array<double, 5>, 4> e11Values{};
    std::array<std::array<double, 5>, 4> e22Values{};
    std::array<std::array<double, 5>, 4> g12Values{};
    std::array<std::array<double, 5>, 4> k11Values{};
    std::array<std::array<double, 5>, 4> k22Values{};
    std::array<std::array<double, 5>, 4> g13Values{};
    std::array<std::array<double, 5>, 4> g23Values{};
    std::array<std::array<double, 5>, 4> k12Values{};
    for (std::size_t nodeIndex = 0U; nodeIndex < nodes.size(); ++nodeIndex) {
        const double x = nodes[nodeIndex].coordinates[0];
        const double y = nodes[nodeIndex].coordinates[1];
        e11Values[nodeIndex][0] = magnitude * x;
        e22Values[nodeIndex][1] = magnitude * y;
        g12Values[nodeIndex][0] = 0.5 * magnitude * y;
        g12Values[nodeIndex][1] = 0.5 * magnitude * x;
        k11Values[nodeIndex][4] = magnitude * x;
        k22Values[nodeIndex][3] = -magnitude * y;
        g13Values[nodeIndex][2] = magnitude * x;
        g23Values[nodeIndex][2] = magnitude * y;
        k12Values[nodeIndex][2] = -0.5 * magnitude * x * y;
        k12Values[nodeIndex][3] = -0.5 * magnitude * x;
        k12Values[nodeIndex][4] = 0.5 * magnitude * y;
    }

    const std::array<fesa::Vector, 8> fields{
        physicalField(e11Values), physicalField(e22Values),
        physicalField(g12Values), physicalField(k11Values),
        physicalField(k22Values), physicalField(g13Values),
        physicalField(g23Values), physicalField(k12Values)};
    expectStrain(shell, fields[0], gauss, -gauss, gauss, {magnitude, 0.0, 0.0, 0.0, 0.0});
    expectStrain(shell, fields[1], gauss, -gauss, gauss, {0.0, magnitude, 0.0, 0.0, 0.0});
    expectStrain(shell, fields[2], gauss, -gauss, gauss, {0.0, 0.0, magnitude, 0.0, 0.0});
    expectStrain(shell, fields[3], gauss, -gauss, gauss, {gauss * magnitude, 0.0, 0.0, 0.0, 0.0});
    expectStrain(shell, fields[4], gauss, -gauss, gauss, {0.0, gauss * magnitude, 0.0, 0.0, 0.0});
    expectStrain(shell, fields[5], gauss, -gauss, gauss, {0.0, 0.0, 0.0, magnitude, 0.0});
    expectStrain(shell, fields[6], gauss, -gauss, gauss, {0.0, 0.0, 0.0, 0.0, magnitude});
    expectStrain(shell, fields[7], gauss, -gauss, gauss, {0.0, 0.0, gauss * magnitude, 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 shellCandidate = fesa::Mitc4Shell::create(
        nodePointers(nodes), directors(), section(0.1), material());
    ASSERT_TRUE(shellCandidate.hasValue());
    const auto stiffnessCandidate = shellCandidate.value().stiffness();
    ASSERT_TRUE(stiffnessCandidate.hasValue());
    const auto& stiffness = stiffnessCandidate.value();

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

// MITC4-KERNEL-006
TEST(Mitc4ShellDrilling, FailsNonfiniteReferenceAndStabilizesEachPureDrillCoordinate) {
    const auto nodes = planarNodes();
    const auto shellCandidate = fesa::Mitc4Shell::create(
        nodePointers(nodes), directors(), section(), material());
    ASSERT_TRUE(shellCandidate.hasValue());
    const auto stiffnessCandidate = shellCandidate.value().stiffness();
    ASSERT_TRUE(stiffnessCandidate.hasValue());
    const auto& stiffness = stiffnessCandidate.value();

    for (std::size_t nodeIndex = 0U; nodeIndex < 4U; ++nodeIndex) {
        fesa::Vector pureDrill{24U};
        pureDrill[6U * nodeIndex + 5U] = 1.0;
        EXPECT_DOUBLE_EQ(
            stiffness.physicalGlobal24.multiply(pureDrill).norm(), 0.0);
        const auto drillAction = stiffness.drillingGlobal24.multiply(pureDrill);
        EXPECT_DOUBLE_EQ(drillAction[6U * nodeIndex + 5U], stiffness.drillingStiffness);
        EXPECT_GT(quadraticEnergy(stiffness.drillingGlobal24, pureDrill), 0.0);
    }

    const std::array<fesa::Node, 4> extremeNodes{
        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 extremeShell = fesa::Mitc4Shell::create(
        nodePointers(extremeNodes), directors(), section(1.0), material(1.0e300));
    ASSERT_TRUE(extremeShell.hasValue());
    const auto failure = extremeShell.value().stiffness();
    ASSERT_FALSE(failure.hasValue());
    ASSERT_EQ(failure.status().diagnostics().size(), 1U);
    EXPECT_EQ(failure.status().diagnostics()[0].code, "invalid-shell-stiffness");
    const auto repeatedFailure = extremeShell.value().stiffness();
    ASSERT_FALSE(repeatedFailure.hasValue());
    ASSERT_EQ(repeatedFailure.status().diagnostics().size(), 1U);
    EXPECT_EQ(
        repeatedFailure.status().diagnostics()[0].code,
        failure.status().diagnostics()[0].code);
    EXPECT_EQ(
        repeatedFailure.status().diagnostics()[0].message,
        failure.status().diagnostics()[0].message);
}

// MITC4-KERNEL-007
TEST(Mitc4ShellDrilling, ExcludesPureDrillFromPhysicalRecoveryAndEnergy) {
    const auto nodes = planarNodes();
    const auto shellCandidate = fesa::Mitc4Shell::create(
        nodePointers(nodes), directors(), section(), material());
    ASSERT_TRUE(shellCandidate.hasValue());
    const auto& shell = shellCandidate.value();
    const auto stiffnessCandidate = shell.stiffness();
    ASSERT_TRUE(stiffnessCandidate.hasValue());

    for (std::size_t nodeIndex = 0U; nodeIndex < nodes.size(); ++nodeIndex) {
        fesa::Vector pureDrill{24U};
        pureDrill[6U * nodeIndex + 5U] = 1.0;
        EXPECT_GT(
            stiffnessCandidate.value().stabilizedGlobal24.multiply(pureDrill).norm(),
            0.0);

        const auto recoveryCandidate = shell.recoverPhysical(pureDrill);
        ASSERT_TRUE(recoveryCandidate.hasValue());
        const auto& recovery = recoveryCandidate.value();
        EXPECT_DOUBLE_EQ(recovery.strainEnergy, 0.0);
        for (const auto& point : recovery.points) {
            for (double value : point.generalizedStrain) {
                EXPECT_DOUBLE_EQ(value, 0.0);
            }
            for (double value : point.sectionResultant) {
                EXPECT_DOUBLE_EQ(value, 0.0);
            }
            for (const auto& stress : point.inPlaneStress) {
                for (double value : stress) {
                    EXPECT_DOUBLE_EQ(value, 0.0);
                }
            }
        }
    }
}

// MITC4-PHYSREC-001
TEST(Mitc4ShellPhysicalRecovery, RecoversHandFieldAtFixedLocationsAndSectionPositions) {
    const auto nodes = planarNodes();
    const auto shellCandidate = fesa::Mitc4Shell::create(
        nodePointers(nodes), directors(), section(), material());
    ASSERT_TRUE(shellCandidate.hasValue());
    const auto& shell = shellCandidate.value();

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

    const auto recoveryCandidate = shell.recoverPhysical(globalField);
    ASSERT_TRUE(recoveryCandidate.hasValue());
    const auto& recovery = recoveryCandidate.value();
    const double gauss = 1.0 / std::sqrt(3.0);
    const std::array<std::array<double, 2>, 4> expectedCoordinates{
        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> expectedResultant{
        22.4, -6.4, 19.2, 22.4, -6.4, 8.0, 32.0, -24.0};
    constexpr std::array<std::array<double, 3>, 3> expectedStress{
        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(), expectedCoordinates.size());
    for (std::size_t pointIndex = 0U;
         pointIndex < recovery.points.size(); ++pointIndex) {
        const auto& point = recovery.points[pointIndex];
        EXPECT_EQ(point.naturalCoordinates, expectedCoordinates[pointIndex]);
        expectOrthonormalRightHanded(point.localFrame);
        expectVectorNear(point.localFrame.e1, {1.0, 0.0, 0.0});
        expectVectorNear(point.localFrame.e2, {0.0, 1.0, 0.0});
        expectVectorNear(point.localFrame.e3, {0.0, 0.0, 1.0});
        for (std::size_t component = 0U;
             component < generalized.size(); ++component) {
            EXPECT_NEAR(
                point.generalizedStrain[component], generalized[component],
                1.0e-12);
            EXPECT_NEAR(
                point.sectionResultant[component], expectedResultant[component],
                1.0e-12);
        }
        for (std::size_t position = 0U;
             position < expectedStress.size(); ++position) {
            for (std::size_t component = 0U;
                 component < expectedStress[position].size(); ++component) {
                EXPECT_NEAR(
                    point.inPlaneStress[position][component],
                    expectedStress[position][component], 1.0e-12);
            }
        }
    }
    EXPECT_NEAR(recovery.strainEnergy, 72.16, 1.0e-12);
}