Discretisation

Mesh

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

class underworld3.discretisation.Mesh

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) – Deprecated and ignored (DeprecationWarning). Mesh coordinate units always come from the model’s reference quantities (model.set_reference_quantities); query them via mesh.units.

• 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)

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)

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)

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()

returns the break down of boundary labels from each processor

clone_dm_hierarchy()

Clone the dm hierarchy on the mesh

extract_region(label_name, label_value=None)

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

extract_surface(label_name, label_value=None, verbose=False)

Extract the codimension-1 surface marked by label_name as a mesh.

The third submesh flavour (alongside extract_region(), which filters cells of the same dimension): a surface submesh is a real Mesh for the parent’s codim-1 boundary stratum, sharing exact vertex positions with the parent. On a 3D SphericalShell, shell.extract_surface("Upper") returns a 2-manifold embedded in 3-space (dim = parent.dim - 1, cdim = parent.cdim).

Mechanism: DMPlexCreateSubmesh on the face label produces a cd-1 DM but retains an upward-DAG phantom stratum (one point per parent volume cell) that breaks closure-based navigation; DMPlexFilter on (depth, dim-1) strips it, leaving a clean standalone manifold. The two subpoint IS’s compose into a single surface→parent point map.

Parent ↔ submesh DOF transfer reuses restrict() / prolongate() (the same KDTree coordinate-match-at-1e-10 path as extract_region); surface vertices are an exact subset of the parent’s, so it is bit-exact.

Parameters:

• label_name (str) – Name of the parent boundary label whose marked faces become the cells of the surface submesh (e.g. "Upper").

• label_value (int, optional) – Stratum value within the label. If None, resolved from self.boundaries[label_name].value.

Returns:

A surface mesh with parent set to this mesh and the standard submesh lineage (subpoint_is, registration with _registered_submeshes).

Return type:

Mesh

Raises:

ValueError – If label_name is missing from the parent or its face stratum is empty (loud-fail contract — no degenerate-mesh fallback).

sync_coordinates_from_parent()

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')

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')

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)

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.

boundary_normal(boundary)

Outward unit normal of a single boundary, tracking deformation.

Assembles the EXACT, outward, area-weighted PETSc facet normals (dm.computeCellGeometryFVM) from ONLY this boundary’s facets onto its P1 vertices. Because each boundary is assembled independently, a vertex shared by two boundaries (a sharp corner) is NOT averaged across the discontinuity — each boundary keeps its own normal. On a smooth boundary (e.g. a free surface) the result is the smooth deformed normal. Cached per boundary; rebuilt lazily after a deform.

Returns a sympy Matrix (row) of the P1 normal-field components, for use as the constraint direction in Nitsche/penalty BCs.

Parameters:

boundary (str or enum) – Boundary label name (or a mesh.boundaries enum member).

cell_size()

Local, per-cell characteristic mesh size as a scalar field symbol.

Returns the .sym of a cell-constant (degree-0, discontinuous) scalar MeshVariable holding each cell’s characteristic length (the volume**(1/dim) equivalent radius, i.e. self._radii). Unlike the single global scalar from get_min_radius() (the smallest cell anywhere), this varies cell to cell, so a stabilisation that scales as \(1/h\) — e.g. the Nitsche free-slip penalty \(\gamma\mu/h\) — is correctly scaled on every facet of a non-uniform or adaptively-refined mesh rather than using the global minimum (which over-penalises coarse cells and drifts as refinement changes the global min).

On a boundary integral the kernel sees the value of the cell adjacent to the facet. The field is cached and rebuilt lazily; its data is refreshed when the mesh deforms or is adapted (see deform()), so it tracks a moving / re-refined mesh — a stale size on a deformed mesh would re-introduce the mis-scaling.

On a uniform mesh every cell is the same size, so this reduces to the global get_min_radius value everywhere and existing behaviour is preserved to tolerance.

Returns:

The cell-size field symbol, for use in JIT-compiled residuals.

Return type:

sympy scalar

property Gamma_P1

Deprecated — use Gamma in integrands and BCs, or boundary_normal() for a per-boundary P1 normal field.

