feat(elements): implement certified full-integration Hex8 kernel
This commit is contained in:
KOKO\Mimi
@@ -0,0 +1,325 @@
|
||||
#include "fesa/elements/hex8_element.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <set>
|
||||
#include <utility>
|
||||
|
||||
#include "hex8_geometry.h"
|
||||
|
||||
namespace fesa {
|
||||
namespace {
|
||||
using hex8_internal::Tensor;
|
||||
using hex8_internal::Triple;
|
||||
/// @brief Checks every tensor component before recovery or inverse use.
|
||||
bool Finite(const Tensor& value) {
|
||||
for (const auto& row : value)
|
||||
for (double entry : row)
|
||||
if (!std::isfinite(entry)) return false;
|
||||
return true;
|
||||
}
|
||||
/// @brief Packs actual tensor shear, never engineering shear, in schema order.
|
||||
std::array<double, 6> Pack(const Tensor& value) {
|
||||
return {value[0][0], value[1][1], value[2][2],
|
||||
value[0][1], value[0][2], value[1][2]};
|
||||
}
|
||||
/// @brief Computes a Euclidean contraction in deterministic component order.
|
||||
double Dot(const Triple& a, const Triple& b) {
|
||||
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
|
||||
}
|
||||
/// @brief Gives the finite-control failure without accepting a geometry guess.
|
||||
const char* CertificateMessage(hex8_internal::Certificate certificate) {
|
||||
switch (certificate) {
|
||||
case hex8_internal::Certificate::kNonpositive:
|
||||
return "Whole-cell Jacobian has a nonpositive witness.";
|
||||
case hex8_internal::Certificate::kUncertain:
|
||||
return "Whole-cell Jacobian positivity is unresolved at depth 10 or cell "
|
||||
"budget 65536.";
|
||||
case hex8_internal::Certificate::kNonfinite:
|
||||
return "Whole-cell Jacobian interval arithmetic is nonfinite.";
|
||||
case hex8_internal::Certificate::kPositive:
|
||||
break;
|
||||
}
|
||||
return "Invalid Jacobian certificate.";
|
||||
}
|
||||
} // namespace
|
||||
|
||||
Hex8Element::Hex8Element(EntityIndex index, ElementDofLayout layout,
|
||||
Coordinates coordinates, double lambda, double mu,
|
||||
std::optional<double> density, SourceLocation location)
|
||||
: element_index_{index},
|
||||
layout_{std::move(layout)},
|
||||
coordinates_{coordinates},
|
||||
lambda_{lambda},
|
||||
mu_{mu},
|
||||
density_{density},
|
||||
location_{std::move(location)} {}
|
||||
|
||||
Status Hex8Element::Failure(const char* message) const {
|
||||
return Status::Failure(FailureCategory::kModel,
|
||||
{{Severity::kError, "invalid-hex8", location_, "C3D8",
|
||||
layout_.source_id.instance_name + ":" +
|
||||
layout_.source_id.source_label_text,
|
||||
message}});
|
||||
}
|
||||
|
||||
Result<Hex8Element> Hex8Element::Create(EntityIndex index,
|
||||
ElementDofLayout layout,
|
||||
Coordinates coordinates,
|
||||
const LinearElasticMaterial& material,
|
||||
std::optional<double> density,
|
||||
SourceLocation location) {
|
||||
const double e = material.YoungsModulus();
|
||||
const double nu = material.PoissonsRatio();
|
||||
Hex8Element candidate{index, std::move(layout), coordinates, 0.0,
|
||||
0.0, density, std::move(location)};
|
||||
const std::vector<DofComponent> components{
|
||||
DofComponent::kUx, DofComponent::kUy, DofComponent::kUz};
|
||||
const std::set<EntityIndex> nodes(candidate.layout_.node_indices.begin(),
|
||||
candidate.layout_.node_indices.end());
|
||||
if (candidate.layout_.node_indices.size() != 8U || nodes.size() != 8U ||
|
||||
candidate.layout_.components_per_node != components)
|
||||
return Result<Hex8Element>::Failure(candidate.Failure(
|
||||
"Hex8 requires eight distinct nodes with exactly UX, UY, UZ."));
|
||||
if (!std::isfinite(e) || e <= 0.0 || !std::isfinite(nu) || nu <= -1.0 ||
|
||||
nu >= 0.5)
|
||||
return Result<Hex8Element>::Failure(
|
||||
candidate.Failure("Hex8 requires finite E>0 and -1<nu<0.5."));
|
||||
candidate.mu_ = e / (2.0 * (1.0 + nu));
|
||||
candidate.lambda_ = e * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
|
||||
if (!std::isfinite(candidate.mu_) || candidate.mu_ <= 0.0 ||
|
||||
!std::isfinite(candidate.lambda_) ||
|
||||
!std::isfinite(candidate.lambda_ + 2.0 * candidate.mu_))
|
||||
return Result<Hex8Element>::Failure(
|
||||
candidate.Failure("Hex8 Lame moduli must remain finite with "
|
||||
"representable positive shear modulus."));
|
||||
if (density && (!std::isfinite(*density) || *density <= 0.0))
|
||||
return Result<Hex8Element>::Failure(candidate.Failure(
|
||||
"Hex8 reference density must be finite and positive."));
|
||||
const auto certificate = hex8_internal::Certify(coordinates);
|
||||
if (certificate != hex8_internal::Certificate::kPositive)
|
||||
return Result<Hex8Element>::Failure(
|
||||
candidate.Failure(CertificateMessage(certificate)));
|
||||
const double gauss = 1.0 / std::sqrt(3.0);
|
||||
for (std::size_t g = 0; g < 8; ++g) {
|
||||
auto& point = candidate.points_[g];
|
||||
point.natural = {g % 2U == 0U ? -gauss : gauss,
|
||||
(g / 2U) % 2U == 0U ? -gauss : gauss,
|
||||
g < 4U ? -gauss : gauss};
|
||||
std::array<Triple, 8> derivatives{};
|
||||
hex8_internal::Shape(point.natural, point.shape, derivatives);
|
||||
const auto jacobian = hex8_internal::Jacobian(coordinates, derivatives);
|
||||
point.weight = hex8_internal::Determinant(jacobian);
|
||||
if (!Finite(jacobian) || !std::isfinite(point.weight) ||
|
||||
point.weight <= 0.0)
|
||||
return Result<Hex8Element>::Failure(candidate.Failure(
|
||||
"Reference quadrature Jacobian must be finite and positive."));
|
||||
const auto inverse = hex8_internal::Inverse(jacobian, point.weight);
|
||||
if (!Finite(inverse))
|
||||
return Result<Hex8Element>::Failure(
|
||||
candidate.Failure("Reference Jacobian inverse is nonfinite."));
|
||||
point.position = coordinates[0];
|
||||
for (std::size_t a = 0; a < 8; ++a)
|
||||
for (std::size_t i = 0; i < 3; ++i) {
|
||||
point.position[i] +=
|
||||
point.shape[a] * (coordinates[a][i] - coordinates[0][i]);
|
||||
for (std::size_t d = 0; d < 3; ++d)
|
||||
point.gradients[a][i] += inverse[d][i] * derivatives[a][d];
|
||||
if (!std::isfinite(point.gradients[a][i]) ||
|
||||
!std::isfinite(point.position[i]))
|
||||
return Result<Hex8Element>::Failure(candidate.Failure(
|
||||
"Reference gradient or point coordinate is nonfinite."));
|
||||
}
|
||||
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 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] +
|
||||
candidate.mu_ * ga[k] * gb[i] +
|
||||
(i == k ? candidate.mu_ * dot : 0.0));
|
||||
}
|
||||
}
|
||||
}
|
||||
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)))
|
||||
return Result<Hex8Element>::Failure(
|
||||
candidate.Failure("Hex8 reference stiffness is nonfinite."));
|
||||
for (std::size_t i = 0; i < 24; ++i)
|
||||
if (candidate.stiffness_(i, i) <= 0.0)
|
||||
return Result<Hex8Element>::Failure(
|
||||
candidate.Failure("Hex8 reference stiffness has an unrepresentable "
|
||||
"positive diagonal."));
|
||||
return Result<Hex8Element>::Success(std::move(candidate));
|
||||
}
|
||||
|
||||
const ElementDofLayout& Hex8Element::DofLayout() const noexcept {
|
||||
return layout_;
|
||||
}
|
||||
Result<ElementStiffnessContribution> Hex8Element::ComputeStiffness() const {
|
||||
return Result<ElementStiffnessContribution>::Success({layout_, stiffness_});
|
||||
}
|
||||
Result<ElementStiffnessContribution> Hex8Element::ComputeMass() const {
|
||||
if (!density_)
|
||||
return Result<ElementStiffnessContribution>::Failure(
|
||||
Failure("Hex8 mass requires reference density."));
|
||||
Matrix mass{24U, 24U};
|
||||
for (const auto& point : points_)
|
||||
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;
|
||||
if (!std::isfinite(value))
|
||||
return Result<ElementStiffnessContribution>::Failure(
|
||||
Failure("Hex8 reference mass is nonfinite."));
|
||||
for (std::size_t i = 0; i < 3; ++i)
|
||||
mass(3U * a + i, 3U * b + i) += value;
|
||||
}
|
||||
for (std::size_t i = 0; i < 24; ++i)
|
||||
for (std::size_t j = 0; j < 24; ++j)
|
||||
if (!std::isfinite(mass(i, j)))
|
||||
return Result<ElementStiffnessContribution>::Failure(
|
||||
Failure("Hex8 accumulated reference mass is nonfinite."));
|
||||
for (std::size_t i = 0; i < 24; ++i)
|
||||
if (mass(i, i) <= 0.0)
|
||||
return Result<ElementStiffnessContribution>::Failure(Failure(
|
||||
"Hex8 reference mass has an unrepresentable positive diagonal."));
|
||||
return Result<ElementStiffnessContribution>::Success(
|
||||
{layout_, std::move(mass)});
|
||||
}
|
||||
Result<ElementResultBundle> Hex8Element::Recover(const Vector& u) const {
|
||||
auto response = EvaluateLinear(u);
|
||||
if (!response.HasValue())
|
||||
return Result<ElementResultBundle>::Failure(response.GetStatus());
|
||||
SolidElementResultRows rows;
|
||||
rows.rows.assign(response.Value().points.begin(),
|
||||
response.Value().points.end());
|
||||
rows.physical_strain_energy = response.Value().strain_energy;
|
||||
return Result<ElementResultBundle>::Success(
|
||||
{layout_.source_id, std::move(rows)});
|
||||
}
|
||||
Result<Hex8Response> Hex8Element::EvaluateLinear(const Vector& u) const {
|
||||
return Evaluate(u, false);
|
||||
}
|
||||
Result<Hex8Response> Hex8Element::EvaluateNonlinear(const Vector& u) const {
|
||||
return Evaluate(u, true);
|
||||
}
|
||||
Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
|
||||
bool nonlinear) const {
|
||||
if (u.Size() != 24U)
|
||||
return Result<Hex8Response>::Failure(
|
||||
Failure("Hex8 trial requires 24 translations."));
|
||||
const auto certificate = hex8_internal::Certify(coordinates_, &u);
|
||||
if (certificate != hex8_internal::Certificate::kPositive)
|
||||
return Result<Hex8Response>::Failure(
|
||||
Failure(CertificateMessage(certificate)));
|
||||
Hex8Response response;
|
||||
if (!nonlinear) {
|
||||
response.material_tangent = stiffness_;
|
||||
response.tangent = stiffness_;
|
||||
response.internal_force = stiffness_.Multiply(u);
|
||||
response.strain_energy = 0.5 * u.Dot(response.internal_force);
|
||||
}
|
||||
for (std::size_t g = 0; g < 8; ++g) {
|
||||
const auto& point = points_[g];
|
||||
Tensor gradient{};
|
||||
for (std::size_t a = 0; a < 8; ++a)
|
||||
for (std::size_t i = 0; i < 3; ++i)
|
||||
for (std::size_t j = 0; j < 3; ++j)
|
||||
gradient[i][j] += (u[3U * a + i] - u[i]) * point.gradients[a][j];
|
||||
Tensor f = gradient;
|
||||
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)
|
||||
for (std::size_t j = 0; j < 3; ++j) {
|
||||
if (nonlinear) {
|
||||
for (std::size_t k = 0; k < 3; ++k)
|
||||
strain[i][j] += 0.5 * f[k][i] * f[k][j];
|
||||
if (i == j) strain[i][j] -= 0.5;
|
||||
} else
|
||||
strain[i][j] = 0.5 * (gradient[i][j] + gradient[j][i]);
|
||||
}
|
||||
const double trace = strain[0][0] + strain[1][1] + strain[2][2];
|
||||
Tensor s{};
|
||||
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);
|
||||
Tensor stress = s;
|
||||
if (nonlinear) {
|
||||
const double det_f = hex8_internal::Determinant(f);
|
||||
if (!std::isfinite(det_f) || det_f <= 0.0)
|
||||
return Result<Hex8Response>::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)
|
||||
response.internal_force[3U * a + i] += Dot(p[i], ga) * point.weight;
|
||||
for (std::size_t b = 0; b < 8; ++b) {
|
||||
const auto& gb = point.gradients[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 * (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));
|
||||
if (i == k)
|
||||
response.geometric_tangent(row, column) +=
|
||||
geometric * point.weight;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if (!Finite(f) || !Finite(strain) || !Finite(s) || !Finite(stress))
|
||||
return Result<Hex8Response>::Failure(
|
||||
Failure("Hex8 trial tensor is nonfinite."));
|
||||
response.points[g] = {element_index_, static_cast<int>(g + 1U),
|
||||
point.natural, point.position,
|
||||
Pack(stress), Pack(strain)};
|
||||
}
|
||||
if (!std::isfinite(response.strain_energy))
|
||||
return Result<Hex8Response>::Failure(
|
||||
Failure("Hex8 trial energy is nonfinite."));
|
||||
for (std::size_t i = 0; i < 24; ++i) {
|
||||
if (!std::isfinite(response.internal_force[i]))
|
||||
return Result<Hex8Response>::Failure(
|
||||
Failure("Hex8 trial force is nonfinite."));
|
||||
for (std::size_t j = 0; j < 24; ++j) {
|
||||
if (nonlinear)
|
||||
response.tangent(i, j) =
|
||||
response.material_tangent(i, j) + response.geometric_tangent(i, j);
|
||||
if (!std::isfinite(response.tangent(i, j)) ||
|
||||
!std::isfinite(response.material_tangent(i, j)) ||
|
||||
!std::isfinite(response.geometric_tangent(i, j)))
|
||||
return Result<Hex8Response>::Failure(
|
||||
Failure("Hex8 trial tangent is nonfinite."));
|
||||
}
|
||||
}
|
||||
return Result<Hex8Response>::Success(std::move(response));
|
||||
}
|
||||
} // namespace fesa
|
||||
@@ -0,0 +1,239 @@
|
||||
#include "hex8_geometry.h"
|
||||
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <limits>
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
namespace fesa::hex8_internal {
|
||||
namespace {
|
||||
/// @brief Encloses exact-real arithmetic on the supplied binary64 inputs.
|
||||
struct Interval {
|
||||
double lower{0.0};
|
||||
double upper{0.0};
|
||||
};
|
||||
using Coefficients = std::array<Interval, 27>;
|
||||
|
||||
/// @brief Rounds each primitive outward by one representable binary64 value.
|
||||
/// @note /fp:precise (repository default) forbids reassociation; nextafter also
|
||||
/// encloses subnormal rounding. Overflow/nonfinite bounds fail certification.
|
||||
Interval Outward(double lower, double upper) {
|
||||
if (!std::isfinite(lower) || !std::isfinite(upper)) {
|
||||
const double invalid = std::numeric_limits<double>::quiet_NaN();
|
||||
return {invalid, invalid};
|
||||
}
|
||||
return {std::nextafter(lower, -std::numeric_limits<double>::infinity()),
|
||||
std::nextafter(upper, std::numeric_limits<double>::infinity())};
|
||||
}
|
||||
/// @brief Recognizes an exact zero, for exact identity operations only.
|
||||
bool Zero(Interval a) { return a.lower == 0.0 && a.upper == 0.0; }
|
||||
/// @brief Encloses addition, including the current X+u construction.
|
||||
Interval Add(Interval a, Interval b) {
|
||||
if (Zero(a)) return b;
|
||||
if (Zero(b)) return a;
|
||||
return Outward(a.lower + b.lower, a.upper + b.upper);
|
||||
}
|
||||
/// @brief Encloses subtraction without assuming correlated interval operands.
|
||||
Interval Subtract(Interval a, Interval b) {
|
||||
if (Zero(b)) return a;
|
||||
if (a.lower == a.upper && b.lower == b.upper && a.lower == b.lower) return {};
|
||||
return Outward(a.lower - b.upper, a.upper - b.lower);
|
||||
}
|
||||
/// @brief Encloses multiplication across all endpoint sign combinations.
|
||||
Interval Multiply(Interval a, Interval b) {
|
||||
if (Zero(a) || Zero(b)) return {};
|
||||
const std::array<double, 4> products{a.lower * b.lower, a.lower * b.upper,
|
||||
a.upper * b.lower, a.upper * b.upper};
|
||||
for (double product : products)
|
||||
if (!std::isfinite(product)) {
|
||||
const double invalid = std::numeric_limits<double>::quiet_NaN();
|
||||
return {invalid, invalid};
|
||||
}
|
||||
return Outward(*std::min_element(products.begin(), products.end()),
|
||||
*std::max_element(products.begin(), products.end()));
|
||||
}
|
||||
/// @brief Multiplies by an exact binary64 constant with outward rounding.
|
||||
Interval Scale(Interval a, double value) { return Multiply(a, {value, value}); }
|
||||
/// @brief Detects an invalid enclosure before any positivity decision.
|
||||
bool Finite(Interval a) {
|
||||
return std::isfinite(a.lower) && std::isfinite(a.upper) && a.lower <= a.upper;
|
||||
}
|
||||
/// @brief Preserves xi-fastest tensor ordering.
|
||||
std::size_t Index(std::size_t i, std::size_t j, std::size_t k) {
|
||||
return i + 3U * j + 9U * k;
|
||||
}
|
||||
/// @brief Evaluates determinant samples with every primitive enclosed.
|
||||
Interval Sample(const std::array<std::array<Interval, 3>, 8>& x,
|
||||
const Triple& natural) {
|
||||
std::array<std::array<Interval, 3>, 3> jacobian{};
|
||||
for (std::size_t a = 0; a < 8; ++a) {
|
||||
for (std::size_t d = 0; d < 3; ++d) {
|
||||
Interval derivative{0.125 * kNodeSigns[a][d], 0.125 * kNodeSigns[a][d]};
|
||||
for (std::size_t axis = 0; axis < 3; ++axis) {
|
||||
if (axis != d) {
|
||||
const auto factor =
|
||||
Add({1.0, 1.0},
|
||||
Scale({natural[axis], natural[axis]}, kNodeSigns[a][axis]));
|
||||
derivative = Multiply(derivative, factor);
|
||||
}
|
||||
}
|
||||
for (std::size_t i = 0; i < 3; ++i)
|
||||
jacobian[i][d] = Add(jacobian[i][d],
|
||||
Multiply(Subtract(x[a][i], x[0][i]), derivative));
|
||||
}
|
||||
}
|
||||
const auto& j = jacobian;
|
||||
return Add(Subtract(Multiply(j[0][0], Subtract(Multiply(j[1][1], j[2][2]),
|
||||
Multiply(j[1][2], j[2][1]))),
|
||||
Multiply(j[0][1], Subtract(Multiply(j[1][0], j[2][2]),
|
||||
Multiply(j[1][2], j[2][0])))),
|
||||
Multiply(j[0][2], Subtract(Multiply(j[1][0], j[2][1]),
|
||||
Multiply(j[1][1], j[2][0]))));
|
||||
}
|
||||
/// @brief Applies the degree-two nodal-to-Bernstein transform on one axis.
|
||||
void Transform(Coefficients& coefficients, std::size_t axis) {
|
||||
const auto original = coefficients;
|
||||
const std::size_t stride = axis == 0U ? 1U : (axis == 1U ? 3U : 9U);
|
||||
for (std::size_t k = 0; k < 3; ++k)
|
||||
for (std::size_t j = 0; j < 3; ++j)
|
||||
for (std::size_t i = 0; i < 3; ++i) {
|
||||
const std::array<std::size_t, 3> coordinate{i, j, k};
|
||||
if (coordinate[axis] != 0U) continue;
|
||||
const auto start = Index(i, j, k);
|
||||
coefficients[start + stride] = Subtract(
|
||||
Scale(original[start + stride], 2.0),
|
||||
Scale(Add(original[start], original[start + 2U * stride]), 0.5));
|
||||
}
|
||||
}
|
||||
/// @brief Uses interval de Casteljau midpoint subdivision along one axis.
|
||||
std::pair<Coefficients, Coefficients> Split(const Coefficients& original,
|
||||
std::size_t axis) {
|
||||
Coefficients left = original;
|
||||
Coefficients right = original;
|
||||
const std::size_t stride = axis == 0U ? 1U : (axis == 1U ? 3U : 9U);
|
||||
for (std::size_t k = 0; k < 3; ++k)
|
||||
for (std::size_t j = 0; j < 3; ++j)
|
||||
for (std::size_t i = 0; i < 3; ++i) {
|
||||
const std::array<std::size_t, 3> coordinate{i, j, k};
|
||||
if (coordinate[axis] != 0U) continue;
|
||||
const auto start = Index(i, j, k);
|
||||
const auto a =
|
||||
Scale(Add(original[start], original[start + stride]), 0.5);
|
||||
const auto b = Scale(
|
||||
Add(original[start + stride], original[start + 2U * stride]), 0.5);
|
||||
const auto middle = Scale(Add(a, b), 0.5);
|
||||
left[start + stride] = a;
|
||||
left[start + 2U * stride] = middle;
|
||||
right[start] = middle;
|
||||
right[start + stride] = b;
|
||||
}
|
||||
return {left, right};
|
||||
}
|
||||
} // namespace
|
||||
|
||||
Certificate Certify(const Hex8Element::Coordinates& coordinates,
|
||||
const Vector* displacement) {
|
||||
std::array<std::array<Interval, 3>, 8> x{};
|
||||
for (std::size_t a = 0; a < 8; ++a)
|
||||
for (std::size_t i = 0; i < 3; ++i) {
|
||||
const auto value = coordinates[a][i];
|
||||
if (!std::isfinite(value)) return Certificate::kNonfinite;
|
||||
x[a][i] = {value, value};
|
||||
if (displacement) {
|
||||
const auto u = (*displacement)[3U * a + i];
|
||||
if (!std::isfinite(u)) return Certificate::kNonfinite;
|
||||
x[a][i] = Add(x[a][i], {u, u});
|
||||
}
|
||||
if (!Finite(x[a][i])) return Certificate::kNonfinite;
|
||||
}
|
||||
Coefficients coefficients{};
|
||||
for (std::size_t k = 0; k < 3; ++k)
|
||||
for (std::size_t j = 0; j < 3; ++j)
|
||||
for (std::size_t i = 0; i < 3; ++i) {
|
||||
const auto sample = Sample(
|
||||
x, {static_cast<double>(i) - 1.0, static_cast<double>(j) - 1.0,
|
||||
static_cast<double>(k) - 1.0});
|
||||
if (!Finite(sample)) return Certificate::kNonfinite;
|
||||
if (sample.upper <= 0.0) return Certificate::kNonpositive;
|
||||
coefficients[Index(i, j, k)] = sample;
|
||||
}
|
||||
for (std::size_t axis = 0; axis < 3; ++axis) Transform(coefficients, axis);
|
||||
struct Cell {
|
||||
Coefficients coefficients;
|
||||
unsigned depth;
|
||||
};
|
||||
std::vector<Cell> pending{{coefficients, 0U}};
|
||||
std::size_t processed = 0U;
|
||||
while (!pending.empty()) {
|
||||
if (++processed > 65536U) return Certificate::kUncertain;
|
||||
Cell cell = std::move(pending.back());
|
||||
pending.pop_back();
|
||||
bool positive = true;
|
||||
for (const auto c : cell.coefficients) {
|
||||
if (!Finite(c)) return Certificate::kNonfinite;
|
||||
positive = positive && c.lower > 0.0;
|
||||
}
|
||||
if (positive) continue;
|
||||
// Corner coefficients are exact polynomial values enclosed by intervals.
|
||||
for (std::size_t k = 0; k < 3; k += 2)
|
||||
for (std::size_t j = 0; j < 3; j += 2)
|
||||
for (std::size_t i = 0; i < 3; i += 2)
|
||||
if (cell.coefficients[Index(i, j, k)].upper <= 0.0)
|
||||
return Certificate::kNonpositive;
|
||||
if (cell.depth >= 10U) return Certificate::kUncertain;
|
||||
const auto xi = Split(cell.coefficients, 0U);
|
||||
for (const auto& xchild : {xi.first, xi.second}) {
|
||||
const auto eta = Split(xchild, 1U);
|
||||
for (const auto& ychild : {eta.first, eta.second}) {
|
||||
const auto zeta = Split(ychild, 2U);
|
||||
pending.push_back({zeta.second, cell.depth + 1U});
|
||||
pending.push_back({zeta.first, cell.depth + 1U});
|
||||
}
|
||||
}
|
||||
}
|
||||
return Certificate::kPositive;
|
||||
}
|
||||
|
||||
void Shape(const Triple& natural, std::array<double, 8>& values,
|
||||
std::array<Triple, 8>& derivatives) {
|
||||
for (std::size_t a = 0; a < 8; ++a) {
|
||||
values[a] = 0.125;
|
||||
for (std::size_t d = 0; d < 3; ++d)
|
||||
values[a] *= 1.0 + kNodeSigns[a][d] * natural[d];
|
||||
for (std::size_t d = 0; d < 3; ++d) {
|
||||
derivatives[a][d] = 0.125 * kNodeSigns[a][d];
|
||||
for (std::size_t axis = 0; axis < 3; ++axis)
|
||||
if (axis != d)
|
||||
derivatives[a][d] *= 1.0 + kNodeSigns[a][axis] * natural[axis];
|
||||
}
|
||||
}
|
||||
}
|
||||
Tensor Jacobian(const Hex8Element::Coordinates& coordinates,
|
||||
const std::array<Triple, 8>& derivatives) {
|
||||
Tensor j{};
|
||||
for (std::size_t a = 0; a < 8; ++a)
|
||||
for (std::size_t i = 0; i < 3; ++i)
|
||||
for (std::size_t d = 0; d < 3; ++d)
|
||||
j[i][d] += (coordinates[a][i] - coordinates[0][i]) * derivatives[a][d];
|
||||
return j;
|
||||
}
|
||||
double Determinant(const Tensor& j) {
|
||||
return j[0][0] * (j[1][1] * j[2][2] - j[1][2] * j[2][1]) -
|
||||
j[0][1] * (j[1][0] * j[2][2] - j[1][2] * j[2][0]) +
|
||||
j[0][2] * (j[1][0] * j[2][1] - j[1][1] * j[2][0]);
|
||||
}
|
||||
Tensor Inverse(const Tensor& j, double determinant) {
|
||||
Tensor inverse{};
|
||||
for (std::size_t i = 0; i < 3; ++i)
|
||||
for (std::size_t d = 0; d < 3; ++d) {
|
||||
const auto r1 = (d + 1U) % 3U;
|
||||
const auto r2 = (d + 2U) % 3U;
|
||||
const auto c1 = (i + 1U) % 3U;
|
||||
const auto c2 = (i + 2U) % 3U;
|
||||
inverse[i][d] =
|
||||
(j[r1][c1] * j[r2][c2] - j[r1][c2] * j[r2][c1]) / determinant;
|
||||
}
|
||||
return inverse;
|
||||
}
|
||||
} // namespace fesa::hex8_internal
|
||||
@@ -0,0 +1,32 @@
|
||||
#ifndef FESA_ELEMENTS_HEX8_GEOMETRY_H_
|
||||
#define FESA_ELEMENTS_HEX8_GEOMETRY_H_
|
||||
|
||||
#include <array>
|
||||
|
||||
#include "fesa/elements/hex8_element.h"
|
||||
|
||||
namespace fesa::hex8_internal {
|
||||
using Triple = std::array<double, 3>;
|
||||
using Tensor = std::array<Triple, 3>;
|
||||
inline constexpr int kNodeSigns[8][3] = {{-1, -1, -1}, {1, -1, -1}, {1, 1, -1},
|
||||
{-1, 1, -1}, {-1, -1, 1}, {1, -1, 1},
|
||||
{1, 1, 1}, {-1, 1, 1}};
|
||||
/// @brief Distinguishes a witness from an exhausted conservative certificate.
|
||||
enum class Certificate { kPositive, kNonpositive, kUncertain, kNonfinite };
|
||||
/// @brief Certifies the closed natural cube with interval Bernstein bounds.
|
||||
/// @note Current addition rounding is enclosed when displacement is supplied.
|
||||
Certificate Certify(const Hex8Element::Coordinates& coordinates,
|
||||
const Vector* displacement = nullptr);
|
||||
/// @brief Evaluates trilinear shapes and natural derivatives in source order.
|
||||
void Shape(const Triple& natural, std::array<double, 8>& values,
|
||||
std::array<Triple, 8>& derivatives);
|
||||
/// @brief Computes a three-column map Jacobian using translated coordinates.
|
||||
Tensor Jacobian(const Hex8Element::Coordinates& coordinates,
|
||||
const std::array<Triple, 8>& derivatives);
|
||||
/// @brief Computes the signed determinant; no absolute value or clamp.
|
||||
double Determinant(const Tensor& matrix);
|
||||
/// @brief Computes the inverse for an already validated positive determinant.
|
||||
Tensor Inverse(const Tensor& matrix, double determinant);
|
||||
} // namespace fesa::hex8_internal
|
||||
|
||||
#endif // FESA_ELEMENTS_HEX8_GEOMETRY_H_
|
||||
@@ -356,10 +356,16 @@ Status ResultRecovery::Recover(const AnalysisModel& model,
|
||||
return results_internal::BeamResultRecovery(
|
||||
{domain.SourcePath(), 0U}, recovered.Value().source_id,
|
||||
element_index, expected_nodes, payload, candidate);
|
||||
} else {
|
||||
} else if constexpr (std::is_same_v<Payload, ShellElementResultRows>) {
|
||||
return results_internal::ShellResultRecovery(
|
||||
{domain.SourcePath(), 0U}, recovered.Value().source_id,
|
||||
element_index, payload, candidate);
|
||||
} else {
|
||||
return RecoveryFailure(
|
||||
"unsupported-solid-result-recovery",
|
||||
{domain.SourcePath(), 0U},
|
||||
recovered.Value().source_id.source_label_text,
|
||||
"Solid results require the solid procedure recovery owner.");
|
||||
}
|
||||
},
|
||||
recovered.Value().payload);
|
||||
|
||||
Reference in New Issue
Block a user