Martin Sandve Alnæs wrote:
> 2008/9/9 Garth N. Wells <[EMAIL PROTECTED]>:
>>
>> 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()) {}
>
> 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.
>
Oh, so it's meant to do nothing. Now I follow.
Garth
> --
> Martin
_______________________________________________
DOLFIN-dev mailing list
[email protected]
http://www.fenics.org/mailman/listinfo/dolfin-dev