This global field point-evaluates Gamma at every mesh vertex, but the underlying petsc_n only exists inside surface-integral kernels; off-kernel evaluation falls back to a coordinate-based direction and, on internal boundaries, averages oppositely-oriented facet normals into sub-unit vectors. Retained unchanged for back-compat (mesh-smoothing internals).

Automatically updated when the mesh deforms.

property boundary_normals

Declared analytic boundary normals (Enum or mapping), or None.

Assigning stamps the declaration against the mesh’s current geometry: any later coordinate change (deform, adapt, direct coordinate writes) marks it stale, and canonical_normal() then refuses to serve it. Re-assign after a deformation to re-declare normals that are valid for the new geometry (e.g. a radial normal after a radius-preserving remesh).

canonical_normal(boundary_name)

Analytic outward-pointing normal for a boundary, or None if no analytic normal was declared for that boundary.

Sourced from the mesh factory’s boundary_normals Enum: for axis-aligned box boundaries this is a constant sympy Matrix, for annulus / spherical-shell radial boundaries it is the analytic radial unit vector, and so on.

The primary caller is code that needs a partition-safe normal on an internal boundary — see :issue:`327`. On an internal boundary at a partition seam, PETSc’s per-quadrature petsc_n[] (surfacing as Gamma / Gamma_N) is derived from support[0] of the DMPlex facet closure, which is partition-dependent; different ranks disagree on which cell is “support[0]” for the one seam facet, and the outward normal of that facet flips sign. A signed integral of Gamma[k] is then wrong by O(seam-facets / total-facets). The analytic normal returned here is partition-independent and does not touch petsc_n, so it sidesteps the defect entirely for the mesh classes that know their internal-boundary geometry (BoxInternalBoundary, AnnulusInternalBoundary, SphericalShellInternalBoundary).

Parameters:

boundary_name (str) – Name of the boundary label to look up (case-sensitive; must match one of boundaries).

Returns:

Row matrix of length cdim giving the outward-pointing normal, or None if this mesh factory did not declare an analytic normal for boundary_name.

Return type:

sympy.Matrix or None

Raises:

RuntimeError – If the mesh coordinates have changed since the normals were declared (deform / adapt): the declaration describes the original geometry and may no longer match the deformed surface. Re-assign boundary_normals to re-declare normals valid for the current geometry.

See also

Gamma

the raw per-quadrature normal — use for external boundaries; may be partition-dependent on internal seams.

Gamma_N

normalised Gamma.

Gamma_P1

projected P1 normals, useful for curved external boundaries.

property return_coords_to_bounds

Callable coords -> coords returning out-of-domain points to the mesh.

Used by the semi-Lagrangian trace-back (systems.ddt) and by swarm advection to pull a departure point / particle that has left the domain back inside.

Two implementations, selected by the state of the geometry:

• analytic (default) — the closure installed by the mesh constructor (radial on an annulus/sphere, face clamps on a box). Cheap and exact while the boundary keeps the shape it was written for.

• general facet restore (once deform() has moved the geometry) — the nearest point on the mesh’s CURRENT boundary facets, with an outward-normal side test (_facet_return_coords_to_bounds()).

The switch matters for a moving free surface: the analytic closure is captured at construction with the ORIGINAL radii/extent, so on a deformed mesh it tests the wrong surface — where the surface has moved inward a point beyond the true boundary is not detected as outside, is never restored, and falls through to the evaluator’s RBF/Shepard fallback, which averages distant DOFs (hot material appearing in a cold boundary layer). The general restore follows the deformed boundary, so it stays correct.

Assigning to this attribute overrides both (the setter replaces the analytic closure and is honoured until the geometry deforms).

property bounding_surfaces

Mapping {label: BoundingSurface} of this mesh’s registered bounding-surface objects (tangent-slip + restore).

This is a NEW collection, separate from and additional to boundaries (the persisted gmsh/DMPlex label Enum, left untouched). Populated by analytic-geometry constructors (Annulus, SphericalShell, CubedSphere, box meshes); user-extendable via register_tangent_slip_provider().

register_tangent_slip_provider(label, surface)

Install a BoundingSurface object for a boundary label (separate from the persisted mesh.boundaries labelling).

Lets a user declare a custom analytic surface (e.g. an ellipsoid) the constructors don’t know about, or replace one (free-surface release).

