Function test_plus_minus_vector

python/test/unit/fem/test_element_integrals.py:401–448

Test that ('+') and ('-') match up with the correct DOFs for DG functions.

(cell_type, pm1, pm2, dtype)

@parametrize_cell_types
@parametrize_dtypes
def test_plus_minus_vector(cell_type, pm1, pm2, dtype):
    """Test that ('+') and ('-') match up with the correct DOFs for DG functions."""
    xtype = np.real(dtype(0)).dtype
    results = []
    orders = []
    spaces = []
    for count in range(3):
        for agree in [True, False]:
            # Two cell mesh with randomly numbered points
            mesh, order = two_unit_cells(cell_type, xtype, agree, return_order=True)
            if cell_type in [CellType.interval, CellType.triangle, CellType.tetrahedron]:
                V = functionspace(mesh, ("DG", 1))
            else:
                V = functionspace(mesh, ("DQ", 1))

            # Assemble vectors with combinations of + and - for a few
            # different numberings
            f = Function(V, dtype=dtype)
            f.interpolate(lambda x: x[0] - 2 * x[1])
            v = ufl.TestFunction(V)
            a = form(ufl.inner(f(pm1), v(pm2)) * ufl.dS, dtype=dtype)
            result = assemble_vector(a)
            spaces.append(V)
            results.append(result.array)
            orders.append(order)

    # Check that the above vectors all have the same values as the first
    # one, but permuted due to differently ordered dofs
    dofmap0 = spaces[0].mesh.geometry.dofmaps[0]
    for result, space in zip(results[1:], spaces[1:]):
        # Get the data relating to two results
        dofmap1 = space.mesh.geometry.dofmaps[0]

        # For each cell
        for cell in range(2):
            # For each point in cell 0 in the first mesh
            for dof0, point0 in zip(spaces[0].dofmap.cell_dofs(cell), dofmap0[cell]):
                # Find the point in the cell 0 in the second mesh
                for dof1, point1 in zip(space.dofmap.cell_dofs(cell), dofmap1[cell]):
                    if np.allclose(
                        spaces[0].mesh.geometry.x[point0], space.mesh.geometry.x[point1]
                    ):
                        break
                else:
                    # If no matching point found, fail
                    assert False

                assert np.isclose(results[0][dof0], result[dof1], atol=1.0e-6)


@pytest.mark.skip_in_parallel