plotting.meshing — mesh and field inspection

tempura/plotting/meshing.py is the inspection toolkit. It is read-only: it consumes a finalized pescado problem plus Tempura's region shapes and renders 2D plane cuts of the mesh and of solved scalar fields, for debugging, documentation, and regression figures. It is a secondary helper module, not part of the core modeling contract.

Plane cuts

The central idea is the PlaneSpec: a description of one 2D slice through the 3D problem — which axis is fixed, at what value, and which two axes are displayed. make_standard_plane_specs(...) builds the usual XY/XZ/YZ trio for a device, nearest_axis_value(...) snaps a requested cut to the nearest realized mesh plane, and transpose_plane_spec(...) / make_xy_emphasis_axes(...) arrange the panels.

Rendering

plot_problem_regions_with_mesh(...) draws the solver-owned regions with the control-volume edges of the discretization overlaid, so you can see exactly how the mesh resolves each layer. plot_scalar_field_on_planes(...) renders a solved field (for example a superposed gate potential) on the same plane specs, and the contour/footprint helpers add region outlines and gate footprints.

The worked examples (Triple Quantum Dot, Minimal Kitaev Chain) use these helpers end to end.

API

plotting

Plotting helpers for Tempura.

PlaneSpec dataclass

Describe one 2D slice through a finalized 3D problem.

Attributes:

Name	Type	Description
title	str	Panel title.
plot_region	Box	2D plotting region in the displayed axes.
plane_axis	int	Fixed axis in 3D coordinates.
plane_value	float	Position of the fixed slicing plane.
axis_indices	tuple[int, int]	Indices of the displayed axes in the 3D coordinate system.
axis_labels	tuple[str, str]	Labels for the displayed axes.

Source code in tempura/plotting/meshing.py

@dataclass(frozen=True)
class PlaneSpec:
    """Describe one 2D slice through a finalized 3D problem.

    Attributes:
        title: Panel title.
        plot_region: 2D plotting region in the displayed axes.
        plane_axis: Fixed axis in 3D coordinates.
        plane_value: Position of the fixed slicing plane.
        axis_indices: Indices of the displayed axes in the 3D coordinate system.
        axis_labels: Labels for the displayed axes.
    """

    title: str
    plot_region: shapes.Box
    plane_axis: int
    plane_value: float
    axis_indices: tuple[int, int]
    axis_labels: tuple[str, str]

    def __post_init__(self) -> None:
        """Validate that the slice description is internally consistent."""
        bbox = np.asarray(self.plot_region.bbox, dtype=float)
        if bbox.shape != (2, 2):
            raise ValueError(
                "PlaneSpec.plot_region must be a 2D shape with bbox shape (2, 2); "
                f"got {bbox.shape}."
            )
        if len(self.axis_indices) != 2 or len(set(self.axis_indices)) != 2:
            raise ValueError(
                "PlaneSpec.axis_indices must contain exactly two distinct axes; "
                f"got {self.axis_indices}."
            )
        if self.plane_axis in self.axis_indices:
            raise ValueError(
                "PlaneSpec.plane_axis must differ from the plotted axis indices; "
                f"got plane_axis={self.plane_axis}, axis_indices={self.axis_indices}."
            )
        if len(self.axis_labels) != 2:
            raise ValueError(
                "PlaneSpec.axis_labels must contain exactly two strings; "
                f"got {self.axis_labels}."
            )

__post_init__()

Validate that the slice description is internally consistent.

Source code in tempura/plotting/meshing.py

def __post_init__(self) -> None:
    """Validate that the slice description is internally consistent."""
    bbox = np.asarray(self.plot_region.bbox, dtype=float)
    if bbox.shape != (2, 2):
        raise ValueError(
            "PlaneSpec.plot_region must be a 2D shape with bbox shape (2, 2); "
            f"got {bbox.shape}."
        )
    if len(self.axis_indices) != 2 or len(set(self.axis_indices)) != 2:
        raise ValueError(
            "PlaneSpec.axis_indices must contain exactly two distinct axes; "
            f"got {self.axis_indices}."
        )
    if self.plane_axis in self.axis_indices:
        raise ValueError(
            "PlaneSpec.plane_axis must differ from the plotted axis indices; "
            f"got plane_axis={self.plane_axis}, axis_indices={self.axis_indices}."
        )
    if len(self.axis_labels) != 2:
        raise ValueError(
            "PlaneSpec.axis_labels must contain exactly two strings; "
            f"got {self.axis_labels}."
        )

add_gate_footprint_contours_xy(ax, x_values, y_values, footprints, *, coordinate_scale=1.0, color='black', linewidth=0.8, alpha=0.7)

Overlay direct XY gate-footprint contours on one existing axis.

Parameters:

Name	Type	Description	Default
ax	Axes	Axis on which to draw the contours.	required
x_values	ndarray	Plane x coordinates used by the footprint masks.	required
y_values	ndarray	Plane y coordinates used by the footprint masks.	required
footprints	Mapping[str, ndarray]	Mapping of gate name to (ny, nx) footprint mask.	required
coordinate_scale	float	Multiplicative scale applied to the displayed coordinates.	1.0
color	str	Contour color.	'black'
linewidth	float	Contour line width.	0.8
alpha	float	Contour alpha.	0.7

Returns:

Type	Description
None	None. Contours are added directly to ax.

Source code in tempura/plotting/meshing.py

def add_gate_footprint_contours_xy(
    ax: plt.Axes,
    x_values: np.ndarray,
    y_values: np.ndarray,
    footprints: Mapping[str, np.ndarray],
    *,
    coordinate_scale: float = 1.0,
    color: str = "black",
    linewidth: float = 0.8,
    alpha: float = 0.7,
) -> None:
    """Overlay direct XY gate-footprint contours on one existing axis.

    Args:
        ax: Axis on which to draw the contours.
        x_values: Plane x coordinates used by the footprint masks.
        y_values: Plane y coordinates used by the footprint masks.
        footprints: Mapping of gate name to ``(ny, nx)`` footprint mask.
        coordinate_scale: Multiplicative scale applied to the displayed
            coordinates.
        color: Contour color.
        linewidth: Contour line width.
        alpha: Contour alpha.

    Returns:
        ``None``. Contours are added directly to ``ax``.
    """
    if not footprints:
        return
    scale = float(coordinate_scale)
    x_vals = np.asarray(x_values, dtype=float) * scale
    y_vals = np.asarray(y_values, dtype=float) * scale
    for mask in footprints.values():
        arr = np.asarray(mask, dtype=bool)
        if not np.any(arr):
            continue
        ax.contour(
            x_vals,
            y_vals,
            arr.astype(float),
            levels=[0.5],
            colors=[color],
            linewidths=linewidth,
            alpha=alpha,
            zorder=2,
        )

