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())