Compare physical slice refinement with shared scalar rendering

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)