add_plane_cut_guides(ax, *, x=None, y=None, coordinate_scale=1.0, color='black', linewidth=1.0, linestyle='--')

Overlay one or two cut-guide lines on an existing 2D plane plot.

Parameters:

Name	Type	Description	Default
ax	Axes	Axis on which to draw the guide lines.	required
x	float | None	Optional vertical guide x coordinate.	None
y	float | None	Optional horizontal guide y coordinate.	None
coordinate_scale	float	Multiplicative scale applied to the displayed coordinates.	1.0
color	str	Line color.	'black'
linewidth	float	Line width.	1.0
linestyle	str	Matplotlib line style.	'--'

Returns:

Type	Description
None	None. Guide lines are added directly to ax.

Source code in tempura/plotting/meshing.py

def add_plane_cut_guides(
    ax: plt.Axes,
    *,
    x: float | None = None,
    y: float | None = None,
    coordinate_scale: float = 1.0,
    color: str = "black",
    linewidth: float = 1.0,
    linestyle: str = "--",
) -> None:
    """Overlay one or two cut-guide lines on an existing 2D plane plot.

    Args:
        ax: Axis on which to draw the guide lines.
        x: Optional vertical guide x coordinate.
        y: Optional horizontal guide y coordinate.
        coordinate_scale: Multiplicative scale applied to the displayed
            coordinates.
        color: Line color.
        linewidth: Line width.
        linestyle: Matplotlib line style.

    Returns:
        ``None``. Guide lines are added directly to ``ax``.
    """
    scale = float(coordinate_scale)
    if x is not None:
        ax.axvline(
            float(x) * scale, color=color, linewidth=linewidth, linestyle=linestyle
        )
    if y is not None:
        ax.axhline(
            float(y) * scale, color=color, linewidth=linewidth, linestyle=linestyle
        )

add_region_contours_on_planes(problem, axes, plane_specs, region_shapes, *, region_display_names=None, coordinate_scale=1.0, linewidth=1.0, alpha=0.9)

Overlay region-outline contours on existing plane plots.

Parameters:

Name	Type	Description	Default
problem	Any	Finalized Pescado problem with a coordinates array.	required
axes	Sequence[Axes]	Matplotlib axes matching plane_specs.	required
plane_specs	Sequence[PlaneSpec]	Plane definitions describing each plotted cut.	required
region_shapes	RegionShapeMap	Mapping of region name to solver-owned shape.	required
region_display_names	RegionDisplayMap | None	Optional mapping that merges several solver regions under one display label.	None
coordinate_scale	float	Multiplicative scale applied to displayed coordinates.	1.0
linewidth	float	Contour line width.	1.0
alpha	float	Contour alpha.	0.9

Returns:

Type	Description
None	None. Contours are added directly to the provided axes.

Source code in tempura/plotting/meshing.py

def add_region_contours_on_planes(
    problem: Any,
    axes: Sequence[plt.Axes],
    plane_specs: Sequence[PlaneSpec],
    region_shapes: RegionShapeMap,
    *,
    region_display_names: RegionDisplayMap | None = None,
    coordinate_scale: float = 1.0,
    linewidth: float = 1.0,
    alpha: float = 0.9,
) -> None:
    """Overlay region-outline contours on existing plane plots.

    Args:
        problem: Finalized Pescado problem with a ``coordinates`` array.
        axes: Matplotlib axes matching ``plane_specs``.
        plane_specs: Plane definitions describing each plotted cut.
        region_shapes: Mapping of region name to solver-owned shape.
        region_display_names: Optional mapping that merges several solver
            regions under one display label.
        coordinate_scale: Multiplicative scale applied to displayed
            coordinates.
        linewidth: Contour line width.
        alpha: Contour alpha.

    Returns:
        ``None``. Contours are added directly to the provided axes.
    """
    if not plane_specs or not region_shapes:
        return

    coordinates = np.asarray(problem.coordinates, dtype=float)
    volume_bounds = _volume_bounds_lookups(coordinates)
    display_names, _ = _display_order(list(region_shapes), region_display_names)
    grouped_names: dict[str, list[str]] = {}
    for region_name in region_shapes:
        display_name = (
            str(region_display_names[region_name])
            if region_display_names is not None and region_name in region_display_names
            else str(region_name)
        )
        grouped_names.setdefault(display_name, []).append(region_name)

    cmap = plt.cm.get_cmap("tab10")
    colors = cmap(np.linspace(0.0, 1.0, max(len(display_names), 1)))

    for ax, spec in zip(axes, plane_specs, strict=False):
        plane_sites = _plane_sites(
            coordinates,
            spec=spec,
            bounds_lookup=volume_bounds[spec.plane_axis],
        )
        if plane_sites.size == 0:
            continue
        x_edges, y_edges = _region_sample_edges(
            plane_sites,
            spec=spec,
            volume_bounds=volume_bounds,
        )
        x_centers = 0.5 * (x_edges[:-1] + x_edges[1:])
        y_centers = 0.5 * (y_edges[:-1] + y_edges[1:])
        xx, yy = np.meshgrid(x_centers, y_centers, indexing="xy")

        sample_points = np.zeros((xx.size, 3), dtype=float)
        sample_points[:, spec.plane_axis] = float(spec.plane_value)
        sample_points[:, spec.axis_indices[0]] = xx.ravel()
        sample_points[:, spec.axis_indices[1]] = yy.ravel()

        for color_index, display_name in enumerate(display_names):
            union = np.zeros(sample_points.shape[0], dtype=bool)
            for region_name in grouped_names.get(display_name, []):
                union |= np.asarray(region_shapes[region_name](sample_points), dtype=bool)
            if not np.any(union):
                continue
            ax.contour(
                xx * float(coordinate_scale),
                yy * float(coordinate_scale),
                union.reshape(len(y_centers), len(x_centers)).astype(float),
                levels=[0.5],
                colors=[colors[color_index]],
                linewidths=linewidth,
                alpha=alpha,
                zorder=2,
            )

