From 9bac57a1df7bb9c6cca949b62eb40c5639d88ded Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Sat, 19 Sep 2026 13:46:11 -0400 Subject: [PATCH] Add physical-plane FEM model slice diagnostics and plots --- experiments/MODEL_SLICES.md | 117 ++++++++++ experiments/model_slice_diagnostics.hpp | 258 +++++++++++++++++++++ experiments/model_slice_geometry.hpp | 194 ++++++++++++++++ experiments/model_slice_geometry_tests.cpp | 13 ++ experiments/plateau_diagnostic.cpp | 34 ++- experiments/plot_model_slices.py | 170 ++++++++++++++ experiments/test_plot_model_slices.py | 103 ++++++++ 7 files changed, 888 insertions(+), 1 deletion(-) create mode 100644 experiments/MODEL_SLICES.md create mode 100644 experiments/model_slice_diagnostics.hpp create mode 100644 experiments/model_slice_geometry.hpp create mode 100644 experiments/model_slice_geometry_tests.cpp create mode 100644 experiments/plot_model_slices.py create mode 100644 experiments/test_plot_model_slices.py diff --git a/experiments/MODEL_SLICES.md b/experiments/MODEL_SLICES.md new file mode 100644 index 0000000..d705323 --- /dev/null +++ b/experiments/MODEL_SLICES.md @@ -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. diff --git a/experiments/model_slice_diagnostics.hpp b/experiments/model_slice_diagnostics.hpp new file mode 100644 index 0000000..b8a2ff1 --- /dev/null +++ b/experiments/model_slice_diagnostics.hpp @@ -0,0 +1,258 @@ +#pragma once + +#include "model_slice_geometry.hpp" +#include "polytrope_physical_state.hpp" + +#include +#include +#include +#include +#include + +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 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, 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 lattice(static_cast(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 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 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 segmentCounts{}; + std::array coplanar{}; + std::array 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 lattice(static_cast(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 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{0, 1, 3}, std::array{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(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 planes = options.plane == "both" ? std::vector{"xz", "xy"} : std::vector{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 diff --git a/experiments/model_slice_geometry.hpp b/experiments/model_slice_geometry.hpp new file mode 100644 index 0000000..408e47e --- /dev/null +++ b/experiments/model_slice_geometry.hpp @@ -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 +#include +#include +#include +#include +#include +#include +#include + +namespace experiment::model_slices { +using Vector = std::array; +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; +using Segment = std::array; + +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 +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& 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::vector EdgeIntersections(const std::array& points, const Plane& plane, + double tolerance, Intersection&& intersect) { + std::vector 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 +std::vector CutTetrahedron(const std::array& 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 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 +std::vector 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, 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 GeometrySelfChecks() { + std::vector 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 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 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 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 diff --git a/experiments/model_slice_geometry_tests.cpp b/experiments/model_slice_geometry_tests.cpp new file mode 100644 index 0000000..44004fb --- /dev/null +++ b/experiments/model_slice_geometry_tests.cpp @@ -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 + +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; + } +} diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index 6cf9b07..70262a5 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -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}; diff --git a/experiments/plot_model_slices.py b/experiments/plot_model_slices.py new file mode 100644 index 0000000..3415ea0 --- /dev/null +++ b/experiments/plot_model_slices.py @@ -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() diff --git a/experiments/test_plot_model_slices.py b/experiments/test_plot_model_slices.py new file mode 100644 index 0000000..2c5f396 --- /dev/null +++ b/experiments/test_plot_model_slices.py @@ -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()