tboudreaux/MeanField
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add physical-plane FEM model slice diagnostics and plots

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-19 13:46:11 -04:00
parent 7793fd9111
commit 9bac57a1df
7 changed files with 888 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

117
experiments/MODEL_SLICES.md Normal file
View File

@@ -0,0 +1,117 @@
# Physical slices of the deformed FEM model
This optional experiment displays the independent density, enthalpy and actual
material boundary of a prepared state. It changes no solver equations, fields,
orders, positivity handling or public library API.
## Use
Add `--model-slices` to a normal `plateau_diagnostic` invocation. It writes
`initial_slices` and `final_slices`, each with `xz` and `xy` subdirectories.
Optional settings:
- `--slice-plane both|xz|xy|yz` (default `both`, meaning xz and xy).
- `--slice-offset VALUE` (default zero). This is the physical coordinate normal
to the plane: y for xz, z for xy, x for yz. `both` applies the same offset to
both planes.
- `--slice-subdivisions N` (default 6, allowed 2–32). Begin with 4–8 and compare
two resolutions for a plot whose small features matter. Work grows with the
number of reference subcells; 32 is an explicit upper bound, not a suggested
default.
Existing harmonic snapshots can be observed through the checked density replay
path, using the matching binary, mesh and physical specifications:
```sh
plateau_diagnostic --index 3 --angular-momentum 0 \
--surface-prescription harmonics --surface-degree 4 \
--mesh /absolute/path/to/saved/input.smesh \
--density-snapshot /absolute/path/to/saved \
--model-slices-only --max-newton 0 --slice-subdivisions 6 \
--output /absolute/path/to/new_observation
```
`--model-slices-only` writes `frozen_slices` and the replay's normal physical
metrics/provenance, skipping the expensive density-region observations. It does
not perform Newton iterations. `--model-slices` can also accompany an existing
frozen density diagnostic. Replay checks still apply; a size match alone does
not establish compatibility of a saved state.
Render one plane with Python, NumPy and Matplotlib:
```sh
python experiments/plot_model_slices.py /absolute/path/to/new_observation/frozen_slices/xz \
--reference-outline --title 'n=3 rotating trial; unconverged' \
--output /absolute/path/to/new_figures/n3_trial_xz
```
The helper writes PNG and PDF and refuses to overwrite either. Both axes use
the same physical length scale. The solid curve is the deformed material
boundary; the optional dashed curve is the undeformed boundary intersected
with the same plane. No deformation magnification is applied. A two-sided
color scale centered at zero retains the sign when a field has both signs;
the displayed ticks are the actual values, and each color half may have a
different range. All-positive or all-negative data use a linear color scale.
No density or enthalpy values are clipped, replaced by an EOS evaluation, or
logged. Extrema in panel titles are sampled extrema.
Final-state plots use the surrounding run's convergence flag when available;
replayed or initial states are labeled as not assessed by the image. Give a
descriptive title, including the physical parameters and any known failure,
when plotting a stalled or otherwise unconverged state. A plausible picture
does not establish virial balance or an unbordered equilibrium.
## Geometry and field contract
Only stellar hexahedra are visited. Each is split into a configurable reference
lattice and consistent tetrahedra within that element. Intersections solve
`mapped_position[normal] = offset` on the actual curved map along each crossed
edge, using a safeguarded bracketed root calculation. This is a physical-plane
cut, not a deformed image of a reference-plane selection. The mapped Jacobian
must remain valid at every evaluation. Exterior compactification is not used.
Density and enthalpy are evaluated directly from their existing FEM fields at
the resolved reference coordinates. Output triangles keep their owning element
and separate vertex records. The plotter uses this supplied connectivity and
never triangulates a global point cloud; consequently it does not connect
disconnected material across vacuum or average the two DG traces on an
interface. A plane exactly coincident with an internal element face uses the
trace from its positive-normal side.
The material boundary is found from stellar/nonstellar faces and evaluated from
the stellar side. Its curve is computed separately from the volume coloring.
The undeformed outline is independently intersected with the physical plane;
projecting reference positions of the deformed cut would be a different curve.
An exactly coplanar surface patch is reported as unsupported instead of being
silently drawn as a unique intersection curve.
Curved cuts are displayed using straight triangles and boundary chords. Field
colors are piecewise linear on those triangles, not exact values throughout
their interiors. Root tolerances, maximum physical-plane residuals, sampled
centroid/midpoint chord errors, counts, bounds and time are recorded in each
plane's `metadata.csv`. The chord errors and bounds are sampled diagnostics,
not rigorous global enclosures. Tangencies or multiple crossings contained
between subdivision vertices can be missed; compare subdivisions for image
stability. This visualization does not certify pointwise positivity or exclude
small unsampled features.
## Checks
Every export runs pure geometric controls: unit-cube section area, integration
of a signed affine field, a curved shear for which the reference-plane cut is
wrong, surface segment intersection, and deterministic one-sided treatment of
a coincident internal face. The same controls can run without FEM:
```sh
c++ -std=c++20 -O0 experiments/model_slice_geometry_tests.cpp -o /tmp/model_slice_checks
/tmp/model_slice_checks
python -m unittest discover -s experiments -p test_plot_model_slices.py -v
```
The Python controls cover signed independent fields, supplied connectivity for
disconnected material, equal axes, PNG/PDF rendering, and rejection of a
cross-element triangle, an off-plane field/outline, and shape exaggeration.
They are geometry and plotting controls, not a substitute for a coupled-model
validation run. The exporter checks that its density, enthalpy and displacement
source buffers remain bitwise unchanged.