restore_to_surface(coords, label)

Snap coords onto the named bounding surface (delegates to the surface object’s restore).

tangent_project(coords, label, reference)

Tangent-slide coords (displacement measured from reference) on the named bounding surface (delegates to the surface object).

boundary_slip(slip_spec=True, reference_coords=None, boundary_labels=None)

Build (is_pinned, project) for tangent slip on this mesh’s bounding surfaces — the orchestrator the metric movers call.

See docs/developer/design/boundary-slip-strategy.md. The mesh classifies which vertices slip vs pin (the cross-surface concern); each surface object owns its tangent-project + restore.

A vertex slips iff it lies on exactly one slip surface that has a registered analytic BoundingSurface; non-boundary, junction (≥2 surfaces), unregistered-surface, or degenerate-normal vertices are pinned (the step-1 safe default — facet restore is a follow-up).

Parameters:

• slip_spec – See _resolve_slip_spec(). Default True (all surfaces).

• reference_coords (ndarray, optional) – Fixed reference vertex positions (local-chart vertex order) that displacements are measured from. Defaults to mesh.X.coords.

• boundary_labels (iterable of str, optional) – Boundary labels defining the boundary (is_bnd). Defaults to all geometric boundary labels; pass a mover’s pinned set for parity.

Returns:

(is_pinned, project) – is_pinned is the per-vertex pinned mask (local-chart order); project(Y) slides+restores the slip vertices of Y in place and returns it.

Return type:

(ndarray[bool], callable)

project_to_slip_surface(coords, slip_spec=True, reference_coords=None, boundary_labels=None)

In-place convenience over boundary_slip(): slide+restore the slip vertices of coords (a full local-chart vertex array) and return it. For callers that just want coordinates snapped back (a checkpoint reload, a diagnostic, the free-surface module).

update_lvec(swarm_sync=True)

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.

swarm_sync=False skips the swarm-dependency hook below. It MUST be passed at call sites that only a SUBSET of ranks reach (e.g. the refresh calls inside petsc_interpolate — ranks with zero interior points skip that function entirely): the hook performs collective reductions, so running it on a subset of ranks deadlocks. Such sites rely on an earlier all-ranks update_lvec() for freshness, exactly as they already do for the collective globalToLocal below.

property lvec: Vec

Returns a local Petsc vector containing the flattened array of all the mesh variables.

register_remesh_hook(op)

Register an operator’s on_remesh(ctx) callback.

Called by the adapt op (smooth_mesh_interior, 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)

Drop an operator’s on_remesh registration. Idempotent.

ephemeral_coords()

Trusted scheme-internal trial coordinate moves, restored on exit.

For schemes (e.g. RK4 surface stages) that probe trial meshes to compute velocities/increments and must NOT commit a transfer. The coordinates are snapshotted on enter and restored on exit, so the intermediate meshes are genuinely ephemeral — only the final committed move (via deform()) updates fields + history.

deform(new_coords, *, dt=None, zero=None, verbose=False)

Move the mesh to new_coords, transferring all fields + history.

The public, foolproof way to impose an arbitrary node displacement (free surface, prescribed mesh motion). Wraps the remesh transaction so that REMAP variables are re-interpolated onto the new layout and every registered on_remesh hook fires — in particular the SemiLagrangian DDt’s coherent ALE (carry its history stack + a one-step v_mesh = Δx/dt correction consumed by the next solve).

Parameters:

• new_coords (ndarray) – Target vertex coordinates (shape of mesh.X.coords).

• dt (float, optional) – Timestep for the ALE v_mesh = Δx/dt correction. Required when SemiLagrangian history is present and the move is advective (a free-surface step); omit for a pure geometric re-mesh.

• zero (list of MeshVariable, optional) – Variables to zero after the move (e.g. V, P for a cold restart).

• verbose (bool)

Returns:

True if the mesh moved (geometry changed), False for a no-op.

Return type:

bool

access(*writeable_vars)

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

Deprecated alias — use Gamma.

Retained for back-compatibility: returns Gamma / |Gamma|, which is numerically identical to Gamma (the quadrature-point normal is already unit length).

Returns:

Row matrix of normalised boundary normal components.

Return type:

sympy.Matrix

