BoundaryMesh.coordinates () gives random results

[Marie Rognes]

I'm tracking down some weirdness when extracting coordinates from a BoundaryMesh.
This seems to be one of the pecularities:

Using

  bdry = BoundaryMesh(mesh)
  coordinates = bdry.coordinates()

seems to work ok, combining the two don't.

---

from dolfin import *

mesh = UnitCube(1, 1, 1)

# This works
#bdry = BoundaryMesh(mesh)
#coordinates = bdry.coordinates()

# This gives funny and random results.
coordinates = BoundaryMesh(mesh).coordinates()
for p in coordinates:
    print p

--

gives on first run:

[ 6.92038619e-310 6.92038619e-310 0.00000000e+000]
[ 1. 1. 0.]
[ 1. 1. 1.]
[ 0. 0. 0.]
[ 1. 0. 1.]
[ 0. 0. 1.]
[ 0. 1. 0.]
[ 0. 1. 1.]

and on second run:

[ 6.93937331e-310 6.93937422e-310 nan]
[ 1.00000000e+000 0.00000000e+000 4.27920928e-317]
[ 6.42285340e-323 2.49208846e-316 1.91505457e-316]
[ 1.48800162e-316 1.92286160e-316 0.00000000e+000]
[ 1.52283403e-316 0.00000000e+000 0.00000000e+000]
[ 0. 0. 0.]
[ 1.48155347e-316 9.44288127e-312 2.58883486e-312]
[ 6.02760088e-322 1.00000000e+000 1.00000000e+000]

[Anders Logg (logg)]

Looks like a SWIG layer Array issue.

--
Anders

On Tue, Mar 27, 2012 at 08:19:43PM -0000, Marie Rognes wrote:
> Public bug reported:
>
>
> I'm tracking down some weirdness when extracting coordinates from a BoundaryMesh.
> This seems to be one of the pecularities:
>
> Using
>
> bdry = BoundaryMesh(mesh)
> coordinates = bdry.coordinates()
>
> seems to work ok, combining the two don't.
>
> ---
>
> from dolfin import *
>
> mesh = UnitCube(1, 1, 1)
>
> # This works
> #bdry = BoundaryMesh(mesh)
> #coordinates = bdry.coordinates()
>
> # This gives funny and random results.
> coordinates = BoundaryMesh(mesh).coordinates()
> for p in coordinates:
> print p
>

[Johan Hake (johan-hake)]

This is because the mesh has gone out of scope and destroyed the coordinatedata which was wrapped to a NumPy array. The same happens if you do:

  coord = UnitCube(1, 1, 1).coordinates()

We can solve this by making adding a deepcopy argument to the coordinates (and the cells) methods, which defaults to True. This might break code as some (hmpfr me!) uses coordinates to move a mesh. Changing the coordinates when a deepcopy is done wont change the coordinates on the mesh.

Another way, is to add a reference of the mesh object to the NumPy array. Then the Mesh object wont go out of scope until the coordinate array does. The latter is more tricky to implement, but would be the best way to do it.

[Garth Wells (garth-wells)]

On 28 March 2012 14:03, Johan Hake <email address hidden> wrote:
> This is because the mesh has gone out of scope and destroyed the
> coordinatedata which was wrapped to a NumPy array. The same happens if
> you do:
>
>  coord = UnitCube(1, 1, 1).coordinates()
>
> We can solve this by making adding a deepcopy argument to the
> coordinates (and the cells) methods, which defaults to True. This might
> break code as some (hmpfr me!) uses coordinates to move a mesh. Changing
> the coordinates when a deepcopy is done wont change the  coordinates on
> the mesh.
>
> Another way, is to add a reference of the mesh object to the NumPy
> array. Then the Mesh object wont go out of scope until the coordinate
> array does. The latter is more tricky to implement, but would be the
> best way to do it.
>

We've seen before that keeping references to large objects will at
some point cause a memory leak. I suggest that we

1. Have one function that returns a copy of the coordinates and one
that provides view

or

2. Make the coordinate data a boost::shared_array<double> or a
boost::shared_ptr<std::vector<double> > object.

Garth

> ** Changed in: dolfin
>       Status: Confirmed => In Progress
>
> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/966553
>
> Title:
>  BoundaryMesh.coordinates () gives random results
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/966553/+subscriptions

[Johan Hake (johan-hake)]

On 03/28/2012 09:41 AM, Garth Wells wrote:
> On 28 March 2012 14:03, Johan Hake<email address hidden> wrote:
>> This is because the mesh has gone out of scope and destroyed the
>> coordinatedata which was wrapped to a NumPy array. The same happens if
>> you do:
>>
>> coord = UnitCube(1, 1, 1).coordinates()
>>
>> We can solve this by making adding a deepcopy argument to the
>> coordinates (and the cells) methods, which defaults to True. This might
>> break code as some (hmpfr me!) uses coordinates to move a mesh. Changing
>> the coordinates when a deepcopy is done wont change the coordinates on
>> the mesh.
>>
>> Another way, is to add a reference of the mesh object to the NumPy
>> array. Then the Mesh object wont go out of scope until the coordinate
>> array does. The latter is more tricky to implement, but would be the
>> best way to do it.
>>
>
> We've seen before that keeping references to large objects will at
> some point cause a memory leak.

The problem with the previous reference storage, was that it was a
global storage. My suggestions is a storage connected to an object. If
an object returns data it owes it should make sure it does not destroy
it for the lifetime of that data. We already do this for for example
down_casting of objects.

> I suggest that we
>
> 1. Have one function that returns a copy of the coordinates and one
> that provides view

I tend to agree, but the question is what should be default. The most
conservative would be to let a copy be default, but that would most
probably break user code, as most people (hmpfr me) uses coordinates to
move meshes.

> or
>
> 2. Make the coordinate data a boost::shared_array<double> or a
> boost::shared_ptr<std::vector<double> > object.

This would introduce a somewhat middle solution, which is pretty
cumbersome to implement and which need to be changed anyhow, when we
revamp the mesh interface.

I vote for adding a reference to the returned object, which happens to
be the solution I just implemented and are now testing ;)

Johan

>
> Garth
>
>> ** Changed in: dolfin
>> Status: Confirmed => In Progress
>>
>> --
>> You received this bug notification because you are a member of DOLFIN
>> Core Team, which is subscribed to DOLFIN.
>> https://bugs.launchpad.net/bugs/966553
>>
>> Title:
>> BoundaryMesh.coordinates () gives random results
>>
>> To manage notifications about this bug go to:
>> https://bugs.launchpad.net/dolfin/+bug/966553/+subscriptions
>

[Anders Logg (logg)]

On Wed, Mar 28, 2012 at 08:06:02AM -0000, Johan Hake wrote:
> On 03/28/2012 09:41 AM, Garth Wells wrote:
> > On 28 March 2012 14:03, Johan Hake<email address hidden> wrote:
> >> This is because the mesh has gone out of scope and destroyed the
> >> coordinatedata which was wrapped to a NumPy array. The same happens if
> >> you do:
> >>
> >> coord = UnitCube(1, 1, 1).coordinates()
> >>
> >> We can solve this by making adding a deepcopy argument to the
> >> coordinates (and the cells) methods, which defaults to True. This might
> >> break code as some (hmpfr me!) uses coordinates to move a mesh. Changing
> >> the coordinates when a deepcopy is done wont change the coordinates on
> >> the mesh.
> >>
> >> Another way, is to add a reference of the mesh object to the NumPy
> >> array. Then the Mesh object wont go out of scope until the coordinate
> >> array does. The latter is more tricky to implement, but would be the
> >> best way to do it.
> >>
> >
> > We've seen before that keeping references to large objects will at
> > some point cause a memory leak.
>
> The problem with the previous reference storage, was that it was a
> global storage. My suggestions is a storage connected to an object. If
> an object returns data it owes it should make sure it does not destroy
> it for the lifetime of that data. We already do this for for example
> down_casting of objects.
>
> > I suggest that we
> >
> > 1. Have one function that returns a copy of the coordinates and one
> > that provides view
>
> I tend to agree, but the question is what should be default. The most
> conservative would be to let a copy be default, but that would most
> probably break user code, as most people (hmpfr me) uses coordinates to
> move meshes.

