baram2584/FESADev
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
b37fddd4e6d037c48f8c2ddab85e577f17988809
FESADev/tests/unit/elements/hex8_high_low_test.cpp
T

771 lines
35 KiB
C++
Raw Blame History

#include "fesa/elements/hex8_element.h"
#include <gtest/gtest.h>
#include <array>
#include <atomic>
#include <cfenv>
#include <cmath>
#include <cstdint>
#include <cstring>
#include <future>
#include <limits>
#include <random>
#include <thread>
#include <utility>
#include <vector>
#include "fesa/materials/isotropic_linear_elastic_material.h"
#include "../../../src/fesa/math/high_low_arithmetic.h"
#include "../../../src/fesa/elements/hex8_geometry.h"
#include "../../../src/fesa/elements/hex8_pair_arithmetic.h"
#include "../../../src/fesa/elements/hex8_pair_simd.h"
#include "../math/high_low_arithmetic_original.h"
namespace {
using fesa::math_internal::HighLow;
using fesa::math_internal::Subtract;
std::uint64_t PairBits(double value) {
std::uint64_t bits = 0;
std::memcpy(&bits, &value, sizeof(bits));
return bits;
}
void ExpectOriginalProduct(HighLow left, HighLow right) {
const auto actual = fesa::hex8_internal::MultiplyPair(left, right);
const auto expected = fesa::high_low_original::Multiply(
{left.high, left.low}, {right.high, right.low});
EXPECT_EQ(PairBits(actual.high), PairBits(expected.high));
EXPECT_EQ(PairBits(actual.low), PairBits(expected.low));
const auto actual_sum = fesa::hex8_internal::AddPair(left, right);
const auto expected_sum = fesa::high_low_original::Add(
{left.high, left.low}, {right.high, right.low});
EXPECT_EQ(PairBits(actual_sum.high), PairBits(expected_sum.high));
EXPECT_EQ(PairBits(actual_sum.low), PairBits(expected_sum.low));
}
TEST(Hex8HighLow, ScalarCrossProductShortcutMatchesOriginalBits) {
ASSERT_TRUE(fesa::math_internal::ArithmeticEnvironmentSupported());
const double minimum = std::numeric_limits<double>::denorm_min();
const double maximum = std::numeric_limits<double>::max();
const double infinity = std::numeric_limits<double>::infinity();
const double nan = std::numeric_limits<double>::quiet_NaN();
std::vector<HighLow> pairs{{}, {-0.0, 0.0}, {0.0, -0.0}, {-0.0, -0.0},
{1, -0.0}, {-1, -0.0}, {1, 1}, {0, minimum}, {infinity, 0},
{0, infinity}, {nan, 0}, {0, nan}, {maximum}, {-maximum},
{minimum}, {-minimum}, {std::numeric_limits<double>::min()},
{0.5}, {std::nextafter(0.5, 1.0)},
{1, std::ldexp(1.0, -53)},
{std::nextafter(1.0, 2.0), -std::ldexp(1.0, -53)}};
for (const auto left : pairs)
for (const auto right : pairs) ExpectOriginalProduct(left, right);
for (int exponent = -1074; exponent <= 1023; ++exponent) {
const double value = std::ldexp(1.0, exponent);
for (double sign : {-1.0, 1.0})
for (const auto right : pairs) {
ExpectOriginalProduct(HighLow{sign * value}, right);
ExpectOriginalProduct(right, HighLow{sign * value});
}
}
std::mt19937_64 random(0x48e81234U);
for (int sample = 0; sample < 12000; ++sample) {
double left = 0;
double right = 0;
const auto left_bits = random();
const auto right_bits = random();
std::memcpy(&left, &left_bits, sizeof(left));
std::memcpy(&right, &right_bits, sizeof(right));
const auto pair = fesa::math_internal::TwoSum(left, right);
ExpectOriginalProduct(pair, HighLow{});
ExpectOriginalProduct(HighLow{}, pair);
ExpectOriginalProduct(pair, HighLow{right});
ExpectOriginalProduct(HighLow{right}, pair);
ExpectOriginalProduct(pair, pair);
}
}
fesa::hex8_internal::PairPacket FrozenMaterialPacket(
const fesa::hex8_internal::MaterialPacketInput& input) {
using Original = fesa::high_low_original::HighLow;
const auto original = [](HighLow value) {
return Original{value.high, value.low};
};
fesa::hex8_internal::PairPacket result;
for (std::size_t lane = 0; lane < 4; ++lane) {
const auto first = fesa::high_low_original::Multiply(
original(input.deviatoric_left[lane]),
original(input.deviatoric_right[lane]));
const auto second = fesa::high_low_original::Multiply(
original(input.shear_left[lane]), original(input.shear_right[lane]));
const auto third = fesa::high_low_original::Multiply(
original(input.tensor_left[lane]), original(input.gradient_dot[lane]));
const auto value = fesa::high_low_original::Add(
fesa::high_low_original::Add(first, second), third);
const auto accumulated = fesa::high_low_original::Add(
original(input.accumulated[lane]),
fesa::high_low_original::Multiply(value, original(input.weight)));
result[lane] = {accumulated.high, accumulated.low};
}
return result;
}
bool ExpectFrozenMaterialPacket(
const fesa::hex8_internal::MaterialPacketInput& input) {
const auto expected = FrozenMaterialPacket(input);
fesa::hex8_internal::PairPacket scalar, selected;
fesa::hex8_internal::EvaluateMaterialPacketScalar(input, &scalar);
const auto evaluator = fesa::hex8_internal::SelectMaterialPacketEvaluator();
evaluator(input, &selected);
for (std::size_t lane = 0; lane < 4; ++lane) {
SCOPED_TRACE(lane);
EXPECT_EQ(PairBits(scalar[lane].high), PairBits(expected[lane].high));
EXPECT_EQ(PairBits(scalar[lane].low), PairBits(expected[lane].low));
EXPECT_EQ(PairBits(selected[lane].high), PairBits(expected[lane].high));
EXPECT_EQ(PairBits(selected[lane].low), PairBits(expected[lane].low));
}
return evaluator != fesa::hex8_internal::EvaluateMaterialPacketScalar;
}
TEST(Hex8HighLow, FourColumnMaterialPacketMatchesFrozenScalarBits) {
ASSERT_TRUE(fesa::math_internal::ArithmeticEnvironmentSupported());
using fesa::hex8_internal::MaterialPacketInput;
RecordProperty("private_packet_backend",
fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
fesa::hex8_internal::EvaluateMaterialPacketScalar ? "scalar" : "avx2");
std::size_t packet_checks = 0, avx2_evaluator_checks = 0;
const double tiny = std::numeric_limits<double>::denorm_min();
const double maximum = std::numeric_limits<double>::max();
const std::vector<HighLow> boundary{{}, {-0.0, 0}, {0, -0.0},
{1}, {-1}, {1, -0.0}, {1, std::ldexp(1.0, -54)},
{-1, std::ldexp(1.0, -54)}, {tiny}, {-tiny},
{std::numeric_limits<double>::min()}, {maximum}, {-maximum},
{1, tiny}, {1, -tiny}, {std::nextafter(1.0, 2.0)},
{1, 1}, {0, tiny}, {std::numeric_limits<double>::infinity(), 0},
{std::numeric_limits<double>::quiet_NaN(), 0},
{std::numeric_limits<double>::signaling_NaN(), 0},
{0, std::numeric_limits<double>::signaling_NaN()}};
for (std::size_t first = 0; first < boundary.size(); ++first)
for (std::size_t second = 0; second < boundary.size(); ++second) {
MaterialPacketInput input;
input.weight = {0.125};
for (std::size_t lane = 0; lane < 4; ++lane) {
input.deviatoric_left[lane] = boundary[first];
input.deviatoric_right[lane] = boundary[second];
input.shear_left[lane] = {2, std::ldexp(1.0, -54)};
input.shear_right[lane] = {-0.5, std::ldexp(1.0, -56)};
input.tensor_left[lane] = {0.75};
input.gradient_dot[lane] = {-0.25};
input.accumulated[lane] = lane == 0 ? HighLow{} : HighLow{0.125};
}
avx2_evaluator_checks += ExpectFrozenMaterialPacket(input) ? 1U : 0U;
++packet_checks;
}
for (int exponent = -1074; exponent <= 1023; ++exponent) {
MaterialPacketInput input;
input.weight = {1};
for (std::size_t lane = 0; lane < 4; ++lane) {
input.deviatoric_left[lane] = {
std::ldexp(lane % 2U == 0 ? 1.0 : -1.0, exponent)};
input.deviatoric_right[lane] = {lane < 2 ? 0.5 : 2.0};
}
avx2_evaluator_checks += ExpectFrozenMaterialPacket(input) ? 1U : 0U;
++packet_checks;
}
std::mt19937_64 random(0x4c0116U);
for (std::size_t sample = 0; sample < 2048; ++sample) {
MaterialPacketInput input;
const auto pair = [&]() {
const double high = std::ldexp(
1.0 + static_cast<double>(random() % 256U) / 256.0,
static_cast<int>(random() % 201U) - 100);
const double signed_high = (random() & 1U) != 0 ? high : -high;
const double low = sample % 3U == 0 ? 0.0 : std::ldexp(high, -55);
return HighLow{signed_high, (random() & 1U) != 0 ? low : -low};
};
input.weight = sample % 2U == 0 ? HighLow{0.125} : pair();
for (std::size_t lane = 0; lane < 4; ++lane) {
input.deviatoric_left[lane] = pair();
input.deviatoric_right[lane] = pair();
input.shear_left[lane] = pair();
input.shear_right[lane] = pair();
input.tensor_left[lane] = pair();
input.gradient_dot[lane] = pair();
input.accumulated[lane] = sample % 4U == 0 ? HighLow{} : pair();
}
avx2_evaluator_checks += ExpectFrozenMaterialPacket(input) ? 1U : 0U;
++packet_checks;
}
RecordProperty("packet_oracle_checks", static_cast<int>(packet_checks));
RecordProperty("avx2_evaluator_packet_checks",
static_cast<int>(avx2_evaluator_checks));
}
TEST(Hex8HighLow, PacketSpecialLaneDoesNotChangeOtherColumns) {
fesa::hex8_internal::MaterialPacketInput input;
input.weight = {0.125};
for (std::size_t lane = 0; lane < 4; ++lane) {
input.deviatoric_left[lane] = {1, std::ldexp(1.0, -54)};
input.deviatoric_right[lane] = {2, std::ldexp(1.0, -55)};
input.shear_left[lane] = {0.5};
input.shear_right[lane] = {-0.25};
input.tensor_left[lane] = {0.75};
input.gradient_dot[lane] = {0.125};
}
const auto reference = FrozenMaterialPacket(input);
const std::vector<HighLow> special{{}, {-0.0, 0}, {1},
{std::numeric_limits<double>::max()},
{std::numeric_limits<double>::denorm_min()}, {1, 1},
{std::numeric_limits<double>::quiet_NaN(), 0}};
for (std::size_t lane = 0; lane < 4; ++lane)
for (const auto value : special) {
auto changed = input;
changed.deviatoric_left[lane] = value;
ExpectFrozenMaterialPacket(changed);
fesa::hex8_internal::PairPacket actual;
fesa::hex8_internal::SelectMaterialPacketEvaluator()(changed, &actual);
for (std::size_t other = 0; other < 4; ++other)
if (other != lane) {
EXPECT_EQ(PairBits(actual[other].high), PairBits(reference[other].high));
EXPECT_EQ(PairBits(actual[other].low), PairBits(reference[other].low));
}
}
}
TEST(Hex8HighLow, UniformDyadicPacketsExerciseEachLowZeroShape) {
std::size_t checks = 0, avx2_checks = 0;
for (int shape = 0; shape < 4; ++shape)
for (int low_weight = 0; low_weight < 2; ++low_weight)
for (int zero_accumulator = 0; zero_accumulator < 2; ++zero_accumulator) {
fesa::hex8_internal::MaterialPacketInput input;
input.weight = {0.5, low_weight != 0 ? std::ldexp(1.0, -57) : 0.0};
for (std::size_t lane = 0; lane < 4; ++lane) {
const bool left_low = (shape & 1) != 0;
const bool right_low = (shape & 2) != 0;
// Every lane has the same normalized nonzero dyadic inputs. Component
// LBE sums are >= -116 and all Add intermediates are bounded by 4.
// Thus each Multiply low-shape is uniform and its component products
// bypass zero/unsupported-tail fallback. No mixed-zero Add fallback
// is possible; the initial zero accumulator is uniform across lanes.
input.deviatoric_left[lane] = {1, left_low ? std::ldexp(1.0, -54) : 0};
input.deviatoric_right[lane] = {2, right_low ? std::ldexp(1.0, -54) : 0};
input.shear_left[lane] = {0.5, left_low ? std::ldexp(1.0, -56) : 0};
input.shear_right[lane] = {0.25, right_low ? std::ldexp(1.0, -57) : 0};
input.tensor_left[lane] = {0.75, left_low ? std::ldexp(1.0, -56) : 0};
input.gradient_dot[lane] = {0.125, right_low ? std::ldexp(1.0, -58) : 0};
input.accumulated[lane] = zero_accumulator != 0 ? HighLow{}
: HighLow{0.0625};
}
avx2_checks += ExpectFrozenMaterialPacket(input) ? 1U : 0U;
++checks;
}
RecordProperty("uniform_low_shape_checks", static_cast<int>(checks));
RecordProperty("avx2_uniform_low_shape_checks", static_cast<int>(avx2_checks));
}
TEST(Hex8HighLow, OverflowLaneIsNeutralizedBeforeResidualFma) {
const auto evaluator = fesa::hex8_internal::SelectMaterialPacketEvaluator();
if (evaluator == fesa::hex8_internal::EvaluateMaterialPacketScalar)
GTEST_SKIP() << "AVX2 masked environment is not selected; routing unverified.";
std::fenv_t original;
ASSERT_EQ(std::fegetenv(&original), 0);
for (std::size_t bad = 0; bad < 4; ++bad) {
fesa::hex8_internal::MaterialPacketInput input;
input.weight = {0.5};
for (std::size_t lane = 0; lane < 4; ++lane) {
input.deviatoric_left[lane] = {
lane == bad ? std::numeric_limits<double>::max() : 1};
input.deviatoric_right[lane] = {2};
input.shear_left[lane] = {0.5};
input.shear_right[lane] = {0.25};
input.tensor_left[lane] = {0.75};
input.gradient_dot[lane] = {0.125};
input.accumulated[lane] = {0.0625};
}
// All operands are finite/nonzero and max*2 has LBE sum 972, so Product
// reaches its vector multiply; only the post-product finite mask rejects it.
const auto expected = FrozenMaterialPacket(input);
EXPECT_EQ(std::feclearexcept(FE_ALL_EXCEPT), 0);
fesa::hex8_internal::PairPacket output;
evaluator(input, &output);
const bool extra_invalid = (std::fetestexcept(FE_INVALID) & FE_INVALID) != 0;
EXPECT_EQ(std::fesetenv(&original), 0);
EXPECT_FALSE(extra_invalid);
for (std::size_t lane = 0; lane < 4; ++lane) {
EXPECT_EQ(PairBits(output[lane].high), PairBits(expected[lane].high));
EXPECT_EQ(PairBits(output[lane].low), PairBits(expected[lane].low));
}
}
RecordProperty("avx2_product_overflow_positions", 4);
}
TEST(Hex8HighLow, PacketBackendSelectionRecoversAfterEnvironmentRestoration) {
const auto initial = fesa::hex8_internal::SelectMaterialPacketEvaluator();
const int rounding = std::fegetround();
ASSERT_EQ(std::fesetround(FE_UPWARD), 0);
EXPECT_EQ(fesa::hex8_internal::SelectMaterialPacketEvaluator(),
fesa::hex8_internal::EvaluateMaterialPacketScalar);
EXPECT_EQ(std::fesetround(rounding), 0);
#if defined(_M_X64) || defined(__SSE2__)
const auto control = _mm_getcsr();
_mm_setcsr(control | 0x8040U);
EXPECT_EQ(fesa::hex8_internal::SelectMaterialPacketEvaluator(),
fesa::hex8_internal::EvaluateMaterialPacketScalar);
_mm_setcsr(control);
_mm_setcsr(control & ~0x80U);
const auto unmasked = fesa::hex8_internal::SelectMaterialPacketEvaluator();
_mm_setcsr(control);
EXPECT_EQ(unmasked, fesa::hex8_internal::EvaluateMaterialPacketScalar);
#endif
EXPECT_EQ(fesa::hex8_internal::SelectMaterialPacketEvaluator(), initial);
}
TEST(Hex8HighLow, PacketDoesNotClearPriorMaskedStickyFlags) {
std::fenv_t original;
ASSERT_EQ(std::fegetenv(&original), 0);
const auto evaluator = fesa::hex8_internal::SelectMaterialPacketEvaluator();
ASSERT_EQ(std::feraiseexcept(FE_INEXACT), 0);
fesa::hex8_internal::MaterialPacketInput input;
input.weight = {1};
for (std::size_t lane = 0; lane < 4; ++lane) {
input.deviatoric_left[lane] = {1, std::ldexp(1.0, -54)};
input.deviatoric_right[lane] = {2, std::ldexp(1.0, -54)};
}
fesa::hex8_internal::PairPacket output;
evaluator(input, &output);
const bool preserved = (std::fetestexcept(FE_INEXACT) & FE_INEXACT) != 0;
EXPECT_EQ(std::fesetenv(&original), 0);
EXPECT_TRUE(preserved);
}
fesa::Hex8Element UnitBrick(double offset = 0, double nu = 0) {
fesa::ElementDofLayout layout;
layout.node_indices = {0, 1, 2, 3, 4, 5, 6, 7};
layout.components_per_node = {fesa::DofComponent::kUx,
fesa::DofComponent::kUy,
fesa::DofComponent::kUz};
fesa::Hex8Element::Coordinates coordinates{{{0, 0, 0}, {1, 0, 0},
{1, 1, 0}, {0, 1, 0},
{0, 0, 1}, {1, 0, 1},
{1, 1, 1}, {0, 1, 1}}};
for (auto& x : coordinates) x[0] += offset;
fesa::IsotropicLinearElasticMaterial material("test", 12, nu, {});
return fesa::Hex8Element::Create(0, layout, coordinates, material, 2).Value();
}
void ExpectSamePairResponse(const fesa::Hex8HighLowResponse& actual,
const fesa::Hex8HighLowResponse& expected) {
EXPECT_EQ(PairBits(actual.high.strain_energy),
PairBits(expected.high.strain_energy));
EXPECT_EQ(PairBits(actual.strain_energy_low),
PairBits(expected.strain_energy_low));
ASSERT_EQ(actual.high.tangent.Rows(), 24U);
ASSERT_EQ(actual.tangent_low.Rows(), 24U);
for (std::size_t i = 0; i < 24; ++i) {
EXPECT_EQ(PairBits(actual.high.internal_force[i]),
PairBits(expected.high.internal_force[i]));
EXPECT_EQ(PairBits(actual.internal_force_low[i]),
PairBits(expected.internal_force_low[i]));
for (std::size_t j = 0; j < 24; ++j) {
EXPECT_EQ(PairBits(actual.high.tangent(i, j)),
PairBits(expected.high.tangent(i, j)));
EXPECT_EQ(PairBits(actual.tangent_low(i, j)),
PairBits(expected.tangent_low(i, j)));
}
}
for (std::size_t point = 0; point < 8; ++point) {
EXPECT_EQ(actual.high.points[point].integration_point,
expected.high.points[point].integration_point);
for (std::size_t i = 0; i < 6; ++i) {
EXPECT_EQ(PairBits(actual.high.points[point].cauchy_stress[i]),
PairBits(expected.high.points[point].cauchy_stress[i]));
EXPECT_EQ(PairBits(actual.high.points[point].strain[i]),
PairBits(expected.high.points[point].strain[i]));
}
}
}
TEST(Hex8HighLow, ConcurrentFirstLinearTrialsAcrossCopiesAndMoveKeepBits) {
const auto element = UnitBrick();
auto copied = element;
const auto moved = std::move(copied);
fesa::Vector high(24), low_a(24), low_b(24);
for (std::size_t i = 0; i < 24; ++i) high[i] = 1;
for (auto a : {1, 2, 5, 6}) {
low_a[3U * a] = std::ldexp(1.0, -60);
low_b[3U * a] = -std::ldexp(1.0, -59);
}
const auto oracle = UnitBrick();
const auto expected_a = oracle.EvaluateLinearHighLow(high, low_a);
const auto expected_b = oracle.EvaluateLinearHighLow(high, low_b);
ASSERT_TRUE(expected_a.HasValue());
ASSERT_TRUE(expected_b.HasValue());
std::promise<void> start;
const auto signal = start.get_future().share();
std::atomic<std::size_t> arrived{0U};
std::array<std::future<fesa::Result<fesa::Hex8HighLowResponse>>, 4> calls;
for (std::size_t i = 0; i < calls.size(); ++i)
calls[i] = std::async(std::launch::async, [&, i] {
arrived.fetch_add(1U, std::memory_order_release);
signal.wait();
return (i % 2U == 0U ? element : moved)
.EvaluateLinearHighLow(high, i < 2U ? low_a : low_b);
});
while (arrived.load(std::memory_order_acquire) != calls.size())
std::this_thread::yield();
start.set_value();
for (std::size_t i = 0; i < calls.size(); ++i) {
const auto response = calls[i].get();
ASSERT_TRUE(response.HasValue());
ExpectSamePairResponse(response.Value(),
i < 2U ? expected_a.Value() : expected_b.Value());
}
}
TEST(Hex8HighLow, FirstLinearPreparationRecoversAfterTrialAndEnvironmentFailure) {
ASSERT_TRUE(fesa::math_internal::ArithmeticEnvironmentSupported());
const auto element = UnitBrick();
fesa::Vector high(24), low(24);
EXPECT_FALSE(element.EvaluateLinearHighLow(fesa::Vector(23), low).HasValue());
high[0] = std::numeric_limits<double>::infinity();
EXPECT_FALSE(element.EvaluateLinearHighLow(high, low).HasValue());
high[0] = 0;
struct RestoreEnvironment {
int rounding = std::fegetround();
#if defined(_M_X64) || defined(__SSE2__)
unsigned csr = _mm_getcsr();
#endif
~RestoreEnvironment() {
std::fesetround(rounding);
#if defined(_M_X64) || defined(__SSE2__)
_mm_setcsr(csr);
#endif
}
} restore;
ASSERT_EQ(std::fesetround(FE_DOWNWARD), 0);
const auto bad_environment = element.EvaluateLinearHighLow(high, low);
ASSERT_FALSE(bad_environment.HasValue());
EXPECT_EQ(bad_environment.GetStatus().Diagnostics()[0].message,
"Hex8 high+low requires nearest rounding and gradual underflow.");
ASSERT_EQ(std::fesetround(FE_TONEAREST), 0);
#if defined(_M_X64) || defined(__SSE2__)
_mm_setcsr(restore.csr | 0x8040U);
EXPECT_FALSE(element.EvaluateLinearHighLow(high, low).HasValue());
_mm_setcsr(restore.csr);
#endif
const auto recovered = element.EvaluateLinearHighLow(high, low);
ASSERT_TRUE(recovered.HasValue());
high[3] = -2;
EXPECT_FALSE(element.EvaluateLinearHighLow(high, low).HasValue());
high[3] = 0;
const auto after_failure = element.EvaluateLinearHighLow(high, low);
ASSERT_TRUE(after_failure.HasValue());
ExpectSamePairResponse(after_failure.Value(), recovered.Value());
}
TEST(Hex8HighLow, LowOnlyAffineStrainSurvivesUniformHighTranslation) {
auto brick = UnitBrick();
fesa::Vector high(24), low(24);
const double strain = std::ldexp(1.0, -60);
for (std::size_t a = 0; a < 8; ++a) high[3 * a] = 1;
for (auto a : {1, 2, 5, 6}) low[3 * a] = strain;
const auto response = brick.EvaluateLinearHighLow(high, low);
ASSERT_TRUE(response.HasValue());
for (std::size_t a = 0; a < 8; ++a) {
const bool positive = a == 1 || a == 2 || a == 5 || a == 6;
EXPECT_NEAR(response.Value().high.internal_force[3 * a],
(positive ? 3 : -3) * strain, 1e-13 * strain);
}
EXPECT_NEAR(response.Value().high.strain_energy, 6 * strain * strain,
1e-13 * strain * strain);
for (const auto& row : response.Value().high.points)
EXPECT_NEAR(row.strain[0], strain, 1e-13 * strain);
}
TEST(Hex8HighLow, NonlinearLowGradientProducesSecondOrderForce) {
auto brick = UnitBrick();
fesa::Vector high(24), low(24);
const double strain = std::ldexp(1.0, -60);
for (std::size_t a = 0; a < 8; ++a) high[3 * a] = 1;
for (auto a : {1, 2, 5, 6}) low[3 * a] = strain;
const auto linear = brick.EvaluateLinearHighLow(high, low);
const auto nonlinear = brick.EvaluateNonlinearHighLow(high, low);
ASSERT_TRUE(linear.HasValue());
ASSERT_TRUE(nonlinear.HasValue());
const auto& l = linear.Value();
const auto& n = nonlinear.Value();
const auto increment = Subtract(
HighLow{n.high.internal_force[3], n.internal_force_low[3]},
HighLow{l.high.internal_force[3], l.internal_force_low[3]});
EXPECT_NEAR(increment.high, 4.5 * strain * strain,
1e-12 * strain * strain);
}
TEST(Hex8HighLow, VolumetricAndTransverseTermsRetainLowStrain) {
auto brick = UnitBrick(0, 0.25);
fesa::Vector high(24), low(24);
const double strain = std::ldexp(1.0, -60);
for (std::size_t a = 0; a < 8; ++a) high[3*a] = 1;
for (auto a : {1, 2, 5, 6}) low[3*a] = strain;
const auto linear = brick.EvaluateLinearHighLow(high, low);
const auto nonlinear = brick.EvaluateNonlinearHighLow(high, low);
ASSERT_TRUE(linear.HasValue()); ASSERT_TRUE(nonlinear.HasValue());
// Uniform axial strain: lambda=mu=4.8, each nodal face integral is 1/4.
EXPECT_NEAR(linear.Value().high.internal_force[3], 3.6*strain, 1e-13*strain);
EXPECT_NEAR(linear.Value().high.internal_force[7], 1.2*strain, 1e-13*strain);
EXPECT_NEAR(linear.Value().high.strain_energy, 7.2*strain*strain,
1e-13*strain*strain);
const auto quadratic = Subtract(
{nonlinear.Value().high.internal_force[3], nonlinear.Value().internal_force_low[3]},
{linear.Value().high.internal_force[3], linear.Value().internal_force_low[3]});
EXPECT_NEAR(quadratic.high, 5.4*strain*strain, 1e-12*strain*strain);
}
TEST(Hex8HighLow, RejectsInvalidLowShapeAndOverlappingComponents) {
auto brick = UnitBrick();
fesa::Vector high(24), low(24);
EXPECT_FALSE(brick.EvaluateLinearHighLow(high, fesa::Vector(23)).HasValue());
high[0] = 1;
low[0] = 1;
EXPECT_FALSE(brick.EvaluateNonlinearHighLow(high, low).HasValue());
}
TEST(Hex8HighLow, LinearAndNonlinearZeroStateHaveSameFullHessian) {
auto brick = UnitBrick();
fesa::Vector zero(24);
const auto linear = brick.EvaluateLinearHighLow(zero, zero);
const auto nonlinear = brick.EvaluateNonlinearHighLow(zero, zero);
ASSERT_TRUE(linear.HasValue());
ASSERT_TRUE(nonlinear.HasValue());
for (std::size_t i = 0; i < 24; ++i)
for (std::size_t j = 0; j < 24; ++j) {
EXPECT_EQ(linear.Value().high.tangent(i, j),
nonlinear.Value().high.tangent(i, j));
EXPECT_EQ(linear.Value().tangent_low(i, j),
nonlinear.Value().tangent_low(i, j));
}
}
TEST(Hex8HighLow, CachedLinearHessianIsIndependentOfTrialAndPriorFailure) {
auto brick = UnitBrick(0, 0.25);
fesa::Vector zero(24), high(24), low(24);
high[3] = 0.1;
low[3] = std::ldexp(1.0, -60);
const auto reference = brick.EvaluateLinearHighLow(zero, zero);
ASSERT_TRUE(reference.HasValue());
EXPECT_FALSE(brick.EvaluateLinearHighLow(high, fesa::Vector(23)).HasValue());
const auto strained = brick.EvaluateLinearHighLow(high, low);
ASSERT_TRUE(strained.HasValue());
for (std::size_t i = 0; i < 24; ++i)
for (std::size_t j = 0; j < 24; ++j) {
EXPECT_EQ(reference.Value().high.tangent(i,j), strained.Value().high.tangent(i,j));
EXPECT_EQ(reference.Value().tangent_low(i,j), strained.Value().tangent_low(i,j));
}
const auto nonlinear_zero = brick.EvaluateNonlinearHighLow(zero, zero);
ASSERT_TRUE(nonlinear_zero.HasValue());
for (std::size_t i = 0; i < 24; ++i)
for (std::size_t j = 0; j < 24; ++j) {
EXPECT_EQ(reference.Value().high.tangent(i,j), nonlinear_zero.Value().high.tangent(i,j));
EXPECT_EQ(reference.Value().tangent_low(i,j), nonlinear_zero.Value().tangent_low(i,j));
}
}
TEST(Hex8HighLow, CoupledBricksTransmitLowStrainAcrossSharedFace) {
const auto left = UnitBrick(), right = UnitBrick(1);
fesa::Vector high(24), left_low(24), right_low(24);
const double strain = std::ldexp(1.0, -60);
for (std::size_t a = 0; a < 8; ++a) {
high[3 * a] = 1;
const bool positive = a == 1 || a == 2 || a == 5 || a == 6;
left_low[3 * a] = positive ? strain : 0;
right_low[3 * a] = positive ? 2 * strain : strain;
}
const auto l = left.EvaluateLinearHighLow(high, left_low);
const auto r = right.EvaluateLinearHighLow(high, right_low);
ASSERT_TRUE(l.HasValue());
ASSERT_TRUE(r.HasValue());
const std::array<std::size_t, 4> face_left{1, 2, 5, 6};
const std::array<std::size_t, 4> face_right{0, 3, 4, 7};
for (std::size_t a = 0; a < 4; ++a) {
const auto i = 3 * face_left[a], j = 3 * face_right[a];
const auto assembled = fesa::math_internal::Add(
{l.Value().high.internal_force[i], l.Value().internal_force_low[i]},
{r.Value().high.internal_force[j], r.Value().internal_force_low[j]});
EXPECT_NEAR(assembled.high, 0, 1e-13 * strain);
}
EXPECT_NEAR(l.Value().high.strain_energy + r.Value().high.strain_energy,
12 * strain * strain, 1e-13 * strain * strain);
}
TEST(Hex8HighLow, NodeZeroSignedDirectionMatchesForceAndEnergyDerivatives) {
using namespace fesa::math_internal;
auto brick = UnitBrick();
fesa::Vector high(24), low(24);
for (std::size_t a = 0; a < 8; ++a) {
high[3 * a] = 1;
low[3 * a] = (a == 1 || a == 2 || a == 5 || a == 6)
? std::ldexp(1.0, -60) : 0;
}
const auto baseline = brick.EvaluateNonlinearHighLow(high, low);
ASSERT_TRUE(baseline.HasValue());
const auto& base = baseline.Value();
for (int exponent : {-65, -70, -75}) {
const double step = std::ldexp(1.0, exponent);
auto plus_high = high, minus_high = high;
auto plus_low = low, minus_low = low;
const auto plus = Add({high[0], low[0]}, HighLow{step});
const auto minus = Subtract({high[0], low[0]}, HighLow{step});
plus_high[0] = plus.high; plus_low[0] = plus.low;
minus_high[0] = minus.high; minus_low[0] = minus.low;
const auto positive = brick.EvaluateNonlinearHighLow(plus_high, plus_low);
const auto negative = brick.EvaluateNonlinearHighLow(minus_high, minus_low);
ASSERT_TRUE(positive.HasValue());
ASSERT_TRUE(negative.HasValue());
for (std::size_t i = 0; i < 24; ++i) {
const auto derivative = Divide(Subtract(
{positive.Value().high.internal_force[i], positive.Value().internal_force_low[i]},
{negative.Value().high.internal_force[i], negative.Value().internal_force_low[i]}),
HighLow{2 * step});
const auto error = Subtract(derivative,
{base.high.tangent(i, 0), base.tangent_low(i, 0)});
EXPECT_LE(std::abs(error.high), 1e-25);
}
const auto gradient = Divide(Subtract(
{positive.Value().high.strain_energy, positive.Value().strain_energy_low},
{negative.Value().high.strain_energy, negative.Value().strain_energy_low}),
HighLow{2 * step});
const auto energy_error = Subtract(gradient,
{base.high.internal_force[0], base.internal_force_low[0]});
// Centered energy differences retain the StVK cubic term: O(step^2).
// For this unit E=12 brick its directional coefficient is bounded by 32.
EXPECT_LE(std::abs(energy_error.high), 32 * step * step + 1e-44);
}
}
TEST(Hex8HighLow, GeometryCertificateIncludesLowOnlyCurrentCellWidth) {
fesa::Hex8Element::Coordinates coordinates{{{0, 0, 0}, {1, 0, 0},
{1, 1, 0}, {0, 1, 0},
{0, 0, 1}, {1, 0, 1},
{1, 1, 1}, {0, 1, 1}}};
fesa::Vector high(24), positive_low(24), negative_low(24);
for (std::size_t a = 0; a < 8; ++a) {
high[3 * a] = -coordinates[a][0];
positive_low[3 * a] = std::ldexp(coordinates[a][0], -60);
negative_low[3 * a] = -positive_low[3 * a];
}
using fesa::hex8_internal::Certificate;
EXPECT_EQ(fesa::hex8_internal::Certify(coordinates, &high, &positive_low),
Certificate::kPositive);
EXPECT_EQ(fesa::hex8_internal::Certify(coordinates, &high, &negative_low),
Certificate::kNonpositive);
EXPECT_NE(fesa::hex8_internal::Certify(coordinates, &high),
Certificate::kPositive);
}
TEST(Hex8HighLow, SignedZeroMaterialCoefficientHasIdenticalResponse) {
auto positive = UnitBrick(0, 0.0), negative = UnitBrick(0, -0.0);
fesa::Vector high(24), low(24);
high[3] = 0.01;
for (bool nonlinear : {false, true}) {
const auto p = nonlinear ? positive.EvaluateNonlinearHighLow(high, low)
: positive.EvaluateLinearHighLow(high, low);
const auto n = nonlinear ? negative.EvaluateNonlinearHighLow(high, low)
: negative.EvaluateLinearHighLow(high, low);
ASSERT_TRUE(p.HasValue());
ASSERT_TRUE(n.HasValue());
EXPECT_EQ(p.Value().high.strain_energy, n.Value().high.strain_energy);
EXPECT_EQ(p.Value().strain_energy_low, n.Value().strain_energy_low);
for (std::size_t i = 0; i < 24; ++i) {
EXPECT_EQ(p.Value().high.internal_force[i], n.Value().high.internal_force[i]);
EXPECT_EQ(p.Value().internal_force_low[i], n.Value().internal_force_low[i]);
for (std::size_t j = 0; j < 24; ++j) {
EXPECT_EQ(p.Value().high.tangent(i,j), n.Value().high.tangent(i,j));
EXPECT_EQ(p.Value().tangent_low(i,j), n.Value().tangent_low(i,j));
}
}
}
}
TEST(Hex8HighLow, FiniteRotationRetainsTinyShearForceCovariance) {
using namespace fesa::math_internal;
auto brick = UnitBrick();
fesa::Hex8Element::Coordinates x{{{0, 0, 0}, {1, 0, 0},
{1, 1, 0}, {0, 1, 0},
{0, 0, 1}, {1, 0, 1},
{1, 1, 1}, {0, 1, 1}}};
fesa::Vector zero(24), rotation(24), shear(24), rotated_shear(24);
const double gamma = std::ldexp(1.0, -60);
for (std::size_t a = 0; a < 8; ++a) {
// Q is an exact quarter turn; Q*(gamma*y,0,0)=(0,gamma*y,0).
rotation[3*a] = HighLow{-x[a][1] - x[a][0]}.high;
rotation[3*a+1] = x[a][0] - x[a][1];
shear[3*a] = gamma*x[a][1];
rotated_shear[3*a+1] = gamma*x[a][1];
}
fesa::Vector turned_high(24), turned_low(24);
for (std::size_t j = 0; j < 24; ++j) {
const auto full = Add(HighLow{rotation[j]}, HighLow{rotated_shear[j]});
turned_high[j] = full.high; turned_low[j] = full.low;
}
const auto origin = brick.EvaluateNonlinearHighLow(zero, zero);
const auto deformed = brick.EvaluateNonlinearHighLow(shear, zero);
const auto turned = brick.EvaluateNonlinearHighLow(rotation, zero);
const auto turned_deformed = brick.EvaluateNonlinearHighLow(turned_high, turned_low);
ASSERT_TRUE(origin.HasValue()); ASSERT_TRUE(deformed.HasValue());
ASSERT_TRUE(turned.HasValue()); ASSERT_TRUE(turned_deformed.HasValue());
// The fixed F64 reference gradients leave coefficient-level rotation error.
// Subtract the same rigid baseline; do not project that force to zero.
for (std::size_t a = 0; a < 8; ++a) {
std::array<HighLow, 3> baseline{}, rotated{};
for (std::size_t i = 0; i < 3; ++i) {
const auto j = 3*a+i;
baseline[i] = Subtract(
{deformed.Value().high.internal_force[j], deformed.Value().internal_force_low[j]},
{origin.Value().high.internal_force[j], origin.Value().internal_force_low[j]});
rotated[i] = Subtract(
{turned_deformed.Value().high.internal_force[j], turned_deformed.Value().internal_force_low[j]},
{turned.Value().high.internal_force[j], turned.Value().internal_force_low[j]});
}
EXPECT_LE(std::abs(Add(rotated[0], baseline[1]).high), 1e-13*gamma);
EXPECT_LE(std::abs(Subtract(rotated[1], baseline[0]).high), 1e-13*gamma);
EXPECT_LE(std::abs(Subtract(rotated[2], baseline[2]).high), 1e-13*gamma);
}
}
TEST(Hex8HighLow,
MeanVolumeLowOnlyNonuniformTraceRetainsMeanAndCrossPointHessian) {
auto brick = UnitBrick(0, .25);
fesa::Vector high(24), low(24), zero(24);
const double c = std::ldexp(1.0, -60);
for (std::size_t n = 0; n < 8; ++n) high[3 * n] = 1;
low[6] = c;
low[18] = c;
const double coefficient = 7 * 4.8 / 18 + 8.0 / 8;
const auto response = brick.EvaluateLinearHighLow(high, low);
const auto high_only = brick.EvaluateLinearHighLow(high, zero);
ASSERT_TRUE(response.HasValue());
ASSERT_TRUE(high_only.HasValue());
EXPECT_EQ(high_only.Value().high.strain_energy, 0);
EXPECT_NEAR(response.Value().high.strain_energy, coefficient * c * c,
1e-13 * c * c);
EXPECT_NEAR(response.Value().high.internal_force[6] +
response.Value().high.internal_force[18],
2 * coefficient * c, 1e-13 * c);
for (const auto& point : response.Value().high.points) {
const auto& s = point.cauchy_stress;
EXPECT_NEAR((s[0] + s[1] + s[2]) / 3, 4 * c, 1e-13 * c);
EXPECT_NEAR(point.strain[0], c * point.reference_coordinates[1], 1e-13 * c);
}
const auto& k = response.Value().high.tangent;
EXPECT_NEAR(k(0, 0), 41.0 / 18, 2e-14);
EXPECT_NEAR(k(0, 1), 19.0 / 30, 2e-14);
const auto nonlinear = brick.EvaluateNonlinearHighLow(high, low);
ASSERT_TRUE(nonlinear.HasValue());
const auto cubic_energy = Subtract({nonlinear.Value().high.strain_energy,
nonlinear.Value().strain_energy_low},
{response.Value().high.strain_energy,
response.Value().strain_energy_low});
EXPECT_NEAR(cubic_energy.high, (5 * 4.8 / 18 + 8.0 / 6) * c * c * c,
1e-12 * c * c * c);
}
} // namespace
Reference in New Issue View Git Blame Copy Permalink