dolfin team mailing list archive

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

2008/9/9 Garth N. Wells <gnw20@xxxxxxxxx>:
>
>
> 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()) {}

When shared_V is deleted, its destructor will decrement its reference
counter. If it reaches 0, it would usually call "delete p;" where p is
the pointer it holds. But in this case, instead of calling "delete
p;", it will call upon the given deleter objects () operator like
"user_defined_deleter(p);", where user_defined_deleter is the object
of type NoDeleter created above.

This user-defined deleter object can also hold pointers to external
reference counters, and decrement those instead. Or simply do nothing
like here.

--
Martin