dolfin team mailing list archive

[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

Attachment: signature.asc
Description: Digital signature