Skip to content

Visualization

The voids.visualization sub-package provides network rendering via Plotly and PyVista, plus reusable scalar/vector field plotting and ParaView export helpers for map-based validation workflows.

PyVista dependency

PyVista is installed as a core dependency of voids and is available by default when you install the package.

Plotly

voids.visualization.plotly

plot_network_plotly 

plot_network_plotly(
    net,
    *,
    point_scalars=None,
    cell_scalars=None,
    point_sizes=None,
    throat_sizes=None,
    point_size=None,
    line_width=None,
    line_opacity=0.4,
    size_scale=1.0,
    point_size_limits=None,
    throat_size_limits=None,
    max_throats=1000,
    title=None,
    show_colorbar=True,
    layout_kwargs=None,
)

Create an interactive Plotly visualization of a pore-throat network.

Parameters:

Name Type Description Default
net Network

Network to render.

 required
point_scalars str | ndarray | None

Pore field name or explicit pore-valued array with shape (Np,).

 None
cell_scalars str | ndarray | None

Throat field name or explicit throat-valued array with shape (Nt,).

 None
point_sizes str | ndarray | bool | None

Pore/throat characteristic size field name, explicit size array, None for automatic size-field detection, or False to disable size-driven rendering. Automatically detected size fields follow the priority diameter_equivalent -> diameter_inscribed -> radius_inscribed -> area.

 None
throat_sizes str | ndarray | bool | None

Pore/throat characteristic size field name, explicit size array, None for automatic size-field detection, or False to disable size-driven rendering. Automatically detected size fields follow the priority diameter_equivalent -> diameter_inscribed -> radius_inscribed -> area.

 None
point_size float | None

Constant marker size for pores when explicit size rendering is disabled. When size-driven rendering is active, this acts as the reference marker size for median-sized pores.

 None
line_width float | None

Constant line width for throats when explicit size rendering is disabled. When size-driven rendering is active, this acts as the reference width for median-sized throats.

 None
line_opacity float

Opacity applied to throat lines.

 0.4
size_scale float

Multiplicative factor applied to size-driven pore markers and throat widths.

 1.0
point_size_limits tuple[float | None, float | None] | None

Optional (min_px, max_px) limits for size-driven rendering in screen-space pixels. When omitted, conservative default clipping is applied for readability. Set to (None, None) to disable clipping and preserve the full relative dynamic range.

 None
throat_size_limits tuple[float | None, float | None] | None

Optional (min_px, max_px) limits for size-driven rendering in screen-space pixels. When omitted, conservative default clipping is applied for readability. Set to (None, None) to disable clipping and preserve the full relative dynamic range.

 None
max_throats int | None

Maximum number of throats to draw. Large networks are downsampled for responsiveness.

 1000
title str | None

Figure title.

 None
show_colorbar bool

If True, display a colorbar for pore scalars.

 True
layout_kwargs dict[str, Any] | None

Optional Plotly layout overrides.

 None

Returns:

Type Description
Figure

Interactive 3D figure.
Notes

If only pore scalars are given, each throat is colored by the arithmetic mean of its endpoint pore values:

s_throat = 0.5 * (s_i + s_j)

The throat colormap is normalized with the same scalar bounds as the pore markers so that equal numerical values map to equal colors across pores and throats.

