Permuting the dofs of DG spaces

[mscroggs]

The DG fem pipeline test (#634) is failing for tetrahedrons. I believe it only passes for other shapes because the Unit*Meshses are ordered nicely.

These tests pass if order_simplex is called on the mesh

For DG spaces, all the dofs are interior dofs, so the dof_permutation is moving all the dofs, including those associated with vertices. This may be the issue, or it may be related to #31 and how jumps etc are computed.

[chrisrichardson]

Probably we need to go back to the generated code and see what assumptions there are in that. I suppose it expects to receive two neighbouring cells, and then you specify which facets (in each cell) are the join. It may be that there is an implicit expectation that the dofs are aligned. I see this as a particular problem for quad/hex elements, since there is no possibility of the orientation always being consistent.

[mscroggs]

We have closed in on the issue by looking at the dg test term-by-term.

Taking u, v = TrialFunction(V), TestFunction(V), then the term u("-") * v("+") * dS does not agree if the mesh is ordered differently. This is one component of the dot(k("+") * jump(u, n), jump(v, n)) * dS term in the DG test.

Strangely, this term does agree if u and v are Functions.

The term inner(k("+") * avg(grad(v)), jump(u, n)) * dS agrees on differently ordered meshes, so it appears this problem only shows up when we have the product of jumps.

[mscroggs]

Have reduced this to the smallest test I can find that fails (

dolfinx/python/test/unit/fem/test_dof_permuting.py

Line 621 in 409d215

def test_plus_and_minus(cell_type, space_type, order):

):

def test_plus_and_minus(cell_type, space_type, order):
    """Test that f('+') and f('-') are computed correctly"""
    # Mesh:
    #   1
    #  /| <- cellA
    # 3-2
    # |/ <-- cellB
    # 0
    points = np.array([[0., -2.], [1., 2.], [1., 0.], [0., 0.]])

    cellB = (3, 2, 0)

    for cellA in [(3, 1, 2), (2, 3, 1)]:
        cells = np.array([cellA, cellB])
        mesh = Mesh(MPI.comm_world, cell_type, points, cells,
                    [], GhostMode.none)
        mesh.geometry.coord_mapping = fem.create_coordinate_map(mesh)

        V = FunctionSpace(mesh, (space_type, order))

        f = Function(V)
        f.interpolate(lambda x: x[0])

        a = f("+") * f("+") * dS
        b = f("-") * f("+") * dS

        v1 = fem.assemble_scalar(a)
        v2 = fem.assemble_scalar(b)

        print(v1, v2)

        assert np.isclose(v1, v2)

[mscroggs]

This issue is fixed in #702