# Redesign Function class (again)

I thought we finally reached a good design after last year's redesign

of the Function/

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.

## Blueprint information

- Status:
- Complete

- Approver:
- None

- Priority:
- Medium

- Drafter:
- None

- Direction:
- Needs approval

- Assignee:
- Anders Logg

- Definition:
- Approved

- Series goal:
- Accepted for 1.0.x

- Implementation:
- Implemented

- Milestone target:
- 0.9.4

- Started by
- Anders Logg

- Completed by
- Anders Logg

### Related branches

### Related bugs

### Sprints

### Whiteboard

See also http://

AL 2009-09-28: I have begun working on this. New classes NewCoefficient and Expression added.