mesh

The dream.mesh module provides mesh-generation utilities, buffer/sponge layer infrastructure, and the condition containers shared by all dream solvers. It wraps Netgen’s meshing API so that common aeroacoustic mesh layouts can be constructed with a few function calls.

Predefined mesh generators

Several helper functions produce ready-to-use Mesh objects for standard geometries:

get_cylinder_mesh() — unstructured mesh around a circular cylinder, optionally with a body-fitted boundary layer, a graded transition region and an outer sponge annulus.
get_cylinder_omesh() — fully structured O-grid ring mesh with a fixed number of elements in the polar and radial directions.
get_2d_naca_occ_profile() and get_3d_naca_occ_profile() — NACA 4-digit airfoil profiles as OCC shapes that can be embedded in any surrounding geometry before meshing.

Structured mesh generators

get_rectangular_mesh() and get_structured_cylinder_mesh() build fully structured meshes from explicit nodal coordinate arrays. The caller supplies one array per spatial direction for each domain region; the function takes their union as the global set of mesh nodes and connects them into quad or triangle elements:

import numpy as np
from dream.mesh import get_rectangular_mesh, get_nodal_points

nx = np.linspace(0, 4, 41)
ny = get_nodal_points(21, distribution='tanh', beta=3) - 0.5   # cluster near walls

domains    = [("channel", (nx, ny))]
boundaries = [("inflow",  (np.array([0.0]), ny)),
              ("outflow", (np.array([4.0]), ny)),
              ("wall",    (nx, np.array([-0.5]))),
              ("wall",    (nx, np.array([0.5])))]

mesh = get_rectangular_mesh(domains, boundaries)

The polar variant get_structured_cylinder_mesh() works the same way but accepts radial and angular coordinates \((r, \varphi)\) instead of Cartesian ones.

get_nodal_points() returns 1-D node distributions in \([0,1]\) with several clustering options ('uniform', 'cosine', 'polynomial', 'tanh', 'exponential') that can be used to concentrate grid points near boundaries.

Buffer and sponge layer infrastructure

dream uses e.g. buffer layers to implement non-reflecting far-field conditions. Two mechanisms are available:

Grid deformation (GridDeformation): the computational mesh is stretched inside the buffer region so that outgoing waves encounter progressively coarser resolution and are damped numerically. The deformation is described by a GridMapping, which maps a BufferCoord in the computational domain to a new position in the physical domain:

from dream.mesh import BufferCoord, GridMapping, GridDeformation

x   = BufferCoord.x(x0=3.0, xn=5.0)          # buffer extends from x=3 to x=5
map = GridMapping.exponential(scale=5, coordinate=x)
deformation = GridDeformation(x=map, order=3)

Four mapping types are available: none() (identity), linear(), exponential(), and tangential(). Polar and spherical buffer coordinates are supported via polar() and spherical().

Sponge layer (SpongeLayer): a volumetric penalty term \(\sigma(\vec{x})\,(\vec{U} - \vec{U}_\infty)\) is added to the right-hand side inside the sponge region, where \(\sigma\) is the sponge weight function provided by SpongeFunction:

from dream.mesh import SpongeFunction

sigma = SpongeFunction.polynomial(weight=2.0, x=x, order=3)
sponge = SpongeLayer(function=sigma, target_state={"rho": 1.0, "u": (1.0, 0.0)})

The p-type variant PSpongeLayer additionally reduces the local polynomial order from a high value at the inner edge of the layer down to a low value at the outer edge, introducing extra numerical dissipation alongside the explicit damping term.

Condition containers

BoundaryConditions and DomainConditions map mesh regions to Condition instances using NGSolve’s pipe-separated region-name syntax:

from dream.mesh import BoundaryConditions, DomainConditions, Periodic, Initial

bcs = BoundaryConditions(mesh, options=[Periodic])
bcs["inflow|outflow"] = Periodic()

dcs = DomainConditions(mesh, options=[SpongeLayer, GridDeformation])
dcs["sponge"] = SpongeLayer(function=sigma, target_state={"rho": 1.0})

Both containers warn when a region pattern does not match any mesh region, and they flag multiple conditions set on the same region.

