How to pass numpy array to compiled expression member variable?

Question

Trying to pass a numpy array to a compiled expression member variable, 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 init_v(
        const Array<double>& vv)
    {
        v = &vv;
        std::cout << "v = " << v->str(1) << std::endl;
    }
    void eval(Array<double> &expr, const Array<double> &pos) const
    {
        std::cout << "v = " << v->str(1) << std::endl;
    }
};
}'''
my_expr = dolfin.Expression(
    cppcode=cpp_code,
    element=fe)
my_expr.init_v(numpy.zeros(2))
my_expr.eval(numpy.zeros(1), numpy.zeros(2))

This compiles and runs fine but of course, since the member pointer v is const, even after assigning a vector 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 still empty when executing the eval function. Now, if I make the member pointer not const:

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

This won't compile with the error:

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

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 init_v(
        Array<double>& vv)
    {
        v = &vv;
        std::cout << "v = " << v->str(1) << std::endl;
    }
    void eval(Array<double> &expr, const Array<double> &pos) const
    {
        std::cout << "v = " << v->str(1) << std::endl;
    }
};
}'''
my_expr = dolfin.Expression(
    cppcode=cpp_code,
    element=fe)
my_expr.init_v(numpy.zeros(2))
my_expr.eval(numpy.zeros(1), numpy.zeros(2))

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

my_expr.init_v(numpy.zeros(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

Comments

Any hint, anyone? Thanks! Martin

---

I found a proposed solution in this thread:

https://answers.launchpad.net/dolfin/+question/227358

(see Response #4)

however, this code does not appear to run in the current version of fenics (?)

When I run this code, I get the error:

... PYTHON_wrap.cxx: In function ‘PyObject* 
        _wrap_K__coeffs_set(PyObject*, PyObject*)’:
... PYTHON_wrap.cxx:6423:58: error: 
    ‘PyArrayObject {aka struct tagPyArrayObject}’ 
     has no member named ‘data’
       double *data = reinterpret_cast<double*>((*xa).data); 
        for (int i=0;i<m;++i) arg2[i] = &data[i*n]; }  else {
                                                      ^

Answer 1

Does anybody on this forum know how to do this using Instant?

Answer 2

Have you tried to cast away the const? The compiler will react on const, but it may be cast away. It is of course a hack of a solution.

Answer 3

I found a solution that involves implementing a linked list in C++ and then adding the array elements one at a time to the list. This doesn't seem very elegant, especially with all the work that people have done on more sophisticated structures using Instant. Unfortunately, I'm not really a C++ programmer and I couldn't figure out how to implement the stuff in Ch 14 of the Fenics Book. Oh well, this works for my purposes and it's pretty fast too.

First, I have this C++ file:

namespace dolfin {

struct node {
    double x;
    double v;
    node *next;
};


class MyFun : public Expression
{
    private:
        node *root;

    public:
        MyFun(): Expression()
        {
            root = new node;
            root->next = 0;
        };
        void eval(Array<double>& values, const Array<double>& x) const
        {
            double vv = findval(x[0]);
            values[0] = vv;
        };

        void update(double _x, double _v)
        {
            node *newNode = new node;
            newNode->x = _x;
            newNode->v = _v;
            newNode->next = root;
            root = newNode;
        };

        double findval(double _x) const
        {
            node *p = root->next;
            node *q = root;
            while (true)
            {
            // Assume that root node has biggest x-value.
            // Traverse down the list until p.x is too small.
                if (p->x <= _x) break;
                else
                {
                    q = p;
                    p = p->next;
                }
            }
            // Linear interpolation
            double interp = p->v + (_x - p->x) * (p->v - q->v) / (p->x - q->x);
            return interp;

        };

};
};

As you can see, in my particular application I need to interpolate over elements of the Numpy Array. This could could obviously be modified to do other things.

Then on the Python side I call C++ as:

header_file = open("ll.cpp")
code = "\n".join(header_file.readlines())
ub = Expression(cppcode=code, degree=2)
for j in range(0,5):  ub.update(vx[j],vv[j])

This then allows me to interpolate the values (vx,vy) over the boundary as a compiled expression:

bc_ux = DirichletBC(W.sub(0).sub(0), ub, bed);