From eabaf6c89399d7187c14e2eb60955e66631f70c1 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Wed, 14 Oct 2026 16:56:11 +0900 Subject: [PATCH] perf(hex8): use MKL nodal mass and preserve rejected diagnostic priority --- src/fesa/elements/hex8_element.cpp | 105 +++++++++++++++++++++++------ 1 file changed, 84 insertions(+), 21 deletions(-) diff --git a/src/fesa/elements/hex8_element.cpp b/src/fesa/elements/hex8_element.cpp index 9e758bc..cc3eb6a 100644 --- a/src/fesa/elements/hex8_element.cpp +++ b/src/fesa/elements/hex8_element.cpp @@ -180,6 +180,7 @@ Result Hex8Element::Create(EntityIndex index, // Zero Hessian of the same reference-relative mean-volume energy. std::array integrated_gradients{}; double volume = 0.0; + bool backend_inputs_failed = false; for (const auto& point : candidate.points_) { auto gradients = point.gradients; gradients[0] = {}; @@ -192,22 +193,34 @@ Result Hex8Element::Create(EntityIndex index, Matrix deviatoric{0U, 0U}; Matrix shear{0U, 0U}; Matrix gram{0U, 0U}; + bool point_products_failed = false; 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.")); + backend_inputs_failed = true; + point_products_failed = true; } 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) += - 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 (point_products_failed) { + // Diagnostic replay only: this candidate can never succeed. + const auto& ga = gradients[a]; + const auto& gb = gradients[b]; + 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(ga, gb) : 0.0)); + } else { + candidate.stiffness_(3U * a + i, 3U * b + k) += + 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)); + } } } } @@ -218,8 +231,15 @@ Result Hex8Element::Create(EntityIndex index, try { bulk = Outer(integrated_gradients, candidate.kappa_ / volume); } catch (const std::invalid_argument&) { - return Result::Failure( - candidate.Failure("Hex8 reference stiffness is nonfinite.")); + backend_inputs_failed = true; + bulk = Matrix{24U, 24U}; + 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) + bulk(3U * a + i, 3U * b + k) = + (candidate.kappa_ / volume) * integrated_gradients[a][i] * + integrated_gradients[b][k]; } for (std::size_t i = 0; i < 24; ++i) for (std::size_t j = 0; j < 24; ++j) @@ -234,6 +254,9 @@ Result Hex8Element::Create(EntityIndex index, return Result::Failure( candidate.Failure("Hex8 reference stiffness has an unrepresentable " "positive diagonal.")); + if (backend_inputs_failed) + return Result::Failure(candidate.Failure( + "Hex8 dense MKL operands failed finite validation.")); const Vector zero{24U}; const auto paired = candidate.EvaluateHighLow(zero, zero, false); if (paired.HasValue()) { @@ -262,17 +285,30 @@ Result Hex8Element::ComputeMass() const { return Result::Failure( Failure("Hex8 mass requires reference density.")); Matrix consistent{24U, 24U}; - for (const auto& point : points_) + for (const auto& point : points_) { + Matrix left{8U, 1U}; + Matrix right{1U, 8U}; + for (std::size_t a = 0; a < 8; ++a) { + left(a, 0) = *density_ * point.shape[a]; + right(0, a) = point.shape[a]; + } + Matrix products{0U, 0U}; + try { + products = left.Multiply(right); + } catch (const std::invalid_argument&) { + return Result::Failure( + Failure("Hex8 reference mass is nonfinite.")); + } for (std::size_t a = 0; a < 8; ++a) for (std::size_t b = 0; b < 8; ++b) { - const double value = - *density_ * point.shape[a] * point.shape[b] * point.weight; + const double value = products(a, b) * point.weight; if (!std::isfinite(value)) return Result::Failure( Failure("Hex8 reference mass is nonfinite.")); for (std::size_t i = 0; i < 3; ++i) consistent(3U * a + i, 3U * b + i) += value; } + } Matrix mass{24U, 24U}; for (std::size_t i = 0; i < 24; ++i) for (std::size_t j = 0; j < 24; ++j) { @@ -387,6 +423,7 @@ Result Hex8Element::Evaluate(const Vector& u, Failure("Hex8 mean volume is nonfinite.")); response.strain_energy = (0.5 * kappa_ * mean_trace) * (volume * mean_trace); std::array integrated_fg{}; + bool backend_inputs_failed = false; for (std::size_t g = 0; g < 8; ++g) { const auto& point = points_[g]; const auto& f = deformations[g]; @@ -428,6 +465,7 @@ Result Hex8Element::Evaluate(const Vector& u, Matrix shear{0U, 0U}; Matrix gram{0U, 0U}; Matrix geometric{0U, 0U}; + bool point_products_failed = false; if (!Finite(f) || !Finite(strain) || !Finite(s) || !Finite(p) || !Finite(fft) || !Finite(stress)) return Result::Failure( @@ -443,8 +481,8 @@ Result Hex8Element::Evaluate(const Vector& u, gram = GradientProduct(gradients[g]); geometric = GradientProduct(gradients[g], &s); } catch (const std::invalid_argument&) { - return Result::Failure( - Failure("Hex8 trial tangent is nonfinite.")); + backend_inputs_failed = true; + point_products_failed = true; } } for (std::size_t a = 0; a < 8; ++a) { @@ -464,18 +502,34 @@ Result Hex8Element::Evaluate(const Vector& u, } if (nonlinear) for (std::size_t b = 0; b < 8; ++b) { + const auto& gb = gradients[g][b]; + double diagnostic_geometric = 0.0; + if (point_products_failed) + for (std::size_t j = 0; j < 3; ++j) + for (std::size_t l = 0; l < 3; ++l) + diagnostic_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 * - (deviatoric(row, column) + - shear(3U * b + i, 3U * a + k) + - mu_ * fft[i][k] * gram(a, b)); + if (point_products_failed) { + // Reproduce original diagnostic priority, never acceptance. + 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)); + } else { + response.material_tangent(row, column) += + point.weight * + (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(a, b) * point.weight; + (point_products_failed ? diagnostic_geometric + : geometric(a, b)) * point.weight; } } } @@ -490,8 +544,14 @@ Result Hex8Element::Evaluate(const Vector& u, try { bulk = Outer(integrated_fg, kappa_ / volume); } catch (const std::invalid_argument&) { - return Result::Failure( - Failure("Hex8 trial tangent is nonfinite.")); + backend_inputs_failed = true; + bulk = Matrix{24U, 24U}; + 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) + bulk(3U * a + i, 3U * b + k) = + (kappa_ / volume) * integrated_fg[a][i] * integrated_fg[b][k]; } for (std::size_t i = 0; i < 24; ++i) for (std::size_t j = 0; j < 24; ++j) @@ -515,6 +575,9 @@ Result Hex8Element::Evaluate(const Vector& u, Failure("Hex8 trial tangent is nonfinite.")); } } + if (backend_inputs_failed) + return Result::Failure( + Failure("Hex8 dense MKL operands failed finite validation.")); return Result::Success(std::move(response)); } } // namespace fesa