fix: evaluate linear Hex8 response relative to translation

src/fesa/elements/hex8_element.cpp

@@ -230,8 +230,14 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
   if (!nonlinear) {
     response.material_tangent = stiffness_;
     response.tangent = stiffness_;
-    response.internal_force = stiffness_.Multiply(u);
-    response.strain_energy = 0.5 * u.Dot(response.internal_force);
+    // K0 annihilates uniform translation. Use the same relative field as
+    // strain recovery to avoid cancellation of a large rigid displacement.
+    Vector relative{24U};
+    for (std::size_t a = 0; a < 8; ++a)
+      for (std::size_t i = 0; i < 3; ++i)
+        relative[3U * a + i] = u[3U * a + i] - u[i];
+    response.internal_force = stiffness_.Multiply(relative);
+    response.strain_energy = 0.5 * relative.Dot(response.internal_force);
   }
   for (std::size_t g = 0; g < 8; ++g) {
     const auto& point = points_[g];

tests/unit/elements/hex8_element_test.cpp

@@ -107,6 +107,55 @@ TEST(Hex8Element, AffineTensorShearAndEightDistinctRecoveryPoints) {
                 .rows.size(),
             8U);
 }
+TEST(Hex8Element, UniformTranslationHasExactlyZeroLinearForceAndEnergy) {
+  auto element = Make(Brick(), 2.1e11, 0.3);
+  ASSERT_TRUE(element.HasValue());
+  fesa::Vector u{24U};
+  for (std::size_t a = 0; a < 8; ++a) {
+    u[3U * a] = 0.125;
+    u[3U * a + 1U] = -0.25;
+    u[3U * a + 2U] = 0.5;
+  }
+  const auto response = element.Value().EvaluateLinear(u);
+  ASSERT_TRUE(response.HasValue());
+  EXPECT_EQ(response.Value().strain_energy, 0.0);
+  for (std::size_t i = 0; i < 24; ++i)
+    EXPECT_EQ(response.Value().internal_force[i], 0.0);
+  for (const auto& point : response.Value().points)
+    for (std::size_t i = 0; i < 6; ++i) {
+      EXPECT_EQ(point.strain[i], 0.0);
+      EXPECT_EQ(point.cauchy_stress[i], 0.0);
+    }
+}
+TEST(Hex8Element, LinearResponseIsUnchangedByExactlyRepresentableTranslation) {
+  auto element = Make();
+  ASSERT_TRUE(element.HasValue());
+  fesa::Vector u{24U};
+  const auto x = Brick();
+  for (std::size_t a = 0; a < 8; ++a) {
+    u[3U * a] = x[a][0] / 64.0 + x[a][1] / 128.0;
+    u[3U * a + 1U] = x[a][2] / 256.0;
+  }
+  const auto original = element.Value().EvaluateLinear(u);
+  ASSERT_TRUE(original.HasValue());
+  for (std::size_t a = 0; a < 8; ++a) {
+    u[3U * a] += 1024.0;
+    u[3U * a + 1U] -= 512.0;
+    u[3U * a + 2U] += 256.0;
+  }
+  const auto translated = element.Value().EvaluateLinear(u);
+  ASSERT_TRUE(translated.HasValue());
+  EXPECT_EQ(translated.Value().strain_energy, original.Value().strain_energy);
+  for (std::size_t i = 0; i < 24; ++i)
+    EXPECT_EQ(translated.Value().internal_force[i],
+              original.Value().internal_force[i]);
+  for (std::size_t g = 0; g < 8; ++g) {
+    EXPECT_EQ(translated.Value().points[g].strain,
+              original.Value().points[g].strain);
+    EXPECT_EQ(translated.Value().points[g].cauchy_stress,
+              original.Value().points[g].cauchy_stress);
+  }
+}
 TEST(Hex8Element, FiniteRigidRotationHasZeroGreenEnergyAndForce) {
   auto element = Make();
   ASSERT_TRUE(element.HasValue());