evaluate components all-at-once/one-by-one for expressions defined over vector elements

Question

Hi,

Consider the following code:

mesh = dolfin.UnitSquareMesh(1,1)
dX = dolfin.dx(mesh)
ve = dolfin.VectorElement(
    family="Quadrature",
    cell=mesh.ufl_cell(),
    degree=1,
    quad_scheme="default")

class MyExpr2(dolfin.Expression):
    def __init__(self, **kwargs):
        pass

    def value_shape(self):
        return (2,)

    def eval(self, Expr, X):
        print "    X = " + str(X)
        Expr[0] = 0
        Expr[1] = 1
        print "    Expr = " + str(Expr)

myExpr = MyExpr(element=ve)
dolfin.assemble(dolfin.dot(myExpr,dolfin.Constant((1,1)))*dX, form_compiler_parameters={'quadrature_degree': 1})

Which returns

X = [ 0.66666667  0.33333333]
Expr = [ 0.  1.]
X = [ 0.66666667  0.33333333]
Expr = [ 0.  1.]
X = [ 0.33333333  0.66666667]
Expr = [ 0.  1.]
X = [ 0.33333333  0.66666667]
Expr = [ 0.  1.]

Clearly, the expression is evaluated twice at each integration point. I assume that components are evaluated one by one. Is there a way to evaluate them all at once? Thanks!

Answer 1

The output is produced when the expression is interpolated into a FunctionSpace.

For the interpolation, the expression is evaluated at every dof for each component separately. In your code this leads to excessive evaluations, but in general, this is necessary because the dofs for each component may differ.

To avoid the excessive evaluations you need a more specialized interpolation function. Ideally, you would loop over quadrature points instead of looping over elements and dofs. The LagrangeInterpolator class appears to do this. For example, consider

Q = FunctionSpace(mesh, ve)
f = Function(Q)
LagrangeInterpolator().interpolate(f, myExpr)
dolfin.assemble(dolfin.dot(f,dolfin.Constant((1,1)))*dX, form_compiler_parameters={'quadrature_degree': 1})

Comments

Good solution but eventually requires additional memory for storing dofmap of Q, dofs of f. In the end construction of Q and f may be much higher penalty.

Answer 2

This is know problem in FFC https://bitbucket.org/fenics-project/ffc/issues/52.