#include "fesa/elements/mitc4_shell.hpp"

#include <algorithm>
#include <cmath>
#include <limits>
#include <stdexcept>
#include <string>
#include <utility>

namespace fesa {
namespace {

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

constexpr std::size_t kNodeCount = 4U;
constexpr std::size_t kPhysicalDofsPerNode = 5U;
constexpr std::size_t kGlobalDofsPerNode = 6U;
constexpr std::size_t kPhysicalDofCount = 20U;
constexpr std::size_t kGlobalDofCount = 24U;
constexpr double kShearCorrection = 5.0 / 6.0;
constexpr double kFrameTolerance = 1.0e-12;

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

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

Vector3 scale(double factor, const Vector3& value) {
    return {factor * value[0], factor * value[1], factor * value[2]};
}

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::hypot(value[0], value[1], value[2]);
}

bool isFinite(const Vector3& value) {
    return std::all_of(value.begin(), value.end(), [](double component) {
        return std::isfinite(component);
    });
}

Vector3 normalized(const Vector3& value) {
    return scale(1.0 / norm(value), value);
}

Vector3 weightedSum(
    const std::array<double, kNodeCount>& weights,
    const std::array<Vector3, kNodeCount>& values) {
    Vector3 result{};
    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        result = add(result, scale(weights[node], values[node]));
    }
    return result;
}

Vector3 derivativeSum(
    const std::array<double, kNodeCount>& derivatives,
    const std::array<Vector3, kNodeCount>& values) {
    std::array<Vector3, kNodeCount> relative{};
    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        relative[node] = subtract(values[node], values[0]);
    }
    return weightedSum(derivatives, relative);
}

bool sameCoordinates(const Vector3& left, const Vector3& right) {
    return left == right;
}

std::array<Vector3, kNodeCount> nodalTangentsA(
    const std::array<Vector3, kNodeCount>& directors) {
    const std::array<Vector3, 3> globalAxes{
        Vector3{1.0, 0.0, 0.0},
        Vector3{0.0, 1.0, 0.0},
        Vector3{0.0, 0.0, 1.0}};
    std::array<Vector3, kNodeCount> tangents{};
    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        std::size_t selected = 0U;
        double alignment = std::abs(dot(globalAxes[0], directors[node]));
        for (std::size_t axis = 1U; axis < globalAxes.size(); ++axis) {
            const double candidate = std::abs(dot(globalAxes[axis], directors[node]));
            if (candidate < alignment) {
                selected = axis;
                alignment = candidate;
            }
        }
        tangents[node] = normalized(subtract(
            globalAxes[selected],
            scale(dot(globalAxes[selected], directors[node]), directors[node])));
    }
    return tangents;
}

std::array<Vector3, kNodeCount> nodalTangentsB(
    const std::array<Vector3, kNodeCount>& directors,
    const std::array<Vector3, kNodeCount>& tangentA) {
    std::array<Vector3, kNodeCount> tangents{};
    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        tangents[node] = cross(directors[node], tangentA[node]);
    }
    return tangents;
}

std::string elementIdentity(const std::array<const Node*, kNodeCount>& nodes) {
    std::string identity;
    for (const Node* node : nodes) {
        if (node == nullptr) {
            continue;
        }
        if (!identity.empty()) {
            identity += "-";
        }
        identity += node->sourceId.source_label_text;
    }
    return identity;
}

Result<Mitc4Shell> modelFailure(
    std::string code,
    const SourceLocation& location,
    const std::string& identity,
    std::string message) {
    return Result<Mitc4Shell>::Failure(Status::Failure(
        FailureCategory::kModel,
        {{Severity::kError,
          std::move(code),
          location,
          "*ELEMENT",
          identity,
          std::move(message)}}));
}

std::array<std::array<double, 3>, 3> covariantStrainColumn(
    const std::array<Vector3, 3>& covariant,
    const std::array<Vector3, 3>& derivatives) {
    std::array<std::array<double, 3>, 3> strain{};
    for (std::size_t first = 0U; first < 3U; ++first) {
        for (std::size_t second = 0U; second < 3U; ++second) {
            strain[first][second] = 0.5 * (
                dot(covariant[first], derivatives[second]) +
                dot(covariant[second], derivatives[first]));
        }
    }
    // Thickness stretch is excluded from the five-component shell law.
    strain[2U][2U] = 0.0;
    return strain;
}