I think we should try hard to avoid that. Accessing coordinates() by
reference (not copy) has been the default behavior and we would break
user code, possibly also some code examples from the book.

> > or
> >
> > 2. Make the coordinate data a boost::shared_array<double> or a
> > boost::shared_ptr<std::vector<double> > object.
>
> This would introduce a somewhat middle solution, which is pretty
> cumbersome to implement and which need to be changed anyhow, when we
> revamp the mesh interface.

Is it up for revamping???

--
Anders

> I vote for adding a reference to the returned object, which happens to
> be the solution I just implemented and are now testing ;)

> Johan
>
> >
> > Garth
> >
> >> ** Changed in: dolfin
> >> Status: Confirmed => In Progress
> >>
> >>
> >> Title:
> >> BoundaryMesh.coordinates () gives random results
> >>
> >> To manage notifications about this bug go to:
> >> https://bugs.launchpad.net/dolfin/+bug/966553/+subscriptions
> >
>

[Johan Hake (johan-hake)]

On 03/28/2012 10:26 AM, Anders Logg wrote:
> On Wed, Mar 28, 2012 at 08:06:02AM -0000, Johan Hake wrote:
>> On 03/28/2012 09:41 AM, Garth Wells wrote:
>>> On 28 March 2012 14:03, Johan Hake<email address hidden> wrote:
>>>> This is because the mesh has gone out of scope and destroyed the
>>>> coordinatedata which was wrapped to a NumPy array. The same happens if
>>>> you do:
>>>>
>>>> coord = UnitCube(1, 1, 1).coordinates()
>>>>
>>>> We can solve this by making adding a deepcopy argument to the
>>>> coordinates (and the cells) methods, which defaults to True. This might
>>>> break code as some (hmpfr me!) uses coordinates to move a mesh. Changing
>>>> the coordinates when a deepcopy is done wont change the coordinates on
>>>> the mesh.
>>>>
>>>> Another way, is to add a reference of the mesh object to the NumPy
>>>> array. Then the Mesh object wont go out of scope until the coordinate
>>>> array does. The latter is more tricky to implement, but would be the
>>>> best way to do it.
>>>>
>>>
>>> We've seen before that keeping references to large objects will at
>>> some point cause a memory leak.
>>
>> The problem with the previous reference storage, was that it was a
>> global storage. My suggestions is a storage connected to an object. If
>> an object returns data it owes it should make sure it does not destroy
>> it for the lifetime of that data. We already do this for for example
>> down_casting of objects.
>>
>>> I suggest that we
>>>
>>> 1. Have one function that returns a copy of the coordinates and one
>>> that provides view
>>
>> I tend to agree, but the question is what should be default. The most
>> conservative would be to let a copy be default, but that would most
>> probably break user code, as most people (hmpfr me) uses coordinates to
>> move meshes.
>
> I think we should try hard to avoid that. Accessing coordinates() by
> reference (not copy) has been the default behavior and we would break
> user code, possibly also some code examples from the book.
>
>>> or
>>>
>>> 2. Make the coordinate data a boost::shared_array<double> or a
>>> boost::shared_ptr<std::vector<double> > object.
>>
>> This would introduce a somewhat middle solution, which is pretty
>> cumbersome to implement and which need to be changed anyhow, when we
>> revamp the mesh interface.
>
> Is it up for revamping???

:)

You keep on talking about this. How on earth are we going to support
coarsening, refinement, and parallel such if we do not make the storage
of coordinates, and connectivities more flexible. And if that will be
the case, providing a numpy array view to these entities need to be
rethinked.

Johan

> --
> Anders
>
>
>> I vote for adding a reference to the returned object, which happens to
>> be the solution I just implemented and are now testing ;)
>
>
>> Johan
>>
>>>
>>> Garth
>>>
>>>> ** Changed in: dolfin
>>>> Status: Confirmed => In Progress
>>>>
>>>>
>>>> Title:
>>>> BoundaryMesh.coordinates () gives random results
>>>>
>>>> To manage notifications about this bug go to:
>>>> https://bugs.launchpad.net/dolfin/+bug/966553/+subscriptions
>>>
>>
>

[Anders Logg (logg)]

On Wed, Mar 28, 2012 at 08:42:15AM -0000, Johan Hake wrote:
> On 03/28/2012 10:26 AM, Anders Logg wrote:
> > On Wed, Mar 28, 2012 at 08:06:02AM -0000, Johan Hake wrote:
> >> On 03/28/2012 09:41 AM, Garth Wells wrote:
> >>> On 28 March 2012 14:03, Johan Hake<email address hidden> wrote:
> >>>> This is because the mesh has gone out of scope and destroyed the
> >>>> coordinatedata which was wrapped to a NumPy array. The same happens if
> >>>> you do:
> >>>>
> >>>> coord = UnitCube(1, 1, 1).coordinates()
> >>>>
> >>>> We can solve this by making adding a deepcopy argument to the
> >>>> coordinates (and the cells) methods, which defaults to True. This might
> >>>> break code as some (hmpfr me!) uses coordinates to move a mesh. Changing
> >>>> the coordinates when a deepcopy is done wont change the coordinates on
> >>>> the mesh.
> >>>>
> >>>> Another way, is to add a reference of the mesh object to the NumPy
> >>>> array. Then the Mesh object wont go out of scope until the coordinate
> >>>> array does. The latter is more tricky to implement, but would be the
> >>>> best way to do it.
> >>>>
> >>>
> >>> We've seen before that keeping references to large objects will at
> >>> some point cause a memory leak.
> >>
> >> The problem with the previous reference storage, was that it was a
> >> global storage. My suggestions is a storage connected to an object. If
> >> an object returns data it owes it should make sure it does not destroy
> >> it for the lifetime of that data. We already do this for for example
> >> down_casting of objects.
> >>
> >>> I suggest that we
> >>>
> >>> 1. Have one function that returns a copy of the coordinates and one
> >>> that provides view
> >>
> >> I tend to agree, but the question is what should be default. The most
> >> conservative would be to let a copy be default, but that would most
> >> probably break user code, as most people (hmpfr me) uses coordinates to
> >> move meshes.
> >
> > I think we should try hard to avoid that. Accessing coordinates() by
> > reference (not copy) has been the default behavior and we would break
> > user code, possibly also some code examples from the book.
> >
> >>> or
> >>>
> >>> 2. Make the coordinate data a boost::shared_array<double> or a
> >>> boost::shared_ptr<std::vector<double> > object.
> >>
> >> This would introduce a somewhat middle solution, which is pretty
> >> cumbersome to implement and which need to be changed anyhow, when we
> >> revamp the mesh interface.
> >
> > Is it up for revamping???
>
> :)
>
> You keep on talking about this. How on earth are we going to support
> coarsening, refinement, and parallel such if we do not make the storage
> of coordinates, and connectivities more flexible. And if that will be
> the case, providing a numpy array view to these entities need to be
> rethinked.

My view is that we should be able to keep the current "static" data
structure and copy to dynamic data structures during refinement etc.

At any rate, I think we should definitely not change the interface
(only the implementation).

--
Anders

> Johan
>
> >
> >
> >> I vote for adding a reference to the returned object, which happens to
> >> be the solution I just...

[Garth Wells (garth-wells)]

On 28 March 2012 16:42, Johan Hake <email address hidden> wrote:
> On 03/28/2012 10:26 AM, Anders Logg wrote:
>> On Wed, Mar 28, 2012 at 08:06:02AM -0000, Johan Hake wrote:
>>> On 03/28/2012 09:41 AM, Garth Wells wrote:
>>>> On 28 March 2012 14:03, Johan Hake<email address hidden>   wrote:
>>>>> This is because the mesh has gone out of scope and destroyed the
>>>>> coordinatedata which was wrapped to a NumPy array. The same happens if
>>>>> you do:
>>>>>
>>>>>    coord = UnitCube(1, 1, 1).coordinates()
>>>>>
>>>>> We can solve this by making adding a deepcopy argument to the
>>>>> coordinates (and the cells) methods, which defaults to True. This might
>>>>> break code as some (hmpfr me!) uses coordinates to move a mesh. Changing
>>>>> the coordinates when a deepcopy is done wont change the  coordinates on
>>>>> the mesh.
>>>>>
>>>>> Another way, is to add a reference of the mesh object to the NumPy
>>>>> array. Then the Mesh object wont go out of scope until the coordinate
>>>>> array does. The latter is more tricky to implement, but would be the
>>>>> best way to do it.
>>>>>
>>>>
>>>> We've seen before that keeping references to large objects will at
>>>> some point cause a memory leak.
>>>
>>> The problem with the previous reference storage, was that it was a
>>> global storage. My suggestions is a storage connected to an object. If
>>> an object returns data it owes it should make sure it does not destroy
>>> it for the lifetime of that data. We already do this for for example
>>> down_casting of objects.