Source code in src/voids/visualization/plotly.py 
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
def plot_network_plotly(
    net: Network,
    *,
    point_scalars: str | np.ndarray | None = None,
    cell_scalars: str | np.ndarray | None = None,
    point_sizes: str | np.ndarray | bool | None = None,
    throat_sizes: str | np.ndarray | bool | None = None,
    point_size: float | None = None,
    line_width: float | None = None,
    line_opacity: float = 0.4,
    size_scale: float = 1.0,
    point_size_limits: tuple[float | None, float | None] | None = None,
    throat_size_limits: tuple[float | None, float | None] | None = None,
    max_throats: int | None = 1000,
    title: str | None = None,
    show_colorbar: bool = True,
    layout_kwargs: dict[str, Any] | None = None,
) -> go.Figure:
    """Create an interactive Plotly visualization of a pore-throat network.

    Parameters
    ----------
    net :
        Network to render.
    point_scalars :
        Pore field name or explicit pore-valued array with shape ``(Np,)``.
    cell_scalars :
        Throat field name or explicit throat-valued array with shape ``(Nt,)``.
    point_sizes, throat_sizes :
        Pore/throat characteristic size field name, explicit size array, ``None`` for
        automatic size-field detection, or ``False`` to disable size-driven rendering.
        Automatically detected size fields follow the priority
        ``diameter_equivalent -> diameter_inscribed -> radius_inscribed -> area``.
    point_size :
        Constant marker size for pores when explicit size rendering is disabled. When
        size-driven rendering is active, this acts as the reference marker size for
        median-sized pores.
    line_width :
        Constant line width for throats when explicit size rendering is disabled. When
        size-driven rendering is active, this acts as the reference width for
        median-sized throats.
    line_opacity :
        Opacity applied to throat lines.
    size_scale :
        Multiplicative factor applied to size-driven pore markers and throat widths.
    point_size_limits, throat_size_limits :
        Optional ``(min_px, max_px)`` limits for size-driven rendering in screen-space
        pixels. When omitted, conservative default clipping is applied for readability.
        Set to ``(None, None)`` to disable clipping and preserve the full relative
        dynamic range.
    max_throats :
        Maximum number of throats to draw. Large networks are downsampled for responsiveness.
    title :
        Figure title.
    show_colorbar :
        If ``True``, display a colorbar for pore scalars.
    layout_kwargs :
        Optional Plotly layout overrides.

    Returns
    -------
    plotly.graph_objects.Figure
        Interactive 3D figure.

    Notes
    -----
    If only pore scalars are given, each throat is colored by the arithmetic mean of
    its endpoint pore values:

    ``s_throat = 0.5 * (s_i + s_j)``

    The throat colormap is normalized with the same scalar bounds as the pore markers
    so that equal numerical values map to equal colors across pores and throats.
    """

    point_values, point_label = _resolve_scalars(
        point_scalars, store=net.pore, expected_shape=(net.Np,), prefix="pore"
    )
    cell_values, cell_label = _resolve_scalars(
        cell_scalars, store=net.throat, expected_shape=(net.Nt,), prefix="throat"
    )
    point_size_values, point_size_label = resolve_size_values(
        point_sizes, store=net.pore, expected_shape=(net.Np,), prefix="pore"
    )
    throat_size_values, throat_size_label = resolve_size_values(
        throat_sizes, store=net.throat, expected_shape=(net.Nt,), prefix="throat"
    )

    coords = np.asarray(net.pore_coords, dtype=float)
    x, y, z = coords[:, 0], coords[:, 1], coords[:, 2]
    sampled = _sample_indices(net.Nt, max_throats)
    point_vmin, point_vmax = _scalar_bounds(point_values)
    cell_vmin, cell_vmax = _scalar_bounds(cell_values)
    point_size_ref = float(
        point_size if point_size is not None else (6.0 if net.Np <= 2000 else 4.0)
    )
    line_width_ref = float(2.0 if line_width is None else line_width)
    use_variable_point_sizes = point_size_values is not None
    use_variable_throat_sizes = throat_size_values is not None
    point_min_size: float | None
    point_max_size: float | None
    if point_size_limits is None:
        point_min_size = max(2.0, 0.5 * point_size_ref)
        point_max_size = max(18.0, 4.0 * point_size_ref)
    else:
        point_min_size, point_max_size = point_size_limits
    throat_min_size: float | None
    throat_max_size: float | None
    if throat_size_limits is None:
        throat_min_size = 0.75
        throat_max_size = max(10.0, 4.0 * line_width_ref)
    else:
        throat_min_size, throat_max_size = throat_size_limits

    if use_variable_point_sizes:
        assert point_size_values is not None
        marker_size: float | np.ndarray = scale_sizes_to_pixels(
            point_size_values,
            reference=point_size_ref,
            scale=size_scale,
            min_size=point_min_size,
            max_size=point_max_size,
        )
    else:
        marker_size = point_size_ref

    if use_variable_throat_sizes:
        assert throat_size_values is not None
        sampled_line_widths = scale_sizes_to_pixels(
            throat_size_values[sampled],
            reference=line_width_ref,
            scale=size_scale,
            min_size=throat_min_size,
            max_size=throat_max_size,
        )
    else:
        sampled_line_widths = np.full(sampled.shape, line_width_ref, dtype=float)

    if point_values is not None:
        marker: dict[str, Any] = {
            "size": marker_size,
            "color": point_values,
            "colorscale": "Viridis",
            "showscale": show_colorbar,
        }
        if use_variable_point_sizes:
            marker["sizemode"] = "diameter"
        if point_vmin is not None and point_vmax is not None:
            marker["cmin"] = point_vmin
            marker["cmax"] = point_vmax
        if show_colorbar:
            marker["colorbar"] = {"title": point_label or "pore scalar"}
    else:
        marker = {
            "size": marker_size,
            "color": "royalblue",
            "showscale": False,
        }
        if use_variable_point_sizes:
            marker["sizemode"] = "diameter"

    pore_text = []
    for idx in range(net.Np):
        hover_lines = [f"Pore {idx}"]
        if point_values is not None:
            hover_lines.append(f"{point_label or 'value'}={point_values[idx]:.3e}")
        if use_variable_point_sizes and point_size_label is not None:
            assert point_size_values is not None
            hover_lines.append(f"{point_size_label}={point_size_values[idx]:.3e}")
        pore_text.append("<br>".join(hover_lines))

    traces: list[Any] = [
        go.Scatter3d(
            x=x,
            y=y,
            z=z,
            mode="markers",
            marker=marker,
            name="Pores",
            text=pore_text,
            hoverinfo="text",
        )
    ]

    if cell_values is not None:
        throat_values = np.asarray(cell_values[sampled], dtype=float)
        throat_label = cell_label or "throat scalar"
        color_vmin, color_vmax = cell_vmin, cell_vmax
    elif point_values is not None:
        conns = net.throat_conns[sampled]
        throat_values = 0.5 * (point_values[conns[:, 0]] + point_values[conns[:, 1]])
        throat_label = f"avg({point_label or 'pore scalar'})"
        color_vmin, color_vmax = point_vmin, point_vmax
    else:
        throat_values = None
        throat_label = None
        color_vmin, color_vmax = None, None

    for local_idx, throat_idx in enumerate(sampled):
        i, j = net.throat_conns[throat_idx]
        hover_lines = [f"Throat {int(throat_idx)}"]
        if throat_values is None:
            color = _rgb_with_opacity("rgb(100,100,100)", line_opacity)
        else:
            if color_vmin is not None and color_vmax is not None and color_vmax > color_vmin:
                norm = float((throat_values[local_idx] - color_vmin) / (color_vmax - color_vmin))
            else:
                norm = 0.5
            color = _rgb_with_opacity(sample_colorscale("Viridis", [norm])[0], line_opacity)
            hover_lines.append(f"{throat_label}={float(throat_values[local_idx]):.3e}")
        if use_variable_throat_sizes and throat_size_label is not None:
            assert throat_size_values is not None
            hover_lines.append(f"{throat_size_label}={float(throat_size_values[throat_idx]):.3e}")
        traces.append(
            go.Scatter3d(
                x=[x[i], x[j]],
                y=[y[i], y[j]],
                z=[z[i], z[j]],
                mode="lines",
                line={"color": color, "width": float(sampled_line_widths[local_idx])},
                name="Throats",
                showlegend=False,
                text="<br>".join(hover_lines),
                hoverinfo="text",
            )
        )

    if title is None:
        title = "Pore network"
        if sampled.size != net.Nt:
            title += f" (showing {sampled.size} of {net.Nt} throats)"

    figure = go.Figure(data=traces)
    layout: dict[str, Any] = {
        "title": title,
        "scene": {
            "xaxis_title": "X",
            "yaxis_title": "Y",
            "zaxis_title": "Z",
            "aspectmode": "data",
        },
        "width": 900,
        "height": 700,
        "hovermode": "closest",
    }
    if layout_kwargs:
        layout.update(layout_kwargs)
    figure.update_layout(**layout)
    return figure

