Extracting Vertex DoFs of a CG2 (Vector) Function in Parallel

Question

Let's say I have a function with data in a CG2 Vector Function Space and would like to extract the values of the solution corresponding to the vertex dofs.

For simplicity, let's just define an Expression and interpolate the expression on both spaces:

from dolfin import *
mesh = UnitSquareMesh(32,32) 
V = VectorFunctionSpace( mesh, "CG", 2 )
V_CG1 = VectorFunctionSpace( mesh, "CG", 1 )
vel_exp = Expression( ('1.0', '-1.0'), degree=1 )
vel = interpolate( vel_exp , V)
vel2 = interpolate( vel_exp , V_CG1)

Ideally, I'd like to be able to extract the vertex dofs from vel corresponding to the same values stored in vel2.

I was able to do this in serial by creating a map of the CG1 dof indices to the vertex CG2 indices:

import numpy as np

dofs = V.dofmap().dofs(mesh, 0)
dofs_CG1 = V_CG1.dofmap().dofs(mesh, 0)
n_dofs
Vdofs_vert = np.zeros( n_dofs, dtype = 'int' )
for i in xrange( n_dofs ):
        j = dofs_CG1[i]
        Vdofs_vert[j] = dofs[i]

And then vel.vector().array()[ Vdofs_vert ] would give the desired result. However this doesn't work in parallel since the global dof indices are returned with the dof maps.

Is there a simple way to do this in parallel? I would like to avoid using interpolate from CG2 to CG1 since it's slow and in my view, unnecessary.

Thanks!

Comments

EDIT: Below the solution is incorrect, still not entirely sure what dofs(mesh, 0) returns in parallel.

I THINK I found a solution. If anyone can comment, please let me know!

This can be done using get_local() but taking the array of vertex dofs in as an argument:

Here's a small example:

import dolfin as dfn
import numpy as np

# define mesh and spaces
nx = 8
mesh = dfn.UnitSquareMesh(nx,nx)
V_CG1 = dfn.VectorFunctionSpace( mesh, "CG", 1 )
V = dfn.VectorFunctionSpace( mesh, "CG", 2 )

# setup dof maps to list vertex dofs
dofs = V.dofmap().dofs(mesh, 0)
dofs_CG1 = V_CG1.dofmap().dofs(mesh, 0)

# define some vector expression and interpolate it on each space
vel_exp = dfn.Expression( ('x[0]', 'x[1]'), degree=1 )
vel = dfn.interpolate( vel_exp , V)
vel2 = dfn.interpolate( vel_exp , V_CG1)

# create numpy arrays of appropriate size to store the dof values
test_vel = np.zeros( dofs.size )
test_CG1 = np.zeros( dofs_CG1.size )

# get the local values (this includes shared dofs)
vel.vector().get_local(test_vel,dofs)
vel2.vector().get_local(test_CG1,dofs_CG1)

# take the difference of the result to check if we get zero
print test_CG1-test_vel

And I get zero for various expressions and mesh sizes.