Domains that are not rectangles

Download notebookView on GitHubRun when the site was built.

The base level can be several boxes of cells instead of the whole rectangle. Their union is the domain: an L, a T, a rectangle with a notch. The edges inside the bounding rectangle are a boundary of their own, "inner", which takes a condition like the faces do.

import numpy as np

from patchwork import AMRHierarchy, Dirichlet, Neumann, NodalPoisson, plot_patches

FACES = ("x_lower", "x_upper", "y_lower", "y_upper")

An L

base_boxes lists the boxes of base cells: here the left half of the square and the upper half of its right half. A patch must lie inside the domain, so a patch around the inner corner is made of two boxes side by side; a patch may lie over several base boxes.

def ell(base, levels):
    half = base // 2
    mesh = AMRHierarchy(
        base_shape=(base, base),
        bounds=((-1.0, 1.0), (-1.0, 1.0)),
        base_boxes=[((0, half), (0, base)), ((half, base), (half, base))],
    )
    for level in range(1, levels + 1):
        centre, width = base * 2**level // 2, base // 4
        mesh.add_patch(
            level=level,
            box=((centre - width, centre), (centre - width, centre + width)),
        )
        mesh.add_patch(
            level=level, box=((centre, centre + width), (centre, centre + width))
        )
    return mesh


fig, ax = plot_patches(ell(16, 2))

The corner of the L

The potential $r^{2/3}\sin(2\theta/3)$ is harmonic on the L and zero on its two inner edges, but its gradient is unbounded at the inner corner. A uniform mesh converges slowly there; nested patches at the corner do much better for the same number of unknowns. The inner edges are grounded, the default for "inner".

def corner(x, y):
    theta = np.arctan2(y, x) % (2 * np.pi)
    return np.hypot(x, y) ** (2 / 3) * np.sin(2 * theta / 3)


for base, levels in ((32, 0), (64, 0), (128, 0), (32, 2), (32, 4)):
    mesh = ell(base, levels)
    solver = NodalPoisson(mesh, faces={face: Dirichlet(corner) for face in FACES})
    nodes = solver.nodes
    result = solver.solve(np.zeros(nodes.num_nodes), potential="potential")
    on_coarse = np.all(np.abs(nodes.coordinates / (2 / 32) % 1) < 1e-9, axis=1)
    error = np.abs(result - corner(*nodes.coordinates.T))[on_coarse].max()
    print(
        f"{base:3d} x {base:<3d} cells, {levels} levels: {solver.free.size:6d} unknowns, error {error:.1e}"
    )

 32 x 32  cells, 0 levels:    705 unknowns, error 1.3e-02
 64 x 64  cells, 0 levels:   2945 unknowns, error 5.2e-03
128 x 128 cells, 0 levels:  12033 unknowns, error 1.9e-03
 32 x 32  cells, 2 levels:    961 unknowns, error 1.9e-03
 32 x 32  cells, 4 levels:   1217 unknowns, error 8.0e-04

fig, ax = plot_patches(mesh, field="potential", show_grid=False)

A T-shaped conductor

A bar with a stem, the foot of the stem at 1 V and the ends of the bar grounded. No current leaves through the other walls: "y_upper" and "inner" (the stem’s walls and the bar’s underside) take Neumann(0.0).

mesh = AMRHierarchy(
    base_shape=(48, 32),
    bounds=((-1.5, 1.5), (-1.0, 1.0)),
    base_boxes=[((0, 48), (16, 32)), ((16, 32), (0, 16))],
)
# at the two inner corners: in the bar across the corner, and in the stem beside it
for level in (1, 2):
    fine = 2**level
    for corner, inward in ((16, 1), (32, -1)):
        x, y = corner * fine, 16 * fine
        mesh.add_patch(level=level, box=((x - 8, x + 8), (y, y + 8)))
        stem = (x, x + 8) if inward == 1 else (x - 8, x)
        mesh.add_patch(level=level, box=(stem, (y - 8, y)))

solver = NodalPoisson(
    mesh,
    faces={
        "y_lower": Dirichlet(1.0),
        "x_lower": Dirichlet(0.0),
        "x_upper": Dirichlet(0.0),
        "y_upper": Neumann(0.0),
        "inner": Neumann(0.0),
    },
)
solver.solve(np.zeros(solver.nodes.num_nodes), potential="potential")
fig, ax = plot_patches(mesh, field="potential", show_grid=False)

On a uniform lattice (to_level, sample_patches) the nodes outside the domain are zero. The cell-centered solve_poisson takes such domains too, with its boundary value on the inner edges as well.