`get_nodal_points`(n[, distribution])	Generate 1D nodal points in [0, 1] with various clustering distributions.
`get_cylinder_mesh`([radius, sponge_layer, ...])	Generates an unstructured mesh around a circular cylinder for aeroacoustic simulations.
`get_cylinder_omesh`(ri, ro, n_polar, n_radial)	Generates a ring mesh with a given inner and outer radius.
`get_rectangular_mesh`(domains, boundaries[, ...])	Generates a structured rectangular mesh from explicit nodal coordinate arrays.
`get_structured_cylinder_mesh`(domains, boundaries)	Generates a structured cylinder mesh based on given radial and angular coordinates.
`get_2d_naca_occ_profile`(number[, AoA, LE, ...])	Returns a 2D NACA airfoil profile with leading edge at position (0, 0).
`get_3d_naca_occ_profile`(number, depth[, ...])	Returns a 3D NACA airfoil profile with leading edge at position (0, 0, 0).
`get_chord_naca_4digit_series_coordinates`(number)	Returns a 2D NACA airfoil profile with leading edge at position (0, 0).
`BufferCoord`(x, x0, xn, shift)	One-dimensional coordinate used in buffer layers.
`SpongeFunction`	Defines some predefined sponge function in a buffer layer
`GridMapping`(x, map)	Mapping used for mesh deformation purposes.
`Condition`()	Base class for a named condition that can be set on a mesh region.
`Periodic`()	Marks a boundary region as periodic.
`Initial`([fields, bonus_int_order])	Sets the initial condition on a domain region.
`Perturbation`([fields])	Adds a perturbation to the initial condition on a domain region.
`Buffer`([function, order])	Base class for a domain buffer condition that contributes an auxiliary term to the weak form.
`GridDeformation`([x, y, z, dim, order])	Specifies a mesh deformation mapping in a domain region.
`SpongeLayer`([target_state, function, order])	Adds a sponge (damping) term in a domain region to absorb outgoing waves.
`PSpongeLayer`([high_order, low_order, ...])	Polynomial sponge layer that additionally reduces the local polynomial order.
`Conditions`(regions, mesh, options)	Container that maps mesh regions to `Condition` instances.
`BoundaryConditions`(mesh, options)	Manages conditions on mesh boundary (edge/face) regions.
`DomainConditions`(mesh, options)	Manages conditions on mesh domain (material) regions.

get_nodal_points(n, distribution='uniform', **kwargs)

Generate 1D nodal points in [0, 1] with various clustering distributions.

Parameters

nint

Number of points.

distributionstr, optional

Type of nodal distribution. Options: ‘uniform’ (default), ‘cosine’, ‘polynomial’, ‘tanh’, ‘exponential’.

**kwargsdict, optional

Additional parameters for specific distributions:

polynomial: p (default=2)
tanh: beta (default=2.5)
exponential: a (default=4)

Returns

xndarray: Array of nodal points in [0, 1].

get_cylinder_mesh(radius: float = 0.5, sponge_layer: bool = False, boundary_layer_levels: int = 5, boundary_layer_thickness: float = 0.0, transition_layer_levels: int = 5, transition_layer_growth: float = 1.4, transition_radial_factor: float = 6, farfield_radial_factor: float = 50, sponge_radial_factor: float = 60, wake_maxh: float = 2, farfield_maxh: float = 4, sponge_maxh: float = 4, bnd: tuple[str, str, str] = ('inflow', 'outflow', 'cylinder'), mat: tuple[str, str] = ('sound', 'sponge'), curve_layers: bool = False, grading: float = 0.3)

Generates an unstructured mesh around a circular cylinder for aeroacoustic simulations.

The mesh is built from up to four concentric regions:

Boundary layer (optional): a stack of \(n\) thin structured rings of uniform thickness boundary_layer_thickness / boundary_layer_levels fitted to the cylinder surface. Enabled when boundary_layer_thickness > 0.
Transition region: a graded layer between the cylinder (or boundary layer) and the farfield, whose element size grows with exponent transition_layer_growth out to radius transition_radial_factor * radius.
Farfield region: the acoustically active region up to radius farfield_radial_factor * radius, resolved with element size farfield_maxh. A downstream wake patch of width wake_maxh is cut into this region.
Sponge layer (optional): an annular absorbing region between farfield_radial_factor * radius and sponge_radial_factor * radius. Enabled when sponge_layer=True.

Parameters:

radius (float, optional) – Radius of the cylinder, defaults to 0.5
sponge_layer (bool, optional) – Include an outer sponge layer, defaults to False
boundary_layer_levels (int, optional) – Number of elements in the radial boundary-layer stack, defaults to 5
boundary_layer_thickness (float, optional) – Total thickness of the boundary layer (0 disables it), defaults to 0
transition_layer_levels (int, optional) – Number of graded rings in the transition region, defaults to 5
transition_layer_growth (float, optional) – Growth exponent of the transition region, defaults to 1.4
transition_radial_factor (float, optional) – Outer radius of the transition region as a multiple of radius, defaults to 6
farfield_radial_factor (float, optional) – Outer radius of the farfield region as a multiple of radius, defaults to 50
sponge_radial_factor (float, optional) – Outer radius of the sponge layer as a multiple of radius, defaults to 60
wake_maxh (float, optional) – Maximum element size in the downstream wake patch, defaults to 2
farfield_maxh (float, optional) – Maximum element size in the farfield region, defaults to 4
sponge_maxh (float, optional) – Maximum element size in the sponge layer, defaults to 4
bnd (tuple[str, str, str], optional) – Names of the inflow, outflow and cylinder boundaries, defaults to (‘inflow’, ‘outflow’, ‘cylinder’)
mat (tuple[str, str], optional) – Names of the acoustic and sponge domain materials, defaults to (‘sound’, ‘sponge’)
curve_layers (bool, optional) – Curve the boundary edges to follow the cylindrical geometry, defaults to False
grading (float, optional) – Netgen grading parameter controlling element size variation, defaults to 0.3

Returns:

The generated mesh

Return type:

ngs.Mesh

get_cylinder_omesh(ri: float, ro: float, n_polar: int, n_radial: int, geom: float = 1, bnd: tuple[str, str, str] = ('cylinder', 'left', 'right'), dom: str = 'default') → Mesh

Generates a ring mesh with a given inner and outer radius.

Parameters:

ri (float) – Inner radius of the ring
ro (float) – Outer radius of the ring
n_polar (int) – Number of elements in the polar direction
n_radial (int) – Number of elements in the radial direction
geom (float, optional) – Geometric factor for the radial direction, defaults to 1

Returns:

Ring mesh

Return type:

ngs.Mesh

get_rectangular_mesh(domains: list[tuple[str, tuple[ndarray, ndarray]]], boundaries: list[tuple[str, tuple[ndarray, ndarray]]], quads: bool = True, periodic_x: bool = False, periodic_y: bool = False) → Mesh

Generates a structured rectangular mesh from explicit nodal coordinate arrays.

Mesh points are placed at every combination of the union of all x and y coordinates supplied across all domains. Elements are then created for each domain by selecting the mesh points that fall within the domain’s x and y extents and subdividing the resulting rectangular patches into quads (default) or pairs of triangles.

Both domains and boundaries are sequences of (name, (x_coords, y_coords)) pairs:

domains: x_coords and y_coords are 1-D arrays spanning the extent of the domain. Several domains may share edge coordinates; the union of all coordinate arrays determines the global mesh nodes.
boundaries: a boundary is identified by the range of its x and y coordinate array. To select a vertical edge at \(x = x_0\), pass x_coords = np.array([x_0]) and y_coords spanning the full y extent of that edge.

Parameters:

domains (list[tuple[str, tuple[np.ndarray, np.ndarray]]]) – Sequence of (name, (x_coords, y_coords)) pairs defining the domain regions
boundaries (list[tuple[str, tuple[np.ndarray, np.ndarray]]]) – Sequence of (name, (x_coords, y_coords)) pairs defining the boundary edges
quads (bool, optional) – Use quadrilateral elements; if False, each quad is split into two triangles, defaults to True
periodic_x (bool, optional) – Identify left and right boundary nodes as periodic, defaults to False
periodic_y (bool, optional) – Identify bottom and top boundary nodes as periodic, defaults to False

Returns:

The generated structured mesh

Return type:

ngs.Mesh

get_structured_cylinder_mesh(domains: dict[str, tuple[ndarray, ndarray]], boundaries: dict[str, tuple[ndarray, ndarray]], close_angular: bool = True, curve_all=False, quads: bool = True) → Mesh

Generates a structured cylinder mesh based on given radial and angular coordinates.

The mesh is constructed by defining domains and boundaries in (r, phi) coordinates. Phi is assumed to be periodic, i.e., the first and last angular coordinates are connected, therefore pass only unique angular coordinates.

Parameters:

domains (dict[str, tuple[np.ndarray, np.ndarray]]) – Dictionary of domain names and their extents in (r, phi) coordinates
boundaries (dict[str, tuple[np.ndarray, np.ndarray]]) – Dictionary of boundary names and their extents in (r, phi) coordinates
close_angular (bool, optional) – Whether the phi direction is closed, defaults to True
curve_all (bool, optional) – Whether to curve also the non-named boundaries, defaults to False
quads (bool, optional) – Whether to use quadrilateral elements, defaults to True

Returns:

The generated structured cylinder mesh

Return type:

ngs.Mesh

get_2d_naca_occ_profile(number: str | int, AoA=0, LE=(0, 0), chord: float = 1.0, n: int = 600) → TopoDS_Shape

Returns a 2D NACA airfoil profile with leading edge at position (0, 0).

Parameters:

number (str | int) – NACA digit number
AOA (float, optional) – Angle of attack in degrees, defaults to 0
LE (tuple, optional) – Leading edge coordinates, defaults to (0, 0)
chord (float, optional) – Chord length, defaults to 1
n (int, optional) – Number of coordinates, defaults to 600

Returns:

NACA airfoil profile

Return type:

occ.TopoDS_Shape

get_3d_naca_occ_profile(number: str | int, depth: float, AoA=0, LE=(0, 0), scale: int = 1, n: int = 600, periodic: bool = True) → TopoDS_Shape

Returns a 3D NACA airfoil profile with leading edge at position (0, 0, 0).

Parameters:

number (str | int) – NACA digit number
depth (float) – Depth of the airfoil
AOA (float, optional) – Angle of attack in degrees, defaults to 0
LE (tuple, optional) – Leading edge coordinates, defaults to (0, 0)
scale (int, optional) – Chord scale factor, defaults to 1
n (int, optional) – Number of coordinates, defaults to 600
periodic (bool, optional) – If True, the airfoil is periodic in the Z direction, defaults to True

Returns:

NACA airfoil profile as occ.Shape

Return type:

occ.TopoDS_Shape

get_chord_naca_4digit_series_coordinates(number: str | int, LE=(0, 0), chord: float = 1.0, n: int = 600, nodal_distribution: str = 'uniform') → list[tuple[float, float, 0]]

Returns a 2D NACA airfoil profile with leading edge at position (0, 0).

Parameters:

number (str | int) – NACA digit number
LE (tuple, optional) – Leading edge coordinates, defaults to (0, 0)
chord (float, optional) – Chord length, defaults to 1
n (int, optional) – Number of coordinates, defaults to 600

Returns:

NACA airfoil points

Return type:

list[tuple[float, float, 0]]

class BufferCoord(x: CoefficientFunction, x0: float | CoefficientFunction, xn: float | CoefficientFunction, shift: tuple[float])

One-dimensional coordinate used in buffer layers.

The physical coordinate is truncated at the starting and end point.

If the coordinate system underwent a translation, it can be set by the shift parameter. This is useful e.g. when dealing with radial or spherical coordinates.

\[\begin{split}\tilde{x} = \begin{cases} 0 & x <= x_0 \\ \frac{x - x_0}{x_n - x_0} & x_0 < x < x_n \\ 1 & x > x_n \end{cases}\end{split}\]

class SpongeFunction: Defines some predefined sponge function in a buffer layer

class GridMapping(x: CoefficientFunction, map: Callable[[float | CoefficientFunction], float | CoefficientFunction])

Mapping used for mesh deformation purposes.

One-dimensional buffer coordinates are mapped from computational to physical buffer coordinates.

classmethod none(coordinate: CoefficientFunction)

Returns a zero grid mapping.

This is mainly used as consistency between mappings in different coordinates. It can be seen as the python equivalent None.

classmethod linear(scale: float, coordinate: BufferCoord)

Returns a linear grid mapping.

The thickness of the buffer layer is scaled by the factor ‘scale’.

\[f(x) = scale * (x - x_0) + x_0\]

classmethod exponential(scale: float, coordinate: BufferCoord)

Returns an exponential grid mapping.

The thickness of the buffer layer is scaled by the factor ‘scale’.

The constants c_0 and c_1 are determined by a fixpoint iteration.

\[f(x) = c_0 * (1 - \exp^{c_1 (x - x_0)}) + x_0\]

classmethod tangential(scale: float, coordinate: BufferCoord)

Returns a tangential grid mapping.

The thickness of the buffer layer is scaled by the factor ‘scale’.

The constants c_0 and c_1 are determined by a fixpoint iteration.