The underlying problem is the same, and I think a poor design. I would
not expect that a simple array object would keep a possibly very large
object in scope.

>>>
>>>> I suggest that we
>>>>
>>>> 1. Have  one function that returns a copy of the coordinates and one
>>>> that provides view
>>>
>>> I tend to agree, but the question is what should be default. The most
>>> conservative would be to let a copy be default, but that would most
>>> probably break user code, as most people (hmpfr me) uses coordinates to
>>> move meshes.
>>
>> I think we should try hard to avoid that. Accessing coordinates() by
>> reference (not copy) has been the default behavior and we would break
>> user code, possibly also some code examples from the book.
>>
>>>> or
>>>>
>>>> 2. Make the coordinate data a boost::shared_array<double>   or a
>>>> boost::shared_ptr<std::vector<double>   >   object.
>>>
>>> This would introduce a somewhat middle solution, which is pretty
>>> cumbersome to implement and which need to be changed anyhow, when we
>>> revamp the mesh interface.
>>

It should be easy to implement - I have already changed some plain
arrays to std::vector. Relying on the contiguous storage that a
std::vector provides, it's easy to interface with functions looking to
accept a plain pointer. The question is how it should work in the SWIG
layer.

>> Is it up for revamping???
>
> :)
>
> You keep on talking about this. How on earth are we going to support
> coarsening, refinement, and parallel such if we do not make the storage
> of coordinates, and connectivities more flexible. And if that will be
> the case, providing a numpy array view to these e...

[Johan Hake (johan-hake)]

On 03/28/2012 12:31 PM, Garth Wells wrote:
> On 28 March 2012 16:42, Johan Hake<email address hidden> wrote:
>> On 03/28/2012 10:26 AM, Anders Logg wrote:
>>> On Wed, Mar 28, 2012 at 08:06:02AM -0000, Johan Hake wrote:
>>>> On 03/28/2012 09:41 AM, Garth Wells wrote:
>>>>> On 28 March 2012 14:03, Johan Hake<email address hidden> wrote:
>>>>>> This is because the mesh has gone out of scope and destroyed the
>>>>>> coordinatedata which was wrapped to a NumPy array. The same happens if
>>>>>> you do:
>>>>>>
>>>>>> coord = UnitCube(1, 1, 1).coordinates()
>>>>>>
>>>>>> We can solve this by making adding a deepcopy argument to the
>>>>>> coordinates (and the cells) methods, which defaults to True. This might
>>>>>> break code as some (hmpfr me!) uses coordinates to move a mesh. Changing
>>>>>> the coordinates when a deepcopy is done wont change the coordinates on
>>>>>> the mesh.
>>>>>>
>>>>>> Another way, is to add a reference of the mesh object to the NumPy
>>>>>> array. Then the Mesh object wont go out of scope until the coordinate
>>>>>> array does. The latter is more tricky to implement, but would be the
>>>>>> best way to do it.
>>>>>>
>>>>>
>>>>> We've seen before that keeping references to large objects will at
>>>>> some point cause a memory leak.
>>>>
>>>> The problem with the previous reference storage, was that it was a
>>>> global storage. My suggestions is a storage connected to an object. If
>>>> an object returns data it owes it should make sure it does not destroy
>>>> it for the lifetime of that data. We already do this for for example
>>>> down_casting of objects.
>
> The underlying problem is the same, and I think a poor design. I would
> not expect that a simple array object would keep a possibly very large
> object in scope.

The whole slicing interface of NumPy is built on this design.

   a = np.array([some big list])
   b = a[::2]

b is here a view into a. Changing values in b will change values in a.
If a is destroyed b is still working, as it knows it came from a and
therefor keeps a reference to it. NumPy provides the interface, we need
the functionality, why not use it?

The coordinates of a mesh is not a small array. Nor are other examples
which are natural to talk about in this regard.

   MeshFunction.array, GenericFoo.data, Mesh.cells, Mesh.coordinates

>>>>> I suggest that we
>>>>>
>>>>> 1. Have one function that returns a copy of the coordinates and one
>>>>> that provides view
>>>>
>>>> I tend to agree, but the question is what should be default. The most
>>>> conservative would be to let a copy be default, but that would most
>>>> probably break user code, as most people (hmpfr me) uses coordinates to
>>>> move meshes.
>>>
>>> I think we should try hard to avoid that. Accessing coordinates() by
>>> reference (not copy) has been the default behavior and we would break
>>> user code, possibly also some code examples from the book.
>>>
>>>>> or
>>>>>
>>>>> 2. Make the coordinate data a boost::shared_array<double> or a
>>>>> boost::shared_ptr<std::vector<double> > object.
>>>>
>>>> This would introduce a somewhat middle solution, which is pretty
>>>> cumbersome to implement and wh...

[Johan Hake (johan-hake)]

On 03/28/2012 12:20 PM, Anders Logg wrote:
> On Wed, Mar 28, 2012 at 08:42:15AM -0000, Johan Hake wrote:
>> On 03/28/2012 10:26 AM, Anders Logg wrote:
>>> On Wed, Mar 28, 2012 at 08:06:02AM -0000, Johan Hake wrote:
>>>> On 03/28/2012 09:41 AM, Garth Wells wrote:
>>>>> On 28 March 2012 14:03, Johan Hake<email address hidden> wrote:
>>>>>> This is because the mesh has gone out of scope and destroyed the
>>>>>> coordinatedata which was wrapped to a NumPy array. The same happens if
>>>>>> you do:
>>>>>>
>>>>>> coord = UnitCube(1, 1, 1).coordinates()
>>>>>>
>>>>>> We can solve this by making adding a deepcopy argument to the
>>>>>> coordinates (and the cells) methods, which defaults to True. This might
>>>>>> break code as some (hmpfr me!) uses coordinates to move a mesh. Changing
>>>>>> the coordinates when a deepcopy is done wont change the coordinates on
>>>>>> the mesh.
>>>>>>
>>>>>> Another way, is to add a reference of the mesh object to the NumPy
>>>>>> array. Then the Mesh object wont go out of scope until the coordinate
>>>>>> array does. The latter is more tricky to implement, but would be the
>>>>>> best way to do it.
>>>>>>
>>>>>
>>>>> We've seen before that keeping references to large objects will at
>>>>> some point cause a memory leak.
>>>>
>>>> The problem with the previous reference storage, was that it was a
>>>> global storage. My suggestions is a storage connected to an object. If
>>>> an object returns data it owes it should make sure it does not destroy
>>>> it for the lifetime of that data. We already do this for for example
>>>> down_casting of objects.
>>>>
>>>>> I suggest that we
>>>>>
>>>>> 1. Have one function that returns a copy of the coordinates and one
>>>>> that provides view
>>>>
>>>> I tend to agree, but the question is what should be default. The most
>>>> conservative would be to let a copy be default, but that would most
>>>> probably break user code, as most people (hmpfr me) uses coordinates to
>>>> move meshes.
>>>
>>> I think we should try hard to avoid that. Accessing coordinates() by
>>> reference (not copy) has been the default behavior and we would break
>>> user code, possibly also some code examples from the book.
>>>
>>>>> or
>>>>>
>>>>> 2. Make the coordinate data a boost::shared_array<double> or a
>>>>> boost::shared_ptr<std::vector<double> > object.
>>>>
>>>> This would introduce a somewhat middle solution, which is pretty
>>>> cumbersome to implement and which need to be changed anyhow, when we
>>>> revamp the mesh interface.
>>>
>>> Is it up for revamping???
>>
>> :)
>>
>> You keep on talking about this. How on earth are we going to support
>> coarsening, refinement, and parallel such if we do not make the storage
>> of coordinates, and connectivities more flexible. And if that will be
>> the case, providing a numpy array view to these entities need to be
>> rethinked.
>
> My view is that we should be able to keep the current "static" data
> structure and copy to dynamic data structures during refinement etc.