property Gamma: MutableDenseMatrix

The boundary unit normal, as a symbolic row matrix.

This is the single user-facing normal symbol: use it in boundary integrands (uw.maths.BdIntegral) and natural boundary conditions on any boundary, external or internal.

The consumer that compiles the integrand resolves the symbol per boundary: on an external boundary the components map to PETSc’s exact per-quadrature outward unit normal (petsc_n[]); on an internal boundary — where the discrete facet normal has no well-defined orientation (issue #327) — they are substituted with the mesh factory’s declared analytic normal (see canonical_normal() and _resolve_boundary_normals()).

Returns:

Row matrix of boundary normal components.

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, petsc_reload=False)

Write mesh and selected variables for timestep output.

This is the standard mesh output method. It always writes:

• one mesh HDF5 file, shared across timesteps unless meshUpdates=True

• one HDF5 file per mesh variable

• raw coordinate/value datasets under /fields for coordinate-based reload with MeshVariable.read_timestep()

The optional payloads are controlled explicitly:

• create_xdmf=True writes ParaView/XDMF output. Variable files also receive /vertex_fields or /cell_fields compatibility groups, and rank 0 writes the companion .xdmf file.

• petsc_reload=True writes PETSc DMPlex section/vector metadata into the same per-variable HDF5 files. These files can then be loaded with MeshVariable.read_checkpoint() for PETSc-native same-mesh reload.

Common choices are:

• visualisation/remap only: create_xdmf=True, petsc_reload=False

• PETSc-native reload only: create_xdmf=False, petsc_reload=True

• unified visualisation/remap and PETSc reload: create_xdmf=True, petsc_reload=True

With both flags enabled, the same variable HDF5 file can be used by MeshVariable.read_timestep() for coordinate/KDTree remapping and by MeshVariable.read_checkpoint() for exact PETSc-native reload.

Parameters:

• filename (str) – Output filename base. Files are written as <filename>.mesh.<index>.h5 and <filename>.mesh.<variable>.<index>.h5.

• index (int) – Timestep/output index used in generated filenames.

• outputPath (str | None) – Directory where output files are written.

• meshVars (list | None) – Mesh variables to write.

• swarmVars (list | None) – Swarm variables to write as proxy fields.

• meshUpdates (bool) – If False, reuse <filename>.mesh.00000.h5 when it already exists. If True, write an indexed mesh file for this timestep.

• create_xdmf (bool) – Write ParaView/XDMF-compatible datasets and companion XDMF file.

• petsc_reload (bool) – Write PETSc DMPlex section/vector metadata for reload with MeshVariable.read_checkpoint().

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

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, create_xdmf=False)

Compatibility wrapper for PETSc DMPlex reload output.

This method is retained for existing callers. New code should use write_timestep(..., petsc_reload=True) so all mesh-variable output goes through the standard timestep writer. By default this compatibility method writes PETSc DMPlex section/vector metadata required for exact parallel reload and does not write XDMF or vertex-field visualisation datasets. Use create_xdmf=True to route through the unified timestep-style output path.

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.

• create_xdmf (bool) – If True, route through write_timestep() and write XDMF, vertex/cell compatibility groups, coordinate/KDTree remap data, and PETSc reload metadata. The output uses the timestep filename convention <base>.mesh.<variable>.<index>.h5. This mode does not support unique_id=True or separate_variable_files=False.

snapshot_payload()

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)

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, petsc_format=None)

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).

• petsc_format (bool | None) – If True, force PETSc DMPlex HDF5 checkpoint/restart topology. If False, force PETSc HDF5_VIZ topology only. If None, use PETSc’s default HDF5 layout, which includes the restart-style topology and labels as well as visualization topology for XDMF.

vtk(filename)

Save mesh to the specified file

Parameters:

filename (str)

generate_xdmf(filename)

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)

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)

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)

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)

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()

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)

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)

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)

Register swarm as dependent on this mesh for coordinate change notifications

unregister_swarm(swarm)

Unregister swarm (called during swarm cleanup)

register_surface(surface)

Register surface as dependent on this mesh for adaptation notifications.

unregister_surface(surface)

Unregister surface (called during surface cleanup).

OT_adapt(*args, **kwargs)

