GenericFunction::restrict() missing in Python

[Jan Blechta]

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 that restrict 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);

[Sigvald]

I just noticed that the declaration for restrict is:

    virtual void restrict(double* w,
                          const FiniteElement& element,
                          const Cell& dolfin_cell,
                          const double* coordinate_dofs,
                          const ufc::cell& ufc_cell) const override;

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

    virtual void restrict(std::vector<double>& w,
                          const FiniteElement& element,
                          const Cell& dolfin_cell,
                          const std::vector<double>& coordinate_dofs,
                          const ufc::cell& ufc_cell) const override;

I don't know how many places this must be changed.

[Jan Blechta]

1. restrict is in a tight assembly loop. We have to be sure to not spoil performance by the datatype change.
2. 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.

[Sigvald]

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.

[Sigvald]

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).

[Jan Blechta]

Cool. I would suggest few fixes.

1. This restrict wrapper should be done for GenericFunction, it should work for Expression too.

2. 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 do ufc::cell ufc_cell; cell.get_cell_data(ufc_cell); in the wrapper.

3. 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.

4. 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.

[Sigvald]

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 the coordinate_dofs argument as well? Or is there any reason why it might not always be cell.get_vertex_coordinates()?

[Jan Blechta]

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 interpolation ufc::finite_element::evaluate_dofs() which in turns calls GenericFunction::eval.) But I agree, from Python high-level perspective, just restrict() using the provided cell. It will prevent ambiguity what it does,