dolfin team mailing list archive

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

2008/9/8 Garth N. Wells <gnw20@xxxxxxxxx>:
>
>
> Anders Logg wrote:
>> On Sun, Sep 07, 2008 at 03:38:44PM +0100, Garth N. Wells wrote:
>>>
>>> Anders Logg wrote:
>>>> On Sun, Sep 07, 2008 at 03:29:14PM +0200, Anders Logg wrote:
>>>>> On Sun, Sep 07, 2008 at 03:11:51PM +0200, Martin Sandve Alnæs wrote:
>>>>>> 2008/9/7 Anders Logg <logg@xxxxxxxxx>:
>>>>>>> On Sun, Sep 07, 2008 at 08:27:41AM +0100, Garth N. Wells wrote:
>>>>>>>> Martin Sandve Alnæs wrote:
>>>>>>>>> 2008/9/6 Anders Logg <logg@xxxxxxxxx>:
>>>>>>>>>> On Sat, Sep 06, 2008 at 04:22:09PM +0200, Martin Sandve Alnæs wrote:
>>>>>>>>>>> 2008/9/6 Garth N. Wells <gnw20@xxxxxxxxx>:
>>>>>>>>>>>> Dag Lindbo wrote:
>>>>>>>>>>>>> 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?
>>>>>>>>>>>>>>
>>>>>>>>>>>>> One major issue which I just want to reiterate is ownership of data. As
>>>>>>>>>>>>> it stands, the DiscreteFunction may or may not be responsible for e.g.
>>>>>>>>>>>>> the dof vector x, depending on whether local_vector is a NULL pointer or
>>>>>>>>>>>>> not. Take a look at the thread "Ownership" from Garth on 06/26/2008.
>>>>>>>>>>>>>
>>>>>>>>>>>> Yes, this is a big problem and has caused me a few headaches with bugs.
>>>>>>>>>>>> For example, passing a user-defined Function to a function to convert it
>>>>>>>>>>>> to a DiscreteFunction via a projection onto a finite element basis
>>>>>>>>>>>> causes a problem because the FiniteElement which the projected Function
>>>>>>>>>>>> points to goes out of scope once the function is exited.
>>>>>>>>>>>>
>>>>>>>>>>>>> A problem related to this is initialization of the DiscreteFunction. We
>>>>>>>>>>>>> had a bug previously where the LinearPDE class maintained ownership of
>>>>>>>>>>>>> the solution vector. The only way to prevent this was to break the
>>>>>>>>>>>>> encapsulation of DiscreteFunction by making it a friend of LinearPDE (as
>>>>>>>>>>>>> with XMLFile for the same reasons). Here is some of the code that
>>>>>>>>>>>>> handles this initializaton today (L101 in LinearPDE.cpp):
>>>>>>>>>>>>>
>>>>>>>>>>>>>   u.init(mesh, *x, a, 1);
>>>>>>>>>>>>>   DiscreteFunction& uu = dynamic_cast<DiscreteFunction&>(*u.f);
>>>>>>>>>>>>>   uu.local_vector = x;
>>>>>>>>>>>>>
>>>>>>>>>>>>> This ain't poetry in my opinion :)
>>>>>>>>>>>>>
>>>>>>>>>>>> Indeed, this isn't nice, and there is something similar in XMLFile.cpp.
>>>>>>>>>>>>
>>>>>>>>>>>> Garth
>>>>>>>>>>> We should start to use std::tr1::shared_ptr. There is some support for it
>>>>>>>>>>> with python in swig 1.3.35, which is part of the upcoming Ubuntu Intrepid
>>>>>>>>>> The main issue is how we want to initialize Functions, and if one
>>>>>>>>>> should allow to set members.
>>>>>>>>>>
>>>>>>>>>> For simplicity, say that a Function is defined only by a Vector.
>>>>>>>>>> Then we have a few different situations to consider:
>>>>>>>>>>
>>>>>>>>>> 1. Function creates the Vector
>>>>>>>>>>
>>>>>>>>>>   Function u;
>>>>>>>>>>   Vector& x = u.vector();
>>>>>>>>>>
>>>>>>>>>> 2. Function gets the Vector
>>>>>>>>>>
>>>>>>>>>>   Vector x;
>>>>>>>>>>   Function u(x);
>>>>>>>>>>
>>>>>>>>>> 3. Function gets initialized with a Vector
>>>>>>>>>>
>>>>>>>>>>   Function u;
>>>>>>>>>>   Vector x;
>>>>>>>>>>   u.init(x);
>>>>>>>>>>
>>>>>>>>>> Do we want to support all of 1-3? Things become considerable easier if
>>>>>>>>>> we can make some simplifying assumptions.
>>>>>>>>>>
>>>>>>>>>> How visible would a shared_ptr be in the interface?
>>>>>>>>> A shared_ptr must be visible to the user every single place
>>>>>>>>> a pointer is passed around, otherwise the reference count
>>>>>>>>> won't be correct and we'll just have more problems.
>>>>>>>>>
>>>>>>>> So, in pseudo code, would it look something link this?
>>>>>>>>
>>>>>>>>    class DiscreteFunction
>>>>>>>>    {
>>>>>>>>      private:
>>>>>>>>
>>>>>>>>        shared_ptr<GenericVector> x;
>>>>>>>>
>>>>>>>>      public:
>>>>>>>>
>>>>>>>>        DiscreteFunction() : x(new Vector) {}
>>>>>>>>
>>>>>>>>        DiscreteFunction(shared_ptr<GenericVector> x)
>>>>>>>>        { x(x); }
>>>>>>>>
>>>>>>>>        shared_ptr<GenericVector> vec()
>>>>>>>>        {return x;}
>>>>>>>>    }
>>>>>>>> ?
>>>>>>>>
>>>>>>>> Garth
>>>>>>> What would the user code look like if we use shared_ptr for examples
>>>>>>> 1-3 above?
>>>>>>>
>>>>>>>>> 1. Function creates the Vector
>>>>>>>>>
>>>>>>>>>   Function u;
>>>>>>>>>   Vector& x = u.vector();
>>>>>> Function u;
>>>>>> Vector& x = u.vector(); // Storing this Vector& for later access is unsafe.
>>>>>> or
>>>>>> Function u;
>>>>>> shared_ptr<Vector> x = u.vector(); // Allows keeping the Vector around
>>>>>> after u is destroyed.
>>>>>>
>>>>>>
>>>>>>>>> 2. Function gets the Vector
>>>>>>>>>
>>>>>>>>>   Vector x;
>>>>>>>>>   Function u(x);
>>>>>> Vector x;
>>>>>> Function u(x); // Copy vector.
>>>>>>
>>>>>> shared_ptr<Vector> x = new Vector();
>>>>>> Function u(x); // Copy vector pointer, x or u may be deleted without
>>>>>> the other getting in trouble.
>>>>> I don't think the first option is what one might expect, and I don't
>>>>> think the second example looks very nice.
>>>>>
>>>>> We initialize Functions with a Mesh all the time and it would then be
>>>>> either very expensive to copy the mesh every time we create a Function
>>>>> from it (and one usually creates many functions on the same mesh), or
>>>>> we would have to write "shared_ptr" and "new" every time we used a
>>>>> Mesh.
>>>>>
>>>>> Isn't there another option? I don't like the all the flags we have now
>>>>> like is_view, local_vector, etc, but this looks worse.
>>>> Is there a way to increase the count for a shared_ptr?
>>>>
>>>> If there is, say a member named increase_ref(), we could do
>>>>
>>>>   class DiscreteFunction
>>>>   {
>>>>   public:
>>>>
>>>>       DiscreteFunction() : x(new Vector) {}
>>>>
>>>>       DiscreteFunction(GenericVector& x) : x(x);
>>>>       {
>>>>           x.increase_ref(1);
>>>>       }
>>>>
>>>>       DiscreteFunction(shared_ptr<GenericVector> x) : x(x) {}
>>>>       { x(x); }
>>>>
>>>>       shared_ptr<GenericVector> vec()
>>>>       {return x;}
>>>>
>>>>   private:
>>>>
>>>>       shared_ptr<GenericVector> x;
>>>>
>>>>   };
>>>>
>>>> Then it would be possible to do
>>>>
>>>>   Vector x;
>>>>   Function u(x);
>>>>
>>>> as we can now and the Function u would know that someone else is
>>>> responsible for deleting the data.
>>>>
>>>> Then if one writes code where the Vector goes out of scope, one must
>>>> use a shared_ptr, but not otherwise. We would not force everyone to
>>>> use shared_ptr all the time.
>>>>
>>> I think that having two options would lead to confusion.
>>>
>>> On a another note, I don't really like the above code anyway. We could
>>> avoid this in many cases by adding the constructor
>>>
>>>    Function::Function(Mesh& mesh, Form& form, uint i = 1);
>>>
>>> and letting Function create the Vector. Function already has a member
>>> function for accessing the underlying Vector.
>>>
>>> Garth
>>
>> That's what I suggested earlier. It would simplify many things if we
>> always know who owns what. We could assume that the Function never
>> owns the Mesh but always owns the Vector (for discrete functions).
>>
>> If someone has a Vector and wants to create a Function from it, then
>> just do
>>
>>   Vector x;
>>   Function u;
>>   u.vector() = x; // Copy data
>>
>
> OK, and if someone really wants a Function to point to x, they can do
>
>     shared_ptr<Vector> x(new Vector);
>     Function u;
>     u.vector() = x; // Share data
>
> If Function uses shared_ptr internally for its Vector, it doesn't have
> to worry about keeping track of ownership.
>
> Garth

You can't overload a function on its return type.

For this to work:
>>   u.vector() = x; // Copy data
vector() must return a reference to the Vector.

For this to work:
>     u.vector() = x; // Share data
vector() must return a reference to the shared_ptr.

What's wrong with just using a set-function?

shared_ptr<Vector> v = u.vector();

shared_ptr<Vector> shared_ptr_x;
u.setVector(shared_ptr_x); // Share data

Vector vector_x;
u.setVector(vector_x); // Copy data


Also, if u.vector returns a shared_ptr, you can't assign it to a
reference, but you should be able to dereference it:

Vector & v = *u.vector();

--
Martin