Utilities

Utility functions and infrastructure for Underworld3.

This module provides supporting utilities including I/O helpers, geometry tools, array callbacks, and mathematical mixins.

I/O Utilities

Xdmf, generateXdmf, generate_uw_Xdmf

XDMF file generation for ParaView visualization.

swarm_h5, swarm_xdmf

Swarm data I/O in HDF5 and XDMF formats.

read_medit_ascii, create_dmplex_from_medit

Import meshes from Medit format.

Geometry Tools

Various geometric helper functions for mesh operations.

Development Utilities

CaptureStdout

Context manager for capturing stdout.

h5_scan

HDF5 file inspection.

mem_footprint

Memory usage tracking.

auditor, postHog

Analytics and debugging tools.

Notes

Units are handled by the enhanced variable system. See underworld3.discretisation.enhanced_variables for details.

See also

underworld3.function

Expression evaluation with unit support.

underworld3.discretisation

Mesh I/O functions.

I/O Utilities

XDMF Generation

underworld3.utilities.generateXdmf(hdfFilename, xdmfFilename=None)

underworld3.utilities.generate_uw_Xdmf(hdfFilename, xdmfFilename=None)

Swarm I/O

underworld3.utilities.swarm_h5(swarm, fileName, timestep, fields=None, outputPath='', compression=False, compressionType='gzip')

Function to save swarm fields to h5 file for checkpointing purposes.

swarm : The swarm that contains the coordinate data fileName : Name of file the file where the swarm data is stored. fields : UW swarm variables that contain the data to save. Should be passed as a list. timestep : timestep of the model to save outputPath : Folder to save the data, can be left undefined to save in the current working directory compression : Whether to compress the h5 files [bool] compressionType : The type of compression to use. ‘gzip’ and ‘lzf’ are the supported types, with ‘gzip’ as the default.

>>>
example usage:
without compression:
swarm_h5(swarm=swarm, fileName='swarm', fields=[material, strainRate], timestep=0)

with compression: swarm_h5(swarm=swarm, fileName=’swarm’, fields=[material, strainRate], timestep=0, compression=True, compressionType=’gzip’)

---

underworld3.utilities.swarm_xdmf(fileName, timestep, fields=None, outputPath='', time=None)

Function to combine h5 files to a single xdmf to visualise the swarm in paraview. Should be run after the ‘swarm_h5’ function.

fileName : Name of file the file where the swarm data is stored. timestep : timestep of the model to save fields : UW swarm variables that contain the data to save. Should be passed as a list. outputPath : Folder to save the data, can be left undefined to save in the current working directory. This directory should also contain the h5 files time : Time of the model, can be left undefined and model output will be sequential but not contain the time data.

>>>
example usage:

swarm_xdmf(fileName=’swarm’, fields=[material, strainRate], timestep=0)

---

Mesh Import

underworld3.utilities.read_medit_ascii(filename, mesh_data_name)

Function to read a medit mesh file in ascii format. Return mesh data and their indices if found.

underworld3.utilities.create_dmplex_from_medit(medit_file, print_info=False)

Reads a medit (.mesh) file and returns a distributed DMPlex with labels.

Development Utilities

class underworld3.utilities.CaptureStdout

Bases: UserString, redirect_stdout

Captures stdout (e.g., from print()) as a variable.

Based on contextlib.redirect_stdout, but saves the user the trouble of defining and reading from an IO stream. Useful for testing the output of functions that are supposed to print certain output.

Citation: https://stackoverflow.com/questions/16571150/how-to-capture-stdout-output-from-a-python-function-call

__init__(seq='', split=False, *args, **kwargs)

start()

stop()

underworld3.utilities.h5_scan(filename)

underworld3.utilities.mem_footprint()

Returns resident set size in Mb for this process

Geometry Utilities

Functions for geometric computations useful for mesh and swarm operations.

Distance Calculations

underworld3.utilities.distance_pointcloud_linesegment(p, a, b)

p - numpy array of points a, b - line-segment points (numpy 1xdim arrays)

returns:

numpy array of distances from each of the points to the nearest point within the triangle (0 if in the plane, within the triangle)

underworld3.utilities.distance_pointcloud_polyline(p, vertices)

Compute distance from points to a polyline (multiple connected segments).

