My first function - Follow up

As said in my previous post, I have been working on extending the implementation of interpolate to allow for an Expression as input. Currently it can also be used as in the Python dolfin interface, see here. Let's see how to use this new function in fem-fenics.

The Poisson equation

This example can be found in the FEniCS Book, it is the very first. The problem at hand is the Poisson equation with Dirichlet boundary conditions:

- Δu = f in Ω
u = u₀ on ∂Ω

We will solve this problem on the unit square, with f constant and equal to -6 and u₀ = 1 + x² + 2y². It can be verified that the exact solution is u_ex = 1 + x² + 2y². With the following ufl file:

element = FiniteElement("Lagrange", triangle, 1)

u = TrialFunction(element)

v = TestFunction(element)

f = Coefficient(element)

a = inner(grad(u), grad(v))*dx

L = f*v*dx

and Octave script:

pkg load fem-fenics msh

import_ufl_Problem ('Poisson')

# Create mesh and define function space

x = y = linspace (0, 1, 20);

mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4));

V = FunctionSpace('Poisson', mesh);

func = @(x,y) 1.0 + x^2 + 2*y^2;

# Define boundary condition

bc = DirichletBC(V, func, 1:4);

f = Constant ('f', -6.0);

# Define exact solution

u_e = Expression ('u_ex', func);

a = BilinearForm ('Poisson', V, V);

L = LinearForm ('Poisson', V, f);

# Compute solution

[A, b] = assemble_system (a, L, bc);

sol = A \ b;

u = Function ('u', V, sol);

# Save solution

save (u, 'poisson');

# Interpolate and save the exact solution

u_e_int = interpolate (u_e, V);

save (u_e_int, 'exact');

it is possible to compute the numerical solution, interpolate the analytical one on the same function space and then compare them. Using a visualisation tool like Paraview, one can verify that the computed solution and the interpolation of the exact one are practically the same. This is due to the fact that the Finite Elements Method with triangle elements on a rectangular domain can exactly represent a second order polynomial, as the solution of the problem at hand.

Here you can see a good solution poorly post-processed in Paraview to the Poisson problem solved in the example.