PyVista

voids.visualization.pyvista

network_to_pyvista_polydata 

network_to_pyvista_polydata(
    net,
    *,
    point_scalars=None,
    cell_scalars=None,
    include_all_numeric_fields=False,
)

Convert a network to pyvista.PolyData.

Parameters:

Name Type Description Default
net Network

Network to convert.

 required
point_scalars str | ndarray | None

Pore/throat scalar field name or explicit array.

 None
cell_scalars str | ndarray | None

Pore/throat scalar field name or explicit array.

 None
include_all_numeric_fields bool

If True, attach every 1D numeric pore/throat array whose length matches Np or Nt.

 False

Returns:

Type Description
PolyData

PolyData with pores as points and throats as line cells.

Raises:

Type Description
KeyError

If a requested scalar field name is missing.
ValueError

If an explicit scalar array has the wrong shape.
Source code in src/voids/visualization/pyvista.py 
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
def network_to_pyvista_polydata(
    net: Network,
    *,
    point_scalars: str | np.ndarray | None = None,
    cell_scalars: str | np.ndarray | None = None,
    include_all_numeric_fields: bool = False,
) -> pv.PolyData:
    """Convert a network to ``pyvista.PolyData``.

    Parameters
    ----------
    net :
        Network to convert.
    point_scalars, cell_scalars :
        Pore/throat scalar field name or explicit array.
    include_all_numeric_fields :
        If ``True``, attach every 1D numeric pore/throat array whose length matches
        ``Np`` or ``Nt``.

    Returns
    -------
    pyvista.PolyData
        PolyData with pores as points and throats as line cells.

    Raises
    ------
    KeyError
        If a requested scalar field name is missing.
    ValueError
        If an explicit scalar array has the wrong shape.
    """

    points = np.asarray(net.pore_coords, dtype=float)
    line_cells = _line_cells_from_conns(net.throat_conns)
    poly: pv.PolyData = pv.PolyData(points, lines=line_cells)

    poly.point_data["pore.id"] = np.arange(net.Np, dtype=np.int64)
    poly.cell_data["throat.id"] = np.arange(net.Nt, dtype=np.int64)

    if include_all_numeric_fields:
        for k, v in net.pore.items():
            a = np.asarray(v)
            if a.ndim == 1 and a.shape[0] == net.Np and np.issubdtype(a.dtype, np.number):
                poly.point_data[f"pore.{k}"] = a
        for k, v in net.throat.items():
            a = np.asarray(v)
            if a.ndim == 1 and a.shape[0] == net.Nt and np.issubdtype(a.dtype, np.number):
                poly.cell_data[f"throat.{k}"] = a

    if isinstance(point_scalars, str):
        if point_scalars not in net.pore:
            raise KeyError(f"Missing pore field '{point_scalars}'")
        poly.point_data["pore.scalar"] = np.asarray(net.pore[point_scalars])
        poly.set_active_scalars("pore.scalar", preference="point")
    elif point_scalars is not None:
        arr = np.asarray(point_scalars)
        if arr.shape != (net.Np,):
            raise ValueError("point_scalars array must have shape (Np,)")
        poly.point_data["pore.scalar"] = arr
        poly.set_active_scalars("pore.scalar", preference="point")

    if isinstance(cell_scalars, str):
        if cell_scalars not in net.throat:
            raise KeyError(f"Missing throat field '{cell_scalars}'")
        poly.cell_data["throat.scalar"] = np.asarray(net.throat[cell_scalars])
    elif cell_scalars is not None:
        arr = np.asarray(cell_scalars)
        if arr.shape != (net.Nt,):
            raise ValueError("cell_scalars array must have shape (Nt,)")
        poly.cell_data["throat.scalar"] = arr

    return poly

plot_network_pyvista 

plot_network_pyvista(
    net,
    *,
    point_scalars=None,
    cell_scalars=None,
    point_sizes=None,
    throat_sizes=None,
    show_points=True,
    show_lines=True,
    line_width=None,
    point_size=None,
    render_tubes=False,
    tube_radius=None,
    size_scale=1.0,
    off_screen=False,
    screenshot=None,
    show_axes=True,
    notebook=None,
    **add_mesh_kwargs,
)

Render a pore network with PyVista.

Parameters:

Name Type Description Default
net Network

Network to render.

 required