Parameters:

• p (numpy array (n, dim)) – Query points

• vertices (numpy array (m, dim)) – Ordered vertices defining the polyline (m-1 segments)

Returns:

Distance from each point to the nearest point on the polyline

Return type:

numpy array (n,)

underworld3.utilities.distance_pointcloud_triangle(p, a, b, c)

p - numpy array of points in 3D a, b, c - triangle points (numpy 1x3 arrays)

returns:

numpy array of distances from each of the points to the nearest point within the triangle (0 if in the plane, within the triangle)

Signed Distance Functions

underworld3.utilities.signed_distance_pointcloud_linesegment_2d(p, a, b)

Compute signed distance from 2D points to a line segment.

The sign is determined by which side of the line the point is on: - Positive: to the “left” of the segment (when looking from a to b) - Negative: to the “right” of the segment

Parameters:

• p (numpy array (n, 2)) – Query points in 2D

• a (numpy array (2,) or (1, 2)) – Start point of line segment

• b (numpy array (2,) or (1, 2)) – End point of line segment

Returns:

Signed distances from each point to the line segment

Return type:

numpy array (n,)

underworld3.utilities.signed_distance_pointcloud_polyline_2d(p, vertices)

Compute signed distance from 2D points to a polyline.

The sign is determined by the nearest segment: - Positive: to the “left” of the polyline (consistent orientation) - Negative: to the “right” of the polyline

Parameters:

• p (numpy array (n, 2)) – Query points in 2D

• vertices (numpy array (m, 2)) – Ordered vertices defining the polyline (m-1 segments)

Returns:

Signed distance from each point to the polyline

Return type:

numpy array (n,)

Point-in-Simplex Tests

underworld3.utilities.points_in_simplex2D(p, a, b, c)

p - numpy array of points in 3D a, b, c - triangle points (numpy 1x2 arrays)

returns:

numpy array of truth values for each of the points in p (is this point in the triangle)

underworld3.utilities.points_in_simplex3D(p, a, b, c, d)

Normals

underworld3.utilities.linesegment_normals_2d(points)

Compute unit normals for line segments defined by ordered points (2D).

For a segment from points[i] to points[i+1], the normal points to the “left” side (90° counterclockwise rotation of the tangent).

Parameters:

points (numpy array (n, 2)) – Ordered vertices defining n-1 line segments

Returns:

• segment_normals (numpy array (n-1, 2)) – Unit normal for each segment

• vertex_normals (numpy array (n, 2)) – Unit normal at each vertex (averaged from adjacent segments)

Array Utilities

UnitAwareArray

Array with attached unit information for dimensional quantities.

class underworld3.utilities.UnitAwareArray

Bases: NDArray_With_Callback

A numpy ndarray subclass that combines callback functionality with unit awareness.

Extends NDArray_With_Callback to provide automatic unit tracking, compatibility checking, and integration with the UW3 unit conversion system.

Operations preserve dimensional consistency: compatible units are added, incompatible units raise errors, and multiplication combines units.

Global Reduction Operations (MPI-aware, unit-preserving):

• global_max() -> UWQuantity (same units as array)

• global_min() -> UWQuantity (same units as array)

• global_sum() -> UWQuantity (same units as array)

• global_mean() -> UWQuantity (same units as array)

• global_size() -> int (count, no units)

• global_norm() -> UWQuantity (same units as array)

• global_rms() -> UWQuantity (same units as array)

• global_var() -> UWQuantity (units squared)

• global_std() -> UWQuantity (same units as array)

Tensor arrays (ndim > 2) raise NotImplementedError for global reductions. Use component-wise operations or slice the array for tensors.

Inherits the callback mechanism from NDArray_With_Callback, enabling it to work with any storage backend (PETSc, pyvista, etc.) by changing only the callback. Provides consistent .units, .magnitude, and unit-aware arithmetic regardless of underlying storage.

static __new__(cls, input_array=None, units=None, owner=None, callback=None, unit_checking=True, auto_convert=True, **kwargs)

Create new UnitAwareArray instance.

Parameters:

• input_array (array-like, optional) – Input data to create array from

• units (str or UWQuantity, optional) – Units for this array (e.g., “m”, “m/s”, “kg”)

• owner (object, optional) – The object that owns this array (stored as weak reference)

