Function assignment failing with Mixed elements with real element

[Mike Welland]

I have a mixed element space with two CG element spaces and a real element space (R). I'm trying to initialize the first space (CG) with an expression that depends on x[0]. Example below will plot the function which is desired and then the result of the assignment.

It works in 1-D (Use the interval mesh), and if you remove the R function space from the mixed elements, but it fails with 2-D. Any ideas? Workarounds? I tried creating two other functions and initializing the whole mixed function with assign(V,[x0, x1, xR]) but got the exact same.

from dolfin import *
# create a mixed function space X = X0 x X1
mesh = RectangleMesh(0,0,1,1,10, 10)
#mesh = IntervalMesh(10,0,1)

V1 = FunctionSpace(mesh, 'CG', 1)
R = FunctionSpace(mesh,'R',0)
V = MixedFunctionSpace([V1, V1,R])
f = Function(V)

# create any nonzero functions in X0 and X1
x0 = interpolate(Expression("x[0]"),V1)
assign( f.sub(0), x0)

plot(x0)
plot(f.split()[0])
interactive()

Thanks!

[Simon Funke]

This is a duplicate of issue #193

[Jan Blechta]

Duplicate of #160.

[Mike Welland]

Here is a workaround from: http://fenicsproject.org/qa/2740/assign-giving-wrong-results-assigning-subfunctions-function

from dolfin import *
# create a mixed function space X = X0 x X1
mesh = RectangleMesh(0,0,1,1,10, 10)
#mesh = IntervalMesh(10,0,1)

V1 = FunctionSpace(mesh, 'CG', 1)
R = FunctionSpace(mesh,'R',0)
V = MixedFunctionSpace([V1, V1,R])
f = Function(V)

# create any nonzero functions in X0 and X1
x0 = interpolate(Expression("x[0]"),V1)
assign( f.sub(0), x0)

dofs = V.sub(0).dofmap().collapse(mesh)[1].values()
f.vector()[dofs] = x0.vector()

plot(x0)
plot(f.split()[0])
interactive()

[Johan Hake]

This workaround is nice and could be used in the FunctionAssigner, however it does not work in parallel. That is because the map returned by the collapsed dof map include ghost dofs. It is easy to sort these out by:

min_dof, max_dof = V.dofmap().ownership_range()
dofs = np.array(V.sub(0).dofmap().collapse(mesh)[1].values())
dofs = dofs[(min_dof <= dofs) * (dofs < max_dof)]

However this reveals another short coming pointed out previously as dofs are randomly distributed to different processors for each mesh entity (in this case vertices).