test: verify translated warped Hex8 linear variations

tests/unit/elements/hex8_element_test.cpp

@@ -537,6 +537,46 @@ TEST(Hex8Element, EnergyForceAndForceTangentDirectionalStepSweep) {
                            trial.Value().geometric_tangent(i, j));
     }
 }
+TEST(Hex8Element, TranslatedWarpedLinearResponseRetainsEnergyAndTangentVariation) {
+  auto element = Make(Warped());
+  ASSERT_TRUE(element.HasValue());
+  fesa::Vector u{24U}, direction{24U};
+  const std::array<double, 3> translation{1024.0, -512.0, 256.0};
+  for (std::size_t i = 0; i < 24; ++i) {
+    u[i] = translation[i % 3U] + .05 * std::sin(.7 * static_cast<double>(i));
+    direction[i] = std::cos(1.1 * static_cast<double>(i));
+  }
+  const auto trial = element.Value().EvaluateLinear(u);
+  ASSERT_TRUE(trial.HasValue());
+  const double derivative = direction.Dot(trial.Value().internal_force);
+  const auto tangent_direction = trial.Value().tangent.Multiply(direction);
+  ASSERT_NE(derivative, 0.0);
+  ASSERT_GT(tangent_direction.Norm(), 0.0);
+  std::array<double, 6> energy_error{}, force_error{};
+  const std::array<double, 6> steps{1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7};
+  for (std::size_t s = 0; s < steps.size(); ++s) {
+    auto plus = u, minus = u;
+    plus.Axpy(steps[s], direction);
+    minus.Axpy(-steps[s], direction);
+    const auto rp = element.Value().EvaluateLinear(plus);
+    const auto rm = element.Value().EvaluateLinear(minus);
+    ASSERT_TRUE(rp.HasValue());
+    ASSERT_TRUE(rm.HasValue());
+    energy_error[s] =
+        std::abs((rp.Value().strain_energy - rm.Value().strain_energy) /
+                     (2.0 * steps[s]) - derivative) / std::abs(derivative);
+    auto difference = rp.Value().internal_force;
+    difference.Axpy(-1.0, rm.Value().internal_force);
+    difference.Scale(1.0 / (2.0 * steps[s]));
+    difference.Axpy(-1.0, tangent_direction);
+    force_error[s] = difference.Norm() / tangent_direction.Norm();
+    std::cout << "translated-linear h=" << steps[s]
+              << " energy->f=" << energy_error[s]
+              << " f->K0=" << force_error[s] << '\n';
+  }
+  EXPECT_LT(*std::min_element(energy_error.begin(), energy_error.end()), 1e-7);
+  EXPECT_LT(*std::min_element(force_error.begin(), force_error.end()), 1e-7);
+}
 TEST(Hex8Element,
      InvalidInputsAndFailedTrialPreservePriorResponseAndReference) {
   EXPECT_FALSE(Make(Brick(), 0.0).HasValue());