std::array<std::array<double, 3>, 3> reconstructCartesianStrain(
    const std::array<std::array<double, 3>, 3>& covariantStrain,
    const std::array<Vector3, 3>& reciprocal) {
    std::array<std::array<double, 3>, 3> tensor{};
    for (std::size_t row = 0U; row < 3U; ++row) {
        for (std::size_t column = 0U; column < 3U; ++column) {
            for (std::size_t first = 0U; first < 3U; ++first) {
                for (std::size_t second = 0U; second < 3U; ++second) {
                    tensor[row][column] +=
                        covariantStrain[first][second] *
                        reciprocal[first][row] * reciprocal[second][column];
                }
            }
        }
    }
    return tensor;
}

double frameComponent(
    const Vector3& left,
    const std::array<std::array<double, 3>, 3>& tensor,
    const Vector3& right) {
    double value = 0.0;
    for (std::size_t row = 0U; row < 3U; ++row) {
        for (std::size_t column = 0U; column < 3U; ++column) {
            value += left[row] * tensor[row][column] * right[column];
        }
    }
    return value;
}

std::array<double, 5> localEngineeringComponents(
    const std::array<std::array<double, 3>, 3>& tensor,
    const Mitc4LocalFrame& frame) {
    return {
        frameComponent(frame.e1, tensor, frame.e1),
        frameComponent(frame.e2, tensor, frame.e2),
        2.0 * frameComponent(frame.e1, tensor, frame.e2),
        2.0 * frameComponent(frame.e1, tensor, frame.e3),
        2.0 * frameComponent(frame.e2, tensor, frame.e3)};
}

Matrix scaledMatrix(const Matrix& source, double factor) {
    Matrix result{source.Rows(), source.Columns()};
    for (std::size_t row = 0U; row < source.Rows(); ++row) {
        for (std::size_t column = 0U; column < source.Columns(); ++column) {
            result(row, column) = factor * source(row, column);
        }
    }
    return result;
}

Matrix congruence(const Matrix& local, const Matrix& transformation) {
    if (local.Rows() != local.Columns() ||
        local.Rows() != transformation.Rows()) {
        throw std::invalid_argument{"MITC4 congruence dimensions are incompatible."};
    }
    Matrix result{transformation.Columns(), transformation.Columns()};
    for (std::size_t row = 0U; row < result.Rows(); ++row) {
        for (std::size_t column = row; column < result.Columns(); ++column) {
            double value = 0.0;
            for (std::size_t localRow = 0U; localRow < local.Rows(); ++localRow) {
                for (std::size_t localColumn = 0U;
                     localColumn < local.Columns(); ++localColumn) {
                    value += transformation(localRow, row) *
                             local(localRow, localColumn) *
                             transformation(localColumn, column);
                }
            }
            result(row, column) = value;
            result(column, row) = value;
        }
    }
    return result;
}

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

bool isFinite(const Vector& vector) {
    for (std::size_t index = 0U; index < vector.Size(); ++index) {
        if (!std::isfinite(vector[index])) {
            return false;
        }
    }
    return true;
}

Result<Mitc4Stiffness> stiffnessFailure(
    const SourceLocation& location,
    const std::string& identity,
    std::string message) {
    return Result<Mitc4Stiffness>::Failure(Status::Failure(
        FailureCategory::kModel,
        {{Severity::kError,
          "invalid-shell-stiffness",
          location,
          "*ELEMENT",
          identity,
          std::move(message)}}));
}

Result<Mitc4PhysicalRecovery> recoveryFailure(
    const SourceLocation& location,
    const std::string& identity,
    std::string message) {
    return Result<Mitc4PhysicalRecovery>::Failure(Status::Failure(
        FailureCategory::kModel,
        {{Severity::kError,
          "invalid-shell-recovery",
          location,
          "*ELEMENT",
          identity,
          std::move(message)}}));
}

} // namespace

