UnitCubeMesh creates duplicate vertices

[Tormod Landet]

DOLFIN will, with the python pybind11 wrappers, create duplicated and wrong nodes for UnitCubeMesh, see commit cde22d1 and the CI test results

The previous unit tests only checked that UnitCubeMesh created the right number of vertices, but UnitCubeMesh(1, 1, 1) creates 8 vertices (right) at only four different locations (obviously wrong)

I have not yet looked into the pybind11 wrapper code to figure out how the wrappers could cause something like this. Could it be that the positions are reported wrongly and the underlying geometry is correct?

[Jan Blechta]

That's not ok. But I can't reproduce it. What is the reproducible way to get this behaviour?

[Tormod Landet]

Very strange, I will test with the latest docker image to see if the error is there. I got the above results with a build of DOLFIN from yesterday. Are you using pybind? I am, but it should not matter?

[Tormod Landet]

I tested now with a Docker image that was built this week based on the FEniCS dev_env image and with dolfin is installed with pybind11 instead of SWIG.

Same error as above

The Dockerfile is here and I also pushed it to Docker Hub:

docker pull trlandet/fenics-dev:py3_CI

It basically builds DOLFIN with -DDOLFIN_ENABLE_PYTHON=OFF and then runs pip installs in the python/ directory. Some extra packages like h5py are also installed in the Docker recipe, but I assume this is not relevant

Can anyone reproduce the bug?

[Tormod Landet]

CircleCI now replicates the bug

[Jan Blechta]

Apparently dangling pointer in

    const double* Point::coordinates() const
    { return _x.data(); }

py::class_<dolfin::Point>(m, "Point")

      .def("array", [](dolfin::Point& self)
           { return Eigen::Map<Eigen::Vector3d>(self.coordinates()); })

Note that Python array() does not use C++ array(). Pybind wrapper apparently creates a Python object referencing to a memory in temporary Point object.

>>> from dolfin import *
>>> mesh = UnitCubeMesh(1, 1, 1)
>>> for v in vertices(mesh): print(v.point()[:])
... 
[ 0.  0.  0.]
[ 1.  0.  0.]
[ 0.  1.  0.]
[ 1.  1.  0.]
[ 0.  0.  1.]
[ 1.  0.  1.]
[ 0.  1.  1.]
[ 1.  1.  1.]
>>> for v in vertices(mesh): print(v.point().array())
... 
[ 0.  0.  0.]
[ 0.  0.  0.]
[ 1.  1.  0.]
[ 1.  1.  0.]
[ 1.  0.  1.]
[ 1.  0.  1.]
[ 1.  1.  1.]
[ 1.  1.  1.]

