fix(analysis): preserve solid static trial force recovery

src/fesa/analysis/solid_static_procedure.cpp

@@ -352,7 +352,7 @@ Status SolidStaticProcedureExecutor::ExecuteInto(
     auto recovery =
         Assemble(solids, dofs, frame.displacement, false, parallel_for_);
     if (!recovery.HasValue()) return recovery.GetStatus();
-    frame.internal_force = stiffness.Value().Multiply(frame.displacement);
+    frame.internal_force = std::move(recovery.Value().force);
     frame.external_force = std::move(load.Value());
     frame.solid_results = std::move(recovery.Value().points);
     frame.physical_strain_energy = recovery.Value().energy;

tests/unit/analysis/solid_static_procedure_test.cpp

@@ -566,6 +566,94 @@ TEST(SolidStaticProcedure,
   ASSERT_EQ(frame.solid_results.size(), 8U);
 }
 
+TEST(SolidStaticProcedure, UniformPrescribedTranslationHasExactlyZeroRecovery) {
+  // A rigid translation has no strain, stress, energy, internal force or RF.
+  // Binary-exact offsets independently expose cancellation in global K*u.
+  const double translation[3] = {2.75, -3.25, 1.5};
+  std::ostringstream boundary;
+  for (int node = 1; node <= 8; ++node)
+    for (int component = 0; component < 3; ++component)
+      boundary << node << ',' << component + 1 << ',' << component + 1 << ','
+               << translation[component] << '\n';
+  auto domain = Map(Deck(boundary.str()));
+  ASSERT_TRUE(domain.HasValue());
+  auto result = ExecuteSolid(domain.Value());
+  ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
+  const auto& frame = result.Value().state.CommittedState();
+  ASSERT_EQ(frame.displacement.Size(), 24U);
+  for (std::size_t i = 0; i < 24U; ++i) {
+    EXPECT_DOUBLE_EQ(frame.displacement[i], translation[i % 3U]);
+    EXPECT_DOUBLE_EQ(frame.internal_force[i], 0.0);
+    EXPECT_DOUBLE_EQ(frame.residual[i], 0.0);
+    EXPECT_DOUBLE_EQ(frame.reaction[i], 0.0);
+  }
+  EXPECT_DOUBLE_EQ(frame.physical_strain_energy, 0.0);
+  ASSERT_EQ(frame.solid_results.size(), 8U);
+  for (const auto& point : frame.solid_results) {
+    for (double strain : point.strain) EXPECT_DOUBLE_EQ(strain, 0.0);
+    for (double stress : point.cauchy_stress) EXPECT_DOUBLE_EQ(stress, 0.0);
+  }
+}
+
+TEST(SolidStaticProcedure, MixedStrainAndTranslationPreserveAnalyticRecovery) {
+  const double translation[3] = {2.75, -3.25, 1.5};
+  // Binary-exact affine strain: ex=1/32, engineering shear=1/16.
+  // Exact constant stress is [4.5,1.5,1.5,3,0,0] for E=120, nu=.25.
+  std::ostringstream translated_boundary;
+  translated_boundary.precision(17);
+  for (int a = 0; a < 8; ++a) {
+    const double u[3] = {translation[0] + kX[a][0] / 32 + kX[a][1] / 16,
+                         translation[1], translation[2]};
+    for (int component = 0; component < 3; ++component)
+      translated_boundary << a + 1 << ',' << component + 1 << ','
+                          << component + 1 << ',' << u[component] << '\n';
+  }
+  auto original_domain = Map(Deck(AffineBoundary(1.0 / 32, 1.0 / 16)));
+  auto translated_domain = Map(Deck(translated_boundary.str()));
+  ASSERT_TRUE(original_domain.HasValue());
+  ASSERT_TRUE(translated_domain.HasValue());
+  auto original = ExecuteSolid(original_domain.Value());
+  auto translated = ExecuteSolid(translated_domain.Value());
+  ASSERT_TRUE(original.HasValue()) << Message(original.GetStatus());
+  ASSERT_TRUE(translated.HasValue()) << Message(translated.GetStatus());
+  const auto& a = original.Value().state.CommittedState();
+  const auto& b = translated.Value().state.CommittedState();
+  const double stress[6] = {4.5, 1.5, 1.5, 3.0, 0.0, 0.0};
+  for (const auto* frame : {&a, &b}) {
+    EXPECT_NEAR(frame->physical_strain_energy, 0.1640625, 1e-13);
+    for (const auto& point : frame->solid_results) {
+      EXPECT_NEAR(point.strain[0], 1.0 / 32, 1e-15);
+      EXPECT_NEAR(point.strain[3], 1.0 / 32, 1e-15);
+      for (std::size_t c = 0; c < 6U; ++c)
+        EXPECT_NEAR(point.cauchy_stress[c], stress[c], 1e-13);
+    }
+    for (std::size_t node = 0; node < 8U; ++node) {
+      // Independent constant-stress surface traction integrals: each incident
+      // unit face contributes its outward normal times sigma/4 at this node.
+      const double nx = (2 * kX[node][0] - 1) / 4;
+      const double ny = (2 * kX[node][1] - 1) / 4;
+      const double nz = (2 * kX[node][2] - 1) / 4;
+      const double exact_force[3] = {4.5 * nx + 3 * ny, 3 * nx + 1.5 * ny,
+                                     1.5 * nz};
+      for (std::size_t c = 0; c < 3U; ++c) {
+        EXPECT_NEAR(frame->internal_force[3 * node + c], exact_force[c], 1e-13);
+        EXPECT_NEAR(frame->reaction[3 * node + c], exact_force[c], 1e-13);
+      }
+    }
+  }
+  EXPECT_DOUBLE_EQ(a.physical_strain_energy, b.physical_strain_energy);
+  for (std::size_t i = 0; i < 24U; ++i) {
+    EXPECT_DOUBLE_EQ(a.internal_force[i], b.internal_force[i]);
+    EXPECT_DOUBLE_EQ(a.reaction[i], b.reaction[i]);
+    EXPECT_DOUBLE_EQ(a.residual[i], b.residual[i]);
+  }
+  for (std::size_t g = 0; g < 8U; ++g) {
+    EXPECT_EQ(a.solid_results[g].strain, b.solid_results[g].strain);
+    EXPECT_EQ(a.solid_results[g].cauchy_stress,
+              b.solid_results[g].cauchy_stress);
+  }
+}
+
 TEST(SolidStaticProcedure, PureLargeRigidRotationEndsAtZeroStrainEnergy) {
   auto domain = Map(Deck(AffineBoundary(0.0, 0.0, 1.2), "", true));
   ASSERT_TRUE(domain.HasValue());