Result<Mitc4Shell> Mitc4Shell::create(
    std::array<const Node*, 4> nodes,
    std::array<std::array<double, 3>, 4> initialDirectors,
    const ShellSection& section,
    const LinearElasticMaterial& material) {
    const std::string identity = elementIdentity(nodes);
    if (std::any_of(nodes.begin(), nodes.end(), [](const Node* value) {
            return value == nullptr;
        })) {
        return modelFailure(
            "invalid-shell-geometry", section.location, identity,
            "MITC4 creation requires four valid node references.");
    }

    std::array<Vector3, kNodeCount> coordinates{};
    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        coordinates[node] = nodes[node]->coordinates;
        if (!isFinite(coordinates[node])) {
            return modelFailure(
                "invalid-shell-geometry", nodes[node]->location, identity,
                "MITC4 node coordinates must be finite.");
        }
        for (std::size_t previous = 0U; previous < node; ++previous) {
            if (sameCoordinates(coordinates[node], coordinates[previous])) {
                return modelFailure(
                    "invalid-shell-geometry", nodes[node]->location, identity,
                    "MITC4 nodes must have distinct coordinates.");
            }
        }
    }

    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        const double directorNorm = norm(initialDirectors[node]);
        if (!isFinite(initialDirectors[node]) || !std::isfinite(directorNorm) ||
            std::abs(directorNorm - 1.0) > kFrameTolerance) {
            return modelFailure(
                "invalid-shell-director", nodes[node]->location, identity,
                "MITC4 initial directors must be finite unit vectors.");
        }
    }

    if (!std::isfinite(section.thickness) || !(section.thickness > 0.0)) {
        return modelFailure(
            "invalid-shell-section", section.location, identity,
            "MITC4 shell thickness must be finite and positive.");
    }
    if (!std::isfinite(material.youngsModulus) ||
        !(material.youngsModulus > 0.0) ||
        !std::isfinite(material.poissonRatio) ||
        !(material.poissonRatio > -1.0) || !(material.poissonRatio < 0.5)) {
        return modelFailure(
            "invalid-shell-material", material.location, identity,
            "MITC4 isotropic material requires finite E>0 and -1<nu<0.5.");
    }

    const auto centerShape = shapeFunctions(0.0, 0.0);
    const Vector3 centerXi = derivativeSum(centerShape.xiDerivatives, coordinates);
    const Vector3 centerEta = derivativeSum(centerShape.etaDerivatives, coordinates);
    const Vector3 centerArea = cross(centerXi, centerEta);
    const double centerMeasure = norm(centerArea);
    if (!isFinite(centerArea) || !std::isfinite(centerMeasure) ||
        !(centerMeasure > 0.0)) {
        return modelFailure(
            "invalid-shell-geometry", nodes[0]->location, identity,
            "MITC4 center surface basis must be finite and nonzero.");
    }
    const Vector3 normalCandidate = scale(1.0 / centerMeasure, centerArea);
    if (std::any_of(
            initialDirectors.begin(), initialDirectors.end(),
            [&normalCandidate](const Vector3& director) {
                return !(dot(normalCandidate, director) > 0.0);
            })) {
        return modelFailure(
            "invalid-shell-director", nodes[0]->location, identity,
            "MITC4 directors must follow the source-order positive face.");
    }

    const auto tangentA = nodalTangentsA(initialDirectors);
    const auto tangentB = nodalTangentsB(initialDirectors, tangentA);
    Mitc4Shell shell{
        coordinates,
        initialDirectors,
        tangentA,
        tangentB,
        normalCandidate,
        section.thickness,
        material.youngsModulus,
        material.poissonRatio,
        nodes[0]->location,
        identity};

    const double shearModulus =
        material.youngsModulus / (2.0 * (1.0 + material.poissonRatio));
    const double planeStressFactor = material.youngsModulus /
                                     (1.0 - material.poissonRatio * material.poissonRatio);
    const double thicknessCubed = section.thickness * section.thickness *
                                  section.thickness;
    const std::array<double, 6> derived{
        shearModulus,
        planeStressFactor,
        kShearCorrection * shearModulus,
        planeStressFactor * section.thickness,
        planeStressFactor * thicknessCubed / 12.0,
        kShearCorrection * shearModulus * section.thickness};
    if (std::any_of(derived.begin(), derived.end(), [](double value) {
            return !std::isfinite(value) || !(value > 0.0);
        })) {
        return modelFailure(
            "invalid-shell-material", material.location, identity,
            "Derived MITC4 constitutive coefficients must be finite and positive.");
    }

    GeometryData geometry{};
    if (!shell.evaluateGeometry(0.0, 0.0, 0.0, geometry)) {
        return modelFailure(
            "invalid-shell-jacobian", nodes[0]->location, identity,
            "MITC4 center frame or Jacobian is invalid.");
    }
    for (const auto& point : volumeQuadrature()) {
        if (!shell.evaluateGeometry(
                point.naturalCoordinates[0],
                point.naturalCoordinates[1],
                point.naturalCoordinates[2],
                geometry)) {
            return modelFailure(
                "invalid-shell-jacobian", nodes[0]->location, identity,
                "MITC4 quadrature frame or Jacobian is invalid.");
        }
    }
    constexpr std::array<Vector3, 4> tyingPoints{
        Vector3{0.0, -1.0, 0.0},
        Vector3{0.0, 1.0, 0.0},
        Vector3{-1.0, 0.0, 0.0},
        Vector3{1.0, 0.0, 0.0}};
    for (const auto& point : tyingPoints) {
        if (!shell.evaluateGeometry(point[0], point[1], point[2], geometry)) {
            return modelFailure(
                "invalid-shell-jacobian", nodes[0]->location, identity,
                "MITC4 tying-point frame or Jacobian is invalid.");
        }
    }

    return Result<Mitc4Shell>::Success(std::move(shell));
}