\[f(x) = c_0 * (\tan{c_1 (x - x_0)}) + x_0\]

polar_to_cartesian() → tuple[GridMapping, GridMapping]: Automates the transformation from a polar mapping to a cartesian one

class Condition

Base class for a named condition that can be set on a mesh region.

Conditions are set on a region (boundary or domain) via a Conditions container, e.g. BoundaryConditions or DomainConditions. Instances compare by identity, such that the same condition object can be associated with several regions while remaining distinguishable from another condition of the same type set on a different region.

class Periodic: Marks a boundary region as periodic.

class Initial(fields: ngsdict | None = None, bonus_int_order: int = 0)

Sets the initial condition on a domain region.

Holds the initial fields, which are projected onto the solution space at the start of a TransientRoutine, PseudoTimeSteppingRoutine, or IMEXTimeRoutine.

property fields: ngsdict: Returns the fields of the initial condition

class Perturbation(fields: ngsdict | None = None)

Adds a perturbation to the initial condition on a domain region.

Holds the perturbation fields, which are added on top of the initial condition, e.g. to seed a transient simulation with a localized disturbance.

property fields: ngsdict: Returns the fields of the initial condition

class Buffer(function: CoefficientFunction | None = None, order: int = 0)

Base class for a domain buffer condition that contributes an auxiliary term to the weak form.

Subclasses hold a function (the buffer weight or deformation field) and a polynomial order used to project it onto a suitable FE space. Concrete subclasses are GridDeformation and SpongeLayer.

class GridDeformation(x: GridMapping = None, y: GridMapping = None, z: GridMapping = None, dim: int = 2, order: int = 0)

Specifies a mesh deformation mapping in a domain region.

The deformation is defined by up to three GridMapping objects, one per spatial coordinate. Each mapping translates a BufferCoord in the computational domain to a deformed position in the physical domain; the difference gives the deformation vector that is projected onto a VectorH1 space and applied to the mesh via get_grid_deformation_function().

class SpongeLayer(target_state: ngsdict = None, function: CoefficientFunction = None, order=0)

Adds a sponge (damping) term in a domain region to absorb outgoing waves.

The sponge drives the solution toward a prescribed target_state with a weight given by function. The weight is projected onto an L2 space via get_sponge_layer_function() and added as a volumetric penalty term in the weak form.

property target_state: ngsdict: Returns the fields of the target state

class PSpongeLayer(high_order: int = 0, low_order: int = 0, target_state: ngsdict = None, function: CoefficientFunction = None, order: int = 0)

Polynomial sponge layer that additionally reduces the local polynomial order.

In addition to the standard damping of SpongeLayer, a PSpongeLayer gradually lowers the polynomial order from high_order (at the inner edge of the layer) down to low_order (at the outer edge). This p-type coarsening reduces the resolution inside the sponge and thereby introduces additional numerical dissipation, which enhances absorption of outgoing waves.

class Conditions(regions: list[str], mesh: Mesh, options: list[Condition])

Container that maps mesh regions to Condition instances.

Regions are identified by their string names as returned by GetBoundaries() or GetMaterials(). A condition is associated with one or more regions by assigning it via conditions['region_a|region_b'] = SomeCondition(). The pipe-separated pattern syntax follows the NGSolve convention.

Assigning the same Condition object to multiple calls merges all matching regions under a single condition instance, which can be retrieved together through items().

class BoundaryConditions(mesh: Mesh, options: list[Condition])

Manages conditions on mesh boundary (edge/face) regions.

Initialised from the list of boundary names returned by GetBoundaries(). Provides helpers to retrieve the periodic boundaries and the non-periodic domain boundaries as NGSolve region patterns.

get_domain_boundaries() → str

Returns a list or pattern of the domain boundaries!

The domain boundaries are deduced by the current set boundary conditions, while periodic boundaries are neglected!

class DomainConditions(mesh: Mesh, options: list[Condition])

Manages conditions on mesh domain (material) regions.

Initialised from the list of material names returned by GetMaterials(). Provides helpers to assemble the grid-deformation field, the sponge-layer weight, and the p-sponge weight as NGSolve GridFunction objects, as well as utilities to reduce the polynomial order of an L2 or FacetFESpace elementwise inside p-sponge regions.

mesh

Predefined mesh generators

Structured mesh generators

Buffer and sponge layer infrastructure

Condition containers

Parameters

Returns

Examples