point_scalars str | ndarray | None

Pore/throat scalar field name or explicit array.

 None
cell_scalars str | ndarray | None

Pore/throat scalar field name or explicit array.

 None
point_sizes str | ndarray | bool | None

Pore/throat characteristic size field name, explicit size array, None for automatic size-field detection, or False to disable size-driven rendering. Automatically detected size fields follow the priority diameter_equivalent -> diameter_inscribed -> radius_inscribed -> area.

 None
throat_sizes str | ndarray | bool | None

Pore/throat characteristic size field name, explicit size array, None for automatic size-field detection, or False to disable size-driven rendering. Automatically detected size fields follow the priority diameter_equivalent -> diameter_inscribed -> radius_inscribed -> area.

 None
show_points bool

Toggle pore and throat rendering.

 True
show_lines bool

Toggle pore and throat rendering.

 True
line_width float | None

Width used for line rendering when throat sizes are not rendered with tubes.

 None
point_size float | None

Marker size used for pores when size-driven sphere rendering is disabled.

 None
render_tubes bool

If True, convert throats from lines to tubes.

 False
tube_radius float | None

Optional tube radius when render_tubes is enabled.

 None
size_scale float

Multiplicative factor applied to point diameters and throat radii in world units.

 1.0
off_screen bool

If True, create an off-screen plotter for headless rendering.

 False
screenshot str | None

Optional screenshot output path. When provided, the plot is rendered and saved.

 None
show_axes bool

If True, display orientation axes.

 True
notebook bool | None

Optional PyVista notebook flag. Defaults to False when omitted.

 None
**add_mesh_kwargs Any

Additional keyword arguments forwarded to :meth:pyvista.Plotter.add_mesh.

 {}

Returns:

Type Description
tuple

Pair (plotter, polydata).
Notes

For throat rendering, scalar selection follows this priority:

  1. explicit throat/cell scalar data
  2. pore/point scalar data, reused on the line representation

This allows pore-defined pressure fields to color both pores and throats in a consistent network visualization.

When tube rendering is requested (via render_tubes, tube_radius, or automatically detected variable throat sizes) but the PyVista tube filter is unavailable or raises an exception, rendering falls back to line geometry with render_lines_as_tubes=True so that a tube-like appearance is preserved. Variable throat radii cannot be accurately represented in this fallback mode; a :class:UserWarning is emitted to notify callers.

