Part 1: Base case

Consider the Dirichlet problem

\[ \begin{aligned} -\Delta u &= 4, & \quad &\text{in $\Omega$} \\ u &= 0, & &\text{on $\partial\Omega$} \end{aligned} \]

on the unit disk \(\Omega=\{(x,y)\;:\; x^2+y^2\le 1\}\)

The exact solution to this problem is

\[ \boxed{u = 1 - x^2 - y^2} \]

A function implementing quadratic Lagrange interpolation over a reference triangle is implemented in the file interpolate.py which is used by the function poisson2 from poisson.py to create a finite element that locally evaluates the following variational Poisson problem

\[ \sum_{j=1}^{n} \mathbf{u}_{j} \int_{\Omega} \nabla \phi_{j} \cdot \nabla \phi_{i}=\int_{\Omega} \phi_{i} f+\int_{\partial \Omega_{N}} \phi_{i} g_{N}-\sum_{j=n+1}^{n+n_{\partial}} \mathbf{u}_{j} \int_{\Omega} \nabla \phi_{j} \cdot \nabla \phi_{i} \]

A model is constructed by loading a pre-defined mesh from the file mesh.m228 as follows:

model =  anabel.MeshGroup.read(f"../dat/circle_iso/mesh.m228", cell="triangle6")
elem = poisson(*[lagrange_t6()]*3, f=f)
U = model.compose(elem=elem, solver="sparse", verbose=True)

This builds a function U which accepts quadrature locations and weights as input. Given a particular quadrature rule, a solution may be computed and plotted as follows:

quad = anon.quad.simplex.Simplex.load(f"../dat/quadrature/gauss{g:02}.m")
model.plot(U(quad.points, quad.weights))

Nonlinear Source

\[ \begin{aligned} -\Delta u &= \frac{\pi^2}{4}\left(\cos\frac{\pi r}{2} + \operatorname{sinc} \frac{\pi r}{2}\right), & \quad &\text{in $\Omega$} \\ u &= 0, & &\text{on $\partial\Omega$} \end{aligned} \]

where \(\Omega\) is again the unit disk. The exact solution is \(u(x,y)=\cos \frac{\pi r}{2}\), where \(r=\sqrt{x^2+y^2}\).

Contributions from Errors

\[ \iint_T \nabla(u_{FE} - u_{exact}) \cdot \nabla(u_{FE} - u_{exact})\,dx \, dy, \qquad (T\in\mathcal{T}) \]

The \(H^1\) semi-norm error in the finite element solution is the square root of the sum of all the errors shown here.

Appendix

Source Code of Interest