Poisson-type problem on a polygon#

A Dirichlet problem on the polygon of Figure 1 in Ferreira and Grigorieff (2006): a non-convex nine-sided shape with oblique edges, discretised on a tensor-product grid whose lines are not fitted to it by hand. The coefficients vary, and the operator has a mixed term.

Problem#

-\partial_x(a\, \partial_x u + b\, \partial_y u) - \partial_y(b\, \partial_x u + c\, \partial_y u) = f \text{ in } P, \qquad u = u_{\text{exact}} \text{ on } \partial P,

with $$P$$ the polygon, $a = 2 + \tfrac12 \sin(x+y)$ , $c = 2 + \tfrac3{10} \cos x$ and the mixed coefficient $b = \tfrac12 + \tfrac15 \sin(xy)$ . The manufactured solution is $u_{\text{exact}} = \sin x \cos(0.7\, y) + 0.1\, xy$ , and $$f$$ is the left-hand side worked out by hand from it, so the error is checkable.

polygon(V) takes the vertices in order, either orientation, and the domain takes it as the marker :P.

using Bramble
using Bramble: polygon_spaces, lift # public but unexported
using Bramble: D₋ₓ, D₋ᵧ, D̃ₓ, D̃ᵧ, Dcₓ, Dcᵧ
using Random

V = [(0.0, 2.0), (2.0, 0.0), (4.0, 0.0), (6.0, 2.0), (6.0, 3.0), (4.0, 3.0), (4.0, 5.0),
    (2.0, 5.0), (0.0, 3.0)]
Ω = domain(interval(0.0, 6.0) × interval(0.0, 5.0), :P => polygon(V))

Domain {2D, Float64}:

  Set:
    x: [0.0, 6.0]
    y: [0.0, 5.0]

  Markers:
    1 marker (1 polygon)
    Labels: :P, :P_boundary, :P_exterior

The grid and the spaces#

The grid is random, requested at 72 by 60 points, which the closure enlarges. The seed makes it the same on every build, and the errors of a random grid scatter with the draw.

polygon_spaces(Ωₕ, :P, g) returns the trial space U and the test space W. Their unknowns are the interior nodes. The nodes on the polygon carry the Dirichlet data g exactly, and lift(U, u) puts a solution back on the whole grid.

Random.seed!(1)
Ωₕ = mesh(Ω, (72, 60), (false, false))
xs, ys = collect(points(Ωₕ(1))), collect(points(Ωₕ(2)))

sol(p) = sin(p[1]) * cos(0.7p[2]) + 0.1 * p[1] * p[2]
U, W = polygon_spaces(Ωₕ, :P, sol)

(GhostGridSpace{2D, Float64, 20404 dofs}, GhostGridSpace{2D, Float64, 20404 dofs})

Coefficients, operator and solution#

$$a$$ is sampled at the half points between neighbours in $$x$$ and $$c$$ between neighbours in $$y$$ , where the conservative form $-\tilde D_x(a\, D_{-x} u)$ reads them. $$b$$ is sampled at the grid points, and the mixed term is Dcₓ(bₕ * Dcᵧ(u)) + Dcᵧ(bₕ * Dcₓ(u)), central differences. The load is avgₕ, the mean of $$f$$ per dual cell.

a(p) = 2 + 0.5sin(p[1] + p[2])
b(p) = 0.5 + 0.2sin(p[1] * p[2])
c(p) = 2 + 0.3cos(p[1])

function load(p)
    x, y = p
    ux, uy = cos(x) * cos(0.7y) + 0.1y, -0.7sin(x) * sin(0.7y) + 0.1x
    uxx, uyy = -sin(x) * cos(0.7y), -0.49sin(x) * cos(0.7y)
    uxy = -0.7cos(x) * sin(0.7y) + 0.1
    ax, bx, by = 0.5cos(x + y), 0.2y * cos(x * y), 0.2x * cos(x * y)
    return -(ax * ux + a(p) * uxx) - c(p) * uyy - (bx * uy + by * ux + 2b(p) * uxy)
end

S = gridspace(Ωₕ)
aₕ, cₕ, fₕ = element(S), element(S), element(S)
nx, ny = length(xs), length(ys)
A3, C3 = reshape(aₕ.data, nx, ny), reshape(cₕ.data, nx, ny)
for l in 1:ny, j in 1:nx

    A3[j, l] = j > 1 ? a(((xs[j - 1] + xs[j]) / 2, ys[l])) : 0.0
    C3[j, l] = l > 1 ? c((xs[j], (ys[l - 1] + ys[l]) / 2)) : 0.0
end
bₕ = Rₕ(S, b)
avgₕ!(fₕ, load)

diffusion(u) = -D̃ₓ(aₕ * D₋ₓ(u)) - D̃ᵧ(cₕ * D₋ᵧ(u))
operator(u) = diffusion(u) - Dcₓ(bₕ * Dcᵧ(u)) - Dcᵧ(bₕ * Dcₓ(u))

The forms take U as trial and W as test space. The Dirichlet data are already in U, so assemble takes no dirichlet keyword:

A, F = assemble(form(U, W, (u, v) -> innerₕ(operator(u), v)), form(W, v -> innerₕ(fₕ, v)))
uh = lift(U, A \ F).data

The solution on the fine grid. Nodes outside the polygon carry no value and are left blank, and the red curve is the polygon; drag the plot to tilt it:

uₚ = element(S, ifelse.(BitVector(markers(Ωₕ)[:P_exterior]), NaN, uh))

The grid closure#

With the polygon as a marker, mesh closes the grid lines: every line that meets an edge of the polygon gets a node on it, so the boundary of $$P$$ is made of grid nodes, and the nodes fall into three labels, :P (interior), :P_boundary and :P_exterior. The closure may add lines, so the grid has more points than the requested count. The labels count the nodes of each kind:

m = markers(Ωₕ)
n_int, n_bnd, n_ext = count.(BitVector.((m[:P], m[:P_boundary], m[:P_exterior])))

(20404, 547, 22705)

n_int + n_bnd + n_ext == length(xs) * length(ys)

true

The figure draws a coarse grid of the same polygon, so its lines stay readable; the problem above is solved on a much finer one. Both are random, with the same seed:

Random.seed!(1)
Ωc = mesh(Ω, (8, 7), (false, false))
xc, yc = collect(points(Ωc(1))), collect(points(Ωc(2)))
length(xc), length(yc)

(22, 21)

(false, false) draws the interior lines at random, as on the other pages. The coarse grid after the closure, with the polygon over it:

The grid after the closure, with the polygon (red) over it.

Error#

The error is e, the difference between the solution and the exact one sampled at the nodes of U. norminf is its largest absolute value, and snorm₁ₕ its discrete $$H^1$$ seminorm.

e = element(U, A \ F) - Rₕ(U, sol)
err_max, err_h1 = norminf(e), snorm₁ₕ(e)
println("max error ", round(err_max; sigdigits = 3), ", H¹ₕ error ", round(err_h1; sigdigits = 3))

max error 0.00173, H¹ₕ error 0.00624

Both are small on this one grid. The errors of a single random grid scatter by a factor of two or three between draws, so no order is read from one: the suite's convergence test measures the orders over six levels. The error is zero at the boundary nodes, where the Dirichlet data sit.