dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Wed, Sep 02, 2009 at 10:03:57PM +0100, Garth N. Wells wrote:
>
>
> 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.

Yes.

> > 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!

Agree.

> Would be good to add the above description to Blueprints. Things tend
> get lost over time on the mailing list.

Just added. I really like Launchpad. Seems very useful.

--
Anders

Attachment: signature.asc
Description: Digital signature