feat: rebuild mitc4 stiffness drilling

PLAN.md

@@ -13,7 +13,7 @@ Every new agent session must read this file together with `PROGRESS.md` before p
 - If an item becomes obsolete, move it to `PROGRESS.md` with a short reason instead of silently deleting it.
 
 ## Current Objective
-Continue the new Phase 1 rebaseline plan in `phases/1-linear-static-mitc4-rebaseline`, starting with P1R-10 MITC4 stiffness, internal force, six-DOF transform, and drilling stabilization. The old `phases/1-linear-static-mitc4` path is historical and superseded by the paper-based MITC4 formulation reset.
+Continue the new Phase 1 rebaseline plan in `phases/1-linear-static-mitc4-rebaseline`, starting with P1R-11 MITC4 patch, locking-sensitivity, and benchmark tests. The old `phases/1-linear-static-mitc4` path is historical and superseded by the paper-based MITC4 formulation reset.
 
 ## Required Reading For New Agents
 1. `AGENTS.md`
@@ -36,7 +36,7 @@ Continue the new Phase 1 rebaseline plan in `phases/1-linear-static-mitc4-rebase
 ## Active Phase Files
 - Active phase directory: `phases/1-linear-static-mitc4-rebaseline`
 - Execute with: `python scripts/execute.py 1-linear-static-mitc4-rebaseline`
-- Step numbering is zero-based. `step0.md` is complete and recorded in `phases/1-linear-static-mitc4-rebaseline/step0-audit.md`; `step1.md` is complete and created the `quad_02_phase1.inp` normalized reference path; `step2.md` is complete and revalidated core harness guardrails; `step3.md` is complete and revalidated the Phase 1 parser/domain subset; `step4.md` is complete and strengthened validation/singular diagnostics; `step5.md` is complete and revalidated the DofManager/reaction foundation; `step6.md` is complete and revalidated the minimum result model plus displacement CSV comparator; `step7.md` is complete and revalidated MITC4 natural coordinates, tying points, center directors, and integration bases; `step8.md` is complete and revalidated degenerated-continuum displacement, direct covariant strain rows, and MITC shear tying rows; `step9.md` is complete and revalidated plane-stress material, convected-to-local transform, and `2 x 2 x 2` material integration scaffolding; `step15.md` is the independent evaluator closeout.
+- Step numbering is zero-based. `step0.md` is complete and recorded in `phases/1-linear-static-mitc4-rebaseline/step0-audit.md`; `step1.md` is complete and created the `quad_02_phase1.inp` normalized reference path; `step2.md` is complete and revalidated core harness guardrails; `step3.md` is complete and revalidated the Phase 1 parser/domain subset; `step4.md` is complete and strengthened validation/singular diagnostics; `step5.md` is complete and revalidated the DofManager/reaction foundation; `step6.md` is complete and revalidated the minimum result model plus displacement CSV comparator; `step7.md` is complete and revalidated MITC4 natural coordinates, tying points, center directors, and integration bases; `step8.md` is complete and revalidated degenerated-continuum displacement, direct covariant strain rows, and MITC shear tying rows; `step9.md` is complete and revalidated plane-stress material, convected-to-local transform, and `2 x 2 x 2` material integration scaffolding; `step10.md` is complete and revalidated MITC4 stiffness, internal force, six-DOF transform, and drilling stabilization; `step15.md` is the independent evaluator closeout.
 - Every step file contains a sprint contract with objective, required reading, scope, allowed files, explicit non-goals, tests to write first, reference artifacts, acceptance command, evaluator checklist, and handoff requirements.
 - Historical phase directory: `phases/1-linear-static-mitc4`
 - Historical phase status: blocked/superseded. Do not resume the old P1-15/P1-16 path unless the user explicitly requests recovery of that exact phase.
@@ -76,7 +76,7 @@ Each gate should be satisfied before moving to the next implementation band unle
 | G1 - Build and validation | satisfied | Build system, test framework, and `scripts/validate_workspace.py` run real checks. | Validation command output |
 | G2 - Parser and domain | satisfied | Parser subset revalidated in step 3; validation and singular diagnostics revalidated in step 4. | Parser acceptance/rejection tests, validation negative tests, and validation output |
 | G3 - DOF/math/results infrastructure | partial | Core aliases, DOF mapping, validation harness, model diagnostic context, DofManager, sparse-connectivity inputs, full-vector reaction formula, result model metadata, and displacement CSV comparator were revalidated in steps 2, 5, and 6; assembly remains for step 12. | P1R-02, P1R-05, and P1R-06 validation output |
-| G4 - MITC4 element readiness | partial | MITC4 formulation was rewritten from local papers; Steps 7 through 9 rebuilt geometry/director/local-basis scaffolding, displacement interpolation, direct covariant strain rows, MITC shear tying rows, plane-stress material, convected-to-local transform, and `2 x 2 x 2` material integration scaffolding, while stiffness/drilling and patch benchmarks remain for steps 10 and 11. | P1R-07 through P1R-09 validation output; future element tests |
+| G4 - MITC4 element readiness | partial | MITC4 formulation was rewritten from local papers; Steps 7 through 10 rebuilt geometry/director/local-basis scaffolding, displacement interpolation, direct covariant strain rows, MITC shear tying rows, plane-stress material, convected-to-local transform, `2 x 2 x 2` material integration scaffolding, stiffness/internal force, six-DOF transform, and drilling stabilization, while patch/benchmark coverage remains for step 11. | P1R-07 through P1R-10 validation output; future element tests |
 | G5 - End-to-end solver | reopened | Linear static path must be revalidated through steps 13 and 14 after the MITC4 rebuild and `quad_02` compatibility path. | Future integration/reference regression output |
 
 ## Phase 1 Implementation Milestones
@@ -91,7 +91,7 @@ All milestones are intended to become one or more self-contained sprint contract
 | P1R-07 | completed | MITC4 generator | Implement MITC4 geometry, node order, tying points, directors, and local bases. | none | Shape/basis/diagnostic tests |
 | P1R-08 | completed | MITC4 generator | Implement degenerated-continuum displacement, covariant strain rows, and MITC shear tying. | P1R-07 | Finite-difference and tying interpolation tests |
 | P1R-09 | completed | MITC4 generator | Implement material matrix, transform, and `2 x 2 x 2` integration scaffolding. | P1R-08 | Material/integration tests |
-| P1R-10 | pending | MITC4 generator | Assemble MITC4 stiffness/internal force with six-DOF transform and drilling stabilization. | P1R-09, P1R-05 | Symmetry, rigid body, drilling sensitivity tests |
+| P1R-10 | completed | MITC4 generator | Assemble MITC4 stiffness/internal force with six-DOF transform and drilling stabilization. | P1R-09, P1R-05 | Symmetry, rigid body, drilling sensitivity tests |
 | P1R-11 | pending | verification generator | Add MITC4 patch, locking-sensitivity, and benchmark tests. | P1R-10 | Membrane/bending/shear/twist/locking tests |
 | P1R-12 | pending | assembly generator | Rebuild assembly, solver adapter boundary, constrained solve, and full-vector RF recovery. | P1R-05, P1R-10 | Assembly and full-vector reaction tests |
 | P1R-13 | pending | analysis generator | Rebuild linear static workflow from input to U/RF result fields. | P1R-03, P1R-04, P1R-06, P1R-12 | End-to-end linear static tests |

PROGRESS.md

@@ -13,10 +13,37 @@ Every new agent session must read this file together with `PLAN.md` before plann
 - Do not remove history unless the user explicitly asks for archival cleanup.
 
 ## Current Status
-Phase 1 has a new rebaseline phase definition in `phases/1-linear-static-mitc4-rebaseline`. Steps 0 through 9 are complete. `quad_02_phase1.inp` is now the normalized Phase 1-compatible input path for the stored `quad_02` S4 reference pair, while the original `quad_02.inp` remains preserved unsupported provenance. Core numeric aliases, DOF mapping, validation harness, model diagnostic context, the Phase 1 parser/domain subset, validation/singular diagnostics, DofManager/reaction foundation, minimum result model metadata, displacement CSV comparator foundation, MITC4 geometry/director scaffolding, MITC4 displacement/strain/tying row scaffolding, and MITC4 material/transform/integration scaffolding have been revalidated. The old `phases/1-linear-static-mitc4` path is historical and superseded after the MITC4 formulation reset.
+Phase 1 has a new rebaseline phase definition in `phases/1-linear-static-mitc4-rebaseline`. Steps 0 through 10 are complete. `quad_02_phase1.inp` is now the normalized Phase 1-compatible input path for the stored `quad_02` S4 reference pair, while the original `quad_02.inp` remains preserved unsupported provenance. Core numeric aliases, DOF mapping, validation harness, model diagnostic context, the Phase 1 parser/domain subset, validation/singular diagnostics, DofManager/reaction foundation, minimum result model metadata, displacement CSV comparator foundation, MITC4 geometry/director scaffolding, MITC4 displacement/strain/tying row scaffolding, MITC4 material/transform/integration scaffolding, and MITC4 stiffness/drilling/internal-force scaffolding have been revalidated. The old `phases/1-linear-static-mitc4` path is historical and superseded after the MITC4 formulation reset.
 
 ## Completed Work
 
+### 2026-05-04 - P1R-10 MITC4 stiffness and drilling completed
+Author: Codex
+
+Changed files:
+- `include/fesa/fesa.hpp`
+- `tests/test_main.cpp`
+- `phases/1-linear-static-mitc4-rebaseline/index.json`
+- `PLAN.md`
+- `PROGRESS.md`
+
+Summary:
+- Added tests and implementation for MITC4 element stiffness accumulation using tied strain rows, transformed material matrices, and `2 x 2 x 2` integration samples.
+- Added a stiffness result API that exposes local physical stiffness, local drilling-stabilized stiffness, global stiffness, integration-point count, drilling reference diagonal, and drilling stiffness metadata.
+- Replaced the old default drilling scale with the documented `drilling_stiffness_scale = 1.0e-3` and applied it to the minimum positive physical local stiffness diagonal before global transformation.
+- Added diagnostics for invalid drilling scale and missing positive drilling reference diagonal.
+- Added tests for stiffness symmetry, drilling scale sensitivity, rigid translation zero energy, documented drilling-only energy, and `f_int_e = K_e u_e` consistency.
+- Routed `MITC4ElementKernel::stiffness` through the rebuilt formulation path so later assembly uses the Step 7 through Step 10 MITC4 helpers.
+
+Verification:
+- First ran `python scripts/validate_workspace.py` after adding Step 10 tests; it failed as expected because the MITC4 stiffness/drilling/internal-force APIs did not exist yet.
+- After implementation, `python scripts/validate_workspace.py` configured CMake, built `fesa_core` and `fesa_tests`, and ran CTest successfully.
+- CTest result: 1 test executable passed.
+
+Follow-up:
+- Continue with P1R-11 MITC4 patch, locking-sensitivity, and benchmark tests.
+- Step 10 intentionally did not tune the drilling scale to a reference case, add stress/resultant output, introduce reduced integration/hourglass logic, or rebuild sparse/global assembly beyond routing the existing kernel to the new element formulation.
+
 ### 2026-05-04 - P1R-09 MITC4 material integration scaffolding completed
 Author: Codex
 

include/fesa/fesa.hpp

@@ -1885,168 +1885,45 @@ inline LocalBasis computeLocalBasis(const std::array<Vec3, 4>& coordinates) {
   return {frame.v1, frame.v2, frame.vn};
 }
 
-struct NaturalDerivatives {
-  ShapeData shape;
-  std::array<Real, 4> dx{};
-  std::array<Real, 4> dy{};
-  Real det_j = 0.0;
-};
-
-inline NaturalDerivatives naturalDerivatives(const std::array<std::array<Real, 2>, 4>& xy, Real r, Real s) {
-  NaturalDerivatives data;
-  data.shape = shapeFunctions(r, s);
-  Real j00 = 0.0;
-  Real j01 = 0.0;
-  Real j10 = 0.0;
-  Real j11 = 0.0;
-  for (std::size_t i = 0; i < 4; ++i) {
-    j00 += data.shape.dr[i] * xy[i][0];
-    j01 += data.shape.ds[i] * xy[i][0];
-    j10 += data.shape.dr[i] * xy[i][1];
-    j11 += data.shape.ds[i] * xy[i][1];
-  }
-  data.det_j = j00 * j11 - j01 * j10;
-  if (std::fabs(data.det_j) < 1.0e-14) {
-    throw std::runtime_error("invalid shell element jacobian");
-  }
-  const Real inv00 = j11 / data.det_j;
-  const Real inv01 = -j01 / data.det_j;
-  const Real inv10 = -j10 / data.det_j;
-  const Real inv11 = j00 / data.det_j;
-  for (std::size_t i = 0; i < 4; ++i) {
-    data.dx[i] = inv00 * data.shape.dr[i] + inv01 * data.shape.ds[i];
-    data.dy[i] = inv10 * data.shape.dr[i] + inv11 * data.shape.ds[i];
-  }
-  return data;
-}
-
 struct ElementStiffnessOptions {
-  Real drilling_stiffness_scale = 1.0e-6;
+  Real drilling_stiffness_scale = 1.0e-3;
 };
 
-class MITC4ElementKernel {
+struct MITC4DrillingStabilizationResult {
- public:
+  DenseMatrix local_with_drilling;
-  DenseMatrix stiffness(const std::array<Vec3, 4>& coordinates, Real elastic_modulus, Real poisson_ratio, Real thickness,
+  Real reference_diagonal = 0.0;
-                        ElementStiffnessOptions options = {}) const {
+  Real drilling_stiffness = 0.0;
-    const LocalBasis basis = computeLocalBasis(coordinates);
+  Real drilling_stiffness_scale = 0.0;
-    std::array<std::array<Real, 2>, 4> xy{};
+  std::vector<Diagnostic> diagnostics;
-    for (std::size_t i = 0; i < 4; ++i) {
-      const Vec3 relative = coordinates[i] - coordinates[0];
-      xy[i] = {dot(relative, basis.e1), dot(relative, basis.e2)};
-    }
 
-    DenseMatrix local(24, 24);
+  bool ok() const {
-    const Real membrane_factor = elastic_modulus * thickness / (1.0 - poisson_ratio * poisson_ratio);
+    return !hasError(diagnostics);
-    const Real bending_factor = elastic_modulus * thickness * thickness * thickness / (12.0 * (1.0 - poisson_ratio * poisson_ratio));
+  }
-    const Real shear_factor = (5.0 / 6.0) * elastic_modulus * thickness / (2.0 * (1.0 + poisson_ratio));
+};
-    const Real gauss = 1.0 / std::sqrt(3.0);
-    const std::array<Real, 2> points = {-gauss, gauss};
 
-    for (Real r : points) {
+struct MITC4ElementStiffnessResult {
-      for (Real s : points) {
+  DenseMatrix local_without_drilling;
-        const NaturalDerivatives d = naturalDerivatives(xy, r, s);
+  DenseMatrix local_with_drilling;
-        DenseMatrix b(8, 24);
+  DenseMatrix global;
-        addMembraneB(b, d);
+  std::size_t integration_point_count = 0;
-        addBendingB(b, d);
+  Real drilling_reference_diagonal = 0.0;
-        addMitcShearB(b, xy, r, s);
+  Real drilling_stiffness = 0.0;
-        accumulateBtDB(local, b, membrane_factor, bending_factor, shear_factor, poisson_ratio, std::fabs(d.det_j));
+  Real drilling_stiffness_scale = 0.0;
-        addDrilling(local, d, elastic_modulus * thickness * options.drilling_stiffness_scale, std::fabs(d.det_j));
+  std::vector<Diagnostic> diagnostics;
-      }
-    }
 
-    return transformToGlobal(local, basis);
+  bool ok() const {
+    return !hasError(diagnostics);
   }
+};
 
- private:
+inline DenseMatrix mitc4LocalDofTransform(const MITC4Geometry& geometry) {
-  static void addMembraneB(DenseMatrix& b, const NaturalDerivatives& d) {
-    for (LocalIndex i = 0; i < 4; ++i) {
-      const LocalIndex c = 6 * i;
-      b(0, c + 0) = d.dx[static_cast<std::size_t>(i)];
-      b(1, c + 1) = d.dy[static_cast<std::size_t>(i)];
-      b(2, c + 0) = d.dy[static_cast<std::size_t>(i)];
-      b(2, c + 1) = d.dx[static_cast<std::size_t>(i)];
-    }
-  }
-
-  static void addBendingB(DenseMatrix& b, const NaturalDerivatives& d) {
-    for (LocalIndex i = 0; i < 4; ++i) {
-      const LocalIndex c = 6 * i;
-      b(3, c + 4) = -d.dx[static_cast<std::size_t>(i)];
-      b(4, c + 3) = d.dy[static_cast<std::size_t>(i)];
-      b(5, c + 3) = d.dx[static_cast<std::size_t>(i)];
-      b(5, c + 4) = -d.dy[static_cast<std::size_t>(i)];
-    }
-  }
-
-  static void addStandardShearRow(DenseMatrix& b, LocalIndex row, const NaturalDerivatives& d, bool gamma_xz, Real scale) {
-    for (LocalIndex i = 0; i < 4; ++i) {
-      const LocalIndex c = 6 * i;
-      if (gamma_xz) {
-        b(row, c + 2) += scale * d.dx[static_cast<std::size_t>(i)];
-        b(row, c + 4) += scale * d.shape.n[static_cast<std::size_t>(i)];
-      } else {
-        b(row, c + 2) += scale * d.dy[static_cast<std::size_t>(i)];
-        b(row, c + 3) -= scale * d.shape.n[static_cast<std::size_t>(i)];
-      }
-    }
-  }
-
-  static void addMitcShearB(DenseMatrix& b, const std::array<std::array<Real, 2>, 4>& xy, Real r, Real s) {
-    addStandardShearRow(b, 6, naturalDerivatives(xy, 0.0, -1.0), true, 0.5 * (1.0 - s));
-    addStandardShearRow(b, 6, naturalDerivatives(xy, 0.0, 1.0), true, 0.5 * (1.0 + s));
-    addStandardShearRow(b, 7, naturalDerivatives(xy, -1.0, 0.0), false, 0.5 * (1.0 - r));
-    addStandardShearRow(b, 7, naturalDerivatives(xy, 1.0, 0.0), false, 0.5 * (1.0 + r));
-  }
-
-  static void accumulateBtDB(DenseMatrix& k, const DenseMatrix& b, Real membrane_factor, Real bending_factor, Real shear_factor, Real poisson_ratio, Real det_j) {
-    std::array<std::array<Real, 8>, 8> d{};
-    d[0][0] = membrane_factor;
-    d[0][1] = poisson_ratio * membrane_factor;
-    d[1][0] = poisson_ratio * membrane_factor;
-    d[1][1] = membrane_factor;
-    d[2][2] = membrane_factor * (1.0 - poisson_ratio) / 2.0;
-    d[3][3] = bending_factor;
-    d[3][4] = poisson_ratio * bending_factor;
-    d[4][3] = poisson_ratio * bending_factor;
-    d[4][4] = bending_factor;
-    d[5][5] = bending_factor * (1.0 - poisson_ratio) / 2.0;
-    d[6][6] = shear_factor;
-    d[7][7] = shear_factor;
-    for (LocalIndex i = 0; i < 24; ++i) {
-      for (LocalIndex j = 0; j < 24; ++j) {
-        Real value = 0.0;
-        for (LocalIndex row = 0; row < 8; ++row) {
-          for (LocalIndex col = 0; col < 8; ++col) {
-            value += b(row, i) * d[static_cast<std::size_t>(row)][static_cast<std::size_t>(col)] * b(col, j);
-          }
-        }
-        k.add(i, j, value * det_j);
-      }
-    }
-  }
-
-  static void addDrilling(DenseMatrix& k, const NaturalDerivatives& d, Real drill_factor, Real det_j) {
-    std::array<Real, 24> b{};
-    for (LocalIndex i = 0; i < 4; ++i) {
-      const LocalIndex c = 6 * i;
-      b[static_cast<std::size_t>(c + 0)] = -0.5 * d.dy[static_cast<std::size_t>(i)];
-      b[static_cast<std::size_t>(c + 1)] = 0.5 * d.dx[static_cast<std::size_t>(i)];
-      b[static_cast<std::size_t>(c + 5)] = -d.shape.n[static_cast<std::size_t>(i)];
-    }
-    for (LocalIndex i = 0; i < 24; ++i) {
-      for (LocalIndex j = 0; j < 24; ++j) {
-        k.add(i, j, b[static_cast<std::size_t>(i)] * drill_factor * b[static_cast<std::size_t>(j)] * det_j);
-      }
-    }
-  }
-
-  static DenseMatrix transformToGlobal(const DenseMatrix& local, const LocalBasis& basis) {
   DenseMatrix transform(24, 24);
-    const std::array<Vec3, 3> axes = {basis.e1, basis.e2, basis.e3};
   for (LocalIndex node = 0; node < 4; ++node) {
+    const LocalIndex base = 6 * node;
+    const auto& frame = geometry.nodal_frames[static_cast<std::size_t>(node)];
+    const std::array<Vec3, 3> axes = {frame.v1, frame.v2, frame.vn};
     for (LocalIndex local_axis = 0; local_axis < 3; ++local_axis) {
       const Vec3 axis = axes[static_cast<std::size_t>(local_axis)];
-        const LocalIndex base = 6 * node;
       transform(base + local_axis, base + 0) = axis.x;
       transform(base + local_axis, base + 1) = axis.y;
       transform(base + local_axis, base + 2) = axis.z;
@@ -2055,13 +1932,95 @@ class MITC4ElementKernel {
       transform(base + 3 + local_axis, base + 5) = axis.z;
     }
   }
+  return transform;
+}
+
+inline MITC4StrainRows mitc4TransformStrainRowsToLocalDofs(const MITC4StrainRows& global_rows,
+                                                           const DenseMatrix& local_dof_transform) {
+  MITC4StrainRows local_rows;
+  local_rows.diagnostics = global_rows.diagnostics;
+  if (hasError(local_rows.diagnostics)) {
+    return local_rows;
+  }
+  for (std::size_t component = 0; component < 6; ++component) {
+    for (LocalIndex local_dof = 0; local_dof < 24; ++local_dof) {
+      Real value = 0.0;
+      for (LocalIndex global_dof = 0; global_dof < 24; ++global_dof) {
+        value += global_rows.rows[component][static_cast<std::size_t>(global_dof)] *
+                 local_dof_transform(local_dof, global_dof);
+      }
+      local_rows.rows[component][static_cast<std::size_t>(local_dof)] = value;
+    }
+  }
+  return local_rows;
+}
+
+inline void accumulateMITC4BtDB(DenseMatrix& stiffness, const MITC4StrainRows& rows,
+                                const MITC4MaterialMatrix& material, Real factor) {
+  for (LocalIndex i = 0; i < 24; ++i) {
+    for (LocalIndex j = 0; j < 24; ++j) {
+      Real value = 0.0;
+      for (std::size_t a = 0; a < 6; ++a) {
+        for (std::size_t b = 0; b < 6; ++b) {
+          value += rows.rows[a][static_cast<std::size_t>(i)] * material[a][b] *
+                   rows.rows[b][static_cast<std::size_t>(j)];
+        }
+      }
+      stiffness.add(i, j, value * factor);
+    }
+  }
+}
+
+inline Real mitc4MinimumPositivePhysicalLocalDiagonal(const DenseMatrix& local_without_drilling,
+                                                     Real tolerance = 1.0e-12) {
+  Real minimum = std::numeric_limits<Real>::infinity();
+  for (LocalIndex node = 0; node < 4; ++node) {
+    for (LocalIndex local_dof = 0; local_dof < 5; ++local_dof) {
+      const Real diagonal = local_without_drilling(6 * node + local_dof, 6 * node + local_dof);
+      if (isFinite(diagonal) && diagonal > tolerance) {
+        minimum = std::min(minimum, diagonal);
+      }
+    }
+  }
+  return minimum;
+}
+
+inline MITC4DrillingStabilizationResult mitc4ApplyDrillingStabilization(const DenseMatrix& local_without_drilling,
+                                                                        Real drilling_stiffness_scale,
+                                                                        Real tolerance = 1.0e-12) {
+  MITC4DrillingStabilizationResult result;
+  result.local_with_drilling = local_without_drilling;
+  result.drilling_stiffness_scale = drilling_stiffness_scale;
+  if (!isFinite(drilling_stiffness_scale) || drilling_stiffness_scale < 0.0) {
+    result.diagnostics.push_back(
+        mitc4Diagnostic("FESA-MITC4-DRILLING-SCALE", "MITC4 drilling stiffness scale must be non-negative and finite"));
+    return result;
+  }
+
+  result.reference_diagonal = mitc4MinimumPositivePhysicalLocalDiagonal(local_without_drilling, tolerance);
+  if (!isFinite(result.reference_diagonal)) {
+    result.diagnostics.push_back(mitc4Diagnostic("FESA-MITC4-DRILLING-REFERENCE",
+                                                 "MITC4 drilling stiffness reference diagonal is not positive"));
+    return result;
+  }
+
+  result.drilling_stiffness = drilling_stiffness_scale * result.reference_diagonal;
+  for (LocalIndex node = 0; node < 4; ++node) {
+    const LocalIndex gamma = 6 * node + 5;
+    result.local_with_drilling.add(gamma, gamma, result.drilling_stiffness);
+  }
+  return result;
+}
+
+inline DenseMatrix mitc4TransformLocalStiffnessToGlobal(const DenseMatrix& local_stiffness,
+                                                        const DenseMatrix& local_dof_transform) {
   DenseMatrix global(24, 24);
   for (LocalIndex i = 0; i < 24; ++i) {
     for (LocalIndex j = 0; j < 24; ++j) {
       Real value = 0.0;
       for (LocalIndex a = 0; a < 24; ++a) {
         for (LocalIndex b = 0; b < 24; ++b) {
-            value += transform(a, i) * local(a, b) * transform(b, j);
+          value += local_dof_transform(a, i) * local_stiffness(a, b) * local_dof_transform(b, j);
         }
       }
       global(i, j) = value;
@@ -2069,6 +2028,75 @@ class MITC4ElementKernel {
   }
   return global;
 }
+
+inline MITC4ElementStiffnessResult mitc4ElementStiffness(const std::array<Vec3, 4>& coordinates,
+                                                        Real elastic_modulus, Real poisson_ratio, Real thickness,
+                                                        ElementStiffnessOptions options = {}) {
+  MITC4ElementStiffnessResult result;
+  result.local_without_drilling = DenseMatrix(24, 24);
+  result.local_with_drilling = DenseMatrix(24, 24);
+  result.global = DenseMatrix(24, 24);
+  result.drilling_stiffness_scale = options.drilling_stiffness_scale;
+
+  const MITC4Geometry geometry = buildMITC4Geometry(coordinates, thickness);
+  appendDiagnostics(result.diagnostics, geometry.diagnostics);
+  if (hasError(result.diagnostics)) {
+    return result;
+  }
+
+  const MITC4MaterialIntegrationData integration =
+      mitc4BuildMaterialIntegrationData(geometry, elastic_modulus, poisson_ratio);
+  appendDiagnostics(result.diagnostics, integration.diagnostics);
+  if (hasError(result.diagnostics)) {
+    return result;
+  }
+  result.integration_point_count = integration.samples.size();
+
+  const DenseMatrix local_dof_transform = mitc4LocalDofTransform(geometry);
+  for (const MITC4MaterialIntegrationSample& sample : integration.samples) {
+    const MITC4StrainRows local_rows =
+        mitc4TransformStrainRowsToLocalDofs(sample.strain_rows, local_dof_transform);
+    appendDiagnostics(result.diagnostics, local_rows.diagnostics);
+    if (hasError(result.diagnostics)) {
+      return result;
+    }
+    const Real factor = std::fabs(sample.basis.jacobian) * sample.point.weight;
+    accumulateMITC4BtDB(result.local_without_drilling, local_rows, sample.convected_material, factor);
+  }
+
+  const auto drilling = mitc4ApplyDrillingStabilization(result.local_without_drilling, options.drilling_stiffness_scale);
+  appendDiagnostics(result.diagnostics, drilling.diagnostics);
+  result.local_with_drilling = drilling.local_with_drilling;
+  result.drilling_reference_diagonal = drilling.reference_diagonal;
+  result.drilling_stiffness = drilling.drilling_stiffness;
+  result.drilling_stiffness_scale = drilling.drilling_stiffness_scale;
+  if (hasError(result.diagnostics)) {
+    return result;
+  }
+
+  result.global = mitc4TransformLocalStiffnessToGlobal(result.local_with_drilling, local_dof_transform);
+  return result;
+}
+
+inline std::vector<Real> mitc4ElementInternalForce(const MITC4ElementStiffnessResult& stiffness,
+                                                  const std::vector<Real>& element_displacement) {
+  if (stiffness.global.rows() != static_cast<LocalIndex>(element_displacement.size())) {
+    throw std::runtime_error("MITC4 internal force size mismatch");
+  }
+  return stiffness.global.multiply(element_displacement);
+}
+
+class MITC4ElementKernel {
+ public:
+  DenseMatrix stiffness(const std::array<Vec3, 4>& coordinates, Real elastic_modulus, Real poisson_ratio, Real thickness,
+                        ElementStiffnessOptions options = {}) const {
+    const auto result = mitc4ElementStiffness(coordinates, elastic_modulus, poisson_ratio, thickness, options);
+    if (!result.ok()) {
+      throw std::runtime_error(result.diagnostics.empty() ? "invalid MITC4 stiffness"
+                                                          : result.diagnostics.front().message);
+    }
+    return result.global;
+  }
 };
 
 struct AssemblyResult {

phases/1-linear-static-mitc4-rebaseline/index.json

@@ -12,7 +12,7 @@
     { "step": 7, "name": "mitc4-geometry-directors", "status": "completed" },
     { "step": 8, "name": "mitc4-covariant-strain-tying", "status": "completed" },
     { "step": 9, "name": "mitc4-material-integration", "status": "completed" },
-    { "step": 10, "name": "mitc4-stiffness-drilling", "status": "pending" },
+    { "step": 10, "name": "mitc4-stiffness-drilling", "status": "completed" },
     { "step": 11, "name": "mitc4-patch-benchmark-tests", "status": "pending" },
     { "step": 12, "name": "assembly-sparse-solver-path", "status": "pending" },
     { "step": 13, "name": "linear-static-workflow", "status": "pending" },

tests/test_main.cpp

@@ -147,6 +147,19 @@ std::array<fesa::Real, 24> zeroElementDofs() {
   return values;
 }
 
+std::vector<fesa::Real> zeroElementVector() {
+  return std::vector<fesa::Real>(24, 0.0);
+}
+
+fesa::Real quadraticEnergy(const fesa::DenseMatrix& stiffness, const std::vector<fesa::Real>& values) {
+  const auto internal = stiffness.multiply(values);
+  fesa::Real energy = 0.0;
+  for (std::size_t i = 0; i < values.size(); ++i) {
+    energy += values[i] * internal[i];
+  }
+  return energy;
+}
+
 fesa::Domain singleElementValidationDomain() {
   fesa::Domain domain;
   domain.nodes[1] = {1, {0, 0, 0}};
@@ -1219,6 +1232,106 @@ FESA_TEST(mitc4_material_integration_samples_carry_tied_rows_basis_and_transform
   }
 }
 
+FESA_TEST(mitc4_element_stiffness_result_is_symmetric_and_uses_2x2x2_samples) {
+  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
+  const auto result = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2);
+  FESA_CHECK(result.ok());
+  FESA_CHECK(result.integration_point_count == 8);
+  FESA_CHECK(result.local_without_drilling.rows() == 24);
+  FESA_CHECK(result.local_with_drilling.rows() == 24);
+  FESA_CHECK(result.global.rows() == 24);
+  FESA_CHECK(result.global.cols() == 24);
+  for (fesa::LocalIndex i = 0; i < 24; ++i) {
+    for (fesa::LocalIndex j = 0; j < 24; ++j) {
+      FESA_CHECK_NEAR(result.local_without_drilling(i, j), result.local_without_drilling(j, i), 1.0e-8);
+      FESA_CHECK_NEAR(result.local_with_drilling(i, j), result.local_with_drilling(j, i), 1.0e-8);
+      FESA_CHECK_NEAR(result.global(i, j), result.global(j, i), 1.0e-8);
+    }
+  }
+}
+
+FESA_TEST(mitc4_drilling_stiffness_uses_minimum_positive_physical_local_diagonal) {
+  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
+  const auto result = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2);
+  FESA_CHECK(result.ok());
+  fesa::Real expected_reference = std::numeric_limits<fesa::Real>::infinity();
+  for (fesa::LocalIndex node = 0; node < 4; ++node) {
+    for (fesa::LocalIndex local_dof = 0; local_dof < 5; ++local_dof) {
+      const fesa::Real diagonal = result.local_without_drilling(6 * node + local_dof, 6 * node + local_dof);
+      if (std::isfinite(diagonal) && diagonal > 0.0) {
+        expected_reference = std::min(expected_reference, diagonal);
+      }
+    }
+  }
+  FESA_CHECK(std::isfinite(expected_reference));
+  FESA_CHECK_NEAR(result.drilling_reference_diagonal, expected_reference, 1.0e-10);
+  FESA_CHECK_NEAR(result.drilling_stiffness_scale, 1.0e-3, 1.0e-15);
+  FESA_CHECK_NEAR(result.drilling_stiffness, 1.0e-3 * expected_reference, 1.0e-10);
+  for (fesa::LocalIndex node = 0; node < 4; ++node) {
+    const fesa::LocalIndex gamma = 6 * node + 5;
+    FESA_CHECK_NEAR(result.local_with_drilling(gamma, gamma) - result.local_without_drilling(gamma, gamma),
+                    result.drilling_stiffness, 1.0e-12);
+  }
+
+  fesa::ElementStiffnessOptions options;
+  options.drilling_stiffness_scale = 2.0e-3;
+  const auto scaled = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2, options);
+  FESA_CHECK(scaled.ok());
+  FESA_CHECK_NEAR(scaled.drilling_reference_diagonal, result.drilling_reference_diagonal, 1.0e-10);
+  FESA_CHECK_NEAR(scaled.drilling_stiffness, 2.0 * result.drilling_stiffness, 1.0e-10);
+}
+
+FESA_TEST(mitc4_drilling_reports_invalid_scale_and_missing_reference_diagonal) {
+  fesa::DenseMatrix zero_physical(24, 24);
+  auto missing_reference = fesa::mitc4ApplyDrillingStabilization(zero_physical, 1.0e-3);
+  FESA_CHECK(!missing_reference.ok());
+  FESA_CHECK(fesa::containsDiagnostic(missing_reference.diagnostics, "FESA-MITC4-DRILLING-REFERENCE"));
+
+  auto invalid_scale = fesa::mitc4ApplyDrillingStabilization(zero_physical, -1.0);
+  FESA_CHECK(!invalid_scale.ok());
+  FESA_CHECK(fesa::containsDiagnostic(invalid_scale.diagnostics, "FESA-MITC4-DRILLING-SCALE"));
+}
+
+FESA_TEST(mitc4_rigid_translation_energy_is_zero_and_drilling_energy_is_documented) {
+  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
+  const auto result = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2);
+  FESA_CHECK(result.ok());
+
+  auto translate_x = zeroElementVector();
+  auto translate_y = zeroElementVector();
+  auto translate_z = zeroElementVector();
+  for (int node = 0; node < 4; ++node) {
+    translate_x[static_cast<std::size_t>(6 * node + 0)] = 1.0;
+    translate_y[static_cast<std::size_t>(6 * node + 1)] = 1.0;
+    translate_z[static_cast<std::size_t>(6 * node + 2)] = 1.0;
+  }
+  FESA_CHECK(std::fabs(quadraticEnergy(result.global, translate_x)) < 1.0e-8);
+  FESA_CHECK(std::fabs(quadraticEnergy(result.global, translate_y)) < 1.0e-8);
+  FESA_CHECK(std::fabs(quadraticEnergy(result.global, translate_z)) < 1.0e-8);
+
+  auto drilling = zeroElementVector();
+  for (int node = 0; node < 4; ++node) {
+    drilling[static_cast<std::size_t>(6 * node + 5)] = 1.0;
+  }
+  FESA_CHECK_NEAR(quadraticEnergy(result.global, drilling), 4.0 * result.drilling_stiffness, 1.0e-8);
+}
+
+FESA_TEST(mitc4_internal_force_is_stiffness_times_displacement_for_linear_phase1) {
+  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
+  const auto result = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2);
+  FESA_CHECK(result.ok());
+  std::vector<fesa::Real> displacement(24, 0.0);
+  for (std::size_t i = 0; i < displacement.size(); ++i) {
+    displacement[i] = 0.001 * static_cast<fesa::Real>(i + 1);
+  }
+  const auto expected = result.global.multiply(displacement);
+  const auto actual = fesa::mitc4ElementInternalForce(result, displacement);
+  FESA_CHECK(actual.size() == expected.size());
+  for (std::size_t i = 0; i < expected.size(); ++i) {
+    FESA_CHECK_NEAR(actual[i], expected[i], 1.0e-12);
+  }
+}
+
 FESA_TEST(mitc4_shape_functions_and_stiffness_baseline) {
   auto shape = fesa::shapeFunctions(0.25, -0.5);
   const fesa::Real sum = shape.n[0] + shape.n[1] + shape.n[2] + shape.n[3];