expand_plane_grid(plane_specs_by_key, grid, *, title_overrides=None)

Expand a named subplot grid into an ordered plane-spec sequence.

Parameters:

Name	Type	Description	Default
plane_specs_by_key	Mapping[str, PlaneSpec]	Mapping from plane key to :class:PlaneSpec.	required
grid	Sequence[Sequence[str]]	Row-major grid of plane keys.	required
title_overrides	Mapping[tuple[int, int], str] | None	Optional subplot-title overrides keyed by (row_index, col_index).	None

Returns:

Type	Description
list[PlaneSpec]	Tuple of (ordered_specs, ncols) where ordered_specs follows the
int	subplot grid order.

Source code in tempura/plotting/meshing.py

def expand_plane_grid(
    plane_specs_by_key: Mapping[str, PlaneSpec],
    grid: Sequence[Sequence[str]],
    *,
    title_overrides: Mapping[tuple[int, int], str] | None = None,
) -> tuple[list[PlaneSpec], int]:
    """Expand a named subplot grid into an ordered plane-spec sequence.

    Args:
        plane_specs_by_key: Mapping from plane key to :class:`PlaneSpec`.
        grid: Row-major grid of plane keys.
        title_overrides: Optional subplot-title overrides keyed by
            ``(row_index, col_index)``.

    Returns:
        Tuple of ``(ordered_specs, ncols)`` where ``ordered_specs`` follows the
        subplot grid order.
    """
    ordered_specs: list[PlaneSpec] = []
    ncols = 0
    for row_idx, row in enumerate(grid):
        ncols = max(ncols, len(row))
        for col_idx, key in enumerate(row):
            try:
                spec = plane_specs_by_key[key]
            except KeyError as exc:
                raise ValueError(f"Unknown plane key {key!r} in grid.") from exc
            title = (
                None
                if title_overrides is None
                else title_overrides.get((row_idx, col_idx))
            )
            ordered_specs.append(spec if title is None else replace(spec, title=title))
    if not ordered_specs:
        raise ValueError("grid must contain at least one plane key.")
    return ordered_specs, ncols

make_standard_plane_specs(device, *, xy_title, xy_plane_z, xz_title, xz_plane_y, yz_title=None, yz_plane_x=None, z_padding=1.0)

Build the standard XY/XZ/YZ plotting cuts for one device.

Parameters:

Name	Type	Description	Default
device	Any	Finalized device whose bounds define the default plotting regions.	required
xy_title	str	Title for the XY panel.	required
xy_plane_z	float	Fixed z coordinate for the XY cut.	required
xz_title	str	Title for the XZ panel.	required
xz_plane_y	float	Fixed y coordinate for the XZ cut.	required
yz_title	str | None	Optional title for the YZ panel.	None
yz_plane_x	float | None	Optional fixed x coordinate for the YZ cut.	None
z_padding	float	Extra z extent added around the realized layer stack.	1.0

Returns:

Type	Description
dict[str, PlaneSpec]	Mapping containing "XY", "XZ", and optionally "YZ"
dict[str, PlaneSpec]	class:PlaneSpec objects.

Source code in tempura/plotting/meshing.py

def make_standard_plane_specs(
    device: Any,
    *,
    xy_title: str,
    xy_plane_z: float,
    xz_title: str,
    xz_plane_y: float,
    yz_title: str | None = None,
    yz_plane_x: float | None = None,
    z_padding: float = 1.0,
) -> dict[str, PlaneSpec]:
    """Build the standard XY/XZ/YZ plotting cuts for one device.

    Args:
        device: Finalized device whose bounds define the default plotting
            regions.
        xy_title: Title for the XY panel.
        xy_plane_z: Fixed z coordinate for the XY cut.
        xz_title: Title for the XZ panel.
        xz_plane_y: Fixed y coordinate for the XZ cut.
        yz_title: Optional title for the YZ panel.
        yz_plane_x: Optional fixed x coordinate for the YZ cut.
        z_padding: Extra z extent added around the realized layer stack.

    Returns:
        Mapping containing ``"XY"``, ``"XZ"``, and optionally ``"YZ"``
        :class:`PlaneSpec` objects.
    """
    z_values = []
    for layer in device.layers.values():
        bbox = np.asarray(layer.shape.bbox, dtype=float)
        z_values.extend([float(bbox[0, 2]), float(bbox[1, 2])])
    if not z_values:
        raise ValueError("The device does not contain any finalized layers.")

    z_min = min(z_values) - float(z_padding) / 2.0
    z_size = (max(z_values) - min(z_values)) + float(z_padding)

    specs = {
        "XY": PlaneSpec(
            title=xy_title,
            plot_region=shapes.Box(
                lower_left=np.array(
                    [-device.length / 2.0, -device.width / 2.0], dtype=float
                ),
                size=np.array([device.length, device.width], dtype=float),
            ),
            plane_axis=2,
            plane_value=float(xy_plane_z),
            axis_indices=(0, 1),
            axis_labels=("x", "y"),
        ),
        "XZ": PlaneSpec(
            title=xz_title,
            plot_region=shapes.Box(
                lower_left=np.array([-device.length / 2.0, z_min], dtype=float),
                size=np.array([device.length, z_size], dtype=float),
            ),
            plane_axis=1,
            plane_value=float(xz_plane_y),
            axis_indices=(0, 2),
            axis_labels=("x", "z"),
        ),
    }
    if yz_title is not None and yz_plane_x is not None:
        specs["YZ"] = PlaneSpec(
            title=yz_title,
            plot_region=shapes.Box(
                lower_left=np.array([-device.width / 2.0, z_min], dtype=float),
                size=np.array([device.width, z_size], dtype=float),
            ),
            plane_axis=0,
            plane_value=float(yz_plane_x),
            axis_indices=(1, 2),
            axis_labels=("y", "z"),
        )
    return specs

