dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Sep 23, 2008 at 12:35:48PM +0200, Marie Rognes wrote:
> Anders Logg wrote:
> > On Wed, Sep 17, 2008 at 01:23:17PM +0200, Evan Lezar wrote:
> >   
> >> On Tue, Sep 16, 2008 at 2:18 PM, Anders Logg <logg@xxxxxxxxx> wrote:
> >>
> >>     On Tue, Sep 16, 2008 at 01:10:35PM +0100, Garth N. Wells wrote:
> >>     >
> >>     >
> >>     > Anders Logg wrote:
> >>     > > On Mon, Sep 15, 2008 at 08:09:43PM +0100, Garth N. Wells wrote:
> >>     > >>
> >>     > >> Anders Logg wrote:
> >>     > >>> On Mon, Sep 15, 2008 at 06:31:56PM +0200, Anders Logg wrote:
> >>     > >>>> On Mon, Sep 15, 2008 at 04:42:32PM +0100, Garth N. Wells wrote:
> >>     > >>>>> Anders Logg wrote:
> >>     > >>>>>> On Mon, Sep 15, 2008 at 04:32:38PM +0100, Garth N. Wells wrote:
> >>     > >>>>>>> Anders Logg wrote:
> >>     > >>>>>>>> On Mon, Sep 15, 2008 at 03:47:59PM +0100, Garth N. Wells wrote:
> >>     > >>>>>>>>> Anders Logg wrote:
> >>     > >>>>>>>>>> On Mon, Sep 15, 2008 at 03:12:58PM +0100, Garth N. Wells
> >>     wrote:
> >>     > >>>>>>>>>>> Anders Logg wrote:
> >>     > >>>>>>>>>>>> On Mon, Sep 15, 2008 at 03:06:55PM +0200, Johan Hake wrote:
> >>     > >>>>>>>>>>>>> On Monday 15 September 2008 14:29:08 Garth N. Wells wrote:
> >>     > >>>>>>>>>>>>>> Could a Python expert take a look at site-packges/dolfin/
> >>     function.py?
> >>     > >>>>>>>>>>>>>> The code directly following the comment
> >>     > >>>>>>>>>>>>>>
> >>     > >>>>>>>>>>>>>>   # Special case, Function(element, mesh, x), need to
> >>     create simple form
> >>     > >>>>>>>>>>>>>> to get arguments
> >>     > >>>>>>>>>>>>>>
> >>     > >>>>>>>>>>>>>> need to be updated but I don't understand it well.
> >>     > >>>>>>>>>>>>> The first special case is for initializing a Function with
> >>     a given Vector, by
> >>     > >>>>>>>>>>>>> constructing a dofmap from the handed element.
> >>     > >>>>>>>>>>>>>
> >>     > >>>>>>>>>>>>> As constructing a Function from a vector is removed from
> >>     the cpp interface,
> >>     > >>>>>>>>>>>>> and we have not, (or have we?) figured out how to wrap a
> >>     shared_ptr in swig,
> >>     > >>>>>>>>>>>>> we should probably just remove the first case for now.
> >>     > >>>>>>>>>>>>>
> >>     > >>>>>>>>>>>>> Johan
> >>     > >>>>>>>>>>>> The question is how we want to create discrete Functions in
> >>     Python.
> >>     > >>>>>>>>>>>> Previously, this was done by
> >>     > >>>>>>>>>>>>
> >>     > >>>>>>>>>>>>   u = Function(element, mesh, Vector())
> >>     > >>>>>>>>>>>>
> >>     > >>>>>>>>>>>> but now the third argument is not needed anymore. If we
> >>     remove it,
> >>     > >>>>>>>>>>>> we get
> >>     > >>>>>>>>>>>>
> >>     > >>>>>>>>>>>>   u = Function(element, mesh)
> >>     > >>>>>>>>>>>>
> >>     > >>>>>>>>>>>> but that doesn't work since that is the way to initialize a
> >>     > >>>>>>>>>>>> user-defined function (something overloading eval()).
> >>     > >>>>>>>>>>>>
> >>     > >>>>>>>>>>>> We could put in a flag and make "discrete" the default. Then
> >>     all
> >>     > >>>>>>>>>>>> user-defined functions need to set the flag to "user".
> >>     > >>>>>>>>>>>>
> >>     > >>>>>>>>>>>> Suggestions? This is a good time to worry about how we want
> >>     to design
> >>     > >>>>>>>>>>>> the Function interface.
> >>     > >>>>>>>>>>>>
> >>     > >>>>>>>>>>> Sounds ok to me. This is basically what Vector() was doing,
> >>     and a flag
> >>     > >>>>>>>>>>> would be more descriptive.
> >>     > >>>>>>>>>>>
> >>     > >>>>>>>>>>> Garth
> >>     > >>>>>>>>>> Maybe we could first try to think seriously about reducing the
> >>     number
> >>     > >>>>>>>>>> of different constructors in Function. There are 14 now! See
> >>     below.
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>>> I guess we need the following two basic constructors (empty
> >>     and copy):
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>>>   /// Create empty function (read data from file)
> >>     > >>>>>>>>>>   Function();
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>>>   /// Copy constructor
> >>     > >>>>>>>>>>   Function(const Function& f);
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>>> Then we have one for reading from file, which seems ok:
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>>>   /// Create function from data file
> >>     > >>>>>>>>>>   explicit Function(const std::string filename);
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>>> And then the following set of constructors for constants:
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>>>   /// Create constant scalar function from given value
> >>     > >>>>>>>>>>   Function(Mesh& mesh, real value);
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>> This one is useful.
> >>     > >>>>>>>>>
> >>     > >>>>>>>>>>   /// Create constant vector function from given size and
> >>     value
> >>     > >>>>>>>>>>   Function(Mesh& mesh, uint size, real value);
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>> We could get rid of this one and use the below constructor.
> >>     > >>>>>>>>>
> >>     > >>>>>>>>>>   /// Create constant vector function from given size and
> >>     values
> >>     > >>>>>>>>>>   Function(Mesh& mesh, const Array<real>& values);
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>> This one is useful.
> >>     > >>>>>>>>>
> >>     > >>>>>>>>>>   /// Create constant tensor function from given shape and
> >>     values
> >>     > >>>>>>>>>>   Function(Mesh& mesh, const Array<uint>& shape, const Array
> >>     <real>& values);
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>> This is the most generic of the constant functions, so I guess
> >>     we need it.
> >>     > >>>>>>>>>
> >>     > >>>>>>>>>> And then there's this constructor which is needed for w.split
> >>     (u, p):
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>>>   /// Create discrete function from sub function
> >>     > >>>>>>>>>>   explicit Function(SubFunction sub_function);
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>>> But then there's the following mess of constructors:
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>> Some of these constructors are necessary to support the
> >>     PyDOLFIN
> >>     > >>>>>>>>> interface. Can we get around this somehow to avoid duplication?
> >>     > >>>>>>>>>
> >>     > >>>>>>>>>>   /// Create function from given ufc::function
> >>     > >>>>>>>>>>   Function(Mesh& mesh, const ufc::function& function, uint
> >>     size);
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>>>   /// Create discrete function for argument function i of form
> >>     > >>>>>>>>>>   Function(Mesh& mesh, Form& form, uint i = 1);
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>>>   /// Create discrete function for argument function i of form
> >>     > >>>>>>>>>>   Function(Mesh& mesh, DofMap& dof_map, const ufc::form& form,
> >>     uint i = 1);
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>>>   /// Create discrete function for argument function i of form
> >>     (data may be shared)
> >>     > >>>>>>>>>>   Function(std::tr1::shared_ptr<Mesh> mesh,
> >>     > >>>>>>>>>>            std::tr1::shared_ptr<GenericVector> x,
> >>     > >>>>>>>>>>            std::tr1::shared_ptr<DofMap> dof_map, const
> >>     ufc::form& form, uint i = 1);
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>>>   /// Create discrete function based on signatures
> >>     > >>>>>>>>>>   Function(std::tr1::shared_ptr<Mesh> mesh,
> >>     > >>>>>>>>>>            const std::string finite_element_signature,
> >>     > >>>>>>>>>>            const std::string dof_map_signature);
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>>>   /// Create user-defined function (evaluation operator must
> >>     be overloaded)
> >>     > >>>>>>>>>>   explicit Function(Mesh& mesh);
> >>     > >>>>>>>>>>
> >>     > >>>>>>>>> We need this one.
> >>     > >>>>>>>>>
> >>     > >>>>>>>>> Garth
> >>     > >>>>>>>> If we just consider discrete functions for a while, the question
> >>     is
> >>     > >>>>>>>> how these may be most conveniently (and naturally) defined in
> >>     C++ and
> >>     > >>>>>>>> Python.
> >>     > >>>>>>>>
> >>     > >>>>>>>> In C++, one only has a dolfin::Form, for example
> >>     PoissonBilinearForm,
> >>     > >>>>>>>> and then it's simple to create a discrete Function by
> >>     > >>>>>>>>
> >>     > >>>>>>>>   Function u(mesh, form);
> >>     > >>>>>>>>
> >>     > >>>>>>>> This will extract the element and dof map for the second
> >>     argument of
> >>     > >>>>>>>> the form (the trial function) which is normally what is needed.
> >>     > >>>>>>>>
> >>     > >>>>>>>> In Python, one does not have a dolfin::Form, but instead one has
> >>     a
> >>     > >>>>>>>> FiniteElement, and then the simplest thing to do is
> >>     > >>>>>>>>
> >>     > >>>>>>>>   u = Function(element, mesh)
> >>     > >>>>>>>>
> >>     > >>>>>>>> The element is the first argument for practical reasons (see
> >>     > >>>>>>>> function.py) but maybe it shouldn't. I'd like to change this so
> >>     that
> >>     > >>>>>>>> the mesh is always first. All Functions require a Mesh and then
> >>     it's
> >>     > >>>>>>>> natural to put this first.
> >>     > >>>>>>>>
> >>     > >>>>>>>> So then we would have
> >>     > >>>>>>>>
> >>     > >>>>>>>>   C++:    Function u(mesh, form);
> >>     > >>>>>>>>   Python: u = Function(mesh, element)
> >>     > >>>>>>>>
> >>     > >>>>>>>> On the other hand, we've been discussing adding a FunctionSpace
> >>     class,
> >>     > >>>>>>>> and then it might be natural to just have
> >>     > >>>>>>>>
> >>     > >>>>>>>>   C++:    Function u(V);
> >>     > >>>>>>>>   Python: u = Function(V)
> >>     > >>>>>>>>
> >>     > >>>>>>>> This would create a discrete Function. Constant Functions and
> >>     > >>>>>>>> user-defined Functions may be created without reference to a
> >>     > >>>>>>>> FunctionSpace. This would solve the problem of overloading
> >>     > >>>>>>>> constructors. It would be very clear that whenever a
> >>     FunctionSpace is
> >>     > >>>>>>>> involved, it is a discrete Function.
> >>     > >>>>>>>>
> >>     > >>>>>>> Agree. It's not appropriate to initialise a discrete function
> >>     with a
> >>     > >>>>>>> form. It seems that using a FunctionSpace will simplify the
> >>     interface
> >>     > >>>>>>> and provide uniformity across the C++ and Python interfaces, so
> >>     let's
> >>     > >>>>>>> get FunctionSpace (or something similar with another name) in
> >>     place and
> >>     > >>>>>>>   then remove some the Function constructors.
> >>     > >>>>>> I think FunctionSpace is a good name.
> >>     > >>>>>>
> >>     > >>>>>>> Function spaces are not only associated with DiscreteFunctions.
> >>     We
> >>     > >>>>>>> usually interpolate user defined function in the finite element
> >>     space,
> >>     > >>>>>>> so perhaps there is some scope to unify discrete and user defined
> >>     > >>>>>>> functions?
> >>     > >>>>>>>
> >>     > >>>>>>> Garth
> >>     > >>>>>> Perhaps, but it would require storing an extra vector of values
> >>     for
> >>     > >>>>>> user-defined functions unnecessarily.
> >>     > >>>>>>
> >>     > >>>>> I'm not suggesting that we store a vector of values - just pointing
> >>     out
> >>     > >>>>> a function space is also associated with user-defined functions.
> >>     > >>>>>
> >>     > >>>>> Garth
> >>     > >>>> Yes, I forgot for a minute that user-defined functions also need to
> >>     be
> >>     > >>>> associated with a particular function space when used in a form.
> >>     > >>>>
> >>     > >>>> Then it seems we still have the problem of differentiating between
> >>     > >>>> user-defined functions and discrete functions when overloading
> >>     > >>>> constructors.
> >>     > >>> I have a suggestion that would solve this as well as another problem
> >>     > >>> we've discussed a few times before regarding Functions
> >>     > >>> (time-stepping). It's rather drastic but it might work.
> >>     > >>>
> >>     > >>> * Remove GenericFunction, DiscreteFunction, UserFunction etc and just
> >>     > >>>   have one single class named Function.
> >>     > >>>
> >>     > >> Sounds good to me. I often found the current design a bit complicated.
> >>     A
> >>     > >> simpler design will make parallisation easier.
> >>     > >>
> >>     > >>> * Always require a FunctionSpace in the inialization of Function:
> >>     > >>>
> >>     > >>>   u = Function(V)
> >>     > >>>
> >>     > >>> * Let a Function consist of two things, a function space V and
> >>     > >>>   degrees of freedom U:
> >>     > >>>
> >>     > >>>   u_h = \sum_i U_i \phi_i
> >>     > >>>
> >>     > >>>   where {\phi_i} are the basis functions for V.
> >>     > >>>
> >>     > >> As long as it's general enough to permit 'quadrature functions'.
> >>     > >>
> >>     > >>> * Keep a shared_ptr to V, and a pointer to a Vector U.
> >>     > >>>
> >>     > >> It might be useful to make the vector a shared pointer (need to check
> >>     > >> what a smart pointer returns if it doesn't point to anything) because
> >>     it
> >>     > >> would be useful with sub-functions.
> >>     > >>
> >>     > >>> * The Vector pointer is initialized to 0 by default and as long as it
> >>     > >>>   remains 0, the Function behaves like a user-defined Function, that
> >>     > >>>   is, the eval() function is used.
> >>     > >>>
> >>     > >>> * Whenever someone calls the vector() member function, the Vector U
> >>     is
> >>     > >>>   initialized to a Vector of correct size (determined by the DofMap
> >>     in
> >>     > >>>   the FunctionSpace V). From this point on, the Function behaves like
> >>     > >>>   a discrete Function, that is, the values in U are used, not eval().
> >>     > >>>
> >>     > >> Sounds good.
> >>     > >>
> >>     > >>> * This would make the behavior dynamic. For example, one may set u0
> >>     to
> >>     > >>>   an initial value in time-stepping and then do u0 = u1 and u0 would
> >>     > >>>   change behavior from user-defined to discrete.
> >>     > >>>
> >>     > >>> * Constant functions are handled by a new class named simply
> >>     Constant,
> >>     > >>>   which is a subclass of Function that overloads eval() and returns
> >>     > >>>   the constant value.
> >>     > >>>
> >>     > >> Perhaps ConstantFunction would be a better name?
> >>     > >
> >>     > > Maybe, but it would also be natural to couple this to Constant in UFL
> >>     > > in the same way as we couple Function to Function.
> >>     > >
> >>     > > Another reason to use Constant is that ConstantFunction implies that
> >>     > > there may be other subclasses of Function than Constant and it would
> >>     > > be good to avoid building a hierarchy (again). The class Constant just
> >>     > > happens to be implemented as a subclass of Function, but it should be
> >>     > > thought of as something different.
> >>     > >
> >>     > >> No matter what happens, it seems that FunctionSpace is the next step.
> >>     > >> Related to this, is it the intention that a FiniteElement owns a
> >>     > >> ufc::finite_elememt and provides wrapper, like DofMap?
> >>     > >
> >>     > > Yes, that was my plan. It only needs to overload the subset of
> >>     > > functionality of ufc::finite_element that we need now.
> >>     > >
> >>     >
> >>     > OK. I can put FiniteElement together pretty quickly if you haven't
> >>     > started already.
> >>     >
> >>     > Garth
> >>
> >>     I haven't (not more than the empty template that's currently there).
> >>
> >>
> >>
> >> Hi
> >>
> >> What is the correct way to initialize a  discrete function in python.  I see
> >> that the way I was using in the electromagnetic demo no longer works.  f =
> >> Function(element, mesh, vector)
> >>
> >> Evan
> >>     
> >
> > That is the correct way, but the vector will be ignored since the
> > Function from now on will create and own its vector of degrees of
> > freedom.
> >
> > So you can do something like
> >
> >   f = Function(element, mesh, Vector())
> >   x = f.vector()
> >
> > The vector you send in will be ignored and then you need to pick it
> > out afterwards.
> >
> > Once the new design is in place it will be
> >
> >   f = Function(V) # V is a FunctionSpace
> >   x = f.vector()
> >
> >   
> 
> 
> Follow-up question: So, what is the right way to initialize a Function 
> with a vector of computed values?

You can't. You have to create the function first, then set the values
in the vector:

  f = Function(element, mesh, Vector())
  x = f.vector()

  solve(A, x, b)

-- 
Anders

Attachment: signature.asc
Description: Digital signature