Mitc4ShapeFunctions Mitc4Shell::shapeFunctions(
    double xi,
    double eta) noexcept {
    constexpr std::array<double, kNodeCount> xiSigns{-1.0, 1.0, 1.0, -1.0};
    constexpr std::array<double, kNodeCount> etaSigns{-1.0, -1.0, 1.0, 1.0};
    Mitc4ShapeFunctions shape{};
    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        shape.values[node] = 0.25 * (1.0 + xiSigns[node] * xi) *
                             (1.0 + etaSigns[node] * eta);
        shape.xiDerivatives[node] =
            0.25 * xiSigns[node] * (1.0 + etaSigns[node] * eta);
        shape.etaDerivatives[node] =
            0.25 * etaSigns[node] * (1.0 + xiSigns[node] * xi);
    }
    return shape;
}

Mitc4TyingWeights Mitc4Shell::tyingWeights(double xi, double eta) noexcept {
    return {
        {(1.0 - eta) * 0.5, (1.0 + eta) * 0.5},
        {(1.0 - xi) * 0.5, (1.0 + xi) * 0.5}};
}

const std::array<Mitc4QuadraturePoint, 8>&
Mitc4Shell::volumeQuadrature() noexcept {
    static const std::array<Mitc4QuadraturePoint, 8> points = [] {
        const double gauss = 1.0 / std::sqrt(3.0);
        return std::array<Mitc4QuadraturePoint, 8>{
            Mitc4QuadraturePoint{{-gauss, -gauss, -gauss}, 1.0},
            {{-gauss, -gauss, gauss}, 1.0},
            {{gauss, -gauss, -gauss}, 1.0},
            {{gauss, -gauss, gauss}, 1.0},
            {{gauss, gauss, -gauss}, 1.0},
            {{gauss, gauss, gauss}, 1.0},
            {{-gauss, gauss, -gauss}, 1.0},
            {{-gauss, gauss, gauss}, 1.0}};
    }();
    return points;
}

Mitc4LocalFrame Mitc4Shell::localFrame(double xi, double eta) const {
    GeometryData geometry{};
    if (!evaluateGeometry(xi, eta, 0.0, geometry)) {
        throw std::invalid_argument{"MITC4 local frame is invalid at the requested point."};
    }
    return geometry.frame;
}

Matrix Mitc4Shell::physicalTransformation20() const {
    Matrix transformation{kPhysicalDofCount, kGlobalDofCount};
    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        const std::size_t physicalOffset = node * kPhysicalDofsPerNode;
        const std::size_t globalOffset = node * kGlobalDofsPerNode;
        for (std::size_t component = 0U; component < 3U; ++component) {
            transformation(physicalOffset + component, globalOffset + component) = 1.0;
            transformation(physicalOffset + 3U, globalOffset + 3U + component) =
                tangentA_[node][component];
            transformation(physicalOffset + 4U, globalOffset + 3U + component) =
                tangentB_[node][component];
        }
    }
    return transformation;
}

Matrix Mitc4Shell::drillingTransformation4() const {
    Matrix transformation{kNodeCount, kGlobalDofCount};
    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        const std::size_t globalOffset = node * kGlobalDofsPerNode;
        for (std::size_t component = 0U; component < 3U; ++component) {
            transformation(node, globalOffset + 3U + component) =
                directors_[node][component];
        }
    }
    return transformation;
}

Matrix Mitc4Shell::directStrainDisplacement20(
    double xi,
    double eta,
    double zeta) const {
    return strainDisplacement(xi, eta, zeta, nullptr);
}

Matrix Mitc4Shell::covariantTyingShearSamples20() const {
    constexpr std::array<Vector3, 4> points{
        Vector3{0.0, -1.0, 0.0},
        Vector3{0.0, 1.0, 0.0},
        Vector3{-1.0, 0.0, 0.0},
        Vector3{1.0, 0.0, 0.0}};
    Matrix samples{4U, kPhysicalDofCount};
    for (std::size_t point = 0U; point < points.size(); ++point) {
        GeometryData geometry{};
        if (!evaluateGeometry(points[point][0], points[point][1], 0.0, geometry)) {
            throw std::logic_error{"Validated MITC4 tying geometry became invalid."};
        }
        const auto derivatives = basisDerivatives(
            points[point][0], points[point][1], 0.0);
        const std::size_t first = point < 2U ? 0U : 1U;
        for (std::size_t dof = 0U; dof < kPhysicalDofCount; ++dof) {
            const auto strain = covariantStrainColumn(
                geometry.covariant, derivatives[dof]);
            samples(point, dof) = strain[first][2U];
        }
    }
    return samples;
}