make_xy_emphasis_axes(*, figsize=(12.0, 9.8), layout='stack', width_ratios=(1.0, 1.0), height_ratios=(1.0, 1.0, 1.0))

Return three axes for XY, XZ, and YZ cuts in a compact layout.

Parameters:

Name	Type	Description	Default
figsize	tuple[float, float]	Matplotlib figure size.	(12.0, 9.8)
layout	str	Either "stack" or "left_span".	'stack'
width_ratios	tuple[float, float]	Width ratios used by "left_span".	(1.0, 1.0)
height_ratios	tuple[float, float, float]	Height ratios used by "stack".	(1.0, 1.0, 1.0)

Returns:

Type	Description
Figure	Tuple of (figure, axes) where axes is ordered as XY,
ndarray	XZ, YZ.

Source code in tempura/plotting/meshing.py

def make_xy_emphasis_axes(
    *,
    figsize: tuple[float, float] = (12.0, 9.8),
    layout: str = "stack",
    width_ratios: tuple[float, float] = (1.0, 1.0),
    height_ratios: tuple[float, float, float] = (1.0, 1.0, 1.0),
) -> tuple[plt.Figure, np.ndarray]:
    """Return three axes for XY, XZ, and YZ cuts in a compact layout.

    Args:
        figsize: Matplotlib figure size.
        layout: Either ``"stack"`` or ``"left_span"``.
        width_ratios: Width ratios used by ``"left_span"``.
        height_ratios: Height ratios used by ``"stack"``.

    Returns:
        Tuple of ``(figure, axes)`` where ``axes`` is ordered as ``XY``,
        ``XZ``, ``YZ``.
    """
    fig = plt.figure(figsize=figsize, constrained_layout=True)
    if layout == "stack":
        gs = fig.add_gridspec(3, 1, height_ratios=height_ratios)
        xy_ax = fig.add_subplot(gs[0, 0])
        xz_ax = fig.add_subplot(gs[1, 0])
        yz_ax = fig.add_subplot(gs[2, 0])
    elif layout == "left_span":
        gs = fig.add_gridspec(2, 2, width_ratios=width_ratios)
        xy_ax = fig.add_subplot(gs[:, 0])
        xz_ax = fig.add_subplot(gs[0, 1])
        yz_ax = fig.add_subplot(gs[1, 1])
    else:
        raise ValueError(f"Unknown layout {layout!r}; expected 'stack' or 'left_span'.")
    return fig, np.asarray([xy_ax, xz_ax, yz_ax], dtype=object)

nearest_axis_value(coordinates, axis, target)

Snap a requested cut to the nearest finalized mesh plane.

Parameters:

Name	Type	Description	Default
coordinates	ndarray	Full solver coordinate array with one row per mesh site.	required
axis	int	Cartesian axis to inspect.	required
target	float	Requested coordinate along axis.	required

Returns:

Type	Description
float	Realized mesh coordinate along axis closest to target.

Source code in tempura/plotting/meshing.py

def nearest_axis_value(coordinates: np.ndarray, axis: int, target: float) -> float:
    """Snap a requested cut to the nearest finalized mesh plane.

    Args:
        coordinates: Full solver coordinate array with one row per mesh site.
        axis: Cartesian axis to inspect.
        target: Requested coordinate along ``axis``.

    Returns:
        Realized mesh coordinate along ``axis`` closest to ``target``.
    """
    values = np.unique(np.round(np.asarray(coordinates, dtype=float)[:, axis], 12))
    return float(values[np.argmin(np.abs(values - target))])

plot_problem_regions_with_mesh(problem, region_shapes, plane_specs, *, axes=None, region_names=None, region_display_names=None, ncols=None, figsize=(13.0, 11.0), suptitle='Pescado region cuts with overlaid mesh nodes', fill_mode='regions', background_alpha=0.88, tile_edgecolors='white', tile_linewidth=0.6, show_volume_edges=False, show_mesh_points=True, mesh_point_size=12.0, mesh_facecolors='white', mesh_edgecolors='black', mesh_linewidth=0.35, coordinate_scale=1.0, axis_unit=None)

Plot finalized region slices with overlaid mesh nodes.

Parameters:

Name	Type	Description	Default
problem	Any	Finalized Pescado problem with a coordinates array.	required
region_shapes	RegionShapeMap	Mapping returned by build_problem(...). These are the solver-owned region shapes after Tempura's lower-face trimming.	required
plane_specs	Sequence[PlaneSpec]	Sequence of PlaneSpec describing the cuts to draw.	required
region_names	Sequence[str] | None	Optional region order for coloring. Defaults to every region except "Boundary" in insertion order.	None
region_display_names	RegionDisplayMap | None	Optional mapping from solver-owned region names to display labels. Regions that share a display label are drawn with the same color and collapsed into one legend entry.	None
ncols	int | None	Number of subplot columns. Defaults to 2 when more than one plane is requested and 1 otherwise.	None
figsize	tuple[float, float]	Matplotlib figure size.	(13.0, 11.0)
suptitle	str	Figure title.	'Pescado region cuts with overlaid mesh nodes'
fill_mode	str	"regions" to render the material cross-section from the finalized region shapes, or "control_volumes" to render the selected mesh control volumes directly.	'regions'
background_alpha	float	Alpha used for the ownership tiles under the markers.	0.88
tile_edgecolors	str	Edge color for the projected mesh control volumes.	'white'
tile_linewidth	float	Edge line width for the projected mesh control volumes.	0.6
show_volume_edges	bool	Whether to overlay the projected control-volume boundaries on top of the filled region cross-section.	False
show_mesh_points	bool	Whether to overlay the selected mesh sites as points.	True
mesh_point_size	float	Scatter marker size for finalized mesh nodes.	12.0
mesh_facecolors	str	Scatter face color.	'white'
mesh_edgecolors	str	Scatter edge color.	'black'
mesh_linewidth	float	Scatter edge line width.	0.35
coordinate_scale	float	Multiplicative scale applied to displayed coordinates.	1.0
axis_unit	str | None	Optional axis-unit label appended to x/y/z labels.	None

Returns:

Type	Description
tuple[Figure, ndarray]	Tuple of (figure, axes).

Notes

