Solution of the 2D Poisson problem using our own FEM code¶
In [1]:
from import_hack import *
from methodsnm.mesh_2d import *
from methodsnm.visualize import *
mesh = StructuredRectangleMesh(10, 10)
DrawMesh2D(mesh)
On the mesh, we define a finite element space:
In [2]:
from methodsnm.fes import *
fes = P1_Triangle_Space(mesh)
On the mesh, with the given finite element space, we define the variational formulation: $$ \int_{\Omega} \nabla u \cdot \nabla v dx + \int_{\Omega} u v dx = \int_{\Omega} f v dx $$ Replacing
- $u \leadsto u_h = \sum_j u_j \phi_j $ and
- $v \leadsto v_h = \phi_i$ we obtain a linear system:
In [3]:
from methodsnm.forms import *
from methodsnm.formint import *
from numpy import pi, cos
blf = BilinearForm(fes)
c = GlobalFunction(lambda x: 1, mesh = mesh)
blf += LaplaceIntegral(c)
blf += MassIntegral(c)
blf.assemble()
lf = LinearForm(fes)
f = GlobalFunction(lambda x: (1 + pi**2) * cos(pi*x[0]), mesh = mesh)
lf += SourceIntegral(f)
lf.assemble()
blf.matrix
In [4]:
print(blf.matrix, "\n", lf.vector)
Next, we solve the linear system. And we use a direct solver from scipy for that (for now):
In [5]:
uh = FEFunction(fes)
from scipy.sparse.linalg import spsolve
uh.vector = spsolve(blf.matrix, lf.vector)
Finally, we can visualize the solution:
In [6]:
DrawFunction2D(uh)
This problem now contained several simplifications, that can be removed step by step:
- qualitative evaluation only (no quantitative evaluation, no convergence study)
- boundary conditions are natural. Different boundary conditions make it more difficult
- 1D $\leadsto$ 2D
- P1 (low order) discretization vs. P2 and higher order discretizations
In [7]:
from methodsnm.forms import compute_difference_L2
uex = GlobalFunction(lambda x: cos(pi*x[0]), mesh = mesh)
l2diff = compute_difference_L2(uh, uex, mesh, intorder = 6)
print("l2diff =", l2diff)