Poisson problem on a disk#

A Dirichlet problem on the unit disk, which is not a polygon and has no grid line along its boundary, solved with the Ferreira and Grigorieff (2006) scheme on a tensor-product grid. The disk is given as a level set \(\{\varphi < 0\}\); mesh adds grid lines until every line that meets the circle has a node on it, and the nodes on the circle are the boundary of the discrete domain.

Problem#

with the manufactured solution \(u_{\text{exact}} = \sin(2x + \tfrac12)\cos(\tfrac32 y - \tfrac15) + e^{0.6x - 0.4y}\). Its Laplacian worked out by hand gives \(f = \tfrac{25}{4} \sin(2x + \tfrac12)\cos(\tfrac32 y - \tfrac15) - 0.52\, e^{0.6x - 0.4y}\), so the error is checkable.

The disk and its grid#

levelset(φ) wraps a function of the point (x, y) that is negative inside, zero on the curve and positive outside. Giving it to the domain as the marker :D gives the mesh three labels, :D (interior), :D_boundary and :D_exterior, as for a polygon. The box around the disk must contain the circle strictly.

using Bramble
using Bramble: polygon_spaces, lift # public but unexported
using Random

φ(p) = hypot(p[1], p[2]) - 1
Ω = domain(interval(-1.25, 1.25) × interval(-1.25, 1.25), :D => levelset(φ))

Domain {2D, Float64}:

  Set:
    x: [-1.25, 1.25]
    y: [-1.25, 1.25]

  Markers:
    1 marker (1 levelset)
    Labels: :D, :D_boundary, :D_exterior

(n, n) is the number of points requested per axis, and they are seeds. For a disk the crossings of the grid lines with the circle close into finitely many lines, so mesh keeps every seed line and adds the crossings of each line with the circle as lines of the other direction until no line is new. The grid has more lines than requested, about \(1.8n\) for the uniform seeds that (true, true) asks for:

Ωₕ = mesh(Ω, (20, 20), (true, true))
xs, ys = collect(points(Ωₕ(1))), collect(points(Ωₕ(2)))
length(xs), length(ys)

(36, 36)

The boundary nodes lie on the circle to round-off, which is what the scheme asks of a domain (the boundary is made of grid nodes), and every node has one of the three labels:

m = markers(Ωₕ)
on_circle = [abs(φ((xs[i], ys[j]))) for i in eachindex(xs), j in eachindex(ys)]
max_off = maximum(on_circle[reshape(BitVector(m[:D_boundary]), length(xs), length(ys))])
sum(count.(BitVector.((m[:D], m[:D_boundary], m[:D_exterior])))) == length(xs) * length(ys)

true

max_off < 1.0e-12

true

Solving#

polygon_spaces(Ωₕ, :D, g) returns the trial space U and the test space W. Their unknowns are the interior nodes, and the Dirichlet data g sit at the boundary nodes. The load is avgₕ, the mean of \(f\) over each dual cell, and lift(U, u) puts a solution back on the whole grid.

u_exact(p) = sin(2p[1] + 0.5) * cos(1.5p[2] - 0.2) + exp(0.6p[1] - 0.4p[2])
f(p) = 6.25 * sin(2p[1] + 0.5) * cos(1.5p[2] - 0.2) - 0.52 * exp(0.6p[1] - 0.4p[2])

function solve_on(Ωₕ)
    U, W = polygon_spaces(Ωₕ, :D, u_exact)
    fₕ = element(gridspace(Ωₕ))
    avgₕ!(fₕ, f)
    A, F = assemble(form(U, W, (u, v) -> innerₕ(-Δₕ(u), v)), form(W, v -> innerₕ(fₕ, v)))
    return U, A \ F
end

The solution on the grid requested at 80 points per axis. Nodes outside the disk carry no value and are left blank; drag the plot to tilt it. The red curve is the polygon inscribed in the circle, whose vertices are the boundary nodes:

Ω₈₀ = mesh(Ω, (80, 80), (true, true))
uₕ = element(gridspace(Ω₈₀))
uₕ.data .= lift(solve_on(Ω₈₀)...).data
uₕ.data[BitVector(markers(Ω₈₀)[:D_exterior])] .= NaN
x₈₀, y₈₀ = collect(points(Ω₈₀(1))), collect(points(Ω₈₀(2)))
on_∂D = reshape(BitVector(markers(Ω₈₀)[:D_boundary]), length(x₈₀), length(y₈₀))
inscribed = sort([(x₈₀[i], y₈₀[j]) for i in eachindex(x₈₀), j in eachindex(y₈₀) if on_∂D[i, j]];
    by = v -> atan(v[2], v[1]))

The grid closure#

A coarse grid of the same disk, eight seeds per axis, keeps the lines readable. Each line that meets the circle has a node there, and those nodes are the vertices of the polygon inscribed in the circle:

Ωc = mesh(Ω, (8, 8), (true, true))
xc, yc = collect(points(Ωc(1))), collect(points(Ωc(2)))
length(xc), length(yc)

(14, 14)

Error#

norminf is the largest absolute error over the interior nodes, and snorm₁ₕ the discrete \(H^1\) seminorm of the error. h is the largest spacing between lines inside the disk's extent.

function errors(Ωₕ, U, u)
    xs, ys = collect(points(Ωₕ(1))), collect(points(Ωₕ(2)))
    nx, ny = length(xs), length(ys)
    e = element(U, u) - Rₕ(U, u_exact)
    E = reshape(lift(U, u).data, nx, ny) .- [u_exact((x, y)) for x in xs, y in ys]
    spacing(z) = maximum(diff(filter(t -> -1 <= t <= 1, z)))
    return (h = max(spacing(xs), spacing(ys)), max = norminf(e), h1 = snorm₁ₕ(e), E = E,
        boundary = reshape(BitVector(markers(Ωₕ)[:D_boundary]), nx, ny))
end

Convergence#

Four grids, requested at 20, 40, 80 and 160 points per axis. The order is the slope of the least-squares line through \(\log e\) against \(\log h\):

levels = (20, 40, 80, 160)
runs = map(levels) do n
    Ωₙ = mesh(Ω, (n, n), (true, true))
    return errors(Ωₙ, solve_on(Ωₙ)...)
end
h = [r.h for r in runs]
emax = [r.max for r in runs]
eh1 = [r.h1 for r in runs]
slope(e) = ([ones(length(h)) log.(h)] \ log.(e))[2]
println("order: max ", round(slope(emax); digits = 2), ", H¹ₕ ", round(slope(eh1); digits = 2))

order: max 2.04, H¹ₕ 1.97

Both orders are about two, as for a polygon: the errors fall by about a factor of four per halving of h, with the max error scattering more from level to level. The error is zero at the boundary nodes, the vertices of the inscribed polygon, where the data sit.

The errors against h on a log-log scale, with a slope-two reference line:

Seeds drawn at random#

The seeds need not be uniform. With (false, false) the interior points are drawn at random, and mesh closes the grid on them just the same, so the lines stay at the draws and the circle's crossings are added. The error of one draw scatters with it, so no order is read from a single one:

Random.seed!(1)
Ωr = mesh(Ω, (80, 80), (false, false))
rr = errors(Ωr, solve_on(Ωr)...)
println("random seeds, n = 80: max error ", round(rr.max; sigdigits = 3), ", H¹ₕ error ",
    round(rr.h1; sigdigits = 3))

random seeds, n = 80: max error 5.53e-5, H¹ₕ error 0.000336

Where the method stops#

A disk is the simplest curve the grid fits. When the crossings do not close, the lines are fitted anew, and some curves cannot be fitted at all: a curve that a horizontal or vertical line meets in more than two points (a star) throws an ArgumentError naming the line. The mesh internals page and the mesh docstring give the class of curves, the errors and the reasons.