I thought we finally reached a good design after last year's redesign
of the Function/FunctionSpace classes, but I think we need to redo it.
I had a long discussion with Marie and Johan yesterday and here are
some of our conclusions.
First a couple of motivating examples and then a simple and elegant
solution. :-)
Consider first the function g(x, y) = x*(1 - x)*y*(1 - y) on a 1 x 1 mesh
of the unit square. This function is zero along the boundary and in
particular in each vertex. Then consider the following code:
g = "x[0]*(1 - x[0])*x[1]*(1 - x[1])"
mesh = UnitSquare(1, 1)
V1 = FunctionSpace(mesh, "CG", 1)
V2 = FunctionSpace(mesh, "CG", 2)
# Create function f in V1
f1 = Function(V1, g)
# Evaluate f at x = (0.5, 0,5)
x = array((0.5, 0.5))
values = array((0.0,))
f1.eval(values, x)
# Project f to piecewise quadratics
f2 = project(f1)
Now, the question is, what is values[0] and what is f2?
One would think that values[0] = f2 = 0. This would be natural since
f1 is zero in each vertex and thus should be zero everywhere. But what
happens is that values[0] will be 1/16 and f2 will be a bubble
function on the unit square.
This is because f1 is not really a function in V1. It only behaves
that way when used as a coefficient in a form in which case it will be
interpolated on each cell.
Adding a call to f1.interpolate() in the above code fixes this.
So our suggestion is to make sure that f1 is really a function in V1
and at the same time make sure that:
1. A Function always has a FunctionSpace
2. A Function always has a Vector
We can do this by breaking out the functionality for user-defined
functions or expressions like g = "x[0]*(1 - x[0])*x[1]*(1 - x[1])"
into a separate class. Then we can have a clean and simple
implementation of the Function class free from all complex logic wrt
having or not having a vector, creating a vector when it doesn't exist
etc.
I suggest calling this new class Expression and keeping it very
simple. An Expression is something like "x[0]*(1 - x[0])*x[1]*(1 -
x[1])" or more specifically something having an eval operator. Both
Functions and Expressions can be used as coefficients in forms, so the
following design seems appropriate:
class Expression : public ufc::function, public Coefficient
{
public:
/// Function evaluation (simple version)
virtual void eval(double* values, const double* x) const;
/// Function evaluation (alternate version)
virtual void eval(double* values, const Data& data) const;
/// Implementation of ufc::function interface
void evaluate(double* values, const double* coordinates, const ufc::cell& cell) const;
/// Implementation of Coefficient interface
void interpolate(double* coefficients, const ufc::cell& ufc_cell, uint cell_index, int local_facet=-1) const;
};
class Function : public Coefficient
{
public:
...
/// Implementation of Coefficient interface
void interpolate(double* coefficients, const ufc::cell& ufc_cell, uint cell_index, int local_facet=-1) const;
};
We will need to change the name of the current Coefficient class and
introduce the new Coefficient base class.
I believe this will also take care of all earlier reported problems on
mysterious behavior of the Function class.