For each requested plane, Tempura first selects the mesh sites whose control volumes intersect that plane. The optional marker overlay shows those selected sites explicitly. The colored background can either show the projected control volumes directly or a sampled cross-section of the finalized region shapes.

Source code in tempura/plotting/meshing.py

def plot_problem_regions_with_mesh(
    problem: Any,
    region_shapes: RegionShapeMap,
    plane_specs: Sequence[PlaneSpec],
    *,
    axes: Sequence[plt.Axes] | None = None,
    region_names: Sequence[str] | None = None,
    region_display_names: RegionDisplayMap | None = None,
    ncols: int | None = None,
    figsize: tuple[float, float] = (13.0, 11.0),
    suptitle: str = "Pescado region cuts with overlaid mesh nodes",
    fill_mode: str = "regions",
    background_alpha: float = 0.88,
    tile_edgecolors: str = "white",
    tile_linewidth: float = 0.6,
    show_volume_edges: bool = False,
    show_mesh_points: bool = True,
    mesh_point_size: float = 12.0,
    mesh_facecolors: str = "white",
    mesh_edgecolors: str = "black",
    mesh_linewidth: float = 0.35,
    coordinate_scale: float = 1.0,
    axis_unit: str | None = None,
) -> tuple[plt.Figure, np.ndarray]:
    """Plot finalized region slices with overlaid mesh nodes.

    Args:
        problem: Finalized Pescado problem with a ``coordinates`` array.
        region_shapes: Mapping returned by ``build_problem(...)``. These are the
            solver-owned region shapes after Tempura's lower-face trimming.
        plane_specs: Sequence of ``PlaneSpec`` describing the cuts to draw.
        region_names: Optional region order for coloring. Defaults to every
            region except ``"Boundary"`` in insertion order.
        region_display_names: Optional mapping from solver-owned region names to
            display labels. Regions that share a display label are drawn with
            the same color and collapsed into one legend entry.
        ncols: Number of subplot columns. Defaults to ``2`` when more than one
            plane is requested and ``1`` otherwise.
        figsize: Matplotlib figure size.
        suptitle: Figure title.
        fill_mode: ``"regions"`` to render the material cross-section from the
            finalized region shapes, or ``"control_volumes"`` to render the
            selected mesh control volumes directly.
        background_alpha: Alpha used for the ownership tiles under the markers.
        tile_edgecolors: Edge color for the projected mesh control volumes.
        tile_linewidth: Edge line width for the projected mesh control volumes.
        show_volume_edges: Whether to overlay the projected control-volume
            boundaries on top of the filled region cross-section.
        show_mesh_points: Whether to overlay the selected mesh sites as points.
        mesh_point_size: Scatter marker size for finalized mesh nodes.
        mesh_facecolors: Scatter face color.
        mesh_edgecolors: Scatter edge color.
        mesh_linewidth: Scatter edge line width.
        coordinate_scale: Multiplicative scale applied to displayed coordinates.
        axis_unit: Optional axis-unit label appended to x/y/z labels.

    Returns:
        Tuple of ``(figure, axes)``.

    Notes:
        For each requested plane, Tempura first selects the mesh sites whose
        control volumes intersect that plane. The optional marker overlay shows
        those selected sites explicitly. The colored background can either show
        the projected control volumes directly or a sampled cross-section of
        the finalized region shapes.
    """
    if not plane_specs:
        raise ValueError("plane_specs must not be empty.")
    if fill_mode not in {"regions", "control_volumes"}:
        raise ValueError(
            "fill_mode must be either 'regions' or 'control_volumes'; "
            f"got {fill_mode!r}."
        )
    if ncols is None:
        ncols = 2 if len(plane_specs) > 1 else 1
    if ncols <= 0:
        raise ValueError(f"ncols must be positive; got {ncols}.")

    coordinates = np.asarray(problem.coordinates, dtype=float)
    if coordinates.ndim != 2 or coordinates.shape[1] != 3:
        raise ValueError(
            "problem.coordinates must have shape (npoints, 3); "
            f"got {coordinates.shape}."
        )

    if region_names is None:
        region_names = [name for name in region_shapes if name != "Boundary"]
    else:
        region_names = list(region_names)
    if not region_names:
        raise ValueError("No region names were provided for plotting.")
    display_names, display_lookup = _display_order(region_names, region_display_names)

    volume_bounds = _volume_bounds_lookups(coordinates)
    if axes is None:
        nrows = int(np.ceil(len(plane_specs) / ncols))
        fig, axes_obj = plt.subplots(
            nrows, ncols, figsize=figsize, constrained_layout=True
        )
        axes_array = np.atleast_1d(axes_obj).ravel()
        created_axes = True
    else:
        axes_array = np.atleast_1d(np.asarray(list(axes), dtype=object)).ravel()
        if len(axes_array) < len(plane_specs):
            raise ValueError(
                f"axes must contain at least {len(plane_specs)} axes; got {len(axes_array)}."
            )
        fig = axes_array[0].figure
        created_axes = False

    for ax, spec in zip(axes_array, plane_specs, strict=False):
        plot_bbox = np.asarray(spec.plot_region.bbox, dtype=float) * float(
            coordinate_scale
        )
        plane_sites = _plane_sites(
            coordinates,
            spec=spec,
            bounds_lookup=volume_bounds[spec.plane_axis],
        )
        plane_points = np.asarray(plane_sites[:, spec.axis_indices], dtype=float) * float(
            coordinate_scale
        )
        labels = _label_plane_sites(
            plane_sites,
            region_shapes,
            region_names,
            region_display_names=region_display_names,
            display_lookup=display_lookup,
        )
        if fill_mode == "regions":
            _region_fill(
                ax,
                spec=spec,
                region_shapes=region_shapes,
                region_names=region_names,
                region_display_names=region_display_names,
                display_lookup=display_lookup,
                display_names=display_names,
                plane_sites=plane_sites,
                volume_bounds=volume_bounds,
                background_alpha=background_alpha,
                coordinate_scale=coordinate_scale,
            )
            if show_volume_edges:
                volume_edges = _volume_edges(
                    plane_sites,
                    spec=spec,
                    coordinate_scale=coordinate_scale,
                    volume_bounds=volume_bounds,
                    tile_edgecolors=tile_edgecolors,
                    tile_linewidth=tile_linewidth,
                )
                if volume_edges is not None:
                    ax.add_collection(volume_edges)
        else:
            ownership_tiles = _ownership_tiles(
                plane_sites,
                labels,
                spec=spec,
                coordinate_scale=coordinate_scale,
                background_alpha=background_alpha,
                volume_bounds=volume_bounds,
                tile_edgecolors=tile_edgecolors,
                tile_linewidth=tile_linewidth,
            )
            if ownership_tiles is not None:
                ax.add_collection(ownership_tiles)
        if show_mesh_points:
            ax.scatter(
                plane_points[:, 0],
                plane_points[:, 1],
                s=mesh_point_size,
                facecolors=mesh_facecolors,
                edgecolors=mesh_edgecolors,
                linewidths=mesh_linewidth,
                zorder=3,
            )
        ax.set_title(spec.title)
        ax.set_xlabel(_axis_label(spec.axis_labels[0], axis_unit))
        ax.set_ylabel(_axis_label(spec.axis_labels[1], axis_unit))
        ax.set_xlim(float(plot_bbox[0, 0]), float(plot_bbox[1, 0]))
        ax.set_ylim(float(plot_bbox[0, 1]), float(plot_bbox[1, 1]))
        ax.set_aspect("equal")

    if created_axes:
        for ax in axes_array[len(plane_specs) :]:
            ax.axis("off")

    fig.suptitle(suptitle, fontsize=16)
    fig.legend(
        handles=_legend_handles(display_names),
        loc="outside lower center",
        ncol=min(4, len(display_names)),
        frameon=False,
    )
    return fig, axes_array