• callback (callable, optional) – Initial callback function to register

• unit_checking (bool, optional) – If True, enforce unit compatibility in operations (default True)

• auto_convert (bool, optional) – If True, automatically convert compatible units (default True)

• **kwargs (dict) – Additional arguments passed to NDArray_With_Callback

property units

Get the units of this array.

property has_units

Check if this array has units.

property dimensionality

Get the dimensionality of this array.

property unit_checking

Check if unit compatibility checking is enabled.

property auto_convert

Check if automatic unit conversion is enabled.

property magnitude

Get the numerical values without units (like Pint’s .magnitude).

This returns a plain numpy array view of the data, stripping units. Useful when you need raw numerical values for dimensionless calculations.

Returns:

Plain numpy array without unit tracking

Return type:

np.ndarray

Examples

>>> coords = mesh.X.coords  # UnitAwareArray with units="km"
>>> x, y = coords[:, 0].magnitude, coords[:, 1].magnitude  # Plain arrays
>>> temperature.array[:, 0, 0] = 300 + 2.6 * y  # Works - no units

to(target_units)

Convert this array to different units.

Provides a unified interface matching Pint’s .to() pattern.

Parameters:

target_units (str) – Target units to convert to

Returns:

New array with converted values and target units

Return type:

UnitAwareArray

Examples

>>> coords = UnitAwareArray([1, 2, 3], units='km')
>>> coords_m = coords.to('m')  # Convert to meters
>>> print(coords_m)
[1000. 2000. 3000.] [meter]

max(axis=None, out=None, keepdims=False, initial=None, where=True)

Return maximum with units preserved.

min(axis=None, out=None, keepdims=False, initial=None, where=True)

Return minimum with units preserved.

mean(axis=None, dtype=None, out=None, keepdims=False, where=True)

Return mean with units preserved.

sum(axis=None, dtype=None, out=None, keepdims=False, initial=None, where=True)

Return sum with units preserved.

std(axis=None, dtype=None, out=None, ddof=0, keepdims=False, where=True)

Return standard deviation with units preserved.

var(axis=None, dtype=None, out=None, ddof=0, keepdims=False, where=True)

Return variance with units squared.

global_max(axis=None, out=None, keepdims=False)

Return maximum across all MPI ranks with units preserved.

For scalar results (axis=None), performs MPI reduction. For array results, performs component-wise maximum. For tensors (ndim > 2), raises NotImplementedError.

Parameters:

• axis (None or int or tuple of ints, optional) – Axis along which to operate (default: None = reduce all dimensions)

• out (ndarray, optional) – Alternative output array

• keepdims (bool, optional) – Keep reduced dimensions as size 1 (default: False)

Returns:

Global maximum with units preserved

Return type:

UWQuantity or ndarray

Raises:

NotImplementedError – If called on tensor data (ndim > 2)

global_min(axis=None, out=None, keepdims=False)

Return minimum across all MPI ranks with units preserved.

For scalar results (axis=None), performs MPI reduction. For array results, performs component-wise minimum. For tensors (ndim > 2), raises NotImplementedError.

Parameters:

• axis (None or int or tuple of ints, optional) – Axis along which to operate (default: None = reduce all dimensions)

• out (ndarray, optional) – Alternative output array

• keepdims (bool, optional) – Keep reduced dimensions as size 1 (default: False)

Returns:

Global minimum with units preserved

Return type:

UWQuantity or ndarray

Raises:

NotImplementedError – If called on tensor data (ndim > 2)

global_sum(axis=None, dtype=None, out=None, keepdims=False)

Return sum across all MPI ranks with units preserved.

For scalar results (axis=None), performs MPI reduction. For array results, performs component-wise sum. For tensors (ndim > 2), raises NotImplementedError.

Parameters:

• axis (None or int or tuple of ints, optional) – Axis along which to operate (default: None = reduce all dimensions)

• dtype (data-type, optional) – Type of returned array

• out (ndarray, optional) – Alternative output array

• keepdims (bool, optional) – Keep reduced dimensions as size 1 (default: False)

Returns:

Global sum with units preserved

Return type:

UWQuantity or ndarray

Raises:

NotImplementedError – If called on tensor data (ndim > 2)

global_mean(axis=None, dtype=None, out=None, keepdims=False)