Ok, but that is a somewhat premature discussion I guess.

> At any rate, I think we should definitely not change the interface
> (only the implementation).

Ok.

Jo...

[Anders Logg (logg)]

On Wed, Mar 28, 2012 at 11:31:35AM -0000, Johan Hake wrote:
> On 03/28/2012 12:31 PM, Garth Wells wrote:
> > On 28 March 2012 16:42, Johan Hake<email address hidden> wrote:
> >> On 03/28/2012 10:26 AM, Anders Logg wrote:
> >>> On Wed, Mar 28, 2012 at 08:06:02AM -0000, Johan Hake wrote:
> >>>> On 03/28/2012 09:41 AM, Garth Wells wrote:
> >>>>> On 28 March 2012 14:03, Johan Hake<email address hidden> wrote:
> >>>>>> This is because the mesh has gone out of scope and destroyed the
> >>>>>> coordinatedata which was wrapped to a NumPy array. The same happens if
> >>>>>> you do:
> >>>>>>
> >>>>>> coord = UnitCube(1, 1, 1).coordinates()
> >>>>>>
> >>>>>> We can solve this by making adding a deepcopy argument to the
> >>>>>> coordinates (and the cells) methods, which defaults to True. This might
> >>>>>> break code as some (hmpfr me!) uses coordinates to move a mesh. Changing
> >>>>>> the coordinates when a deepcopy is done wont change the coordinates on
> >>>>>> the mesh.
> >>>>>>
> >>>>>> Another way, is to add a reference of the mesh object to the NumPy
> >>>>>> array. Then the Mesh object wont go out of scope until the coordinate
> >>>>>> array does. The latter is more tricky to implement, but would be the
> >>>>>> best way to do it.
> >>>>>>
> >>>>>
> >>>>> We've seen before that keeping references to large objects will at
> >>>>> some point cause a memory leak.
> >>>>
> >>>> The problem with the previous reference storage, was that it was a
> >>>> global storage. My suggestions is a storage connected to an object. If
> >>>> an object returns data it owes it should make sure it does not destroy
> >>>> it for the lifetime of that data. We already do this for for example
> >>>> down_casting of objects.
> >
> > The underlying problem is the same, and I think a poor design. I would
> > not expect that a simple array object would keep a possibly very large
> > object in scope.
>
> The whole slicing interface of NumPy is built on this design.
>
> a = np.array([some big list])
> b = a[::2]
>
> b is here a view into a. Changing values in b will change values in a.
> If a is destroyed b is still working, as it knows it came from a and
> therefor keeps a reference to it. NumPy provides the interface, we need
> the functionality, why not use it?
>
> The coordinates of a mesh is not a small array. Nor are other examples
> which are natural to talk about in this regard.
>
> MeshFunction.array, GenericFoo.data, Mesh.cells, Mesh.coordinates

I think Garth meant that the Mesh is a big data structure which would
then be attached to the array.

I don't see that as a big problem. The coordinates array is itself of
considerable size compared to the Mesh, and if someone calls
Mesh.coordinates() they will likely assume that the mesh is still
alive. If not, they could call Mesh.coordinates().copy().

--
Anders

> >>>>> I suggest that we
> >>>>>
> >>>>> 1. Have one function that returns a copy of the coordinates and one
> >>>>> that provides view
> >>>>
> >>>> I tend to agree, but the question is what should be default. The most
> >>>> conservative would be to let a copy be default, but that would most
> >>>> probably break user code, as ...

[Garth Wells (garth-wells)]

On 28 March 2012 19:31, Johan Hake <email address hidden> wrote:
> On 03/28/2012 12:31 PM, Garth Wells wrote:
>> On 28 March 2012 16:42, Johan Hake<email address hidden>  wrote:
>>> On 03/28/2012 10:26 AM, Anders Logg wrote:
>>>> On Wed, Mar 28, 2012 at 08:06:02AM -0000, Johan Hake wrote:
>>>>> On 03/28/2012 09:41 AM, Garth Wells wrote:
>>>>>> On 28 March 2012 14:03, Johan Hake<email address hidden>     wrote:
>>>>>>> This is because the mesh has gone out of scope and destroyed the
>>>>>>> coordinatedata which was wrapped to a NumPy array. The same happens if
>>>>>>> you do:
>>>>>>>
>>>>>>>     coord = UnitCube(1, 1, 1).coordinates()
>>>>>>>
>>>>>>> We can solve this by making adding a deepcopy argument to the
>>>>>>> coordinates (and the cells) methods, which defaults to True. This might
>>>>>>> break code as some (hmpfr me!) uses coordinates to move a mesh. Changing
>>>>>>> the coordinates when a deepcopy is done wont change the  coordinates on
>>>>>>> the mesh.
>>>>>>>
>>>>>>> Another way, is to add a reference of the mesh object to the NumPy
>>>>>>> array. Then the Mesh object wont go out of scope until the coordinate
>>>>>>> array does. The latter is more tricky to implement, but would be the
>>>>>>> best way to do it.
>>>>>>>
>>>>>>
>>>>>> We've seen before that keeping references to large objects will at
>>>>>> some point cause a memory leak.
>>>>>
>>>>> The problem with the previous reference storage, was that it was a
>>>>> global storage. My suggestions is a storage connected to an object. If
>>>>> an object returns data it owes it should make sure it does not destroy
>>>>> it for the lifetime of that data. We already do this for for example
>>>>> down_casting of objects.
>>
>> The underlying problem is the same, and I think a poor design. I would
>> not expect that a simple array object would keep a possibly very large
>> object in scope.
>
> The whole slicing interface of NumPy is built on this design.
>
>   a = np.array([some big list])
>   b = a[::2]
>
> b is here a view into a. Changing values in b will change values in a.
> If a is destroyed b is still working, as it knows it came from a and
> therefor keeps a reference to it. NumPy provides the interface, we need
> the functionality, why not use it?
>

The above argument doesn't make sense. How is it relevant? A DOLFIN
Vector is essentially an array, with some small and constant in size
extra data. A Mesh comes with numerous arrays, which would all be
kept in scope for the sake of one array.

> The coordinates of a mesh is not a small array. Nor are other examples
> which are natural to talk about in this regard.
>
>   MeshFunction.array, GenericFoo.data, Mesh.cells, Mesh.coordinates
>

Whether coordinates is small or not is irrelevant. The problem is that
the other objects associated with a Mesh are of the same order in
size.

