How to pass numpy array to compiled expression member variable?

Coming back to an old problem I never got to solve properly… I need to pass a numpy array to a compiled expression, to no avail. Consider the following code:

import numpy
import dolfin
mesh = dolfin.UnitSquareMesh(1,1)
fe = dolfin.FiniteElement(
    family="Quadrature",
    cell=mesh.ufl_cell(),
    degree=1,
    quad_scheme="default")
cpp_code = '''\
namespace dolfin
{
class MyExpr : public Expression
{
    const Array<double>* v;
public:
    MyExpr(): Expression() {}
    void print_v() const
    {
        std::cout << "v = " << v->str(1) << std::endl;
    }
    void init_v(
        const Array<double>& v_)
    {
        print_v();
        v = &v_;
        print_v();
    }
    void eval(
              Array<double> &expr,
        const Array<double> &pos) const
    {
        print_v();
    }
};
}'''
my_expr = dolfin.Expression(
    cppcode=cpp_code,
    element=fe)
my_expr.init_v(numpy.ones(2))
my_expr.print_v()
my_expr.eval(numpy.ones(1), numpy.ones(3))

This compiles and runs fine but of course, since the member pointer v is const, even after assigning an array to it in the init_v function (actually I’m not even sure how/why I can do that since it is const, but anyway), it is empty when executing the print_v function, and equal to first parameter of the function when executing the eval function (so weird, right?). Now, if I make the member pointer not const:

cpp_code = '''\
namespace dolfin
{
class MyExpr : public Expression
{
    Array<double>* v;
public:
    MyExpr(): Expression() {}
    void print_v() const
    {
        std::cout << "v = " << v->str(1) << std::endl;
    }
    void init_v(
        const Array<double>& v_)
    {
        print_v();
        v = &v_;
        print_v();
    }
    void eval(
              Array<double> &expr,
        const Array<double> &pos) const
    {
        print_v();
    }
};
}'''
my_expr = dolfin.Expression(
    cppcode=cpp_code,
    element=fe)
my_expr.init_v(numpy.ones(2))
my_expr.print_v()
my_expr.eval(numpy.ones(1), numpy.ones(3))

This won’t compile with the error:

error: invalid conversion from ‘const dolfin::Array<double>*’ to ‘dolfin::Array<double>*’ [-fpermissive]
     v = &v_;
       ^

This makes sense. Now, if I remove the other const as well:

cpp_code = '''\
namespace dolfin
{
class MyExpr : public Expression
{
    Array<double>* v;
public:
    MyExpr(): Expression() {}
    void print_v() const
    {
        std::cout << "v = " << v->str(1) << std::endl;
    }
    void init_v(
        Array<double>& v_)
    {
        print_v();
        v = &v_;
        print_v();
    }
    void eval(
              Array<double> &expr,
        const Array<double> &pos) const
    {
        print_v();
    }
};
}'''
my_expr = dolfin.Expression(
    cppcode=cpp_code,
    element=fe)
my_expr.init_v(numpy.ones(2))
my_expr.print_v()
my_expr.eval(numpy.ones(1), numpy.ones(3))

This compiles fine, but at runtime I’m getting the following error:

my_expr.init_v(numpy.ones(2))
TypeError: in method 'MyExpr_init_v', argument 2 of type 'dolfin::Array< double > &'

So apparently, the compiled module recognizes a numpy array as const dolfin::Array< double > & , but not as dolfin::Array< double > &. Any idea on how to make this work? Thanks! Martin

Hey Martin,

Here’s a tester code I tried to pass a numpy array to a dolfin expression in C++ . It should help :

import dolfin as d
import numpy as np

code = """
    
    #include <pybind11/pybind11.h>
    #include <pybind11/eigen.h>
    #include <dolfin/function/Expression.h>
    
    class test_exp : public dolfin::Expression {
      public:
        
        Eigen::VectorXd arr;
        
        test_exp() : dolfin::Expression() { }

        void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x) const {
            
            values[0] = arr.sum();
        }

    };
    
    PYBIND11_MODULE(SIGNATURE, m) {
        pybind11::class_<test_exp, std::shared_ptr<test_exp>, dolfin::Expression>
        (m, "test_exp")
        .def(pybind11::init<>())
        .def_readwrite("arr", &test_exp::arr)
        ;
    }
    """

exp = d.CompiledExpression(d.compile_cpp_code(code).test_exp(), degree=1)
exp.arr = np.array([1., 2., 4., 8.], dtype=float)
print(exp([0.,0.,0.]))   # Prints 15

Note : Initialize the exp.arr object in python to avoid seg_faults.
Ref: test_jit.py
(Tested on 2019.1.0)

3 Likes

OMG thank you!! (test_jit.py is such a useful source of solution!) Just wondering: is there a similar procedure with dolfin::Array object?

Hey Martin,

I haven’t personally tried using the dolfin::Array class as an interface between Python and C++, but I think it should technically be possible. I can maybe suggest you this -

Interfacing NumPy objects in pybind is possible by exposing pointers to the data buffer. Maybe you can pass along the buffer as reference and construct a dolfin::Array object within C++ code and use it? (may require verifying data type compatibility and conversion)

If you do use the interface between Python and C++, I think using Eigen classes with Numpy is the easiest as in the earlier answer. You can read a little more on it here. Keep in mind that there is a proposal in FEniCS to adapt to this (Ref), although it hasn’t yet been sorted.

1 Like

All clear, thank you so much!

1 Like