Return mean across all MPI ranks with units preserved.

Computes the true global mean by summing all values across ranks and dividing by total count. For tensors (ndim > 2), raises NotImplementedError.

Parameters:

• axis (None or int or tuple of ints, optional) – Axis along which to operate (default: None = reduce all dimensions)

• dtype (data-type, optional) – Type of returned array

• out (ndarray, optional) – Alternative output array

• keepdims (bool, optional) – Keep reduced dimensions as size 1 (default: False)

Returns:

Global mean with units preserved

Return type:

UWQuantity or ndarray

Raises:

NotImplementedError – If called on tensor data (ndim > 2)

global_var(axis=None, dtype=None, ddof=0, keepdims=False)

Return variance across all MPI ranks with units squared preserved.

Uses parallel variance algorithm (Welford/Chan) for numerical stability. For tensors (ndim > 2), raises NotImplementedError.

Parameters:

• axis (None or int or tuple of ints, optional) – Axis along which to operate (default: None = reduce all dimensions)

• dtype (data-type, optional) – Type of returned array

• ddof (int, optional) – Delta degrees of freedom (default: 0)

• keepdims (bool, optional) – Keep reduced dimensions as size 1 (default: False)

Returns:

Global variance with units squared

Return type:

UWQuantity or ndarray

Raises:

NotImplementedError – If called on tensor data (ndim > 2)

global_std(axis=None, dtype=None, ddof=0, keepdims=False)

Return standard deviation across all MPI ranks with units preserved.

Computed as square root of global variance. For tensors (ndim > 2), raises NotImplementedError.

Parameters:

• axis (None or int or tuple of ints, optional) – Axis along which to operate (default: None = reduce all dimensions)

• dtype (data-type, optional) – Type of returned array

• ddof (int, optional) – Delta degrees of freedom (default: 0)

• keepdims (bool, optional) – Keep reduced dimensions as size 1 (default: False)

Returns:

Global standard deviation with units preserved

Return type:

UWQuantity or ndarray

Raises:

NotImplementedError – If called on tensor data (ndim > 2)

global_norm(ord=None)

Return norm across all MPI ranks.

For scalars (ndim=1), computes sqrt(sum of squares). For vectors, computes vector norm. For tensors (ndim > 2), raises NotImplementedError.

Parameters:

ord ({non-zero int, inf, -inf, 'fro', 'nuc'}, optional) – Order of the norm (default: None = 2-norm)

Returns:

Global norm with units preserved

Return type:

UWQuantity or float

Raises:

NotImplementedError – If called on tensor data (ndim > 2)

global_size()

Return total number of elements across all MPI ranks.

Useful for computing global statistics that require total element count, such as RMS or normalized quantities.

Returns:

Total number of elements summed across all MPI ranks

Return type:

int

Examples

>>> coords = mesh.X.coords  # Shape: (N_local, 2)
>>> total_points = coords.global_size()  # Sum of N_local across all ranks
>>> rms = coords.global_norm() / np.sqrt(total_points)

global_rms()

Return root mean square across all MPI ranks with units preserved.

Computes RMS = sqrt(sum of squares / total count) across all ranks. For tensors (ndim > 2), raises NotImplementedError.

The RMS is computed as: RMS = global_norm() / sqrt(global_size())

Returns:

Global RMS with units preserved

Return type:

UWQuantity or float

Raises:

NotImplementedError – If called on tensor data (ndim > 2)

Examples

>>> coords = mesh.X.coords  # UnitAwareArray with units="km"
>>> rms_coord = coords.global_rms()  # Returns UWQuantity in km
>>> print(f"RMS coordinate: {rms_coord}")

copy(order='C')

Return a copy of the array with preserved units.

astype(dtype, order='K', casting='unsafe', subok=True, copy=True)

Convert array type while preserving units.

view(dtype=None, type=None)

Return a view of the array with preserved units.

reshape(*shape, order='C')

Return a reshaped array with preserved units.

flatten(order='C')

Return a flattened array with preserved units.

squeeze(axis=None)

Return a squeezed array with preserved units.

transpose(*axes)

Return a transposed array with preserved units.

NDArray_With_Callback

Array wrapper that triggers callbacks on modification.

class underworld3.utilities.NDArray_With_Callback