258
experiments/model_slice_diagnostics.hpp Normal file
View File

@@ -0,0 +1,258 @@
#pragma once
#include "model_slice_geometry.hpp"
#include "polytrope_physical_state.hpp"
#include <chrono>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <map>
namespace experiment::model_slices {
struct Options final {
std::string plane{"both"};
double offset{0.};
int subdivisions{6};
void Validate() const {
if (subdivisions < 2 || subdivisions > 32)
throw std::invalid_argument("Slice subdivisions must be between 2 and 32.");
(void)Plane::Parse(plane == "both" ? "xz" : plane, offset);
}
};
namespace detail {
inline std::ofstream Open(const std::filesystem::path& path) {
if (std::filesystem::exists(path)) throw std::invalid_argument("Refusing to overwrite slice output: " + path.string());
std::ofstream stream(path);
stream.exceptions(std::ios::failbit | std::ios::badbit);
stream << std::setprecision(17);
return stream;
}
inline bool Identical(const mfem::Vector& a, const mfem::Vector& b) {
return a.Size() == b.Size() && (a.Size() == 0 ||
std::memcmp(a.GetData(), b.GetData(), std::size_t(a.Size()) * sizeof(double)) == 0);
}
inline mfem::IntegrationPoint IntegrationPoint(const Vector& reference) {
mfem::IntegrationPoint point;
point.Set3(reference[0], reference[1], reference[2]);
return point;
}
class ElementMap final {
public:
ElementMap(polytrope_validation::PhysicalState& state, int element, bool deformed = true)
: state(state), deformed(deformed), transformation(*state.finiteElements.mesh->GetElementTransformation(element)) {}
Point operator()(const Vector& reference) {
const auto point = IntegrationPoint(reference);
Vector physical{};
if (deformed) {
if (state.mapping.EvaluatePoint(transformation, point, mapped) != mean_field::mapping::MappingStatus::valid || mapped.compactified)
throw std::domain_error("Model slices require valid, noncompactified stellar geometry.");
for (int d = 0; d < 3; ++d) physical[d] = mapped.physical_position(d);
} else {
transformation.Transform(point, referencePosition);
for (int d = 0; d < 3; ++d) physical[d] = referencePosition(d);
}
for (double value : physical) if (!std::isfinite(value)) throw std::domain_error("Nonfinite slice geometry.");
return {reference, physical};
}
private:
polytrope_validation::PhysicalState& state;
bool deformed;
mfem::ElementTransformation& transformation;
mean_field::mapping::MappingPointContext mapped;
mfem::Vector referencePosition = mfem::Vector(3);
};
class Intersections final {
public:
Intersections(ElementMap& map, Plane plane, double tolerance) : map(map), plane(std::move(plane)), tolerance(tolerance) {}
Point operator()(const Point& a, const Point& b) {
const Vector& first = a.reference < b.reference ? a.reference : b.reference;
const Vector& last = a.reference < b.reference ? b.reference : a.reference;
const std::array<double, 6> key{first[0], first[1], first[2], last[0], last[1], last[2]};
if (const auto found = cache.find(key); found != cache.end()) return found->second;
const Point root = Intersect(a, b, plane, tolerance, map);
cache.emplace(key, root);
return root;
}
private:
ElementMap& map;
Plane plane;
double tolerance;
std::map<std::array<double, 6>, Point> cache;
};
struct Bounds final {
Vector lower{INFINITY, INFINITY, INFINITY}, upper{-INFINITY, -INFINITY, -INFINITY};
void Add(const Vector& value) {
for (int d = 0; d < 3; ++d) { lower[d] = std::min(lower[d], value[d]); upper[d] = std::max(upper[d], value[d]); }
}
};
inline void ExportPlane(polytrope_validation::PhysicalState& state, const Options& options, const Plane& plane,
const std::filesystem::path& directory) {
const auto start = std::chrono::steady_clock::now();
if (!std::filesystem::create_directory(directory)) throw std::invalid_argument("Slice directory already exists.");
auto vertices = Open(directory / "vertices.csv"), triangles = Open(directory / "triangles.csv");
vertices << "vertex,element,attribute,xi,eta,zeta,x,y,z,rho,h,plane_error\n";
triangles << "triangle,element,v0,v1,v2\n";
auto& mesh = *state.finiteElements.mesh;
const int n = options.subdivisions, width = n + 1;
const auto index = [width](int i, int j, int k) { return i + width * (j + width * k); };
Bounds bounds;
std::size_t vertexCount = 0, triangleCount = 0;
double maximumPlaneError = 0., maximumChordError = 0., minimumDensity = INFINITY, minimumEnthalpy = INFINITY;
// Coordinate-scaled tolerance is independent of density units and order.
double scale = 0.;
for (int element = 0; element < mesh.GetNE(); ++element) if (state.isStellar(element)) {
ElementMap map(state, element);
for (int corner = 0; corner < 8; ++corner) {
const auto point = map({double(corner & 1), double((corner >> 1) & 1), double((corner >> 2) & 1)});
scale = std::max(scale, std::hypot(point.physical[0], point.physical[1], point.physical[2]));
}
}
if (!(scale > 0) || !std::isfinite(scale)) throw std::domain_error("Invalid stellar slice scale.");
const double tolerance = 2e-12 * std::max(scale, std::abs(plane.offset));
for (int element = 0; element < mesh.GetNE(); ++element) if (state.isStellar(element)) {
if (mesh.GetElementBaseGeometry(element) != mfem::Geometry::CUBE)
throw std::invalid_argument("Model slices currently require hexahedral stellar elements.");
ElementMap map(state, element);
Intersections intersections(map, plane, tolerance);
std::vector<Point> lattice(static_cast<std::size_t>(width * width * width));
for (int k = 0; k <= n; ++k) for (int j = 0; j <= n; ++j) for (int i = 0; i <= n; ++i) {
auto& point = lattice[index(i, j, k)];
point = map({double(i) / n, double(j) / n, double(k) / n});
bounds.Add(point.physical);
}
for (int k = 0; k < n; ++k) for (int j = 0; j < n; ++j) for (int i = 0; i < n; ++i) {
std::array<Point, 8> cube;
for (int c = 0; c < 8; ++c) cube[c] = lattice[index(i + (c & 1), j + ((c >> 1) & 1), k + ((c >> 2) & 1))];
for (const auto& indices : cubeTetrahedra) {
std::array<Point, 4> tetrahedron;
for (int c = 0; c < 4; ++c) tetrahedron[c] = cube[indices[c]];
for (const auto& triangle : CutTetrahedron(tetrahedron, plane, tolerance, intersections)) {
triangles << triangleCount++ << ',' << element << ',' << vertexCount << ',' << vertexCount + 1 << ',' << vertexCount + 2 << '\n';
Vector referenceCenter{}, physicalCenter{};
for (const auto& point : triangle) {
const auto ip = IntegrationPoint(point.reference);
const double rho = state.density.GetValue(element, ip), h = state.enthalpy.GetValue(element, ip);
if (!std::isfinite(rho) || !std::isfinite(h)) throw std::domain_error("Nonfinite slice field value.");
const double error = std::abs(plane.Distance(point));
if (error > tolerance) throw std::logic_error("A slice vertex is off its physical plane.");
maximumPlaneError = std::max(maximumPlaneError, error);
minimumDensity = std::min(minimumDensity, rho); minimumEnthalpy = std::min(minimumEnthalpy, h);
vertices << vertexCount++ << ',' << element << ',' << mesh.GetAttribute(element);
for (double value : point.reference) vertices << ',' << value;
for (double value : point.physical) vertices << ',' << value;
vertices << ',' << rho << ',' << h << ',' << error << '\n';
for (int d = 0; d < 3; ++d) {
referenceCenter[d] += point.reference[d] / 3.; physicalCenter[d] += point.physical[d] / 3.;
}
bounds.Add(point.physical);
}
maximumChordError = std::max(maximumChordError, Distance(map(referenceCenter).physical, physicalCenter));
}
}
}
}
vertices.flush(); triangles.flush();
std::array<std::size_t, 2> segmentCounts{};
std::array<int, 2> coplanar{};
std::array<double, 2> surfaceError{}, surfaceChordError{};
for (int reference = 0; reference < 2; ++reference) {
auto surface = Open(directory / (reference ? "reference_surface.csv" : "surface.csv"));
surface << "segment,face,element,x0,y0,z0,x1,y1,z1,maximum_plane_error\n";
for (int face = 0; face < mesh.GetNumFaces(); ++face) {
int left = -1, right = -1;
mesh.GetFaceElements(face, &left, &right);
const bool leftStar = left >= 0 && state.isStellar(left), rightStar = right >= 0 && state.isStellar(right);
if (leftStar == rightStar) continue;
const int element = leftStar ? left : right;
auto* faceMap = mesh.GetFaceElementTransformations(face);
if (faceMap->GetGeometryType() != mfem::Geometry::SQUARE)
throw std::invalid_argument("Model slices currently require quadrilateral material faces.");
ElementMap map(state, element, !reference);
Intersections intersections(map, plane, tolerance);
std::vector<Point> lattice(static_cast<std::size_t>(width * width));
for (int j = 0; j <= n; ++j) for (int i = 0; i <= n; ++i) {
mfem::IntegrationPoint facePoint, elementPoint;
facePoint.Set2(double(i) / n, double(j) / n);
if (leftStar) faceMap->Loc1.Transform(facePoint, elementPoint); else faceMap->Loc2.Transform(facePoint, elementPoint);
lattice[i + width * j] = map({elementPoint.x, elementPoint.y, elementPoint.z});
}
for (int j = 0; j < n; ++j) for (int i = 0; i < n; ++i) {
const std::array<Point, 4> square{lattice[i + width * j], lattice[i + 1 + width * j],
lattice[i + width * (j + 1)], lattice[i + 1 + width * (j + 1)]};
for (const auto& indices : {std::array<int, 3>{0, 1, 3}, std::array<int, 3>{0, 3, 2}}) {
const Triangle triangle{square[indices[0]], square[indices[1]], square[indices[2]]};
for (const auto& segment : CutSurfaceTriangle(triangle, plane, tolerance, intersections, coplanar[reference])) {
const double error = std::max(std::abs(plane.Distance(segment[0])), std::abs(plane.Distance(segment[1])));
if (error > tolerance) throw std::logic_error("A surface intersection is off its physical plane.");
surfaceError[reference] = std::max(surfaceError[reference], error);
surfaceChordError[reference] = std::max(surfaceChordError[reference],
Distance(map(Blend(segment[0].reference, segment[1].reference, .5)).physical,
Blend(segment[0].physical, segment[1].physical, .5)));
surface << segmentCounts[reference]++ << ',' << face << ',' << element;
for (const auto& point : segment) { for (double value : point.physical) surface << ',' << value; bounds.Add(point.physical); }
surface << ',' << error << '\n';
}
}
}
}
}
auto metadata = Open(directory / "metadata.csv");
metadata << "metric,value\nplane," << plane.name << "\nnormal_axis," << plane.normal
<< "\nhorizontal_axis," << plane.horizontal << "\nvertical_axis," << plane.vertical
<< "\noffset," << plane.offset << "\nsubdivisions," << n << "\nroot_tolerance," << tolerance
<< "\ndeformation_scale,1\ntriangle_count," << triangleCount << "\nvertex_count," << vertexCount
<< "\nmaximum_plane_error," << maximumPlaneError << "\nmaximum_sampled_chord_error," << maximumChordError
<< "\nminimum_sampled_density," << minimumDensity << "\nminimum_sampled_enthalpy," << minimumEnthalpy
<< "\nsurface_segments," << segmentCounts[0] << "\nreference_surface_segments," << segmentCounts[1]
<< "\nsurface_coplanar_triangles," << coplanar[0] << "\nreference_surface_coplanar_triangles," << coplanar[1]
<< "\nsurface_maximum_plane_error," << surfaceError[0] << "\nreference_surface_maximum_plane_error," << surfaceError[1]
<< "\nsurface_maximum_sampled_chord_error," << surfaceChordError[0]
<< "\nreference_surface_maximum_sampled_chord_error," << surfaceChordError[1]
<< "\nu_min," << bounds.lower[plane.horizontal] << "\nu_max," << bounds.upper[plane.horizontal]
<< "\nv_min," << bounds.lower[plane.vertical] << "\nv_max," << bounds.upper[plane.vertical]
<< "\nseconds," << std::chrono::duration<double>(std::chrono::steady_clock::now() - start).count() << '\n';
if (triangleCount && !segmentCounts[0]) throw std::runtime_error("Material slice has no resolved outer boundary; refine or move the plane.");
if (coplanar[0] || coplanar[1]) throw std::runtime_error("Slice is coplanar with a surface patch; choose a different plane offset.");
}
} // namespace detail
inline void Record(polytrope_validation::PhysicalState& state, const Options& options,
const std::filesystem::path& output, const std::string& label) {
options.Validate();
int ranks = 0;
MPI_Comm_size(state.finiteElements.mesh->GetComm(), &ranks);
if (ranks != 1 || state.finiteElements.mesh->Dimension() != 3 || !state.finiteElements.mesh->Conforming())
throw std::invalid_argument("Model slices require a conforming 3D mesh on one MPI rank.");
const mfem::Vector density(state.density), enthalpy(state.enthalpy), displacement(state.displacement);
const auto directory = output / (label + "_slices");
if (!std::filesystem::create_directory(directory)) throw std::invalid_argument("Slice output directory already exists.");
auto checks = detail::Open(directory / "geometry_self_checks.csv");
checks << "case,observed,expected,tolerance,passed\n";
for (const auto& check : GeometrySelfChecks())
checks << check.name << ',' << check.observed << ',' << check.expected << ',' << check.tolerance << ",1\n";
const std::vector<std::string> planes = options.plane == "both" ? std::vector<std::string>{"xz", "xy"} : std::vector<std::string>{options.plane};
for (const auto& name : planes) {
std::cout << label << " physical " << name << " slice, subdivisions=" << options.subdivisions << std::endl;
detail::ExportPlane(state, options, Plane::Parse(name, options.offset), directory / name);
}
if (!detail::Identical(density, state.density) || !detail::Identical(enthalpy, state.enthalpy) || !detail::Identical(displacement, state.displacement))
throw std::logic_error("Model slice export changed a source field.");
auto scope = detail::Open(directory / "scope.txt");
scope << "source_fields_bitwise_identical=1\ndeformation_scale=1\n"
<< "Vertices evaluate the independent FEM density and enthalpy on the actual deformed physical plane.\n"
<< "Triangles remain local to their material element; exterior fields are never sampled.\n"
<< "A coincident internal face uses the positive-normal-side DG trace.\n"
<< "Surface and undeformed reference outlines are separately intersected with the same physical plane.\n"
<< "Interior colors and boundary chords approximate curved cuts at the recorded subdivision count.\n"
<< "Sampled chord errors and bounds are diagnostics, not rigorous geometric enclosures.\n"
<< "Refine subdivisions to assess image stability; subcell tangencies or multiple edge crossings can be missed.\n"
<< "The slice does not certify positivity, equilibrium, virial balance or the absence of unsampled features.\n";
}
} // namespace experiment::model_slices