RETIRED (2026-07). The optimal-transport reset adapt was superseded by the variational MMPDE mover; calling this raises RuntimeError.

Use uw.meshing.follow_metric(mesh, field, refinement=...) (the two-knob gradient-following adapter) or uw.meshing.smooth_mesh_interior(mesh, metric=..., method="mmpde") for fixed-topology node redistribution; use adapt() when more resolution than a fixed node budget can provide is needed.

redistribute_nodes(metric, *, verbose=False, **kwargs)

Move this mesh’s nodes so cell sizes follow metric, with the topology fixed.

Node redistribution concentrates the existing node budget where the metric demands resolution: vertex count, vertex ids, DOF layout and the parallel partition are all preserved — only coordinates move. Mesh variables (and solver / transport history) are transferred onto the moved nodes automatically. Contrast adapt(), which refines (returns a child mesh with more cells), and remesh(), which regenerates the mesh in place.

This method is how each mesh type controls whether (and how) it can be modified: the base implementation supports 2D (triangle) and 3D (tetrahedral) simplex meshes, where it drives the Huang–Kamenski variational MMPDE mover (non-folding by construction, parallel-safe; scalar metric → isotropic equidistribution, tensor metric → anisotropic clustering and alignment). Quadrilateral / hexahedral meshes and constrained manifolds raise NotImplementedError — no mover is implemented for them yet.

Parameters:

• metric (sympy expression, MeshVariable, or sympy Matrix) – Target density \(\rho(x)\) (scalar, larger ⇒ finer cells) or a full \(d \times d\) SPD metric tensor (anisotropic: small across a feature, base along it).

• verbose (bool, default False) – Print mover progress.

• **kwargs – Forwarded to underworld3.meshing.smooth_mesh_interior() — e.g. pinned_labels, slip_surfaces, skip_threshold, method_kwargs (the mover tunables).

Examples

>>> x, y = mesh.CoordinateSystem.X
>>> rho = 1 + 8 * sympy.exp(-(((x - 0.5)**2 + (y - 0.5)**2) / 0.05))
>>> mesh.redistribute_nodes(rho)

See also

underworld3.meshing.node_redistribution

The free-function spelling of this operation.

underworld3.meshing.follow_metric

Two-knob adapter that builds the metric from a field gradient.

adapt

Add resolution (topology change, returns a child mesh).

relax(metric=None, *, verbose=False, **kwargs)

Improve this mesh’s element shapes without changing its size distribution or its topology.

The companion to adapt(). Refinement chooses where new nodes go from combinatorics — which edge the tagging rule nominated (bisection), or the cell centroid (Alfeld) — never from geometry, so a refined mesh carries needles and slivers that reflect the base mesh’s arbitrary choices rather than anything about the problem. Relaxation moves those nodes to where the geometry wants them, keeping the resolution the refinement installed:

child = mesh.adapt(metric, max_levels=3)
child.relax()

Implemented with the same MMPDE mover as redistribute_nodes(), in its ideal reference frame (reference="ideal"): each cell’s reference element is a regular simplex scaled to that cell’s own current volume. Two consequences, and they are the whole point:

• the size term starts and stays at its optimum, so the graded spacing is preserved rather than re-derived;

• “distortion” is measured against equilateral, not against the mesh as supplied — so a distorted mesh is no longer its own optimum, which is exactly why redistribute_nodes(metric) cannot do this job (its reference IS the mesh it was handed, so under a uniform metric it moves nothing at all).

Non-folding by construction, and the parallel partition, vertex count and DOF layout are unchanged.

Moving the nodes can upset the custom-P geometric-MG transfers, which are built by geometric point location rather than from the refinement relation: the local-support barycentric builder can be left with a coarse DOF that has no fine image. The FMG build now retries with the global-support rbf builder before giving up, so the hierarchy survives (measured in 3D: GAMG at 23 iterations before the retry, pc=mg at 2 after). See #424.

Two valid placements, neither dominant. adapt(metric, ..., relax=True) relaxes once at the end — the recommended default. relax="per-generation" relaxes inside the refinement loop, so each generation marks from already-relaxed coordinates. Both beat no relaxation; which is better depends on which property you weight, and we deliberately do not rank them:

property	unrelaxed	at end	per generation
P1 interpolation error	2.75e-2	2.56e-2	2.38e-2
99th pct max angle	112 deg	110 deg	96 deg
on-fault size spread	2.80	2.31	2.31
far-field closure halo	42.8	27.7	43.9

Relaxing at the end keeps the cell count identical to the unrelaxed mesh and cuts the far-field halo most; relaxing per generation gives the cleanest element shapes and the lowest error, but spends ~3% more cells to do it.

In 3D it improves mesh QUALITY but not interpolation error — two different things, and worth keeping apart. On an adapted 3D mesh it halves the near-degenerate population (cells with q < 0.1: 3.6% -> 1.8%), lifts median quality 0.32 -> 0.39 and pulls the 99th percentile dihedral angle back from 153 to 146 degrees; but the interpolation error of an isotropic feature is unchanged (+0.5%). That is consistent: relaxation holds the size distribution, and in 2D the error gain came from cells ALIGNING onto the feature, which an isotropic metric gives no reason to do in 3D. Use it in 3D for conditioning and element quality, not expecting an accuracy win.

Parameters:

• metric (sympy expression, MeshVariable, or sympy Matrix, optional) –

  Usually omitted, and omitting it is what makes this a shape guarantee. None (default) relaxes under a uniform metric — pure shape repair at fixed size.

  Passing a metric switches to the ideal-metric frame, which re-grades as well as reshapes, and it will trade element shape away to chase the size field. Measured on a 4-level graded box, the 99th-percentile max angle went 117.9 -> 113.8 degrees with no metric but 117.9 -> 127.4 with one. Pass a metric when you want the sizes corrected too, and accept that shape is no longer the objective.

• verbose (bool, default False) – Print mover progress.

• **kwargs – Forwarded to redistribute_nodes() — e.g. pinned_labels, slip_surfaces, method_kwargs (mover tunables such as n_outer).

See also

adapt

Add resolution (topology change, returns a child mesh).

redistribute_nodes

Move nodes to follow a metric (changes the size distribution; the reference is the mesh as supplied).

adapt(metric_field, max_levels=None, node_budget=None, builder=None, adapter=None, engine=None, verbose=False, relax=False, relax_kwargs=None)

Nested adapt-on-top: return a refined child mesh.

Locally refine the static base finest where the metric demands resolution, on top of the existing uniform hierarchy, and return a new child mesh (child.parent is self). The base mesh is not modified — this is adapt / re-adapt, not node movement: each call re-marks from the static base finest, so successive adapts are non-cumulative (cf. remesh(), which regenerates the mesh in place via MMG and may redistribute).

The default call needs no engine choice — mesh.adapt(metric, max_levels=...) refines with the graded newest-vertex-bisection engine, on 2D triangle and 3D tetrahedral meshes, serial and parallel (the refined mesh is partition-independent: the same mesh at any communicator size). engine= remains available as an advanced / internal selector (the algorithm names live here, not in the everyday call):

• "nvb" (default) — newest-vertex bisection, a graded engine with a bounded conforming closure: a marked cell adds O(1) cells locally, so successive levels grade (a level+1 ring around a finer core) and DOFs concentrate near the feature. Runs in parallel in both dimensions via the native uwnvb DMPlexTransform driver (in-place, co-partitioned with the parent, bit-confluent serial↔parallel; in 3D the per-cell refinement state is seeded identically on every rank from geometry). When the compiled extension is absent it falls back to the serial cell-list engines (NVBMesh / TaggedBisectionMesh; NotImplementedError at np>1). Bisects 1→2 (volume halves), so one isotropic-equivalent max_levels is run as dim bisection generations.

• "sbr" — PETSc skeleton-based (longest-edge) bisection, 2D only (PETSc’s SBR transform cannot handle tetrahedra). Each pass refines marked cells isotropically (1→4). Its conforming closure is unbounded for region marking, so it produces a uniform-finest patch, not a graded mesh (a marked cell drains the longest-edge path to the patch edge). Robust and fine for the MG hierarchy.

The child owns a custom-P geometric-MG hierarchy with one level per refinement step — [base L0 … base finest, refine-1, …, refine-n] — so every solver built on it drives geometric multigrid on the refined operator with no per-solver setup. (Each max_levels pass adds its own MG level; the transfers between consecutive levels each span a single refinement.)

