GenericFunction::restrict() missing in Python
Two parts are missing to use restrict from Python:
- better pybind11 interface for
ufc::cell
, GenericFunction::restrict
bindings which would take lists or NumPy arrays; now pybind11 thinks incorrectly thatrestrict
takes floats (instead of float ptr/arrays) as coefficients and coordinates.
See a report by Sigvald Marholm https://www.allanswered.com/post/xzjpj/functionrestrict-changed-in-fenis-2018-1-0/
Partial fix for the first part:
diff --git a/python/src/function.cpp b/python/src/function.cpp
index 755f5b440..6aab98303 100644
--- a/python/src/function.cpp
+++ b/python/src/function.cpp
@@ -57,6 +57,7 @@ namespace dolfin_wrappers
// ufc::cell
py::class_<ufc::cell, std::shared_ptr<ufc::cell>>(m, "ufc_cell")
+ .def(py::init<>())
.def_readonly("cell_shape", &ufc::cell::cell_shape)
.def_readonly("topological_dimension", &ufc::cell::topological_dimension)
.def_readonly("geometric_dimension", &ufc::cell::geometric_dimension)
diff --git a/python/src/mesh.cpp b/python/src/mesh.cpp
index d8322a737..f549d6fc0 100644
--- a/python/src/mesh.cpp
+++ b/python/src/mesh.cpp
@@ -318,6 +318,8 @@ namespace dolfin_wrappers
std::vector<double> x;
self.get_vertex_coordinates(x);
return x; }, "Get cell vertex coordinates")
+ .def("get_cell_data", &dolfin::Cell::get_cell_data)
+ .def("get_cell_topology", &dolfin::Cell::get_cell_topology)
.def("orientation", (std::size_t (dolfin::Cell::*)() const) &dolfin::Cell::orientation)
.def("orientation", (std::size_t (dolfin::Cell::*)(const dolfin::Point&) const) &dolfin::Cell::orientation);
Comments (7)
-
-
reporter restrict
is in a tight assembly loop. We have to be sure to not spoil performance by the datatype change.- some C++ code might rely on the current signature.
By the above reasons it is preferable that correct mapping is done in the pybind11 wrapper (and any changes of
restrict
signature are to be done in DOLFINx if neeed). There might be a simple way how to achieve this in pybind11 - Eigen array wrapper to the pointers. Eigen arrays map natively to NumPy arrays in pybind11, so this could even give us backwards compatibility with SWIG/DOLFIN. -
I see. Well performance was about the only reason I could think of (if you're willing to go through all of dolfin to take care of the signature). Changes must happen in the pybind11 interface then.
-
The following (in python/src/function.cpp) works:
// dolfin::Function py::class_<dolfin::Function, std::shared_ptr<dolfin::Function>, dolfin::GenericFunction, dolfin::Hierarchical<dolfin::Function>> (m, "Function", "A finite element function") // ... .def("restrict", [](const dolfin::Function& self, Eigen::Ref<Eigen::VectorXd> w, const dolfin::FiniteElement& element, const dolfin::Cell& cell, const Eigen::Ref<Eigen::VectorXd> coordinate_dofs, const dolfin::Cell& _cell) { ufc::cell ufc_cell; _cell.get_cell_data(ufc_cell); self.restrict(w.data(), element, cell, coordinate_dofs.data(), ufc_cell); })
but I'm not quite sure about the ufc_cell argument. I expected it to work out-of-the-box using a ufc::cell argument with the fix you suggested but it didn't.
Also, Eigen vectors accepts numpy arrays but not lists.
cell.get_vertex_coordinates()
in python returns a list, not a numpy array so now it must be converted to numpy array (this wasn't necessary in 2017.2.0). -
reporter Cool. I would suggest few fixes.
-
This
restrict
wrapper should be done forGenericFunction
, it should work forExpression
too. -
Do we really need the second cell argument? It's just a heritage from old times when
ufc::cell
was exposed in DOLFIN/SWIG and the ufc_cell was of this datatype on Python side. I would suggest to remove it and just doufc::cell ufc_cell; cell.get_cell_data(ufc_cell);
in the wrapper. -
Regarding data types, you could implement more
restrict
wrappers to convert to NumPy if you really have a desire for that - it's trivial to implement the conversion in Python. But it is preferable to have just one version with nice pybind11 datatype checking at runtime. -
Could you add a test to
python/test/unit/function/test_{function,expression}.py
?
If you're not up to it, I will do it when I have time.
-
-
Now that I've started (on my first ever change in the FEniCS source code) it would be fun trying to complete it.
I kept the
ufc::cell
for consistency with the C++restrict()
method but I'll be glad to get rid of it. I'm wondering if it's possible to get rid of thecoordinate_dofs
argument as well? Or is there any reason why it might not always becell.get_vertex_coordinates()
? -
reporter Or is there any reason why it might not always be cell.get_vertex_coordinates()?
Good point. (In the internals it happens that if mesh is non-matching, the call is dispatched to
restrict_as_ufc_function
which just drops the cell (because mesh is non-matching) and uses the coordinates and UFC element to call interpolationufc::finite_element::evaluate_dofs()
which in turns callsGenericFunction::eval
.) But I agree, from Python high-level perspective, justrestrict()
using the provided cell. It will prevent ambiguity what it does, - Log in to comment
I just noticed that the declaration for restrict is:
Perhaps the plain old C-style dynamic arrays is the reason why pybind11 thinks it is floats? Is there any good reason to use this instead of vectors? How about
I don't know how many places this must be changed.