From 18b36c1c228c535558d15923a5e0fcc0b55be486 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Sat, 10 Oct 2026 23:32:40 +0900 Subject: [PATCH] feat(hex8): use reference-weighted mean-volume energy --- include/fesa/elements/hex8_element.h | 10 +- src/fesa/elements/hex8_element.cpp | 199 ++++++++------- src/fesa/elements/hex8_high_low.cpp | 66 ++++- .../solid_frequency_procedure_test.cpp | 7 +- .../solid_hht_dynamic_procedure_test.cpp | 13 +- .../solid_hht_temporal_verification_test.cpp | 121 ++++++--- tests/unit/elements/hex8_element_test.cpp | 240 +++++++++++++++++- tests/unit/elements/hex8_high_low_test.cpp | 36 +++ 8 files changed, 537 insertions(+), 155 deletions(-) diff --git a/include/fesa/elements/hex8_element.h b/include/fesa/elements/hex8_element.h index a2fd838..9bd0f98 100644 --- a/include/fesa/elements/hex8_element.h +++ b/include/fesa/elements/hex8_element.h @@ -30,7 +30,7 @@ struct Hex8HighLowResponse { double strain_energy_low{0.0}; }; -/// @brief Implements full-integration C3D8 linear and TL StVK mechanics. +/// @brief Implements C3D8 full deviatoric and mean-volume Green energy. /// @note Reference geometry is owned and immutable; no trial is committed. class Hex8Element final : public Element { public: @@ -53,14 +53,14 @@ class Hex8Element final : public Element { /// @brief Returns f=K0*u and infinitesimal tensor strain at eight points. Result EvaluateLinear(const Vector& displacement) const; - /// @brief Returns one certified StVK trial and its analytic energy Hessian. + /// @brief Returns one certified mean-Green TL trial and its energy Hessian. Result EvaluateNonlinear(const Vector& displacement) const; /// @brief Evaluates linear energy/force/Hessian on the complete paired state. /// @note Uses the same reference-relative gradient operator as the paired TL /// path. Inputs must have 24 finite normalized high+low translations. Result EvaluateLinearHighLow( const Vector& high, const Vector& low) const; - /// @brief Evaluates certified TL StVK mechanics without collapsing low state. + /// @brief Evaluates mean-Green TL mechanics without collapsing low state. /// @note Failure preserves immutable reference data and all caller state. Result EvaluateNonlinearHighLow( const Vector& high, const Vector& low) const; @@ -81,7 +81,7 @@ class Hex8Element final : public Element { }; /// @brief Constructs an owned candidate before reference validation. Hex8Element(EntityIndex element_index, ElementDofLayout layout, - Coordinates coordinates, double lambda, double mu, + Coordinates coordinates, double kappa, double mu, std::optional density, SourceLocation location); /// @brief Evaluates without mutating the cached reference or previous values. Result Evaluate(const Vector& displacement, @@ -95,7 +95,7 @@ class Hex8Element final : public Element { EntityIndex element_index_; ElementDofLayout layout_; Coordinates coordinates_; - double lambda_; + double kappa_; double mu_; std::optional density_; SourceLocation location_; diff --git a/src/fesa/elements/hex8_element.cpp b/src/fesa/elements/hex8_element.cpp index e7f711f..8e282e3 100644 --- a/src/fesa/elements/hex8_element.cpp +++ b/src/fesa/elements/hex8_element.cpp @@ -61,12 +61,12 @@ const char* CertificateMessage(hex8_internal::Certificate certificate) { } // namespace Hex8Element::Hex8Element(EntityIndex index, ElementDofLayout layout, - Coordinates coordinates, double lambda, double mu, + Coordinates coordinates, double kappa, double mu, std::optional density, SourceLocation location) : element_index_{index}, layout_{std::move(layout)}, coordinates_{coordinates}, - lambda_{lambda}, + kappa_{kappa}, mu_{mu}, density_{density}, location_{std::move(location)} {} @@ -102,23 +102,12 @@ Result Hex8Element::Create(EntityIndex index, return Result::Failure( candidate.Failure("Hex8 requires finite E>0 and -1::Failure( - candidate.Failure("Hex8 Lame moduli must remain finite with " - "representable positive shear modulus.")); - // Validate the rounded moduli actually used by the kernel, not merely the - // positive exact-real E/nu formula: cancellation near nu=-1 can erase bulk. - const double bulk = candidate.lambda_ + (2.0 * candidate.mu_) / 3.0; - const double volumetric_stress = - 3.0 * candidate.lambda_ + 2.0 * candidate.mu_; - if (!std::isfinite(bulk) || bulk <= 0.0 || - !std::isfinite(volumetric_stress) || volumetric_stress <= 0.0) - return Result::Failure( - candidate.Failure("Hex8 rounded Lame moduli require finite " - "representable positive bulk stiffness.")); + candidate.Failure("Hex8 shear and direct bulk moduli must remain " + "finite and representably positive.")); if (density && (!std::isfinite(*density) || *density <= 0.0)) return Result::Failure(candidate.Failure( "Hex8 reference density must be finite and positive.")); @@ -156,20 +145,43 @@ Result Hex8Element::Create(EntityIndex index, return Result::Failure(candidate.Failure( "Reference gradient or point coordinate is nonfinite.")); } + } + // Zero Hessian of the same reference-relative mean-volume energy. + std::array integrated_gradients{}; + double volume = 0.0; + for (const auto& point : candidate.points_) { + auto gradients = point.gradients; + gradients[0] = {}; + for (std::size_t a = 1; a < 8; ++a) + for (std::size_t i = 0; i < 3; ++i) gradients[0][i] -= gradients[a][i]; + volume += point.weight; + 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]; for (std::size_t a = 0; a < 8; ++a) for (std::size_t b = 0; b < 8; ++b) { - const auto& ga = point.gradients[a]; - const auto& gb = point.gradients[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 * (candidate.lambda_ * ga[i] * gb[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)); } } } + 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]; 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))) @@ -278,48 +290,33 @@ Result Hex8Element::Evaluate(const Vector& u, if (!nonlinear) { response.material_tangent = stiffness_; response.tangent = stiffness_; - // K0 annihilates uniform translation. Use the same relative field as - // strain recovery to avoid cancellation of a large rigid displacement. - // Retain low product bits and cancelled sum terms in bending-dominated - // rows; the stored stiffness and its tangent are unchanged. - for (std::size_t i = 0; i < 24; ++i) { - double sum = 0.0; - double correction = 0.0; - for (std::size_t j = 0; j < 24; ++j) { - AccumulateProduct(stiffness_(i, j), relative[j], sum, correction); - AccumulateProduct(stiffness_(i, j), relative_correction[j], sum, - correction); - } - response.internal_force[i] = sum + correction; - } - double energy = 0.0; - double energy_correction = 0.0; - for (std::size_t i = 0; i < 24; ++i) { - AccumulateProduct(relative[i], response.internal_force[i], energy, - energy_correction); - AccumulateProduct(relative_correction[i], response.internal_force[i], - energy, energy_correction); - } - response.strain_energy = 0.5 * (energy + energy_correction); } + std::array deformations{}, strains{}; + std::array, 8> gradients{}; + double volume = 0.0, trace_integral = 0.0; Vector force_correction{24U}; for (std::size_t g = 0; g < 8; ++g) { const auto& point = points_[g]; + gradients[g] = point.gradients; + Triple sum{}; + for (std::size_t a = 1; a < 8; ++a) + for (std::size_t j = 0; j < 3; ++j) sum[j] += gradients[g][a][j]; + for (std::size_t j = 0; j < 3; ++j) gradients[g][0][j] = -sum[j]; Tensor gradient{}; Tensor gradient_correction{}; - for (std::size_t a = 0; a < 8; ++a) + for (std::size_t a = 1; a < 8; ++a) for (std::size_t i = 0; i < 3; ++i) for (std::size_t j = 0; j < 3; ++j) { AccumulateProduct(relative[3U * a + i], point.gradients[a][j], gradient[i][j], gradient_correction[i][j]); AccumulateProduct(relative_correction[3U * a + i], - point.gradients[a][j], - gradient[i][j], gradient_correction[i][j]); + point.gradients[a][j], gradient[i][j], + gradient_correction[i][j]); } for (std::size_t i = 0; i < 3; ++i) for (std::size_t j = 0; j < 3; ++j) gradient[i][j] += gradient_correction[i][j]; - Tensor f = gradient; + Tensor f = nonlinear ? gradient : Tensor{}; for (std::size_t i = 0; i < 3; ++i) f[i][i] += 1.0; Tensor strain{}; for (std::size_t i = 0; i < 3; ++i) @@ -330,51 +327,76 @@ Result Hex8Element::Evaluate(const Vector& u, for (std::size_t k = 0; k < 3; ++k) strain[i][j] += 0.5 * gradient[k][i] * gradient[k][j]; } + deformations[g] = f; + strains[g] = strain; + volume += point.weight; + trace_integral += + point.weight * (strain[0][0] + strain[1][1] + strain[2][2]); + if (!Finite(gradient) || !Finite(f) || !Finite(strain)) + return Result::Failure( + Failure("Hex8 trial tensor is nonfinite.")); + } + const double mean_trace = trace_integral / volume; + if (!std::isfinite(volume) || volume <= 0 || !std::isfinite(mean_trace)) + return Result::Failure( + Failure("Hex8 mean volume is nonfinite.")); + response.strain_energy = (0.5 * kappa_ * mean_trace) * (volume * mean_trace); + std::array integrated_fg{}; + for (std::size_t g = 0; g < 8; ++g) { + const auto& point = points_[g]; + const auto& f = deformations[g]; + const auto& strain = strains[g]; const double trace = strain[0][0] + strain[1][1] + strain[2][2]; - Tensor s{}; + Tensor s{}, p{}, fft{}; + double density_energy = 0.0; for (std::size_t i = 0; i < 3; ++i) - for (std::size_t j = 0; j < 3; ++j) - s[i][j] = 2.0 * mu_ * strain[i][j] + (i == j ? lambda_ * trace : 0.0); + for (std::size_t j = 0; j < 3; ++j) { + const double dev = strain[i][j] - (i == j ? trace / 3.0 : 0.0); + s[i][j] = 2.0 * mu_ * dev + (i == j ? kappa_ * mean_trace : 0.0); + density_energy += mu_ * dev * dev; + } + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) { + for (std::size_t k = 0; k < 3; ++k) { + p[i][j] += f[i][k] * s[k][j]; + fft[i][j] += f[i][k] * f[j][k]; + } + } + response.strain_energy += density_energy * point.weight; Tensor stress = s; if (nonlinear) { const double det_f = hex8_internal::Determinant(f); if (!std::isfinite(det_f) || det_f <= 0.0) return Result::Failure(Failure( "Hex8 trial deformation determinant must be finite and positive.")); - Tensor p{}; - Tensor fft{}; - double energy = 0.5 * lambda_ * trace * trace; - for (std::size_t i = 0; i < 3; ++i) - for (std::size_t j = 0; j < 3; ++j) { - energy += mu_ * strain[i][j] * strain[i][j]; - for (std::size_t k = 0; k < 3; ++k) { - p[i][j] += f[i][k] * s[k][j]; - fft[i][j] += f[i][k] * f[j][k]; - } - } - response.strain_energy += energy * point.weight; stress = {}; for (std::size_t i = 0; i < 3; ++i) for (std::size_t j = 0; j < 3; ++j) for (std::size_t k = 0; k < 3; ++k) stress[i][j] += p[i][k] * f[j][k] / det_f; - for (std::size_t a = 0; a < 8; ++a) { - const auto& ga = point.gradients[a]; - for (std::size_t i = 0; i < 3; ++i) { - double contraction = 0.0; - double contraction_correction = 0.0; - for (std::size_t j = 0; j < 3; ++j) - AccumulateProduct(p[i][j], ga[j], contraction, - contraction_correction); - AccumulateProduct(contraction, point.weight, - response.internal_force[3U * a + i], - force_correction[3U * a + i]); - AccumulateProduct(contraction_correction, point.weight, - response.internal_force[3U * a + i], - force_correction[3U * a + i]); - } + } + for (std::size_t a = 0; a < 8; ++a) { + for (std::size_t i = 0; i < 3; ++i) + integrated_fg[a][i] += point.weight * Dot(f[i], gradients[g][a]); + } + for (std::size_t a = 0; a < 8; ++a) { + const auto& ga = gradients[g][a]; + for (std::size_t i = 0; i < 3; ++i) { + double contraction = 0.0; + double contraction_correction = 0.0; + for (std::size_t j = 0; j < 3; ++j) + AccumulateProduct(p[i][j], ga[j], contraction, + contraction_correction); + AccumulateProduct(contraction, point.weight, + response.internal_force[3U * a + i], + force_correction[3U * a + i]); + AccumulateProduct(contraction_correction, point.weight, + response.internal_force[3U * a + i], + force_correction[3U * a + i]); + } + if (nonlinear) for (std::size_t b = 0; b < 8; ++b) { - const auto& gb = point.gradients[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) @@ -384,26 +406,33 @@ Result Hex8Element::Evaluate(const Vector& u, const auto row = 3U * a + i; const auto column = 3U * b + k; response.material_tangent(row, column) += - point.weight * (lambda_ * Dot(f[i], ga) * Dot(f[k], gb) + - mu_ * Dot(f[i], gb) * Dot(f[k], ga) + - mu_ * fft[i][k] * Dot(ga, gb)); + 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)); if (i == k) response.geometric_tangent(row, column) += geometric * point.weight; } } - } } - if (!Finite(f) || !Finite(strain) || !Finite(s) || !Finite(stress)) + if (!Finite(f) || !Finite(strain) || !Finite(s) || !Finite(p) || + !Finite(fft) || !Finite(stress)) return Result::Failure( Failure("Hex8 trial tensor is nonfinite.")); response.points[g] = {element_index_, static_cast(g + 1U), point.natural, point.position, Pack(stress), Pack(strain)}; } + for (std::size_t i = 0; i < 24; ++i) + response.internal_force[i] += force_correction[i]; if (nonlinear) - for (std::size_t i = 0; i < 24; ++i) - response.internal_force[i] += force_correction[i]; + 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 (!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 c6577be..e18fe8c 100644 --- a/src/fesa/elements/hex8_high_low.cpp +++ b/src/fesa/elements/hex8_high_low.cpp @@ -86,7 +86,10 @@ Result Hex8Element::EvaluateHighLow( Hex8HighLowResponse response; std::array force{}; std::array, 24> material{}, geometric{}; - HighLow energy; + HighLow energy, volume, trace_integral; + std::array deformations{}, strains{}; + std::array, 8> point_gradients{}; + std::array integrated_fg{}; for (std::size_t g = 0; g < 8; ++g) { const auto& point = points_[g]; std::array gradients{}; @@ -104,7 +107,7 @@ Result Hex8Element::EvaluateHighLow( for (std::size_t j = 0; j < 3; ++j) gradients[0][j] = math_internal::Negate(gradient_sum[j]); Tensor f{}; - Tensor strain{}, s{}, p{}, fft{}; + Tensor strain{}; for (std::size_t i = 0; i < 3; ++i) for (std::size_t j = 0; j < 3; ++j) { f[i][j] = nonlinear ? h[i][j] : HighLow{}; @@ -114,13 +117,48 @@ Result Hex8Element::EvaluateHighLow( for (std::size_t k = 0; k < 3; ++k) strain[i][j] += HighLow{0.5} * h[k][i] * h[k][j]; } + if (!Finite(h) || !Finite(f) || !Finite(strain)) + return Result::Failure(Failure( + "Hex8 high+low strain arithmetic is outside its supported range.")); + point_gradients[g] = gradients; + deformations[g] = f; + strains[g] = strain; + volume += HighLow{point.weight}; + trace_integral += + HighLow{point.weight} * (strain[0][0] + strain[1][1] + strain[2][2]); + } + if (!IsNormalized(volume) || math_internal::Compare(volume, HighLow{}) != 1 || + !IsNormalized(trace_integral)) + return Result::Failure( + Failure("Hex8 high+low reference volume or trace is unsupported.")); + const HighLow mean_trace = trace_integral / volume; + if (!IsNormalized(mean_trace)) + return Result::Failure( + Failure("Hex8 high+low weighted mean is unsupported.")); + const HighLow twice_mu = HighLow{2} * HighLow{mu_}; + const HighLow deviatoric_coefficient = + math_internal::Negate(twice_mu) / HighLow{3}; + const HighLow bulk_coefficient = HighLow{kappa_} / volume; + if (!IsNormalized(twice_mu) || !IsNormalized(deviatoric_coefficient) || + !IsNormalized(bulk_coefficient)) + return Result::Failure( + Failure("Hex8 high+low constitutive coefficients are unsupported.")); + energy = HighLow{0.5} * HighLow{kappa_} * volume * mean_trace * mean_trace; + for (std::size_t g = 0; g < 8; ++g) { + const auto& point = points_[g]; + const auto& gradients = point_gradients[g]; + const auto& f = deformations[g]; + const auto& strain = strains[g]; const HighLow trace = strain[0][0] + strain[1][1] + strain[2][2]; - HighLow density_energy = HighLow{0.5} * HighLow{lambda_} * trace * trace; + Tensor s{}, p{}, fft{}; + HighLow density_energy; for (std::size_t i = 0; i < 3; ++i) for (std::size_t j = 0; j < 3; ++j) { - s[i][j] = HighLow{2} * HighLow{mu_} * strain[i][j]; - if (i == j) s[i][j] += HighLow{lambda_} * trace; - density_energy += HighLow{mu_} * strain[i][j] * strain[i][j]; + const HighLow dev = + strain[i][j] - (i == j ? trace / HighLow{3} : HighLow{}); + s[i][j] = twice_mu * dev; + if (i == j) s[i][j] += HighLow{kappa_} * mean_trace; + density_energy += HighLow{mu_} * dev * dev; } energy += density_energy * HighLow{point.weight}; for (std::size_t i = 0; i < 3; ++i) @@ -144,8 +182,8 @@ Result Hex8Element::EvaluateHighLow( stress[i][j] = stress[i][j] / det; } } - if (!Finite(h) || !Finite(f) || !Finite(strain) || !Finite(s) || - !Finite(p) || !Finite(fft) || !Finite(stress)) + if (!Finite(f) || !Finite(strain) || !Finite(s) || !Finite(p) || + !Finite(fft) || !Finite(stress)) return Result::Failure(Failure( "Hex8 high+low tensor arithmetic is outside its supported range.")); response.high.points[g] = {element_index_, static_cast(g + 1U), @@ -160,6 +198,9 @@ Result Hex8Element::EvaluateHighLow( 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[a]); + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t i = 0; i < 3; ++i) + integrated_fg[a][i] += fg[a][i] * HighLow{point.weight}; for (std::size_t a = 0; a < 8; ++a) for (std::size_t b = 0; b < 8; ++b) { const auto dot = Dot(gradients[a], gradients[b]); @@ -172,7 +213,7 @@ Result Hex8Element::EvaluateHighLow( for (std::size_t k = 0; k < 3; ++k) { const auto row = 3U * a + i; const auto column = 3U * b + k; - const auto mat = HighLow{lambda_} * fg[a][i] * fg[b][k] + + const auto mat = deviatoric_coefficient * fg[a][i] * fg[b][k] + HighLow{mu_} * fg[b][i] * fg[a][k] + HighLow{mu_} * fft[i][k] * dot; material[row][column] += mat * HighLow{point.weight}; @@ -181,6 +222,13 @@ Result Hex8Element::EvaluateHighLow( } } } + if (nonlinear || !paired_linear_tangent_ready_) + 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) + material[3 * a + i][3 * b + k] += + bulk_coefficient * integrated_fg[a][i] * integrated_fg[b][k]; if (!IsNormalized(energy)) return Result::Failure( Failure("Hex8 high+low energy is unsupported or nonfinite.")); diff --git a/tests/unit/analysis/solid_frequency_procedure_test.cpp b/tests/unit/analysis/solid_frequency_procedure_test.cpp index 62479b5..508dcc5 100644 --- a/tests/unit/analysis/solid_frequency_procedure_test.cpp +++ b/tests/unit/analysis/solid_frequency_procedure_test.cpp @@ -99,8 +99,11 @@ std::vector DenseOracle(const std::vector& free) { for (std::size_t d = 0; d < 3; ++d) gradient_dot += GradientIntegral(a, b, d, d); k[row * n + col] = - 48.0 * (GradientIntegral(a, b, i, j) + GradientIntegral(a, b, j, i) + - (i == j ? gradient_dot : 0.0)); + 48.0 * + (-(2.0 / 3.0) * GradientIntegral(a, b, i, j) + + GradientIntegral(a, b, j, i) + (i == j ? gradient_dot : 0.0)) + + (80.0 / 48.0) * (48.0 * kSigns[a][i] / (8.0 * (i + 1))) * + (48.0 * kSigns[b][j] / (8.0 * (j + 1))); if (i == j) { m[row * n + col] = 2.0 * 48.0 / 64.0; for (std::size_t d = 0; d < 3; ++d) diff --git a/tests/unit/analysis/solid_hht_dynamic_procedure_test.cpp b/tests/unit/analysis/solid_hht_dynamic_procedure_test.cpp index 78c8fd4..4bdd4bf 100644 --- a/tests/unit/analysis/solid_hht_dynamic_procedure_test.cpp +++ b/tests/unit/analysis/solid_hht_dynamic_procedure_test.cpp @@ -573,7 +573,7 @@ TEST(SolidHhtDynamics, auto result = RunDynamics(OneFree()); ASSERT_TRUE(result.HasValue()) << solid_hht_test::Message(result.GetStatus()); const auto d = solid_hht_test::Brick(); - constexpr double lambda = 4.8; + constexpr double bulk = 8.0; constexpr double mu = 4.8; const double roundoff = 512 * std::numeric_limits::epsilon(); const auto& frames = result.Value().state.AcceptedFrames(); @@ -592,17 +592,20 @@ TEST(SolidHhtDynamics, double stiffness; double stiffness_term_scale; if (c == 0) { - stiffness = (lambda + 2 * mu) * (-sign[0]) * product[1] * product[2] + + stiffness = (4 * mu / 3) * (-sign[0]) * product[1] * product[2] + mu * (-sign[1]) * product[0] * product[2] + mu * (-sign[2]) * product[0] * product[1]; - stiffness_term_scale = (lambda + 2 * mu) * product[1] * product[2] + + stiffness_term_scale = (4 * mu / 3) * product[1] * product[2] + mu * product[0] * product[2] + mu * product[0] * product[1]; } else { const unsigned other = c == 1 ? 2 : 1; - stiffness = -.25 * (lambda * sign[c] + mu * sign[0]) * product[other]; - stiffness_term_scale = .25 * (lambda + mu) * product[other]; + stiffness = + -.25 * ((-2 * mu / 3) * sign[c] + mu * sign[0]) * product[other]; + stiffness_term_scale = .25 * (2 * mu / 3 + mu) * product[other]; } + stiffness -= bulk * sign[c] / 16; + stiffness_term_scale += bulk / 16; const double inertia = c == 0 ? scalar_mass * f.acceleration[0] : 0; const double internal = stiffness * f.displacement[0]; const double applied = b == 0 && c == 0 ? .01 : 0; diff --git a/tests/unit/analysis/solid_hht_temporal_verification_test.cpp b/tests/unit/analysis/solid_hht_temporal_verification_test.cpp index bc9f5ad..b9d29dd 100644 --- a/tests/unit/analysis/solid_hht_temporal_verification_test.cpp +++ b/tests/unit/analysis/solid_hht_temporal_verification_test.cpp @@ -24,7 +24,7 @@ AffineMotion Motion(double time, bool nonlinear) { const double dq = 2 * std::sin(2 * time); const double ddq = 4 * std::cos(2 * time); AffineMotion m; - constexpr double lambda = 4.8; + constexpr double bulk = 8; constexpr double mu = 4.8; if (!nonlinear) { const M3 d{{{.015, .008, -.004}, {.006, -.01, .003}, {.005, .004, .012}}}; @@ -35,8 +35,8 @@ AffineMotion Motion(double time, bool nonlinear) { m.b[i][j] = q * d[i][j]; m.velocity[i][j] = dq * d[i][j]; m.acceleration[i][j] = ddq * d[i][j]; - m.piola[i][j] = - (i == j ? lambda * trace : 0) + mu * q * (d[i][j] + d[j][i]); + m.piola[i][j] = (i == j ? (bulk * trace - 2 * mu * trace / 3) : 0) + + mu * q * (d[i][j] + d[j][i]); } } else { const double s = .004 * q, ds = .004 * dq, dds = .004 * ddq; @@ -46,7 +46,7 @@ AffineMotion Motion(double time, bool nonlinear) { const M3 jr{{{-sn, -c, 0}, {c, -sn, 0}, {0, 0, 0}}}; const M3 jjr{{{-c, sn, 0}, {-sn, -c, 0}, {0, 0, 0}}}; const double eg = .5 * ((1 + s) * (1 + s) - 1); - const double stress = (3 * lambda + 2 * mu) * eg; + const double stress = 3 * bulk * eg; for (unsigned i = 0; i < 3; ++i) for (unsigned j = 0; j < 3; ++j) { m.b[i][j] = (1 + s) * r[i][j] - (i == j ? 1 : 0); @@ -506,7 +506,7 @@ TEST(SolidHhtTemporal, TEST(SolidHhtTemporal, IndependentThreeDofBrickOscillatorHasSecondOrderUVA) { // Node0's three translations are free. Exact polynomial integrals yield - // M=(rho/27)I, Kii=(lambda+4mu)/9, Kij=(lambda+mu)/12. + // M=(rho/27)I, Kii=41/18, Kij=19/30 for mu=4.8 and direct kappa=8. std::array, 3> errors{}; for (unsigned level = 0; level < 3; ++level) { const double h = .02 / static_cast(1U << level); @@ -529,8 +529,8 @@ TEST(SolidHhtTemporal, IndependentThreeDofBrickOscillatorHasSecondOrderUVA) { << solid_hht_test::Message(result.GetStatus()); const auto& f = result.Value().state.CommittedState(); const double m = 2.0 / 27; - const double diagonal = (4.8 + 4 * 4.8) / 9; - const double off = (4.8 + 4.8) / 12; + const double diagonal = 41.0 / 18; + const double off = 19.0 / 30; const std::array k{diagonal + 2 * off, diagonal - off}; for (unsigned c = 0; c < 3; ++c) { std::array truth{}; @@ -563,7 +563,7 @@ TEST(SolidHhtTemporal, IndependentThreeDofBrickOscillatorHasSecondOrderUVA) { TEST(SolidHhtTemporal, IndependentScalarBrickOscillatorHasSecondOrderUVA) { std::array, 3> errors{}; const double m = 2.0 / 27; - const double k = (4.8 + 4 * 4.8) / 9; + const double k = 41.0 / 18; const double w = std::sqrt(k / m); for (unsigned level = 0; level < 3; ++level) { const double h = .02 / static_cast(1U << level); @@ -608,46 +608,95 @@ TEST(SolidHhtTemporal, IndependentScalarBrickOscillatorHasSecondOrderUVA) { } } -TEST(SolidHhtTemporal, - BendingMeshRefinementQuantifiesFullIntegrationShearLocking) { - constexpr double thickness = .125; - constexpr double curvature = .01; - constexpr double lambda = 4.8; - constexpr double mu = 4.8; - const double continuum = .5 * (lambda + 2 * mu) * curvature * curvature * - thickness * thickness * thickness / 12; - double previous_error = 0; - for (unsigned nx : {1U, 2U, 4U, 8U}) { - double numerical = 0; - for (unsigned e = 0; e < nx; ++e) { - auto d = solid_hht_test::Brick(); +namespace { +// Compatible nodal pure-bending field; no equilibrium solve is claimed. +std::array BendingEnergy(unsigned nx, unsigned ny) { + constexpr double thickness = .125, curvature = .01, bulk = 8; + std::array result{}; + auto d = solid_hht_test::Brick(); + for (unsigned ex = 0; ex < nx; ++ex) + for (unsigned ey = 0; ey < ny; ++ey) { + const auto index = ex * ny + ey; fesa::Hex8Element::Coordinates x{}; fesa::ElementDofLayout layout{ - {"I", e + 1, std::to_string(e + 1)}, + {"I", index + 1, std::to_string(index + 1)}, {0, 1, 2, 3, 4, 5, 6, 7}, {fesa::DofComponent::kUx, fesa::DofComponent::kUy, fesa::DofComponent::kUz}}; fesa::Vector u(24); for (unsigned a = 0; a < 8; ++a) { - x[a] = {static_cast(e) / nx + d.nodes[a].coordinates[0] / nx, - thickness * (d.nodes[a].coordinates[1] - .5), + x[a] = {(ex + d.nodes[a].coordinates[0]) / nx, + thickness * ((ey + d.nodes[a].coordinates[1]) / ny - .5), d.nodes[a].coordinates[2]}; u[3 * a] = -curvature * x[a][0] * x[a][1]; u[3 * a + 1] = .5 * curvature * x[a][0] * x[a][0]; } - auto element = fesa::Hex8Element::Create(e, layout, x, d.materials[0], 2); - ASSERT_TRUE(element.HasValue()); + auto element = + fesa::Hex8Element::Create(index, layout, x, d.materials[0], 2); + EXPECT_TRUE(element.HasValue()); + if (!element.HasValue()) return {}; auto response = element.Value().EvaluateLinear(u); - ASSERT_TRUE(response.HasValue()); - numerical += response.Value().strain_energy; + EXPECT_TRUE(response.HasValue()); + if (!response.HasValue()) return {}; + const auto weights = element.Value().ReferenceWeights(); + double volume = 0, trace_integral = 0; + for (unsigned g = 0; g < 8; ++g) { + volume += weights[g]; + const auto& strain = response.Value().points[g].strain; + trace_integral += weights[g] * (strain[0] + strain[1] + strain[2]); + } + const double volumetric = + .5 * bulk * trace_integral * trace_integral / volume; + result[0] += response.Value().strain_energy; + result[1] += response.Value().strain_energy - volumetric; + result[2] += volumetric; } - const double error = numerical / continuum - 1; - const double exact_parasitic = - mu / (lambda + 2 * mu) / std::pow(nx * thickness, 2); - EXPECT_NEAR(error, exact_parasitic, 2e-11); - if (previous_error != 0) EXPECT_NEAR(previous_error / error, 4, 1e-10); - std::cout << "Bending nx=" << nx - << " energy/continuum=" << numerical / continuum << '\n'; - previous_error = error; + return result; +} +} // namespace + +TEST(SolidHhtTemporal, FixedThicknessBendingRetainsDiscreteBulkBias) { + constexpr double thickness = .125, curvature = .01, mu = 4.8, bulk = 8; + constexpr double axial = bulk + 4 * mu / 3; + const double continuum = + axial * curvature * curvature * std::pow(thickness, 3) / 24; + const double limit = 4 * mu / (3 * axial); + double previous_parasitic = 0; + for (unsigned nx : {1U, 2U, 4U, 8U}) { + const auto energy = BendingEnergy(nx, 1); + const double parasitic = mu / (axial * std::pow(nx * thickness, 2)); + EXPECT_NEAR(energy[0] / continuum, limit + parasitic, 2e-11); + EXPECT_NEAR(energy[2] / continuum, 0, 2e-13); + if (previous_parasitic != 0) + EXPECT_NEAR(previous_parasitic / parasitic, 4, 1e-13); + previous_parasitic = parasitic; + std::cout << "Fixed thickness nx=" << nx + << " U/Ucont=" << energy[0] / continuum + << " limiting ratio=" << limit << " shear=" << parasitic << '\n'; + } +} + +TEST(SolidHhtTemporal, + BendingMeshRefinementQuantifiesFullIntegrationShearLocking) { + constexpr double thickness = .125, curvature = .01, mu = 4.8, bulk = 8; + constexpr double axial = bulk + 4 * mu / 3; + const double continuum = + axial * curvature * curvature * std::pow(thickness, 3) / 24; + double previous_error = 0; + for (unsigned n : {1U, 2U, 4U, 8U}) { + const auto energy = BendingEnergy(n, n); + const double shear = mu / (axial * std::pow(n * thickness, 2)); + const double missing_bulk = bulk / (axial * n * n); + const double error = energy[0] / continuum - 1; + EXPECT_NEAR(error, shear - missing_bulk, 2e-11); + // Independent moments separate competing errors, so cancellation cannot + // make the continuum-convergence claim pass accidentally. + EXPECT_NEAR(energy[1] / continuum, 4 * mu / (3 * axial) + shear, 2e-11); + EXPECT_NEAR(energy[2] / continuum, bulk / axial - missing_bulk, 2e-11); + if (previous_error != 0) EXPECT_NEAR(previous_error / error, 4, 1e-10); + previous_error = error; + std::cout << "Bending nx=ny=" << n << " U/Ucont=" << energy[0] / continuum + << " positive shear=" << shear + << " negative bulk=" << missing_bulk << '\n'; } } diff --git a/tests/unit/elements/hex8_element_test.cpp b/tests/unit/elements/hex8_element_test.cpp index 8395968..706bd9a 100644 --- a/tests/unit/elements/hex8_element_test.cpp +++ b/tests/unit/elements/hex8_element_test.cpp @@ -11,6 +11,8 @@ #include #include +#include "../../../src/fesa/math/high_low_arithmetic.h" + namespace { using Coordinates = fesa::Hex8Element::Coordinates; constexpr int kSigns[8][3] = {{-1, -1, -1}, {1, -1, -1}, {1, 1, -1}, @@ -64,11 +66,14 @@ TEST(Hex8Element, ExactBrickStiffnessAndConsistentMass) { scalar *= kSigns[a][d] == kSigns[b][d] ? 2.0 : 1.0; for (std::size_t i = 0; i < 3; ++i) for (std::size_t j = 0; j < 3; ++j) { - double exact = 48.0 * GradientProduct(a, b, i, j) + + double exact = -(2.0 / 3.0) * 48.0 * GradientProduct(a, b, i, j) + 48.0 * GradientProduct(a, b, j, i); if (i == j) for (std::size_t d = 0; d < 3; ++d) exact += 48.0 * GradientProduct(a, b, d, d); + // Mean bulk term is the product of two independent first moments. + exact += (80.0 / 48.0) * (48.0 * kSigns[a][i] / (8.0 * (i + 1))) * + (48.0 * kSigns[b][j] / (8.0 * (j + 1))); EXPECT_NEAR(k(3 * a + i, 3 * b + j), exact, 2e-12); EXPECT_NEAR(m(3 * a + i, 3 * b + j), i == j ? scalar : 0.0, 2e-14); } @@ -542,8 +547,9 @@ TEST(Hex8Element, EnergyForceAndForceTangentDirectionalStepSweep) { ASSERT_TRUE(trial.HasValue()); const double derivative = direction.Dot(trial.Value().internal_force); const auto tangent_direction = trial.Value().tangent.Multiply(direction); - std::array energy_error{}, force_error{}; - const std::array steps{1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7}; + std::array energy_error{}, force_error{}; + const std::array steps{1e-2, 1e-3, 1e-4, 1e-5, + 1e-6, 1e-7, 1e-8, 1e-9}; for (std::size_t s = 0; s < steps.size(); ++s) { auto plus = u, minus = u; plus.Axpy(steps[s], direction); @@ -751,14 +757,14 @@ TEST(Hex8Element, RejectsUnrepresentableShearAndMassInsteadOfReturningZeros) { ASSERT_TRUE(element.HasValue()); EXPECT_FALSE(element.Value().ComputeMass().HasValue()); } -TEST(Hex8Element, RejectsRoundedLossOfPositiveVolumetricModulus) { - // Binary64 Lame arithmetic at this valid nu yields lambda+2mu/3 == 0 - // even though exact E/[3(1-2nu)] is approximately1/9, which is positive. +TEST(Hex8Element, DirectBulkSurvivesNearNegativePoissonBound) { const double nu = std::nextafter(-1.0, 0.0); auto element = Make(Brick(), 1.0, nu); - ASSERT_FALSE(element.HasValue()); - EXPECT_NE(element.GetStatus().Diagnostics()[0].message.find("bulk"), - std::string::npos); + ASSERT_TRUE(element.HasValue()); + // Direct kappa is positive even where a lambda-based difference loses it. + EXPECT_GT(1.0 / (3.0 * (1.0 - 2.0 * nu)), 0.0); + fesa::Vector zero(24); + EXPECT_TRUE(element.Value().EvaluateLinear(zero).HasValue()); } TEST(Hex8Element, OrdinaryNegativePoissonRatioRetainsPositiveVolumetricEnergy) { auto element = Make(Brick(), 120.0, -.5); @@ -838,9 +844,11 @@ Triple Cross(const Triple& a, const Triple& b) { a[0] * b[1] - a[1] * b[0]}; } std::array WarpedQuadratureEnergies(const fesa::Vector& u, - unsigned count) { + unsigned count, + bool nonlinear = false) { const auto rule = LegendreRule(count); std::array result{}; + double volume = 0.0, trace_integral = 0.0; for (const auto& xi : rule) for (const auto& eta : rule) for (const auto& zeta : rule) { @@ -877,19 +885,30 @@ std::array WarpedQuadratureEnergies(const fesa::Vector& u, for (std::size_t j = 0; j < 3; ++j) for (std::size_t d = 0; d < 3; ++d) gradient[i][j] += natural_gradient[i][d] * reciprocal[d][j]; + Tensor strains{}; + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) { + strains[i][j] = .5 * (gradient[i][j] + gradient[j][i]); + if (nonlinear) + for (std::size_t k = 0; k < 3; ++k) + strains[i][j] += .5 * gradient[k][i] * gradient[k][j]; + } double trace = 0.0, strain_squared = 0.0, velocity_squared = 0.0; for (std::size_t i = 0; i < 3; ++i) { - trace += gradient[i][i]; + trace += strains[i][i]; velocity_squared += velocity[i] * velocity[i]; for (std::size_t j = 0; j < 3; ++j) { - const double strain = .5 * (gradient[i][j] + gradient[j][i]); + const double strain = strains[i][j]; strain_squared += strain * strain; } } const double weight = xi[1] * eta[1] * zeta[1] * det; - result[0] += weight * (24.0 * trace * trace + 48.0 * strain_squared); + result[0] += weight * 48.0 * (strain_squared - trace * trace / 3.0); + volume += weight; + trace_integral += weight * trace; result[1] += weight * velocity_squared; // .5*rho*v^2 with rho=2. } + result[0] += 40.0 * trace_integral * trace_integral / volume; return result; } TEST(Hex8Element, WarpedNonaffineQuadratureDifferenceIsQuantifiedAndRefines) { @@ -922,4 +941,199 @@ TEST(Hex8Element, WarpedNonaffineQuadratureDifferenceIsQuantifiedAndRefines) { << " K3-vs-K7=" << std::abs(q3[0] - q7[0]) / q7[0] << " K5-vs-K7=" << std::abs(q5[0] - q7[0]) / q7[0] << '\n'; } +Coordinates MeanVolumeUnitBrick() { + Coordinates x{}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t i = 0; i < 3; ++i) x[a][i] = .5 * (kSigns[a][i] + 1); + return x; +} + +TEST(Hex8Element, + MeanVolumeNonaffinePolynomialEnergyForceHessianAndActualStrain) { + const auto x = MeanVolumeUnitBrick(); + const auto element = Make(x, 12.0, .25); + ASSERT_TRUE(element.HasValue()); + constexpr double mu = 4.8, bulk = 8.0, c = .125; + // Exact monomial integrals of u_x=c*x*y over the unit cube. The other + // translations vanish; this direction varies trace among all eight points. + constexpr double a = 7 * mu / 18 + bulk / 8; + constexpr double b = 5 * mu / 18 + bulk / 6; + // Fixed two-point quadrature has Q[x^4]=7/36, not the continuum 1/5. + constexpr double d = 11 * mu / 108 + bulk / 18; + fesa::Vector u(24); + for (std::size_t n = 0; n < 8; ++n) u[3 * n] = c * x[n][0] * x[n][1]; + for (bool nonlinear : {false, true}) { + const auto response = nonlinear ? element.Value().EvaluateNonlinear(u) + : element.Value().EvaluateLinear(u); + ASSERT_TRUE(response.HasValue()); + const auto& r = response.Value(); + const double energy = + a * c * c + (nonlinear ? b * c * c * c + d * c * c * c * c : 0); + const double force = + 2 * a * c + (nonlinear ? 3 * b * c * c + 4 * d * c * c * c : 0); + const double hessian = 2 * a + (nonlinear ? 6 * b * c + 12 * d * c * c : 0); + EXPECT_NEAR(r.strain_energy, energy, 2e-14); + EXPECT_NEAR(r.internal_force[6] + r.internal_force[18], force, 2e-13); + EXPECT_NEAR(r.tangent(6, 6) + r.tangent(6, 18) + r.tangent(18, 6) + + r.tangent(18, 18), + hessian, 2e-12); + const double mean = c / 2 + (nonlinear ? c * c / 3 : 0); + for (const auto& point : r.points) { + const double px = point.reference_coordinates[0], + py = point.reference_coordinates[1]; + const double t = 1 + c * py; + const double exx = c * py + (nonlinear ? .5 * c * c * py * py : 0); + const double eyy = nonlinear ? .5 * c * c * px * px : 0; + const double exy = .5 * c * px + (nonlinear ? .5 * c * c * px * py : 0); + const double trace = exx + eyy; + const double sxx = 2 * mu * (exx - trace / 3) + bulk * mean; + const double syy = 2 * mu * (eyy - trace / 3) + bulk * mean; + const double szz = -2 * mu * trace / 3 + bulk * mean; + const double sxy = 2 * mu * exy; + EXPECT_NEAR(point.strain[0], exx, 2e-15); + EXPECT_NEAR(point.strain[1], eyy, 2e-15); + EXPECT_NEAR(point.strain[3], exy, 2e-15); + EXPECT_NEAR( + point.cauchy_stress[0], + nonlinear + ? (t * t * sxx + 2 * t * c * px * sxy + c * c * px * px * syy) / t + : sxx, + 2e-13); + EXPECT_NEAR(point.cauchy_stress[1], nonlinear ? syy / t : syy, 2e-13); + EXPECT_NEAR(point.cauchy_stress[2], nonlinear ? szz / t : szz, 2e-13); + EXPECT_NEAR(point.cauchy_stress[3], + nonlinear ? (t * sxy + c * px * syy) / t : sxy, 2e-13); + } + } + const auto k = element.Value().ComputeStiffness().Value().values; + EXPECT_NEAR(k(0, 0), 41.0 / 18, 2e-14); + EXPECT_NEAR(k(0, 1), 19.0 / 30, 2e-14); +} + +TEST(Hex8Element, MeanVolumeWarpedStressUsesReferenceWeightedTrace) { + auto element = Make(Warped()); + ASSERT_TRUE(element.HasValue()); + fesa::Vector u(24); + u[1] = 1; // Independent derivative coefficients separate the two means by + // >1e-8. + const auto response = element.Value().EvaluateLinear(u); + ASSERT_TRUE(response.HasValue()); + const auto weights = element.Value().ReferenceWeights(); + double volume = 0, integral = 0, arithmetic = 0; + for (std::size_t g = 0; g < 8; ++g) { + const auto& e = response.Value().points[g].strain; + const double trace = e[0] + e[1] + e[2]; + volume += weights[g]; + integral += weights[g] * trace; + arithmetic += trace / 8; + } + const double mean = integral / volume; + ASSERT_GT(std::abs(mean - arithmetic), 1e-8); + for (const auto& point : response.Value().points) { + const auto& s = point.cauchy_stress; + EXPECT_NEAR((s[0] + s[1] + s[2]) / 3, 80 * mean, 2e-13); + } + EXPECT_NEAR(response.Value().strain_energy, WarpedQuadratureEnergies(u, 2)[0], + 2e-12); +} + +TEST(Hex8Element, + MeanVolumeDirectBulkSupportsBoundsWithoutLambdaOrThreeBulkGuard) { + EXPECT_TRUE( + Make(MeanVolumeUnitBrick(), 12, std::nextafter(-1.0, 0.0)).HasValue()); + EXPECT_TRUE( + Make(MeanVolumeUnitBrick(), 1, std::nextafter(.5, 0.0)).HasValue()); + const auto huge = Make(MeanVolumeUnitBrick(), 1e308, .4); + ASSERT_TRUE(huge.HasValue()); + const double direct_bulk = 1e308 / (3 * (1 - 2 * .4)); + ASSERT_TRUE(std::isfinite(direct_bulk)); + ASSERT_FALSE(std::isfinite(3 * direct_bulk)); + fesa::Vector u(24); + constexpr double shear = 1e-150; + const auto x = MeanVolumeUnitBrick(); + for (unsigned n = 0; n < 8; ++n) u[3 * n] = shear * x[n][1]; + const auto response = huge.Value().EvaluateLinear(u); + ASSERT_TRUE(response.HasValue()); + const double expected = .5 * (1e308 / (2 * (1 + .4))) * shear * shear; + EXPECT_NEAR(response.Value().strain_energy, expected, 1e-13 * expected); + EXPECT_FALSE(Make(MeanVolumeUnitBrick(), 1e308, .499).HasValue()); + EXPECT_FALSE(Make(MeanVolumeUnitBrick(), + std::numeric_limits::denorm_min(), -.5) + .HasValue()); +} +TEST(Hex8Element, + MeanVolumeWarpedIndependentEnergyHasCoupledNodeZeroDerivatives) { + auto element = Make(Warped()); + ASSERT_TRUE(element.HasValue()); + fesa::Vector u(24), direction(24), zero(24); + for (std::size_t i = 0; i < 24; ++i) u[i] = .03 * std::sin(.7 * i); + direction[0] = 1; + direction[8] = -.375; + direction[22] = .25; + constexpr double h = .01; + for (bool nonlinear : {false, true}) { + const auto response = nonlinear ? element.Value().EvaluateNonlinear(u) + : element.Value().EvaluateLinear(u); + const auto paired = nonlinear + ? element.Value().EvaluateNonlinearHighLow(u, zero) + : element.Value().EvaluateLinearHighLow(u, zero); + ASSERT_TRUE(response.HasValue()); + ASSERT_TRUE(paired.HasValue()); + const double energy = WarpedQuadratureEnergies(u, 2, nonlinear)[0]; + EXPECT_NEAR(response.Value().strain_energy, energy, 2e-12); + EXPECT_NEAR(paired.Value().high.strain_energy, energy, 2e-12); + std::array first{}, second{}; + for (unsigned level = 0; level < 2; ++level) { + const double step = h / (1U << level); + auto plus = u, minus = u; + plus.Axpy(step, direction); + minus.Axpy(-step, direction); + const double ep = WarpedQuadratureEnergies(plus, 2, nonlinear)[0]; + const double em = WarpedQuadratureEnergies(minus, 2, nonlinear)[0]; + first[level] = (ep - em) / (2 * step); + second[level] = (ep - 2 * energy + em) / (step * step); + } + // This energy is quartic in nodal translations; Richardson removes the + // exact centered-difference h² term, independent of either kernel path. + const double df = (4 * first[1] - first[0]) / 3; + const double ddf = (4 * second[1] - second[0]) / 3; + EXPECT_NEAR(direction.Dot(response.Value().internal_force), df, 2e-10); + EXPECT_NEAR(direction.Dot(paired.Value().high.internal_force), df, 2e-10); + EXPECT_NEAR(direction.Dot(response.Value().tangent.Multiply(direction)), + ddf, 2e-8); + EXPECT_NEAR(direction.Dot(paired.Value().high.tangent.Multiply(direction)), + ddf, 2e-8); + } +} + +TEST(Hex8Element, MeanVolumeWarpedLowOnlyStateRetainsWeightedCubicEnergy) { + using namespace fesa::math_internal; + auto element = Make(Warped()); + ASSERT_TRUE(element.HasValue()); + const double c = std::ldexp(1.0, -60); + fesa::Vector high(24), low(24), zero(24), basis(24), minus(24); + for (unsigned a = 0; a < 8; ++a) + for (unsigned i = 0; i < 3; ++i) high[3 * a + i] = 1; + low[1] = c; + basis[1] = 1; + minus[1] = -1; + const auto linear = element.Value().EvaluateLinearHighLow(high, low); + const auto nonlinear = element.Value().EvaluateNonlinearHighLow(high, low); + const auto high_only = element.Value().EvaluateNonlinearHighLow(high, zero); + ASSERT_TRUE(linear.HasValue()); + ASSERT_TRUE(nonlinear.HasValue()); + ASSERT_TRUE(high_only.HasValue()); + EXPECT_EQ(high_only.Value().high.strain_energy, 0); + const double quadratic = WarpedQuadratureEnergies(basis, 2)[0]; + const double cubic = .5 * (WarpedQuadratureEnergies(basis, 2, true)[0] - + WarpedQuadratureEnergies(minus, 2, true)[0]); + EXPECT_NEAR(linear.Value().high.strain_energy, quadratic * c * c, + 1e-12 * c * c); + const auto actual_cubic = Subtract( + {nonlinear.Value().high.strain_energy, + nonlinear.Value().strain_energy_low}, + {linear.Value().high.strain_energy, linear.Value().strain_energy_low}); + EXPECT_NEAR(actual_cubic.high, cubic * c * c * c, 1e-11 * c * c * c); + EXPECT_NE(actual_cubic.high, 0); +} } // namespace diff --git a/tests/unit/elements/hex8_high_low_test.cpp b/tests/unit/elements/hex8_high_low_test.cpp index e837e0c..add510b 100644 --- a/tests/unit/elements/hex8_high_low_test.cpp +++ b/tests/unit/elements/hex8_high_low_test.cpp @@ -295,4 +295,40 @@ TEST(Hex8HighLow, FiniteRotationRetainsTinyShearForceCovariance) { 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