evaluate_dofs for expressions is sub-optimal for vector, tensor or mixed function spaces

Interpolating an Expression is currently sub-optimal for two reasons:

it iterates over the cells and the over the dofs (by restricting it), but for Lagrange-type basis function it would be more convenient to iterate over the dofs;
for vector, tensor, mixed function spaces, it recompute the same function (same cell and coordinates) for each component, when the result is clearly the same.

The first point is DOLFIN-related, but the second is due to FFC.

Here an example:

from fenics import *

code = \
"""\
class MyExpression : public Expression
{
public:
    MyExpression() : Expression(2) {};

    void eval(dolfin::Array<double>&,
              const dolfin::Array<double>&) const
    {
        std::cout << "call" << std::endl;
    }
};
"""

mesh = UnitTriangleMesh()
V = VectorFunctionSpace(mesh, "P", 1)
e = Expression(cppcode = code, element = V.ufl_element())
interpolate(e, V)

From the output you can easily count 6 "call"s (1 cell x 6 dofs). The generated code for the vector space overloads the method ufc::evaluate_dofs as follows:

  /// Evaluate linear functionals for all dofs on the function f
  virtual void evaluate_dofs(double* values,
                             const ufc::function& f,
                             const double* vertex_coordinates,
                             int cell_orientation,
                             const ufc::cell& c) const
  {
    // Declare variables for result of evaluation
    double vals[2];

    // Declare variable for physical coordinates
    double y[2];
    y[0] = vertex_coordinates[0];
    y[1] = vertex_coordinates[1];
    f.evaluate(vals, y, c);
    values[0] = vals[0];
    y[0] = vertex_coordinates[2];
    y[1] = vertex_coordinates[3];
    f.evaluate(vals, y, c);
    values[1] = vals[0];
    y[0] = vertex_coordinates[4];
    y[1] = vertex_coordinates[5];
    f.evaluate(vals, y, c);
    values[2] = vals[0];
    y[0] = vertex_coordinates[0];
    y[1] = vertex_coordinates[1];
    f.evaluate(vals, y, c);
    values[3] = vals[1];
    y[0] = vertex_coordinates[2];
    y[1] = vertex_coordinates[3];
    f.evaluate(vals, y, c);
    values[4] = vals[1];
    y[0] = vertex_coordinates[4];
    y[1] = vertex_coordinates[5];
    f.evaluate(vals, y, c);
    values[5] = vals[1];
  }

can be rewritten as

  /// Evaluate linear functionals for all dofs on the function f
  virtual void evaluate_dofs(double* values,
                             const ufc::function& f,
                             const double* vertex_coordinates,
                             int cell_orientation,
                             const ufc::cell& c) const
  {
    // Declare variables for result of evaluation
    double vals[2];

    // Declare variable for physical coordinates
    double y[2];
    y[0] = vertex_coordinates[0];
    y[1] = vertex_coordinates[1];
    f.evaluate(vals, y, c);
    values[0] = vals[0];
    values[3] = vals[1];
    y[0] = vertex_coordinates[2];
    y[1] = vertex_coordinates[3];
    f.evaluate(vals, y, c);
    values[1] = vals[0];
    values[4] = vals[1];
    y[0] = vertex_coordinates[4];
    y[1] = vertex_coordinates[5];
    f.evaluate(vals, y, c);
    values[2] = vals[0];
    values[5] = vals[1];
  }

Here we save half of the cost of calling eval, but for tensor-valued functions in 3d the cost would be 1/9th.

Comments (22)