Anders Logg wrote:
> 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.
>
Sounds good. It follows what was suggested in recent discussions on the
list - a function should be user-defined or discrete and stay that way.
> 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.
>
Yes, the parallel effort should be the top priority!
Would be good to add the above description to Blueprints. Things tend
get lost over time on the mailing list.
Garth
> --
> Anders
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> DOLFIN-dev mailing list
> [email protected]
> http://www.fenics.org/mailman/listinfo/dolfin-dev
_______________________________________________
DOLFIN-dev mailing list
[email protected]
http://www.fenics.org/mailman/listinfo/dolfin-dev
