From e18271051c699e2f80a7c4ca9055e52077781140 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Wed, 14 Oct 2026 16:49:15 +0900 Subject: [PATCH] perf(hex8): use MKL Hessian products and exact zero pair addition --- src/fesa/elements/hex8_element.cpp | 119 ++++++++++++++++----- src/fesa/elements/hex8_high_low.cpp | 4 +- src/fesa/elements/hex8_pair_arithmetic.h | 20 +++- tests/unit/elements/hex8_high_low_test.cpp | 7 ++ 4 files changed, 116 insertions(+), 34 deletions(-) diff --git a/src/fesa/elements/hex8_element.cpp b/src/fesa/elements/hex8_element.cpp index b715f5d..e477244 100644 --- a/src/fesa/elements/hex8_element.cpp +++ b/src/fesa/elements/hex8_element.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include "hex8_geometry.h" @@ -26,6 +27,36 @@ std::array Pack(const Tensor& value) { double Dot(const Triple& a, const Triple& b) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; } +/// @brief Forms a scaled 24-component outer product through dense MKL. +Matrix Outer(const std::array& values, double coefficient) { + Matrix left{24U, 1U}; + Matrix right{1U, 24U}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t i = 0; i < 3; ++i) { + left(3U * a + i, 0) = coefficient * values[a][i]; + right(0, 3U * a + i) = values[a][i]; + } + return left.Multiply(right); +} +/// @brief Forms node-gradient Gram and stress congruence through dense MKL. +Matrix GradientProduct(const std::array& gradients, + const Tensor* stress = nullptr) { + Matrix left{8U, 3U}; + Matrix transpose{3U, 8U}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t i = 0; i < 3; ++i) { + left(a, i) = gradients[a][i]; + transpose(i, a) = gradients[a][i]; + } + if (stress != nullptr) { + Matrix stress_transpose{3U, 3U}; + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) + stress_transpose(i, j) = (*stress)[j][i]; + left = left.Multiply(stress_transpose); + } + return left.Multiply(transpose); +} /// @brief Retains product and cancellation errors in fixed-order contractions. void AccumulateProduct(double left, double right, double& sum, double& correction) { @@ -158,30 +189,41 @@ Result Hex8Element::Create(EntityIndex index, for (std::size_t a = 0; a < 8; ++a) for (std::size_t i = 0; i < 3; ++i) integrated_gradients[a][i] += point.weight * gradients[a][i]; + Matrix deviatoric{0U, 0U}; + Matrix shear{0U, 0U}; + Matrix gram{0U, 0U}; + try { + deviatoric = Outer(gradients, -(2.0 / 3.0) * candidate.mu_); + shear = Outer(gradients, candidate.mu_); + gram = GradientProduct(gradients); + } catch (const std::invalid_argument&) { + return Result::Failure( + candidate.Failure("Hex8 reference stiffness is nonfinite.")); + } for (std::size_t a = 0; a < 8; ++a) for (std::size_t b = 0; b < 8; ++b) { - const auto& ga = gradients[a]; - const auto& gb = gradients[b]; - const double dot = Dot(ga, gb); for (std::size_t i = 0; i < 3; ++i) for (std::size_t k = 0; k < 3; ++k) { candidate.stiffness_(3U * a + i, 3U * b + k) += - point.weight * (-(2.0 / 3.0) * candidate.mu_ * ga[i] * gb[k] + - candidate.mu_ * ga[k] * gb[i] + - (i == k ? candidate.mu_ * dot : 0.0)); + point.weight * (deviatoric(3U * a + i, 3U * b + k) + + shear(3U * a + k, 3U * b + i) + + (i == k ? candidate.mu_ * gram(a, b) : 0.0)); } } } if (!std::isfinite(volume) || volume <= 0.0) return Result::Failure( candidate.Failure("Hex8 reference volume is nonfinite.")); - for (std::size_t a = 0; a < 8; ++a) - for (std::size_t b = 0; b < 8; ++b) - for (std::size_t i = 0; i < 3; ++i) - for (std::size_t k = 0; k < 3; ++k) - candidate.stiffness_(3U * a + i, 3U * b + k) += - (candidate.kappa_ / volume) * integrated_gradients[a][i] * - integrated_gradients[b][k]; + Matrix bulk{0U, 0U}; + try { + bulk = Outer(integrated_gradients, candidate.kappa_ / volume); + } catch (const std::invalid_argument&) { + return Result::Failure( + candidate.Failure("Hex8 reference stiffness is nonfinite.")); + } + for (std::size_t i = 0; i < 24; ++i) + for (std::size_t j = 0; j < 24; ++j) + candidate.stiffness_(i, j) += bulk(i, j); for (std::size_t i = 0; i < 24; ++i) for (std::size_t j = 0; j < 24; ++j) if (!std::isfinite(candidate.stiffness_(i, j))) @@ -382,6 +424,25 @@ Result Hex8Element::Evaluate(const Vector& u, for (std::size_t i = 0; i < 3; ++i) integrated_fg[a][i] += point.weight * Dot(f[i], gradients[g][a]); } + Matrix deviatoric{0U, 0U}; + Matrix shear{0U, 0U}; + Matrix gram{0U, 0U}; + Matrix geometric{0U, 0U}; + if (nonlinear) { + std::array fg{}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t i = 0; i < 3; ++i) + fg[a][i] = Dot(f[i], gradients[g][a]); + try { + deviatoric = Outer(fg, -(2.0 / 3.0) * mu_); + shear = Outer(fg, mu_); + gram = GradientProduct(gradients[g]); + geometric = GradientProduct(gradients[g], &s); + } catch (const std::invalid_argument&) { + return Result::Failure( + Failure("Hex8 trial tangent is nonfinite.")); + } + } for (std::size_t a = 0; a < 8; ++a) { const auto& ga = gradients[g][a]; for (std::size_t i = 0; i < 3; ++i) { @@ -399,23 +460,18 @@ Result Hex8Element::Evaluate(const Vector& u, } if (nonlinear) for (std::size_t b = 0; b < 8; ++b) { - const auto& gb = gradients[g][b]; - double geometric = 0.0; - for (std::size_t j = 0; j < 3; ++j) - for (std::size_t l = 0; l < 3; ++l) - geometric += ga[j] * s[l][j] * gb[l]; for (std::size_t i = 0; i < 3; ++i) for (std::size_t k = 0; k < 3; ++k) { const auto row = 3U * a + i; const auto column = 3U * b + k; response.material_tangent(row, column) += point.weight * - (-(2.0 / 3.0) * mu_ * Dot(f[i], ga) * Dot(f[k], gb) + - mu_ * Dot(f[i], gb) * Dot(f[k], ga) + - mu_ * fft[i][k] * Dot(ga, gb)); + (deviatoric(row, column) + + shear(3U * b + i, 3U * a + k) + + mu_ * fft[i][k] * gram(a, b)); if (i == k) response.geometric_tangent(row, column) += - geometric * point.weight; + geometric(a, b) * point.weight; } } } @@ -429,13 +485,18 @@ Result Hex8Element::Evaluate(const Vector& u, } for (std::size_t i = 0; i < 24; ++i) response.internal_force[i] += force_correction[i]; - if (nonlinear) - for (std::size_t a = 0; a < 8; ++a) - for (std::size_t b = 0; b < 8; ++b) - for (std::size_t i = 0; i < 3; ++i) - for (std::size_t k = 0; k < 3; ++k) - response.material_tangent(3 * a + i, 3 * b + k) += - (kappa_ / volume) * integrated_fg[a][i] * integrated_fg[b][k]; + if (nonlinear) { + Matrix bulk{0U, 0U}; + try { + bulk = Outer(integrated_fg, kappa_ / volume); + } catch (const std::invalid_argument&) { + return Result::Failure( + Failure("Hex8 trial tangent is nonfinite.")); + } + for (std::size_t i = 0; i < 24; ++i) + for (std::size_t j = 0; j < 24; ++j) + response.material_tangent(i, j) += bulk(i, j); + } if (!std::isfinite(response.strain_energy)) return Result::Failure( Failure("Hex8 trial energy is nonfinite.")); diff --git a/src/fesa/elements/hex8_high_low.cpp b/src/fesa/elements/hex8_high_low.cpp index bba1016..6d5842b 100644 --- a/src/fesa/elements/hex8_high_low.cpp +++ b/src/fesa/elements/hex8_high_low.cpp @@ -14,7 +14,9 @@ using math_internal::HighLow; using Triple = std::array; using Tensor = std::array; -HighLow operator+(HighLow a, HighLow b) { return math_internal::Add(a, b); } +HighLow operator+(HighLow a, HighLow b) { + return hex8_internal::AddPair(a, b); +} HighLow operator-(HighLow a, HighLow b) { return math_internal::Subtract(a, b); } diff --git a/src/fesa/elements/hex8_pair_arithmetic.h b/src/fesa/elements/hex8_pair_arithmetic.h index d640288..2b8fb58 100644 --- a/src/fesa/elements/hex8_pair_arithmetic.h +++ b/src/fesa/elements/hex8_pair_arithmetic.h @@ -4,6 +4,18 @@ #include "../math/high_low_arithmetic.h" namespace fesa::hex8_internal { +/// @brief Preserves original zero-add normalization without redundant TwoSums. +/// @note Valid nonzero terms retain the original addition association. +inline math_internal::HighLow AddPair(math_internal::HighLow left, + math_internal::HighLow right) { + using namespace math_internal; + if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow(); + if (left.high == 0 && left.low == 0) + return {right.high, right.low == 0 ? 0.0 : right.low}; + if (right.high == 0 && right.low == 0) + return {left.high, left.low == 0 ? 0.0 : left.low}; + return Add(left, right); +} /// @brief Omits exact-zero cross products without changing nonzero association. /// @note A retained zero addition performs the original canonicalization and /// invalid-result check; its repeated application is idempotent. @@ -14,14 +26,14 @@ inline math_internal::HighLow MultiplyPair(math_internal::HighLow left, if (left.low != 0 && right.low != 0) return Multiply(left, right); auto product = TwoProduct(left.high, right.high); if (right.low != 0) { - product = Add(product, TwoProduct(left.high, right.low)); + product = AddPair(product, TwoProduct(left.high, right.low)); } else if (left.low != 0) { // The original zero addition precedes the nonzero low-high term. - product = Add(product, HighLow{}); + product = AddPair(product, HighLow{}); } if (left.low != 0) - product = Add(product, TwoProduct(left.low, right.high)); - return Add(product, HighLow{}); + product = AddPair(product, TwoProduct(left.low, right.high)); + return AddPair(product, HighLow{}); } } // namespace fesa::hex8_internal diff --git a/tests/unit/elements/hex8_high_low_test.cpp b/tests/unit/elements/hex8_high_low_test.cpp index a589a98..6ae57cc 100644 --- a/tests/unit/elements/hex8_high_low_test.cpp +++ b/tests/unit/elements/hex8_high_low_test.cpp @@ -31,6 +31,11 @@ void ExpectOriginalProduct(HighLow left, HighLow right) { {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) { @@ -65,6 +70,8 @@ TEST(Hex8HighLow, ScalarCrossProductShortcutMatchesOriginalBits) { 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);