plot_scalar_field_on_planes(problem, scalar_values, plane_specs, *, axes=None, ncols=None, nrows=None, figsize=(13.0, 11.0), suptitle='Scalar field cuts', colorbar_label='Value', cmap='coolwarm', symmetric=True, vmin=None, vmax=None, tile_edgecolors='none', tile_linewidth=0.0, show_volume_edges=False, aspect='equal', colorbar_shrink=1.0, coordinate_scale=1.0, axis_unit=None)

Plot one scalar field on a set of mesh-aligned 2D plane cuts.

Parameters:

Name	Type	Description	Default
problem	Any	Finalized Pescado problem with a coordinates array.	required
scalar_values	ndarray	One scalar value per solver mesh site.	required
plane_specs	Sequence[PlaneSpec]	Sequence of cuts to render.	required
axes	Sequence[Axes] | None	Optional existing axes that should receive the plots.	None
ncols	int | None	Number of subplot columns when axes are created internally.	None
nrows	int | None	Number of subplot rows when axes are created internally.	None
figsize	tuple[float, float]	Matplotlib figure size used when axes are created internally.	(13.0, 11.0)
suptitle	str	Figure title.	'Scalar field cuts'
colorbar_label	str	Label shown on the shared colorbar.	'Value'
cmap	str	Matplotlib colormap name.	'coolwarm'
symmetric	bool	Whether to choose symmetric color limits around zero when vmin and vmax are omitted.	True
vmin	float | None	Optional lower color limit.	None
vmax	float | None	Optional upper color limit.	None
tile_edgecolors	str	Edge color for optional projected control-volume outlines.	'none'
tile_linewidth	float	Line width for optional projected control-volume outlines.	0.0
show_volume_edges	bool	Whether to overlay projected control-volume edges.	False
aspect	str | float	Matplotlib aspect setting applied to each panel.	'equal'
colorbar_shrink	float	Shrink factor passed to the shared colorbar.	1.0
coordinate_scale	float	Multiplicative scale applied to displayed coordinates.	1.0
axis_unit	str | None	Optional axis-unit label appended to x/y/z labels.	None

Returns:

Type	Description
tuple[Figure, ndarray]	Tuple of (figure, axes) containing the rendered scalar-field cuts.

Source code in tempura/plotting/meshing.py