Source code in src/voids/visualization/pyvista.py 
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
def plot_network_pyvista(
    net: Network,
    *,
    point_scalars: str | np.ndarray | None = None,
    cell_scalars: str | np.ndarray | None = None,
    point_sizes: str | np.ndarray | bool | None = None,
    throat_sizes: str | np.ndarray | bool | None = None,
    show_points: bool = True,
    show_lines: bool = True,
    line_width: float | None = None,
    point_size: float | None = None,
    render_tubes: bool = False,
    tube_radius: float | None = None,
    size_scale: float = 1.0,
    off_screen: bool = False,
    screenshot: str | None = None,
    show_axes: bool = True,
    notebook: bool | None = None,
    **add_mesh_kwargs: Any,
) -> tuple[pv.Plotter, pv.PolyData]:
    """Render a pore network with PyVista.

    Parameters
    ----------
    net :
        Network to render.
    point_scalars, cell_scalars :
        Pore/throat scalar field name or explicit array.
    point_sizes, throat_sizes :
        Pore/throat characteristic size field name, explicit size array, ``None`` for
        automatic size-field detection, or ``False`` to disable size-driven rendering.
        Automatically detected size fields follow the priority
        ``diameter_equivalent -> diameter_inscribed -> radius_inscribed -> area``.
    show_points, show_lines :
        Toggle pore and throat rendering.
    line_width :
        Width used for line rendering when throat sizes are not rendered with tubes.
    point_size :
        Marker size used for pores when size-driven sphere rendering is disabled.
    render_tubes :
        If ``True``, convert throats from lines to tubes.
    tube_radius :
        Optional tube radius when ``render_tubes`` is enabled.
    size_scale :
        Multiplicative factor applied to point diameters and throat radii in world units.
    off_screen :
        If ``True``, create an off-screen plotter for headless rendering.
    screenshot :
        Optional screenshot output path. When provided, the plot is rendered and saved.
    show_axes :
        If ``True``, display orientation axes.
    notebook :
        Optional PyVista notebook flag. Defaults to ``False`` when omitted.
    **add_mesh_kwargs :
        Additional keyword arguments forwarded to :meth:`pyvista.Plotter.add_mesh`.

    Returns
    -------
    tuple
        Pair ``(plotter, polydata)``.

    Notes
    -----
    For throat rendering, scalar selection follows this priority:

    1. explicit throat/cell scalar data
    2. pore/point scalar data, reused on the line representation

    This allows pore-defined pressure fields to color both pores and throats in a
    consistent network visualization.

    When tube rendering is requested (via ``render_tubes``, ``tube_radius``, or
    automatically detected variable throat sizes) but the PyVista tube filter is
    unavailable or raises an exception, rendering falls back to line geometry with
    ``render_lines_as_tubes=True`` so that a tube-like appearance is preserved.
    Variable throat radii **cannot** be accurately represented in this fallback mode;
    a :class:`UserWarning` is emitted to notify callers.
    """

    poly = network_to_pyvista_polydata(
        net,
        point_scalars=point_scalars,
        cell_scalars=cell_scalars,
        include_all_numeric_fields=True,
    )

    if notebook is None:
        notebook = False
    pl = pv.Plotter(off_screen=off_screen, notebook=notebook)

    line_scalars_name = "throat.scalar" if "throat.scalar" in poly.cell_data else None
    point_scalars_name = "pore.scalar" if "pore.scalar" in poly.point_data else None
    if line_scalars_name is None and point_scalars_name is not None:
        line_scalars_name = point_scalars_name
    point_size_values, _ = resolve_size_values(
        point_sizes, store=net.pore, expected_shape=(net.Np,), prefix="pore"
    )
    throat_size_values, _ = resolve_size_values(
        throat_sizes, store=net.throat, expected_shape=(net.Nt,), prefix="throat"
    )
    use_variable_point_sizes = point_size_values is not None
    use_variable_throat_sizes = throat_size_values is not None
    point_size_value = float(point_size if point_size is not None else 9.0)
    line_width_value = float(3.0 if line_width is None else line_width)
    if use_variable_point_sizes:
        poly.point_data["pore.render_diameter"] = float(size_scale) * np.asarray(
            point_size_values, dtype=float
        )
    if use_variable_throat_sizes:
        poly.cell_data["throat.render_radius"] = (
            0.5 * float(size_scale) * np.asarray(throat_size_values, dtype=float)
        )

    if show_lines and net.Nt > 0:
        line_mesh = poly
        render_tubes_effective = (
            render_tubes or use_variable_throat_sizes or tube_radius is not None
        )
        if render_tubes_effective:
            kwargs: dict[str, Any] = {}
            if use_variable_throat_sizes:
                kwargs["scalars"] = "throat.render_radius"
                kwargs["absolute"] = True
                kwargs["preference"] = "cell"
            elif tube_radius is not None:
                kwargs["radius"] = float(tube_radius)
            try:
                line_mesh = poly.tube(**kwargs)
            except Exception as exc:
                line_mesh = poly
                msg = (
                    f"PyVista tube filter failed ({type(exc).__name__}: {exc}); "
                    "falling back to line rendering with render_lines_as_tubes=True."
                )
                if use_variable_throat_sizes:
                    msg += (
                        " Variable throat radii cannot be represented accurately "
                        "without the tube filter."
                    )
                warnings.warn(msg, stacklevel=2)
        line_kwargs: dict[str, Any] = {
            "scalars": line_scalars_name,
            "show_scalar_bar": (line_scalars_name is not None),
            **add_mesh_kwargs,
        }
        if line_mesh is poly:
            line_kwargs["line_width"] = line_width_value
            # Use line-tube approximation when tubes were requested but unavailable.
            line_kwargs["render_lines_as_tubes"] = render_tubes_effective
        pl.add_mesh(line_mesh, **line_kwargs)

    if show_points and net.Np > 0:
        if use_variable_point_sizes:
            point_mesh = pv.PolyData(np.asarray(net.pore_coords, dtype=float))
            point_mesh.point_data["pore.render_diameter"] = poly.point_data["pore.render_diameter"]
            if point_scalars_name is not None:
                point_mesh.point_data[point_scalars_name] = poly.point_data[point_scalars_name]
            sphere = pv.Sphere(radius=0.5)
            point_mesh = point_mesh.glyph(
                scale="pore.render_diameter",
                orient=False,
                factor=1.0,
                geom=sphere,
            )
            pl.add_mesh(
                point_mesh,
                scalars=point_scalars_name,
                show_scalar_bar=(point_scalars_name is not None and not show_lines),
            )
        else:
            pl.add_mesh(
                poly,
                style="points",
                point_size=point_size_value,
                render_points_as_spheres=True,
                scalars=point_scalars_name,
                show_scalar_bar=(point_scalars_name is not None and not show_lines),
            )

    if show_axes:
        pl.add_axes()  # type: ignore[call-arg]

    if screenshot is not None:
        pl.show(auto_close=False)
        pl.screenshot(screenshot)
    return pl, poly

Field Plots And Exports

The field helpers are intended for pressure and velocity diagnostics generated by finite-volume, finite-element, and lattice-Boltzmann single-phase workflows. Velocity midplane plots draw velocity magnitude as the scalar background and in-plane quiver arrows as the vector overlay. Structured vector fields can be written as VTU/VTK cell data, and DOLFINx functions can be written as XDMF/HDF5 after interpolation to first-order visualization spaces for ParaView compatibility. For pressure diagnostics, reference_pressure_to_outlet removes arbitrary pressure gauges by shifting the field so the outlet layer has a prescribed reference pressure, while preserving all pressure differences.

voids.visualization.fields

vector_magnitude 

vector_magnitude(vector)

Return the Euclidean magnitude of a dim-first vector field.

Source code in src/voids/visualization/fields.py 
15
16
17
18
19
def vector_magnitude(vector: np.ndarray) -> np.ndarray:
    """Return the Euclidean magnitude of a dim-first vector field."""

    vec = _as_vector_field(vector)
    return np.asarray(np.linalg.norm(vec, axis=0), dtype=float)

reference_pressure_to_outlet 

reference_pressure_to_outlet(
    pressure,
    *,
    flow_axis="x",
    reference_pressure=100000.0,
    pressure_outlet=0.0,
)

Shift a pressure field so the outlet layer has a reference pressure.

In incompressible Darcy, Brinkman, and Stokes solves, pressure is only determined up to an additive constant. This helper preserves the computed pressure differences and adds the constant needed to make the mean outlet layer pressure equal to reference_pressure + pressure_outlet.

The outlet layer is assumed to be the maximum-index face along flow_axis, matching the pressure convention used by the map-based validation notebooks.

