def METHOD_NAME():
    """Check that discrete operator is consistent with normal
    interpolation between non-matching maps on non-affine geometries"""
    points = np.array([[0, 0, 0], [1, 0, 0], [0, 2, 0], [1, 2, 0],
                       [0, 0, 3], [1, 0, 3], [0, 2, 3], [1, 2, 3],
                       [0.5, 0, 0], [0, 1, 0], [0, 0, 1.5], [1, 1, 0],
                       [1, 0, 1.5], [0.5, 2, 0], [0, 2, 1.5], [1, 2, 1.5],
                       [0.5, 0, 3], [0, 1, 3], [1, 1, 3], [0.5, 2, 3],
                       [0.5, 1, 0], [0.5, -0.1, 1.5], [0, 1, 1.5], [1, 1, 1.5],
                       [0.5, 2, 1.5], [0.5, 1, 3], [0.5, 1, 1.5]], dtype=default_real_type)
    cells = np.array([range(len(points))], dtype=np.int32)
    cell_type = CellType.hexahedron
    domain = ufl.Mesh(element("Lagrange", cell_type.name, 2, rank=1))
    mesh = create_mesh(MPI.COMM_WORLD, cells, points, domain)
    W = VectorFunctionSpace(mesh, ("DG", 1))
    V = FunctionSpace(mesh, ("NCE", 4))
    w, v = Function(W), Function(V)
    w.interpolate(lambda x: x)
    v.interpolate(w)
    G = interpolation_matrix(W._cpp_object, V._cpp_object)
    G.assemble()
    # Compute global matrix vector product
    v_vec = Function(V)
    G.mult(w.vector, v_vec.vector)
    v_vec.x.scatter_forward()
    atol = 10 * np.finfo(default_real_type).resolution
    assert np.allclose(v_vec.x.array, v.x.array, atol=atol)
    s = assemble_scalar(form(ufl.inner(w - v, w - v) * ufl.dx))
    assert np.isclose(s, 0, atol=atol)
    G.destroy()