194
experiments/model_slice_geometry.hpp Normal file
View File

@@ -0,0 +1,194 @@
#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

13
experiments/model_slice_geometry_tests.cpp Normal file
View File

@@ -0,0 +1,13 @@
// Standalone, no MFEM or MPI setup: c++ -std=c++20 -O0 this_file.cpp -o /tmp/slice_checks
#include "model_slice_geometry.hpp"
#include <iostream>
int main() {
try {
for (const auto& check : experiment::model_slices::GeometrySelfChecks())
std::cout << check.name << ": " << check.observed << '\n';
} catch (const std::exception& error) {
std::cerr << error.what() << '\n';
return 1;
}
}

34
experiments/plateau_diagnostic.cpp
View File

@@ -38,6 +38,7 @@ import mean_field;
#include "density_order_projection.hpp"
#include "continuation_snapshot.hpp"
#include "harmonic_linear_trials.hpp"
#include "model_slice_diagnostics.hpp"
using namespace mean_field;
namespace {
@@ -53,6 +54,8 @@ struct Options {
bool densityDiagnostic = false;
bool densityRegions = false;
bool densityOrderProjection = false, densityOrderFinal = false;
bool modelSlices = false, modelSlicesOnly = false;
experiment::model_slices::Options slices;
int surfaceDegree = 4, surfaceMaximumOrder = -1;
std::string surfacePrescription = "nodal";
bool coupledTrialProbe = false, coupledSeedControl = false;
@@ -154,6 +157,17 @@ Options Parse(int argc, char **argv) {
o.allowUnconvergedRestart = true;
else if (a == "--benchmark-trials")
o.benchmarkTrials = true;
else if (a == "--model-slices")
o.modelSlices = true;
else if (a == "--model-slices-only") {
o.modelSlices = o.modelSlicesOnly = o.densityDiagnostic = true;
}
else if (a == "--slice-plane")
o.slices.plane = v();
else if (a == "--slice-offset")
o.slices.offset = std::stod(v());
else if (a == "--slice-subdivisions")
o.slices.subdivisions = std::stoi(v());
else if (a == "--density-snapshot")
o.densitySnapshot = v();
else if (a == "--surface-prescription")
@@ -231,6 +245,14 @@ Options Parse(int argc, char **argv) {
throw std::invalid_argument("Snapshot restart requires native harmonic coordinates and no legacy replay");
if (o.benchmarkTrials && !o.preconditionerProbe)
throw std::invalid_argument("--benchmark-trials requires --preconditioner-probe");
o.slices.Validate();
if (o.modelSlicesOnly && (o.densityOrderProjection || o.densityRegions || o.densityOrderFinal))
throw std::invalid_argument("--model-slices-only cannot request another density observation");
if (o.modelSlices && (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe ||
o.positivityProbe || o.positivityCertificate || o.preconditionerProbe || o.nonlinearStepProbe ||
o.reducedSurfaceProbe || o.reducedSurfaceAlgebraCheck || o.weakShapeProbe || o.coupledTrialProbe ||
o.dampedCoupledProbe || !o.probeState.empty() || !o.probeDirection.empty() || !o.replayCandidate.empty()))
throw std::invalid_argument("Model slices require live initial/final diagnostics or density-snapshot replay");
if (o.surfacePrescription != "nodal" && o.surfacePrescription != "harmonics")
throw std::invalid_argument("Surface prescriptions are nodal or harmonics");
if (!o.densitySnapshot.empty() && !o.densityDiagnostic)
@@ -755,6 +777,9 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
<< "\ndisplacement_order=" << field::Displacement::Vector::familyOrder
<< "\ndensity_order_projection=" << o.densityOrderProjection
<< "\ndensity_order_final=" << o.densityOrderFinal << '\n';
metadata << "model_slices=" << o.modelSlices << "\nmodel_slices_only=" << o.modelSlicesOnly
<< "\nslice_plane=" << o.slices.plane << "\nslice_offset=" << o.slices.offset
<< "\nslice_subdivisions=" << o.slices.subdivisions << '\n';
if (o.densityDiagnostic) {
metadata << "density_diagnostic=1\ndensity_region_integral=" << o.densityRegions << "\ndensity_source="
<< (o.densitySnapshot.empty() ? "lane_emden_seed" : o.densitySnapshot.string()) << '\n';
@@ -795,7 +820,12 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
adapter.scatter(closureResidual, closureFull);
SaveVector(o.output / "production_closure_residual_full.txt", closureFull);
SaveVector(o.output / "observed_density_coefficients.txt", physical.density);
if (o.densityOrderProjection)
if (o.modelSlices)
experiment::model_slices::Record(physical, o.slices, o.output, "frozen");
if (o.modelSlicesOnly) {
// Replay keeps its ordinary physical metrics and provenance,
// but a plot request need not integrate all negative regions.
} else if (o.densityOrderProjection)
experiment::plateau_diagnostics::RecordDensityOrderProjection(physical, o.n, scales.K, o.output);
else if (o.densityRegions)
experiment::plateau_diagnostics::RecordDensityNegativeRegions(physical, o.n, scales.K, o.output);
@@ -997,6 +1027,8 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
std::cout << label << " physical diagnostics" << std::endl;
experiment::polytrope_validation::PhysicalState physical(s, f);
auto out = File(o.output / (label + "_physical.csv"));
if (o.modelSlices)
experiment::model_slices::Record(physical, o.slices, o.output, label);
out << "quadrature,metric,value\n";
if (o.surfacePrescription == "harmonics") {
deformation::SphericalHarmonicSurfaceOptions surfaceOptions{.maximumDegree = o.surfaceDegree};

170
experiments/plot_model_slices.py Normal file
View File

@@ -0,0 +1,170 @@
#!/usr/bin/env python3
"""Plot physical-plane FEM cuts without joining material across the vacuum."""
from __future__ import annotations
import argparse
import csv
from dataclasses import dataclass
from pathlib import Path
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from matplotlib.colors import Normalize, TwoSlopeNorm
from matplotlib.tri import Triangulation
import numpy as np
@dataclass
class Slice:
directory: Path
metadata: dict[str, str]
vertices: np.ndarray
triangles: np.ndarray
surface: np.ndarray
reference_surface: np.ndarray
def table(path: Path, names: list[str]) -> np.ndarray:
with path.open(newline="") as stream:
reader = csv.DictReader(stream)
if reader.fieldnames is None or any(name not in reader.fieldnames for name in names):
raise ValueError(f"Missing required columns in {path}")
rows = [[float(row[name]) for name in names] for row in reader]
values = np.asarray(rows, dtype=float).reshape((-1, len(names)))
if not np.isfinite(values).all():
raise ValueError(f"Nonfinite coordinates or fields in {path}")
return values
def load_slice(directory: Path) -> Slice:
with (directory / "metadata.csv").open(newline="") as stream:
metadata = {row["metric"]: row["value"] for row in csv.DictReader(stream)}
axes = {"xz": (0, 2, 1), "xy": (0, 1, 2), "yz": (1, 2, 0)}
plane = metadata["plane"]
if plane not in axes or tuple(int(metadata[name]) for name in ("horizontal_axis", "vertical_axis", "normal_axis")) != axes[plane]:
raise ValueError("Slice orientation metadata is inconsistent")
if float(metadata["deformation_scale"]) != 1:
raise ValueError("This plotter requires actual deformation with scale exactly one")
tolerance, offset = float(metadata["root_tolerance"]), float(metadata["offset"])
if not np.isfinite([tolerance, offset]).all() or tolerance <= 0:
raise ValueError("Invalid physical-plane tolerance or offset")
v = table(directory / "vertices.csv", ["vertex", "element", "x", "y", "z", "rho", "h", "plane_error"])
t = table(directory / "triangles.csv", ["triangle", "element", "v0", "v1", "v2"])
if len(v) != int(metadata["vertex_count"]) or len(t) != int(metadata["triangle_count"]):
raise ValueError("Slice row counts disagree with metadata")
if not np.array_equal(v[:, 0], np.arange(len(v))) or not np.array_equal(t[:, 0], np.arange(len(t))):
raise ValueError("Slice indices must be contiguous")
connectivity = t[:, 2:].astype(int)
if not np.array_equal(connectivity, t[:, 2:]) or (connectivity.size and (connectivity.min() < 0 or connectivity.max() >= len(v))):
raise ValueError("Invalid supplied triangle connectivity")
if len(t) and not np.all(v[connectivity, 1] == t[:, 1, None]):
raise ValueError("A triangle joins different material elements")
normal = axes[plane][2]
if len(v) and (np.max(np.abs(v[:, 2 + normal] - offset)) > tolerance * 1.01 or np.max(v[:, 7]) > tolerance * 1.01):
raise ValueError("Field samples do not lie on the requested physical plane")
surfaces = []
for filename, count in (("surface.csv", "surface_segments"), ("reference_surface.csv", "reference_surface_segments")):
s = table(directory / filename, ["x0", "y0", "z0", "x1", "y1", "z1", "maximum_plane_error"])
if len(s) != int(metadata[count]):
raise ValueError("Surface row count disagrees with metadata")
if len(s) and (np.max(np.abs(s[:, [normal, 3 + normal]] - offset)) > tolerance * 1.01 or np.max(s[:, 6]) > tolerance * 1.01):
raise ValueError("An outline does not intersect the requested physical plane")
surfaces.append(s)
if int(metadata["surface_coplanar_triangles"]) or int(metadata["reference_surface_coplanar_triangles"]):
raise ValueError("A surface patch is coplanar; this is not a resolved curve intersection")
if len(t) and not len(surfaces[0]):
raise ValueError("A nonempty material slice needs its deformed boundary")
return Slice(directory, metadata, v, connectivity, *surfaces)
def color_scale(values: np.ndarray):
"""Keep the sign visible; neither abs, logarithms nor clipping changes data."""
lower, upper = float(values.min()), float(values.max())
if lower < 0 < upper:
return TwoSlopeNorm(vmin=lower, vcenter=0, vmax=upper), "coolwarm", "color scale centered at zero"
if lower == upper:
delta = max(abs(lower) * .01, 1e-12)
return Normalize(lower - delta, upper + delta), "viridis", "linear color scale"
return Normalize(lower, upper), "viridis", "linear color scale"
def render(data: Slice, *, reference_outline: bool = False, title: str | None = None):
m, vertices = data.metadata, data.vertices
u, v = int(m["horizontal_axis"]), int(m["vertical_axis"])
labels = ("x", "y", "z")
figure, axes = plt.subplots(1, 2, figsize=(12.4, 5.6), layout="constrained")
connectivity = None
if len(data.triangles):
# Providing connectivity is essential: Delaunay triangulation of this
# point cloud would bridge holes and average across DG interfaces.
connectivity = Triangulation(vertices[:, 2 + u], vertices[:, 2 + v], triangles=data.triangles)
for ax, column, name in zip(axes, (5, 6), ("Independent density", "Specific enthalpy")):
if connectivity is not None:
values = vertices[:, column]
norm, cmap, note = color_scale(values)
artist = ax.tripcolor(connectivity, values, shading="gouraud", norm=norm, cmap=cmap, rasterized=True)
bar = figure.colorbar(artist, ax=ax, shrink=.82, label="model field units")
ticks = bar.get_ticks()
ticks = ticks[(ticks >= norm.vmin) & (ticks <= norm.vmax)]
bar.set_ticks(np.unique(np.r_[norm.vmin, ticks, norm.vmax]))
ax.set_title(f"{name}\nsampled min {values.min():.5g}; max {values.max():.5g}\n{note}", fontsize=10)
else:
ax.set_title(name)
ax.text(.5, .5, "No material intersection resolved", ha="center", va="center", transform=ax.transAxes)
for surface, color, linestyle, label in (
(data.surface, "#111111", "solid", "Deformed material boundary"),
(data.reference_surface if reference_outline else np.empty((0, 7)), "#656565", "dashed", "Undeformed boundary"),
):
if len(surface):
segments = np.stack((surface[:, [u, v]], surface[:, [3 + u, 3 + v]]), axis=1)
ax.add_collection(LineCollection(segments, colors=color, linewidths=1.05, linestyles=linestyle, label=label))
lower = np.array([float(m["u_min"]), float(m["v_min"])])
upper = np.array([float(m["u_max"]), float(m["v_max"])])
if not np.isfinite(np.r_[lower, upper]).all() or np.any(upper <= lower):
raise ValueError("Invalid slice display bounds")
padding = .04 * max(upper - lower)
ax.set_xlim(lower[0] - padding, upper[0] + padding)
ax.set_ylim(lower[1] - padding, upper[1] + padding)
ax.set_aspect("equal", adjustable="box")
ax.set_xlabel(f"{labels[u]} (model length units)")
ax.set_ylabel(f"{labels[v]} (model length units)")
ax.set_facecolor("#f3f3f3")
if len(data.surface) or (reference_outline and len(data.reference_surface)):
ax.legend(loc="upper right", fontsize=7, framealpha=.9)
description = title or data.directory.parent.name.replace("_", " ")
status = "Equilibrium not assessed by this image"
run_metadata = data.directory.parent.parent / "metadata.txt"
if data.directory.parent.name == "final_slices" and run_metadata.exists():
entries = dict(line.split("=", 1) for line in run_metadata.read_text().splitlines() if "=" in line)
if entries.get("converged") == "0":
status = "UNCONVERGED state"
elif entries.get("converged") == "1":
status = "Solver reported convergence; physical validation is separate"
figure.suptitle(f"{description} — physical {m['plane']} plane, {labels[int(m['normal_axis'])]}={float(m['offset']):.6g}\n"
f"{status}; actual deformation; subdivisions {m['subdivisions']}", fontsize=11)
return figure
def main() -> None:
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("directory", type=Path, help="One plane directory, e.g. final_slices/xz")
parser.add_argument("--output", type=Path, help="Output prefix; defaults to DIRECTORY/model_slice")
parser.add_argument("--reference-outline", action="store_true")
parser.add_argument("--title")
args = parser.parse_args()
output = args.output or args.directory / "model_slice"
paths = [output.with_suffix(suffix) for suffix in (".png", ".pdf")]
if any(path.exists() for path in paths):
raise FileExistsError("Refusing to overwrite an existing slice figure")
figure = render(load_slice(args.directory), reference_outline=args.reference_outline, title=args.title)
output.parent.mkdir(parents=True, exist_ok=True)
for path in paths:
figure.savefig(path, dpi=200)
print(path)
plt.close(figure)
if __name__ == "__main__":
main()

103
experiments/test_plot_model_slices.py Normal file
View File

@@ -0,0 +1,103 @@
#!/usr/bin/env python3
"""Small schema/plot controls; these are not stellar-equilibrium validation."""
import csv
from pathlib import Path
import tempfile
import unittest
import matplotlib.pyplot as plt
import numpy as np
import plot_model_slices as slices
def write_csv(path, header, rows):
with path.open("w", newline="") as stream:
writer = csv.writer(stream)
writer.writerow(header)
writer.writerows(rows)
def fixture(directory):
# Two disconnected material triangles leave a genuine empty gap. Their
# enthalpy is independent of density and one density sample is negative.
coordinates = [(-2, 0), (-1, 0), (-1.5, 1), (1, 0), (2, 0), (1.5, 1)]
rho, h = [-.4, .2, 1, .5, 1, 2], [1, 2, 3, 4, 5, 6]
vertices = [[i, i // 3, u, .2, v, rho[i], h[i], 0] for i, (u, v) in enumerate(coordinates)]
write_csv(directory / "vertices.csv", ["vertex", "element", "x", "y", "z", "rho", "h", "plane_error"], vertices)
write_csv(directory / "triangles.csv", ["triangle", "element", "v0", "v1", "v2"], [[0, 0, 0, 1, 2], [1, 1, 3, 4, 5]])
segments = []
for element in range(2):
for i, j in ((0, 1), (1, 2), (2, 0)):
a, b = coordinates[3 * element + i], coordinates[3 * element + j]
segments.append([a[0], .2, a[1], b[0], .2, b[1], 0])
for name in ("surface.csv", "reference_surface.csv"):
write_csv(directory / name, ["x0", "y0", "z0", "x1", "y1", "z1", "maximum_plane_error"], segments)
metadata = dict(plane="xz", horizontal_axis=0, vertical_axis=2, normal_axis=1, offset=.2,
root_tolerance=1e-11, deformation_scale=1, vertex_count=6, triangle_count=2,
surface_segments=6, reference_surface_segments=6, surface_coplanar_triangles=0,
reference_surface_coplanar_triangles=0, subdivisions=6, u_min=-2, u_max=2, v_min=0, v_max=1)
write_csv(directory / "metadata.csv", ["metric", "value"], metadata.items())
class ModelSlicePlotTests(unittest.TestCase):
def setUp(self):
self.temporary = tempfile.TemporaryDirectory()
self.directory = Path(self.temporary.name)
fixture(self.directory)
def tearDown(self):
plt.close("all")
self.temporary.cleanup()
def test_signed_independent_fields_are_preserved(self):
data = slices.load_slice(self.directory)
self.assertEqual(data.vertices[0, 5], -.4)
self.assertEqual(data.vertices[0, 6], 1.)
norm, _, _ = slices.color_scale(data.vertices[:, 5])
self.assertEqual(norm(0.), .5)
self.assertEqual(norm(-.4), 0.)
def test_no_triangulation_across_vacuum(self):
data = slices.load_slice(self.directory)
np.testing.assert_array_equal(data.triangles, [[0, 1, 2], [3, 4, 5]])
for triangle in data.triangles:
horizontal = data.vertices[triangle, 2]
self.assertTrue(np.all(horizontal < 0) or np.all(horizontal > 0))
figure = slices.render(data, reference_outline=True, title="Analytic geometry/schema control")
self.assertEqual(figure.axes[0].get_aspect(), 1.)
self.assertEqual(figure.axes[1].get_aspect(), 1.)
for suffix in ("png", "pdf"):
path = self.directory / ("control." + suffix)
figure.savefig(path)
self.assertGreater(path.stat().st_size, 1000)
def test_cross_element_triangle_is_rejected(self):
write_csv(self.directory / "triangles.csv", ["triangle", "element", "v0", "v1", "v2"], [[0, 0, 0, 1, 4], [1, 1, 3, 4, 5]])
with self.assertRaisesRegex(ValueError, "different material elements"):
slices.load_slice(self.directory)
def test_wrong_physical_plane_is_rejected(self):
path = self.directory / "vertices.csv"
text = path.read_text().replace("0,0,-2,0.2,0,", "0,0,-2,0.25,0,")
path.write_text(text)
with self.assertRaisesRegex(ValueError, "physical plane"):
slices.load_slice(self.directory)
def test_wrong_reference_plane_is_rejected(self):
path = self.directory / "reference_surface.csv"
values = slices.table(path, ["x0", "y0", "z0", "x1", "y1", "z1", "maximum_plane_error"])
values[0, 1] = 0.
write_csv(path, ["x0", "y0", "z0", "x1", "y1", "z1", "maximum_plane_error"], values)
with self.assertRaisesRegex(ValueError, "outline"):
slices.load_slice(self.directory)
def test_shape_exaggeration_is_rejected(self):
path = self.directory / "metadata.csv"
path.write_text(path.read_text().replace("deformation_scale,1", "deformation_scale,10"))
with self.assertRaisesRegex(ValueError, "scale exactly one"):
slices.load_slice(self.directory)
if __name__ == "__main__":
unittest.main()