Trouble Interpolating piecewise Constant onto function space

Question

I want to project a user defined piecewise constant onto a function space

from dolfin import *

V = FunctionSpace ( mesh, "CG", 1 )

class MyExpression(Expression):
def eval(self, value, x):
if x[0] <= .5:
value[0] = 1.0
else:
value[0] = 0.0

def value_shape(self):
  return (1,)

u_init = MyExpression()
u = interpolate ( u_init, V )

When I run this code I get the following error message

*** Error: Unable to interpolate function into function space.
*** Reason: Rank of function (1) does not match rank of function space (0).
*** Where: This error was encountered inside FunctionSpace.cpp.
*** Process: unknown

Where as if I use an expression such as f = Expression("x[0]") everything runs fine. I don't understand why the two are fundamentally different.

Answer 1

The expression has rank 1 because of

def value_shape(self):
    return (1,)

Rank 1 expressions are vector-valued, even though it's just a length 1 vector in this case. Scalar expression should have rank 0. Since the function space has rank 0, you get an error.

If you replace the code with

def value_shape(self):
    return ()

then the expression gets rank 0 and your code will work. This is also the default for Expression subclasses.