fix(elements): compensate cancelled Hex8 linear force products
This commit is contained in:
@@ -243,7 +243,22 @@ Result<Hex8Response> 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) {
|
||||
|
||||
Reference in New Issue
Block a user