Source code in src/voids/visualization/fields.py 
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
def reference_pressure_to_outlet(
    pressure: np.ndarray,
    *,
    flow_axis: str = "x",
    reference_pressure: float = 1.0e5,
    pressure_outlet: float = 0.0,
) -> np.ndarray:
    """Shift a pressure field so the outlet layer has a reference pressure.

    In incompressible Darcy, Brinkman, and Stokes solves, pressure is only
    determined up to an additive constant. This helper preserves the computed
    pressure differences and adds the constant needed to make the mean outlet
    layer pressure equal to ``reference_pressure + pressure_outlet``.

    The outlet layer is assumed to be the maximum-index face along
    ``flow_axis``, matching the pressure convention used by the map-based
    validation notebooks.
    """

    values = _as_scalar_field(pressure)
    target_outlet_pressure = float(reference_pressure) + float(pressure_outlet)
    if not np.isfinite(target_outlet_pressure):
        raise ValueError("reference and outlet pressures must be finite")
    axis_index = _axis_index(flow_axis, values.ndim)
    outlet_selector: list[slice | int] = [slice(None)] * values.ndim
    outlet_selector[axis_index] = values.shape[axis_index] - 1
    outlet_reference = float(np.mean(values[tuple(outlet_selector)]))
    return values + target_outlet_pressure - outlet_reference

reconstruct_tpfa_cell_velocity 

reconstruct_tpfa_cell_velocity(
    pressure,
    permeability,
    *,
    flow_axis="x",
    viscosity=1.0,
    pressure_inlet=1.0,
    pressure_outlet=0.0,
    cell_size=None,
)

Reconstruct a cell-centered Darcy velocity field from a TPFA pressure solve.

The reconstruction uses the same two-point face fluxes as the TPFA balance: harmonic permeability on interior faces, half-cell Dirichlet transmissibility on inlet/outlet faces, and zero transverse boundary flux. The returned field is a dim-first array with shape (ndim, *pressure.shape).

Source code in src/voids/visualization/fields.py 
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
def reconstruct_tpfa_cell_velocity(
    pressure: np.ndarray,
    permeability: PermeabilityMap | np.ndarray,
    *,
    flow_axis: str = "x",
    viscosity: float = 1.0,
    pressure_inlet: float = 1.0,
    pressure_outlet: float = 0.0,
    cell_size: float | Sequence[float] | None = None,
) -> np.ndarray:
    """Reconstruct a cell-centered Darcy velocity field from a TPFA pressure solve.

    The reconstruction uses the same two-point face fluxes as the TPFA balance:
    harmonic permeability on interior faces, half-cell Dirichlet transmissibility
    on inlet/outlet faces, and zero transverse boundary flux. The returned field
    is a dim-first array with shape ``(ndim, *pressure.shape)``.
    """

    pressure_values = np.asarray(pressure, dtype=float)
    if pressure_values.ndim not in {2, 3}:
        raise ValueError("pressure must be a 2D or 3D cell-centered field")
    if not np.all(np.isfinite(pressure_values)):
        raise ValueError("pressure must contain only finite values")
    permeability_values, size = _permeability_values_and_cell_size(
        permeability,
        pressure_values.shape,
        cell_size=cell_size,
    )
    if viscosity <= 0.0 or not np.isfinite(viscosity):
        raise ValueError("viscosity must be positive and finite")
    flow_axis_index = _axis_index(flow_axis, pressure_values.ndim)

    velocity = np.zeros((pressure_values.ndim, *pressure_values.shape), dtype=float)
    for direction in range(pressure_values.ndim):
        low_face = np.zeros_like(pressure_values, dtype=float)
        high_face = np.zeros_like(pressure_values, dtype=float)
        spacing = float(size[direction])

        low_selector: list[slice | int] = [slice(None)] * pressure_values.ndim
        high_selector: list[slice | int] = [slice(None)] * pressure_values.ndim
        low_selector[direction] = slice(0, -1)
        high_selector[direction] = slice(1, None)
        low_index = tuple(low_selector)
        high_index = tuple(high_selector)

        k_left = permeability_values[low_index]
        k_right = permeability_values[high_index]
        k_face = _harmonic_permeability(k_left, k_right)
        face_velocity = k_face * (pressure_values[low_index] - pressure_values[high_index])
        face_velocity /= float(viscosity) * spacing
        high_face[low_index] = face_velocity
        low_face[high_index] = face_velocity

        if direction == flow_axis_index:
            inlet_selector: list[slice | int] = [slice(None)] * pressure_values.ndim
            outlet_selector: list[slice | int] = [slice(None)] * pressure_values.ndim
            inlet_selector[direction] = 0
            outlet_selector[direction] = pressure_values.shape[direction] - 1
            inlet_index = tuple(inlet_selector)
            outlet_index = tuple(outlet_selector)
            low_face[inlet_index] = (
                permeability_values[inlet_index]
                * (float(pressure_inlet) - pressure_values[inlet_index])
                / (float(viscosity) * (spacing / 2.0))
            )
            high_face[outlet_index] = (
                permeability_values[outlet_index]
                * (pressure_values[outlet_index] - float(pressure_outlet))
                / (float(viscosity) * (spacing / 2.0))
            )

        velocity[direction] = 0.5 * (low_face + high_face)

    return velocity

write_structured_vector_field 

write_structured_vector_field(
    vector,
    grid,
    path,
    *,
    name="velocity",
    extra_cell_data=None,
    file_format=None,
)

Write a regular-grid vector field to a ParaView-readable mesh file.

