diff --git a/src/fesa/elements/hex8_element.cpp b/src/fesa/elements/hex8_element.cpp index 09aab50..cf7facb 100644 --- a/src/fesa/elements/hex8_element.cpp +++ b/src/fesa/elements/hex8_element.cpp @@ -243,7 +243,22 @@ Result Hex8Element::Evaluate(const Vector& u, for (std::size_t a = 0; a < 8; ++a) for (std::size_t i = 0; i < 3; ++i) relative[3U * a + i] = u[3U * a + i] - u[i]; - response.internal_force = stiffness_.Multiply(relative); + // 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) { + const double product = stiffness_(i, j) * relative[j]; + const double next = sum + product; + correction += std::fma(stiffness_(i, j), relative[j], -product) + + (std::abs(sum) >= std::abs(product) + ? (sum - next) + product + : (product - next) + sum); + sum = next; + } + response.internal_force[i] = sum + correction; + } response.strain_energy = 0.5 * relative.Dot(response.internal_force); } for (std::size_t g = 0; g < 8; ++g) {