fix: restore thickness factor in shell bending resultants

src/fesa/elements/nonlinear_mitc4_shell.cpp

@@ -1389,9 +1389,9 @@ Result<NonlinearMitc4Recovery> NonlinearMitc4Shell::EvaluateKinematics(
             3.0 * thickness_sign * gauss * point.strain[component] / thickness_;
         row.section_resultant[component] +=
             0.5 * thickness_ * point.stress[component];
-        row.section_resultant[3U + component] += 0.25 * thickness_ *
+        row.section_resultant[3U + component] +=
-                                                 thickness_sign * gauss *
+            0.25 * thickness_ * thickness_ * thickness_sign * gauss *
-                                                 point.stress[component];
+            point.stress[component];
       }
       for (std::size_t component = 0U; component < 2U; ++component) {
         row.generalized_strain[6U + component] +=

tests/unit/elements/nonlinear_mitc4_shell_test.cpp

@@ -91,6 +91,37 @@ TEST(NonlinearMitc4Shell, UndeformedStateHasZeroPhysicalRecovery) {
   ExpectZeroRecovery(recovery.Value());
 }
 
+TEST(NonlinearMitc4Recovery, BendingResultantsMatchStressThicknessMoments) {
+  const auto nodes = Nodes();
+  auto trial = ZeroTrial();
+  trial.rotation_vector = {{{0.12, -0.05, 0.0},
+                            {-0.02, 0.13, 0.0},
+                            {-0.12, 0.05, 0.0},
+                            {0.02, -0.13, 0.0}}};
+  for (const double thickness : {0.2, 0.5, 2.0}) {
+    SCOPED_TRACE(thickness);
+    auto section = Section();
+    section.thickness = thickness;
+    const auto shell = fesa::NonlinearMitc4Shell::Create(
+        NodePointers(nodes), Frames(), section, Material());
+    ASSERT_TRUE(shell.HasValue());
+    const auto recovery = shell.Value().RecoverState(trial);
+    ASSERT_TRUE(recovery.HasValue());
+    for (const auto& point : recovery.Value().gauss_points) {
+      for (std::size_t component = 0U; component < 3U; ++component) {
+        // Flat reference geometry gives quadratic stress in z. Its odd
+        // thickness moment is exact from independently sampled face stresses.
+        const double moment = thickness * thickness / 12.0 *
+                              (point.stress[2U].components[component] -
+                               point.stress[0U].components[component]);
+        ASSERT_GT(std::abs(moment), 0.0);
+        EXPECT_NEAR(point.section_resultant[3U + component], moment,
+                    1.0e-12 * std::abs(moment));
+      }
+    }
+  }
+}
+
 // NSNR-ELEM-KIN-001A
 TEST(NonlinearMitc4Shell, RejectsInvalidNodeAndInitialFrameContracts) {
   const auto nodes = Nodes();