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
> ?
>
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> 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
