Array utilities¶

Array functions shared across the codebase.

Each function below is a standalone operation on plain NumPy arrays. This module acts a a place to collect common conventions, patterns, and conversions so that they can be used consistently across the code.

Flow angle and coordinate conversions¶

Paired functions convert between the native ember data representation and and other conventions, each the inverse of its partner. There are utilities for conversion between flow angles and velocity components, between Cartesian and polar coordinates, and between polar and pseudo-Cartesian coordinates.

`angles_to_components`(V, Alpha, Beta)	Resolve velocity magnitude into polar components.
`components_to_angles`(Vx, Vr, Vt)	Convert velocity components to velocity magnitude and flow angles.
`cart_to_pol`(xyz, Vxyz[, perm, signs])	Convert Cartesian coordinates and velocities to polar form.
`pol_to_cart`(xrt, Vxrt[, perm, signs])	Convert polar coordinates and velocities to Cartesian form.
`pol_to_pseudocart`(xrt[, inplace])	Convert (x, r, θ) coordinates to pseudo-Cartesian (x, r, rθ).

Vector and matrix operations¶

These functions are intended for use batches of vectors or matrices stacked along leading axes, with the last one or two axes holding the vector or matrix components. Some are thin wrappers around numpy.einsum, but for operations on the hot path of a solver loop they are implemented in Fortran for better performance.

`dot`(a, b)	Dot product of two arrays along the last axis.
`vecnorm`(x)	Calculate the norm of a vector array along the last axis.
`matmat`(A, B)	Matrix-matrix multiplication over trailing dimensions.
`matvec`(A, b[, out])	Matrix-vector multiplication using einsum over trailing dimensions.

Grid construction¶

The below functions are used to build coordinate arrays for the solver. For simple duct geometries, they allow convenient construction of uniform or clustered xrt arrays. The resample() is useful for coarsening or refining the distribution of an existing grid.

`meshgrid3`(xv, rv, tv)	Create 3D coordinate rectangular meshgrid.
`linmesh3`(x, r, t, shape)	Create 3D coordinate meshgrid from ranges and shape.
`cosine_cluster`(n)	Generate cosine-clustered points from 0 to 1, dense at both ends.
`cluster`(ni, ER, dmax)	Generate geometrically spaced points from 0 to 1 with expansion ratio and max spacing.
`cluster_symmetric`(n, ER[, dmax])	Generate geometrically spaced points from 0 to 1, dense at both ends.
`resample`(factor, vector[, i_crit])	Resampled a vector with specified factor, optionally preserving critical points.

Array allocation and buffers¶

Many of these functions are are thin wrappers around their numpy namesakes that fix the dtype and memory layout that weember standardises on: Fortran order and single-precision float32. Calling these rather than the bare NumPy functions signals intent and keeps that convention in one place – use them for any new array creation. There are also functions for memory management: allocation, or buffer creation and reuse.

`zeros`(shape[, dtype])	Zero-filled array in standard layout and dtype.
`array`(x[, dtype])	Copy array data into standard layout and dtype.
`empty`(shape)	Uninitialised array in standard layout and dtype.
`full`(shape, fill_value)	Constant-filled array in standard layout and dtype.
`allocate_or_reuse`(out, shape[, dtype])	Allocate output array if not provided, otherwise reuse existing array.
`bcast_if_needed`(a, shape)	Broadcast `a` to `shape` only if it doesn't already have it.
`carve_view`(buf, *shapes)	Carve one or more zero-copy Fortran-order views from a buffer of any shape.
`rss_bytes`()	Resident and peak-resident memory of this process, in bytes.

Miscellaneous geometry¶

Assorted utilities for working with coordinates and bounding boxes.

`extent`(*args)	Calculate per-component min and max values.
`bounding_box`(xyz)	Calculate bounding box vertices from Cartesian coordinates.
`apply_perm_flip`(array, perm[, flip])	Apply permutation and flipping to array while preserving coordinate dimension.
`unwrap_meridional`(xr_curve, xr_query)	Unwrap meridional coordinates onto conformal distance along a curve.

util.profile(func)[source]¶

util.dot(a, b)[source]¶

Dot product of two arrays along the last axis.

Parameters:

a (Array, shape (..., n)) – First input array.
b (Array, shape (..., n)) – Second input array.

Returns:

prod – Dot product of the input arrays, components summed along the last axis.

