From 6614a5bdf65aa2faee6c50991e6e910af256280b Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Sat, 19 Sep 2026 14:15:35 -0400 Subject: [PATCH] Compare physical slice refinement with shared scalar rendering --- experiments/MODEL_SLICES.md | 48 +++++++ experiments/compare_model_slices.py | 162 +++++++++++++++++++++++ experiments/plot_model_slices.py | 86 +++++++++--- experiments/test_compare_model_slices.py | 124 +++++++++++++++++ experiments/test_plot_model_slices.py | 16 +++ 5 files changed, 419 insertions(+), 17 deletions(-) create mode 100644 experiments/compare_model_slices.py create mode 100644 experiments/test_compare_model_slices.py diff --git a/experiments/MODEL_SLICES.md b/experiments/MODEL_SLICES.md index 44fcafa..5384b15 100644 --- a/experiments/MODEL_SLICES.md +++ b/experiments/MODEL_SLICES.md @@ -56,6 +56,16 @@ 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. +Field values are linearly interpolated within each supplied triangle at the +centers of a 512-by-512 scalar grid, then mapped to colors. This avoids Gouraud +shading's interpolation of RGBA colors, which can change an unchanged affine +field's appearance under triangle refinement. `--pixels N` changes this display +resolution (32–2048). Uncovered pixels remain masked; a pixel covered by more +than one approximate triangle takes the first supplied triangle, with no +cross-element averaging. No global triangulation or global point-location +topology is inferred. Pixels can miss subpixel features and straddle a boundary; +the separately drawn physical outline remains the geometric reference. + 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, @@ -115,6 +125,36 @@ neighboring curved hexahedra and do not replace this export refinement check. Globally consistent face triangulation would remove the diagonal mismatch; it is not implemented in this bounded diagnostic. +The comparison helper enforces a common plane/orientation/offset, increasing +subdivision counts, shared physical axes and shared per-field color scales: + +```sh +python experiments/compare_model_slices.py \ + /absolute/path/to/coarse/final_slices/xz \ + /absolute/path/to/fine/final_slices/xz \ + --source-state /absolute/path/to/original/accepted_state.txt \ + --reference-outline --title 'Same rotating trial; unconverged' \ + --output /absolute/path/outside/source/slice_refinement +``` + +It writes a side-by-side PNG/PDF and JSON containing triangle-area sums, signed +piecewise-linear 2D density/enthalpy integrals, surface bounds, segment-length +sums, exported geometric diagnostics and their absolute/relative changes. Area +and field integrals include any overlapping triangles; segment-length sums +include any duplicate segments. These are measures of the displayed slice, +not exact curved-FEM integrals, 3D mass or virial diagnostics. Near-zero signed +integrals make relative changes sensitive, so also inspect absolute changes. +The two panels use the same scalar raster grid and bounds; its pixel resolution +and overlap ownership are recorded. Quantitative integrals are computed from +the supplied triangles, independently of the display raster. + +`--source-state` is the caller's declaration of the common saved state and its +SHA256 is recorded. The caller must verify that both exports used that state, +mesh, field orders and physical specifications. The current export CSVs do not +bind themselves to a snapshot, and this helper does not establish provenance +from matching plane metadata. Output artifacts must be outside the source +checkout and existing files are not overwritten. + ## Checks Every export runs pure geometric controls: unit-cube section area, integration @@ -126,6 +166,7 @@ a coincident internal face. The same controls can run without FEM: 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 +python -m unittest discover -s experiments -p test_compare_model_slices.py -v ``` The Python controls cover signed independent fields, supplied connectivity for @@ -134,3 +175,10 @@ 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. + +The comparison controls use an affine signed density and independent enthalpy +on a unit square at two subdivisions. They check exact area/integrals/perimeter, +explicit counting of duplicated triangles, shared axes/color normalizations, +and rejection of a mismatched offset or reversed refinement. They do not run a +FEM model. An additional affine control verifies that coarse and fine scalar +rasters agree, and the disconnected-material control verifies masked vacuum. diff --git a/experiments/compare_model_slices.py b/experiments/compare_model_slices.py new file mode 100644 index 0000000..5eeb40e --- /dev/null +++ b/experiments/compare_model_slices.py @@ -0,0 +1,162 @@ +#!/usr/bin/env python3 +"""Compare two visualization resolutions of the same physical FEM state.""" +from __future__ import annotations + +import argparse +import hashlib +import json +from pathlib import Path + +import numpy as np + +from plot_model_slices import Slice, color_scale, load_slice, plt, render + + +def slice_metrics(data: Slice) -> dict: + """Integrate the displayed triangles, not the underlying curved FE slice.""" + u, v = int(data.metadata["horizontal_axis"]), int(data.metadata["vertical_axis"]) + points = data.vertices[data.triangles] + x, y = points[:, :, 2 + u], points[:, :, 2 + v] + area = .5 * np.abs((x[:, 1] - x[:, 0]) * (y[:, 2] - y[:, 0]) - + (x[:, 2] - x[:, 0]) * (y[:, 1] - y[:, 0])) + metrics = { + "triangle_count": len(data.triangles), + "triangle_area_sum": float(np.sum(area)), + "piecewise_linear_density_integral": float(np.sum(area * points[:, :, 5].mean(axis=1))), + "piecewise_linear_enthalpy_integral": float(np.sum(area * points[:, :, 6].mean(axis=1))), + } + metrics["exported_geometric_diagnostics"] = { + key: float(data.metadata[key]) if key in data.metadata else None + for key in ("maximum_plane_error", "maximum_sampled_chord_error", "surface_maximum_plane_error", + "reference_surface_maximum_plane_error", "surface_maximum_sampled_chord_error", + "reference_surface_maximum_sampled_chord_error") + } + for name, segments in (("surface", data.surface), ("reference_surface", data.reference_surface)): + if len(segments): + a, b = segments[:, [u, v]], segments[:, [3 + u, 3 + v]] + endpoints = np.concatenate((a, b)) + values = {"segment_count": len(segments), + "segment_length_sum": float(np.linalg.norm(b - a, axis=1).sum()), + "u_min": float(endpoints[:, 0].min()), "u_max": float(endpoints[:, 0].max()), + "v_min": float(endpoints[:, 1].min()), "v_max": float(endpoints[:, 1].max())} + else: + values = {"segment_count": 0, "segment_length_sum": 0., "u_min": None, "u_max": None, "v_min": None, "v_max": None} + metrics[name] = values + return metrics + + +def differences(coarse: dict, fine: dict) -> dict: + result = {} + for key in coarse: + a, b = coarse[key], fine[key] + if isinstance(a, dict): + result[key] = differences(a, b) + elif a is None or b is None: + result[key] = {"coarse": a, "fine": b, "signed_change": None, "absolute_change": None, + "relative_change": None, "relative_absolute_change": None} + else: + delta = b - a + relative = delta / abs(a) if a != 0 else (0. if delta == 0 else None) + result[key] = {"coarse": a, "fine": b, "signed_change": delta, "absolute_change": abs(delta), + "relative_change": relative, + "relative_absolute_change": abs(relative) if relative is not None else None} + return result + + +def validate_pair(coarse: Slice, fine: Slice) -> None: + if coarse.directory.resolve() == fine.directory.resolve(): + raise ValueError("Coarse and fine exports must be distinct directories") + for key in ("plane", "horizontal_axis", "vertical_axis", "normal_axis"): + if coarse.metadata[key] != fine.metadata[key]: + raise ValueError("The paired exports have different physical planes") + if float(coarse.metadata["offset"]) != float(fine.metadata["offset"]): + raise ValueError("The paired exports have different physical plane offsets") + if int(coarse.metadata["subdivisions"]) >= int(fine.metadata["subdivisions"]): + raise ValueError("Fine subdivisions must exceed coarse subdivisions") + + +def comparison(coarse: Slice, fine: Slice, source_state: Path) -> dict: + validate_pair(coarse, fine) + if not source_state.is_file(): + raise ValueError("--source-state must identify the caller-verified common saved state file") + digest = hashlib.sha256() + with source_state.open("rb") as stream: + for chunk in iter(lambda: stream.read(1024 * 1024), b""): + digest.update(chunk) + return { + "coarse_directory": str(coarse.directory.resolve()), "fine_directory": str(fine.directory.resolve()), + "plane": coarse.metadata["plane"], "offset": float(coarse.metadata["offset"]), + "coarse_subdivisions": int(coarse.metadata["subdivisions"]), "fine_subdivisions": int(fine.metadata["subdivisions"]), + "state_provenance": { + "caller_asserted_common_source": str(source_state.resolve()), "source_sha256": digest.hexdigest(), + "export_binding_verified": False, + "limitation": "Caller must verify that both exports used this same state, mesh, field orders and physical specifications. The CSVs do not bind themselves cryptographically to the source state.", + }, + "changes": differences(slice_metrics(coarse), slice_metrics(fine)), + "scope": [ + "Triangle area and field integrals sum supplied triangles, including any overlaps; they are not union-area measures.", + "Field integrals are signed integrals of piecewise-linear vertex data in 2D; they are not 3D mass, pressure, virial or exact curved-FEM integrals.", + "Surface segment-length sums approximate perimeter and include any duplicated/overlapping segments.", + "Relative change is (fine-coarse)/abs(coarse); it is null for a zero coarse value and a nonzero fine value.", + "Near-zero coarse signed integrals can have large relative changes; inspect absolute changes as well.", + "Shared physical axes and separate shared density/enthalpy color normalizations preserve actual deformation and signed field values.", + "Agreement under one refinement does not certify global geometric conformity, positivity or equilibrium.", + ], + } + + +def render_comparison(coarse: Slice, fine: Slice, *, reference_outline=False, title=None, pixels=512): + validate_pair(coarse, fine) + bounds = (min(float(d.metadata["u_min"]) for d in (coarse, fine)), + max(float(d.metadata["u_max"]) for d in (coarse, fine)), + min(float(d.metadata["v_min"]) for d in (coarse, fine)), + max(float(d.metadata["v_max"]) for d in (coarse, fine))) + normalizations = [color_scale(np.concatenate((coarse.vertices[:, column], fine.vertices[:, column]))) + for column in (5, 6)] if len(coarse.vertices) + len(fine.vertices) else None + figure, axes = plt.subplots(2, 2, figsize=(13.4, 10.6), layout="constrained") + for column, (label, data) in enumerate((("Coarse", coarse), ("Fine", fine))): + render(data, reference_outline=reference_outline, + title=f"{label}: {data.metadata['subdivisions']} subdivisions", + normalizations=normalizations, physical_bounds=bounds, + figure=figure, axes=axes[:, column], heading=False, pixels=pixels) + figure.suptitle(f"{title or 'Same-state slice refinement'} — physical {coarse.metadata['plane']} plane, offset {float(coarse.metadata['offset']):.6g}\n" + f"Common axes/color scales; {pixels}×{pixels} scalar pixels; actual deformation; equilibrium not assessed", fontsize=11) + return figure + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("coarse", type=Path, help="Coarse plane directory, e.g. first/final_slices/xz") + parser.add_argument("fine", type=Path, help="Fine plane directory for the same saved state") + parser.add_argument("--source-state", required=True, type=Path, + help="Common saved state file; caller asserts both exports used this state and matching discretization") + parser.add_argument("--output", required=True, type=Path, help="New JSON/PNG/PDF prefix outside the source checkout") + parser.add_argument("--reference-outline", action="store_true") + parser.add_argument("--pixels", type=int, default=512, help="Common scalar sample grid per axis, 32–2048") + parser.add_argument("--title") + args = parser.parse_args() + source_checkout = Path(__file__).resolve().parents[1] + if args.output.resolve().is_relative_to(source_checkout): + raise ValueError("Comparison artifacts must be written outside the source checkout") + paths = [args.output.with_suffix(suffix) for suffix in (".json", ".png", ".pdf")] + if any(path.exists() for path in paths): + raise FileExistsError("Refusing to overwrite an existing comparison artifact") + coarse, fine = load_slice(args.coarse), load_slice(args.fine) + report = comparison(coarse, fine, args.source_state) + figure = render_comparison(coarse, fine, reference_outline=args.reference_outline, title=args.title, pixels=args.pixels) + report["rasterization"] = {"pixels_per_axis": args.pixels, "sample_location": "pixel centers", + "scalar_interpolation": "barycentric within each supplied triangle before color mapping", + "overlap_owner": "first supplied triangle; no trace averaging", "vacuum": "masked", + "limitation": "Subpixel features may be missed; this raster is not an exact curved-FEM evaluation."} + report["figures"] = [str(path.resolve()) for path in paths[1:]] + args.output.parent.mkdir(parents=True, exist_ok=True) + paths[0].write_text(json.dumps(report, indent=2, allow_nan=False) + "\n") + for path in paths[1:]: + figure.savefig(path, dpi=200) + plt.close(figure) + for path in paths: + print(path) + + +if __name__ == "__main__": + main() diff --git a/experiments/plot_model_slices.py b/experiments/plot_model_slices.py index 3415ea0..31b3fe2 100644 --- a/experiments/plot_model_slices.py +++ b/experiments/plot_model_slices.py @@ -12,7 +12,6 @@ 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 @@ -90,21 +89,69 @@ def color_scale(values: np.ndarray): return Normalize(lower, upper), "viridis", "linear color scale" -def render(data: Slice, *, reference_outline: bool = False, title: str | None = None): +def rasterize_fields(data: Slice, bounds, pixels=512): + """Evaluate each triangle's linear scalar interpolant at pixel centers. + + No global point-location topology is assumed: curved-cell chords can have + gaps or overlaps. In an overlap the first supplied triangle owns the pixel. + This deterministic display rule does not turn overlaps into a union-area + diagnostic and never averages independent element traces. + """ + if not 32 <= pixels <= 2048: + raise ValueError("Raster resolution must be between 32 and 2048 pixels per axis") + u0, u1, v0, v1 = bounds + if not np.isfinite(bounds).all() or u1 <= u0 or v1 <= v0: + raise ValueError("Invalid scalar raster bounds") + du, dv = (u1 - u0) / pixels, (v1 - v0) / pixels + us = u0 + (np.arange(pixels) + .5) * du + vs = v0 + (np.arange(pixels) + .5) * dv + values = np.zeros((2, pixels, pixels)) + covered = np.zeros((pixels, pixels), dtype=bool) + u, v = int(data.metadata["horizontal_axis"]), int(data.metadata["vertical_axis"]) + for indices in data.triangles: + triangle = data.vertices[indices] + x, y = triangle[:, 2 + u], triangle[:, 2 + v] + left = max(0, int(np.ceil((x.min() - u0) / du - .5 - 1e-12))) + right = min(pixels - 1, int(np.floor((x.max() - u0) / du - .5 + 1e-12))) + bottom = max(0, int(np.ceil((y.min() - v0) / dv - .5 - 1e-12))) + top = min(pixels - 1, int(np.floor((y.max() - v0) / dv - .5 + 1e-12))) + if left > right or bottom > top: + continue + dx1, dx2, dy1, dy2 = x[1] - x[0], x[2] - x[0], y[1] - y[0], y[2] - y[0] + determinant = dx1 * dy2 - dy1 * dx2 + if determinant == 0: + continue + dx, dy = us[None, left:right + 1] - x[0], vs[bottom:top + 1, None] - y[0] + a = (dx * dy2 - dy * dx2) / determinant + b = (dy * dx1 - dx * dy1) / determinant + ownership = covered[bottom:top + 1, left:right + 1] + inside = (a >= -2e-13) & (b >= -2e-13) & (a + b <= 1 + 2e-13) & ~ownership + for field, column in enumerate((5, 6)): + scalar = triangle[0, column] + a * (triangle[1, column] - triangle[0, column]) + b * (triangle[2, column] - triangle[0, column]) + values[field, bottom:top + 1, left:right + 1][inside] = scalar[inside] + ownership[inside] = True + return np.ma.array(values, mask=np.broadcast_to(~covered, values.shape)) + + +def render(data: Slice, *, reference_outline: bool = False, title: str | None = None, + normalizations=None, physical_bounds=None, figure=None, axes=None, heading=True, pixels=512): 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: + if figure is None: + figure, axes = plt.subplots(1, 2, figsize=(12.4, 5.6), layout="constrained") + elif axes is None or len(axes) != 2: + raise ValueError("An existing figure requires two field axes") + bounds = physical_bounds if physical_bounds is not None else (float(m["u_min"]), float(m["u_max"]), float(m["v_min"]), float(m["v_max"])) + raster = rasterize_fields(data, bounds, pixels) if len(data.triangles) else None + for field, (ax, column, name) in enumerate(zip(axes, (5, 6), ("Independent density", "Specific enthalpy"))): + if raster 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) + norm, cmap, note = color_scale(values) if normalizations is None else normalizations[field] + # Colormapping follows scalar interpolation. Gouraud shading would + # interpolate RGBA instead, changing an affine field's appearance + # under triangle refinement even with identical scalar norms. + artist = ax.imshow(raster[field], extent=bounds, origin="lower", interpolation="none", 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)] @@ -120,8 +167,8 @@ def render(data: Slice, *, reference_outline: bool = False, title: str | None = 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"])]) + lower = np.array([bounds[0], bounds[2]]) + upper = np.array([bounds[1], bounds[3]]) 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) @@ -142,8 +189,12 @@ def render(data: Slice, *, reference_outline: bool = False, title: str | None = 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) + if heading: + 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']}; {pixels}×{pixels} scalar pixels", fontsize=11) + else: + for ax in axes: + ax.set_title(f"{description}\n{ax.get_title()}", fontsize=10) return figure @@ -152,13 +203,14 @@ def main() -> None: 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("--pixels", type=int, default=512, help="Scalar sample grid per axis, 32–2048 (default 512)") 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) + figure = render(load_slice(args.directory), reference_outline=args.reference_outline, title=args.title, pixels=args.pixels) output.parent.mkdir(parents=True, exist_ok=True) for path in paths: figure.savefig(path, dpi=200) diff --git a/experiments/test_compare_model_slices.py b/experiments/test_compare_model_slices.py new file mode 100644 index 0000000..c86ecf2 --- /dev/null +++ b/experiments/test_compare_model_slices.py @@ -0,0 +1,124 @@ +#!/usr/bin/env python3 +"""Affine same-state controls for physical-slice refinement summaries.""" +import csv +from pathlib import Path +import tempfile +import unittest + +import numpy as np + +import compare_model_slices as compare +from plot_model_slices import load_slice, plt, rasterize_fields + + +def write_csv(path, header, rows): + with path.open("w", newline="") as stream: + writer = csv.writer(stream) + writer.writerow(header) + writer.writerows(rows) + + +def square(directory, subdivisions, offset=.2): + directory.mkdir() + vertices, triangles = [], [] + n = subdivisions + for j in range(n): + for i in range(n): + corners = [(i / n, j / n), ((i + 1) / n, j / n), (i / n, (j + 1) / n), ((i + 1) / n, (j + 1) / n)] + for indices in ((0, 1, 3), (0, 3, 2)): + start = len(vertices) + triangles.append([len(triangles), 0, start, start + 1, start + 2]) + for corner in indices: + x, z = corners[corner] + vertices.append([len(vertices), 0, x, offset, z, x - .5, 1 + x + 2 * z, 0]) + 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"], triangles) + segments = [] + for i in range(n): + a, b = i / n, (i + 1) / n + for x0, z0, x1, z1 in ((a, 0, b, 0), (1, a, 1, b), (b, 1, a, 1), (0, b, 0, a)): + segments.append([x0, offset, z0, x1, offset, z1, 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=offset, + root_tolerance=1e-11, deformation_scale=1, vertex_count=len(vertices), triangle_count=len(triangles), + surface_segments=len(segments), reference_surface_segments=len(segments), surface_coplanar_triangles=0, + reference_surface_coplanar_triangles=0, subdivisions=n, u_min=0, u_max=1, v_min=0, v_max=1) + write_csv(directory / "metadata.csv", ["metric", "value"], metadata.items()) + + +class ComparisonTests(unittest.TestCase): + def setUp(self): + self.temporary = tempfile.TemporaryDirectory() + self.root = Path(self.temporary.name) + square(self.root / "coarse", 2) + square(self.root / "fine", 4) + self.coarse = load_slice(self.root / "coarse") + self.fine = load_slice(self.root / "fine") + self.source = self.root / "synthetic_state.txt" + self.source.write_text("rho=x-0.5; h=1+x+2z; unit square; y=0.2\n") + + def tearDown(self): + plt.close("all") + self.temporary.cleanup() + + def test_affine_integrals_and_perimeter_do_not_change_under_refinement(self): + report = compare.comparison(self.coarse, self.fine, self.source) + for key, expected in (("triangle_area_sum", 1.), ("piecewise_linear_density_integral", 0.), + ("piecewise_linear_enthalpy_integral", 2.5)): + self.assertAlmostEqual(report["changes"][key]["coarse"], expected, places=14) + self.assertAlmostEqual(report["changes"][key]["fine"], expected, places=14) + self.assertLess(report["changes"][key]["absolute_change"], 1e-14) + self.assertEqual(report["changes"]["surface"]["segment_length_sum"]["fine"], 4.) + self.assertEqual(report["changes"]["surface"]["u_min"]["fine"], 0.) + self.assertEqual(report["changes"]["surface"]["v_max"]["fine"], 1.) + self.assertFalse(report["state_provenance"]["export_binding_verified"]) + self.assertEqual(len(report["state_provenance"]["source_sha256"]), 64) + + def test_overlap_is_counted_not_hidden_as_union_area(self): + self.fine.triangles = np.concatenate((self.fine.triangles, self.fine.triangles)) + metrics = compare.slice_metrics(self.fine) + self.assertEqual(metrics["triangle_area_sum"], 2.) + self.assertEqual(metrics["piecewise_linear_enthalpy_integral"], 5.) + + def test_shared_axes_and_signed_color_ranges(self): + figure = compare.render_comparison(self.coarse, self.fine, reference_outline=True) + axes = figure.axes[:4] + for ax in axes: + self.assertEqual(ax.get_xlim(), axes[0].get_xlim()) + self.assertEqual(ax.get_ylim(), axes[0].get_ylim()) + self.assertEqual(ax.get_aspect(), 1.) + self.assertIs(axes[0].images[0].norm, axes[1].images[0].norm) + self.assertIs(axes[2].images[0].norm, axes[3].images[0].norm) + self.assertEqual(axes[0].images[0].norm.vmin, -.5) + self.assertEqual(axes[0].images[0].norm.vmax, .5) + for extension in ("png", "pdf"): + path = self.root / ("comparison." + extension) + figure.savefig(path) + self.assertGreater(path.stat().st_size, 1000) + + def test_affine_scalar_raster_is_refinement_invariant(self): + a = rasterize_fields(self.coarse, (0, 1, 0, 1), 65) + b = rasterize_fields(self.fine, (0, 1, 0, 1), 65) + self.assertFalse(a.mask.any()) + self.assertFalse(b.mask.any()) + np.testing.assert_allclose(a.data, b.data, rtol=0, atol=1e-14) + + def test_different_plane_offset_is_rejected(self): + square(self.root / "wrong_offset", 4, offset=.3) + with self.assertRaisesRegex(ValueError, "offsets"): + compare.comparison(self.coarse, load_slice(self.root / "wrong_offset"), self.source) + + def test_reversed_refinement_is_rejected(self): + with self.assertRaisesRegex(ValueError, "exceed"): + compare.comparison(self.fine, self.coarse, self.source) + + def test_relative_zero_denominator_is_explicit(self): + result = compare.differences({"zero": 0., "negative": -2.}, {"zero": 1., "negative": -3.}) + self.assertIsNone(result["zero"]["relative_change"]) + self.assertEqual(result["negative"]["relative_change"], -.5) + self.assertEqual(result["negative"]["relative_absolute_change"], .5) + + +if __name__ == "__main__": + unittest.main() diff --git a/experiments/test_plot_model_slices.py b/experiments/test_plot_model_slices.py index 2c5f396..3caba75 100644 --- a/experiments/test_plot_model_slices.py +++ b/experiments/test_plot_model_slices.py @@ -64,6 +64,8 @@ class ModelSlicePlotTests(unittest.TestCase): for triangle in data.triangles: horizontal = data.vertices[triangle, 2] self.assertTrue(np.all(horizontal < 0) or np.all(horizontal > 0)) + raster = slices.rasterize_fields(data, (-2, 2, 0, 1), 64) + self.assertTrue(raster.mask[:, :, 24:40].all()) 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.) @@ -77,6 +79,20 @@ class ModelSlicePlotTests(unittest.TestCase): with self.assertRaisesRegex(ValueError, "different material elements"): slices.load_slice(self.directory) + def test_raster_overlap_uses_first_triangle_without_trace_averaging(self): + data = slices.load_slice(self.directory) + original = slices.rasterize_fields(data, (-2, 2, 0, 1), 64) + extra = data.vertices[:3].copy() + extra[:, 0] += len(data.vertices) + extra[:, 1] = 2 + extra[:, 5] = 100 + extra[:, 6] = -100 + data.vertices = np.concatenate((data.vertices, extra)) + data.triangles = np.concatenate((data.triangles, [[6, 7, 8]])) + overlapped = slices.rasterize_fields(data, (-2, 2, 0, 1), 64) + np.testing.assert_array_equal(original.mask, overlapped.mask) + np.testing.assert_array_equal(original.data, overlapped.data) + 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,")