dolfin team mailing list archive

[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

Attachment: signature.asc
Description: Digital signature