dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Sep 23, 2008 at 11:03:16AM -0400, Shawn Walker wrote:
>
>
> On Tue, 23 Sep 2008, Anders Logg wrote:
>
>> 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)
>
> I'm sorry, should this be:
>
> f.vector() = x?

No, it should be x = f.vector(). The Function now owns the vector so
to get hold of the vector, one needs to ask for it.

I imagine that f.vector() = x works in C++ (it will assign the values
of x to the values in the vector of f), but not in Python.

-- 
Anders

Attachment: signature.asc
Description: Digital signature