Elliptic grid generation

Kernels for elliptic grid generation.

ember does not generate grids, but a mesher that smooths them by solving the Thompson-Thames-Mastin equations spends most of its time forming their stencils, and this is where the compiled code lives. The kernels take node numbers and coordinates only, so they serve any block structure a mesher assembles.

elliptic.ttm_stencil(X, stencil_g, stencil_s, PQ, pitch=0.0)[source]

Nine-point stencils of the Thompson-Thames-Mastin equations.

The equations in Thomas-Middlecoff form,

\[\alpha (r_{\xi\xi} + P r_\xi) - 2 \beta r_{\xi\eta} + \gamma (r_{\eta\eta} + Q r_\eta) = 0,\]

with \(\alpha = |r_\eta|^2\), \(\beta = r_\xi \cdot r_\eta\) and \(\gamma = |r_\xi|^2\), are frozen at the coordinates X and discretised by central differences, with the control function term lagged onto the right-hand side.

Parameters:

  • X (array, shape (n_node, 2)) – Node coordinates \((x, y)\).

  • stencil_g (int array, shape (n_reg, 9)) – Node numbers of each stencil, counted from zero, in the order centre, E, NE, N, NW, W, SW, S, SE.

  • stencil_s (int array, shape (n_reg, 9)) – Pitch shift of each stencil node in whole periods of y, for a stencil that straddles a periodic boundary.

  • PQ (array, shape (n_reg, 2)) – Control functions \(P, Q\) at each stencil centre.

  • pitch (float) – Period in y, or zero for none.

Returns:

  • w (array, shape (n_reg, 9)) – Stencil coefficients, normalised to one at the centre.

  • src (array, shape (n_reg, 2)) – Right-hand side of each row, so that sum(w[n] * X[stencil_g[n]]) == src[n] is the discrete equation, with the neighbours’ pitch shifts already moved into the y component.

elliptic.winslow_sweeps(X, stencil_g, stencil_s, pitch, irregular, slide, corner, wall, n_sweep)[source]

Point-Jacobi sweeps of Winslow’s equations, in single precision.

Each sweep moves every regular node to the solution of its nine-point stencil of

\[\alpha r_{\xi\xi} - 2 \beta r_{\xi\eta} + \gamma r_{\eta\eta} = 0,\]

frozen at the coordinates the sweep starts from, and every singular point to a weighted row of the same; then every sliding node in y to a weighted row of the coordinates after that; then draws every node near a corner towards it; then slides every wall node along a closed line to the foot of the perpendicular from the node outside it, no further than half way to either neighbour, its gaps to them kept above a floor.

Parameters:

  • X (array, shape (n_node, 2)) – Node coordinates \((x, y)\).

  • stencil_g (int arrays, shape (n_reg, 9)) – Node numbers of each regular node’s stencil, counted from zero, centre first then E, NE, N, NW, W, SW, S, SE, and the pitch shift of each.

  • stencil_s (int arrays, shape (n_reg, 9)) – Node numbers of each regular node’s stencil, counted from zero, centre first then E, NE, N, NW, W, SW, S, SE, and the pitch shift of each.

  • pitch (float) – Period in y, or zero for none.

  • irregular (tuple) – Singular points (nodes, indptr, indices, weights, const): each node to weights @ X[indices] over its row of the compressed sparse rows, plus const in y.

  • slide (tuple) – Sliding nodes (nodes, indptr, indices, weights, const), the same in y only.

  • corner (tuple) – (nodes, corners, dy, f): each node to c + (X[node] - c) * f, c its corner seen shifted dy in y.

  • wall (tuple) – (nodes, outside, dy, dy_outside, line, s_line, floor, seg, seg_out): the wall nodes in order round a closed ring, the last the first a turn on, and the node outside each, both seen shifted in y; the closed line (n_line, 2) they lie on, its last point its first, and arc lengths along it; the least cumulative arc length of each node past the first; and the samples of the line each node and the node outside it were nearest the sweep before.

  • n_sweep (int) – Number of sweeps.

Returns:

  • X (array, shape (n_node, 2)) – Swept coordinates, single precision.

  • seg, seg_out (int arrays) – Samples nearest each wall node and the node outside it, to start the next call from.

elliptic.polish_cost(X, ids, dy, la, lb, eps, power=1.0, w_size=1.0)[source]

Target-matrix cost of the corners of a set of cells, and its gradient.

At a corner with edges \(a\) along i and \(b\) along k, scaled by the cell’s target lengths into \(T = [a / l_a, b / l_b]\), the cost is the shape term \(\mu = |T|^2 / (2 \det T)\) raised to power, plus w_size times the size term \((\det T + 1 / \det T) / 2\), with \(\det T\) regularised by eps so that a folded corner costs a finite, steep amount – the target-matrix optimisation of Knupp, as a mesher’s polish minimises it.

Parameters:

  • X (array, shape (n_node, 2)) – Node coordinates \((x, y)\).

  • ids (int array, shape (n_cell, 4)) – Node numbers of each cell, counted from zero, in the order \((i, k)\), \((i + 1, k)\), \((i + 1, k + 1)\), \((i, k + 1)\).

  • dy (array, shape (n_cell, 4)) – Shift in y of each of a cell’s nodes, where it straddles a periodic boundary.

  • la (array, shape (n_cell,)) – Target lengths of each cell’s edges along i and k.

  • lb (array, shape (n_cell,)) – Target lengths of each cell’s edges along i and k.

  • eps (float) – Regularisation of the corner Jacobian.

  • power (float) – Power each corner’s shape term is raised to.

  • w_size (float) – Weight of the size term against the shape term.

Returns:

  • f (float) – The cost summed over every corner of every cell.

  • g (array, shape (n_node, 2)) – Its gradient with respect to X.