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);