def plot_scalar_field_on_planes(
    problem: Any,
    scalar_values: np.ndarray,
    plane_specs: Sequence[PlaneSpec],
    *,
    axes: Sequence[plt.Axes] | None = None,
    ncols: int | None = None,
    nrows: int | None = None,
    figsize: tuple[float, float] = (13.0, 11.0),
    suptitle: str = "Scalar field cuts",
    colorbar_label: str = "Value",
    cmap: str = "coolwarm",
    symmetric: bool = True,
    vmin: float | None = None,
    vmax: float | None = None,
    tile_edgecolors: str = "none",
    tile_linewidth: float = 0.0,
    show_volume_edges: bool = False,
    aspect: str | float = "equal",
    colorbar_shrink: float = 1.0,
    coordinate_scale: float = 1.0,
    axis_unit: str | None = None,
) -> tuple[plt.Figure, np.ndarray]:
    """Plot one scalar field on a set of mesh-aligned 2D plane cuts.

    Args:
        problem: Finalized Pescado problem with a ``coordinates`` array.
        scalar_values: One scalar value per solver mesh site.
        plane_specs: Sequence of cuts to render.
        axes: Optional existing axes that should receive the plots.
        ncols: Number of subplot columns when axes are created internally.
        nrows: Number of subplot rows when axes are created internally.
        figsize: Matplotlib figure size used when axes are created internally.
        suptitle: Figure title.
        colorbar_label: Label shown on the shared colorbar.
        cmap: Matplotlib colormap name.
        symmetric: Whether to choose symmetric color limits around zero when
            ``vmin`` and ``vmax`` are omitted.
        vmin: Optional lower color limit.
        vmax: Optional upper color limit.
        tile_edgecolors: Edge color for optional projected control-volume
            outlines.
        tile_linewidth: Line width for optional projected control-volume
            outlines.
        show_volume_edges: Whether to overlay projected control-volume edges.
        aspect: Matplotlib aspect setting applied to each panel.
        colorbar_shrink: Shrink factor passed to the shared colorbar.
        coordinate_scale: Multiplicative scale applied to displayed
            coordinates.
        axis_unit: Optional axis-unit label appended to x/y/z labels.

    Returns:
        Tuple of ``(figure, axes)`` containing the rendered scalar-field cuts.
    """
    if not plane_specs:
        raise ValueError("plane_specs must not be empty.")
    if ncols is None and nrows is None:
        ncols = 2 if len(plane_specs) > 1 else 1
    elif ncols is None:
        if nrows is None or nrows <= 0:
            raise ValueError(f"nrows must be positive; got {nrows}.")
        ncols = int(np.ceil(len(plane_specs) / nrows))
    if ncols <= 0:
        raise ValueError(f"ncols must be positive; got {ncols}.")
    if nrows is not None and nrows <= 0:
        raise ValueError(f"nrows must be positive; got {nrows}.")

    coordinates = np.asarray(problem.coordinates, dtype=float)
    scalar_values = np.asarray(scalar_values, dtype=float)
    if scalar_values.shape != (coordinates.shape[0],):
        raise ValueError(
            "scalar_values must have shape (npoints,); "
            f"got {scalar_values.shape} for coordinates {coordinates.shape}."
        )

    volume_bounds = _volume_bounds_lookups(coordinates)
    panel_payloads: list[tuple[PlaneSpec, np.ndarray, np.ndarray]] = []
    panel_values: list[np.ndarray] = []
    for spec in plane_specs:
        indices = _exact_plane_site_indices(
            coordinates,
            spec=spec,
        )
        if indices.size == 0:
            indices = _plane_site_indices(
                coordinates,
                spec=spec,
                bounds_lookup=volume_bounds[spec.plane_axis],
            )
        sites = coordinates[indices] if indices.size else np.empty((0, 3), dtype=float)
        values = scalar_values[indices] if indices.size else np.empty((0,), dtype=float)
        panel_payloads.append((spec, sites, values))
        if values.size:
            panel_values.append(values)

    if panel_values:
        flattened = np.concatenate(panel_values)
        if symmetric:
            color_limit = float(np.max(np.abs(flattened)))
            if np.isclose(color_limit, 0.0):
                color_limit = 1.0
            resolved_vmin = -color_limit if vmin is None else float(vmin)
            resolved_vmax = color_limit if vmax is None else float(vmax)
        else:
            resolved_vmin = float(np.min(flattened) if vmin is None else vmin)
            resolved_vmax = float(np.max(flattened) if vmax is None else vmax)
            if np.isclose(resolved_vmin, resolved_vmax):
                resolved_vmin -= 1.0
                resolved_vmax += 1.0
    else:
        resolved_vmin = -1.0 if vmin is None else float(vmin)
        resolved_vmax = 1.0 if vmax is None else float(vmax)

    cmap_obj = plt.cm.get_cmap(cmap)
    norm = Normalize(vmin=resolved_vmin, vmax=resolved_vmax)
    if axes is None:
        if nrows is None:
            nrows = int(np.ceil(len(plane_specs) / ncols))
        fig, axes_obj = plt.subplots(
            nrows, ncols, figsize=figsize, constrained_layout=True
        )
        axes_array = np.atleast_1d(axes_obj).ravel()
        created_axes = True
    else:
        axes_array = np.atleast_1d(np.asarray(list(axes), dtype=object)).ravel()
        if len(axes_array) < len(plane_specs):
            raise ValueError(
                f"axes must contain at least {len(plane_specs)} axes; got {len(axes_array)}."
            )
        fig = axes_array[0].figure
        created_axes = False

    for ax, (spec, plane_sites, values) in zip(axes_array, panel_payloads, strict=False):
        plot_bbox = np.asarray(spec.plot_region.bbox, dtype=float) * float(
            coordinate_scale
        )
        plane_grid = _scalar_plane_grid(plane_sites, values, spec=spec)
        if plane_grid is not None:
            x_values, y_values, display_values = plane_grid
            ax.contourf(
                x_values * float(coordinate_scale),
                y_values * float(coordinate_scale),
                display_values,
                levels=np.linspace(resolved_vmin, resolved_vmax, 33),
                cmap=cmap_obj,
                norm=norm,
                antialiased=False,
                zorder=1,
            )
        if show_volume_edges:
            volume_edges = _volume_edges(
                plane_sites,
                spec=spec,
                coordinate_scale=coordinate_scale,
                volume_bounds=volume_bounds,
                tile_edgecolors=tile_edgecolors,
                tile_linewidth=tile_linewidth,
            )
            if volume_edges is not None:
                ax.add_collection(volume_edges)
        ax.set_title(spec.title)
        ax.set_xlabel(_axis_label(spec.axis_labels[0], axis_unit))
        ax.set_ylabel(_axis_label(spec.axis_labels[1], axis_unit))
        ax.set_xlim(float(plot_bbox[0, 0]), float(plot_bbox[1, 0]))
        ax.set_ylim(float(plot_bbox[0, 1]), float(plot_bbox[1, 1]))
        ax.set_aspect(aspect)

    if created_axes:
        for ax in axes_array[len(plane_specs) :]:
            ax.axis("off")

    sm = plt.cm.ScalarMappable(norm=norm, cmap=cmap_obj)
    sm.set_array([])
    fig.suptitle(suptitle, fontsize=16)
    fig.colorbar(
        sm,
        ax=list(axes_array[: len(plane_specs)]),
        label=colorbar_label,
        location="right",
        pad=0.03,
        fraction=0.05,
        shrink=float(colorbar_shrink),
    )
    return fig, axes_array

save_problem_regions_with_mesh(problem, region_shapes, plane_specs, output_path, *, region_names=None, region_display_names=None, ncols=None, figsize=(13.0, 11.0), suptitle='Pescado region cuts with overlaid mesh nodes', dpi=220, fill_mode='regions', background_alpha=0.88, tile_edgecolors='white', tile_linewidth=0.6, show_volume_edges=False, show_mesh_points=True, mesh_point_size=12.0, mesh_facecolors='white', mesh_edgecolors='black', mesh_linewidth=0.35, coordinate_scale=1.0, axis_unit=None)

Render region slices with overlaid mesh nodes and save the figure.

This is a convenience wrapper around plot_problem_regions_with_mesh(...) that also creates the parent directory and closes the Matplotlib figure after saving.

Parameters:

