Coding a split/assign for Vector()

Question

I want to code a 'split' and 'assign' function that would work with Vector(). I was thinking about using the dofmap().dofs for each view of the MixedFunctionSpace. That is,

V1V2 = MixedFunctionSpace([V1, V2])

V1dofs = V1V2.sub(0).dofmap().dofs()
V2dofs = V1V2.sub(1).dofmap().dofs()

To split, I then build two matrices of dimension V1xV1V2, and V2xV1V2 with 1's and 0's given by the V1dofs and V2dofs. This seems to work fine if V1 and V2 are the same FunctionSpace. However, it fails if V1 and V2 are different. And most surprisingly, at least for CG spaces that differ only by their polynomial order, it seems that the split fails b/c of a single dof (the first one).

Looking at FunctionAssigner, it appears that Fenics builds these maps by iterating over all cells; I'd rather avoid that.

So my questions are:
1) how could I get the right map between the dofs of V1V2.sub(i) and Vi? Is there any way to extract that map from FunctionAssigner?
2) Am I correct to assume that what I'm doing will work if V1=V2 (say for CG spaces only)?
3) Would there be another way to split/assign for Vector()?

Answer 1

You cannot assume that V1V2.sub(0) and V1 (and V1V2.sub(1) and V2) have the same dofmaps. In parallel, they will usually not even have the same sizes.
In the most recent versions of dolfin the construction of mixed spaces have changed to avoid this confusion.

You can get maps from V1 and V2 into V1V2 using FunctionAssigner. For example,

v1, v2 = Function(V1), Function(V2)
w = Function(V1V2)
w.vector() = numpy.arange(*w.vector().local_range())

assigner = FunctionAssigner([V1,V2], V1V2)
assigner.assign([v1, v2], w)

map1 = v1.vector().array().astype(int)
map2 = v2.vector().array().astype(int)

Also note that collapsing a subspace preserves sizes but does not in general preserve dof ordering.

Comments

Thanks, Magne. But I'd rather have a solution that does not involve numpy.

My example doesn't seem to require FunctionAssigner to work, in parallel, if both FunctionSpaces are the same. Here is an example which uses the splitting technique I described in the original question (SplitAndAssign):

mesh = dl.UnitSquareMesh(40,40)
V1 = dl.FunctionSpace(mesh, "Lagrange", 1)
V2 = dl.FunctionSpace(mesh, "Lagrange", 1)
V1V2 = V1*V2
splitassign = SplitAndAssign(V1, V2, mesh.mpi_comm())

mpirank = dl.MPI.rank(mesh.mpi_comm())

u = dl.interpolate(dl.Expression(("x[0]*x[1]", "11+x[0]+x[1]")), V1V2)
uu = dl.Function(V1V2)
u1, u2 = u.split(deepcopy=True)
u1v, u2v = splitassign.split(u.vector())
u11 = dl.interpolate(dl.Expression("x[0]*x[1]"), V1)
u22 = dl.interpolate(dl.Expression("11+x[0]+x[1]"), V2)
a,b,c,d= dl.norm(u1.vector()-u1v), dl.norm(u2.vector()-u2v),\
dl.norm(u1.vector()-u11.vector()), dl.norm(u2.vector()-u22.vector())
if mpirank == 0:
    print 'Splitting an interpolated function:', a, b, c, d

Whether in serial or in parallel, I consistently get a=b=c=d=0.
Although as I was saying in the original question, the code breaks down if V2 is quadratic elements.