>>>>>> I suggest that we
>>>>>>
>>>>>> 1. Have  one function that returns a copy of the coordinates and one
>>>>>> that provides view
>>>>>
>>>>> I tend to agree, but the question is what should be default. The most
>>>>> conservative would be to let a copy be default, but that would most
>>>>> probably break user code, as most people (...

[Johan Hake (johan-hake)]

On 03/28/2012 02:02 PM, Garth Wells wrote:
> On 28 March 2012 19:31, Johan Hake<email address hidden> wrote:
>> On 03/28/2012 12:31 PM, Garth Wells wrote:
>>> On 28 March 2012 16:42, Johan Hake<email address hidden> wrote:
>>>> On 03/28/2012 10:26 AM, Anders Logg wrote:
>>>>> On Wed, Mar 28, 2012 at 08:06:02AM -0000, Johan Hake wrote:
>>>>>> On 03/28/2012 09:41 AM, Garth Wells wrote:
>>>>>>> On 28 March 2012 14:03, Johan Hake<email address hidden> wrote:
>>>>>>>> This is because the mesh has gone out of scope and destroyed the
>>>>>>>> coordinatedata which was wrapped to a NumPy array. The same happens if
>>>>>>>> you do:
>>>>>>>>
>>>>>>>> coord = UnitCube(1, 1, 1).coordinates()
>>>>>>>>
>>>>>>>> We can solve this by making adding a deepcopy argument to the
>>>>>>>> coordinates (and the cells) methods, which defaults to True. This might
>>>>>>>> break code as some (hmpfr me!) uses coordinates to move a mesh. Changing
>>>>>>>> the coordinates when a deepcopy is done wont change the coordinates on
>>>>>>>> the mesh.
>>>>>>>>
>>>>>>>> Another way, is to add a reference of the mesh object to the NumPy
>>>>>>>> array. Then the Mesh object wont go out of scope until the coordinate
>>>>>>>> array does. The latter is more tricky to implement, but would be the
>>>>>>>> best way to do it.
>>>>>>>>
>>>>>>>
>>>>>>> We've seen before that keeping references to large objects will at
>>>>>>> some point cause a memory leak.
>>>>>>
>>>>>> The problem with the previous reference storage, was that it was a
>>>>>> global storage. My suggestions is a storage connected to an object. If
>>>>>> an object returns data it owes it should make sure it does not destroy
>>>>>> it for the lifetime of that data. We already do this for for example
>>>>>> down_casting of objects.
>>>
>>> The underlying problem is the same, and I think a poor design. I would
>>> not expect that a simple array object would keep a possibly very large
>>> object in scope.
>>
>> The whole slicing interface of NumPy is built on this design.
>>
>> a = np.array([some big list])
>> b = a[::2]
>>
>> b is here a view into a. Changing values in b will change values in a.
>> If a is destroyed b is still working, as it knows it came from a and
>> therefor keeps a reference to it. NumPy provides the interface, we need
>> the functionality, why not use it?
>>
>
> The above argument doesn't make sense. How is it relevant? A DOLFIN
> Vector is essentially an array, with some small and constant in size
> extra data. A Mesh comes with numerous arrays, which would all be
> kept in scope for the sake of one array.

This has nothing to do with dolfin Vectors, but on cashing large objects
to potentially small NumPy arrays. In the above example that large
object is another NumPy array.

>> The coordinates of a mesh is not a small array. Nor are other examples
>> which are natural to talk about in this regard.
>>
>> MeshFunction.array, GenericFoo.data, Mesh.cells, Mesh.coordinates
>>
>
> Whether coordinates is small or not is irrelevant. The problem is that
> the other objects associated with a Mesh are of the same order in
> size.

Not sure i see any problems here. Remember ...

[Garth Wells (garth-wells)]

On 28 March 2012 19:58, Anders Logg <email address hidden> wrote:
> On Wed, Mar 28, 2012 at 11:31:35AM -0000, Johan Hake wrote:
>> On 03/28/2012 12:31 PM, Garth Wells wrote:
>> > On 28 March 2012 16:42, Johan Hake<email address hidden>  wrote:
>> >> On 03/28/2012 10:26 AM, Anders Logg wrote:
>> >>> On Wed, Mar 28, 2012 at 08:06:02AM -0000, Johan Hake wrote:
>> >>>> On 03/28/2012 09:41 AM, Garth Wells wrote:
>> >>>>> On 28 March 2012 14:03, Johan Hake<email address hidden>     wrote:
>> >>>>>> This is because the mesh has gone out of scope and destroyed the
>> >>>>>> coordinatedata which was wrapped to a NumPy array. The same happens if
>> >>>>>> you do:
>> >>>>>>
>> >>>>>>     coord = UnitCube(1, 1, 1).coordinates()
>> >>>>>>
>> >>>>>> We can solve this by making adding a deepcopy argument to the
>> >>>>>> coordinates (and the cells) methods, which defaults to True. This might
>> >>>>>> break code as some (hmpfr me!) uses coordinates to move a mesh. Changing
>> >>>>>> the coordinates when a deepcopy is done wont change the  coordinates on
>> >>>>>> the mesh.
>> >>>>>>
>> >>>>>> Another way, is to add a reference of the mesh object to the NumPy
>> >>>>>> array. Then the Mesh object wont go out of scope until the coordinate
>> >>>>>> array does. The latter is more tricky to implement, but would be the
>> >>>>>> best way to do it.
>> >>>>>>
>> >>>>>
>> >>>>> We've seen before that keeping references to large objects will at
>> >>>>> some point cause a memory leak.
>> >>>>
>> >>>> The problem with the previous reference storage, was that it was a
>> >>>> global storage. My suggestions is a storage connected to an object. If
>> >>>> an object returns data it owes it should make sure it does not destroy
>> >>>> it for the lifetime of that data. We already do this for for example
>> >>>> down_casting of objects.
>> >
>> > The underlying problem is the same, and I think a poor design. I would
>> > not expect that a simple array object would keep a possibly very large
>> > object in scope.
>>
>> The whole slicing interface of NumPy is built on this design.
>>
>>    a = np.array([some big list])
>>    b = a[::2]
>>
>> b is here a view into a. Changing values in b will change values in a.
>> If a is destroyed b is still working, as it knows it came from a and
>> therefor keeps a reference to it. NumPy provides the interface, we need
>> the functionality, why not use it?
>>
>> The coordinates of a mesh is not a small array. Nor are other examples
>> which are natural to talk about in this regard.
>>
>>    MeshFunction.array, GenericFoo.data, Mesh.cells, Mesh.coordinates
>
> I think Garth meant that the Mesh is a big data structure which would
> then be attached to the array.
>

Yes.

> I don't see that as a big problem. The coordinates array is itself of
> considerable size compared to the Mesh, and if someone calls
> Mesh.coordinates() they will likely assume that the mesh is still
> alive.

 I don't think so. From

   coord = UnitCube(1, 1, 1).coordinates()

I would not expect the Mesh to be in scope.

> If not, they could call Mesh.coordinates().copy().
>

We shouldn't make this assumption on what a user wants, and we
shouldn't assume that the obj...

[Garth Wells (garth-wells)]

On 28 March 2012 20:39, Johan Hake <email address hidden> wrote:
> On 03/28/2012 02:02 PM, Garth Wells wrote:
>> On 28 March 2012 19:31, Johan Hake<email address hidden>  wrote:
>>> On 03/28/2012 12:31 PM, Garth Wells wrote:
>>>> On 28 March 2012 16:42, Johan Hake<email address hidden>    wrote:
>>>>> On 03/28/2012 10:26 AM, Anders Logg wrote:
>>>>>> On Wed, Mar 28, 2012 at 08:06:02AM -0000, Johan Hake wrote:
>>>>>>> On 03/28/2012 09:41 AM, Garth Wells wrote:
>>>>>>>> On 28 March 2012 14:03, Johan Hake<email address hidden>       wrote:
>>>>>>>>> This is because the mesh has gone out of scope and destroyed the
>>>>>>>>> coordinatedata which was wrapped to a NumPy array. The same happens if
>>>>>>>>> you do:
>>>>>>>>>
>>>>>>>>>      coord = UnitCube(1, 1, 1).coordinates()
>>>>>>>>>
>>>>>>>>> We can solve this by making adding a deepcopy argument to the
>>>>>>>>> coordinates (and the cells) methods, which defaults to True. This might
>>>>>>>>> break code as some (hmpfr me!) uses coordinates to move a mesh. Changing
>>>>>>>>> the coordinates when a deepcopy is done wont change the  coordinates on
>>>>>>>>> the mesh.
>>>>>>>>>
>>>>>>>>> Another way, is to add a reference of the mesh object to the NumPy
>>>>>>>>> array. Then the Mesh object wont go out of scope until the coordinate
>>>>>>>>> array does. The latter is more tricky to implement, but would be the
>>>>>>>>> best way to do it.
>>>>>>>>>
>>>>>>>>
>>>>>>>> We've seen before that keeping references to large objects will at
>>>>>>>> some point cause a memory leak.
>>>>>>>
>>>>>>> The problem with the previous reference storage, was that it was a
>>>>>>> global storage. My suggestions is a storage connected to an object. If
>>>>>>> an object returns data it owes it should make sure it does not destroy
>>>>>>> it for the lifetime of that data. We already do this for for example
>>>>>>> down_casting of objects.
>>>>
>>>> The underlying problem is the same, and I think a poor design. I would
>>>> not expect that a simple array object would keep a possibly very large
>>>> object in scope.
>>>
>>> The whole slicing interface of NumPy is built on this design.
>>>
>>>    a = np.array([some big list])
>>>    b = a[::2]
>>>
>>> b is here a view into a. Changing values in b will change values in a.
>>> If a is destroyed b is still working, as it knows it came from a and
>>> therefor keeps a reference to it. NumPy provides the interface, we need
>>> the functionality, why not use it?
>>>
>>
>> The above argument doesn't make sense. How is it relevant? A DOLFIN
>> Vector is essentially an array, with some small and constant in size
>> extra data.  A Mesh comes with numerous arrays, which would all be
>> kept in scope for the sake of one array.
>
> This has nothing to do with dolfin Vectors, but on cashing large objects
> to potentially small NumPy arrays. In the above example that large
> object is another NumPy array.
>
>>> The coordinates of a mesh is not a small array. Nor are other examples
>>> which are natural to talk about in this regard.
>>>
>>>    MeshFunction.array, GenericFoo.data, Mesh.cells, Mesh.coordinates
>>>
>>
>> Whether coordinates is small or not is irrelevant. Th...

[Anders Logg (logg)]

On Wed, Mar 28, 2012 at 01:11:33PM -0000, Garth Wells wrote:

> > I don't see that as a big problem. The coordinates array is itself of
> > considerable size compared to the Mesh, and if someone calls
> > Mesh.coordinates() they will likely assume that the mesh is still
> > alive.
>
> I don't think so. From
>
> coord = UnitCube(1, 1, 1).coordinates()
>
> I would not expect the Mesh to be in scope.

Me neither, I would expect coord to be a copy.

But if I do

  mesh = UnitCube(1, 1, 1)
  coord = mesh.coordinates()

I expect to be able to modify the mesh coordinates by modifying coord.

> > If not, they could call Mesh.coordinates().copy().
> >
>
> We shouldn't make this assumption on what a user wants, and we
> shouldn't assume that the objects stored in a Mesh won't change. Say,
> for example, the Mesh data structure is changed in the future that it
> stores a hierarchy of refinements. Then by accessing the array of
> coordinates, one is, behind the scenes and unintuitively, keeping a
> possibly very hierarchy of objects in scope.
>
> We have had cases in past where we artificially keep objects in scope,
> and this has lead to unanticipated problems. This is why I believe
> that it's a flawed design that should be avoided.

I think this is something of a special case. The mesh.coordinates()
thing has been around for quite some time, it is very convenient and
I think it's being used quite a lot (I use it).

--
Anders

> Garth
>
> >
> >
> >> >>>>> I suggest that we
> >> >>>>>
> >> >>>>> 1. Have  one function that returns a copy of the coordinates and one
> >> >>>>> that provides view
> >> >>>>
> >> >>>> I tend to agree, but the question is what should be default. The most
> >> >>>> conservative would be to let a copy be default, but that would most
> >> >>>> probably break user code, as most people (hmpfr me) uses coordinates to
> >> >>>> move meshes.
> >> >>>
> >> >>> I think we should try hard to avoid that. Accessing coordinates() by
> >> >>> reference (not copy) has been the default behavior and we would break
> >> >>> user code, possibly also some code examples from the book.
> >> >>>
> >> >>>>> or
> >> >>>>>
> >> >>>>> 2. Make the coordinate data a boost::shared_array<double>     or a
> >> >>>>> boost::shared_ptr<std::vector<double>     >     object.
> >> >>>>
> >> >>>> This would introduce a somewhat middle solution, which is pretty
> >> >>>> cumbersome to implement and which need to be changed anyhow, when we
> >> >>>> revamp the mesh interface.
> >> >>>
> >> >
> >> > It should be easy to implement - I have already changed some plain
> >> > arrays to std::vector. Relying on the contiguous storage that a
> >> > std::vector provides, it's easy to interface with functions looking to
> >> > accept a plain pointer. The question is how it should work in the SWIG
> >> > layer.
> >>
> >> Exactly, and here NumPy provides an excellent functionality. For example
> >> if we a provide a shared_ptr version of std::vector we might be able to
> >> do things in the way you suggest, but then we are out of NumPy land.
> >>
> >> A middle way might be to create a dummy Python object which holds a
> >> reference to a shared_ptr vector, which we then cre...

[Garth Wells (garth-wells)]

On 29 March 2012 04:44, Anders Logg <email address hidden> wrote:
> On Wed, Mar 28, 2012 at 01:11:33PM -0000, Garth Wells wrote:
>
>> > I don't see that as a big problem. The coordinates array is itself of
>> > considerable size compared to the Mesh, and if someone calls
>> > Mesh.coordinates() they will likely assume that the mesh is still
>> > alive.
>>
>>  I don't think so. From
>>
>>    coord = UnitCube(1, 1, 1).coordinates()
>>
>> I would not expect the Mesh to be in scope.
>
> Me neither, I would expect coord to be a copy.
>

Agree.

> But if I do
>
>  mesh = UnitCube(1, 1, 1)
>  coord = mesh.coordinates()
>
> I expect to be able to modify the mesh coordinates by modifying coord.
>

Agree.

>> > If not, they could call Mesh.coordinates().copy().
>> >
>>
>> We shouldn't make this assumption on what a user wants, and we
>> shouldn't assume that the objects stored in a Mesh won't change. Say,
>> for example, the Mesh data structure is changed in the future that it
>> stores a hierarchy of refinements. Then by accessing the array of
>> coordinates, one is, behind the scenes and unintuitively, keeping a
>> possibly very hierarchy of objects in scope.
>>
>> We have had cases in past where we artificially keep objects in scope,
>> and this has lead to unanticipated problems. This is why I believe
>> that it's a flawed design that should be avoided.
>
> I think this is something of a special case. The mesh.coordinates()
> thing has been around for quite some time, it is very convenient and
> I think it's being used quite a lot (I use it).
>

I don't see a problem with the current design on the C++ side. What I
see as a flawed design is the proposal that the Python array that
wraps the coordinate array should also hold a reference to the Mesh.
In the case

   coord = UnitCube(1, 1, 1).coordinates()

the Mesh should go out of scope and coord should just be an array.
We've dealt with this issue elsewhere using shared pointers to member
data.

Garth

> --
> Anders
>
>
>> Garth
>>
>> >
>> >
>> >> >>>>> I suggest that we
>> >> >>>>>
>> >> >>>>> 1. Have  one function that returns a copy of the coordinates and one
>> >> >>>>> that provides view
>> >> >>>>
>> >> >>>> I tend to agree, but the question is what should be default. The most
>> >> >>>> conservative would be to let a copy be default, but that would most
>> >> >>>> probably break user code, as most people (hmpfr me) uses coordinates to
>> >> >>>> move meshes.
>> >> >>>
>> >> >>> I think we should try hard to avoid that. Accessing coordinates() by
>> >> >>> reference (not copy) has been the default behavior and we would break
>> >> >>> user code, possibly also some code examples from the book.
>> >> >>>
>> >> >>>>> or
>> >> >>>>>
>> >> >>>>> 2. Make the coordinate data a boost::shared_array<double>     or a
>> >> >>>>> boost::shared_ptr<std::vector<double>     >     object.
>> >> >>>>
>> >> >>>> This would introduce a somewhat middle solution, which is pretty
>> >> >>>> cumbersome to implement and which need to be changed anyhow, when we
>> >> >>>> revamp the mesh interface.
>> >> >>>
>> >> >
>> >> > It should be easy to implement - I have already changed some plain
>> >> > arrays to std::vector. Rel...

[Johan Hake (johan-hake)]

On 03/29/2012 01:53 AM, Garth Wells wrote:
> On 29 March 2012 04:44, Anders Logg<email address hidden> wrote:
>> On Wed, Mar 28, 2012 at 01:11:33PM -0000, Garth Wells wrote:
>>
>>>> I don't see that as a big problem. The coordinates array is itself of
>>>> considerable size compared to the Mesh, and if someone calls
>>>> Mesh.coordinates() they will likely assume that the mesh is still
>>>> alive.
>>>
>>> I don't think so. From
>>>
>>> coord = UnitCube(1, 1, 1).coordinates()
>>>
>>> I would not expect the Mesh to be in scope.
>>
>> Me neither, I would expect coord to be a copy.
>>
>
> Agree.
>
>> But if I do
>>
>> mesh = UnitCube(1, 1, 1)
>> coord = mesh.coordinates()
>>
>> I expect to be able to modify the mesh coordinates by modifying coord.
>>
>
> Agree.
>
>>>> If not, they could call Mesh.coordinates().copy().
>>>>
>>>
>>> We shouldn't make this assumption on what a user wants, and we
>>> shouldn't assume that the objects stored in a Mesh won't change. Say,
>>> for example, the Mesh data structure is changed in the future that it
>>> stores a hierarchy of refinements. Then by accessing the array of
>>> coordinates, one is, behind the scenes and unintuitively, keeping a
>>> possibly very hierarchy of objects in scope.
>>>
>>> We have had cases in past where we artificially keep objects in scope,
>>> and this has lead to unanticipated problems. This is why I believe
>>> that it's a flawed design that should be avoided.
>>
>> I think this is something of a special case. The mesh.coordinates()
>> thing has been around for quite some time, it is very convenient and
>> I think it's being used quite a lot (I use it).
>>
>
> I don't see a problem with the current design on the C++ side. What I
> see as a flawed design is the proposal that the Python array that
> wraps the coordinate array should also hold a reference to the Mesh.
> In the case
>
> coord = UnitCube(1, 1, 1).coordinates()
>
> the Mesh should go out of scope and coord should just be an array.
> We've dealt with this issue elsewhere using shared pointers to member
> data.

See your point Garth. The solution can then be to store the mesh data in
a shared_ptr<std::vector<double> > and shared_ptr<std::vector<uint> >.
We can add:

   shared_coordinates and shared_cells

to the C++ interface and map that to coordinates and cells in the Python
interface. Then we need to add some magic that turns the shared_ptr
stored data into a NumPy array. This can be done but the C++ interface
need to be changed first. What I mean with this can be done is that the
std::vector will be wrapped into a Python object, and a NumPy array will
be constructed wrapping the data of the std::vector. The Python wrapper
will then be attached to the NumPy array, keeping that std::vector
together with its wrapper in scope.

What with the other examples:

   MeshFunction.array, GenericFoo.data

Is it possible to come with a similar solution there? It looks like we
can do it for MeshFunction, and GenericVector (uBLASVector, MTL4Vector)
not sure we can do it for GenericMatrix (uBLASSparseMatrix, MTL4Matrix),
as the underlying sparse matrices do not store their data using shared_ptr.

Johan

> Garth
>...

[Garth Wells (garth-wells)]

On 29 March 2012 15:55, Johan Hake <email address hidden> wrote:
> On 03/29/2012 01:53 AM, Garth Wells wrote:
>> On 29 March 2012 04:44, Anders Logg<email address hidden>  wrote:
>>> On Wed, Mar 28, 2012 at 01:11:33PM -0000, Garth Wells wrote:
>>>
>>>>> I don't see that as a big problem. The coordinates array is itself of
>>>>> considerable size compared to the Mesh, and if someone calls
>>>>> Mesh.coordinates() they will likely assume that the mesh is still
>>>>> alive.
>>>>
>>>>   I don't think so. From
>>>>
>>>>     coord = UnitCube(1, 1, 1).coordinates()
>>>>
>>>> I would not expect the Mesh to be in scope.
>>>
>>> Me neither, I would expect coord to be a copy.
>>>
>>
>> Agree.
>>
>>> But if I do
>>>
>>>   mesh = UnitCube(1, 1, 1)
>>>   coord = mesh.coordinates()
>>>
>>> I expect to be able to modify the mesh coordinates by modifying coord.
>>>
>>
>> Agree.
>>
>>>>> If not, they could call Mesh.coordinates().copy().
>>>>>
>>>>
>>>> We shouldn't make this assumption on what a user wants, and we
>>>> shouldn't assume that the objects stored in a Mesh won't change. Say,
>>>> for example, the Mesh data structure is changed in the future that it
>>>> stores a hierarchy of refinements. Then by accessing the array of
>>>> coordinates, one is, behind the scenes and unintuitively, keeping a
>>>> possibly very hierarchy of objects in scope.
>>>>
>>>> We have had cases in past where we artificially keep objects in scope,
>>>> and this has lead to unanticipated problems. This is why I believe
>>>> that it's a flawed design that should be avoided.
>>>
>>> I think this is something of a special case. The mesh.coordinates()
>>> thing has been around for quite some time, it is very convenient and
>>> I think it's being used quite a lot (I use it).
>>>
>>
>> I don't see a problem with the current design on the C++ side.  What I
>> see as a flawed design is the proposal that the Python array that
>> wraps the coordinate array should also hold a reference to the Mesh.
>> In the case
>>
>>     coord = UnitCube(1, 1, 1).coordinates()
>>
>> the Mesh should go out of scope and coord should just be an array.
>> We've dealt with this issue elsewhere using shared pointers to member
>> data.
>
> See your point Garth. The solution can then be to store the mesh data in
> a shared_ptr<std::vector<double> > and shared_ptr<std::vector<uint> >.
> We can add:
>
>   shared_coordinates and shared_cells
>
> to the C++ interface and map that to coordinates and cells in the Python
> interface. Then we need to add some magic that turns the shared_ptr
> stored data into a NumPy array. This can be done but the C++ interface
> need to be changed first. What I mean with this can be done is that the
> std::vector will be wrapped into a Python object, and a NumPy array will
> be constructed wrapping the data of the std::vector. The Python wrapper
> will then be attached to the NumPy array, keeping that std::vector
> together with its wrapper in scope.
>
> What with the other examples:
>
>   MeshFunction.array, GenericFoo.data
>
> Is it possible to come with a similar solution there? It looks like we
> can do it for MeshFunction, and GenericVector (uBLASVector, MTL4Vector)
> not sure w...

[Johan Hake (johan-hake)]

The solution with shared_ptr of std::vectors is something we need to wait for. I have a tested and running implementation of the mesh attached array, discussed above. If you rather not see that solution pushed into trunk please speak up. Note the proposed solution of attaching the object to the returned numpy array, solves the same problem for data structures of MTL4 and uBLAS objects, which cannot be solved with a shared_ptr of std::vector data discussed above.

I could push it as some of the code can be reused for the solution using shared_ptr of std::vector and then open another bug for the storage of mesh data?

[Garth Wells (garth-wells)]

On 3 April 2012 20:43, Johan Hake <email address hidden> wrote:
> The solution with shared_ptr of std::vectors is something we need to
> wait for. I have a tested and running implementation of the mesh
> attached array, discussed above. If you rather not see that solution
> pushed into trunk please speak up.

Is the the mesh attached to the array? If it is, I think it's already
clear that I don't like it ;).

> Note the proposed solution of
> attaching the object to the returned numpy array, solves the same
> problem for data structures of MTL4 and uBLAS objects, which cannot be
> solved with a shared_ptr of std::vector data discussed above.
>

I'd say that the 'problem' in the above cases is that we expose the
data in the first place. It's a very low level interface, so I'm happy
to leave these case as is.

Garth

> I could push it as some of the code can be reused for the solution using
> shared_ptr of std::vector and then open another bug for the storage of
> mesh data?
>
> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/966553
>
> Title:
>  BoundaryMesh.coordinates () gives random results
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/966553/+subscriptions

[Anders Logg (logg)]

On Tue, Apr 03, 2012 at 12:43:23PM -0000, Johan Hake wrote:
> The solution with shared_ptr of std::vectors is something we need to
> wait for. I have a tested and running implementation of the mesh
> attached array, discussed above. If you rather not see that solution
> pushed into trunk please speak up. Note the proposed solution of
> attaching the object to the returned numpy array, solves the same
> problem for data structures of MTL4 and uBLAS objects, which cannot be
> solved with a shared_ptr of std::vector data discussed above.
>
> I could push it as some of the code can be reused for the solution using
> shared_ptr of std::vector and then open another bug for the storage of
> mesh data?

I think this solution is ok.

--
Anders

[Anders Logg (logg)]

On Tue, Apr 03, 2012 at 01:03:48PM -0000, Garth Wells wrote:
> On 3 April 2012 20:43, Johan Hake <email address hidden> wrote:
> > The solution with shared_ptr of std::vectors is something we need to
> > wait for. I have a tested and running implementation of the mesh
> > attached array, discussed above. If you rather not see that solution
> > pushed into trunk please speak up.
>
> Is the the mesh attached to the array? If it is, I think it's already
> clear that I don't like it ;).
>
> > Note the proposed solution of
> > attaching the object to the returned numpy array, solves the same
> > problem for data structures of MTL4 and uBLAS objects, which cannot be
> > solved with a shared_ptr of std::vector data discussed above.
> >
>
> I'd say that the 'problem' in the above cases is that we expose the
> data in the first place. It's a very low level interface, so I'm happy
> to leave these case as is.

I think we should allow for low-level access of array data, not always
but when it makes sense. It is particularly useful from Python to work
with NumPy arrays as basic containers instead of being forced to work
through an object-oriented interface that hides the data. OO hiding of
data is a good concept but it can easily go too far.

--
Anders

> Garth
>
> > I could push it as some of the code can be reused for the solution using
> > shared_ptr of std::vector and then open another bug for the storage of
> > mesh data?
> >
> >
> > Title:
> >  BoundaryMesh.coordinates () gives random results
> >
> > To manage notifications about this bug go to:
> > https://bugs.launchpad.net/dolfin/+bug/966553/+subscriptions
>

[Johan Hake (johan-hake)]

And we have a conclusion! Not... ;)

