How to get global to local edge dof mapping for triangular P2 elements

I’m trying to understand the ordering of dofs for a solution 2nd order Lagrange elements and have a minimal example with 4 cells.

    Cells         Vertices

 ----- -----    4-----5-----6
| 3 / | \ 4 |   |   / | \   |
| / 1 | 2 \ |   | /   |   \ |
 ----- -----    1-----2-----3

Where the dof ordering per cell is the following (first row vertex dofs, second edge dofs)

7    2    7    5       Vtx dofs          Edge dofs (from mapping coordinates)
2    5   10    3
3    3    3   13    10-----3-----13     -11-- --14-
4    6   11   14     |   / | \   |   12|  8 /| \6  |15
8    4    8   15     | /   |   \ |     | /   4   \ |
9    1   12    6     7-----2-----5      --9-- --1--

I have understood the ordering of vertex dofs (by global vertex number), but I can’t seem to understand the algorithm for how the edge dofs are ordered (global to local (mesh) dof ordering). From the only brief reference in the 2011 UFC user’s manual 1.1 (if still valid), I interpret it as edge dofs are ordered in by the opposite vertex. For example, doing this from the global node numbering gives

  --6-- --8--
7|   / | \   |9    1  4  6  8
 | /2  1  4\ |     2  1  2  9
  --3-- --5--      3  5  7  4

which isn’t right. And no matter what I can’t seem to reconstruct the correct mapping. Is there some more up to date details on this, or maybe some easier way to do this?

For reference here is the code for the minimal example:

import dolfin
from fenics import *

with open('test.xml', 'wt') as f:
   f.write('<?xml version="1.0"?><dolfin xmlns:dolfin="http://fenicsproject.org"><mesh celltype="triangle" dim="2"><vertices size="6"><vertex index="0" x="0" y="0" z="0"/><vertex index="1" x="0.5" y="0" z="0"/><vertex index="2" x="1" y="0" z="0"/><vertex index="3" x="0" y="1" z="0"/><vertex index="4" x="0.5" y="1" z="0"/><vertex index="5" x="1" y="1" z="0"/></vertices><cells size="4"><cell index="0" v0="0" v1="1" v2="4"/><cell index="1" v0="2" v1="4" v2="1"/><cell index="2" v0="0" v1="4" v2="3"/><cell index="3" v0="5" v1="4" v2="2"/></cells></mesh></dolfin>')
mesh = Mesh('test.xml')

E = FiniteElement("P", mesh.ufl_cell(), 2)
V = FunctionSpace(mesh, E)
u = Function(V)
v = TestFunction(V)

a = dot(grad(u), grad(v))*dx
f = Constant(1.0)*v*dx
def boundary(x, on_boundary):
    return on_boundary
bc = DirichletBC(V, Constant(0), boundary)

solve(a - f == 0, u, bc)

h5 = dolfin.HDF5File(MPI.comm_world, 'test.h5', 'w')
h5.write(u,'/u')
h5.close()

print(u.vector().get_local())
print(V.tabulate_dof_coordinates())

The ordering in the UFC users manual refers to the local ordering of dofs on a cell.
If we add the following code to your problem:

from dolfin import *

with open('test.xml', 'wt') as f:
   f.write('<?xml version="1.0"?><dolfin xmlns:dolfin="http://fenicsproject.org"><mesh celltype="triangle" dim="2"><vertices size="6"><vertex index="0" x="0" y="0" z="0"/><vertex index="1" x="0.5" y="0" z="0"/><vertex index="2" x="1" y="0" z="0"/><vertex index="3" x="0" y="1" z="0"/><vertex index="4" x="0.5" y="1" z="0"/><vertex index="5" x="1" y="1" z="0"/></vertices><cells size="4"><cell index="0" v0="0" v1="1" v2="4"/><cell index="1" v0="2" v1="4" v2="1"/><cell index="2" v0="0" v1="4" v2="3"/><cell index="3" v0="5" v1="4" v2="2"/></cells></mesh></dolfin>')
mesh = Mesh('test.xml')

E = FiniteElement("P", mesh.ufl_cell(), 2)
V = FunctionSpace(mesh, E)
dofmap = V.dofmap()
x = V.tabulate_dof_coordinates()
for cell in cells(mesh):
    print(cell.index(), dofmap.cell_dofs(cell.index()), x[dofmap.cell_dofs(cell.index())])

We can sketch it to get the following picture.


You should not think too much about the global ordering, as it is optimized for data locality and will vary for each mesh (and depends on the number of processes).

However, you should consider the local ordering (obtained by dofmap.cell_dofs(i)).
As you see here, I’ve sketched each of the cells, and given them local numbering corresponding to the index of the cell dofs array.
We then observe that all of the local orderings follow the same pattern (local 3 opposite of local 0, local 4 opposite of local 1, local 5 opposite of local 2).

1 Like

Thank you for the detailed walk though, going step by step for each triangle helped and I think I got it right now.