Return type:

Array, shape (…)

util.angles_to_components(V, Alpha, Beta)[source]¶

Resolve velocity magnitude into polar components.

Uses numerically stable trigonometry for all angles including 90 degrees.

Flow angle conventions:

Alpha = 0°: No swirl (Vt = 0)
Beta = 0°: Pure axial flow (Vr = 0)
Beta = ±90°: Pure radial flow (Vx = 0)

Parameters:

V_rel (Array) – Relative velocity magnitude [m/s]
Alpha (Array) – Relative yaw angle (tangential flow direction) [degrees]
Beta (Array) – Pitch angle (radial flow direction) [degrees]

Returns:

Vx (Array) – Axial velocity component [m/s]
Vr (Array) – Radial velocity component [m/s]
Vt (Array) – Tangential velocity component [m/s]

util.components_to_angles(Vx, Vr, Vt)[source]¶

Convert velocity components to velocity magnitude and flow angles.

This is the inverse of angles_to_components(), obeying the conventions:

Alpha = 0°: No swirl (Vt_rel = 0)
Beta = 0°: Pure axial flow (Vr = 0)
Beta = ±90°: Pure radial flow (Vx = 0)

For velocity components with very small magnitudes, angles may be numerically unstable. Zero velocity returns (0, 0, 0).

Parameters:

Vx (Array) – Axial velocity component [m/s]
Vr (Array) – Radial velocity component [m/s]
Vt (Array) – Tangential velocity component [m/s]

Returns:

V (Array) – Velocity magnitude [m/s]
Alpha (Array) – Yaw angle (tangential flow direction) [degrees]
Beta (Array) – Pitch angle (radial flow direction) [degrees]

util.vecnorm(x)[source]¶

Calculate the norm of a vector array along the last axis.

Parameters:: x (Array, shape (..., n)) – Input array with last dimension as vector components.
Returns:: norm – Norm of the input vector array.
Return type:: Array, shape (…)

util.zeros(shape, dtype=<class 'numpy.float32'>)[source]¶

Zero-filled array in standard layout and dtype.

Equivalent to numpy.zeros, but defaulting to dtype=float32 and always order="F".

Parameters:

shape (int or tuple of int) – Shape of the new array.
dtype (numpy dtype, optional) – Data type for the new array. Default is np.float32.

Returns:

Fortran-ordered array of the given shape and dtype, filled with 0.

Return type:

Array

util.array(x, dtype=<class 'numpy.float32'>)[source]¶

Copy array data into standard layout and dtype.

Equivalent to numpy.array, but defaulting to dtype=float32 and always order="F".

Parameters:

x (array_like) – Data to copy into the new array.
dtype (numpy dtype, optional) – Data type for the new array. Default is np.float32.

Returns:

Fortran-ordered array of the given dtype, holding a copy of x.

Return type:

Array

util.empty(shape)[source]¶

Uninitialised array in standard layout and dtype.

Equivalent to numpy.empty, but always dtype=float32 and order="F". As with numpy.empty, contents are arbitrary until written – use this only where every element will be set before being read.

Parameters:: shape (int or tuple of int) – Shape of the new array.
Returns:: Fortran-ordered, float32 array of the given shape with uninitialised contents.
Return type:: Array

util.full(shape, fill_value)[source]¶

Constant-filled array in standard layout and dtype.

Equivalent to numpy.full, but always dtype=float32 and order="F".

Parameters:

shape (int or tuple of int) – Shape of the new array.
fill_value (scalar) – Value to fill every element with.

Returns:

Fortran-ordered, float32 array of the given shape, filled with fill_value.

Return type:

Array

util.allocate_or_reuse(out, shape, dtype=<class 'numpy.float32'>)[source]¶

Allocate output array if not provided, otherwise reuse existing array.

Helper function for functions that accept an optional out parameter. If out is None, allocates a new F-contiguous zero array with the specified shape and dtype. Otherwise returns the provided array.

Parameters:

out (Array or None) – Pre-allocated output array, or None to allocate new array.
shape (tuple) – Shape of array to allocate if out is None.
dtype (numpy dtype, optional) – Data type for new array. Default is np.float32.

Returns:

Either the provided out array or a newly allocated zero array.

Return type:

Array

Examples

>>> def my_function(x, out=None):
...     out = allocate_or_reuse(out, x.shape)
...     # ... compute results into out ...
...     return out

