feat(linear-static-mitc4-shell): step 4 - mitc4-stiffness-drilling

include/fesa/elements/mitc4_shell.hpp

@@ -5,6 +5,7 @@
 #include "fesa/model/model_types.hpp"
 
 #include <array>
+#include <string>
 
 namespace fesa {
 
@@ -30,8 +31,16 @@ struct Mitc4QuadraturePoint {
     double weight;
 };
 
-// Concrete small-rotation MITC4 kinematics and material value kernel. Global
-// equation ownership, drilling stiffness, recovery, and assembly remain outside.
+struct Mitc4Stiffness {
+    Matrix physicalLocal20;
+    Matrix physicalGlobal24;
+    Matrix drillingGlobal24;
+    Matrix stabilizedGlobal24;
+    double drillingStiffness;
+};
+
+// Concrete small-rotation MITC4 kinematics, constitutive, and stiffness kernel.
+// Global equation ownership, physical recovery, and assembly remain outside.
 class Mitc4Shell {
 public:
     static Result<Mitc4Shell> create(
@@ -63,6 +72,7 @@ public:
     [[nodiscard]] Matrix membraneSectionMatrix() const;
     [[nodiscard]] Matrix bendingSectionMatrix() const;
     [[nodiscard]] Matrix transverseShearSectionMatrix() const;
+    [[nodiscard]] Result<Mitc4Stiffness> stiffness() const;
 
 private:
     using Vector3 = std::array<double, 3>;
@@ -82,7 +92,9 @@ private:
         Vector3 normalCandidate,
         double thickness,
         double youngsModulus,
-        double poissonRatio);
+        double poissonRatio,
+        SourceLocation sourceLocation,
+        std::string identity);
 
     bool evaluateGeometry(
         double xi,
@@ -107,6 +119,8 @@ private:
     double thickness_;
     double youngsModulus_;
     double poissonRatio_;
+    SourceLocation sourceLocation_;
+    std::string identity_;
 };
 
 } // namespace fesa

src/fesa/elements/mitc4_shell.cpp

@@ -2,6 +2,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include <limits>
 #include <stdexcept>
 #include <string>
 #include <utility>
@@ -213,6 +214,55 @@ Matrix scaledMatrix(const Matrix& source, double factor) {
     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;
+}
+
+Result<Mitc4Stiffness> stiffnessFailure(
+    const SourceLocation& location,
+    const std::string& identity,
+    std::string message) {
+    return Result<Mitc4Stiffness>::failure(Status::failure(
+        FailureCategory::model,
+        {{Severity::error,
+          "invalid-shell-stiffness",
+          location,
+          "*ELEMENT",
+          identity,
+          std::move(message)}}));
+}
+
 } // namespace
 
 Result<Mitc4Shell> Mitc4Shell::create(
@@ -302,7 +352,9 @@ Result<Mitc4Shell> Mitc4Shell::create(
         normalCandidate,
         section.thickness,
         material.youngsModulus,
-        material.poissonRatio};
+        material.poissonRatio,
+        nodes[0]->location,
+        identity};
 
     const double shearModulus =
         material.youngsModulus / (2.0 * (1.0 + material.poissonRatio));
@@ -519,6 +571,108 @@ Matrix Mitc4Shell::transverseShearSectionMatrix() const {
     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});
+}
+
 Mitc4Shell::Mitc4Shell(
     std::array<Vector3, 4> coordinates,
     std::array<Vector3, 4> directors,
@@ -527,7 +681,9 @@ Mitc4Shell::Mitc4Shell(
     Vector3 normalCandidate,
     double thickness,
     double youngsModulus,
-    double poissonRatio)
+    double poissonRatio,
+    SourceLocation sourceLocation,
+    std::string identity)
     : coordinates_{std::move(coordinates)},
       directors_{std::move(directors)},
       tangentA_{std::move(tangentA)},
@@ -535,7 +691,9 @@ Mitc4Shell::Mitc4Shell(
       normalCandidate_{std::move(normalCandidate)},
       thickness_{thickness},
       youngsModulus_{youngsModulus},
-      poissonRatio_{poissonRatio} {}
+      poissonRatio_{poissonRatio},
+      sourceLocation_{std::move(sourceLocation)},
+      identity_{std::move(identity)} {}
 
 bool Mitc4Shell::evaluateGeometry(
     double xi,

tests/unit/elements/mitc4_shell_test.cpp

@@ -2,12 +2,15 @@
 
 #include <gtest/gtest.h>
 
+#include <algorithm>
 #include <array>
 #include <cmath>
 #include <cstddef>
 #include <cstdint>
 #include <limits>
+#include <stdexcept>
 #include <string>
+#include <vector>
 
 namespace {
 
@@ -124,6 +127,194 @@ bool hasPositiveCholeskyPivots(const fesa::Matrix& matrix) {
     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}),
@@ -381,3 +572,254 @@ TEST(Mitc4ShellKinematics, UsesOneFixedTwoByTwoByTwoQuadratureOrder) {
         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);
+}