Matrix Mitc4Shell::strainDisplacement20(
    double xi,
    double eta,
    double zeta) const {
    const Matrix samples = covariantTyingShearSamples20();
    return strainDisplacement(xi, eta, zeta, &samples);
}

Matrix Mitc4Shell::planeStressConstitutive() const {
    const double factor = youngsModulus_ /
                          (1.0 - poissonRatio_ * poissonRatio_);
    Matrix constitutive{3U, 3U};
    constitutive(0U, 0U) = factor;
    constitutive(0U, 1U) = factor * poissonRatio_;
    constitutive(1U, 0U) = factor * poissonRatio_;
    constitutive(1U, 1U) = factor;
    constitutive(2U, 2U) = factor * (1.0 - poissonRatio_) * 0.5;
    return constitutive;
}

Matrix Mitc4Shell::materialConstitutive5() const {
    Matrix constitutive{5U, 5U};
    const Matrix planeStress = planeStressConstitutive();
    for (std::size_t row = 0U; row < 3U; ++row) {
        for (std::size_t column = 0U; column < 3U; ++column) {
            constitutive(row, column) = planeStress(row, column);
        }
    }
    const double shearModulus =
        youngsModulus_ / (2.0 * (1.0 + poissonRatio_));
    constitutive(3U, 3U) = kShearCorrection * shearModulus;
    constitutive(4U, 4U) = kShearCorrection * shearModulus;
    return constitutive;
}

Matrix Mitc4Shell::membraneSectionMatrix() const {
    return scaledMatrix(planeStressConstitutive(), thickness_);
}

Matrix Mitc4Shell::bendingSectionMatrix() const {
    return scaledMatrix(
        planeStressConstitutive(),
        thickness_ * thickness_ * thickness_ / 12.0);
}

Matrix Mitc4Shell::transverseShearSectionMatrix() const {
    const double shearModulus =
        youngsModulus_ / (2.0 * (1.0 + poissonRatio_));
    Matrix result{2U, 2U};
    result(0U, 0U) = kShearCorrection * shearModulus * thickness_;
    result(1U, 1U) = result(0U, 0U);
    return result;
}

Result<Mitc4Stiffness> Mitc4Shell::stiffness() const {
    const Matrix constitutive = materialConstitutive5();
    const Matrix tyingSamples = covariantTyingShearSamples20();
    Matrix physicalLocal{kPhysicalDofCount, kPhysicalDofCount};
    for (const auto& point : volumeQuadrature()) {
        GeometryData geometry{};
        if (!evaluateGeometry(
                point.naturalCoordinates[0],
                point.naturalCoordinates[1],
                point.naturalCoordinates[2],
                geometry)) {
            return stiffnessFailure(
                sourceLocation_, identity_,
                "Validated MITC4 quadrature geometry became invalid.");
        }
        const Matrix strain = strainDisplacement(
            point.naturalCoordinates[0],
            point.naturalCoordinates[1],
            point.naturalCoordinates[2],
            &tyingSamples);
        for (std::size_t row = 0U; row < kPhysicalDofCount; ++row) {
            for (std::size_t column = row;
                 column < kPhysicalDofCount; ++column) {
                double integrand = 0.0;
                for (std::size_t first = 0U; first < 5U; ++first) {
                    for (std::size_t second = 0U; second < 5U; ++second) {
                        integrand += strain(first, row) *
                                     constitutive(first, second) *
                                     strain(second, column);
                    }
                }
                const double contribution =
                    integrand * geometry.jacobian * point.weight;
                physicalLocal(row, column) += contribution;
                if (row != column) {
                    physicalLocal(column, row) += contribution;
                }
            }
        }
    }

    double drillingReference = (std::numeric_limits<double>::max)();
    bool hasDrillingReference = false;
    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        const std::size_t offset = node * kPhysicalDofsPerNode;
        for (std::size_t rotation = 3U; rotation < 5U; ++rotation) {
            const double diagonal = physicalLocal(offset + rotation, offset + rotation);
            if (std::isfinite(diagonal) && diagonal > 0.0) {
                drillingReference = (std::min)(drillingReference, diagonal);
                hasDrillingReference = true;
            }
        }
    }
    if (!hasDrillingReference) {
        return stiffnessFailure(
            sourceLocation_, identity_,
            "MITC4 drilling stabilization requires a finite positive physical "
            "tangent-rotation diagonal.");
    }
    if (!isFinite(physicalLocal)) {
        return stiffnessFailure(
            sourceLocation_, identity_,
            "MITC4 physical stiffness must contain only finite values.");
    }

    const double drillingStiffness = 1.0e-3 * drillingReference;
    if (!std::isfinite(drillingStiffness) || !(drillingStiffness > 0.0)) {
        return stiffnessFailure(
            sourceLocation_, identity_,
            "MITC4 drilling stiffness must be finite and positive.");
    }

    Matrix physicalGlobal = congruence(
        physicalLocal, physicalTransformation20());
    Matrix drillingLocal{kNodeCount, kNodeCount};
    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        drillingLocal(node, node) = drillingStiffness;
    }
    Matrix drillingGlobal = congruence(
        drillingLocal, drillingTransformation4());
    Matrix stabilizedGlobal{kGlobalDofCount, kGlobalDofCount};
    for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
        for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
            stabilizedGlobal(row, column) =
                physicalGlobal(row, column) + drillingGlobal(row, column);
        }
    }
    if (!isFinite(physicalGlobal) || !isFinite(drillingGlobal) ||
        !isFinite(stabilizedGlobal)) {
        return stiffnessFailure(
            sourceLocation_, identity_,
            "MITC4 transformed stiffness must contain only finite values.");
    }

    return Result<Mitc4Stiffness>::Success(Mitc4Stiffness{
        std::move(physicalLocal),
        std::move(physicalGlobal),
        std::move(drillingGlobal),
        std::move(stabilizedGlobal),
        drillingStiffness});
}