Note

mesh.adapt(metric) with no other keyword raises a TypeError: that call shape used to be the in-place MMG remesher (now remesh()), and a legacy caller would silently discard the returned child. Pass any keyword — e.g. adapt(metric, max_levels=2) — to opt in to the new semantics.

Parameters:

• metric_field (MeshVariable, sympy expression, or callable) – Scalar metric M = 1/h² (target edge length h); larger ⇒ finer. A MeshVariable or sympy/UWexpression is sampled through uw.function.evaluate (so anything it references is interpolated from the base mesh). A callable metric(centroids) -> M is evaluated directly at each refined level’s centroids — use this for a metric built from exact geometry (e.g. Surface.refinement_metric_function()) so a thin feature refines to a clean, uniform-width band instead of a P1-aliased patchy one. Same interface as remesh() / adaptivity.create_metric.

• max_levels (int, default 2) – Maximum refinement depth applied on top of the base finest (bounds the on-rank imbalance). Each level re-marks against the metric.

• node_budget (int or None) – Optional cap on the number of seed cells marked per level (highest- metric first). Caveat: this caps the marked seeds, not the resulting DOFs — SBR’s conforming closure re-refines the whole connected patch from any seed in it, so a seed cap does not bound the added DOFs and cannot, on its own, concentrate the finest level near a feature. To make the finest level hug a feature (a funnel) with bounded per-level growth, shape the metric so element size grows with distance (e.g. Surface.refinement_metric(..., profile="linear"), a wedge), not this budget. A flat-core metric (e.g. a Gaussian) refines the whole core uniformly at every level (no funnel). None ⇒ uncapped.

• builder ({"barycentric", "rbf"}) – Per-level node-prolongation builder for the child’s custom-P hierarchy.

• adapter ({"sbr", "mmg"}) – "sbr" (default) is the nested adapt-on-top path (the refinement engine is then chosen by engine). "mmg" is a deprecated shim that forwards to remesh() (in-place, returns self).

• engine ({"nvb", "sbr"}, optional) – Advanced selector for the nested refinement engine (ignored when adapter="mmg"). Default "nvb" — graded newest-vertex bisection; "sbr" is longest-edge bisection (uniform patch, still the right choice when a uniform-finest MG patch is wanted). See above.

• verbose (bool)

Returns:

The refined child (or self when adapter='mmg').

Return type:

Mesh

Notes

Controlling the grading (funnel toward a feature). Each SBR pass refines every cell that is still coarser than the metric target, so the final grading is whatever the metric specifies. To make the finest level hug a feature (a funnel), use a metric whose target size grows with distance from the feature — a wedge, e.g. Surface.refinement_metric(h_near, h_far, width, profile="linear") with a small width. A flat-core metric (Gaussian) instead refines the whole core uniformly at every level, so every level has the same width.

SBR’s conforming closure re-refines the whole connected patch from any marked seed, so the funnel must come from the metric shape, not from capping seed counts (see node_budget). For a 1-D feature in 2-D the per-level DOFs still grow (the along-feature resolution doubles each level); the wedge keeps the total finite and feature-concentrated.

remesh(metric_field, verbose=False)

Re-mesh (regenerate) the discretization in place from a metric field.

This is the MMG / topology-changing remesher: it regenerates the mesh in place (cells are created/destroyed and the partition may change), automatically transferring all attached MeshVariables, updating Surfaces, and marking Solvers for rebuild on their next solve() call.

Contrast adapt(), which performs nested skeleton-based refinement on top of the static base and returns a refined child (the mesh is not modified in place). remesh is the in-place, redistributing path; prefer adapt when you want a parent/child geometric-MG hierarchy.

This method was formerly called adapt; adapt now performs the nested SBR adapt-on-top.

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)
>>> # Smaller H near fault, larger far away
>>> metric.data[:, 0] = 0.01 + 0.09 * fault.distance_from(mesh.X.coords)
>>> mesh.remesh(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)

Parameters:

• filename (str)

• meshUpdates (bool)

• meshVars (list | None)

• swarmVars (list | None)

• index (int | None)