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