util.bcast_if_needed(a, shape)[source]¶

Broadcast a to shape only if it doesn’t already have it.

A shape-tuple comparison is far cheaper than numpy.broadcast_to() itself, so this is a near-zero-cost no-op on the common path where a is already shape – e.g. a hot-path setter that documents its inputs as “must broadcast to X” for callers that need it, but whose actual callers already pass exactly X-shaped arrays.

Parameters:

a (array_like) – Candidate array; may already have shape or be broadcastable to it.
shape (tuple) – Target shape.

Returns:

a unchanged if a.shape == shape, otherwise np.broadcast_to(a, shape).

Return type:

Array

util.carve_view(buf, *shapes)[source]¶

Carve one or more zero-copy Fortran-order views from a buffer of any shape.

Reinterprets successive spans of a Fortran-contiguous buf as arrays of the requested shapes, packing them end-to-end (view k starts where view k-1 ends) so every returned view aliases distinct storage and all may be held live simultaneously. Used to borrow differently-shaped scratch arrays out of one oversized solver buffer without allocating.

The offsets are computed internally, so a caller carving several coexisting slots cannot accidentally overlap or gap them. buf.reshape(-1, order="F") on such a buffer is itself a free view, so repeated calls stay zero-copy and there is no need to hoist a shared flat view at the call site.

Parameters:

buf (Array) – Fortran-contiguous source buffer to reinterpret. Must be large enough to hold the concatenation of all shapes (sum(prod(shape)) elements).
*shapes (tuple) – One shape per view, in packing order. Passing a single shape returns a single view; passing several returns a list, one per shape.

Returns:

A single F-order view when one shape is given, else a list of them, each aliasing a disjoint span of buf.

Return type:

Array or list of Array

Raises:

ValueError – If the packed views would not fit within buf.

util.rss_bytes()[source]¶

Resident and peak-resident memory of this process, in bytes.

Reads VmRSS/VmHWM from /proc/self/status, so the peak is the kernel’s own high-water mark: it covers transients that came and went between two calls, which sampling VmRSS alone would miss. Where /proc is unavailable (non-Linux) both come back as 0, so callers on a debug-logging path never have to branch.

Returns:: (rss, peak_rss) in bytes.
Return type:: tuple of int

util.meshgrid3(xv, rv, tv)[source]¶

Create 3D coordinate rectangular meshgrid.

This function combines the common pattern of creating a 3D meshgrid and stacking the results into a single coordinate array. It preserves the input dtype unless the inputs have mixed dtypes.

Parameters:

xv (array_like) – X-coordinate vector.
rv (array_like) – R-coordinate vector.
tv (array_like) – Theta-coordinate vector.

Returns:

Coordinate array with [x, r, t] components.

Return type:

Array, shape (len(xv), len(rv), len(tv), 3)

Examples

>>> xv = np.linspace(0, 1, 3)
>>> rv = np.linspace(1, 2, 4)
>>> tv = np.linspace(0, np.pi, 5)
>>> xrt = meshgrid3(xv, rv, tv)
>>> xrt.shape
(3, 4, 5, 3)

util.linmesh3(x, r, t, shape)[source]¶

Create 3D coordinate meshgrid from ranges and shape.

This function creates linearly spaced vectors from coordinate ranges and then generates a 3D meshgrid using meshgrid3(). It combines the common pattern of creating linspace vectors and then meshing them.

Parameters:

x (tuple or array_like) – Axial coordinate range [x_min, x_max].
r (tuple or array_like) – Radial coordinate range [r_min, r_max].
t (tuple or array_like) – Angular coordinate range [t_min, t_max].
shape (tuple) – Shape of the grid (ni, nj, nk).

Returns:

Coordinate array with (x, r, t) components, dtype=float32.

Return type:

Array, shape (ni, nj, nk, 3)

Examples

>>> xrt = linmesh3([0, 1], [1, 2], [0, np.pi], (3, 4, 5))
>>> xrt.shape
(3, 4, 5, 3)
>>> xrt.dtype
dtype('float32')

util.pol_to_pseudocart(xrt, inplace=False)[source]¶

Convert (x, r, θ) coordinates to pseudo-Cartesian (x, r, rθ).

