Meshing¶

Underworld3 Meshing Module

This module provides mesh generation functions organized by geometry type. All functions maintain backward compatibility with the original meshing interface.

Cartesian Meshes¶

StructuredQuadBox¶

underworld3.meshing.StructuredQuadBox(elementRes=(16, 16), minCoords=None, maxCoords=None, degree=1, qdegree=2, filename=None, refinement=None, gmsh_verbosity=0, units=None, verbose=False)[source]¶

Create a structured quadrilateral or hexahedral box mesh.

Generates a mesh with regular rectangular (2D) or brick (3D) elements using transfinite meshing. Provides precise control over element count in each direction.

Parameters:

elementRes (tuple of int, default=(16, 16)) – Number of elements in each direction. Use (nx, ny) for 2D or (nx, ny, nz) for 3D. This tuple also determines the mesh dimensionality.
minCoords (tuple of float, optional) – Minimum corner coordinates. Defaults to (0.0, 0.0) for 2D or (0.0, 0.0, 0.0) for 3D based on elementRes length. Supports plain numbers or UWQuantity objects.
maxCoords (tuple of float, optional) – Maximum corner coordinates. Defaults to (1.0, 1.0) for 2D or (1.0, 1.0, 1.0) for 3D. Supports UWQuantity objects.
degree (int, default=1) – Polynomial degree of finite element basis functions. Use degree=1 for bilinear/trilinear elements, degree=2 for biquadratic/triquadratic.
qdegree (int, default=2) – Quadrature degree for numerical integration.
filename (str, optional) – Path to save the mesh file. If None, auto-generates in .meshes/.
refinement (int, optional) – Number of uniform refinement levels to apply.
gmsh_verbosity (int, default=0) – Gmsh output verbosity level.
units (str, optional) – Deprecated. Mesh coordinates are always in model reference units.
verbose (bool, default=False) – Print diagnostic information during mesh construction.

Returns:

An Underworld mesh object with structured elements and boundaries:

2D boundaries (same as UnstructuredSimplexBox):

Bottom: \(y = y_{min}\) edge
Top: \(y = y_{max}\) edge
Right: \(x = x_{max}\) edge
Left: \(x = x_{min}\) edge

3D boundaries:

Bottom: \(z = z_{min}\) face
Top: \(z = z_{max}\) face
Right: \(x = x_{max}\) face
Left: \(x = x_{min}\) face
Front: \(y = y_{min}\) face
Back: \(y = y_{max}\) face

Return type:

Mesh

See also

UnstructuredSimplexBox: For triangular/tetrahedral meshes.
BoxInternalBoundary: For box meshes with an internal interface.

Examples

Create a 2D structured mesh with 32x32 elements:

>>> import underworld3 as uw
>>> mesh = uw.meshing.StructuredQuadBox(
...     elementRes=(32, 32),
...     minCoords=(0.0, 0.0),
...     maxCoords=(1.0, 1.0)
... )

Create a 3D mesh (note the 3-element tuple):

>>> mesh3d = uw.meshing.StructuredQuadBox(
...     elementRes=(16, 16, 8),
...     maxCoords=(2.0, 2.0, 1.0)
... )

Notes

Structured meshes have predictable element layouts which can be advantageous for:

Consistent interpolation behaviour
Benchmark problems with known analytical solutions
Simpler mesh-to-mesh comparisons in convergence studies

The mesh dimensionality is determined by the length of elementRes: 2-tuple creates a 2D mesh, 3-tuple creates a 3D mesh.

UnstructuredSimplexBox¶

underworld3.meshing.UnstructuredSimplexBox(minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.1, degree=1, qdegree=2, regular=False, filename=None, refinement=None, gmsh_verbosity=0, units=None, verbose=False)[source]¶

Create an unstructured simplex mesh on a rectangular box domain.

Generates a triangular (2D) or tetrahedral (3D) mesh using Gmsh, with named boundary labels for applying boundary conditions.

Parameters:

minCoords (tuple of float) – Minimum corner coordinates (x_min, y_min) for 2D or (x_min, y_min, z_min) for 3D. Supports plain numbers (model units) or UWQuantity objects (auto-converted).
maxCoords (tuple of float) – Maximum corner coordinates (x_max, y_max) for 2D or (x_max, y_max, z_max) for 3D. Supports plain numbers or UWQuantity objects.
cellSize (float) – Target mesh element size. Controls mesh density; smaller values produce finer meshes with more elements.
degree (int, default=1) – Polynomial degree of finite element basis functions. Use degree=1 for linear elements, degree=2 for quadratic.
qdegree (int, default=2) – Quadrature degree for numerical integration. Should typically be at least 2 * degree for accuracy.
regular (bool, default=False) – If True, use transfinite meshing for a more structured layout. Currently only works for 2D meshes.
filename (str, optional) – Path to save the mesh file. If None, generates a unique name in the .meshes/ directory based on mesh parameters.
refinement (int, optional) – Number of uniform refinement levels to apply after mesh generation. Each level approximately quadruples element count.
gmsh_verbosity (int, default=0) – Gmsh output verbosity level. 0 is silent, higher values produce more diagnostic output.
units (str, optional) – Deprecated. Mesh coordinates are always in model reference units. This parameter is retained for backward compatibility.
verbose (bool, default=False) – If True, print additional diagnostic information during mesh construction.

Returns:

An Underworld mesh object with the following boundaries defined:

2D boundaries (accessible via mesh.boundaries):

Bottom: \(y = y_{min}\) edge
Top: \(y = y_{max}\) edge
Right: \(x = x_{max}\) edge
Left: \(x = x_{min}\) edge

3D boundaries:

Bottom: \(z = z_{min}\) face
Top: \(z = z_{max}\) face
Right: \(x = x_{max}\) face
Left: \(x = x_{min}\) face
Front: \(y = y_{min}\) face
Back: \(y = y_{max}\) face

Return type:

Mesh

See also

StructuredQuadBox: For quadrilateral/hexahedral meshes.
BoxInternalBoundary: For box meshes with an internal interface.

Examples

Create a 2D unit square mesh:

>>> import underworld3 as uw
>>> mesh = uw.meshing.UnstructuredSimplexBox(
...     minCoords=(0.0, 0.0),
...     maxCoords=(1.0, 1.0),
...     cellSize=0.1
... )
>>> mesh.dim
2

Create a 3D box with finer resolution:

>>> mesh3d = uw.meshing.UnstructuredSimplexBox(
...     minCoords=(0.0, 0.0, 0.0),
...     maxCoords=(2.0, 1.0, 1.0),
...     cellSize=0.05,
...     degree=2
... )

Access boundary labels for boundary conditions:

>>> mesh.boundaries.Bottom
<boundaries_2D.Bottom: 11>

Notes

Mesh coordinates are always in non-dimensional (scaled) units (set via uw.scaling.get_coefficients()). If UWQuantity objects with physical units are passed, they are automatically converted using uw.scaling.non_dimensionalise().

The regular=True option produces a more structured mesh layout but currently only works for 2D meshes.

BoxInternalBoundary¶

underworld3.meshing.BoxInternalBoundary(elementRes=(8, 8, 8), zelementRes=(4, 4), cellSize=0.1, minCoords=(0, 0, 0), maxCoords=(1, 1, 1), zintCoord=0.5, simplex=False, degree=1, qdegree=2, filename=None, refinement=None, gmsh_verbosity=0, units=None, verbose=False)[source]¶

Create a box mesh with an internal horizontal boundary.

Generates a 2D or 3D mesh with an embedded internal boundary surface, useful for problems with material interfaces, phase boundaries, or layered domains that require flux calculations across the interface.

Parameters:

elementRes (tuple of int, default=(8, 8, 8)) – Number of elements in each direction (nx, ny) for 2D or (nx, ny, nz) for 3D. Used when simplex=False (structured).
zelementRes (tuple of int, default=(4, 4)) – Number of elements (n_below, n_above) in the vertical direction below and above the internal boundary. Allows different resolution in each layer.
cellSize (float, default=0.1) – Target element size for unstructured meshing (simplex=True). Ignored for structured meshes.
minCoords (tuple of float, default=(0, 0, 0)) – Minimum corner coordinates. Length determines dimensionality: 2-tuple for 2D, 3-tuple for 3D.
maxCoords (tuple of float, default=(1, 1, 1)) – Maximum corner coordinates.
zintCoord (float, default=0.5) – Vertical coordinate of the internal boundary surface. In 2D this is the y-coordinate; in 3D the z-coordinate.
simplex (bool, default=False) – If False, create a structured quadrilateral/hexahedral mesh. If True, create an unstructured triangular/tetrahedral mesh.
degree (int, default=1) – Polynomial degree of finite element basis functions.
qdegree (int, default=2) – Quadrature degree for numerical integration.
filename (str, optional) – Path to save the mesh file. If None, auto-generates in .meshes/.
refinement (int, optional) – Number of uniform refinement levels to apply.
gmsh_verbosity (int, default=0) – Gmsh output verbosity level.
units (str, optional) – Coordinate units for unit-aware arrays.
verbose (bool, default=False) – Print diagnostic information during mesh construction.

Returns:

An Underworld mesh object with boundaries including an internal interface:

2D boundaries:

Bottom: \(y = y_{min}\) edge
Top: \(y = y_{max}\) edge
Right: \(x = x_{max}\) edge
Left: \(x = x_{min}\) edge
Internal: \(y = z_{int}\) interface

3D boundaries:

Bottom: \(z = z_{min}\) face
Top: \(z = z_{max}\) face
Right: \(x = x_{max}\) face
Left: \(x = x_{min}\) face
Front: \(y = y_{min}\) face
Back: \(y = y_{max}\) face
Internal: \(z = z_{int}\) interface

Return type:

Mesh

See also

UnstructuredSimplexBox: Box mesh without internal boundary.
AnnulusInternalBoundary: Annular mesh with internal boundary.

Examples

