From ca806100c99913809a31861178e0c3b97493dc08 Mon Sep 17 00:00:00 2001
From: "KOKO\\Mimi" <kjkim2584@gmail.com>
Date: Fri, 9 Oct 2026 02:13:40 +0900
Subject: [PATCH] fix: evaluate linear Hex8 response relative to translation

---
 src/fesa/elements/hex8_element.cpp        | 10 ++++-
 tests/unit/elements/hex8_element_test.cpp | 49 +++++++++++++++++++++++
 2 files changed, 57 insertions(+), 2 deletions(-)

diff --git a/src/fesa/elements/hex8_element.cpp b/src/fesa/elements/hex8_element.cpp
index a42c897..0669033 100644
--- a/src/fesa/elements/hex8_element.cpp
+++ b/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];
diff --git a/tests/unit/elements/hex8_element_test.cpp b/tests/unit/elements/hex8_element_test.cpp
index 3aa7f90..8c512a2 100644
--- a/tests/unit/elements/hex8_element_test.cpp
+++ b/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());