efficient function composition

Question

I would like to have a composition of two functions f(g(x)) in a form which can be achieved by a custom Expression in python. In order to speed up assembly, I tried to implement this in C++ and pass it to Expression. Unfortunately, this fails for reasons unclear to me. Even if I comment out all everything in (B), the code still throws an exception.
Another issue is that the example won't compile with (A) although the header for Array is included.

Thanks for advise, Martin

mesh = UnitSquareMesh(5,5)
V = FunctionSpace(mesh,'CG',1)
V2 = VectorFunctionSpace(mesh,'CG',1)
g = interpolate(Expression(('x[0]+x[1]','-x[0]')),V2)
f = interpolate(Expression('1.+x[0]*x[0]'),V)

code = '''
#include "dolfin/common/Array.h"
namespace dolfin{
class MyFunc : public Expression
{
public:
// (A)
//  Array<double> val(2);
  boost::shared_ptr<Function> f;
  boost::shared_ptr<Function> g;

  MyFunc() : Expression() { }

  void eval(Array<double>& values, const Array<double>& x,
          const ufc::cell& c) const
  {
// (B)
      static Array<double> val(2);
      g->eval(val, x, c);
      f->eval(values, val, c);
  }
};}'''

fg = Expression(code)
fg.f = f
fg.g = g
fg(0.1,0.2)

Answer 1

You have implemented the version of Expression::eval that needs a cell. When you evaluate the Expression at an arbitrary point it does not know what cell it is in. The following code fixes that and the compilation of the Expression

from dolfin import *
import numpy as np

mesh = UnitSquareMesh(5,5)
V = FunctionSpace(mesh,'CG',1)
V2 = VectorFunctionSpace(mesh,'CG',1)
g = interpolate(Expression(('x[0]+x[1]','-x[0]')),V2)
f = interpolate(Expression('1.+x[0]*x[0]'),V)

code = '''
class MyFunc : public Expression
{
public:
  boost::shared_ptr<Function> f;
  boost::shared_ptr<Function> g;

  MyFunc() : Expression() { }

  void eval(Array<double>& values, const Array<double>& x,
            const ufc::cell& c) const
  {
      Array<double> val(2);
      g->eval(val, x, c);
      f->eval(values, val, c);
  }
};
'''

fg = Expression(code)
fg.f = f
fg.g = g

# Find cell that intersects Point(0.1, 0.2)
p = Point(0.1, 0.2)
bt = mesh.bounding_box_tree()
cell_id = bt.compute_first_entity_collision(p)
cell = Cell(mesh, cell_id)
values = np.zeros(1)
x = np.array([p.x(), p.y()])
fg.eval_cell(values, x, cell)
print values

Comments

Ok, I see. Is bounding_box_tree a new method? It does not seem to be present in the nightly snapshot release and I cant find documentation. Btw, since the nightly builds seem to have failed for quite some time, see
https://launchpad.net/~fenics-packages/+archive/fenics-dev/+packages
the available dolfin packages from that PPA are from 20130618 instead of 20131004.

---

Yes bounding_box_tree is a new method but in DOLFIN 1.2 and earlier it is even easier:

cell_id = mesh.intersected_cell(p)