Create a 2D layered domain with an interface at y=0.5:

>>> import underworld3 as uw
>>> mesh = uw.meshing.BoxInternalBoundary(
...     minCoords=(0.0, 0.0),
...     maxCoords=(1.0, 1.0),
...     zintCoord=0.5,
...     elementRes=(16, 16),
...     zelementRes=(8, 8)
... )

Access the internal boundary for flux calculations:

>>> mesh.boundaries.Internal
<boundaries_2D.Internal: 15>

Notes

The internal boundary is useful for:

Calculating heat flux across a thermal boundary layer
Tracking mass flux between mantle and crust
Applying different material properties in each layer

Annulus Meshes¶

Annulus¶

underworld3.meshing.Annulus(radiusOuter=1.0, radiusInner=0.547, cellSize=0.1, cellSizeOuter=None, cellSizeInner=None, centre=False, degree=1, qdegree=2, filename=None, refinement=None, gmsh_verbosity=0, verbose=False)[source]¶

Create a full 2D annulus mesh.

Generates a ring-shaped mesh (or disc if radiusInner=0) using triangular elements. Provides \((r, \theta)\) coordinates with radial boundary normal directions.

Parameters:

radiusOuter (float, default=1.0) – Outer radius of the annulus. Supports UWQuantity objects.
radiusInner (float, default=0.547) – Inner radius of the annulus. Set to 0 for a full disc.
cellSize (float, default=0.1) – Default target mesh element size.
cellSizeOuter (float, optional) – Element size at the outer boundary. Defaults to cellSize.
cellSizeInner (float, optional) – Element size at the inner boundary. Defaults to cellSize. Use different sizes to create graded meshes.
centre (bool, default=False) – If True and radiusInner=0, mark the centre point as a boundary for applying point constraints.
degree (int, default=1) – Polynomial degree of finite element basis functions.
qdegree (int, default=2) – Quadrature degree for numerical integration.
filename (str, optional) – Path to save the mesh file.
refinement (int, optional) – Number of uniform refinement levels to apply. Each level approximately quadruples element count.
gmsh_verbosity (int, default=0) – Gmsh output verbosity level.
verbose (bool, default=False) – Print diagnostic information.

Returns:

A 2D mesh with boundaries:

Lower: Inner boundary at \(r = r_{inner}\)
Upper: Outer boundary at \(r = r_{outer}\)
Centre: Centre point (if \(r_{inner} = 0\))

The mesh includes a refinement callback that snaps boundary nodes back to true circular geometry after refinement.

Return type:

Mesh

See also

QuarterAnnulus: Partial annulus with radial edges.
AnnulusInternalBoundary: Annulus with an internal interface.
SphericalShell: 3D equivalent for spherical geometries.

Examples

Create a standard annulus for mantle convection:

>>> import underworld3 as uw
>>> mesh = uw.meshing.Annulus(
...     radiusOuter=1.0,
...     radiusInner=0.547,
...     cellSize=0.05
... )

Create a graded mesh (finer at inner boundary):

>>> mesh = uw.meshing.Annulus(
...     radiusOuter=1.0,
...     radiusInner=0.5,
...     cellSizeOuter=0.1,
...     cellSizeInner=0.02
... )

Notes

The inner radius default of 0.547 corresponds approximately to the Earth’s core-mantle boundary radius ratio (3480/6371).

The mesh coordinate system provides unit vectors via mesh.CoordinateSystem:

unit_e_0: radial direction \((r)\)
unit_e_1: tangential direction \((\theta)\)

For free-slip boundary conditions on vector problems (e.g., Stokes), use a penalty on the normal velocity component:

Gamma_N = mesh.Gamma  # discrete normal, or
Gamma_N = mesh.CoordinateSystem.unit_e_0  # analytic radial
stokes.add_natural_bc(
    penalty * Gamma_N.dot(v_soln.sym) * Gamma_N, "Upper"
)

QuarterAnnulus¶

underworld3.meshing.QuarterAnnulus(radiusOuter=1.0, radiusInner=0.547, angle=45, cellSize=0.1, centre=False, degree=1, qdegree=2, filename=None, gmsh_verbosity=0, verbose=False)[source]¶

Create a quarter-annulus (wedge) mesh in 2D.

Generates a pie-slice shaped mesh bounded by inner and outer circular arcs and two radial edges. Provides \((r, \theta)\) coordinates for convenient representation of radial problems.

Parameters:

radiusOuter (float, default=1.0) – Outer radius of the annulus. Supports UWQuantity objects for automatic unit conversion.
radiusInner (float, default=0.547) – Inner radius of the annulus. Set to 0 for a disc sector extending to the centre.
angle (float, default=45) – Angular extent of the wedge in degrees. The wedge spans from -angle to +angle (total sweep = 2 * angle).
cellSize (float, default=0.1) – Target mesh element size.
centre (bool, default=False) – If True and radiusInner=0, mark the centre point as a boundary for applying point constraints.
degree (int, default=1) – Polynomial degree of finite element basis functions.
qdegree (int, default=2) – Quadrature degree for numerical integration.
filename (str, optional) – Path to save the mesh file.
gmsh_verbosity (int, default=0) – Gmsh output verbosity level.
verbose (bool, default=False) – Print diagnostic information.

Returns:

A 2D mesh with boundaries:

Lower: Inner arc at \(r = r_{inner}\)
Upper: Outer arc at \(r = r_{outer}\)
Left: Left radial edge
Right: Right radial edge
Centre: Centre point (if \(r_{inner} = 0\))

Return type:

Mesh

See also

Annulus: Full 360-degree annulus.
SegmentofAnnulus: Partial annulus with arbitrary angle.

Examples

Create a quarter-disc (no inner hole):

>>> import underworld3 as uw
>>> mesh = uw.meshing.QuarterAnnulus(
...     radiusOuter=1.0,
...     radiusInner=0.0,
...     angle=45,
...     cellSize=0.05
... )

Notes

The mesh coordinate system provides unit vectors via mesh.CoordinateSystem:

unit_e_0: radial direction \((r)\)
unit_e_1: tangential direction \((\theta)\)

SegmentofAnnulus¶

underworld3.meshing.SegmentofAnnulus(radiusOuter=1.0, radiusInner=0.547, angleExtent=45, cellSize=0.1, centre=False, degree=1, qdegree=2, filename=None, refinement=None, gmsh_verbosity=0, verbose=False)[source]¶

Generates a segment of an annulus using Gmsh. This function creates a 2D mesh of an annular segment defined by outer and inner radii, and the extent of the angle. The mesh can be customized with various parameters like cell size, element degree, and verbosity.

Parameters:

radiusOuter (float, optional) – The outer radius of the annular segment. Default is 1.0.
radiusInner (float, optional) – The inner radius of the annular segment. Default is 0.547.
angleExtent (float, optional) – The angular extent of the segment in degrees. Default is 45.
cellSize (float, optional) – The target size for the mesh elements. Default is 0.1.
centre (bool, optional) – If True, the segment will be centered at the origin. Default is False.
degree (int, optional) – The polynomial degree of the finite elements. Default is 1.
qdegree (int, optional) – The quadrature degree for integration. Default is 2.
filename (str, optional) – The name of the file where the mesh will be saved. Default is None.
refinement (optional) – Refinement level for the mesh. Default is None.
gmsh_verbosity (int, optional) – Gmsh output verbosity (0=quiet). Default is 0.
verbose (bool, optional) – If True, print additional information. Default is False.

Returns:

The generated annular segment mesh.

Return type:

Mesh

Examples

>>> mesh = uw.meshing.SegmentofAnnulus(
...     radiusOuter=2.0,
...     radiusInner=1.0,
...     angleExtent=90.0,
...     cellSize=0.05,
...     centre=True,
... )

AnnulusWithSpokes¶

underworld3.meshing.AnnulusWithSpokes(radiusOuter=1.0, radiusInner=0.547, cellSizeOuter=0.1, cellSizeInner=None, centre=False, spokes=3, degree=1, qdegree=2, filename=None, refinement=None, gmsh_verbosity=0, verbose=False)[source]¶

Create an annulus mesh divided into sectors by radial spokes.

Generates a ring-shaped mesh with radial internal boundaries (spokes) that divide the annulus into equal angular sectors. Useful for problems requiring sector-wise analysis or periodic boundary conditions.

Parameters:

radiusOuter (float, default=1.0) – Outer radius of the annulus.
radiusInner (float, default=0.547) – Inner radius of the annulus.
cellSizeOuter (float, default=0.1) – Element size at the outer boundary.
cellSizeInner (float, optional) – Element size at the inner boundary. Defaults to cellSizeOuter.
centre (bool, default=False) – Mark the centre point as a boundary (if radiusInner=0).
spokes (int, default=3) – Number of radial spokes dividing the annulus. Creates spokes equal angular sectors.
degree (int, default=1) – Polynomial degree of finite element basis functions.
qdegree (int, default=2) – Quadrature degree for numerical integration.
filename (str, optional) – Path to save the mesh file.
refinement (int, optional) – Number of uniform refinement levels to apply.
gmsh_verbosity (int, default=0) – Gmsh output verbosity level.
verbose (bool, default=False) – Print diagnostic information.

Returns:

A 2D mesh with boundaries:

Lower: Inner boundary at \(r = r_{inner}\)
LowerPlus: Inner boundary + spokes (for BCs)
Upper: Outer boundary at \(r = r_{outer}\)
UpperPlus: Outer boundary + spokes (for BCs)
Spokes: Radial spoke boundaries

Return type:

Mesh

See also

Annulus: Standard annulus without spokes.

Examples

Create a 6-sector annulus:

>>> import underworld3 as uw
>>> mesh = uw.meshing.AnnulusWithSpokes(
...     radiusOuter=1.0,
...     radiusInner=0.5,
...     spokes=6,
...     cellSizeOuter=0.05
... )

Notes