Result<Mitc4PhysicalRecovery> Mitc4Shell::recoverPhysical(
    const Vector& globalElementDisplacement24) const {
    if (globalElementDisplacement24.Size() != kGlobalDofCount) {
        return recoveryFailure(
            sourceLocation_, identity_,
            "MITC4 physical recovery requires exactly 24 global element DOFs.");
    }
    if (!isFinite(globalElementDisplacement24)) {
        return recoveryFailure(
            sourceLocation_, identity_,
            "MITC4 physical recovery displacement must be finite.");
    }

    const Vector physicalDisplacement =
        physicalTransformation20().Multiply(globalElementDisplacement24);
    const Matrix tyingSamples = covariantTyingShearSamples20();
    const Matrix constitutive = materialConstitutive5();
    const Matrix planeStress = planeStressConstitutive();
    const double gauss = 1.0 / std::sqrt(3.0);
    const std::array<std::array<double, 2>, 4> surfacePoints{
        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, 2> thicknessPoints{-1.0, 1.0};
    constexpr std::array<double, 3> sectionPositions{-1.0, 0.0, 1.0};

    Mitc4PhysicalRecovery recovery{};
    for (std::size_t surface = 0U; surface < surfacePoints.size(); ++surface) {
        auto& point = recovery.points[surface];
        point.naturalCoordinates = surfacePoints[surface];
        GeometryData midsurfaceGeometry{};
        if (!evaluateGeometry(
                point.naturalCoordinates[0], point.naturalCoordinates[1], 0.0,
                midsurfaceGeometry)) {
            return recoveryFailure(
                sourceLocation_, identity_,
                "MITC4 midsurface recovery geometry is invalid.");
        }
        point.localFrame = midsurfaceGeometry.frame;

        for (double thicknessSign : thicknessPoints) {
            const double zeta = thicknessSign * gauss;
            const Matrix strainMatrix = strainDisplacement(
                point.naturalCoordinates[0], point.naturalCoordinates[1], zeta,
                &tyingSamples);
            const Vector strain = strainMatrix.Multiply(physicalDisplacement);
            const Vector stress = constitutive.Multiply(strain);
            GeometryData geometry{};
            if (!evaluateGeometry(
                    point.naturalCoordinates[0], point.naturalCoordinates[1],
                    zeta, geometry)) {
                return recoveryFailure(
                    sourceLocation_, identity_,
                    "Validated MITC4 recovery geometry became invalid.");
            }

            for (std::size_t component = 0U; component < 3U; ++component) {
                point.generalizedStrain[component] += 0.5 * strain[component];
                point.generalizedStrain[3U + component] +=
                    3.0 * zeta * strain[component] / thickness_;
                point.sectionResultant[component] +=
                    0.5 * thickness_ * stress[component];
                point.sectionResultant[3U + component] +=
                    0.25 * thickness_ * thickness_ * zeta * stress[component];
            }
            for (std::size_t component = 0U; component < 2U; ++component) {
                point.generalizedStrain[6U + component] +=
                    0.5 * strain[3U + component];
                point.sectionResultant[6U + component] +=
                    0.5 * thickness_ * stress[3U + component];
            }
            recovery.strainEnergy +=
                0.5 * strain.Dot(stress) * geometry.jacobian;
        }

        for (std::size_t position = 0U;
             position < sectionPositions.size(); ++position) {
            GeometryData sectionGeometry{};
            if (!evaluateGeometry(
                    point.naturalCoordinates[0], point.naturalCoordinates[1],
                    sectionPositions[position], sectionGeometry)) {
                return recoveryFailure(
                    sourceLocation_, identity_,
                    "MITC4 section-position recovery geometry is invalid.");
            }
            const Vector strain = strainDisplacement(
                point.naturalCoordinates[0], point.naturalCoordinates[1],
                sectionPositions[position], &tyingSamples)
                                      .Multiply(physicalDisplacement);
            Vector inPlaneStrain{3U};
            for (std::size_t component = 0U; component < 3U; ++component) {
                inPlaneStrain[component] = strain[component];
            }
            const Vector stress = planeStress.Multiply(inPlaneStrain);
            for (std::size_t component = 0U; component < 3U; ++component) {
                point.inPlaneStress[position][component] = stress[component];
            }
        }
    }

    if (!std::isfinite(recovery.strainEnergy)) {
        return recoveryFailure(
            sourceLocation_, identity_,
            "MITC4 physical strain energy must be finite.");
    }
    for (const auto& point : recovery.points) {
        const auto finite = [](const auto& values) {
            return std::all_of(values.begin(), values.end(), [](double value) {
                return std::isfinite(value);
            });
        };
        if (!finite(point.generalizedStrain) ||
            !finite(point.sectionResultant) ||
            std::any_of(
                point.inPlaneStress.begin(), point.inPlaneStress.end(),
                [&finite](const auto& stress) { return !finite(stress); })) {
            return recoveryFailure(
                sourceLocation_, identity_,
                "MITC4 physical recovery values must be finite.");
        }
    }

    return Result<Mitc4PhysicalRecovery>::Success(std::move(recovery));
}

Mitc4Shell::Mitc4Shell(
    std::array<Vector3, 4> coordinates,
    std::array<Vector3, 4> directors,
    std::array<Vector3, 4> tangentA,
    std::array<Vector3, 4> tangentB,
    Vector3 normalCandidate,
    double thickness,
    double youngsModulus,
    double poissonRatio,
    SourceLocation sourceLocation,
    std::string identity)
    : coordinates_{std::move(coordinates)},
      directors_{std::move(directors)},
      tangentA_{std::move(tangentA)},
      tangentB_{std::move(tangentB)},
      normalCandidate_{std::move(normalCandidate)},
      thickness_{thickness},
      youngsModulus_{youngsModulus},
      poissonRatio_{poissonRatio},
      sourceLocation_{std::move(sourceLocation)},
      identity_{std::move(identity)} {}

bool Mitc4Shell::evaluateGeometry(
    double xi,
    double eta,
    double zeta,
    GeometryData& result) const noexcept {
    const auto shape = shapeFunctions(xi, eta);
    const Vector3 midsurfaceXi = derivativeSum(shape.xiDerivatives, coordinates_);
    const Vector3 midsurfaceEta = derivativeSum(shape.etaDerivatives, coordinates_);
    const Vector3 directorXi = derivativeSum(shape.xiDerivatives, directors_);
    const Vector3 directorEta = derivativeSum(shape.etaDerivatives, directors_);
    const Vector3 directorValue = weightedSum(shape.values, directors_);
    const double halfThickness = 0.5 * thickness_;

    result.covariant[0] = add(
        midsurfaceXi, scale(halfThickness * zeta, directorXi));
    result.covariant[1] = add(
        midsurfaceEta, scale(halfThickness * zeta, directorEta));
    result.covariant[2] = scale(halfThickness, directorValue);
    result.jacobian = dot(
        result.covariant[0], cross(result.covariant[1], result.covariant[2]));
    if (!isFinite(result.covariant[0]) || !isFinite(result.covariant[1]) ||
        !isFinite(result.covariant[2]) || !std::isfinite(result.jacobian) ||
        !(result.jacobian > 0.0)) {
        return false;
    }

    result.reciprocal[0] = scale(
        1.0 / result.jacobian,
        cross(result.covariant[1], result.covariant[2]));
    result.reciprocal[1] = scale(
        1.0 / result.jacobian,
        cross(result.covariant[2], result.covariant[0]));
    result.reciprocal[2] = scale(
        1.0 / result.jacobian,
        cross(result.covariant[0], result.covariant[1]));

    const Vector3 area = cross(midsurfaceXi, midsurfaceEta);
    const double directorNorm = norm(directorValue);
    if (!isFinite(area) || !isFinite(directorValue) ||
        !std::isfinite(directorNorm) || !(directorNorm > 0.0) ||
        !(dot(area, normalCandidate_) > 0.0)) {
        return false;
    }
    result.frame.e3 = scale(1.0 / directorNorm, directorValue);
    if (!(dot(area, result.frame.e3) > 0.0)) {
        return false;
    }
    const Vector3 e1Candidate = subtract(
        midsurfaceXi,
        scale(dot(midsurfaceXi, result.frame.e3), result.frame.e3));
    const double e1Norm = norm(e1Candidate);
    if (!isFinite(e1Candidate) || !std::isfinite(e1Norm) || !(e1Norm > 0.0)) {
        return false;
    }
    result.frame.e1 = scale(1.0 / e1Norm, e1Candidate);
    result.frame.e2 = cross(result.frame.e3, result.frame.e1);
    return isFinite(result.reciprocal[0]) && isFinite(result.reciprocal[1]) &&
           isFinite(result.reciprocal[2]) && isFinite(result.frame.e2);
}

std::array<std::array<Vector3, 3>, 20> Mitc4Shell::basisDerivatives(
    double xi,
    double eta,
    double zeta) const noexcept {
    const auto shape = shapeFunctions(xi, eta);
    const double halfThickness = 0.5 * thickness_;
    std::array<std::array<Vector3, 3>, kPhysicalDofCount> derivatives{};
    constexpr std::array<Vector3, 3> globalAxes{
        Vector3{1.0, 0.0, 0.0},
        Vector3{0.0, 1.0, 0.0},
        Vector3{0.0, 0.0, 1.0}};

    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        const std::size_t offset = node * kPhysicalDofsPerNode;
        for (std::size_t component = 0U; component < 3U; ++component) {
            derivatives[offset + component][0] =
                scale(shape.xiDerivatives[node], globalAxes[component]);
            derivatives[offset + component][1] =
                scale(shape.etaDerivatives[node], globalAxes[component]);
        }

        const Vector3 alphaDirection = scale(-halfThickness, tangentB_[node]);
        derivatives[offset + 3U][0] =
            scale(zeta * shape.xiDerivatives[node], alphaDirection);
        derivatives[offset + 3U][1] =
            scale(zeta * shape.etaDerivatives[node], alphaDirection);
        derivatives[offset + 3U][2] =
            scale(shape.values[node], alphaDirection);

        const Vector3 betaDirection = scale(halfThickness, tangentA_[node]);
        derivatives[offset + 4U][0] =
            scale(zeta * shape.xiDerivatives[node], betaDirection);
        derivatives[offset + 4U][1] =
            scale(zeta * shape.etaDerivatives[node], betaDirection);
        derivatives[offset + 4U][2] =
            scale(shape.values[node], betaDirection);
    }
    return derivatives;
}

