Entries of array 'parent_vertex_indices' change values randomly when leaving the function where it was computed

[Maximilian Albert]

The attached script illustrates a bug where the entries of a certain array can randomly change their values when leaving the subroutine where the array was computed.

In the attached example I'm creating a UnitCubeMesh and extract the vertex indices of the submesh corresponding to the upper half. Inside the helper function where this computation happens the indices are correct, but when the resulting array is passed back to the caller, the entries of the array end up with random values.

Here is some sample output of the script:

Indices computed in subroutine: [32 33 37 53 49 48 36 52 34 38 54 50 35 39 55 51 41 57 40 56 42 58 43 59 45 61 44 60 46 62 47 63]

Indices in main: [54044368, 140140727019384, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2308706344189829152, 2314885530818453536, 2314885530818453536, 2314898793677463584, 2314885530818453536, 2318280822725025824, 2314885530818453536, 3184080258900959264, 2314885530818453536, 3467807035425300512, 2314885530818447916, 2314885530818453536, 3467807035425300512, 2314885530818453548, 2314885530818453536, 2314885530818456624, 2314885530818453536, 2314885530819244064, 2314885530818453536]

Note that the bug does not happen when the array is copied before returning it (as illustrated by calling the helper function with 'return_copy=True').

Also, interestingly it doesn't happen if the two meshes are plotted before the array is returned (set 'plot_meshes=True' to illustrate this).

I'm using the development version of dolfin. I tested it briefly with version 1.2.0 (with the old syntax for getting the 'parent_vertex_indices' array) and the bug doesn't seem to be present there.

[Johan Hake]

The parent_vertex_indices is a view into the c++ strored values. When the SubMesh is destroyed upon leaving the function the data the view points to is also destroyed and that is the random values you see. This is of course potentially dangerous, but you as a user needs to take measure to prevent this being dangerous, by forexample as you suggest, return a copy.

[Maximilian Albert]

Thanks for the quick reply, Johan! I suspected something along those lines, but missed the fact that the returned array actually belongs to the submesh and is thus destroyed with it. This might be obvious to someone more familiar with the mesh data structure, but I don't think it is clear from the documentation - especially since it seems to be inconsistent with, say, f.vector().array() (where f is a dolfin.Function), which if I am not mistaken does return a copy of the vector. It took me a good few hours to debug it this morning (partly because the issue was hidden under some other layers of code) so I was wondering whether it's worth always returning a copy to prevent something like this happening accidentally. On the other hand, there might be efficiency reasons to keep it the way it is. But I thought since this "bug" is not present in dolfin version 1.2 it might have slipped through. In any case, if you decide to keep it the way it is, would you mind adding a comment to the documentation that this returns a view into the underlying C++ array rather than a copy? Many thanks!

[Johan Hake]

You are correct that Foo.array is not consistently returning a view of some underlying data. Making it a bit more confusing is that for most places we add an array function where we would like to return a numpy array, and then it is easy to add notion about view and copy in the doc string. However for MesData.array it is different because the method is already there in the C++ interface and it is this method that returns a reference to the underlying data. The return array is then automatically converted to a numpy array via a SWIG typemap.

For us to change the docstring for that method we would need to make changes to the documentation in the C++ header file, which might or might not make sense as in C++ it is pretty clear that the array is a view as it is a reference.

This is a good example where mixing C++ and Python not always works as it should.

[Maximilian Albert]

Thanks again for the explanation, Johan. (And sorry for the delay in replying.) I see now that it's tricky to document it properly. Would it be possible to always return a copy then, or are changes in the code subject to similar limitations as in the documentation due to the automatic SWIG wrapping?

[Johan Hake]

Not automatically. One would need to implement this in the Python layer, then one can fix both documentation and return a copy.

Not sure we want to do that as the Mesh::data structure might change after the next release.

[Maximilian Albert]

I see. No worries then, and thanks for the clarifications!