value shape in subclass vs cpp expression

Question

Trying to turn some expression from python subclass into JIT cpp.

First the basics:

import dolfin

mesh = dolfin.UnitSquareMesh(1,1)
dX = dolfin.dx(mesh)

fe = dolfin.FiniteElement(
    family="Quadrature",
    cell=mesh.ufl_cell(),
    degree=1,
    quad_scheme="default")

Now the python sublass:

class PyExpr(dolfin.Expression):
    def eval(self, values, position):
        print "position = "+str(position)
        values[0] = 1.
        print "values = "+str(values)
pyExpr = PyExpr(element=fe)
dolfin.assemble(pyExpr * dX)

This works fine. Now if I try to use a JIT expression without specifying the dimension

cppExprCode='''
namespace dolfin
{

class CppExpr : public Expression
{
public:

    CppExpr(): Expression()
    {
    }

    void eval(Array<double>& values, const Array<double>& position) const
    {
        std::cout << "position = " << position << std::endl;
        values[0] = 1.;
        std::cout << "values = " << values << std::endl;
    }
};

}'''
cppExpr = dolfin.Expression(cppExprCode, element=fe)
dolfin.assemble(cppExpr * dX)

It works ok but the dimension of the expression is set to three. Now if I try to set it to one:

cppExprCode='''
namespace dolfin
{

class CppExpr : public Expression
{
public:

    CppExpr(): Expression(1)
    {
    }

    void eval(Array<double>& values, const Array<double>& position) const
    {
        std::cout << "position = " << position << std::endl;
        values[0] = 1.;
        std::cout << "values = " << values << std::endl;
    }
};

}'''
cppExpr = dolfin.Expression(cppExprCode, element=fe)
dolfin.assemble(cppExpr * dX)

I'm getting the error:

ValueError: The value shape of the passed 'element', is not equal to the value shape of the compiled expression.

What would be the correct way to do so? Thanks!

Comments

Hi, what makes you say that the dimension of the expression is set to three in the first compiled example?

---

Thanks for your help! If you run the code, you get

position = [ 0.66666667  0.33333333]
values = [ 1.]
position = [ 0.33333333  0.66666667]
values = [ 1.]

for the python subclass expression, but you get

position = <Point x = 0.666667 y = 0.333333 z = 0>
values = <Point x = 1 y = 0 z = 0>
position = <Point x = 0.333333 y = 0.666667 z = 0>
values = <Point x = 1 y = 0 z = 0>

for the cpp jit expression.

---

Actually, writing

std::cout << "values.size() = " << values.size() << std::endl;
std::cout << "values = " << values << std::endl;

gives

values.size() = 1
values = <Point x = 1 y = 0 z = 0>
values.size() = 1
values = <Point x = 1 y = 0 z = 0>

which I find a little strange. Anyway, the dimension is indeed 1. Thanks!