dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Friday 03 July 2009 10:43:17 Garth N. Wells wrote:
> Johan Hake wrote:
> > On Thursday 02 July 2009 23:22:35 Garth N. Wells wrote:
> >> Anders Logg wrote:
> >>> On Thu, Jul 02, 2009 at 01:42:41PM +0200, Johan Hake wrote:
> >>>> On Thursday 02 July 2009 13:24:28 Garth N. Wells wrote:
> >>>>> Johan Hake wrote:
> >>>>>> On Thursday 02 July 2009 13:07:47 Garth N. Wells wrote:
> >>>>>>> Marie Rognes wrote:
> >>>>>>>> Garth N. Wells wrote:
> >>>>>>>>> Marie Rognes wrote:
> >>>>>>>>>> The following code gives r = 0.0. It is not supposed to be.
> >>>>>>>>>>
> >>>>>>>>>> The problem seems to be that f's vector is still all zeros at
> >>>>>>>>>> the call to interpolate. Could this be easily fixed?
> >>>>>>>>>
> >>>>>>>>> This example should have led to an error message since f is not a
> >>>>>>>>> discrete function. I'll take a look.
> >>>>>>>>
> >>>>>>>> Ok, thanks!
> >>>>>>>>
> >>>>>>>> However,
> >>>>>>>>
> >>>>>>>> (a) Why is f not a discrete function? (It is defined on a finite
> >>>>>>>> element space?)
> >>>>>>>
> >>>>>>> On second thought, it may be a discrete function. I think that this
> >>>>>>> is defined in the Python interface and not the C++ interface, so
> >>>>>>> I'll take a look.
> >>>>>>
> >>>>>> A user defined function is not a discrete function untill you either
> >>>>>> call interpolate() or vector, also in python. The problem with the
> >>>>>> later is that you then create a vector which is initialized to 0.
> >>>>>>
> >>>>>> I think this has been discussed before, but should we populate the
> >>>>>> vector using f.interpolate() when vector is called on a userdefined
> >>>>>> function?
> >>>>>
> >>>>> Or perhaps Function::vector() should throw an error if the vector has
> >>>>> not already been allocated.
> >>>>
> >>>> I vote for this.
> >>>>
> >>>> The error message can include information about the user might want to
> >>>> call interpolate?
> >>>>
> >>>> Johan
> >>>
> >>> Sounds good.
> >>
> >> Unfortunately this won't work because we often do
> >>
> >>      Function u(V);
> >>      solve(A, u.vector(), b);
> >
> > I see.
>
> The best short-term solution is a clear comment above Function::vector()
> in Function.h that calling this function will create a new and empty
> vector, thereby  and that it is the users responsibility to deal with
> this situation
>
> The real problem is a deeper design issue. For example, say a user
> creates an object
>
>      MyFunction my_function(V);
>
> and provides an eval function, and then calls vector(). Now,
> if the user calls eval(....) on my_function the eval function will be
> entered. However, if the user does
>
>      Function& function = my_function;
>
> (or just passes my_function to a function which expects a Function
> object) and then calls eval, the entries in the dof vector (which are
> zero) will be used to evaluate the function rather than the eval
> function. The behaviour of the object is inconsistent due to the
> call-back function eval(..) inside Function. The root of the problem is
> the old one of aFfunction object that can change state between
> user-defined and discrete.
>
> I don't see a simple solution. We could add a switch to indicate that a
> Function is user-defined, in which case a vector can never be created
> and a Function cannot change state,
>
>     Function f(V, Function::user);
>     Function pi_f0 = f;   // Interpolate f in V
>     Function pi_f1(W);
>     pi_f1.interpolate(f); // Interpolate f (defined on V) in W
>
>     GenericVector& x = f.vect();  // Throw an error
>
> Or, we add a new class like 'UserFunction' which is derived from
> (Generic)Function and cannot create a vector, and let Function create a
> vector upon construction. Now that we have interpolation, the need for
> allowing Functions to change state is diminished.

Yes I think this might be a solution. The problem is that in both cases we 
cannot prevent evil things to happen. A used can just derive his user defined 
class from Function instead of UserFunction, or he might forget to call the 
constructor of Function using the user flag.

Another solution might be to use a Function::discrete flag (Maybe 
Function::init_dofs is better?). Whenever this is used in the constructor 
init() is called. We throw an error when a user tries to call 
Function::vector(), when the vector is not initialized.

It will create more code when initializing the obligatory solution function:

  Function f(V,Function::discrete);

instead of just:

  Function f(V);

but it might be more intuitive and definitely more explicit. 

This is soooooo easy in python ;)

  f = Function(V)

has always an intialized vector.

Cheers!

Johan

> Garth
>
> > It would have worked in python though, as u would have been a discrete
> > Function when created. However this is now implemented as a call to
> > vector() in the constructor of DiscreteFunction, so that would also need
> > some attention.
> >
> > Johan
> >
> >> Garth
> >>
> >>> Just to check: this only occurs (in Python) when a user defines a
> >>> Function using a C++ expression or overloads eval(), right?
> >>>
> >>>
> >>>
> >>> -----------------------------------------------------------------------
> >>>-
> >>>
> >>> _______________________________________________
> >>> DOLFIN-dev mailing list
> >>> DOLFIN-dev@xxxxxxxxxx
> >>> http://www.fenics.org/mailman/listinfo/dolfin-dev
> >>
> >> _______________________________________________
> >> DOLFIN-dev mailing list
> >> DOLFIN-dev@xxxxxxxxxx
> >> http://www.fenics.org/mailman/listinfo/dolfin-dev