The LowerPlus and UpperPlus boundaries are composite labels that include both the curved boundaries and the spoke edges, useful for applying no-slip conditions that include radial walls.

AnnulusInternalBoundary¶

underworld3.meshing.AnnulusInternalBoundary(radiusOuter=1.5, radiusInternal=1.0, radiusInner=0.547, cellSize=0.1, cellSize_Outer=None, cellSize_Inner=None, cellSize_Internal=None, centre=False, degree=1, qdegree=2, filename=None, gmsh_verbosity=0, verbose=False)[source]¶

Create an annulus mesh with a circular internal boundary.

Generates an annular mesh with an embedded internal boundary at a specified radius, useful for tracking flux across interfaces or defining different material layers.

Parameters:

radiusOuter (float, default=1.5) – Outer radius of the annulus.
radiusInternal (float, default=1.0) – Radius of the internal boundary surface.
radiusInner (float, default=0.547) – Inner radius of the annulus.
cellSize (float, default=0.1) – Default target mesh element size.
cellSize_Outer (float, optional) – Element size at the outer boundary.
cellSize_Inner (float, optional) – Element size at the inner boundary.
cellSize_Internal (float, optional) – Element size at the internal boundary.
centre (bool, default=False) – Mark the centre point as a boundary (if radiusInner=0).
degree (int, default=1) – Polynomial degree of finite element basis functions.
qdegree (int, default=2) – Quadrature degree for numerical integration.
filename (str, optional) – Path to save the mesh file.
gmsh_verbosity (int, default=0) – Gmsh output verbosity level.
verbose (bool, default=False) – Print diagnostic information.

Returns:

A 2D mesh with boundaries:

Lower: Inner boundary at \(r = r_{inner}\)
Internal: Internal boundary at \(r = r_{internal}\)
Upper: Outer boundary at \(r = r_{outer}\)
Centre: Centre point (if radiusInner=0)

Return type:

Mesh

See also

Annulus: Annulus without internal boundary.
BoxInternalBoundary: Cartesian box with internal boundary.
SphericalShellInternalBoundary: 3D spherical equivalent.

Examples

Create an annulus with an internal tracking surface:

>>> import underworld3 as uw
>>> mesh = uw.meshing.AnnulusInternalBoundary(
...     radiusOuter=1.0,
...     radiusInternal=0.7,
...     radiusInner=0.5,
...     cellSize=0.03
... )

Notes

The internal boundary is useful for calculating radial heat or mass flux across a specified radius (e.g., core-mantle boundary, thermal boundary layer).

DiscInternalBoundaries¶

underworld3.meshing.DiscInternalBoundaries(radiusUpper=1.5, radiusInternal=1.0, radiusLower=0.547, cellSize=0.1, cellSize_Upper=None, cellSize_Lower=None, cellSize_Internal=None, cellSize_Centre=None, degree=1, qdegree=2, filename=None, gmsh_verbosity=0, verbose=False)[source]¶

Create a disc mesh with multiple concentric internal boundaries.

Generates a full disc (extending to centre) with two embedded circular internal boundaries at specified radii. Useful for multi-layer problems or tracking flux across multiple interfaces.

Parameters:

radiusUpper (float, default=1.5) – Outer radius of the disc.
radiusInternal (float, default=1.0) – Radius of the middle internal boundary.
radiusLower (float, default=0.547) – Radius of the inner internal boundary.
cellSize (float, default=0.1) – Default target mesh element size.
cellSize_Upper (float, optional) – Element size at the outer boundary.
cellSize_Lower (float, optional) – Element size at the lower internal boundary.
cellSize_Internal (float, optional) – Element size at the middle internal boundary.
cellSize_Centre (float, optional) – Element size at the disc centre.
degree (int, default=1) – Polynomial degree of finite element basis functions.
qdegree (int, default=2) – Quadrature degree for numerical integration.
filename (str, optional) – Path to save the mesh file.
gmsh_verbosity (int, default=0) – Gmsh output verbosity level.
verbose (bool, default=False) – Print diagnostic information.

Returns:

A 2D mesh with boundaries:

Lower: Inner internal boundary at \(r = r_{lower}\)
Internal: Middle internal boundary at \(r = r_{internal}\)
Upper: Outer boundary at \(r = r_{upper}\)
Centre: Centre point

Return type:

Mesh

See also

AnnulusInternalBoundary: Annulus with one internal boundary.

Examples

Create a disc with two tracking surfaces:

>>> import underworld3 as uw
>>> mesh = uw.meshing.DiscInternalBoundaries(
...     radiusUpper=1.0,
...     radiusInternal=0.7,
...     radiusLower=0.4,
...     cellSize=0.03
... )

Notes

Unlike AnnulusInternalBoundary, this mesh extends to the centre (no inner hole). All internal boundaries are embedded as mesh edges, allowing flux integration across each interface.

Spherical Meshes¶

SphericalShell¶

underworld3.meshing.SphericalShell(radiusOuter=1.0, radiusInner=0.547, cellSize=0.1, degree=1, qdegree=2, filename=None, refinement=None, gmsh_verbosity=0, verbose=False)[source]¶

Create a 3D spherical shell mesh.

Generates a tetrahedral mesh on a full spherical shell (or solid sphere if radiusInner=0) using Gmsh’s OpenCASCADE geometry kernel. Provides \((r, \theta, \phi)\) coordinates for convenient representation of radial problems.

Parameters:

radiusOuter (float, default=1.0) – Outer radius of the spherical shell.
radiusInner (float, default=0.547) – Inner radius of the spherical shell. Set to 0 for a solid sphere extending to the centre.
cellSize (float, default=0.1) – Target mesh element size.
degree (int, default=1) – Polynomial degree of finite element basis functions.
qdegree (int, default=2) – Quadrature degree for numerical integration.
filename (str, optional) – Path to save the mesh file.
refinement (int, optional) – Number of uniform refinement levels to apply. Each level approximately octruples element count.
gmsh_verbosity (int, default=0) – Gmsh output verbosity level.
verbose (bool, default=False) – Print diagnostic information.

Returns:

A 3D mesh with boundaries:

Lower: Inner surface at \(r = r_{inner}\)
Upper: Outer surface at \(r = r_{outer}\)
Centre: Centre point (if radiusInner=0)

The mesh includes a refinement callback that snaps boundary nodes back to true spherical geometry after refinement.

Return type:

Mesh

See also

Annulus: 2D equivalent.
SphericalShellInternalBoundary: Spherical shell with internal surface.
CubedSphere: Structured spherical shell mesh.

Examples

Create a spherical shell for mantle convection:

>>> import underworld3 as uw
>>> mesh = uw.meshing.SphericalShell(
...     radiusOuter=1.0,
...     radiusInner=0.547,
...     cellSize=0.1
... )

Create a solid sphere:

>>> mesh = uw.meshing.SphericalShell(
...     radiusOuter=1.0,
...     radiusInner=0.0,
...     cellSize=0.1
... )

Notes

The inner radius default of 0.547 corresponds approximately to the Earth’s core-mantle boundary radius ratio (3480/6371).

The mesh coordinate system provides unit vectors via mesh.CoordinateSystem:

unit_e_0: radial direction \((r)\)
unit_e_1: colatitude direction \((\theta)\)
unit_e_2: longitude direction \((\phi)\)

For free-slip boundary conditions on vector problems (e.g., Stokes), use a penalty on the normal velocity component:

Gamma_N = mesh.Gamma  # discrete normal, or
Gamma_N = mesh.CoordinateSystem.unit_e_0  # analytic radial
stokes.add_natural_bc(
    penalty * Gamma_N.dot(v_soln.sym) * Gamma_N, "Upper"
)

CubedSphere¶

underworld3.meshing.CubedSphere(radiusOuter=1.0, radiusInner=0.547, numElements=5, degree=1, qdegree=2, simplex=False, filename=None, refinement=None, gmsh_verbosity=0, verbose=False)[source]¶

Create a cubed-sphere mesh with structured elements.

Generates a spherical shell mesh using the cubed-sphere projection, which maps six cube faces onto the sphere. Produces hexahedral elements (or tetrahedra if simplex=True) with uniform resolution.

Parameters:

radiusOuter (float, default=1.0) – Outer radius of the spherical shell.
radiusInner (float, default=0.547) – Inner radius of the spherical shell.
numElements (int, default=5) – Number of elements along each edge of each cube face. Total elements approximately 6 * numElements^2 per radial layer.
degree (int, default=1) – Polynomial degree of finite element basis functions.
qdegree (int, default=2) – Quadrature degree for numerical integration.
simplex (bool, default=False) – If False, use hexahedral elements. If True, subdivide into tetrahedra (simplex elements).
filename (str, optional) – Path to save the mesh file.
refinement (int, optional) – Number of uniform refinement levels to apply.
gmsh_verbosity (int, default=0) – Gmsh output verbosity level.
verbose (bool, default=False) – Print diagnostic information.

Returns:

A 3D mesh with boundaries:

Lower: Inner surface at \(r = r_{inner}\)
Upper: Outer surface at \(r = r_{outer}\)

Return type:

Mesh

See also

SphericalShell: Unstructured tetrahedral spherical mesh.
SegmentofSphere: Partial spherical shell.

Examples

Create a cubed-sphere mesh with 10 elements per edge:

>>> import underworld3 as uw
>>> mesh = uw.meshing.CubedSphere(
...     radiusOuter=1.0,
...     radiusInner=0.55,
...     numElements=10
... )

Notes

The cubed-sphere projection avoids polar singularities present in latitude-longitude grids, providing quasi-uniform element sizes across the sphere. This is particularly advantageous for global geodynamics simulations.

The mesh coordinate system provides unit vectors via mesh.CoordinateSystem:

unit_e_0: radial direction \((r)\)
unit_e_1: colatitude direction \((\theta)\)
unit_e_2: longitude direction \((\phi)\)

SegmentofSphere¶

