fix(elements): retain Hex8 Green strain precision at small deformation

src/fesa/elements/hex8_element.cpp

@@ -258,12 +258,11 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
     Tensor strain{};
     for (std::size_t i = 0; i < 3; ++i)
       for (std::size_t j = 0; j < 3; ++j) {
-        if (nonlinear) {
+        strain[i][j] = 0.5 * (gradient[i][j] + gradient[j][i]);
+        // E = sym(H) + H^T*H/2 avoids subtracting I at tiny strain.
+        if (nonlinear)
           for (std::size_t k = 0; k < 3; ++k)
-            strain[i][j] += 0.5 * f[k][i] * f[k][j];
+            strain[i][j] += 0.5 * gradient[k][i] * gradient[k][j];
-          if (i == j) strain[i][j] -= 0.5;
-        } else
-          strain[i][j] = 0.5 * (gradient[i][j] + gradient[j][i]);
       }
     const double trace = strain[0][0] + strain[1][1] + strain[2][2];
     Tensor s{};

tests/unit/elements/hex8_element_test.cpp

@@ -156,6 +156,49 @@ TEST(Hex8Element, LinearResponseIsUnchangedByExactlyRepresentableTranslation) {
               original.Value().points[g].cauchy_stress);
   }
 }
+TEST(Hex8Element, TinySignedAffineStretchRetainsStvkStressForceAndEnergy) {
+  constexpr double youngs_modulus = 2.1e11;
+  constexpr double poisson_ratio = 0.3;
+  const double lambda = youngs_modulus * poisson_ratio /
+                        ((1.0 + poisson_ratio) * (1.0 - 2.0 * poisson_ratio));
+  const double mu = youngs_modulus / (2.0 * (1.0 + poisson_ratio));
+  const double axial_modulus = lambda + 2.0 * mu;
+  auto element = Make(Brick(), youngs_modulus, poisson_ratio);
+  ASSERT_TRUE(element.HasValue());
+  const auto x = Brick();
+  for (const double stretch : {1e-10, -1e-10, 1e-8, -1e-8}) {
+    SCOPED_TRACE(stretch);
+    fesa::Vector u{24U};
+    for (std::size_t a = 0; a < 8; ++a)
+      u[3U * a] = stretch * x[a][0];
+    const auto response = element.Value().EvaluateNonlinear(u);
+    ASSERT_TRUE(response.HasValue());
+    // Independent homogeneous diagonal F and reference face-traction integrals.
+    const double green = stretch + 0.5 * stretch * stretch;
+    const double sxx = axial_modulus * green;
+    const double syy = lambda * green;
+    const double expected_energy = 24.0 * axial_modulus * green * green;
+    EXPECT_NEAR(response.Value().strain_energy, expected_energy,
+                1e-12 * expected_energy);
+    for (const auto& point : response.Value().points) {
+      EXPECT_NEAR(point.strain[0], green, 1e-12 * std::abs(green));
+      EXPECT_NEAR(point.cauchy_stress[0], (1.0 + stretch) * sxx,
+                  1e-12 * std::abs(sxx));
+      EXPECT_NEAR(point.cauchy_stress[1], syy / (1.0 + stretch),
+                  1e-12 * std::abs(syy));
+      EXPECT_NEAR(point.cauchy_stress[2], syy / (1.0 + stretch),
+                  1e-12 * std::abs(syy));
+    }
+    for (std::size_t a = 0; a < 8; ++a) {
+      const std::array<double, 3> expected{
+          6.0 * kSigns[a][0] * (1.0 + stretch) * sxx,
+          3.0 * kSigns[a][1] * syy, 2.0 * kSigns[a][2] * syy};
+      for (std::size_t i = 0; i < 3; ++i)
+        EXPECT_NEAR(response.Value().internal_force[3U * a + i], expected[i],
+                    1e-12 * std::abs(expected[i]));
+    }
+  }
+}
 TEST(Hex8Element, FiniteRigidRotationHasZeroGreenEnergyAndForce) {
   auto element = Make();
   ASSERT_TRUE(element.HasValue());