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 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) {}
      };

...

      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.

--
Martin