KOKO\Mimi
245 lines
9.9 KiB
C++
245 lines
9.9 KiB
C++
#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, const Vector* displacement_low) {
|
|
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 (displacement_low) {
|
|
const auto low = (*displacement_low)[3U * a + i];
|
|
if (!std::isfinite(low)) return Certificate::kNonfinite;
|
|
x[a][i] = Add(x[a][i], {low, low});
|
|
}
|
|
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
|