Interpolation on mixed functionspace - shapes not matching

Question

I have encountered a weird problem with shapes of the mixed function space not matching the "original" function space. The problem only arises, when i apply periodic boundary conditions. I have setup a small MWE where

1) Meshes and spaces are created.
2) A dummy problem is solved.
3) A constant is added to the obtained solution.

It is in step (3) that my problem arises. When i try to add a constant (interpolated on the relevant sub space) i get the following error

ValueError: operands could not be broadcast together with shapes (441,) (420,)

From the error message it is clear that length of the solution vector (u1) is different from the length of the interpolated expression. However, i can't figure out why this happens. Can anyone enlighten me and/or suggest a fix to the problem?

from fenics import *

class PeriodicDomain(SubDomain):
    def __init__(self, ly):
        self.ly = ly
        SubDomain.__init__(self)

    def inside(self, x, on_boundary):
        # return True if on bottom boundary
        return near(x[1], 0) and on_boundary

    def map(self, x, y):
        # return Map top boundary to bottom boundary
        y[0] = x[0]
        y[1] = x[1] - self.ly

# Create mesh and function spaces.
mesh = UnitSquareMesh(10, 10)
el = FiniteElement("CG", triangle, 2)
V = FunctionSpace(mesh, el, constrained_domain=PeriodicDomain(1.0))
W = FunctionSpace(mesh, MixedElement([el, el]), constrained_domain=PeriodicDomain(1.0))
# Solve dummy problem.
u1, u2 = split(TrialFunction(W))
v1, v2 = split(TestFunction(W))
a = inner(grad(u1), grad(v1))*dx() + inner(grad(u2), grad(v2))*dx()
L = v1*dx()
u = Function(W)
solve(a == L, u)
# Try to add something to the solution.
u1, u2 = u.split(True)
u1.vector()[:] = u1.vector().array() + interpolate(Expression('1.0', el), V).vector().array()

Answer 1

On the master branch of dolfin the u.split(True) doesn't work, but this works, though:

# Split the function
fa = FunctionAssigner([V, V], W)
u1, u2 = Function(V), Function(V)
fa.assign([u1, u2], u)

# Test sizes
u3 = interpolate(Expression('1.0', el), V)
print u1.vector().array().shape
print u2.vector().array().shape
print u3.vector().array().shape

Comments

Sorry, i forgot to note that i am on the 1.7.0dev branch. Here

u1, u2 = u.split(True)

has worked for me so far, but in this special case it apparently causes problems. Replacing the line with your suggestion, e.g.

fa = FunctionAssigner([V, V], W)
u1, u2 = Function(V), Function(V)
fa.assign([u1, u2], u)

fixed the problem.