Matrix Mitc4Shell::strainDisplacement(
    double xi,
    double eta,
    double zeta,
    const Matrix* tyingSamples) const {
    GeometryData geometry{};
    if (!evaluateGeometry(xi, eta, zeta, geometry)) {
        throw std::invalid_argument{"MITC4 strain geometry is invalid at the requested point."};
    }
    const auto derivatives = basisDerivatives(xi, eta, zeta);
    const Mitc4TyingWeights weights = tyingWeights(xi, eta);
    Matrix result{5U, kPhysicalDofCount};
    for (std::size_t dof = 0U; dof < kPhysicalDofCount; ++dof) {
        auto covariant = covariantStrainColumn(
            geometry.covariant, derivatives[dof]);
        if (tyingSamples != nullptr) {
            covariant[0U][2U] =
                weights.xiZeta[0] * (*tyingSamples)(0U, dof) +
                weights.xiZeta[1] * (*tyingSamples)(1U, dof);
            covariant[2U][0U] = covariant[0U][2U];
            covariant[1U][2U] =
                weights.etaZeta[0] * (*tyingSamples)(2U, dof) +
                weights.etaZeta[1] * (*tyingSamples)(3U, dof);
            covariant[2U][1U] = covariant[1U][2U];
        }
        const auto engineering = localEngineeringComponents(
            reconstructCartesianStrain(covariant, geometry.reciprocal),
            geometry.frame);
        for (std::size_t component = 0U; component < engineering.size(); ++component) {
            result(component, dof) = engineering[component];
        }
    }
    return result;
}

} // namespace fesa