On 04/03/2012 03:03 PM, Garth Wells wrote:
> On 3 April 2012 20:43, Johan Hake<email address hidden> wrote:
>> The solution with shared_ptr of std::vectors is something we need to
>> wait for. I have a tested and running implementation of the mesh
>> attached array, discussed above. If you rather not see that solution
>> pushed into trunk please speak up.
>
> Is the the mesh attached to the array? If it is, I think it's already
> clear that I don't like it ;).

Sure, and therefore I will open another bug addressing this issue. But
as mentioned that is a larger refactoring.

>> Note the proposed solution of
>> attaching the object to the returned numpy array, solves the same
>> problem for data structures of MTL4 and uBLAS objects, which cannot be
>> solved with a shared_ptr of std::vector data discussed above.
>>
>
> I'd say that the 'problem' in the above cases is that we expose the
> data in the first place. It's a very low level interface, so I'm happy
> to leave these case as is.

I disagree. Exposing these data types is a nice way of getting access to
low lever data structure for free. With the proposed solution the
original data will not go out of scope and the numpy array attached
object only add a small overhead compared to the numpy array as they
share the actual data. This is not the case for the Mesh, which
potentially comes with a lot of overhead, which I agree is suboptimal.

Johan

> Garth
>
>> I could push it as some of the code can be reused for the solution using
>> shared_ptr of std::vector and then open another bug for the storage of
>> mesh data?
>>
>> --
>> You received this bug notification because you are a member of DOLFIN
>> Core Team, which is subscribed to DOLFIN.
>> https://bugs.launchpad.net/bugs/966553
>>
>> Title:
>> BoundaryMesh.coordinates () gives random results
>>
>> To manage notifications about this bug go to:
>> https://bugs.launchpad.net/dolfin/+bug/966553/+subscriptions
>