Bases: ndarray

A numpy ndarray subclass that triggers callbacks when array data is modified.

This class maintains full numpy array compatibility while providing reactive programming capabilities for scientific computing applications.

Callback Function Signature:

def callback(array: NDArray_With_Callback, change_info: dict) -> None:
    pass

The change_info dictionary contains:

• operation (str): Operation name (‘setitem’, ‘iadd’, ‘fill’, etc.)

• indices (tuple/slice/None): Location of change (for setitem operations)

• old_value (array-like/None): Previous values (when available)

• new_value (array-like): New values being assigned

• array_shape (tuple): Current shape of the array

• array_dtype (np.dtype): Data type of the array

Features:

• Multiple callbacks: add_callback(), remove_callback(), clear_callbacks()

• Enable/disable: enable_callbacks(), disable_callbacks()

• Delayed execution: delay_callback(), delay_callbacks_global()

• MPI synchronization: Automatic barriers in parallel contexts

• Weak references: Owner tracking without circular dependencies

• Global reductions: MPI-aware global_max(), global_min(), global_sum(), etc.

Global Reduction Operations (MPI-aware):

• global_max(axis=None): Maximum value across all MPI ranks

• global_min(axis=None): Minimum value across all MPI ranks

• global_sum(axis=None): Sum of all values across all MPI ranks

• global_mean(axis=None): True mean (global sum / global count)

• global_size(): Total number of elements across all ranks

• global_norm(ord=2): 2-norm (Euclidean) across all ranks

• global_rms(): Root mean square across all ranks

These methods use MPI collective operations (allreduce). All ranks must call these methods (they are collective operations). Subclasses like UnitAwareArray override these to preserve units.

static __new__(cls, input_array=None, owner=None, callback=None, disable_inplace_operators=False)

Create new NDArray_With_Callback instance.

Parameters:

• input_array (array-like, optional) – Input data to create array from (defaults to empty array if None)

• owner (object, optional) – The object that owns this array (stored as weak reference)

• callback (callable, optional) – Initial callback function to register

• disable_inplace_operators (bool, optional) – If True, in-place operators (+=, -=, *=, /=, etc.) will raise RuntimeError for parallel safety. Default is False for backward compatibility.

property mask

Public mask property for numpy.ma compatibility.

Matplotlib’s quiver and other plotting functions access .mask directly. This aliases to _mask which returns np.ma.nomask (no masking).

filled(fill_value=None)

Return array with masked values filled.

For numpy.ma compatibility. Since we have no mask, this just returns a copy of the data (as numpy array to avoid further masked array operations).

Parameters:

fill_value (scalar, optional) – Value used to fill masked entries. Ignored since we have no mask.

Returns:

A copy of the data as a plain numpy array.

Return type:

ndarray

set_callback(callback)

Set a single callback function (replaces any existing callbacks).

Parameters:

callback (callable) – Function with signature: callback(array, change_info) - array: the NDArray_With_Callback instance - change_info: dict with operation details

add_callback(callback)

Add an additional callback function.

Parameters:

callback (callable) – Function to add to callback list

remove_callback(callback)

Remove a specific callback function.

Parameters:

callback (callable) – Function to remove from callback list

clear_callbacks()

Remove all registered callbacks.

enable_callbacks()

Enable callback triggering.

disable_callbacks()

Disable callback triggering (useful for batch operations).

property owner

Get the owner object (may be None if owner was garbage collected).

delay_callback(context_info=None)

Context manager to delay callback execution until context exit.

During the context, all callbacks from this array (and any other arrays using delay_callback) will be accumulated and executed when the outermost context exits.

Parameters:

context_info (str, optional) – Optional information about the context (for debugging)

Example

>>> with arr.delay_callback("batch update"):
...     arr[0] = 1
...     arr[1] = 2
...     arr[2] = 3
# All callbacks fire here at context exit

static delay_callbacks_global(context_info=None)

Static method to create a global delay context for all NDArray_With_Callback instances.

This is useful when you don’t have a specific array instance but want to delay callbacks from multiple arrays.

Example

>>> with NDArray_With_Callback.delay_callbacks_global("mesh update"):
...     mesh.data[0] = new_pos
...     swarm.data += displacement
# All callbacks from all arrays fire here

