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@fenics.org > >>>>> http://www.fenics.org/mailman/listinfo/dolfin-dev > >>>> > >>>> _______________________________________________ > >>>> DOLFIN-dev mailing list > >>>> DOLFIN-dev@fenics.org > >>>> http://www.fenics.org/mailman/listinfo/dolfin-dev _______________________________________________ DOLFIN-dev mailing list DOLFIN-dev@fenics.org http://www.fenics.org/mailman/listinfo/dolfin-dev