dolfin team mailing list archive

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

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?
>
> --
> Anders

As Garth says. This allows the user to decide whether to shoot away
his foot or not. :-/

Also, this won't break the python interface since we
can continue to use references as before.

--
Martin