(Another memory bug in Pybind11. It's easy to write memory-corrupting code in Pybind. We have to be more careful about shallow copies and write more tests.)

[Jan Blechta]

What did we agree before? That (in Python) Point.array() is a view or a copy? In C++ it is a copy contrary to Point::coordinates().

[Tormod Landet]

Thanks for debugging! I have not had time to look at this yet

I believe we agreed that array_view() is a view and array() is a copy. For EigenMatrix we also have data() and data_view() where data is a copy and the view is a writeable view

[Prof Garth Wells]

Can we just avoid exposing Point.coordinates()?

[Jan Blechta]

Fix here

diff --git a/python/src/geometry.cpp b/python/src/geometry.cpp
index 5af9fe7c0..f4b3d7002 100644
--- a/python/src/geometry.cpp
+++ b/python/src/geometry.cpp
@@ -121,8 +121,9 @@ namespace dolfin_wrappers
       .def(py::self == py::self)
       .def(py::self * float())
       .def(py::self / float())
-      .def("array", [](dolfin::Point& self)
-           { return Eigen::Map<Eigen::Vector3d>(self.coordinates()); })
+      .def("array",
+           [](dolfin::Point& self) { return Eigen::Vector3d(self.coordinates()); },
+           "Return copy of coordinate array")
       .def("norm", &dolfin::Point::norm)
       .def("x", &dolfin::Point::x)
       .def("y", &dolfin::Point::y)

Test here

p = Point(1, 2, 3)
p.array()[:] *= 0
assert all(p.array() != 0)
p[:] *= 0
assert all(p.array() == 0)

[Tormod Landet]

Maybe, but then indexing point should be better documented

Edit: we should just expose .coordinates() directly (returns a Python list) and then convert std::array to Eigen::Vector3d later in dolfin::Point which will make it return a numpy array Edit2: sorry, I see that coordinates returns double*

[Jan Blechta]

The essence of the bug is

let_us_take_reference_to_plain_ptr = Eigen::Map<Eigen::Vector3d>(self.coordinates());

[Jan Blechta]

So the problem is not using self.coordinates() but using it wrongly.

[Prof Garth Wells]

The point is do we need more than one way to do the same thing?

[Jan Blechta]

Maybe, but then indexing point should be better documented.

I've added docstring noting it is a copy.

p[:] boils down do getitem (gets a copy if used on rhs) and setitem (sets if used on lhs).

[Jan Blechta]

The point is do we need more than one way to do the same thing?

We are not having more than one way. Point.array() is a copy in both C++ and Python. Then there is Point.coordinates() in C++ and Point.__(gs)etitem__ in Python allowing some shallow manipulations.

[Prof Garth Wells]

@blechta We do have a whole bunch of ways to access and set coordinates of a Point. Why do we need so many?

[Jan Blechta]

Let me remind you that the current state of both C++ and Python interface to Point is a result of very long recent discussion. I was for not adding any array() and only having coordinates() on both sides, but I was overvoted. See https://bitbucket.org/fenics-project/dolfin/pull-requests/328/throw-error-for-index-out-of-range-and-add/activity

If you want to change things, I suggest to open a PR. This is about fixing an obvious bug which I proposed a fix and test, see above.

[Jan Blechta]

Btw, we only have one way of setting a coordinate, which is Point.__setitem__. Getting is possible through p.x() (stupid), p[0] (can get, can set), p.array() ((documentedly) explicit copy of array). Point.coordinates() does not exist in Pybind11 DOLFIN.

[Prof Garth Wells]

@blechta In the discussion of the PR you reference there is no clear decision (not sure why it was merged in view of the discussion), and it was premised on SWIG. We have a chance now to try and do things 'right' with pybind11, and to make it simpler.

Point.array maps to Point::coordinates, and the C++ code is not nice in returning a plain pointer. If we were to change it to an Eigen array we could set and get values in a NumPy array, and deprecate the other methods.

[Jan Blechta]

I agree, but I note that this

Point.array maps to Point::coordinates

is an implementational detail. Python Point.array() is currently buggy; with the proposed fix (including docstring) it does semantically the same thing as SWIG DOLFIN and C++ DOLFIN. Feel free to change it further.

we could set and get values in a NumPy array

Fair enough, we would not need to handcode setitem and getitem. But let's do it quickly if there is a will. C++ and Python array() was added recently as a copy and now we say,"oh actually we don't like copy that much, let's make it a view."

[Prof Garth Wells]

I think we need to fix it. I objected to PR #328 because it was not clear if it was for setting, getting or both. This boils down of course to copies versus views. Let's not propagate or be dogmatic about fixing a recent design mistake.

[Jan Blechta]

Ok. Let's merge my fix anyway. It fixes the bug for the case that the design mistake is not addressed soon.

[Prof Garth Wells]

@blechta Can you register an issue so we don't forget to re-visit and fix the design?

[Jan Blechta]

Yes.

[Jan Blechta]

Fix in https://bitbucket.org/fenics-project/dolfin/branch/jan/fix-issue-954.

[Jan Blechta]

The interface design issue filed as #959.

[Jan Blechta]

Fix in next.

[Jan Blechta]

df8698e1ee2d33bb70e9d861244e137b95cb2030