ufl binding

This week I started my work on the ufl function: it is now possible to write ufl code on-the-go, directly in your m-files. You can see below how the Poisson.ufl file of the homonymous example provided with fem-fenics (on the left) can be translated to a snippet of Octave code:

# Copyright (C) 2005-2009 Anders Logg element = FiniteElement("Lagrange", triangle, 1) u = TrialFunction(element) v = TestFunction(element) f = Coefficient(element) g = Coefficient(element) a = inner(grad(u), grad(v))*dx L = f*v*dx + g*v*ds

# Copyright (C) 2005-2009 Anders Logg ufl start Poisson ufl element = FiniteElement("Lagrange", triangle, 1) ufl ufl u = TrialFunction(element) ufl v = TestFunction(element) ufl f = Coefficient(element) ufl g = Coefficient(element) ufl ufl a = inner(grad(u), grad(v))*dx ufl L = f*v*dx + g*v*ds ufl end

How to use ufl

Basically, you just need to prepend what you would have written in your .ufl file with ufl. As you can see, anyway, there are also two new instructions. fem-fenics still needs to store your code in a separate file, which is then compiled using ffc, the FEniCS form compiler, but now ufl takes care of the process.

Your code should begin with the start command, and optionally with the name you want to assign to the file: in this example, we choose to open a new Poisson.ufl file. Be aware that ufl will not overwrite an existing file so, if you plan to use your script for several runs, my suggestion is to keep your working directory clean and tidy with a delete ('Poisson.ufl') after the snippet above.

When you are fine with your ufl code, the end command will tell ufl that it can compile and provide you with your freshly built problem. You can also specify options like BilinearForm (it is not the only one available, find a comprehensive list in the help message, in Octave), in case you wrote just part of the problem in your last lines.

What now?

A lot of commitment was devoted to this function. This is not due to intrinsic difficulties: a sketch of the function's code has been around for a while and the current implementation has not consistently slid away from it. The goal was to obtain a robust piece of code, since it will be the cornerstone of a new paradigm in fem-fenics usage. At least each and every example provided with the package needs to be modified to take advantage of this change, and this will be my next task.

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.

My first function

Lately I have coded my first function for fem-fenics: it is interpolate, which wraps the homonymous method of the dolfin::Function class. This allows to interpolate a Function, G, on the FunctionSpace of Function F, even if they are not defined on the same mesh, with a call like this:

res = interpolate (F, G)

I am working on extending it to allow for an Expression as input. With this function it is possible to make a quantitative comparison between the results of different discretisation approaches or to check the accuracy of a method, comparing the computed solution and an analytically obtained one.

The implementation

I provide here a comprehensive overview of the code for interpolate. First of all, the number of input arguments is obtained and an octave_value is declared, in order to hold the output. Then there is a check ensuring that exactly two arguments are provided and no more than one output value is asked for.

After verifying these preliminary conditions, there are some instructions checking that the function type is loaded and, if necessary, registering it. This way Octave is aware of it and can store it in an octave_value.

Eventually, the real computation is performed. After checking that the inputs are of the function type, with a static_cast the actual objects are extracted from the arguments:

const function & u0 = static_cast<const function&> (args(0).get_rep ());

const function & u1 = static_cast<const function&> (args(1).get_rep ());

Here comes the tricky part. The classes in fem-fenics are designed to hold constant dolfin objects, but dolfin::Function::interpolate is not a constant method. In order to be able to call it, a local dolfin::Function is constructed, used to perform the interpolation, then fed to the function constructor and assigned to the return value:

boost::shared_ptr<dolfin::Function> output (new dolfin::Function (u0.get_fun ()));

const dolfin::Function & input = u1.get_fun ();

output->interpolate (input);

std::string name = u1.get_str ();

retval = new function (name, output);