UnitCubeMesh creates duplicate vertices
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?
Comments (29)
-
-
reporter 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?
-
reporter 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 thepython/
directory. Some extra packages likeh5py
are also installed in the Docker recipe, but I assume this is not relevantCan anyone reproduce the bug?
-
reporter -
reporter CircleCI now replicates the bug
-
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 temporaryPoint
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.)
-
What did we agree before? That (in Python)
Point.array()
is a view or a copy? In C++ it is a copy contrary toPoint::coordinates()
. -
reporter 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
-
Can we just avoid exposing
Point.coordinates()
? -
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)
-
reporter Maybe, but then indexing point should be better documented
Edit: we should just expose
.coordinates()
directly (returns a Python list) and then convertstd::array
toEigen::Vector3d
later indolfin::Point
which will make it return a numpy array Edit2: sorry, I see that coordinates returnsdouble*
-
The essence of the bug is
let_us_take_reference_to_plain_ptr = Eigen::Map<Eigen::Vector3d>(self.coordinates());
-
So the problem is not using
self.coordinates()
but using it wrongly. -
The point is do we need more than one way to do the same thing?
-
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). -
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 isPoint.coordinates()
in C++ andPoint.__(gs)etitem__
in Python allowing some shallow manipulations. -
@blechta We do have a whole bunch of ways to access and set coordinates of a
Point
. Why do we need so many? -
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 anyarray()
and only havingcoordinates()
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/activityIf 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.
-
Btw, we only have one way of setting a coordinate, which is
Point.__setitem__
. Getting is possible throughp.x()
(stupid),p[0]
(can get, can set),p.array()
((documentedly) explicit copy of array).Point.coordinates()
does not exist in Pybind11 DOLFIN. -
@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 toPoint::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. -
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." -
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.
-
Ok. Let's merge my fix anyway. It fixes the bug for the case that the design mistake is not addressed soon.
-
@blechta Can you register an issue so we don't forget to re-visit and fix the design?
-
Yes.
-
- changed milestone to 2017.2
-
assigned issue to
-
The interface design issue filed as #959.
-
Fix in next.
-
- changed status to resolved
- Log in to comment
That's not ok. But I can't reproduce it. What is the reproducible way to get this behaviour?