How to interpret output of vertex_to_dof_map and dof_to_vertex_map for vector function spaces?

Question

I am confused by the way vertex_to_dof_map and dof_to_vertex_map work for vector function spaces.

Consider the short snippet below.

from dolfin import *
mesh = UnitSquareMesh(5,5)
Q = VectorFunctionSpace(mesh, "CG", 1, dim=2)
q = Function(Q)
v2d=vertex_to_dof_map(Q)
d2v=dof_to_vertex_map(Q)

This is a 2D vector function so I expect to have 2 degrees of freedom for each vertex. When I run dof_to_vertex_map(Q) I get a list of numbers from 0 to some integer which is twice the amount of vertices (minus 1). I would think the zero'th entry of d2v would be the vertex corresponding to the zero'th degree of freedom. Vice versa for v2d. This means in the dof to vertex map the vertices should be repeated so we can know which two degrees of freedom correspond to which vertex. This is not the case. I don't know how to interpret what the values returned by the maps mean.

How should I understand what the integers in these 2 mappings mean?

Answer 1

As the documentation says, "For mixed FunctionSpaces vertex index is offset with the number of dofs per vertex" (this also applies for VectorFunctionSpaces). So, the correct way to obtain the map between (mesh) vertex indices and dofs would be something like this:

from dolfin import *

mesh = UnitSquareMesh(5,5)
Q = VectorFunctionSpace(mesh, "CG", 1, dim=2)

v2d=vertex_to_dof_map(Q)
d2v=dof_to_vertex_map(Q)

v2d = v2d.reshape((-1, mesh.geometry().dim()))
d2v = d2v[xrange(0, len(d2v), 2)]/2

# Test
v = VertexFunction("size_t", mesh)
q = Function(Q)

v.array()[10] = 1.0
v.array()[25] = 1.0
q.vector()[v2d[10]] = 1.0
q.vector()[v2d[25]] = 1.0


plot(v)
plot(q)          # press m to visualize the mesh
interactive()

Comments

Perfect. Well explained and you even inadvertently answered a second question I had. It appears q.vector() is organized by degrees of freedom that have 1 or more values depending on the space it originates from.

If I had a mixed function space with three degrees of freedom per vertex then the v2d would be a list of triplets and q.vector requires replacing triplets of values. Is that correct?

---

Yes. Modifying the above example we can get this:

from dolfin import *

mesh = UnitSquareMesh(5,5)
P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V1 = VectorElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, P1*V1)

v2d=vertex_to_dof_map(W)
d2v=dof_to_vertex_map(W)

v2d = v2d.reshape((-1, 3))
d2v = d2v[xrange(0, len(d2v), 3)]/3

# Test
v = VertexFunction("size_t", mesh)
U = Function(W)

v.array()[10] = 1.0
v.array()[25] = 1.0
U.vector()[v2d[10]] = 1.0
U.vector()[v2d[25]] = 1.0

(p, q) = U.split()

plot(v)
plot(p)          # press m to visualize the mesh
plot(q)          # press m to visualize the mesh
interactive()