underworld3.meshing.SegmentofSphere(radiusOuter=1.0, radiusInner=0.547, longitudeExtent=90.0, latitudeExtent=90.0, cellSize=0.1, degree=1, qdegree=2, filename=None, refinement=None, gmsh_verbosity=0, verbose=False, centroid=(0.0, 0.0, 0.0))[source]¶

Generates a segment of a sphere using Gmsh. This function creates a 3D mesh of a spherical segment defined by outer and inner radii, and the extent in longitude and latitude. The mesh can be customized in terms of size, polynomial degree, and verbosity.

Parameters:

radiusOuter (float, optional) – The outer radius of the spherical segment. Default is 1.0.
radiusInner (float, optional) – The inner radius of the spherical segment. Default is 0.547.
longitudeExtent (float, optional) – The angular extent in longitude (degrees). Default is 90.0.
latitudeExtent (float, optional) – The angular extent in latitude (degrees). Default is 90.0.
cellSize (float, optional) – The target size for the mesh elements. Default is 0.1.
degree (int, optional) – The polynomial degree of the finite elements. Default is 1.
qdegree (int, optional) – The quadrature degree for integration. Default is 2.
filename (str, optional) – The name of the file where the mesh will be saved. Default is None.
refinement (optional) – Refinement level for the mesh. Default is None.
gmsh_verbosity (int, optional) – Gmsh output verbosity (0=quiet). Default is 0.
verbose (bool, optional) – If True, print additional information. Default is False.
centroid (tuple of float, optional) – The coordinates of the centroid of the sphere segment. Default is (0.0, 0.0, 0.0).

Returns:

The generated spherical segment mesh.

Return type:

Mesh

Examples

>>> mesh = uw.meshing.SegmentofSphere(
...     radiusOuter=2.0,
...     radiusInner=1.0,
...     longitudeExtent=120.0,
...     latitudeExtent=60.0,
...     cellSize=0.05,
... )

SphericalShellInternalBoundary¶

underworld3.meshing.SphericalShellInternalBoundary(radiusOuter=1.0, radiusInternal=0.8, radiusInner=0.547, cellSize=0.1, degree=1, qdegree=2, filename=None, refinement=None, gmsh_verbosity=0, verbose=False)[source]¶

Generates a spherical shell with an internal boundary using Gmsh. The function creates a 3D mesh of a spherical shell defined by outer, internal, and inner radii. Mesh size, polynomial degree, and Gmsh verbosity can be customized.

Parameters:

radiusOuter (float, optional) – The outer radius of the spherical shell. Default is 1.0.
radiusInternal (float, optional) – The radius of the internal boundary. Default is 0.8.
radiusInner (float, optional) – The inner radius of the spherical shell. Default is 0.547.
cellSize (float, optional) – The target size for the mesh elements. Default is 0.1.
degree (int, optional) – The polynomial degree of the finite elements. Default is 1.
qdegree (int, optional) – The quadrature degree for integration. Default is 2.
filename (str, optional) – The name of the file where the mesh will be saved. Default is None.
refinement (optional) – Refinement level for the mesh. Default is None.
gmsh_verbosity (int, optional) – Gmsh output verbosity (0=quiet). Default is 0.
verbose (bool, optional) – If True, print additional information. Default is False.

Returns:

The generated spherical shell mesh with internal boundary.

Return type:

Mesh

Examples

>>> mesh = uw.meshing.SphericalShellInternalBoundary(
...     radiusOuter=2.0,
...     radiusInternal=1.5,
...     radiusInner=1.0,
...     cellSize=0.05,
... )

Geographic Meshes¶

Regional meshes defined using longitude/latitude coordinates.

RegionalSphericalBox¶

underworld3.meshing.RegionalSphericalBox(radiusOuter=1.0, radiusInner=0.547, SWcorner=[-45, -45], NEcorner=[45, 45], numElementsLon=5, numElementsLat=5, numElementsDepth=5, degree=1, qdegree=2, simplex=False, filename=None, refinement=None, coarsening=None, gmsh_verbosity=0, verbose=False)[source]¶

Create a regional spherical box mesh (cubed-sphere section).

Generates a 3D structured mesh for a regional section of a spherical shell, using a cubed-sphere projection. The domain is defined by corner coordinates in degrees (longitude, latitude) and radial bounds.

Parameters:

radiusOuter (float, default=1.0) – Outer radius of the spherical shell.
radiusInner (float, default=0.547) – Inner radius of the spherical shell.
SWcorner (list of float, default=[-45, -45]) – Southwest corner as [longitude, latitude] in degrees.
NEcorner (list of float, default=[+45, +45]) – Northeast corner as [longitude, latitude] in degrees.
numElementsLon (int, default=5) – Number of elements in the longitude direction.
numElementsLat (int, default=5) – Number of elements in the latitude direction.
numElementsDepth (int, default=5) – Number of elements in the radial (depth) direction.
degree (int, default=1) – Polynomial degree of finite elements.
qdegree (int, default=2) – Quadrature degree for numerical integration.
simplex (bool, default=False) – If True, use tetrahedral elements; if False, use hexahedral.
filename (str, optional) – Path to save the mesh file.
refinement (int, optional) – Number of uniform refinement levels to apply.
coarsening (int, optional) – Number of coarsening levels to apply.
gmsh_verbosity (int, default=0) – Gmsh output verbosity level.
verbose (bool, default=False) – Print diagnostic information.

Returns:

A 3D mesh with boundaries:

Lower: Inner surface at \(r = r_{inner}\)
Upper: Outer surface at \(r = r_{outer}\)
North: Northern boundary at \(\phi = \phi_{max}\)
South: Southern boundary at \(\phi = \phi_{min}\)
East: Eastern boundary at \(\lambda = \lambda_{max}\)
West: Western boundary at \(\lambda = \lambda_{min}\)

The mesh uses a SPHERICAL coordinate system and includes a refinement callback that snaps boundary nodes to true spherical geometry.

Return type:

Mesh

See also

CubedSphere: Full cubed-sphere mesh.
RegionalGeographicBox: Geographic mesh with ellipsoidal geometry.
SphericalShell: Unstructured spherical shell.

Examples

Create a regional mesh for the Australian region:

>>> import underworld3 as uw
>>> mesh = uw.meshing.RegionalSphericalBox(
...     radiusOuter=1.0,
...     radiusInner=0.9,
...     SWcorner=[110, -45],
...     NEcorner=[155, -10],
...     numElementsLon=10,
...     numElementsLat=8,
...     numElementsDepth=5
... )

Notes

This mesh uses a cubed-sphere projection, which provides more uniform element sizes than a latitude-longitude grid. The structured mesh is suitable for regional mantle convection models where boundary-aligned elements are beneficial.

The coordinate system provides unit vectors via mesh.CoordinateSystem:

unit_e_0: radial direction \((r)\)
unit_e_1: colatitude direction \((\theta)\)
unit_e_2: longitude direction \((\phi)\)

RegionalGeographicBox¶

underworld3.meshing.RegionalGeographicBox(lon_range=(135.0, 140.0), lat_range=(-35.0, -30.0), depth_range=(0.0, 400.0), ellipsoid='WGS84', numElements=(10, 10, 10), degree=1, qdegree=2, simplex=True, filename=None, refinement=None, coarsening=None, gmsh_verbosity=0, verbose=False)[source]¶

Create a regional geographic mesh with ellipsoidal geometry.

This function creates a structured 3D mesh in geographic coordinates (longitude, latitude, depth) on an ellipsoidal planet. The mesh uses geodetic latitude (perpendicular to ellipsoid surface) and measures depth below the reference ellipsoid surface.

Parameters:

lon_range (tuple of float, optional) – Longitude range in degrees East (lon_min, lon_max). Default: (135.0, 140.0) for southeastern Australia.
lat_range (tuple of float, optional) – Latitude range in degrees North (lat_min, lat_max). Geodetic latitude (perpendicular to ellipsoid). Default: (-35.0, -30.0) for southeastern Australia.
depth_range (tuple of float, optional) – Depth range in km below ellipsoid surface (depth_min, depth_max). Positive downward. depth_min=0 means surface. Default: (0.0, 400.0) for 0-400 km depth.
ellipsoid (str, tuple, or bool, optional) – Ellipsoid specification: - str: Name from ELLIPSOIDS dict (‘WGS84’, ‘Mars’, ‘Moon’, ‘Venus’, ‘sphere’) - tuple: (semi_major_axis_km, semi_minor_axis_km) for custom ellipsoid - True: Use WGS84 (default) - False or ‘sphere’: Use perfect sphere with Earth mean radius Default: ‘WGS84’
numElements (tuple of int, optional) – Number of elements in (lon, lat, depth) directions. Default: (10, 10, 10)
degree (int, optional) – Polynomial degree for finite elements (1=linear, 2=quadratic). Default: 1
qdegree (int, optional) – Quadrature degree for numerical integration. Default: 2
simplex (bool, optional) – If True, use tetrahedral elements. If False, use hexahedral elements. Default: True
filename (str, optional) – Path to save generated mesh file. If None, uses automatic naming. Default: None
refinement (int, optional) – Number of uniform refinement steps to apply. Default: None
coarsening (int, optional) – Number of coarsening steps to apply. Default: None
gmsh_verbosity (int, optional) – Gmsh output verbosity level (0=quiet, 5=very verbose). Default: 0
verbose (bool, optional) – If True, print mesh generation details. Default: False

Returns:

Underworld3 mesh object with GEOGRAPHIC coordinate system. Access geographic coordinates via mesh.geo: - mesh.geo.lon, mesh.geo.lat, mesh.geo.depth (data arrays) - mesh.geo[:] for symbolic coordinates (λ_lon, λ_lat, λ_d) - mesh.geo.unit_east, mesh.geo.unit_north, mesh.geo.unit_down (basis vectors)

Return type:

Mesh

Examples

Create mesh for southeastern Australia, 0-400 km depth:

mesh = uw.meshing.RegionalGeographicBox(
    lon_range=(135, 140),
    lat_range=(-35, -30),
    depth_range=(0, 400),
    ellipsoid='WGS84',
    numElements=(20, 20, 10),
)

Access geographic coordinates:

lon = mesh.geo.lon
lat = mesh.geo.lat
depth = mesh.geo.depth

Notes

Uses geodetic latitude (GPS/map standard), not geocentric latitude. Depth is measured from reference ellipsoid surface, not from center. mesh.R provides spherical coordinates for backward compatibility. mesh.geo provides geographic coordinates with ellipsoid geometry. Right-handed coordinate system: WE x SN = down.

Segmented Meshes¶

Meshes with internal segment boundaries for multi-region problems.

SegmentedSphericalShell¶

underworld3.meshing.SegmentedSphericalShell(radiusOuter=1.0, radiusInner=0.547, cellSize=0.1, numSegments=6, degree=1, qdegree=2, filename=None, refinement=None, coordinatesNative=False, gmsh_verbosity=0, verbose=False)[source]¶

Create a 3D spherical shell mesh using orange-peel segmentation.

Generates a spherical shell mesh by dividing the sphere into longitudinal wedge segments, similar to an orange peel. Each segment is a 3D wedge extending from the inner to outer radius.

Parameters:

radiusOuter (float, default=1.0) – Outer radius of the shell.
radiusInner (float, default=0.547) – Inner radius of the shell.
cellSize (float, default=0.1) – Target mesh element size.
numSegments (int, default=6) – Number of longitudinal segments.
degree (int, default=1) – Polynomial degree of finite elements.
qdegree (int, default=2) – Quadrature degree for numerical integration.
filename (str, optional) – Path to save the mesh file.
refinement (int, optional) – Number of uniform refinement levels to apply.
coordinatesNative (bool, default=False) – If True, use native \((r, \theta, \phi)\) coordinates with periodicity in \(\phi\). If False, use Cartesian coordinates with spherical coordinate system overlay.
gmsh_verbosity (int, default=0) – Gmsh output verbosity level.
verbose (bool, default=False) – Print diagnostic information.

Returns:

A 3D mesh with boundaries:

Lower: Inner surface at \(r = r_{inner}\)
LowerPlus: Inner surface + slice boundaries
Upper: Outer surface at \(r = r_{outer}\)
UpperPlus: Outer surface + slice boundaries
Slices: Radial slice boundaries between segments

The mesh provides boundary normals for free-slip conditions on the radial boundaries.

Return type:

Mesh

See also

SegmentedSphericalBall: Solid ball version (no inner boundary).
SphericalShell: Standard unstructured spherical shell.
CubedSphere: Alternative structured spherical mesh.

Warning

When coordinatesNative=True, this mesh uses PETSc periodicity features that may be unreliable. Consider using coordinatesNative=False or alternative mesh types for production work.

Notes

The LowerPlus and UpperPlus boundaries combine the radial surfaces with the slice boundaries, useful for applying no-slip conditions that include the radial walls.

The coordinate system provides unit vectors via mesh.CoordinateSystem:

unit_e_0: radial direction \((r)\)
unit_e_1: colatitude direction \((\theta)\)
unit_e_2: longitude direction \((\phi)\)

SegmentedSphericalSurface2D¶

underworld3.meshing.SegmentedSphericalSurface2D(radius=1.0, cellSize=0.05, numSegments=6, degree=1, qdegree=2, filename=None, gmsh_verbosity=0, verbose=False)[source]¶

Create a 2D spherical surface mesh using orange-peel segmentation.

Generates a mesh on the surface of a sphere by dividing it into longitudinal segments (like an orange peel), with coordinates converted to \((\lambda, \phi)\) (longitude, latitude).

Parameters:

radius (float, default=1.0) – Radius of the sphere.
cellSize (float, default=0.05) – Target mesh element size.
numSegments (int, default=6) – Number of longitudinal segments (orange-peel divisions).
degree (int, default=1) – Polynomial degree of finite elements.
qdegree (int, default=2) – Quadrature degree for numerical integration.
filename (str, optional) – Path to save the mesh file.
gmsh_verbosity (int, default=0) – Gmsh output verbosity level.
verbose (bool, default=False) – Print diagnostic information.

Returns:

A 2D surface mesh with point boundaries:

NPole: North pole point
SPole: South pole point
Poles: Both poles combined

Uses SPHERE_SURFACE_NATIVE coordinate system with \((\lambda, \phi)\) coordinates.

Return type:

Mesh

See also

SegmentedSphericalShell: 3D segmented spherical shell.
SphericalShell: Standard unstructured spherical shell.

Warning

This mesh uses PETSc periodicity features that may be unreliable. Consider using alternative mesh types for production work.

Notes

The mesh is constructed by creating curved triangular segments from pole to pole, then converting the Cartesian coordinates to longitude- latitude representation. Periodicity is set up to handle the \(\pm\pi\) longitude wrap-around.

SegmentedSphericalBall¶

underworld3.meshing.SegmentedSphericalBall(radius=1.0, cellSize=0.1, numSegments=6, degree=1, qdegree=2, filename=None, refinement=None, coordinatesNative=False, verbosity=0, gmsh_verbosity=0, verbose=False)[source]¶

Create a 3D solid spherical ball mesh using orange-peel segmentation.

Generates a solid sphere mesh (no inner cavity) by dividing it into longitudinal wedge segments extending from the centre to the surface.

Parameters:

radius (float, default=1.0) – Radius of the sphere.
cellSize (float, default=0.1) – Target mesh element size.
numSegments (int, default=6) – Number of longitudinal segments.
degree (int, default=1) – Polynomial degree of finite elements.
qdegree (int, default=2) – Quadrature degree for numerical integration.
filename (str, optional) – Path to save the mesh file.
refinement (int, optional) – Number of uniform refinement levels to apply.
coordinatesNative (bool, default=False) – If True, use native \((r, \theta, \phi)\) coordinates with periodicity in \(\phi\). If False, use Cartesian coordinates with spherical coordinate system overlay.
verbosity (int, default=0) – Deprecated. Use verbose instead.
gmsh_verbosity (int, default=0) – Gmsh output verbosity level.
verbose (bool, default=False) – Print diagnostic information.

Returns:

A 3D mesh with boundaries:

Upper: Outer surface at \(r = radius\)
UpperPlus: Outer surface + slice boundaries
Centre: Central point
Slices: Radial slice boundaries between segments

The mesh provides boundary normals for free-slip conditions on the outer surface.

Return type:

Mesh

See also

SegmentedSphericalShell: Shell version with inner boundary.
SphericalShell: Standard unstructured spherical shell.

Warning

When coordinatesNative=True, this mesh uses PETSc periodicity features that may be unreliable. Consider using coordinatesNative=False or alternative mesh types for production work.

Notes

This mesh is useful for modelling solid bodies (e.g., planetary cores, asteroids) where elements must extend to the centre. The segmented approach avoids the pole singularity issues of latitude-longitude grids.

The coordinate system provides unit vectors via mesh.CoordinateSystem:

unit_e_0: radial direction \((r)\)
unit_e_1: colatitude direction \((\theta)\)
unit_e_2: longitude direction \((\phi)\)

Faults and Internal Boundaries¶

Tools for creating meshes with fault surfaces and internal boundaries.

FaultSurface¶

class underworld3.meshing.FaultSurface[source]¶

Bases: object

A triangulated fault surface with orientation data.

Represents a single fault segment as a 2D surface embedded in 3D space. Can be created from: - A point cloud (requires triangulation via pyvista) - A VTK file (pre-triangulated surface)

Attributes:

name: Identifier for this fault segment points: (N, 3) array of surface points triangles: (M, 3) array of triangle vertex indices (after triangulation) normals: (M, 3) array of face normals (after triangulation) pv_mesh: PyVista PolyData object (if pyvista available) is_triangulated: Whether the surface has been triangulated

Example:

>>> # Create from points and triangulate
>>> points = np.array([[0, 0, 0], [1, 0, 0], [0.5, 1, 0], [0.5, 0.5, 1]])
>>> fault = uw.meshing.FaultSurface("fault1", points)
>>> fault.triangulate()
>>> fault.to_vtk("fault1.vtk")
>>>
>>> # Load from VTK
>>> fault2 = uw.meshing.FaultSurface.from_vtk("fault1.vtk")

__init__(name, points=None)[source]¶

Create a fault surface.

Args:: name: Identifier for this fault segment. points: (N, 3) array of 3D points, or None for an empty fault.

Parameters:

name (str)
points (ndarray)

property points: ndarray | None¶: (N, 3) array of surface points.

property triangles: ndarray | None¶: (M, 3) array of triangle vertex indices.

property normals: ndarray | None¶: (M, 3) array of face normals.

property pv_mesh¶: PyVista PolyData mesh (None if not triangulated or pyvista unavailable).

property is_triangulated: bool¶: Whether the surface has been triangulated.

property n_points: int¶: Number of points in the surface.

property n_triangles: int¶: Number of triangles in the surface.

classmethod from_vtk(filename, name=None)[source]¶

Load fault surface from VTK file.

Args:: filename: Path to VTK file (.vtk or .vtp) name: Name for the fault. If None, uses filename stem.
Returns:: FaultSurface: Loaded fault surface with triangulation and normals
Raises:: FileNotFoundError: If file doesn’t exist ImportError: If pyvista not available

Parameters:

filename (str)
name (str)

Return type:

FaultSurface

triangulate(offset=0.01)[source]¶

Triangulate point cloud using pyvista delaunay_2d.

This creates a triangulated surface from the point cloud by projecting points onto a best-fit plane, performing 2D Delaunay triangulation, and mapping back to 3D.

Args:

offset: Height offset for delaunay_2d (controls curvature tolerance).: Larger values allow more curved surfaces.

Raises:

ImportError: If pyvista not available ValueError: If points too sparse for triangulation (< 3 points) RuntimeError: If triangulation fails

Parameters:: offset (float)
Return type:: None

compute_normals(consistent_normals=True)[source]¶

Recompute face normals for triangulated surface.

Args:: consistent_normals: If True, attempt to make normals consistently oriented

Parameters:: consistent_normals (bool)
Return type:: None

flip_normals()[source]¶

Flip the direction of all face normals.

Return type:: None

to_vtk(filename)[source]¶

Export triangulated surface to VTK file.

Args:: filename: Output path (.vtk or .vtp)
Raises:: RuntimeError: If surface not triangulated ImportError: If pyvista not available

Parameters:: filename (str)
Return type:: None

build_kdtree()[source]¶

Build KDTree for nearest-neighbor queries on face centers.

Returns:: KDTree built from triangle centroids
Raises:: RuntimeError: If surface not triangulated

Return type:: KDTree

property face_centers: ndarray¶: (M, 3) array of triangle centroids.

FaultCollection¶

class underworld3.meshing.FaultCollection[source]¶

Bases: object

Collection of fault surfaces for mesh integration.

Manages multiple FaultSurface objects and provides methods to: - Compute minimum distance from mesh points to any fault - Transfer fault normals to mesh variables via nearest-neighbor - Create rheology functions for fault-weakened zones

Example:

>>> faults = uw.meshing.FaultCollection()
>>> faults.add_from_vtk("fault1.vtk")
>>> faults.add_from_vtk("fault2.vtk")
>>>
>>> # Compute distance field
>>> fault_distance = faults.compute_distance_field(mesh)
>>>
>>> # Transfer normals
>>> fault_normals = faults.transfer_normals(mesh)
>>>
>>> # Create weakness function for rheology
>>> eta_weak = faults.create_weakness_function(
...     fault_distance,
...     fault_width=mesh.get_min_radius() * 5,
...     eta_weak=0.01,
... )

__init__()[source]¶: Create an empty fault collection.

add(fault)[source]¶

Add a fault surface to the collection.

Args:: fault: FaultSurface to add
Raises:: ValueError: If fault with same name already exists

Parameters:: fault (FaultSurface)
Return type:: None

add_from_vtk(filename, name=None)[source]¶

Load and add a fault from VTK file.

Args:: filename: Path to VTK file name: Name for the fault. If None, uses filename stem.
Returns:: The loaded FaultSurface

Parameters:

filename (str)
name (str)

Return type:

FaultSurface

remove(name)[source]¶

Remove and return a fault from the collection.

Args:: name: Name of fault to remove
Returns:: The removed FaultSurface
Raises:: KeyError: If fault not found

Parameters:: name (str)
Return type:: FaultSurface

property names: List[str]¶: List of fault names.

compute_distance_field(mesh, distance_var=None, variable_name='fault_distance')[source]¶

Compute minimum distance from mesh points to any fault surface.

Uses pyvista’s compute_implicit_distance for accurate signed distance computation. The returned field contains the absolute distance to the nearest fault surface at each mesh point.

Parameters:

mesh (Mesh) – The mesh to compute distances on.
distance_var (MeshVariable, optional) – Existing MeshVariable to store results. If None, creates a new variable.
variable_name (str, optional) – Name for new variable if distance_var is None.

Returns:

Scalar variable with distance values (1 component).

Return type:

MeshVariable

Raises:

ValueError – If collection is empty or no faults are triangulated.
ImportError – If pyvista is not available.

transfer_normals(mesh, coords=None, normal_var=None, variable_name='fault_normals')[source]¶

Transfer fault normals to mesh points via nearest-neighbor lookup.

For each mesh point, finds the closest fault face (from any fault in the collection) and copies that face’s normal vector.

Parameters:

mesh (Mesh) – The mesh to transfer normals to.
coords (ndarray, optional) – Coordinates to query. If None, uses mesh.X.coords.
normal_var (MeshVariable, optional) – Existing MeshVariable to store results. If None, creates a new variable.
variable_name (str, optional) – Name for new variable if normal_var is None.

Returns:

Variable with normal vectors (3 components).

Return type:

MeshVariable

Raises:

ValueError – If collection is empty or no faults are triangulated.

create_weakness_function(distance_var, fault_width, eta_weak=0.01, eta_background=1.0)[source]¶

Create Piecewise viscosity function for fault weakness.

Creates a sympy Piecewise expression that gives: - eta_weak when distance < fault_width - eta_background otherwise

This can be used directly with TransverseIsotropicFlowModel.Parameters.eta_1 for creating anisotropic weakness along fault zones.

Args:

distance_var: MeshVariable containing fault distances fault_width: Width of the weak zone around faults eta_weak: Viscosity within fault zone (default 0.01) eta_background: Viscosity outside fault zone (default 1.0)

Returns:

sympy.Piecewise expression for use in constitutive models

Example:

>>> eta_1 = faults.create_weakness_function(
...     fault_distance,
...     fault_width=mesh.get_min_radius() * 5,
...     eta_weak=0.01,
... )
>>> stokes.constitutive_model.Parameters.eta_1 = eta_1

Parameters:

distance_var (MeshVariable)
fault_width (float)
eta_weak (float)
eta_background (float)

Return type:

sympy.Expr

Surfaces¶

Tools for working with surface meshes and embedded boundaries.

Surface¶

class underworld3.meshing.Surface[source]¶

Bases: object

A discretized embedded surface with variable storage.

Represents a surface embedded in a computational mesh using pyvista PolyData for storage. In 3D, surfaces are triangulated meshes; in 2D, they are polylines. Supports per-vertex variables with symbolic access for use in Underworld expressions.

The surface uses lazy evaluation with stale flags: - Discretization is computed when first accessed - Distance field is computed when first accessed - Proxy MeshVariables are updated when .sym is accessed

Attributes:

name: Identifier for this surface mesh: Associated computational mesh (for proxy MeshVariables) vertices: (N, 3) array of vertex positions (via pyvista) normals: (N, 3) array of vertex normals (via pyvista) n_vertices: Number of vertices n_triangles: Number of triangles (3D only)

Example:

>>> # Create from points and discretize
>>> surface = uw.meshing.Surface("fault", mesh, points)
>>> surface.discretize()
>>>
>>> # Add a variable
>>> friction = surface.add_variable("friction")
>>> friction.data[:] = 0.6
>>>
>>> # Use distance field in influence function
>>> eta_weak = surface.influence_function(
...     width=0.05,
...     value_near=0.01,
...     value_far=1.0,
...     profile="gaussian",
... )

__init__(name, mesh=None, control_points=None, symbol=None)[source]¶

Create a surface.

Parameters:

name (str) – Identifier for this surface.
mesh (Mesh, optional) – Computational mesh (required for .sym access and distance field).
control_points (ndarray, optional) – (N, 3) array of 3D points defining the surface. If None, the surface is empty and must be loaded or have points set later.
symbol (str, optional) – Short LaTeX-friendly symbol for math display (e.g., "F"). If None, defaults to first letter of name capitalized.

property dim: int¶

Spatial dimension (2 or 3).

Detected from mesh.dim if available, otherwise from control points shape.

property is_2d: bool¶: True if this is a 2D surface (1D curve in 2D space).

property symbol: str¶

Short LaTeX-friendly symbol for math expressions.

Used in distance field expressions like \(d_F\) instead of the full variable name {surf_fault_distance}.

property control_points: ndarray | None¶

(N, 3) array of control points defining the surface.

Returns coordinates in physical (dimensional) units when the units system is active, matching the mesh.X.coords convention. Internally, points are stored in model (non-dimensional) space.

set_control_points(points)[source]¶

Set control points and mark discretization as stale.

Coordinates are stored internally in model (non-dimensional) space, matching the mesh’s internal coordinate representation (mesh._coords).

Args:

points: (N, 2) or (N, 3) array of points in one of these forms:

Raw numpy array: Assumed to be in model coordinates (same space as mesh._coords). When units are active, this means nondimensional coordinates.

Pint Quantity / UnitAwareArray: Automatically converted to model coordinates via the scaling system.

For 2D points, a z=0 column is appended automatically.

Parameters:: points (ndarray)
Return type:: None

property vertices: ndarray | None¶

(N, 3) array of vertex positions.

Returns coordinates in physical (dimensional) units when the units system is active, matching the mesh.X.coords convention. Internally, vertices are stored in model (non-dimensional) space.

property n_vertices: int¶: Number of vertices in the discretized surface.

property n_triangles: int¶: Number of triangles in the surface.

property normals: ndarray | None¶: (N, 3) array of vertex normals (point normals from pyvista).

property triangles: ndarray | None¶: (M, 3) array of triangle vertex indices.

property face_centers: ndarray | None¶: (M, 3) array of triangle centroids.

property face_normals: ndarray | None¶: (M, 3) array of face normals (cell normals from pyvista).

property pv_mesh¶

PyVista PolyData mesh in the same coordinate space as mesh.X.coords.

Returns a shallow copy whose points have been dimensionalised so that they overlay correctly with vis.mesh_to_pv_mesh(mesh). The internal _pv_mesh stays in nondimensional model space so that distance calculations remain consistent with the solver.

Returns None if the surface has not been discretized.

property is_discretized: bool¶: Whether the surface has been discretized.

discretize(offset=0.01, n_segments=None)[source]¶

Discretize control points into a surface mesh.

For 3D surfaces: Uses pyvista delaunay_2d to create a triangulated mesh. For 2D surfaces: Uses scipy to fit a spline and create a polyline.

Args:: offset: (3D only) Height offset for delaunay_2d (controls curvature tolerance). n_segments: (2D only) Number of line segments, or None for auto.
Raises:: ImportError: If pyvista not available ValueError: If points too sparse for discretization RuntimeError: If discretization fails

Parameters:

offset (float)
n_segments (int)

Return type:

None

deform_vertices(displacement)[source]¶

Deform surface vertices in-place (no re-discretization).

