Discretisation¶

Mesh¶

The computational mesh class that provides the spatial discretisation for finite element calculations.

class underworld3.discretisation.Mesh[source]¶

Bases: Stateful, uw_object

Unstructured mesh with PETSc DMPlex backend.

The Mesh class provides the spatial discretisation for finite element computations. It wraps PETSc’s DMPlex for unstructured mesh management, supporting various cell types (triangles, quadrilaterals, tetrahedra, hexahedra) and coordinate systems.

Parameters:

plex_or_meshfile (PETSc.DMPlex or str) – Either a PETSc DMPlex object or path to a mesh file (gmsh, exodus).
degree (int, optional) – Polynomial degree for the coordinate field (default 1).
simplex (bool, optional) – True for simplicial elements (triangles/tets), False for quads/hexes.
coordinate_system_type (CoordinateSystemType, optional) – Coordinate system for vector calculus (Cartesian, cylindrical, etc.).
qdegree (int, optional) – Quadrature degree for numerical integration (default 2).
boundaries (list of NamedTuple, optional) – Boundary region definitions with names and values.
boundary_normals (dict, optional) – Outward normal vectors for each boundary.
units (str or pint.Unit, optional) – Physical units for mesh coordinates.
verbose (bool, optional) – Print mesh construction information.

Examples

Meshes are typically created via the meshing module:

>>> mesh = uw.meshing.UnstructuredSimplexBox(
...     minCoords=(0, 0), maxCoords=(1, 1), cellSize=0.1
... )
>>> T = mesh.add_variable("T", vtype=uw.VarType.SCALAR)

See also

underworld3.meshing: Mesh generation utilities.
underworld3.discretisation.MeshVariable: Field variables on meshes.

mesh_instances = 0¶

__init__(plex_or_meshfile, degree=1, simplex=True, coordinate_system_type=None, qdegree=2, markVertices=None, useRegions=None, useMultipleTags=None, filename=None, refinement=None, refinement_callback=None, coarsening=None, coarsening_callback=None, return_coords_to_bounds=None, boundaries=None, boundary_normals=None, name=None, units=None, verbose=False, *args, **kwargs)[source]¶

property dim: int¶

Topological dimension of the mesh.

Returns:: The mesh dimension (2 for 2D, 3 for 3D).
Return type:: int

property cdim: int¶

Coordinate dimension (embedding space dimension).

For most meshes, cdim == dim. For surface meshes embedded in 3D (e.g., a 2D spherical shell), dim=2 but cdim=3.

Returns:: The coordinate dimension.
Return type:: int

property element: dict¶

Element type information for the mesh.

Contains details about the finite element discretization including cell type, polynomial degree, and quadrature order.

Returns:: Element information dictionary.
Return type:: dict

Notes

UW3 does not support mixed-element meshes; this applies uniformly to all cells.

property length_scale: float¶

Length scale for non-dimensionalization.

This property is IMMUTABLE after mesh creation to ensure synchronization with all spatial operators (gradient, divergence, curl, etc.).

The length scale is derived from model reference quantities at mesh creation: - Priority 1: domain_depth from model.set_reference_quantities() - Priority 2: length from model.set_reference_quantities() - Default: 1.0 (no scaling)

Returns:: Length scale value for non-dimensionalization
Return type:: float

Examples

>>> model.set_reference_quantities(domain_depth=uw.quantity(100, "km"))
>>> mesh = uw.meshing.UnstructuredSimplexBox(...)
>>> mesh.length_scale
100000.0  # meters

See also

length_units: Units string for length scale

property length_units: str¶

Unit string for the length scale.

Returns:: Units for the length scale (e.g., “meter”, “kilometer”)
Return type:: str

Examples

>>> mesh.length_units
'kilometer'

quality(per_cell=False)[source]¶

Cell-quality diagnostics relevant to FE / solver conditioning.

Bulk volume ratios (min/mean) hide the handful of near-degenerate cells that nonetheless dominate stiffness-matrix conditioning — a Stokes / saddle-point solve line-search-fails on the worst element, not the mean. This reports the tail metrics that actually predict that. For a 2-D simplex (triangle) mesh, per cell:

shape quality q = 4√3·A / Σℓ² (1 = equilateral, → 0 = sliver; folds skew + stretch into one number)
largest interior angle (→ 180° is the conditioning killer)
aspect ratio ℓ_max² / (2A) (longest edge / shortest altitude)
neighbour size-jump (adjacent-cell area ratio — the mesh gradation the solver actually sees)

The conditioning-relevant numbers are the worst cell (q_min, angle_max_deg, aspect_max) and the poor-cell counts, not the means. Non-2-D-simplex meshes get the dimension-agnostic cell-volume-spread subset only.

Parameters:: per_cell (bool, default False) – Also return per-cell arrays (q, angle_deg, aspect, volume) under "per_cell" — for plotting or locating the bad cells.
Returns:: Aggregate + tail stats. Headline scalars (min/max/counts) are MPI-reduced so they are correct in parallel; percentiles and the neighbour size-jump are rank-local estimates (exact in serial — the convention for the mesh-redistribution tooling).
Return type:: dict

Examples

>>> q = mesh.quality()
>>> q["q_min"], q["n_q_lt_0p3"], q["aspect_max"]
>>> mesh.quality(per_cell=True)["per_cell"]["q"]  # to plot

view(level=0)[source]¶

Displays mesh information at different levels.

Parameters:: level (int (0 default)) – The display level. 0, for basic mesh information (variables and boundaries), while level=1 displays detailed mesh information (including PETSc information)

view_parallel()[source]¶: returns the break down of boundary labels from each processor

clone_dm_hierarchy()[source]¶: Clone the dm hierarchy on the mesh

extract_region(label_name, label_value=None)[source]¶

Extract a submesh containing only cells with the given region label.

Uses DMPlexFilter to create a new mesh sharing exact node positions with the parent. The submesh carries a subpoint_is mapping back to the parent for restrict/prolongate operations, and a parent reference.

Boundary labels from the parent survive the filter. For example, an “Internal” boundary on the parent becomes an exterior boundary on the submesh and can be referenced by the same name.

Parameters:

label_name (str) – DM label name identifying the region (e.g., "Inner").
label_value (int, optional) – Stratum value within the label. If None, uses mesh.regions.<label_name>.value when available.

Returns:

A new mesh covering only the specified region.

Return type:

Mesh

Examples

>>> full_mesh = uw.meshing.AnnulusInternalBoundary(...)
>>> rock_mesh = full_mesh.extract_region("Inner")
>>> rock_mesh.parent is full_mesh
True

sync_coordinates_from_parent()[source]¶

Update submesh coordinates from the parent mesh.

Called automatically when the parent mesh deforms. Uses the cached vertex map to copy parent vertex positions to the submesh, then calls _deform_mesh to rebuild geometry.

Raises:: ValueError – If this mesh has no parent.

restrict(parent_var, sub_var, mode='replace')[source]¶

Copy data from a parent-mesh variable to a submesh variable.

Parameters:

parent_var (MeshVariable) – Source variable on the parent mesh.
sub_var (MeshVariable) – Destination variable on this (sub)mesh.
mode (str) – "replace" overwrites submesh values (INSERT_VALUES). "add" adds parent values into submesh (ADD_VALUES).

Raises:

ValueError – If this mesh has no parent, or the variable meshes don’t match.

prolongate(sub_var, parent_var, mode='replace')[source]¶

Copy data from a submesh variable to a parent-mesh variable.

Parameters:

sub_var (MeshVariable) – Source variable on this (sub)mesh.
parent_var (MeshVariable) – Destination variable on the parent mesh.
mode (str) – "replace" overwrites parent values at submesh DOFs. "add" adds submesh values into parent.

Raises:

ValueError – If this mesh has no parent, or the variable meshes don’t match.

nuke_coords_and_rebuild(verbose=False, active_vars=None)[source]¶

Rebuild DM/DS, the kd-tree, mesh sizes, and per-variable DOF coordinate caches after a coordinate change.

