Anders Logg wrote:
> On Tue, Sep 09, 2008 at 04:31:41PM +0200, Martin Sandve Alnæs wrote:
>> 2008/9/9 Garth N. Wells <[EMAIL PROTECTED]>:
>>>
>>> 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 <[EMAIL PROTECTED]>:
>>>>>>>>>> On Mon, Sep 08, 2008 at 11:12:14AM +0200, Martin Sandve Alnæs wrote:
>>>>>>>>>>> 2008/9/8 Johan Hoffman <[EMAIL PROTECTED]>:
>>>>>>>>>>>>> 2008/9/8 Dag Lindbo <[EMAIL PROTECTED]>:
>>>>>>>>>>>>>> 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
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> DOLFIN-dev mailing list
> [email protected]
> http://www.fenics.org/mailman/listinfo/dolfin-dev
_______________________________________________
DOLFIN-dev mailing list
[email protected]
http://www.fenics.org/mailman/listinfo/dolfin-dev
