195 lines
9.8 KiB
C++
195 lines
9.8 KiB
C++
#pragma once
|
|
|
|
// Pure geometry keeps the physical-plane construction independently testable.
|
|
// The caller supplies the actual curved element map, never a displaced image
|
|
// of an already selected reference-plane slice.
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <limits>
|
|
#include <stdexcept>
|
|
#include <string>
|
|
#include <utility>
|
|
#include <vector>
|
|
|
|
namespace experiment::model_slices {
|
|
using Vector = std::array<double, 3>;
|
|
struct Point final { Vector reference{}, physical{}; };
|
|
struct Plane final {
|
|
int horizontal{0}, vertical{2}, normal{1};
|
|
double offset{0.};
|
|
std::string name{"xz"};
|
|
static Plane Parse(const std::string& name, double offset = 0.) {
|
|
if (!std::isfinite(offset)) throw std::invalid_argument("Slice offset must be finite.");
|
|
if (name == "xz") return {0, 2, 1, offset, name};
|
|
if (name == "xy") return {0, 1, 2, offset, name};
|
|
if (name == "yz") return {1, 2, 0, offset, name};
|
|
throw std::invalid_argument("Slice plane must be xz, xy or yz.");
|
|
}
|
|
double Distance(const Point& point) const { return point.physical[normal] - offset; }
|
|
};
|
|
using Triangle = std::array<Point, 3>;
|
|
using Segment = std::array<Point, 2>;
|
|
|
|
inline double Distance(const Vector& a, const Vector& b) {
|
|
return std::hypot(a[0] - b[0], a[1] - b[1], a[2] - b[2]);
|
|
}
|
|
inline Vector Blend(const Vector& a, const Vector& b, double fraction) {
|
|
Vector out{};
|
|
for (int i = 0; i < 3; ++i) out[i] = (1. - fraction) * a[i] + fraction * b[i];
|
|
return out;
|
|
}
|
|
inline double TwiceArea(const Triangle& t, const Plane& plane) {
|
|
const int u = plane.horizontal, v = plane.vertical;
|
|
return (t[1].physical[u] - t[0].physical[u]) * (t[2].physical[v] - t[0].physical[v]) -
|
|
(t[1].physical[v] - t[0].physical[v]) * (t[2].physical[u] - t[0].physical[u]);
|
|
}
|
|
|
|
template <class Evaluate>
|
|
Point Intersect(Point left, Point right, const Plane& plane, double tolerance, Evaluate&& evaluate) {
|
|
double fl = plane.Distance(left), fr = plane.Distance(right);
|
|
if (std::abs(fl) <= tolerance) return left;
|
|
if (std::abs(fr) <= tolerance) return right;
|
|
if ((fl < 0) == (fr < 0)) throw std::logic_error("Slice root is not bracketed.");
|
|
for (int iteration = 0; iteration < 80; ++iteration) {
|
|
// Smooth maps usually need only a few secant evaluations. Periodic
|
|
// bisection guarantees shrinking brackets even near a grazing cut.
|
|
const double fraction = iteration % 4 == 3 ? .5 : std::clamp(fl / (fl - fr), .02, .98);
|
|
const Point candidate = evaluate(Blend(left.reference, right.reference, fraction));
|
|
const double fc = plane.Distance(candidate);
|
|
if (!std::isfinite(fc)) throw std::domain_error("Nonfinite physical-plane root.");
|
|
if (std::abs(fc) <= tolerance) return candidate;
|
|
if ((fc < 0) == (fl < 0)) { left = candidate; fl = fc; }
|
|
else { right = candidate; fr = fc; }
|
|
}
|
|
throw std::runtime_error("Curved physical-plane intersection did not meet its tolerance.");
|
|
}
|
|
|
|
inline void AppendUnique(std::vector<Point>& points, const Point& point, double tolerance) {
|
|
for (const auto& old : points) if (Distance(old.physical, point.physical) <= 4 * tolerance) return;
|
|
points.push_back(point);
|
|
}
|
|
|
|
template <std::size_t Size, class Intersection>
|
|
std::vector<Point> EdgeIntersections(const std::array<Point, Size>& points, const Plane& plane,
|
|
double tolerance, Intersection&& intersect) {
|
|
std::vector<Point> cut;
|
|
for (std::size_t i = 0; i < Size; ++i) {
|
|
const double a = plane.Distance(points[i]);
|
|
if (std::abs(a) <= tolerance) AppendUnique(cut, points[i], tolerance);
|
|
for (std::size_t j = i + 1; j < Size; ++j) {
|
|
const double b = plane.Distance(points[j]);
|
|
if ((a < -tolerance && b > tolerance) || (a > tolerance && b < -tolerance))
|
|
AppendUnique(cut, intersect(points[i], points[j]), tolerance);
|
|
}
|
|
}
|
|
return cut;
|
|
}
|
|
|
|
template <class Intersection>
|
|
std::vector<Triangle> CutTetrahedron(const std::array<Point, 4>& points, const Plane& plane,
|
|
double tolerance, Intersection&& intersect) {
|
|
bool positive = false, nonpositive = false;
|
|
for (const auto& point : points) {
|
|
positive |= plane.Distance(point) > tolerance;
|
|
nonpositive |= plane.Distance(point) <= tolerance;
|
|
}
|
|
// A plane coincident with an internal face takes the trace from its
|
|
// positive-normal side. Rendering both sides would arbitrarily hide one
|
|
// of the two independent DG traces and double the slice area.
|
|
if (!positive || !nonpositive) return {};
|
|
auto cut = EdgeIntersections(points, plane, tolerance, intersect);
|
|
if (cut.size() < 3) return {};
|
|
Vector center{};
|
|
for (const auto& point : cut) for (int d = 0; d < 3; ++d) center[d] += point.physical[d] / cut.size();
|
|
std::sort(cut.begin(), cut.end(), [&](const Point& a, const Point& b) {
|
|
return std::atan2(a.physical[plane.vertical] - center[plane.vertical], a.physical[plane.horizontal] - center[plane.horizontal]) <
|
|
std::atan2(b.physical[plane.vertical] - center[plane.vertical], b.physical[plane.horizontal] - center[plane.horizontal]);
|
|
});
|
|
std::vector<Triangle> triangles;
|
|
for (std::size_t i = 1; i + 1 < cut.size(); ++i) {
|
|
Triangle triangle{cut[0], cut[i], cut[i + 1]};
|
|
if (std::abs(TwiceArea(triangle, plane)) > tolerance * tolerance) triangles.push_back(triangle);
|
|
}
|
|
return triangles;
|
|
}
|
|
|
|
template <class Intersection>
|
|
std::vector<Segment> CutSurfaceTriangle(const Triangle& points, const Plane& plane,
|
|
double tolerance, Intersection&& intersect, int& coplanar) {
|
|
int on = 0;
|
|
bool positive = false;
|
|
for (const auto& point : points) {
|
|
on += std::abs(plane.Distance(point)) <= tolerance;
|
|
positive |= plane.Distance(point) > tolerance;
|
|
}
|
|
if (on == 3) { ++coplanar; return {}; }
|
|
if (on == 2 && !positive) return {};
|
|
auto cut = EdgeIntersections(points, plane, tolerance, intersect);
|
|
if (cut.size() == 2 && Distance(cut[0].physical, cut[1].physical) > 4 * tolerance) return {{cut[0], cut[1]}};
|
|
return {};
|
|
}
|
|
|
|
// The common body diagonal makes face triangulations agree between adjacent
|
|
// reference subcubes. Field connectivity remains local to the owning element.
|
|
inline constexpr std::array<std::array<int, 4>, 6> cubeTetrahedra{{
|
|
{0, 1, 3, 7}, {0, 3, 2, 7}, {0, 2, 6, 7},
|
|
{0, 6, 4, 7}, {0, 4, 5, 7}, {0, 5, 1, 7}
|
|
}};
|
|
|
|
struct Check final { std::string name; double observed, expected, tolerance; };
|
|
inline std::vector<Check> GeometrySelfChecks() {
|
|
std::vector<Check> checks;
|
|
const auto record = [&](std::string name, double observed, double expected, double tolerance) {
|
|
if (!std::isfinite(observed) || std::abs(observed - expected) > tolerance)
|
|
throw std::runtime_error("Model slice geometry self-check failed: " + name);
|
|
checks.push_back({std::move(name), observed, expected, tolerance});
|
|
};
|
|
const auto identity = [](Vector reference) { return Point{reference, reference}; };
|
|
const Plane plane = Plane::Parse("xy", .5);
|
|
auto intersection = [&](Point a, Point b) { return Intersect(a, b, plane, 1e-13, identity); };
|
|
std::array<Point, 8> cube;
|
|
for (int i = 0; i < 8; ++i) cube[i] = identity({double(i & 1), double((i >> 1) & 1), double((i >> 2) & 1)});
|
|
double area = 0., affineIntegral = 0., planeError = 0.;
|
|
for (const auto& indices : cubeTetrahedra) {
|
|
std::array<Point, 4> tet;
|
|
for (int i = 0; i < 4; ++i) tet[i] = cube[indices[i]];
|
|
for (const auto& triangle : CutTetrahedron(tet, plane, 1e-13, intersection)) {
|
|
const double a = .5 * std::abs(TwiceArea(triangle, plane));
|
|
area += a;
|
|
for (const auto& point : triangle) {
|
|
affineIntegral += a / 3. * (point.physical[0] + 2 * point.physical[1] - point.physical[2]);
|
|
planeError = std::max(planeError, std::abs(plane.Distance(point)));
|
|
}
|
|
}
|
|
}
|
|
record("unit_cube_slice_area", area, 1., 1e-12);
|
|
record("affine_signed_field_integral", affineIntegral, 1., 1e-12);
|
|
record("physical_plane_residual", planeError, 0., 1e-12);
|
|
|
|
// This shear bends the preimage of z=0.5. Cutting z_ref=0.5 and
|
|
// displacing afterward gives the wrong result, even though it looks smooth.
|
|
const auto curved = [](Vector reference) {
|
|
Vector physical = reference;
|
|
physical[2] += .3 * reference[0] * reference[0];
|
|
return Point{reference, physical};
|
|
};
|
|
const Point root = Intersect(curved({0, 0, 0}), curved({1, 0, 1}), plane, 1e-13, curved);
|
|
const double expected = (-1. + std::sqrt(1.6)) / .6;
|
|
record("curved_edge_reference_root", root.reference[0], expected, 1e-12);
|
|
record("curved_edge_physical_plane", root.physical[2], .5, 1e-12);
|
|
int coplanar = 0;
|
|
const Triangle surface{identity({0, 0, 0}), identity({1, 0, 0}), identity({0, 0, 1})};
|
|
const auto segments = CutSurfaceTriangle(surface, plane, 1e-13, intersection, coplanar);
|
|
record("surface_intersection_segment_count", double(segments.size()), 1., 0.);
|
|
record("surface_intersection_length", Distance(segments[0][0].physical, segments[0][1].physical), .5, 1e-12);
|
|
const Plane facePlane = Plane::Parse("xy", 0.);
|
|
auto faceIntersection = [&](Point a, Point b) { return Intersect(a, b, facePlane, 1e-13, identity); };
|
|
const std::array<Point, 4> above{identity({0, 0, 0}), identity({1, 0, 0}), identity({0, 1, 0}), identity({0, 0, 1})};
|
|
auto below = above; below[3] = identity({0, 0, -1});
|
|
record("coincident_face_positive_side", double(CutTetrahedron(above, facePlane, 1e-13, faceIntersection).size()), 1., 0.);
|
|
record("coincident_face_negative_side", double(CutTetrahedron(below, facePlane, 1e-13, faceIntersection).size()), 0., 0.);
|
|
return checks;
|
|
}
|
|
} // namespace experiment::model_slices
|