In the pseudo-Cartesian system, the theta coordinate is multiplied by the radial coordinate to give rθ, which behaves like a Cartesian coordinate for distance calculations in cylindrical coordinates.

Parameters:

xrt (array_like) – Input coordinates with shape (…, 3) where the last dimension contains [x, r, θ] coordinates.
inplace (bool, optional) – If True, modify the input array in-place. If False (default), return a copy with the conversion applied.

Returns:

Coordinates with shape (…, 3) containing [x, r, rθ]. If inplace=True, returns the modified input array. If inplace=False, returns a new array.

Return type:

Array

util.extent(*args)[source]¶

Calculate per-component min and max values.

Parameters:: *args (array_like) – Arrays to compute the extent of, must have same number of dimensions, and same trailing dimension if ndim > 1.
Returns:: Extent array where extent[0, …] contains minimum values and extent[1, …] contains maximum values.
Return type:: Array, shape (2, …)

util.bounding_box(xyz)[source]¶

Calculate bounding box vertices from Cartesian coordinates.

Parameters:: xyz (array_like, shape (N, 3)) – Cartesian coordinates [x, y, z] with components on last axis
Returns:: Eight vertices of the bounding box representing all combinations of (min/max x, min/max y, min/max z)
Return type:: Array, shape (8, 3)

util.unwrap_meridional(xr_curve, xr_query)[source]¶

Unwrap meridional coordinates onto conformal distance along a curve.

Returns the conformal (or “blade-to-blade”) meridional coordinate

\[m' = \int \frac{\mathrm{d}m}{r}, \qquad \mathrm{d}m = \sqrt{\mathrm{d}x^2 + \mathrm{d}r^2},\]

integrated along xr_curve from its first point. Paired with \(\theta\) in radians, \(m'\) spans a conformal plane: angles and aspect ratios are preserved, so an aerofoil section drawn on \((m', \theta)\) axes keeps its shape at any radius.

The coordinate is a property of the curve, not of any grid: every query point is referred to one datum by one integration, so points from different blocks — or from different cuts of the same machine — come back on a common scale with nothing to match up afterwards. Where the curve is the same one passed to ember.cut.structured_meridional(), the coordinate and the surface it describes are exactly consistent.

Query points are projected onto the curve and their perpendicular offset discarded, so points lying near it rather than exactly on it — the triangulated output of ember.cut.unstructured(), say, or a blade surface — are handled without special-casing. A point beyond either end clamps to that end.

Parameters:

xr_curve (array_like, shape (n_point, 2)) – Meridional \((x, r)\) polyline defining the curve, in order. Must not touch the axis, where \(m'\) diverges.
xr_query (array_like, shape (..., 2)) – Meridional \((x, r)\) coordinates to evaluate.

Returns:

Conformal distance \(m'\) [-] at each query point, zero at the first point of the curve.

Return type:

Array, shape (…)

Raises:

ValueError – If the curve has fewer than two points, or reaches the axis.

Examples

>>> curve = np.array([[0.0, 2.0], [1.0, 2.0]])  # cylindrical, r = 2
>>> float(unwrap_meridional(curve, np.array([1.0, 2.0])))  # m' = m / r
0.5

util.cart_to_pol(xyz, Vxyz, perm=(0, 1, 2), signs=(1, 1, 1))[source]¶

Convert Cartesian coordinates and velocities to polar form.

Inverts the ember \((x, r, \theta)\) convention described at Coordinate system, i.e. the map implemented by pol_to_cart():

\[y = r \cos\theta, \qquad z = -r \sin\theta\]

This function takes Cartesian data stacked on the last axis, and returns polar data of the same shape, e.g. for setting a block’s raw coordinates and velocities xrt and Vxrt after reading a Cartesian solution.

perm and signs arguments, when supplied, reorient the input before conversion, for source data that doesn’t already follow our Cartesian axis order or sense. perm=(1, 0, 2) swaps \(x\) and \(y\) values at every point. signs[i] = -1 negates the (already permuted) component i. A mesh authored with \(z\) pointing the opposite way to our convention should be read with signs=(1, 1, -1), for example.

Parameters:

xyz (array_like, shape (..., 3)) – Cartesian coordinates (x, y, z) with components on last axis
Vxyz (array_like, shape (..., 3)) – Cartesian velocity components (Vx, Vy, Vz) with components on last axis
perm (tuple of int, optional) – Coordinate permutation (0, 1, 2) -> reordered indices. Default: (0, 1, 2)
signs (tuple of int, optional) – Coordinate signs (-1 or 1 for each axis). Default: (1, 1, 1)

Returns:

xrt (Array, shape (…, 3)) – Polar coordinates (x, r, t) with components on last axis
Vxrt (Array, shape (…, 3)) – Polar velocity components (Vx, Vr, Vt) with components on last axis

util.pol_to_cart(xrt, Vxrt, perm=(0, 1, 2), signs=(1, 1, 1))[source]¶

Convert polar coordinates and velocities to Cartesian form.

The reverse of cart_to_pol(). Implements the ember \((x, r, \theta)\) convention described at Coordinate system:

\[y = r \cos\theta, \qquad z = -r \sin\theta\]

The inputs are typically read from xrt and Vxrt; the returned arrays are Cartesian data of the same shape.

perm and signs reorient the output after conversion, for a target convention that isn’t ember’s own axis order or sense. perm[i] says which converted component supplies output component i, and then signs[i] = -1 negates the (already permuted) component i.

Note that reorientation is applied on the way in by cart_to_pol() but on the way out here, so the two are not simply symmetric for general perm and signs. ember.grid.Grid.align_cart_unstr() detects the correct orientation to align with a destination grid automatically when given Cartesian data.

Parameters:

xrt (array_like, shape (..., 3)) – Polar coordinates (x, r, t) with components on last axis
Vxrt (array_like, shape (..., 3)) – Polar velocity components (Vx, Vr, Vt) with components on last axis
perm (tuple of int, optional) – Coordinate permutation (0, 1, 2) -> reordered indices. Default: (0, 1, 2)
signs (tuple of int, optional) – Coordinate signs (-1 or 1 for each axis). Default: (1, 1, 1)

Returns:

xyz (Array, shape (…, 3)) – Cartesian coordinates (x, y, z) with components on last axis
Vxyz (Array, shape (…, 3)) – Cartesian velocity components (Vx, Vy, Vz) with components on last axis

util.matmat(A, B)[source]¶

Matrix-matrix multiplication over trailing dimensions.

Performs matrix multiplication on stacks of matrices where the matrices are stored in the trailing dimensions. This is optimized for arrays with matrix data in the last two dimensions and arbitrary leading dimensions.

Parameters:

A (Array, shape (..., m, k)) – First input array with matrices in trailing dimensions.
B (Array, shape (..., k, n)) – Second input array with matrices in trailing dimensions.

Returns:

Result of matrix multiplication A @ B for each corresponding pair of matrices in the trailing dimensions. Uses f32 precision and Fortran ordering for optimal performance.

Return type:

Array, shape (…, m, n)

Examples

>>> # Stack of 2x2 matrices
>>> A = np.random.randn(10, 5, 2, 2).astype(np.float32, order='F')
>>> B = np.random.randn(10, 5, 2, 2).astype(np.float32, order='F')
>>> C = matmat(A, B)  # Shape: (10, 5, 2, 2)

>>> # Single matrix multiplication
>>> A = np.eye(3, dtype=np.float32, order='F')
>>> B = np.ones((3, 3), dtype=np.float32, order='F')
>>> C = matmat(A, B)  # C = B

util.rotation_matrices(chi)[source]¶

Build a paired 2x2 rotation matrix and its inverse from a meridional-plane angle.

Both RevolutionPatch’s interface frame and resolve_to_interface() rotate a velocity pair \((V_x, V_r)\) by the same convention, so this is the one place that convention is written down:

\[\begin{split}V_n &= \cos\chi\, V_x + \sin\chi\, V_r \\ V_s &= -\sin\chi\, V_x + \cos\chi\, V_r\end{split}\]

rot_from is the transpose of rot_to – the rotation is orthogonal – so it undoes exactly this and nothing has to be re-derived to invert it.

Parameters:: chi (float or Array) – Angle [rad] of the frame axis from \(+x\), any shape.
Returns:: rot_to, rot_from – rot_to turns \((V_x, V_r)\) into \((V_n, V_s)\); rot_from turns \((V_n, V_s)\) back into \((V_x, V_r)\).
Return type:: Array, shape chi.shape + (2, 2)

util.matvec(A, b, out=None)[source]¶

Matrix-vector multiplication using einsum over trailing dimensions.

Parameters:

A (Array, shape (..., n, m)) – Input matrices with matrix dimensions in the last two axes.
b (Array, shape (..., m)) – Input vectors with vector dimension in the last axis.
out (Array, shape (..., n), optional) – Preallocated array to write the result into, sparing a new allocation.

Returns:

Result of matrix-vector multiplication A @ b for each corresponding matrix and vector in the trailing dimensions. Uses f32 precision and Fortran ordering for optimal performance.

Return type:

Array, shape (…, n)

Notes

When out is given, A and b are checked against two specific hot-path shapes – the per-grid-point 5x5 Jacobian contractions the solver runs every iteration – and dispatched to a specialised Fortran kernel if they match. Any other shape combination falls back to numpy.matmul(..., out=...), which is still correct but without the specialised kernel. Either way out is updated in place and returned by identity.

The two fast-path shapes, with b.shape[-1] fixed at 5 in every case:

Single batch axis: b.shape == (N, 5) and A.shape == (N, 5, 5) – one matrix per vector, no broadcasting.
Three-dimensional grid with one broadcast: b.shape === (ni, nj, nk, 5), and A.shape one of:
- (ni, nj, 1, 5, 5) – one matrix per \((i,j)\), shared across \(k\)
- (1, 1, nk, 5, 5) – one matrix per \(k\), shared across \((i,j)\)
- (ni, 1, 1, 5, 5) – one matrix per \(i\), shared across \((j,k)\)
A parallel set of kernels handles the same three broadcast patterns for b.shape[-1] == 2 instead of 5. Any other broadcast pattern within the 3-grid-axis case, including a full A shape (ni, nj, nk, n, m) with no broadcast axis at all, falls back to numpy.matmul.

Examples

>>> # Stack of 3x3 matrices with 3-element vectors
>>> A = np.random.randn(10, 5, 3, 3).astype(np.float32, order='F')
>>> b = np.random.randn(10, 5, 3).astype(np.float32, order='F')
>>> y = matvec(A, b)  # Shape: (10, 5, 3)

>>> # Single matrix-vector multiplication
>>> A = np.eye(3, dtype=np.float32, order='F')
>>> b = np.array([1, 2, 3], dtype=np.float32, order='F')
>>> y = matvec(A, b)  # y = b

>>> # Write into a preallocated buffer instead of allocating
>>> A = np.random.randn(4, 5, 5).astype(np.float32, order='F')
>>> b = np.random.randn(4, 5).astype(np.float32, order='F')
>>> out = np.empty((4, 5), dtype=np.float32, order='F')
>>> matvec(A, b, out=out) is out
True

util.resample(factor, vector, i_crit=None)[source]¶

Resampled a vector with specified factor, optionally preserving critical points.

Creates a new vector by resampling with a given factor. The new length equals len(vector) * factor (approximately). Critical indices are preserved, and for upsampling, fractional indices are linearly interpolated.

Parameters:

factor (float) – Resampling factor. Values > 1 increase resolution, < 1 decrease resolution. Must be > 0.
vector (array_like) – Input vector to resample.
i_crit (array_like, optional) – Critical indices that must be preserved in the output. Must be sorted and within [0, len(vector)-1]. If None, uses [0, len(vector)-1].

Returns:

resampled_values (Array) – Resampled vector with length approximately len(vector) * factor, with values at critical indices preserved and interpolated elsewhere.
crit_mapping (dict) – Dictionary mapping old critical indices to their positions in the new vector. Format: {old_index: new_index}

Examples

>>> # Increase resolution by factor of 2 with linear interpolation
>>> x = np.array([0, 1, 4, 9, 16])
>>> x_resampled, mapping = resample(2.0, x)
>>> len(x_resampled)  # approximately 10
>>> mapping[0], mapping[4]  # endpoints preserved at new positions

>>> # Decrease resolution preserving critical points
>>> x_resampled, mapping = resample(0.5, x, [0, 2, 4])
>>> len(x_resampled)  # approximately 2-3
>>> mapping  # {0: 0, 2: 1, 4: 2} - critical indices mapped to new positions

util.apply_perm_flip(array, perm, flip=())[source]¶

Apply permutation and flipping to array while preserving coordinate dimension.

This function applies a permutation and optional flipping to the spatial dimensions of an array, while preserving the last dimension (typically coordinates). This is commonly used in patch operations for coordinate transformations.

Parameters:

array (Array) – Input array to transform, typically with shape (…, 3) where the last dimension contains coordinate components.
perm (tuple) – Permutation to apply to spatial dimensions (0, 1, 2).
flip (tuple, optional) – Dimensions to flip after permutation. Default is () (no flipping).

Returns:

Transformed array with permutation and flipping applied to spatial dimensions, last dimension preserved.

Return type:

Array

Examples

>>> coords = np.random.randn(5, 4, 3).astype(np.float32, order='F')
>>> # Swap i and j dimensions
>>> transformed = apply_perm_flip(coords, perm=(1, 0, 2))
>>> transformed.shape == (4, 5, 3)  # i,j swapped, coordinates preserved
True

>>> # Swap dimensions and flip along first dimension
>>> transformed = apply_perm_flip(coords, perm=(1, 0, 2), flip=(0,))

util.cosine_cluster(n)[source]¶

Generate cosine-clustered points from 0 to 1, dense at both ends.

Produces n points on the unit interval following a half-cosine distribution, which clusters nodes near both endpoints (zeta = 0 and zeta = 1) and spaces them most coarsely in the middle:

\[\zeta_k = \tfrac{1}{2}\left(1 - \cos\frac{\pi k}{n - 1}\right), \quad k = 0, \ldots, n - 1.\]

The result is symmetric about 0.5 with exact endpoints 0 and 1.

Parameters:: n (int) – Number of points to generate (must be >= 2).
Returns:: Cosine-clustered vector from 0 to 1, dtype=float32.
Return type:: Array, shape (n,)

Examples

>>> z = cosine_cluster(5)
>>> z[0], z[-1]  # exact endpoints
(0.0, 1.0)
>>> np.allclose(z + z[::-1], 1.0)  # symmetric
True

util.cluster(ni, ER, dmax)[source]¶

Generate geometrically spaced points from 0 to 1 with expansion ratio and max spacing.

Creates a vector from 0 to 1 with geometrically spaced points, where the spacing expands at ratio ER but is capped to maximum size dmax. The algorithm creates an initial geometric dx vector, then iteratively adjusts the scaling to achieve unit total length while respecting dmax after capping.

Parameters:

ni (int) – Number of points to generate (must be >= 2)
ER (float) – Expansion ratio for geometric spacing (must be > 0)
dmax (float) – Maximum allowed spacing between consecutive points (must be > 0)

Returns:

Vector from 0 to 1 with clustered spacing, dtype=float32

Return type:

Array, shape (ni,)

Examples

>>> # Basic clustering with expansion ratio 1.2
>>> x = cluster(10, 1.2, 0.5)
>>> x[0], x[-1]  # Should be (0.0, 1.0)

>>> # Uniform spacing when ER=1.0
>>> x = cluster(5, 1.0, 1.0)
>>> np.allclose(x, np.linspace(0, 1, 5))
True

util.cluster_symmetric(n, ER, dmax=1.0)[source]¶

Generate geometrically spaced points from 0 to 1, dense at both ends.

Where cluster() expands away from a single end, this mirrors a half-width cluster() vector about the centreline, so the spacing grows at expansion ratio ER away from both endpoints and is coarsest in the middle. Unlike cosine_cluster(), which is also symmetric, the growth rate is controlled rather than fixed by the distribution.

Parameters:

n (int) – Number of points to generate. Must be odd and >= 3, so that the two mirrored halves share their midpoint.
ER (float) – Expansion ratio for geometric spacing (must be > 0).
dmax (float) – Maximum allowed spacing in the returned vector (must be > 0).

Returns:

Vector from 0 to 1 with spacing clustered at both ends, dtype=float32.

Return type:

Array, shape (n,)

Examples

>>> z = cluster_symmetric(9, 1.2, 1.0)
>>> z[0], z[-1]  # exact endpoints
(0.0, 1.0)
>>> np.allclose(z + z[::-1], 1.0)  # symmetric about 0.5
True

>>> # Uniform spacing when ER=1.0
>>> np.allclose(cluster_symmetric(5, 1.0, 1.0), np.linspace(0, 1, 5))
True

ember

Navigation

Versions

Related Topics

Array utilities¶

Flow angle and coordinate conversions¶

Vector and matrix operations¶

Grid construction¶

Array allocation and buffers¶

Miscellaneous geometry¶