dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Friday 03 July 2009 15:24:31 Garth N. Wells wrote:
> Johan Hake wrote:
> > 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.
>
> I'm more inclined to default to discrete, and require a user to specify
> that something is user-defined. It would only appear in the class
> definition, e.g.
>
>    class MyFunction : Function
>    {
>       public:
>       MyFunction(FunctionSpace& V) : Function(V, Function::user) {}
>
>       void eval(....)
>
>    };
>
>    MyFunction f(V);

After some back and forth I might go for this too. However should we take it a 
bit further and always create a discrete function whenever a FunctionSpace is 
provided. The vector that is create is then always interpolated. This will 
work for:

   Function f(V);  (1)

if we just remove the error message last in

  Function::eval(double* values, const double* x)

we would be fine. It would then just be a dummy eval, which do nothing for 
Function with no vector.

During construction of f above we can then just initialize a vector with 
zeros. Interpolate it using the dummy eval, and we have a discrete 
and "interpolated" Function.

A Userdefined function intialized with a FunctionSpace is then also intialized 
as a discrete function, by interpolation. A user defined function constructed 
without a FunctionSpace will then by default be indiscretable.

In

  Function::vector

we throw an error when there are no vector to return.

With this design we define a Functions state during construction.

Also:

  Function f;

Would then just be what it looks like, a dummy function that do not do 
anything.

Sorry for the unstructured presentation, it is late...

> > This is soooooo easy in python ;)
> >
> >   f = Function(V)
> >
> > has always an intialized vector.
>
> I don't see how the problem is alleviated in PyDOLFIN. I could supply an
> eval function to a Python class MyFunction and pass the Function to
> DOLFIN. If a vector is created, then the vector is used to evaluate the
> function, but I can still call the provided eval function from the
> Python side.

This is correct. I was aiming at how we can fix the default state of a 
Function during construction.

Johan

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