The output uses the structured map mesh representation in voids.mesh and stores the vector as cell data. Use a .vtu or .vtk suffix for common ParaView workflows.

Source code in src/voids/visualization/fields.py 
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def write_structured_vector_field(
    vector: np.ndarray,
    grid: PorosityMap,
    path: str | Path,
    *,
    name: str = "velocity",
    extra_cell_data: Mapping[str, np.ndarray] | None = None,
    file_format: str | None = None,
) -> Path:
    """Write a regular-grid vector field to a ParaView-readable mesh file.

    The output uses the structured map mesh representation in `voids.mesh` and
    stores the vector as cell data. Use a `.vtu` or `.vtk` suffix for common
    ParaView workflows.
    """

    if not name:
        raise ValueError("name must be a non-empty string")
    vec = _as_vector_field(vector, shape=grid.shape)
    if not np.all(np.isfinite(vec)):
        raise ValueError("vector field must contain only finite values")

    mesh = structured_map_mesh(
        grid,
        extra_cell_data=extra_cell_data,
        require_finite_cell_data=True,
    )
    vector_cell_data = np.moveaxis(vec, 0, -1).reshape((-1, vec.shape[0]), order="C")
    if vector_cell_data.shape[1] == 2:
        vector_cell_data = np.column_stack(
            [vector_cell_data, np.zeros(vector_cell_data.shape[0], dtype=vector_cell_data.dtype)]
        )
    mesh.cell_data[str(name)] = [vector_cell_data]
    return mesh.write(path, file_format=file_format)

write_dolfinx_function_xdmf 

write_dolfinx_function_xdmf(function, path, *, name=None)

Write a DOLFINx function to an XDMF/HDF5 file readable by ParaView.

DOLFINx XDMF output expects a function layout compatible with the mesh geometry. To make high-order and discontinuous fields robustly viewable, the function is first interpolated to a first-order Lagrange visualization space.

Source code in src/voids/visualization/fields.py 
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
def write_dolfinx_function_xdmf(
    function: Any, path: str | Path, *, name: str | None = None
) -> Path:
    """Write a DOLFINx function to an XDMF/HDF5 file readable by ParaView.

    DOLFINx XDMF output expects a function layout compatible with the mesh
    geometry. To make high-order and discontinuous fields robustly viewable, the
    function is first interpolated to a first-order Lagrange visualization space.
    """

    from dolfinx.io import XDMFFile

    destination = Path(path)
    destination.parent.mkdir(parents=True, exist_ok=True)
    export_function = _linear_dolfinx_export_function(function, name=name)
    mesh = export_function.function_space.mesh
    with XDMFFile(mesh.comm, str(destination), "w") as handle:
        handle.write_mesh(mesh)
        handle.write_function(export_function)
    return destination

sample_dolfinx_function_on_grid 

sample_dolfinx_function_on_grid(
    function,
    *,
    shape,
    cell_size,
    origin=None,
    fill_value=np.nan,
)

Sample a DOLFINx scalar or vector function at regular cell centers.

Source code in src/voids/visualization/fields.py 
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
def sample_dolfinx_function_on_grid(
    function: Any,
    *,
    shape: Sequence[int],
    cell_size: float | Sequence[float],
    origin: Sequence[float] | None = None,
    fill_value: float = np.nan,
) -> np.ndarray:
    """Sample a DOLFINx scalar or vector function at regular cell centers."""

    from dolfinx import geometry

    grid_shape = tuple(int(value) for value in shape)
    if len(grid_shape) not in {2, 3}:
        raise ValueError("shape must describe a 2D or 3D grid")
    size = _cell_size_tuple(cell_size, len(grid_shape))
    start = (0.0,) * len(grid_shape) if origin is None else tuple(float(value) for value in origin)
    if len(start) != len(grid_shape):
        raise ValueError("origin dimensionality must match shape")

    axes = [
        start[axis] + (np.arange(grid_shape[axis], dtype=float) + 0.5) * size[axis]
        for axis in range(len(grid_shape))
    ]
    meshgrid = np.meshgrid(*axes, indexing="ij")
    points = np.column_stack([component.reshape(-1, order="C") for component in meshgrid])
    if len(grid_shape) == 2:
        points = np.column_stack([points, np.zeros(points.shape[0], dtype=float)])

    mesh = function.function_space.mesh
    tree = geometry.bb_tree(mesh, mesh.topology.dim)
    candidate_cells = geometry.compute_collisions_points(tree, points)
    colliding_cells = geometry.compute_colliding_cells(mesh, candidate_cells, points)

    valid_indices: list[int] = []
    cells: list[int] = []
    for point_index in range(points.shape[0]):
        links = colliding_cells.links(np.int32(point_index))
        if len(links) > 0:
            valid_indices.append(point_index)
            cells.append(int(links[0]))
    if not valid_indices:
        raise RuntimeError("no grid sample points were found inside the DOLFINx mesh")

    valid_points = points[np.asarray(valid_indices, dtype=np.int64)]
    values = np.asarray(function.eval(valid_points, np.asarray(cells, dtype=np.int32)), dtype=float)
    if values.ndim == 1:
        values = values[:, np.newaxis]

    flat = np.full((points.shape[0], values.shape[1]), float(fill_value), dtype=float)
    flat[np.asarray(valid_indices, dtype=np.int64)] = values
    if values.shape[1] == 1:
        return flat[:, 0].reshape(grid_shape, order="C")
    return np.moveaxis(flat.reshape((*grid_shape, values.shape[1]), order="C"), -1, 0)

plot_scalar_midplanes 