fill(value)

Fill array with scalar value, triggering callback.

sort(axis=-1, kind=None, order=None)

Sort array in-place, triggering callback.

resize(new_shape, refcheck=True)

Resize array in-place, triggering callback.

copy(order='C')

Return a copy of the array.

The copy will have the same callbacks registered but will be independent.

view(dtype=None, type=None)

Return a view of the array.

Views share callbacks with the original array.

sync_data(new_data)

Update array with new data, preserving callbacks and all metadata.

This method efficiently handles both same-size and different-size data updates. For same-size updates, it uses efficient in-place copying. For different sizes, it creates a new array object but preserves all metadata and callbacks.

Parameters:

new_data (array-like) – New data to sync into this array. Can be different size/shape.

Returns:

result – For same-size: returns self (same object) For different-size: returns new object with same metadata

Return type:

NDArray_With_Callback

Notes

• For same-size data: Uses efficient in-place copy (preserves object identity)

• For different sizes: Creates new object but copies all callbacks/metadata

• All callbacks, owner references, and settings are preserved

• Triggers ‘sync_data’ callback after update

Examples

>>> arr = NDArray_With_Callback([1, 2, 3])
>>> result = arr.sync_data([4, 5, 6])  # Same size: returns same object
>>> assert result is arr
>>> result = arr.sync_data([7, 8, 9, 10, 11])  # Different size: new object
>>> assert result is not arr  # Different object
>>> assert len(result._callbacks) == len(arr._callbacks)  # Same callbacks

global_max(axis=None, out=None, keepdims=False)

Return maximum across all MPI ranks.

For scalar results (axis=None), performs MPI reduction. For array results, performs component-wise maximum.

Parameters:

• axis (None or int or tuple of ints, optional) – Axis along which to operate (default: None = reduce all dimensions)

• out (ndarray, optional) – Alternative output array

• keepdims (bool, optional) – Keep reduced dimensions as size 1 (default: False)

Returns:

Global maximum value(s)

Return type:

scalar or ndarray

global_min(axis=None, out=None, keepdims=False)

Return minimum across all MPI ranks.

For scalar results (axis=None), performs MPI reduction. For array results, performs component-wise minimum.

Parameters:

• axis (None or int or tuple of ints, optional) – Axis along which to operate (default: None = reduce all dimensions)

• out (ndarray, optional) – Alternative output array

• keepdims (bool, optional) – Keep reduced dimensions as size 1 (default: False)

Returns:

Global minimum value(s)

Return type:

scalar or ndarray

global_sum(axis=None, dtype=None, out=None, keepdims=False)

Return sum across all MPI ranks.

For scalar results (axis=None), performs MPI reduction. For array results, performs component-wise sum.

Parameters:

• axis (None or int or tuple of ints, optional) – Axis along which to operate (default: None = reduce all dimensions)

• dtype (data-type, optional) – Type of returned array

• out (ndarray, optional) – Alternative output array

• keepdims (bool, optional) – Keep reduced dimensions as size 1 (default: False)

Returns:

Global sum value(s)

Return type:

scalar or ndarray

global_mean(axis=None, dtype=None, out=None, keepdims=False)

Return mean across all MPI ranks.

Computes the true global mean by summing all values across ranks and dividing by total count.

Parameters:

• axis (None or int or tuple of ints, optional) – Axis along which to operate (default: None = reduce all dimensions)

• dtype (data-type, optional) – Type of returned array

• out (ndarray, optional) – Alternative output array

• keepdims (bool, optional) – Keep reduced dimensions as size 1 (default: False)

Returns:

Global mean value(s)

Return type:

scalar or ndarray

global_size()

Return total number of elements across all MPI ranks.

Useful for computing global statistics that require total element count.

Returns:

Total number of elements summed across all MPI ranks

Return type:

int

global_norm(ord=None)

Return 2-norm across all MPI ranks.

Computes sqrt(sum of squares) across all ranks.

Parameters:

ord ({None, 2}, optional) – Order of the norm (only 2-norm supported, default: None = 2-norm)

Returns:

Global 2-norm value

Return type:

float

global_rms()

Return root mean square across all MPI ranks.

Computes RMS = sqrt(sum of squares / total count) across all ranks.

Returns:

Global RMS value

Return type:

float