Mapping degrees of freedom between related meshes

Question

As the title suggests, I would like to map the degrees of freedom from a submesh to the original mesh, and from a boundarymesh to its parent mesh.
If the dof's correspond to vertices (CG1) I can do this by mapping and mapping back the corresponding vertices, using functions like (see second link):

dof_to_vertex_map(Vs)
submesh.data().array('parent_vertex_indices', 0)
boundarymesh.entity_map(0)
vertex_to_dof_map(V)

Now my question is whether there exists a similar piece of information that I can use to relate the degrees of freedom between meshes, not just the vertices? Reason for asking is that my problem requires a CG2 functionspace. Any pointing in the right direction will be helpful.

See also:
interpolate-works-functionspaces-with-degree-polynomials
vertex-mapping-from-submesh-boundarymesh-back-actual-mesh

Answer 1

If you have the corresponding cell objects in both meshes you can use the dofmap of the functions to get the degrees of freedom for each cell (call the dofmap.cell_dofs(cell.index()) method).

You will have to special case the mapping depending on the types of elements. The CG2 elements in 2D will i.e have the cell_dofs() method return six dof numbers, the first three degrees belong to the vertices and the next three correspond to facet centroids.

For CG2 in 2D, as far as I remember, the first facet centroid dof is on the facet facing (not connected to) the first vertex and so on. You can see in the FEniCS book for details on the dof numbering of various elements, or you can just ask dolfin for the coordinates of all dofs and infer the element ordering from that (http://fenicsproject.org/qa/3932/accessing-the-coordinates-of-a-degree-of-freedom)

Comments

Although it is probably possible to find matching cells and map the dof's, I found it eventually a lot easier to use the coordinates (thank you for mentioning them) of the dof's in both meshes, and map them that way. Not pretty, but simple and it works (for any order CG polynomial). I hope the soon-to-be-implemented MeshView will make life easier.

Perhaps the code is useful to someone else:

from dolfin import *
import numpy as np

degree = 3
n = 20

mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, 'CG', degree)
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)
leftboundary = LeftBoundary()
boundaries = FacetFunction('size_t', mesh)
boundaries.set_all(0)
leftboundary.mark(boundaries, 1)
ds = Measure('ds')[boundaries]

boundarymesh = BoundaryMesh(mesh, 'exterior')
left_mesh = SubMesh(boundarymesh, leftboundary)
L = FunctionSpace(left_mesh, 'CG', degree)
bf = Function(L)
# here goes some more code to find the correct bf
# but here a dummy function will do
bf.vector()[:] = np.linspace(0.0, 10.0, n * degree + 1)


# make a map of dofs from the boundarymesh to the original mesh
gdim = left_mesh.geometry().dim()
bmsh_dof_coordinates = L.dofmap().tabulate_all_coordinates(left_mesh).reshape(-1, gdim)
mesh_dof_coordinates = V.dofmap().tabulate_all_coordinates(mesh).reshape(-1, gdim)
bnd_to_glob_map = {}
for bnd_dof_nr, bnd_dof_coords in enumerate(bmsh_dof_coordinates):
    corresponding_dofs = [i for i, coords in enumerate(mesh_dof_coordinates) if np.array_equal(coords, bnd_dof_coords)]
    if len(corresponding_dofs) == 1:
        bnd_to_glob_map[bnd_dof_nr] = corresponding_dofs[0]
    else:
        raise NameError("Degrees of freedom not matching.")
print bnd_to_glob_map

# warp bf to V
bf_to_V = Function(V)
for bnd_dof, glob_dof in bnd_to_glob_map.iteritems():
    bf_to_V.vector()[glob_dof] = bf.vector().array()[bnd_dof]
plot(bf_to_V, interactive=True)