This modifies vertex positions while keeping topology fixed. Normals are automatically recomputed.

Args:: displacement: (n_vertices, 3) array of displacements to add
Note:: This does NOT update control points. If you want topology changes, use set_control_points() instead and call discretize().

Parameters:: displacement (ndarray)
Return type:: None

property distance: EnhancedMeshVariable¶

Signed distance from mesh nodes to surface (lazily computed).

The signed distance is positive on one side of the surface and negative on the other. Use sympy.Abs(surface.distance.sym[0]) for unsigned distance, or use influence_function() which does this automatically.

Returns:

MeshVariable with signed distance values at each mesh node. Access .sym[0] for use in expressions.

Raises:

RuntimeError: If mesh not set or surface not discretized

Example:

>>> # Use signed distance for different properties on each side
>>> d = surface.distance.sym[0]
>>> prop = sympy.Piecewise((upper_value, d > 0), (lower_value, True))
>>>
>>> # Use absolute distance for symmetric influence
>>> mask = sympy.Piecewise((1, sympy.Abs(d) < width), (0, True))

influence_function(width, value_near=1.0, value_far=0.0, profile='step')[source]¶

Create level-set-like influence function based on distance.

Creates a sympy expression that varies from value_near (at the surface) to value_far (far from the surface) based on the chosen profile.

Uses the absolute value of the signed distance field, so the influence is symmetric on both sides of the surface. For asymmetric behavior, access the signed distance directly via surface.distance.sym[0].

Parameters:

width (float) – Characteristic width of the transition zone.
value_near (float or sympy.Expr, optional) – Value at/near the surface. Default is 1.0.
value_far (float or sympy.Expr, optional) – Value far from the surface. Default is 0.0.
profile (str, optional) – Transition profile type. One of "step" (sharp transition), "linear" (linear ramp), "gaussian" (smooth decay), or "smoothstep" (C1-continuous Hermite). Default is "step".

Returns:

Expression that can be used in Underworld equations.

Return type:

sympy.Expr

Examples

Step function for fault zone viscosity:

>>> eta = surface.influence_function(
    ...     width=0.05,
    ...     value_near=0.01,
    ...     value_far=1.0,
    ...     profile="step",
    ... )
    >>>
    >>> # Gaussian decay for smooth transitions
    >>> eta = surface.influence_function(
    ...     width=0.1,
    ...     value_near=friction.sym,  # Variable on surface
    ...     value_far=1.0,
    ...     profile="gaussian",
    ... )

add_variable(name, size=1, proxy_degree=1, units=None, mask_width=None, mask_profile='gaussian')[source]¶

Add a variable on surface vertices.

Creates a SurfaceVariable stored in pyvista point_data with symbolic access via .sym for use in expressions.

Args:

name: Variable name size: Number of components (1 for scalar, 3 for vector) proxy_degree: Degree of proxy MeshVariable for .sym access units: Optional units for this variable (e.g., “Pa”, “m/s”) mask_width: Width for distance-based mask (enables .mask property) mask_profile: Profile for mask function (“step”, “linear”, “gaussian”, “smoothstep”)

Returns:

SurfaceVariable that can be modified via .data and used via .sym

Example:

>>> # Variable with units and mask
>>> friction = surface.add_variable("friction", size=1, mask_width=0.1)
>>> friction.data[:] = 0.6
>>>
>>> # Use in expressions with explicit masking
>>> eta = friction.sym[0] * friction.mask

Parameters:

name (str)
size (int)
proxy_degree (int)
units (str | None)
mask_width (float | None)
mask_profile (str)

Return type:

SurfaceVariable

get_variable(name)[source]¶

Get an existing variable by name.

Args:: name: Variable name
Returns:: SurfaceVariable
Raises:: KeyError: If variable doesn’t exist

Parameters:: name (str)
Return type:: SurfaceVariable

property variables: Dict[str, SurfaceVariable]¶: Dictionary of all variables on this surface.

save(filename)[source]¶

Save surface with all variables to VTK file.

All SurfaceVariable data is automatically included in the VTK file as point_data arrays.

Args:: filename: Output path (.vtk or .vtp)

Parameters:: filename (str)
Return type:: None

refinement_metric(h_near, h_far, width=None, profile='linear', name=None)[source]¶

Create a metric field for mesh adaptation based on distance from this surface.

Returns a MeshVariable containing refinement metric values that can be passed directly to mesh.adapt(). Higher metric values produce finer mesh (smaller elements).

Parameters:

h_near (float or quantity) – Target edge length near the surface (smaller = finer mesh). Accepts a plain float (in nondimensional mesh coordinates) or a uw.quantity (e.g., uw.quantity(3, "km")) which is automatically nondimensionalised.
h_far (float or quantity) – Target edge length far from the surface (larger = coarser mesh). Same unit handling as h_near.
width (float or quantity, optional) – Distance over which to transition from h_near to h_far. If None, defaults to 2 * h_far. Same unit handling as h_near.
profile (str, optional) – Transition profile: “linear”, “smoothstep”, or “gaussian”. Default is “linear”.
name (str, optional) – Name for the metric MeshVariable. Defaults to “{surface_name}_metric”.

Returns:

Scalar MeshVariable containing refinement metric values.

Return type:

MeshVariable

Notes

Metric Tensor Mathematics

For isotropic mesh adaptation, MMG/PETSc uses a metric tensor:

\[M = h^{-2} \cdot I\]

where \(h\) is the target edge length and \(I\) is the identity matrix. This relationship is dimension-independent - the same formula applies in 2D and 3D because the metric defines edge lengths, not areas or volumes.

The adaptation algorithm seeks to make all edges have unit length in the metric space (i.e., \(\mathbf{e}^T M \mathbf{e} = 1\) for edge vector \(\mathbf{e}\)). Higher metric values produce smaller elements.

Refinement Ratio

The refinement ratio is h_far / h_near. For example, if h_near=0.01 and h_far=0.1, the mesh will be ~10× finer near the surface.

Element Count Control

To maintain approximately the same total element count while refining near the surface, the far-field should use similar h to the original mesh’s cell size. The refined region is small, so coarsening the far-field slightly can compensate.

References

Examples

>>> fault = uw.meshing.Surface("fault", mesh, fault_points)
>>> fault.discretize()
>>>
>>> # With plain floats (nondimensional coordinates)
>>> metric = fault.refinement_metric(h_near=0.005, h_far=0.05)
>>> mesh.adapt(metric)
>>>
>>> # With quantities (automatic nondimensionalisation)
>>> metric = fault.refinement_metric(
...     h_near=uw.quantity(3, "km"),
...     h_far=uw.quantity(30, "km"),
...     width=uw.quantity(10, "km"),
... )
>>> mesh.adapt(metric)

classmethod from_trace(name, mesh, trace_points, depth_range, n_depth_layers=None, depth_spacing=None, trace_resolution=None, smoothing=0.0, dip=None, dip_direction='right', symbol=None)[source]¶

Create a surface by extruding a polyline trace to depth.

This is the recommended way to create fault surfaces from map-view trace data. The trace polyline is optionally interpolated to a target resolution, then extruded to depth to create a ruled surface with explicit triangulation.

When dip is specified, the surface follows a parabolic curve from vertical at the surface to the given dip angle at maximum depth. The offset is applied perpendicular to the local trace direction. This produces geologically realistic fault geometry where faults are steep near the surface and flatten at depth.

For geographic meshes, trace points are (lon, lat) in degrees and depths are physical distances below the ellipsoid surface. The ellipsoid parameters are read from the mesh automatically.

For Cartesian meshes, trace points are (x, y) in model coordinates and depth is the z-coordinate (downward positive).

Parameters:

name (str) – Surface identifier.
mesh (Mesh) – Computational mesh. For geographic meshes the ellipsoid is read from mesh.CoordinateSystem.ellipsoid.
trace_points (ndarray, shape (N, 2)) – Surface trace polyline. For geographic meshes: (lon, lat) in degrees. For Cartesian meshes: (x, y) in the same coordinate space as control points (model or physical, depending on whether units are active).
depth_range (tuple of (min_depth, max_depth)) – Depth extent for extrusion (positive downward). Each value can be a plain float (km for geographic, model coords for Cartesian) or a uw.quantity.
n_depth_layers (int, optional) – Number of depth levels (including surface and deepest). If None, computed from depth_spacing. If both are None, defaults to 7 layers.
depth_spacing (float or quantity, optional) – Spacing between depth layers. Alternative to n_depth_layers.
trace_resolution (float or quantity, optional) – Target point spacing along the trace. If None the trace is used as-is. For geographic meshes this is in km; for Cartesian meshes it is in model coordinates.
smoothing (float) – Spline smoothing parameter for splprep (0 = interpolate exactly through points, >0 allows smoothing).
dip (float, optional) –
Dip angle in degrees measured from horizontal. 90 means vertical (no offset), 45 means equal horizontal and vertical extent at maximum depth. When None (default), the surface is extruded vertically.

The surface profile is parabolic: vertical at the surface (depth = depth_range[0]) and reaching dip at the deepest layer.
dip_direction (str) – Direction of the dip offset relative to the trace: "right" (default) or "left" when looking along the trace from first to last point.
symbol (str, optional) – Short LaTeX symbol for expressions (e.g. "F").

Returns:

A fully discretized surface ready for use.

Return type:

Surface

Examples

>>> # Geographic mesh — trace in lon/lat, depths in km
>>> trace = np.column_stack([lon_points, lat_points])
>>> s = Surface.from_trace(
...     "fault_1", mesh, trace,
...     depth_range=(uw.quantity(0, "km"), uw.quantity(30, "km")),
...     depth_spacing=uw.quantity(5, "km"),
...     trace_resolution=uw.quantity(3, "km"),
...     dip=60,  # 60° from horizontal, parabolic profile
... )

classmethod from_vtk(filename, mesh=None, name=None)[source]¶