[Garth Wells (garth-wells)]

On 3 April 2012 21:25, Johan Hake <email address hidden> wrote:
> And we have a conclusion! Not... ;)
>
> On 04/03/2012 03:03 PM, Garth Wells wrote:
>> On 3 April 2012 20:43, Johan Hake<email address hidden>  wrote:
>>> The solution with shared_ptr of std::vectors is something we need to
>>> wait for. I have a tested and running implementation of the mesh
>>> attached array, discussed above. If you rather not see that solution
>>> pushed into trunk please speak up.
>>
>> Is the the mesh attached to the array? If it is, I think it's already
>> clear that I don't like it ;).
>
> Sure, and therefore I will open another bug addressing this issue. But
> as mentioned that is a larger refactoring.
>
>>> Note the proposed solution of
>>> attaching the object to the returned numpy array, solves the same
>>> problem for data structures of MTL4 and uBLAS objects, which cannot be
>>> solved with a shared_ptr of std::vector data discussed above.
>>>
>>
>> I'd say that the 'problem' in the above cases is that we expose the
>> data in the first place. It's a very low level interface, so I'm happy
>> to leave  these case as is.
>
> I disagree. Exposing these data types is a nice way of getting access to
> low lever data structure for free. With the proposed solution the
> original data will not go out of scope and the numpy array attached
> object only add a small overhead compared to the numpy array as they
> share the actual data.

