Redesign Function class (again)

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.

This redesign is probably something that should wait until after the
parallel effort is fully in place, but it would be good to spend some
effort on formulating a plan now.