Load surface from VTK file.

All point_data arrays in the VTK file are automatically wrapped as SurfaceVariables.

Args:: filename: Path to VTK file (.vtk or .vtp) mesh: Computational mesh (required for .sym access) name: Name for the surface. If None, uses filename stem.
Returns:: Surface with all variables from VTK file
Raises:: FileNotFoundError: If file doesn’t exist ImportError: If pyvista not available

Parameters:

filename (str)
mesh (Mesh)
name (str)

Return type:

Surface

compute_normals(consistent_normals=True)[source]¶

Recompute vertex and face normals.

Args:: consistent_normals: If True, attempt to make normals consistently oriented

Parameters:: consistent_normals (bool)
Return type:: None

flip_normals()[source]¶

Flip the direction of all normals.

This directly negates the normal vectors. Note that pyvista’s compute_normals() uses consistent orientation so we can’t rely on flip_faces() to reverse normals - we must negate them explicitly.

Return type:: None

SurfaceCollection¶

class underworld3.meshing.SurfaceCollection[source]¶

Bases: object

Collection of surfaces for combined operations.

Manages multiple Surface objects and provides methods to: - Compute minimum distance from mesh points to any surface - Transfer surface normals to mesh variables via nearest-neighbor - Create combined influence functions

Example:

>>> surfaces = uw.meshing.SurfaceCollection()
>>> surfaces.add(fault1)
>>> surfaces.add(fault2)
>>>
>>> # Compute combined distance field
>>> dist = surfaces.compute_distance_field(mesh)
>>>
>>> # Combined influence function
>>> eta = surfaces.influence_function(
...     width=0.05,
...     value_near=0.01,
...     value_far=1.0,
... )

__init__()[source]¶: Create an empty surface collection.

add(surface)[source]¶

Add a surface to the collection.

Args:: surface: Surface to add
Raises:: ValueError: If surface with same name already exists

Parameters:: surface (Surface)
Return type:: None

add_from_vtk(filename, mesh=None, name=None)[source]¶

Load and add a surface from VTK file.

Args:: filename: Path to VTK file mesh: Computational mesh name: Name for the surface. If None, uses filename stem.
Returns:: The loaded Surface

Parameters:

filename (str)
mesh (Mesh)
name (str)

Return type:

Surface

remove(name)[source]¶

Remove and return a surface from the collection.

Args:: name: Name of surface to remove
Returns:: The removed Surface
Raises:: KeyError: If surface not found

Parameters:: name (str)
Return type:: Surface

property names: List[str]¶: List of surface names.

compute_distance_field(mesh, distance_var=None, variable_name='surface_distance')[source]¶

Compute minimum distance from mesh points to any surface.

Args:: mesh: The mesh to compute distances on distance_var: Optional existing MeshVariable to store results variable_name: Name for new variable if distance_var is None
Returns:: MeshVariable with distance values (scalar)

Parameters:

mesh (Mesh)
distance_var (MeshVariable)
variable_name (str)

Return type:

MeshVariable

transfer_normals(mesh, coords=None, normal_var=None, variable_name='surface_normals')[source]¶

Transfer surface normals to mesh points via nearest-neighbor.

For each mesh point, finds the closest surface face and copies that face’s normal vector.

Args:: mesh: The mesh to transfer normals to coords: Optional coordinates to query. If None, uses mesh.data normal_var: Optional existing MeshVariable variable_name: Name for new variable
Returns:: MeshVariable with normal vectors (3 components)

Parameters:

mesh (Mesh)
coords (np.ndarray)
normal_var (MeshVariable)
variable_name (str)

Return type:

MeshVariable

influence_function(mesh, width, value_near=1.0, value_far=0.0, profile='step')[source]¶

Create combined influence function from all surfaces.

Args:: mesh: Computational mesh width: Characteristic width of the transition zone value_near: Value at/near surfaces value_far: Value far from surfaces profile: Transition profile (step, linear, gaussian, smoothstep)
Returns:: sympy.Expr based on combined distance field

Parameters:

mesh (Mesh)
width (float)
value_near (Union[float, sympy.Expr])
value_far (Union[float, sympy.Expr])
profile (str)

Return type:

sympy.Expr

refinement_metric(mesh, h_near, h_far, width=None, profile='linear', variable_name='fault_metric')[source]¶

Create a combined refinement metric for mesh adaptation.

Computes a single metric field using the minimum unsigned distance across all surfaces in the collection. This creates 2 MeshVariables (distance + metric) regardless of how many surfaces are in the collection, avoiding the O(N²) DM-rebuild cost of computing per-surface metrics in a loop.

Parameters:

mesh (Mesh) – The computational mesh.
h_near (float or quantity) – Target edge length near any surface.
h_far (float or quantity) – Target edge length far from surfaces.
width (float or quantity, optional) – Transition distance. Defaults to 2 * h_far.
profile (str) – "linear", "smoothstep", or "gaussian".
variable_name (str) – Name for the metric MeshVariable.

Returns:

Scalar metric field suitable for mesh.adapt().

Return type:

MeshVariable

compute_nearest_fields(mesh, fault_width=None, normal_var_name='fault_n', id_var_name='fault_id', distance_var_name='d_faults', weight_var_name='fault_w')[source]¶

Compute per-node nearest-surface fields for rheology.

For every mesh node, finds the closest surface vertex across all surfaces and returns:

normal: unit normal of the nearest surface vertex
id: surface identifier for the nearest surface
distance: minimum unsigned distance to any surface
weight (optional): Gaussian influence exp(-0.5*(d/width)²)

This consolidates what would otherwise be manual KDTree code into a single reusable call, and creates only 3–4 MeshVariables total.

Parameters:

mesh (Mesh) – The computational mesh.
fault_width (float or quantity, optional) – Gaussian half-width for the weight field. If None, the weight field is omitted from the returned dict.
normal_var_name (str) – Names for the output MeshVariables.
id_var_name (str) – Names for the output MeshVariables.
distance_var_name (str) – Names for the output MeshVariables.
weight_var_name (str) – Names for the output MeshVariables.

Returns:

Keys: "normal", "id", "distance", and optionally "weight". Values are MeshVariable instances.

Return type:

dict

create_weakness_function(distance_var, fault_width, eta_weak=0.01, eta_background=1.0)[source]¶

Create Piecewise viscosity function for fault weakness.

DEPRECATED: Use influence_function() instead.

Args:: distance_var: MeshVariable containing distances fault_width: Width of the weak zone eta_weak: Viscosity within weak zone eta_background: Viscosity outside weak zone
Returns:: sympy.Piecewise expression

Parameters:

distance_var (MeshVariable)
fault_width (float)
eta_weak (float)
eta_background (float)

Return type:

sympy.Expr

SurfaceVariable¶

Variables defined on surface meshes.

class underworld3.meshing.SurfaceVariable[source]¶

Bases: object

Variable defined on surface vertices, stored in pyvista point_data.

Provides direct access to per-vertex data via .data property, and symbolic access for expressions via .sym property.

The .sym property uses a proxy MeshVariable that is lazily interpolated from surface vertices to mesh nodes when accessed.

Unit awareness is provided via UnitAwareArray wrapper on .data access. Distance-based masking is available via .mask property when mask_width is set.

Attributes:

name: Variable name (also the key in pyvista point_data) surface: Parent Surface object size: Number of components per vertex units: Optional units for this variable (e.g., “Pa”, “m/s”)

Example:

>>> friction = surface.add_variable("friction", size=1, mask_width=0.1)
>>> friction.data[:] = 0.6  # Set values on surface vertices
>>> friction.data[weak_mask] = 0.3  # Set weak zone values
>>>
>>> # Use in expressions with explicit masking
>>> eta = friction.sym[0] * friction.mask

__init__(name, surface, size=1, proxy_degree=1, existing=False, units=None, mask_width=None, mask_profile='gaussian')[source]¶

Create a variable on surface vertices.

Args:: name: Variable name (key in pyvista point_data) surface: Parent Surface object size: Number of components per vertex (1 for scalar, 3 for vector) proxy_degree: Degree of proxy MeshVariable for .sym access existing: If True, wraps existing point_data (for loading from VTK) units: Optional units for this variable (e.g., “Pa”, “m/s”) mask_width: Width for distance-based mask (enables .mask property) mask_profile: Profile for mask function (“step”, “linear”, “gaussian”, “smoothstep”)

Parameters:

name (str)
surface (Surface)
size (int)
proxy_degree (int)
existing (bool)
units (str | None)
mask_width (float | None)
mask_profile (str)

property units: str | None¶: Units for this variable (None if dimensionless).

property has_units: bool¶: Check if this variable has units.

property data¶

Direct access to vertex data with optional unit awareness.

Returns:: UnitAwareArray if units are set, otherwise plain numpy array. Shape is (n_vertices,) for size=1, or (n_vertices, size) otherwise.
Note:: When units are set, modifications via array operations automatically sync to pyvista storage and mark proxy as stale. For raw numpy arrays, call mark_stale() after modifications.

mark_stale()[source]¶

Mark the proxy as stale so it will be recomputed on next .sym access.

Return type:: None

property mask: Expr¶

1 near surface, 0 far away.

Uses the surface’s signed distance field with the configured profile. Must be explicitly applied in expressions: use var.sym * var.mask.

Returns:

sympy.Expr representing the mask (1 near, 0 far)

Raises:

ValueError: If mask_width was not set in add_variable()

Example:

>>> friction = surface.add_variable("friction", mask_width=0.1)
>>> friction.data[:] = 0.6
>>> eta = friction.sym[0] * friction.mask  # Masked value

Type:: Distance-based mask

property sym: MutableDenseMatrix¶

Symbolic representation for use in expressions.

Returns:: sympy.Matrix that can be used in Underworld expressions.
Note:: On first access (or after mark_stale()), this triggers interpolation from surface vertices to mesh nodes. The interpolation uses inverse distance weighting from nearby surface vertices.