Perhaps in this case, but it goes deeper. For example, with a sparse
matrix keeping the data in scope may involve keeping an entire LU
factorisation in scope behind the scenes. It is a very low level
interface, which is why I can live with the user who wants to access
such low-data taking some of the responsibility for managing it.

Garth

> This is not the case for the Mesh, which
> potentially comes with a lot of overhead, which I agree is suboptimal.
>
> Johan
>
>> Garth
>>
>>> I could push it as some of the code can be reused for the solution using
>>> shared_ptr of std::vector and then open another bug for the storage of
>>> mesh data?
>>>
>>> --
>>> You received this bug notification because you are a member of DOLFIN
>>> Core Team, which is subscribed to DOLFIN.
>>> https://bugs.launchpad.net/bugs/966553
>>>
>>> Title:
>>>   BoundaryMesh.coordinates () gives random results
>>>
>>> To manage notifications about this bug go to:
>>> https://bugs.launchpad.net/dolfin/+bug/966553/+subscriptions
>>
>
> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/966553
>
> Title:
>  BoundaryMesh.coordinates () gives random results
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/966553/+subscriptions

[Johan Hake (johan-hake)]

Ok, will push the fix and report a bug once I am back from easter holiday.

Johan

On 04/03/2012 04:51 PM, Garth Wells wrote:
> On 3 April 2012 21:25, Johan Hake<email address hidden> wrote:
>> And we have a conclusion! Not... ;)
>>
>> On 04/03/2012 03:03 PM, Garth Wells wrote:
>>> On 3 April 2012 20:43, Johan Hake<email address hidden> wrote:
>>>> The solution with shared_ptr of std::vectors is something we need to
>>>> wait for. I have a tested and running implementation of the mesh
>>>> attached array, discussed above. If you rather not see that solution
>>>> pushed into trunk please speak up.
>>>
>>> Is the the mesh attached to the array? If it is, I think it's already
>>> clear that I don't like it ;).
>>
>> Sure, and therefore I will open another bug addressing this issue. But
>> as mentioned that is a larger refactoring.
>>
>>>> Note the proposed solution of
>>>> attaching the object to the returned numpy array, solves the same
>>>> problem for data structures of MTL4 and uBLAS objects, which cannot be
>>>> solved with a shared_ptr of std::vector data discussed above.
>>>>
>>>
>>> I'd say that the 'problem' in the above cases is that we expose the
>>> data in the first place. It's a very low level interface, so I'm happy
>>> to leave these case as is.
>>
>> I disagree. Exposing these data types is a nice way of getting access to
>> low lever data structure for free. With the proposed solution the
>> original data will not go out of scope and the numpy array attached
>> object only add a small overhead compared to the numpy array as they
>> share the actual data.
>
> Perhaps in this case, but it goes deeper. For example, with a sparse
> matrix keeping the data in scope may involve keeping an entire LU
> factorisation in scope behind the scenes. It is a very low level
> interface, which is why I can live with the user who wants to access
> such low-data taking some of the responsibility for managing it.
>
> Garth
>
>> This is not the case for the Mesh, which
>> potentially comes with a lot of overhead, which I agree is suboptimal.
>>
>> Johan
>>
>>> Garth
>>>
>>>> I could push it as some of the code can be reused for the solution using
>>>> shared_ptr of std::vector and then open another bug for the storage of
>>>> mesh data?
>>>>
>>>> --
>>>> You received this bug notification because you are a member of DOLFIN
>>>> Core Team, which is subscribed to DOLFIN.
>>>> https://bugs.launchpad.net/bugs/966553
>>>>
>>>> Title:
>>>> BoundaryMesh.coordinates () gives random results
>>>>
>>>> To manage notifications about this bug go to:
>>>> https://bugs.launchpad.net/dolfin/+bug/966553/+subscriptions
>>>
>>
>> --
>> You received this bug notification because you are a member of DOLFIN
>> Core Team, which is subscribed to DOLFIN.
>> https://bugs.launchpad.net/bugs/966553
>>
>> Title:
>> BoundaryMesh.coordinates () gives random results
>>
>> To manage notifications about this bug go to:
>> https://bugs.launchpad.net/dolfin/+bug/966553/+subscriptions
>