plot_scalar_midplanes(
    scalar,
    *,
    title,
    path=None,
    cmap="viridis",
    colorbar_label=None,
    vmin=None,
    vmax=None,
    colorbar_use_offset=True,
)

Plot scalar-field mid-slices and optionally save the figure.

Source code in src/voids/visualization/fields.py 
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
def plot_scalar_midplanes(
    scalar: np.ndarray,
    *,
    title: str,
    path: str | Path | None = None,
    cmap: str = "viridis",
    colorbar_label: str | None = None,
    vmin: float | None = None,
    vmax: float | None = None,
    colorbar_use_offset: bool = True,
) -> Any:
    """Plot scalar-field mid-slices and optionally save the figure."""

    import matplotlib.pyplot as plt

    values = _as_scalar_field(scalar)
    image_vmin = float(np.min(values)) if vmin is None else vmin
    image_vmax = float(np.max(values)) if vmax is None else vmax
    specs = _midplane_specs(values.shape)
    fig, axes = plt.subplots(
        1, len(specs), figsize=(4.2 * len(specs), 3.8), constrained_layout=True
    )
    axes_array = np.atleast_1d(axes)
    for ax, spec in zip(axes_array, specs, strict=True):
        image = _scalar_slice(values, spec.plane_axis).T
        im = ax.imshow(image, origin="lower", cmap=cmap, vmin=image_vmin, vmax=image_vmax)
        ax.set_title(spec.label)
        ax.set_xlabel(_AXIS_NAMES[spec.in_plane_axes[0]])
        ax.set_ylabel(_AXIS_NAMES[spec.in_plane_axes[1]])
        colorbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label=colorbar_label)
        if not colorbar_use_offset:
            from matplotlib.ticker import ScalarFormatter

            formatter = ScalarFormatter(useOffset=False)
            formatter.set_scientific(False)
            colorbar.formatter = formatter
            colorbar.update_ticks()
    fig.suptitle(title)
    if path is not None:
        destination = Path(path)
        destination.parent.mkdir(parents=True, exist_ok=True)
        fig.savefig(destination, dpi=180)
    return fig

plot_vector_midplanes 

plot_vector_midplanes(
    vector,
    *,
    title,
    path=None,
    quiver_stride=3,
    cmap="magma",
    colorbar_label="velocity magnitude",
    vmin=None,
    vmax=None,
)

Plot vector-field magnitude mid-slices with in-plane quiver arrows.

Source code in src/voids/visualization/fields.py 
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
def plot_vector_midplanes(
    vector: np.ndarray,
    *,
    title: str,
    path: str | Path | None = None,
    quiver_stride: int = 3,
    cmap: str = "magma",
    colorbar_label: str = "velocity magnitude",
    vmin: float | None = None,
    vmax: float | None = None,
) -> Any:
    """Plot vector-field magnitude mid-slices with in-plane quiver arrows."""

    import matplotlib.pyplot as plt

    if quiver_stride <= 0:
        raise ValueError("quiver_stride must be positive")
    vec = _as_vector_field(vector)
    magnitude = vector_magnitude(vec)
    finite_magnitude = magnitude[np.isfinite(magnitude)]
    if not finite_magnitude.size:
        raise ValueError("vector field must contain at least one finite magnitude")
    image_vmin = float(np.min(finite_magnitude)) if vmin is None else vmin
    image_vmax = float(np.max(finite_magnitude)) if vmax is None else vmax
    specs = _midplane_specs(magnitude.shape)
    fig, axes = plt.subplots(
        1, len(specs), figsize=(4.5 * len(specs), 3.9), constrained_layout=True
    )
    axes_array = np.atleast_1d(axes)
    for ax, spec in zip(axes_array, specs, strict=True):
        magnitude_slice = _scalar_slice(magnitude, spec.plane_axis)
        im = ax.imshow(
            magnitude_slice.T, origin="lower", cmap=cmap, vmin=image_vmin, vmax=image_vmax
        )
        comp_x = _scalar_slice(vec[spec.in_plane_axes[0]], spec.plane_axis)
        comp_y = _scalar_slice(vec[spec.in_plane_axes[1]], spec.plane_axis)
        x = np.arange(magnitude_slice.shape[0])
        y = np.arange(magnitude_slice.shape[1])
        grid_x, grid_y = np.meshgrid(x, y)
        stride = (slice(None, None, quiver_stride), slice(None, None, quiver_stride))
        u = comp_x.T[stride]
        v = comp_y.T[stride]
        in_plane_magnitude = np.hypot(u, v)
        finite_magnitude = in_plane_magnitude[np.isfinite(in_plane_magnitude)]
        max_in_plane_magnitude = float(np.max(finite_magnitude)) if finite_magnitude.size else 0.0
        if max_in_plane_magnitude > 0.0:
            max_arrow_length = max(1.0, 0.65 * float(quiver_stride))
            ax.quiver(
                grid_x[stride],
                grid_y[stride],
                u,
                v,
                color="white",
                angles="xy",
                scale_units="xy",
                scale=max_in_plane_magnitude / max_arrow_length,
                width=0.003,
            )
        ax.set_title(spec.label)
        ax.set_xlabel(_AXIS_NAMES[spec.in_plane_axes[0]])
        ax.set_ylabel(_AXIS_NAMES[spec.in_plane_axes[1]])
        fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label=colorbar_label)
    fig.suptitle(title)
    if path is not None:
        destination = Path(path)
        destination.parent.mkdir(parents=True, exist_ok=True)
        fig.savefig(destination, dpi=180)
    return fig