active_vars (optional set/list of MeshVariables): restrict the per-variable DOF coordinate-cache recomputation to this set. When None (default) every registered variable is recomputed eagerly, matching the BUGFIX(#130) collective-safe behaviour. Movers that thread their own work-vars through _deform_mesh(..., active_vars=...) skip recomputing the non-mover variables n_outer× during the inner sweep; the wrapper does one final _deform_mesh (or a direct nuke_coords_and_rebuild()) with active_vars=None at sweep exit to bring the full cache back into sync.

Naming note: “nuke and rebuild” historically referred to the DS+DM tear-down/recreate; what was called “refill” of the per-variable cache (line 1890 in the old code) is in fact recomputation of each variable’s DOF coordinates from the new mesh coordinates — the storage is per-variable, the values are derived. The active_vars whitelist controls which of those recomputations runs now versus deferring to the next full rebuild.

property Gamma_P1¶

Projected P1 boundary normals as a sympy Matrix.

Returns the normalised, vertex-averaged PETSc face normals as a smooth P1 field. Preferred over Gamma_N for penalty and Nitsche BCs on curved boundaries — gives consistent orientation and better convergence in 3D.

Automatically updated when the mesh deforms.

update_lvec()[source]¶: This method creates and/or updates the mesh variable local vector. If the local vector is already up to date, this method will do nothing.

property lvec: Vec¶: Returns a local Petsc vector containing the flattened array of all the mesh variables.

register_remesh_hook(op)[source]¶

Called by the adapt op (smooth_mesh_interior, OT_adapt, follow_metric) after the generic per-variable REMAP pass. op must expose an on_remesh(ctx) method; ctx is a RemeshContext with the old/new coords, total displacement, dt, and a scratch dict for stashing things like v_mesh for the next solve.

Stored as a weak reference so operators that go out of scope are cleaned up automatically. Idempotent: registering the same operator twice is a no-op.

unregister_remesh_hook(op)[source]¶: Drop an operator’s on_remesh registration. Idempotent.

access(*writeable_vars)[source]¶

Dummy access manager that provides deferred sync for backward compatibility. Uses NDArray_With_Callback.delay_callbacks_global() internally.

This is a compatibility wrapper that allows existing code using the access() context manager to work with the new direct-access variable interfaces. All variable modifications are deferred and synchronized at context exit.

Parameters:: writeable_vars (MeshVariable) – Variables that will be modified (ignored - all variables are writable with the new interface, this parameter is kept for API compatibility)
Return type:: Context manager that defers variable synchronization until exit

Notes

This method is deprecated. New code should access variable.data or variable.array directly without requiring an access context.

property N: CoordSys3D¶

SymPy coordinate system for symbolic calculus.

The base coordinate system used for gradient, divergence, and curl operations. Access base scalars via mesh.N.x, mesh.N.y, mesh.N.z and base vectors via mesh.N.i, mesh.N.j, mesh.N.k.

Returns:: The SymPy coordinate system object.
Return type:: sympy.vector.CoordSys3D

See also

X: Coordinate system with data access.
r: Tuple of coordinate scalars.

property Gamma_N: MutableDenseMatrix¶

Normalised boundary/surface normal as a row matrix.

Returns Gamma / |Gamma| so that the result is a unit normal regardless of element size. Use this for penalty and Nitsche BCs where mesh-independent scaling is required.

Returns:: Row matrix of normalised boundary normal components.
Return type:: sympy.Matrix

property Gamma: MutableDenseMatrix¶

Raw (un-normalised) boundary coordinate scalars as a row matrix.

The magnitude scales with face edge length (2D) or face area (3D). For a unit normal, use Gamma_N instead.

Returns:: Row matrix of boundary coordinate scalars.
Return type:: sympy.Matrix

property X¶

Coordinate system with symbolic coordinates and data access.

The primary interface for mesh coordinates, providing both symbolic expressions for equations and numerical data for evaluation.

Returns:

Coordinate system object with:

mesh.X[0], mesh.X[1]: Symbolic coordinate functions
mesh.X.coords: Coordinate data array (vertex positions)
mesh.X.units: Coordinate units
x, y = mesh.X: Unpack symbolic coordinates

Return type:

CoordinateSystem

Examples

>>> x, y = mesh.X  # Symbolic coordinates for equations
>>> coords = mesh.X.coords  # Numerical vertex positions

See also

N: SymPy coordinate system for vector calculus.

property CoordinateSystem: CoordinateSystem¶: Alias for X (the coordinate system object).

property t¶

Symbolic time coordinate.

PETSc passes a time value (petsc_t) to all pointwise residual and Jacobian functions. Use mesh.t in expressions to reference this time without forcing JIT recompilation each timestep.

The low-level PETSc solver accepts time=t to set the value of petsc_t for pointwise functions. If not provided, petsc_t defaults to 0. Note: the high-level Python solve() wrappers do not yet pass time= through — set it directly via UW_DMSetTime at the Cython level if needed.

When the scaling system is active, mesh.t carries time units (derived from the model’s time scale) so that dimensional analysis works correctly in expressions.

Examples

>>> omega = 2 * np.pi / period
>>> stokes.add_dirichlet_bc((V0 * sympy.sin(omega * mesh.t), 0.0), "Top")
>>> stokes.solve(time=current_time)   # sets petsc_t before SNES

property nullspace_rotations¶

Symbolic velocity fields for rigid-body rotation null modes.

Returns a list of SymPy Matrix expressions in mesh Cartesian coordinates. Empty for meshes with no rotation nullspace (boxes, wedge segments with walls). Set by mesh factory functions for closed surfaces (annulus, spherical shell, etc.).

Each entry represents a rigid rotation: v = omega x r.

Returns:: Velocity fields for each independent rotation mode.
Return type:: list of sympy.Matrix

Examples

>>> annulus = uw.meshing.Annulus(...)
>>> annulus.nullspace_rotations  # [Matrix([-y, x])]
>>> shell = uw.meshing.SphericalShell(...)
>>> shell.nullspace_rotations   # 3 rotation matrices

property r: Tuple[BaseScalar]¶

Tuple of coordinate scalars \((x, y)\) or \((x, y, z)\).

Returns:: Tuple of SymPy base scalars (N.x, N.y[, N.z]).
Return type:: tuple

See also

rvec: Position vector form.

property rvec: Vector¶

Position vector \(\mathbf{r} = x\hat{i} + y\hat{j} [+ z\hat{k}]\).

Returns:: The position vector in the mesh coordinate system.
Return type:: sympy.vector.Vector

property data: ndarray¶

The array of mesh element vertex coordinates.

Deprecated since version 0.99.0: Use X.coords instead. mesh.data is deprecated in favor of mesh.X.coords (coordinate-system-aware interface).

This is an alias for mesh.points (which is also deprecated).

property points¶

Mesh node coordinates in physical units.

Deprecated since version 0.99.0: Use X.coords instead. mesh.points is deprecated in favor of mesh.X.coords (coordinate-system-aware interface).

When the mesh has coordinate scaling applied (via model units), this property automatically converts from internal model coordinates to physical coordinates for user access.

When the mesh has coordinate units specified, returns a unit-aware array.

Returns:: numpy.ndarray or UnitAwareArray: Node coordinates (with units if specified)

property physical_coordinates¶

Mesh coordinates in physical units.

Returns the mesh coordinate array scaled to physical units using the model’s length scale. Requires the mesh to be associated with a model that has reference quantities set.

Returns:: Coordinates in physical units, or None if no model scaling available
Return type:: UWQuantity or None

Examples

>>> model.set_reference_quantities(domain_length=1000*uw.units.km, ...)
>>> mesh = uw.meshing.StructuredQuadBox(...)
>>> physical_coords = mesh.physical_coordinates  # In kilometers

property physical_bounds¶

Mesh bounds in physical units.

Returns the mesh bounding box scaled to physical units using the model’s length scale.

Returns:: (min_coords, max_coords) in physical units, or None if no model scaling
Return type:: tuple of UWQuantity or None

Examples

>>> physical_min, physical_max = mesh.physical_bounds
>>> print(f"Domain: {physical_min} to {physical_max}")

property physical_extent¶

Mesh spatial extent in physical units.

Returns the mesh size (max - min) in each dimension scaled to physical units.

Returns:: Extent in physical units, or None if no model scaling
Return type:: UWQuantity or None

Examples

>>> extent = mesh.physical_extent
>>> print(f"Domain size: {extent}")

write_timestep(filename, index, outputPath='', meshVars=[], swarmVars=[], meshUpdates=False, create_xdmf=True)[source]¶

Write mesh and selected variables for visualisation output.

This writes: - one mesh HDF5 file (shared/static or per-step, depending on meshUpdates) - one HDF5 file per mesh variable - optional proxy files for swarm variables - optional XDMF file linking all output files

When create_xdmf=True (the default), variable files also include ParaView-compatible groups (/vertex_fields or /cell_fields), and an XDMF file is generated on rank 0.

Parameters:

filename (str)
index (int)
outputPath (str | None)
meshVars (list | None)
swarmVars (list | None)
meshUpdates (bool)
create_xdmf (bool)

petsc_save_checkpoint(index, meshVars=[], outputPath='')[source]¶

Save the mesh and mesh variables to HDF5 with XDMF.

This is a convenience wrapper around write_timestep() that provides the simpler interface used by earlier Underworld3 code. Output uses the same per-variable file layout and XDMF generation (including vertex/cell compatibility groups, field projection, and tensor repacking) as write_timestep().

Parameters:

meshVars (list | None) – List of UW mesh variables to save. If left empty then just the mesh is saved.
index (int) – An index which might correspond to the timestep or output number (for example).
outputPath (str | None) – Path to save the data. If left empty it will save the data in the current working directory.

write_checkpoint(filename, outputPath='', meshUpdates=True, meshVars=[], swarmVars=[], index=0, unique_id=False, separate_variable_files=True)[source]¶

Write PETSc DMPlex checkpoint files for restart/postprocessing.

Checkpoint output stores PETSc DMPlex section/vector metadata required for exact parallel reload. Unlike write_timestep(), this is restart output and does not write XDMF or vertex-field visualisation datasets.

Parameters:

filename (str) – Checkpoint base filename. With outputPath unset, this may include a directory. With outputPath set, it is joined to that directory.
outputPath (str) – Optional output directory, matching the write_timestep() style.
meshUpdates (bool) – If False, write the mesh checkpoint only when it does not already exist. If True, always write the indexed mesh checkpoint.
meshVars (list | None) – Variables to write into checkpoint files.
swarmVars (list | None) – Variables to write into checkpoint files.
index (int | None) – Checkpoint index used in output filenames.
unique_id (bool | None) – Preserve existing unique-rank filename behaviour for checkpoint data.
separate_variable_files (bool) – If True (default), write one file per variable: <base>.<variable>.<index>.h5. If False, write all variables into one file: <base>.checkpoint.<index>.h5.

snapshot_payload()[source]¶

Return a self-contained dict describing this mesh’s state.

The returned dict is consumed by underworld3.checkpoint.snapshot capture. Keys:

name: stable string identifier for the mesh.
mesh_version: current _mesh_version integer.
coords: deformed mesh coordinates (numpy array).
vars: {var.clean_name: gvec_array.copy()} for every mesh variable on this mesh.

v1.2 will additionally populate a topology key with section / DM-topology data sufficient to rebuild the DM on restore.

Return type:: dict

apply_snapshot_payload(payload)[source]¶

Restore this mesh from a payload produced by snapshot_payload().

v1 implementation writes coordinates and per-variable DOFs back in place. The captured DOF arrays must match the current section, which means _mesh_version must equal the captured value — mesh-adapt during the interval would have resized the section and is detected as a v1 refusal here.

v1.2 will replace the _mesh_version refusal with a rebuild-from-payload path: destroy the current DM, rebuild from payload["topology"], allocate vectors, write DOFs, and re-bind MeshVariable / Swarm wrappers. The interface stays the same; only this method’s body changes.

Parameters:: payload (dict)
Return type:: None

write(filename, index=None)[source]¶

Save mesh data to the specified hdf5 file.

Parameters:

filename (str) – The filename for the mesh checkpoint file.
index (int | None) – Not yet implemented. An optional index which might correspond to the timestep (for example).

vtk(filename)[source]¶

Save mesh to the specified file

Parameters:: filename (str)

generate_xdmf(filename)[source]¶

This method generates an xdmf schema for the specified file.

The filename of the generated file will be the same as the hdf5 file but with the xmf extension.

Parameters:: filename (str) – File name of the checkpointed hdf5 file for which the xdmf schema will be written.

property vars¶: A list of variables recorded on the mesh.

property block_vars¶: A list of variables recorded on the mesh.

points_in_domain(points, strict_validation=True)[source]¶

Determine if the given points lie in this domain. Uses a mesh-boundary skeletonization array to determine whether the point is inside the boundary or outside. If close to the boundary, it checks if points are in a cell.

Parameters:

points (array-like) – Coordinate array in any physical unit system (will be auto-converted). Plain numbers are assumed to be in model coordinates.
strict_validation (bool) – Whether to perform strict validation near boundaries

get_closest_cells(coords)[source]¶

This method uses a kd-tree algorithm to find the closest cells to the provided coords. For a regular mesh, this should be exactly the owning cell, but if the mesh is deformed, this is not guaranteed. Note, the nearest point may not be all that close by - use get_closest_local_cells to filter out points that are (probably) not within any local cell.

Parameters:¶

coords:: An array of the coordinates for which we wish to determine the closest cells. This should be a 2-dimensional array of shape (n_coords,dim) in any physical unit system (will be auto-converted). Plain numbers are assumed to be in model coordinates.

Returns:¶

closest_cells:: An array of indices representing the cells closest to the provided coordinates. This will be a 1-dimensional array of shape (n_coords).

Parameters:: coords (ndarray)
Return type:: ndarray

test_if_points_in_cells(points, cells, on_boundary=True, tol=0.0)[source]¶

Determine if the given points lie in the suggested cells. Uses a mesh skeletonization array to determine whether the point is with the convex polygon / polyhedron defined by a cell.

Exact if applied to a linear mesh, approximate otherwise.

Parameters:

points (array-like) – Coordinate array in any physical unit system (will be auto-converted)
cells (array-like) – Cell indices to test
on_boundary (bool, default True) – If True (the default), points exactly on a cell face count as inside the cell (natural for FE evaluation, where the basis at a shared face/vertex is consistent across adjacent cells). If False, points on the closure of a cell are reported as NOT in it (strict-inside semantics — useful when uniqueness matters).
tol (float, default 0.0) – Face-relative tolerance forwarded to _test_if_points_in_cells_internal. When > 0 takes precedence over on_boundary: the test admits points within tol of the face relative to the control-point separation² — used by the parallel evaluation locator for on-face / near-face queries at the mesh-spacing scale.

Returns:

Boolean array indicating if points are in cells

Return type:

numpy.ndarray

get_closest_local_cells(coords, on_boundary=True, tol=0.0)[source]¶

This method uses a kd-tree algorithm to find the closest cells to the provided coords. For a regular mesh, this should be exactly the owning cell, but if the mesh is deformed, this is not guaranteed. Also compares the distance from the cell to the point - if this is larger than the “cell size” then returns -1

Parameters:¶

coords:: An array of the coordinates for which we wish to determine the closest cells. This should be a 2-dimensional array of shape (n_coords,dim) in any physical unit system (will be auto-converted).
on_boundarybool, default True: If True (the default), queries exactly on a cell face are treated as inside that cell (natural for FE-evaluation hints — mesh vertices sit on cell faces by definition). If False, strict-inside semantics; boundary queries return -1.
tolfloat, default 0.0: Face-relative tolerance forwarded to _test_if_points_in_cells_internal. When > 0 takes precedence over on_boundary: the test admits points within tol of the face relative to the control-point separation² — used by the parallel evaluation locator for on-face / near-face queries at the mesh-spacing scale.

Returns:¶

closest_cells:: An array of indices representing the cells closest to the provided coordinates. This will be a 1-dimensional array of shape (n_coords).

Parameters:

coords (ndarray)
on_boundary (bool)
tol (float)

Return type:

ndarray

get_min_radius_old()[source]¶

This method returns the global minimum distance from any cell centroid to a face. It wraps to the PETSc DMPlexGetMinRadius routine. The petsc4py equivalent always returns zero.

Return type:: float

get_min_radius(**kwargs)¶: This method returns the global minimum distance from any cell centroid to a face. It wraps to the PETSc DMPlexGetMinRadius routine. The petsc4py equivalent always returns zero.

get_max_radius(**kwargs)¶: This method returns the global maximum distance from any cell centroid to a face.

get_mean_radius(**kwargs)¶

Global mean of the characteristic cell length scale (volume^(1/dim), i.e. the equivalent radius derived from each cell’s volume — the same quantity averaged by get_min_radius and get_max_radius to obtain global min/max). Parallel-safe via MPI allreduce of the local sum and count.

Together with get_min_radius() / get_max_radius() this is the canonical “mesh length” API. Use this anywhere you need a representative h0 (smoothing-length defaults, diffusion- stability heuristics, problem-scale normalisation) rather than reaching for the rank-local self._radii array, which gives different answers on different MPI ranks and leaks downstream (e.g. into JIT C source via per-rank pointwise-function inputs).

stats(uw_function, uw_meshVariable, basis=None)[source]¶

Returns various norms on the mesh for the provided function.

size
mean
min
max
sum
L2 norm
rms

NOTE: this currently assumes scalar variables !

meshVariable_mask_from_label(label_name, label_value)[source]¶: Extract single label value and make a point mask - note: this produces a mask on the mesh points and assumes a 1st order mesh. Cell labels are not respected in this function.

register_swarm(swarm)[source]¶: Register swarm as dependent on this mesh for coordinate change notifications

unregister_swarm(swarm)[source]¶: Unregister swarm (called during swarm cleanup)

register_surface(surface)[source]¶: Register surface as dependent on this mesh for adaptation notifications.

unregister_surface(surface)[source]¶: Unregister surface (called during surface cleanup).

OT_adapt(field, *, refinement=3.0, coarsening='auto', grad_smoothing_length='auto', metric_choice='front-following', fields_to_remap=None, fields_to_zero=None, skip_threshold=None, reference_coords=None, verbose=False)[source]¶

Adapt the mesh in place so cell sizes track |∇field|, using the validated optimal-transport reset pattern.

Each call resets the mesh to a cached reference (the initial uniform coordinates), FE-remaps field onto that clean canvas, builds a gradient-density metric, runs the OT mover, and FE-remaps the requested fields onto the adapted positions. Resetting every event (rather than composing adaptations across time steps) is what keeps the mover sliver-free over long runs. The “reset” is internal — from the caller’s point of view this just tracks the moving feature.

Topology is preserved (vertex count, DOF maps, rank partition unchanged); only coordinates move. Registered solvers are marked for rebuild via _deform_mesh.

Reference coordinates¶

The reset target is snapshotted lazily on the first call as self._ot_adapt_reference_coords (a copy of the current mesh.X.coords) and reused thereafter.

Warning

If the mesh is deformed by something other than OT_adapt between calls (e.g. a manual mesh._deform_mesh(...) or a resume that loads a deformed snapshot), the cached reference no longer matches the intended pristine state. Use OT_adapt_reset_reference() to re-baseline, or pass an explicit reference_coords for a one-off override.

param field:: Scalar field whose gradient drives refinement (typically T). Always FE-remapped onto the adapted mesh.
type field:: MeshVariable
param refinement:: Cell-size envelope h0/refinement for the densest cells. Validated range 1.5–5; 3 is the Nu sweet spot.
type refinement:: float, default 3.0
param coarsening:: "auto" = budget-conserving refinement**(1/d).
type coarsening:: float or “auto”, default “auto”
param grad_smoothing_length:: Screened-Poisson de-noising length for |∇field| before the metric is built — the most effective sliver lever; without it, production refinement chases sub-cell gradient noise. "auto" (default) ≈ the mesh’s uniform cell size (mean edge length) — the validated setting. None turns it off. A number or Pint length sets L explicitly; user-supplied lengths are unit-aware (non-dimensionalised via the projection), so pass a Pint quantity (or a non-dimensional number) — ≈ h0 is mild, ≈ 2·h0 stronger.
type grad_smoothing_length:: “auto”, None, float, or Pint Quantity, default “auto”
param metric_choice:
type metric_choice:: {“front-following”, “gradient-uniform”}, default “front-following”
param fields_to_remap:: Extra fields to FE-remap onto the adapted positions (field is always remapped). None ⇒ just field.
type fields_to_remap:: list of MeshVariable, optional
param fields_to_zero:: Fields to zero after the adapt (e.g. [V, P] for a cold restart of the flow solve).
type fields_to_zero:: list of MeshVariable, optional
param skip_threshold:: If the mesh is already aligned with the metric (misalignment below this; see mesh_metric_mismatch()), skip the adapt and return False. None ⇒ always adapt.
type skip_threshold:: float, optional
param reference_coords:: One-off override of the reset target (does not update the cache).
type reference_coords:: array, optional
param verbose:
type verbose:: bool, default False
returns:: True if the mesh was adapted, False if the skip_threshold check short-circuited it.
rtype:: bool

Notes

Boundary nodes slide tangentially and stay on the boundary for radial coordinate systems (Annulus / shell), using the projected boundary normal mesh.Gamma_P1. Cartesian boundaries are pinned (the vertex-evaluated normal is degenerate there).

Constrained-manifold meshes (mesh.dim != mesh.cdim, e.g. a 2D spherical surface in 3D) are not supported: the OT mover would have to constrain every node to the surface, not just boundary nodes. See docs/developer/design/ot-adapt-api-proposal.md.

Examples

>>> mesh = uw.meshing.Annulus(radiusOuter=1.0, radiusInner=0.5,
...                           cellSize=1/16, qdegree=3)
>>> T = uw.discretisation.MeshVariable("T", mesh, 1, degree=3)
>>> # ... initialise T ...
>>> mesh.OT_adapt(T, refinement=3.0, fields_to_remap=[T])

See also

OT_adapt_reset_reference: Re-baseline the reset reference coords.
underworld3.meshing.follow_metric: The single-shot anisotropic mover.
adapt: Topology-changing MMG remeshing (different mechanism).

OT_adapt_reset_reference(coords=None)[source]¶

Re-baseline the reference coordinates used by OT_adapt().

coords=None re-snapshots the current mesh.X.coords as the new reset target; passing explicit coords (e.g. the initial uniform mesh loaded from a checkpoint) sets those instead. Use on resume, when the loaded mesh is in a deformed state and the cache would otherwise lazily initialise from it.

adapt(metric_field, verbose=False)[source]¶

Adapt the mesh discretization based on a metric field.

This method refines or coarsens the mesh in place, automatically transferring all attached MeshVariables, updating Surfaces, and marking Solvers for rebuild on their next solve() call.

Parameters:

metric_field (MeshVariable) – A scalar MeshVariable containing metric values (1/h² where h is target edge length). Larger values mean finer mesh (smaller elements). Use Surface.refinement_metric() to create this field from distance.
verbose (bool, optional) – If True, print progress and statistics during adaptation.

Notes

The adaptation uses PETSc’s mesh adaptation with MMG/pragmatic backend.

What happens automatically:

MeshVariables are interpolated to the new mesh
Surfaces recompute their distance fields
Swarms are marked as stale (particle-element associations invalidated)
Solvers are marked for rebuild (happens lazily on next solve())

Examples

>>> # Define metric from fault distance
>>> metric = uw.discretisation.MeshVariable("H", mesh, 1)
>>> with mesh.access(metric):
...     # Smaller H near fault, larger far away
...     metric.data[:, 0] = 0.01 + 0.09 * fault.distance_from(mesh.data)
>>> mesh.adapt(metric, verbose=True)
>>> stokes.solve()  # Solver rebuilds automatically

MeshVariable¶

The primary class for field data on meshes.

underworld3.discretisation.MeshVariable¶: alias of EnhancedMeshVariable

Checkpointing¶

checkpoint_xdmf¶

underworld3.discretisation.checkpoint_xdmf(filename, meshUpdates=True, meshVars=[], swarmVars=[], index=0)[source]¶

Parameters:

filename (str)
meshUpdates (bool)
meshVars (list | None)
swarmVars (list | None)
index (int | None)