Help: calling C++ Tensor Expressions from Python (dolfin 1.5)

Hi all,

I want to assemble

$$ \left( \Theta \nabla u, \nabla v \right) $$

where $\Theta$ is a s.p.d. tensor defined in a c++ class that inherits from Expression.

I have prepared the following code for the 2D case, but I get some run-time errors I can not figure out.

cpp_expression.h

#include <dolfin/function/Expression.h>

namespace dolfin
{

class AnisTensor2D : public Expression
{
public:

  AnisTensor2D();
  void set(const double & t0, const double & t1, const double & a);    
  void eval(Array<double>& values, const Array<double>& x) const;

private:

  int ndim;
  double theta0;
  double theta1;
  double alpha;

  double c00;
  double c01;
  double c11;

};

cpp_expression.cpp

AnisTensor2D::AnisTensor2D() :
Expression(2,2),
ndim(2),
theta0(1.),
theta1(1.),
alpha(0),
c00(1.),
c01(0.),
c11(1.)
{

}

void AnisTensor2D::set(const double & t0, const double & t1, const double & a)
{
    theta0 = t0;
    theta1 = t1;
    alpha = a;

  double sa = sin(alpha);
  double ca = cos(alpha);
  c00 = theta0*sa*sa + theta1*ca*ca;
  c01 = (theta0 - theta1)*sa*ca;
  c11 = theta0*ca*ca + theta1*sa*sa;

}

double AnisTensor2D::Theta0() const {return theta0;}
double AnisTensor2D::Theta1() const {return theta1;}
double AnisTensor2D::Alpha() const {return alpha;}
int AnisTensor2D::NDim() const {return ndim;}

void AnisTensor2D::eval(Array<double>& values, const Array<double>& x) const
{
   values[0] = c00;
   values[1] = c01;
   values[2] = c01;
   values[3] = c11;
}

Python script

import dolfin as dl
import os
import math

if __name__ == "__main__":
    dl.set_log_active(False)
    ndim = 2
    nx = 64
    ny = 64
    mesh = dl.UnitSquareMesh(nx, ny)
    Vh = dl.FunctionSpace(mesh, "CG", 1)

    sdir = "<path>"
    header_file = open(os.path.join(sdir,"cpp_expression.h"), "r")
    code = header_file.read()
    header_file.close()
    cpp_sources = ["cpp_expression.cpp"] 
    hypExpression = dl.compile_extension_module(code=code, source_directory=sdir, sources=cpp_sources, include_dirs=[".",  sdir])

    T = hypExpression.AnisTensor2D()
    T.set(2., .5, math.pi/4)

    u = dl.TrialFunction(Vh)
    v = dl.TestFunction(Vh)
    a = dl.inner(T*dl.grad(u), dl.grad(v))*dl.dx

    A = dl.assemble(a)

Error:

Traceback (most recent call last):
  File "<path>/AnisTensor2.py", line 25, in <module>
    a = dl.inner(T*dl.grad(u), dl.grad(v))*dl.dx
TypeError: unsupported operand type(s) for *: 'AnisTensor2D' and 'Grad'

How can I modify the class AnisTensor2D, so that the expression T*dl.grad(u) is supported?

NOTE: Please do not suggest to rewrite the class in python, because I need to have it accessible from C++ in a different part of the code.

Thank you in advance for your help.

Best,

Umbe