diff --git a/include/fesa/elements/mitc4_shell.hpp b/include/fesa/elements/mitc4_shell.hpp
new file mode 100644
index 0000000..c8593a4
--- /dev/null
+++ b/include/fesa/elements/mitc4_shell.hpp
@@ -0,0 +1,112 @@
+#pragma once
+
+#include "fesa/core/status.hpp"
+#include "fesa/math/matrix.hpp"
+#include "fesa/model/model_types.hpp"
+
+#include <array>
+
+namespace fesa {
+
+struct Mitc4ShapeFunctions {
+    std::array<double, 4> values;
+    std::array<double, 4> xiDerivatives;
+    std::array<double, 4> etaDerivatives;
+};
+
+struct Mitc4LocalFrame {
+    std::array<double, 3> e1;
+    std::array<double, 3> e2;
+    std::array<double, 3> e3;
+};
+
+struct Mitc4TyingWeights {
+    std::array<double, 2> xiZeta;
+    std::array<double, 2> etaZeta;
+};
+
+struct Mitc4QuadraturePoint {
+    std::array<double, 3> naturalCoordinates;
+    double weight;
+};
+
+// Concrete small-rotation MITC4 kinematics and material value kernel. Global
+// equation ownership, drilling stiffness, recovery, and assembly remain outside.
+class Mitc4Shell {
+public:
+    static Result<Mitc4Shell> create(
+        std::array<const Node*, 4> nodes,
+        std::array<std::array<double, 3>, 4> initialDirectors,
+        const ShellSection& section,
+        const LinearElasticMaterial& material);
+
+    static Mitc4ShapeFunctions shapeFunctions(double xi, double eta) noexcept;
+    static Mitc4TyingWeights tyingWeights(double xi, double eta) noexcept;
+    static const std::array<Mitc4QuadraturePoint, 8>&
+    volumeQuadrature() noexcept;
+
+    [[nodiscard]] Mitc4LocalFrame localFrame(double xi, double eta) const;
+    [[nodiscard]] Matrix physicalTransformation20() const;
+    [[nodiscard]] Matrix drillingTransformation4() const;
+    [[nodiscard]] Matrix directStrainDisplacement20(
+        double xi,
+        double eta,
+        double zeta) const;
+    [[nodiscard]] Matrix covariantTyingShearSamples20() const;
+    [[nodiscard]] Matrix strainDisplacement20(
+        double xi,
+        double eta,
+        double zeta) const;
+
+    [[nodiscard]] Matrix planeStressConstitutive() const;
+    [[nodiscard]] Matrix materialConstitutive5() const;
+    [[nodiscard]] Matrix membraneSectionMatrix() const;
+    [[nodiscard]] Matrix bendingSectionMatrix() const;
+    [[nodiscard]] Matrix transverseShearSectionMatrix() const;
+
+private:
+    using Vector3 = std::array<double, 3>;
+
+    struct GeometryData {
+        std::array<Vector3, 3> covariant;
+        std::array<Vector3, 3> reciprocal;
+        Mitc4LocalFrame frame;
+        double jacobian;
+    };
+
+    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);
+
+    bool evaluateGeometry(
+        double xi,
+        double eta,
+        double zeta,
+        GeometryData& result) const noexcept;
+    std::array<std::array<Vector3, 3>, 20> basisDerivatives(
+        double xi,
+        double eta,
+        double zeta) const noexcept;
+    Matrix strainDisplacement(
+        double xi,
+        double eta,
+        double zeta,
+        const Matrix* tyingSamples) const;
+
+    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_;
+};
+
+} // namespace fesa
diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt
index f5bbcb7..1dc761f 100644
--- a/src/fesa/CMakeLists.txt
+++ b/src/fesa/CMakeLists.txt
@@ -13,6 +13,7 @@ add_library(
     core/diagnostic.cpp
     core/status.cpp
     elements/euler_beam_3d.cpp
+    elements/mitc4_shell.cpp
     fem/dof_manager.cpp
     io/abaqus/domain_mapper.cpp
     io/abaqus/input_reader.cpp
diff --git a/src/fesa/elements/mitc4_shell.cpp b/src/fesa/elements/mitc4_shell.cpp
new file mode 100644
index 0000000..2f6c843
--- /dev/null
+++ b/src/fesa/elements/mitc4_shell.cpp
@@ -0,0 +1,675 @@
+#include "fesa/elements/mitc4_shell.hpp"
+
+#include <algorithm>
+#include <cmath>
+#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.sourceLabelText;
+    }
+    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::model,
+        {{Severity::error,
+          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;
+}
+
+} // 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};
+
+    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;
+}
+
+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)
+    : 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} {}
+
+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
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index ede415f..6cf4378 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -13,6 +13,7 @@ add_executable(
     unit/core/source_identity_test.cpp
     unit/core/status_test.cpp
     unit/elements/euler_beam_3d_test.cpp
+    unit/elements/mitc4_shell_test.cpp
     unit/fem/dof_manager_test.cpp
     unit/math/matrix_test.cpp
     unit/math/sparse_matrix_test.cpp
diff --git a/tests/unit/elements/mitc4_shell_test.cpp b/tests/unit/elements/mitc4_shell_test.cpp
new file mode 100644
index 0000000..2591b4d
--- /dev/null
+++ b/tests/unit/elements/mitc4_shell_test.cpp
@@ -0,0 +1,383 @@
+#include "fesa/elements/mitc4_shell.hpp"
+
+#include <gtest/gtest.h>
+
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <string>
+
+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;
+}
+
+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);
+    }
+}