dolfin team mailing list archive

["Martin Sandve Alnæs" <martinal@xxxxxxxxx>]

2008/9/9 Anders Logg <logg@xxxxxxxxx>:
> On Tue, Sep 09, 2008 at 07:37:47PM +0200, Martin Sandve Alnæs wrote:
>> 2008/9/9 Anders Logg <logg@xxxxxxxxx>:
>> > On Tue, Sep 09, 2008 at 03:47:10PM +0100, Garth N. Wells wrote:
>> >>
>> >>
>> >> Anders Logg wrote:
>> >> > On Tue, Sep 09, 2008 at 04:31:41PM +0200, Martin Sandve Alnæs wrote:
>> >> >> 2008/9/9 Garth N. Wells <gnw20@xxxxxxxxx>:
>> >> >>>
>> >> >>> Anders Logg wrote:
>> >> >>>> On Tue, Sep 09, 2008 at 11:30:26AM +0100, Garth N. Wells wrote:
>> >> >>>>> Anders Logg wrote:
>> >> >>>>>> On Tue, Sep 09, 2008 at 09:06:00AM +0200, Anders Logg wrote:
>> >> >>>>>>> On Mon, Sep 08, 2008 at 11:53:09PM +0200, Johan Hake wrote:
>> >> >>>>>>>> On Monday 08 September 2008 21:45:27 Martin Sandve Alnæs wrote:
>> >> >>>>>>>>> 2008/9/8 Anders Logg <logg@xxxxxxxxx>:
>> >> >>>>>>>>>> On Mon, Sep 08, 2008 at 11:12:14AM +0200, Martin Sandve Alnæs wrote:
>> >> >>>>>>>>>>> 2008/9/8 Johan Hoffman <jhoffman@xxxxxxxxxx>:
>> >> >>>>>>>>>>>>> 2008/9/8 Dag Lindbo <dag@xxxxxxxxxx>:
>> >> >>>>>>>>>>>>>> Anders Logg wrote:
>> >> >>>>>>>>>>>>>>> There seems to be a problem (among many) with the current design of
>> >> >>>>>>>>>>>>>>> the Function classes (see thread "evaluating higher order mesh
>> >> >>>>>>>>>>>>>>> function").
>> >> >>>>>>>>>>>>>>>
>> >> >>>>>>>>>>>>>>> In particular, the finite element is missing in DiscreteFunction.
>> >> >>>>>>>>>>>>>>> My suggestion would be to just add it and let a DiscreteFunction
>> >> >>>>>>>>>>>>>>> consist of the following four items which are always available:
>> >> >>>>>>>>>>>>>>>
>> >> >>>>>>>>>>>>>>>   mesh, x, dof_map, finite_element
>> >> >>>>>>>>>>>>>>>
>> >> >>>>>>>>>>>>>>> Is this enough, and what other issues to we need to fix?
>> >> >>>>>>>>>>>>>> I'm not sure I agree that the dof map and finite element should be
>> >> >>>>>>>>>>>>>> owned by the discrete function. There was a great suggestion from
>> >> >>>>>>>>>>>>>> Martin, in a thread "Abstraction idea" from 06/05/2008, to create a
>> >> >>>>>>>>>>>>>> class FunctionSpace where the mesh, element and dof_map(s) are
>> >> >>>>>>>>>>>>>> aggregated. Citing Martin:
>> >> >>>>>>>>>>>>>> U = FunctionSpace(mesh, dofmapset, form, 0) # or something similar
>> >> >>>>>>>>>>>>>> u = Function(U)
>> >> >>>>>>>>>>>>>> v = Function(U)
>> >> >>>>>>>>>>>>>>
>> >> >>>>>>>>>>>>>> This seems a solid approach to me since it would provide a way of
>> >> >>>>>>>>>>>>>> encapsulating the mathematical formulation of the problem, which is
>> >> >>>>>>>>>>>>>> more or less const and likely to be reused by many discrete
>> >> >>>>>>>>>>>>>> functions in a solver.
>> >> >>>>>>>>>>>>>>
>> >> >>>>>>>>>>>>>> It seems to me that there is an obvious risk that a lot of redundant
>> >> >>>>>>>>>>>>>> initialization would occur if all discrete functions should own
>> >> >>>>>>>>>>>>>> their own elements and dof maps. There seems to be consensus that
>> >> >>>>>>>>>>>>>> the mesh should be "global" for efficiency reasons, so why not treat
>> >> >>>>>>>>>>>>>> the function space the same way?
>> >> >>>>>>>>>>>>>>
>> >> >>>>>>>>>>>>>> Is there a problem with an approach where the funciton _always_ owns
>> >> >>>>>>>>>>>>>> the vector and _never_ owns the function space (and mesh)? A very
>> >> >>>>>>>>>>>>>> strict design would avoid shared/smart pointers, provide a
>> >> >>>>>>>>>>>>>> comprehensible user interface and probably help the parallellization
>> >> >>>>>>>>>>>>>> effort.
>> >> >>>>>>>>>>>>>>
>> >> >>>>>>>>>>>>>> /Dag
>> >> >>>>>>>>>>>>> If the Function always owns the vector, there are cases you'll have
>> >> >>>>>>>>>>>>> to make unneccessary copies of a vector, in particular such scenarios
>> >> >>>>>>>>>>>>> may occur when trying to combine dolfin with something else.
>> >> >>>>>>>>>>>>>
>> >> >>>>>>>>>>>>> If the Function never owns the function space, it must always be
>> >> >>>>>>>>>>>>> constructed explicitly by the user. This may not be a bad thing.
>> >> >>>>>>>>>>>>> However, if the Function is loaded from a file, nobody owns the
>> >> >>>>>>>>>>>>> FunctionSpace.
>> >> >>>>>>>>>>>> Conceptually, I agree with Dag (and Martin?) that it is natural to
>> >> >>>>>>>>>>>> have global function spaces. And if the explicit construction of such
>> >> >>>>>>>>>>>> spaces can be made simple, it may not be a bad thing but a natural
>> >> >>>>>>>>>>>> part in setting up the mathematical problem. And I do not really like
>> >> >>>>>>>>>>>> that functions should be initialized from a form, which defines an
>> >> >>>>>>>>>>>> equation.
>> >> >>>>>>>>>>>>
>> >> >>>>>>>>>>>> I think one idea was to not force less mathematically oriented users
>> >> >>>>>>>>>>>> to worry about function spaces. I guess there are (at least) 2 types
>> >> >>>>>>>>>>>> of functions: (i) functions part of the form, and (ii) functions not
>> >> >>>>>>>>>>>> part of the form, but used in pre/postprocessing etc.
>> >> >>>>>>>>>>>>
>> >> >>>>>>>>>>>> For (i) it may be natural to construct the function space from the
>> >> >>>>>>>>>>>> form, and for (ii) it may be convenient in some cases, but it is not
>> >> >>>>>>>>>>>> really obvious that this is the best solution.
>> >> >>>>>>>>>>>>
>> >> >>>>>>>>>>>> Maybe an explicit construction of a function space can come with a
>> >> >>>>>>>>>>>> default, such as a nodal basis of piecewise linears?
>> >> >>>>>>>>>>>>
>> >> >>>>>>>>>>>> /Johan
>> >> >>>>>>>>>>> So:
>> >> >>>>>>>>>>>   FunctionSpace V(mesh);
>> >> >>>>>>>>>>>   Function f(V);
>> >> >>>>>>>>>>> gives a function f on piecewise linears?
>> >> >>>>>>>>>>> That's ok with me.
>> >> >>>>>>>>>>>
>> >> >>>>>>>>>>>
>> >> >>>>>>>>>>> About ownership, I think the only both robust and intuitive solution
>> >> >>>>>>>>>>> is that an object never should store a reference (or regular pointer)
>> >> >>>>>>>>>>> to another object. But as long as we are aware of the cost of doing
>> >> >>>>>>>>>>> this and state it clearly in the documentation, I'm ok with keeping
>> >> >>>>>>>>>>> e.g. the Mesh references like we do.
>> >> >>>>>>>>>>>
>> >> >>>>>>>>>>> Any time object A stores a reference to object B, the user must
>> >> >>>>>>>>>>> take care that the lifetime of B exceeds the lifetime of A. There
>> >> >>>>>>>>>>> are no exceptions to this. This puts some real limitations on the
>> >> >>>>>>>>>>> way the user must structure his program, e.g. he must sometimes
>> >> >>>>>>>>>>> (often?) keep objects around longer than they're explicitly needed.
>> >> >>>>>>>>>>>
>> >> >>>>>>>>>>> This may be a good thing, since it forces the user to think about
>> >> >>>>>>>>>>> dependencies and object lifetimes, and the objects in question
>> >> >>>>>>>>>>> use some memory.
>> >> >>>>>>>>>> I think this is ok. There are many ways to create a segfault in C++.
>> >> >>>>>>>>>> If you program in C++, you will have to think about memory.
>> >> >>>>>>>>>>
>> >> >>>>>>>>>>> But if we use references instead of shared_ptr,
>> >> >>>>>>>>>>> we should never have default values:
>> >> >>>>>>>>>>> - A Function has a reference to a Mesh, which is ok since
>> >> >>>>>>>>>>>   it's always created outside.
>> >> >>>>>>>>>>> - If a DiscreteFunction is to have a reference to a Vector, or a
>> >> >>>>>>>>>>>   Function is to have a reference to a FunctionSpace, it cannot
>> >> >>>>>>>>>>>   create its own without adding memory management code.
>> >> >>>>>>>>>>>
>> >> >>>>>>>>>>> Every place we accept these limitations and requirements of how
>> >> >>>>>>>>>>> the user structures his programs, we can use references and be
>> >> >>>>>>>>>>> done with it. But don't think that the pretty syntax means the user
>> >> >>>>>>>>>>> doesn't have to think about memory management, since all the
>> >> >>>>>>>>>>> responsibility for memory management (object destruction order)
>> >> >>>>>>>>>>> is in fact placed on the user, and errors from the users side will
>> >> >>>>>>>>>>> lead to invalid references we cannot detect and segfaults.
>> >> >>>>>>>>>> I want the syntax to be simple and pretty, but I don't necessarily
>> >> >>>>>>>>>> want to hide the user from problems that are part of the design of
>> >> >>>>>>>>>> C++. It isn't Python or Java. You should be expected to know what
>> >> >>>>>>>>>> you are doing. :-)
>> >> >>>>>>>>> It's not only about knowing what you're doing. It forces very
>> >> >>>>>>>>> hard restrictions on the design/flow of your program, which can
>> >> >>>>>>>>>
>> >> >>>>>>>>> 1) Be a major source of bugs in nontrivial apps (see Garths email), which
>> >> >>>>>>>>> are not locally visible because they depend on the global program flow.
>> >> >>>>>>>>>
>> >> >>>>>>>>> 2) Make it impossible to initialize e.g. Function in e.g a file reader,
>> >> >>>>>>>>> since the caller of the file reader would need to get the objects
>> >> >>>>>>>>> Function depends on. This is not limited to file readers, but is
>> >> >>>>>>>>> a recurring pattern in nontrivial apps.
>> >> >>>>>>>>>
>> >> >>>>>>>>> If we want to use dolfin or want dolfin to be used in apps that
>> >> >>>>>>>>> are more complicated than the traditional "read input, compute
>> >> >>>>>>>>> something, output something" app, these restrictions become
>> >> >>>>>>>>> a larger problem.
>> >> >>>>>>>>>
>> >> >>>>>>>>>> Anyway, I like the idea about having a FunctionSpace class which
>> >> >>>>>>>>>> several Functions may share. The problem we need to solve is
>> >> >>>>>>>>>> reading from file:
>> >> >>>>>>>>>>
>> >> >>>>>>>>>>  FunctionSpace V(mesh);
>> >> >>>>>>>>>>  Function u(V);
>> >> >>>>>>>>>>  file >> u;
>> >> >>>>>>>>>>
>> >> >>>>>>>>>> The last line just fills out the data in both u and V.
>> >> >>>>>>>>>>
>> >> >>>>>>>>>> This will lead to side effects as V might be changed when doing
>> >> >>>>>>>>>>
>> >> >>>>>>>>>>  FunctionSpace V(mesh);
>> >> >>>>>>>>>>  Function u(V);
>> >> >>>>>>>>>>  Function v(V);
>> >> >>>>>>>>>>  file >> u;
>> >> >>>>>>>>>>
>> >> >>>>>>>>>> V will be changed, both for u and v. In fact, the mesh will also be
>> >> >>>>>>>>>> changed.
>> >> >>>>>>>>>>
>> >> >>>>>>>>>> The best thing would be if we could do
>> >> >>>>>>>>>>
>> >> >>>>>>>>>>  file >> (mesh, V, u);
>> >> >>>>>>>>> This is _exactly_ the kind of issue that smart pointers solve.
>> >> >>>>>>>>>
>> >> >>>>>>>>> Btw, I tried to search the swig documentation for shared_ptr, and
>> >> >>>>>>>>> found nothing...
>> >> >>>>>>>>> I don't know what exactly they mean by "shared_ptr support".
>> >> >>>>>>>> It seems to be a set of typemaps that should kick in at the "right places".
>> >> >>>>>>>> They are defined in
>> >> >>>>>>>>
>> >> >>>>>>>>    <Lib/python/boost_shared_ptr.i>
>> >> >>>>>>>>
>> >> >>>>>>>> and used very rudimentary in
>> >> >>>>>>>>
>> >> >>>>>>>>   <Examples/test_suite/li_boost_shared_ptr.i>
>> >> >>>>>>>>
>> >> >>>>>>>> in the source tree of the 1.3.36 release. It seems that it is not specific for
>> >> >>>>>>>> boost::share_ptr but should also figure out tr1::shared_ptr
>> >> >>>>>>>>
>> >> >>>>>>>> I think:
>> >> >>>>>>>>
>> >> >>>>>>>>   %include <boost_shared_ptr.i>
>> >> >>>>>>>>
>> >> >>>>>>>> at the appropriate place should do the trick.
>> >> >>>>>>>>
>> >> >>>>>>>> Johan
>> >> >>>>>>> I was about to start sketching on a FunctionSpace class when I
>> >> >>>>>>> realized that this class might look different for different types of
>> >> >>>>>>> Functions.
>> >> >>>>>>>
>> >> >>>>>>> In particular, some functions (DiscreteFunctions) need to have a
>> >> >>>>>>> dof map and a finite element, while this is not needed for other
>> >> >>>>>>> functions.
>> >> >>>>>>>
>> >> >>>>>>> This means that we might need to duplicate the hierarchy of different
>> >> >>>>>>> function classes in a number of different function space classes.
>> >> >>>>>>>
>> >> >>>>>>> I see two other options:
>> >> >>>>>>>
>> >> >>>>>>> 1. Require that a function space consists of
>> >> >>>>>>>
>> >> >>>>>>>   mesh, finite_element, dof_map
>> >> >>>>>>>
>> >> >>>>>>> and pick a suitable (trivial) dof map for example constant functions.
>> >> >>>>>>>
>> >> >>>>>>> 2. Implement a number of different function space classes but have a
>> >> >>>>>>> single Function class which holds data but where each function call
>> >> >>>>>>> is passed on to the specific type of function call in the
>> >> >>>>>>> corresponding function space class.
>> >> >>>>>> Here are some more thoughts about the shared_ptr issue.
>> >> >>>>>>
>> >> >>>>>> As I understand it, using shared_ptr to store the data in
>> >> >>>>>> DiscreteFunction means that users will need to use shared_ptr.
>> >> >>>>>> In particular, a Function must be initialized with a shared_ptr to a
>> >> >>>>>> Mesh, not a reference to a Mesh. Then in all demos, we would have to
>> >> >>>>>> change from
>> >> >>>>>>
>> >> >>>>>>   Mesh mesh("mesh.xml");
>> >> >>>>>>
>> >> >>>>>> to
>> >> >>>>>>
>> >> >>>>>>   shared_ptr<Mesh> mesh = new Mesh("mesh.xml");
>> >> >>>>>>
>> >> >>>>>> This doesn't look very pretty, and it would make the interface
>> >> >>>>>> inconsistent. One would need to know that some classes likes Mesh
>> >> >>>>>> should "always" be stored in a shared_ptr (since it will most
>> >> >>>>>> certainly be used to initialize a Function), while other classes like
>> >> >>>>>> Matrix, File, KrylovSolver, can be created as usual.
>> >> >>>>>>
>> >> >>>>>> There seem to be two options:
>> >> >>>>>>
>> >> >>>>>> 1. Don't use shared_ptr anywhere (at least not in public interfaces)
>> >> >>>>>>
>> >> >>>>>> If we choose this option, we will need to decide who owns what. For
>> >> >>>>>> example, we could say that a Function always owns the Vector, but
>> >> >>>>>> never the Mesh.
>> >> >>>>>>
>> >> >>>>>> This works out also for reading from files:
>> >> >>>>>>
>> >> >>>>>>   Mesh mesh;
>> >> >>>>>>   Function u(mesh);
>> >> >>>>>>   file >> u;
>> >> >>>>>>
>> >> >>>>>> This would lead to side effects (the mesh will be changed), but that
>> >> >>>>>> might be ok.
>> >> >>>>>>
>> >> >>>>>> 2. Use shared_ptr everywhere
>> >> >>>>>>
>> >> >>>>>> If we choose this option, we need to redefine Mesh so that it is not a
>> >> >>>>>> mesh class, but a shared_ptr for a mesh class (which we rename to
>> >> >>>>>> something else). Sundance does something like this and it gives a nice
>> >> >>>>>> and pretty syntax:
>> >> >>>>>>
>> >> >>>>>>   Mesh mesh = new Mesh();             // no need to delete
>> >> >>>>>>   Function u = new Function(mesh);    // no need to delete
>> >> >>>>>>
>> >> >>>>>> Essentially, it would look like Java. See here:
>> >> >>>>>>
>> >> >>>>>>   http://software.sandia.gov/sundance/basics.html
>> >> >>>>>>
>> >> >>>>>> in the section "Handles and Memory Management".
>> >> >>>>>>
>> >> >>>>>> I think we need to pick one of these options and not a mix since it
>> >> >>>>>> would make the interface inconsistent.
>> >> >>>>>>
>> >> >>>>>> Opinions? 1 or 2?
>> >> >>>>>>
>> >> >>>>> I don't think that mixing is inconsistent. At the moment we use a
>> >> >>>>> mixture of references and pointers, and there is usually a logical
>> >> >>>>> choice. Most users never need to work with the pointers. If we use smart
>> >> >>>>> pointers, they will be used primarily internally in place of where we
>> >> >>>>> currently use plain pointers.
>> >> >>>>>
>> >> >>>>> For the specific case of a Function, part of the issue will be resolved
>> >> >>>>> by requiring that a Function own it's own vector (forgetting about
>> >> >>>>> sub-functions for the time being). Ownership of a finite element, or
>> >> >>>>> possibly a FunctionSpace in the future, remains an issue because it's
>> >> >>>>> natural that these are shared.
>> >> >>>>>
>> >> >>>>> If you look at uBLASVector which now uses a smart pointer, not a single
>> >> >>>>> piece of code outside uBLASVector had to be changed, but an advanced
>> >> >>>>> user can now create a uBLASVector from a uBLAS object and not worry
>> >> >>>>> about it going out of scope,
>> >> >>>>>
>> >> >>>>> Garth
>> >> >>>> I'm all for it if it can be hidden (or visible only by choice), but
>> >> >>>> isn't it the case that the mesh (or function space) will always have
>> >> >>>> to be a shared_ptr?
>> >> >>>>
>> >> >>>> In all applications, one will need to write
>> >> >>>>
>> >> >>>>   shared_ptr<Mesh> mesh = new Mesh(...);
>> >> >>>>
>> >> >>> I think that it's logical for a Function to have a reference to a Mesh.
>> >> >>> I'm not sure exactly the nicest way to deal with a FunctionSpace. If we
>> >> >>> want pretty syntax + plus smart pointer flexibility, a Function could
>> >> >>> have a pointer to the FunctionSpace, in which case it *does not* own the
>> >> >>> FunctionSpace, and a smart pointer in which case it may or may not  own
>> >> >>> the FunctionSpace.
>> >> >>>
>> >> >>> The class Function would look like
>> >> >>>
>> >> >>>   class Function
>> >> >>>   {
>> >> >>>     public:
>> >> >>>
>> >> >>>       Function(FunctionSpace& space) : V(space) {}
>> >> >>>
>> >> >>>       Function(shared_prt<FunctionSpace> space) : V(0), shared_V(space)
>> >> >>>         {}
>> >> >>>
>> >> >>>     private:
>> >> >>>       FunctionSpace* V;
>> >> >>>       shared_ptr<FunctionSpace>  shared_V;
>> >> >>>   }
>> >> >>>
>> >> >>> and a Function could be initialised by
>> >> >>>
>> >> >>>   FunctionSpace V;
>> >> >>>   Function u(V);
>> >> >>>
>> >> >>> in which case it does not own the FunctionSpace, or by
>> >> >>>
>> >> >>>   shared_ptr<FunctionSpace> V(new FunctionSpace);
>> >> >>>   Function u(V);
>> >> >>>
>> >> >>> in which case it may or may not own the FunctionSpace. The former would
>> >> >>> be used in simple programs, and the latter in, for example,
>> >> >>> LinearPDE::solve.
>> >> >>>
>> >> >>> Garth
>> >> >>>
>> >> >>>> ?
>> >> >> Sounds a bit dangerous to me... It does enable sharing objects
>> >> >> if you want to, but keeps the possibility to shoot of your foot.
>> >> >>
>> >> >> Anyway, a little tip: shared_ptr can take a deleter policy as input
>> >> >> to its constructor, so you don't need the extra FunctionSpace* to
>> >> >> implement this approach. Something like this should do the trick:
>> >> >>
>> >> >>       class NoDeleter {
>> >> >>         public:
>> >> >>           void operator() (FunctionSpace *p) {}
>> >> >>       };
>> >> >>
>> >>
>> >> I don't follow this.
>> >>
>> >> >> ...
>> >> >>
>> >> >>       Function(FunctionSpace& space) : shared_V(&space, NoDeleter()) {}
>> >> >>
>> >> >> (not tested, but modified from some working code).
>> >> >>
>> >> >>
>> >> >> And by the way, a thing to look out for is when using shared_ptr
>> >> >> for array pointers. Then you need to use a deleter policy like above,
>> >> >> except with function body { delete [] p; } instead.
>> >> >
>> >> > This sounds like the trick I was asking for earlier: a way to provide
>> >> > two different ways to initialize objects, either passing a reference
>> >> > or a shared_ptr (something like decrease_count).
>> >> >
>> >> > Do we want to do it this way?
>> >> >
>> >>
>> >> It's not ideal in terms of safety, but it does provide what we want  -
>> >> simplicity and flexibility.
>> >>
>> >> Garth
>> >
>> > ok, so yes? Martin?
>> >
>>
>> As Garth says. This allows the user to decide whether to shoot away
>> his foot or not. :-/
>
> Yes, this is what I want too. :-)
>
>> Also, this won't break the python interface since we
>> can continue to use references as before.
>
> Very nice, so let's do this.
>
> The next problem we need to solve is the design of the FunctionSpace
> class(es).
>
> Do we want to implement a hierarchy of FunctionSpaces? There seem to
> be two different options as sketched in a previous post.
>
> --
> Anders

I think... this is basically a placeholder for a set of variables,
so we can just have different constructors, allow null pointers
and use some enums.


Btw, this would be nice:
  f = project(g, V)

--
Martin