Poisson equation on an annulus

\[\begin{split} \mathbb{S}_u= \begin{cases} \Omega = \{(x, y)~:~R_{\text{inner}}^2 < x^2 + y^2 < R_{\text{outer}}^2\} \\ \partial\Omega_{\text{inner}} = \{(x, y)~:~ x^2 + y^2 = R_{\text{inner}}^2 \} \\ \partial\Omega_{\text{outer}} = \{(x, y)~:~ x^2 + y^2 = R_{\text{outer}}^2 \} \\ u_{\text{D}}\vert_{\partial\Omega_{\text{inner}}}=0 \\ u_{\text{D}}\vert_{\partial\Omega_{\text{outer}}}=1 \\ f=0 \\ u_{\text{e}}(r)=\frac{\ln(r/R_{\text{inner}})}{\ln(R_{\text{outer}}/R_{\text{inner}})}~~,~~r=\sqrt{x^2+y^2} \end{cases}\end{split}\]

import numpy as np
from dolfinx.fem import FunctionSpace
from lucifex.mesh import annulus_mesh, mesh_boundary, as_tri_mesh
from lucifex.fem import Function, Constant, as_tri_function
from lucifex.solver import bvp, BoundaryConditions
from lucifex.plt import plot_colormap, plot_mesh, save_figure
from lucifex.pde.poisson import poisson

Rinner = 1.0
Router = 2.0
Nradial = 10
dr = (Router - Rinner) / Nradial
mesh = annulus_mesh(dr, 'triangle')(Rinner, Router)
boundary = mesh_boundary(
    mesh, 
    {
        "inner": lambda x: x[0]**2 + x[1]**2 - Rinner**2,
        "outer": lambda x: x[0]**2 + x[1]**2 - Router**2,
    },
)

fs = FunctionSpace(mesh, ('P', 1))
f = Constant(mesh, 0.0)

bcs = BoundaryConditions(
    ("dirichlet", boundary['inner'], 0.0),
    ("dirichlet", boundary['outer'], 1.0),  
)
u = Function(fs, name='u')
u_solver = bvp(poisson, bcs)(u, f)
u_solver.solve()

fig, ax = plot_colormap(u, title='$u$')
plot_mesh(fig, ax, mesh, color='cyan', linewidth=0.5)
save_figure('u(x,y)_mesh', thumbnail=True)(fig)

def exact_solution(x, y, Rinner, Router):
    r = np.sqrt(x**2 + y**2)
    return np.log(r / Rinner) / np.log(Router / Rinner)

tri_mesh = as_tri_mesh(mesh)
ue = [exact_solution(xi, yi, Rinner, Router) for xi, yi in zip(tri_mesh.x_coordinates, tri_mesh.y_coordinates)]

fig, ax = plot_colormap(
    (tri_mesh.x_coordinates, tri_mesh.y_coordinates, ue), 
    title='$u_{\mathrm{e}}$', 
    cartesian=False, 
    triang=tri_mesh.triangulation,
)
save_figure('ue(x,y)')(fig)

u_tri = as_tri_function(u)
error = np.subtract(ue, u_tri.value)
fig, ax = plot_colormap(
    (tri_mesh.x_coordinates, tri_mesh.y_coordinates, error), 
    title='$u_{\mathrm{e}} - u$', 
    triang=tri_mesh.triangulation,
)
save_figure('error(x,y)')(fig)