Name	Type	Description	Default
problem	Any	Finalized Pescado problem with a coordinates array.	required
region_shapes	RegionShapeMap	Mapping of region name to solver-owned shape.	required
plane_specs	Sequence[PlaneSpec]	Sequence of cuts to render.	required
output_path	str | Path	Destination image path.	required
region_names	Sequence[str] | None	Optional region order for coloring.	None
region_display_names	RegionDisplayMap | None	Optional mapping that merges several solver regions under one display label.	None
ncols	int | None	Number of subplot columns.	None
figsize	tuple[float, float]	Matplotlib figure size.	(13.0, 11.0)
suptitle	str	Figure title.	'Pescado region cuts with overlaid mesh nodes'
dpi	int	Output resolution passed to figure.savefig.	220
fill_mode	str	Either "regions" or "control_volumes".	'regions'
background_alpha	float	Alpha used for the filled region background.	0.88
tile_edgecolors	str	Edge color for projected control volumes.	'white'
tile_linewidth	float	Edge line width for projected control volumes.	0.6
show_volume_edges	bool	Whether to overlay projected control-volume edges.	False
show_mesh_points	bool	Whether to overlay the selected mesh sites as markers.	True
mesh_point_size	float	Scatter marker size.	12.0
mesh_facecolors	str	Scatter face color.	'white'
mesh_edgecolors	str	Scatter edge color.	'black'
mesh_linewidth	float	Scatter edge line width.	0.35
coordinate_scale	float	Multiplicative scale applied to displayed coordinates.	1.0
axis_unit	str | None	Optional axis-unit label appended to x/y/z labels.	None

Returns:

Type	Description
Path	Output path after the figure has been written.

Source code in tempura/plotting/meshing.py

def save_problem_regions_with_mesh(
    problem: Any,
    region_shapes: RegionShapeMap,
    plane_specs: Sequence[PlaneSpec],
    output_path: str | Path,
    *,
    region_names: Sequence[str] | None = None,
    region_display_names: RegionDisplayMap | None = None,
    ncols: int | None = None,
    figsize: tuple[float, float] = (13.0, 11.0),
    suptitle: str = "Pescado region cuts with overlaid mesh nodes",
    dpi: int = 220,
    fill_mode: str = "regions",
    background_alpha: float = 0.88,
    tile_edgecolors: str = "white",
    tile_linewidth: float = 0.6,
    show_volume_edges: bool = False,
    show_mesh_points: bool = True,
    mesh_point_size: float = 12.0,
    mesh_facecolors: str = "white",
    mesh_edgecolors: str = "black",
    mesh_linewidth: float = 0.35,
    coordinate_scale: float = 1.0,
    axis_unit: str | None = None,
) -> Path:
    """Render region slices with overlaid mesh nodes and save the figure.

    This is a convenience wrapper around ``plot_problem_regions_with_mesh(...)``
    that also creates the parent directory and closes the Matplotlib figure
    after saving.

    Args:
        problem: Finalized Pescado problem with a ``coordinates`` array.
        region_shapes: Mapping of region name to solver-owned shape.
        plane_specs: Sequence of cuts to render.
        output_path: Destination image path.
        region_names: Optional region order for coloring.
        region_display_names: Optional mapping that merges several solver
            regions under one display label.
        ncols: Number of subplot columns.
        figsize: Matplotlib figure size.
        suptitle: Figure title.
        dpi: Output resolution passed to ``figure.savefig``.
        fill_mode: Either ``"regions"`` or ``"control_volumes"``.
        background_alpha: Alpha used for the filled region background.
        tile_edgecolors: Edge color for projected control volumes.
        tile_linewidth: Edge line width for projected control volumes.
        show_volume_edges: Whether to overlay projected control-volume edges.
        show_mesh_points: Whether to overlay the selected mesh sites as
            markers.
        mesh_point_size: Scatter marker size.
        mesh_facecolors: Scatter face color.
        mesh_edgecolors: Scatter edge color.
        mesh_linewidth: Scatter edge line width.
        coordinate_scale: Multiplicative scale applied to displayed
            coordinates.
        axis_unit: Optional axis-unit label appended to x/y/z labels.

    Returns:
        Output path after the figure has been written.
    """
    output_path = Path(output_path)
    figure, _ = plot_problem_regions_with_mesh(
        problem,
        region_shapes,
        plane_specs,
        region_names=region_names,
        region_display_names=region_display_names,
        ncols=ncols,
        figsize=figsize,
        suptitle=suptitle,
        fill_mode=fill_mode,
        background_alpha=background_alpha,
        tile_edgecolors=tile_edgecolors,
        tile_linewidth=tile_linewidth,
        show_volume_edges=show_volume_edges,
        show_mesh_points=show_mesh_points,
        mesh_point_size=mesh_point_size,
        mesh_facecolors=mesh_facecolors,
        mesh_edgecolors=mesh_edgecolors,
        mesh_linewidth=mesh_linewidth,
        coordinate_scale=coordinate_scale,
        axis_unit=axis_unit,
    )
    output_path.parent.mkdir(parents=True, exist_ok=True)
    figure.savefig(output_path, dpi=dpi)
    plt.close(figure)
    return output_path

transpose_plane_spec(spec, *, title=None)

Return one plane spec with the displayed axes swapped.

Parameters:

Name	Type	Description	Default
spec	PlaneSpec	Original plane specification.	required
title	str | None	Optional replacement title for the transposed panel.	None

Returns:

Type	Description
PlaneSpec	Copy of spec with its 2D plot region, axis indices, and axis labels
PlaneSpec	reversed.

Source code in tempura/plotting/meshing.py

def transpose_plane_spec(spec: PlaneSpec, *, title: str | None = None) -> PlaneSpec:
    """Return one plane spec with the displayed axes swapped.

    Args:
        spec: Original plane specification.
        title: Optional replacement title for the transposed panel.

    Returns:
        Copy of ``spec`` with its 2D plot region, axis indices, and axis labels
        reversed.
    """
    bbox = np.asarray(spec.plot_region.bbox, dtype=float)
    lower_left = np.asarray(bbox[0, ::-1], dtype=float)
    size = np.asarray((bbox[1] - bbox[0])[::-1], dtype=float)
    return replace(
        spec,
        title=spec.title if title is None else title,
        plot_region=shapes.Box(lower_left=lower_left, size=size),
        axis_indices=tuple(reversed(spec.axis_indices)),
        axis_labels=tuple(reversed(spec.axis_labels)),
    )