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 <[EMAIL PROTECTED]> 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?
Maybe I just misread the question.
- Shawn
_______________________________________________
DOLFIN-dev mailing list
[email protected]
http://www.fenics.org/mailman/listinfo/dolfin-dev