Get velocity field at a boundary

dolfinx
2

You are using a function space with a block size (this is due to the usage of vectorelement). This means that the dofs are not unrolled.
Consider this MWE:

import dolfinx
import numpy as np
from dolfinx import fem
from dolfinx.fem import Function, FunctionSpace, locate_dofs_topological
from mpi4py import MPI
from ufl import VectorElement

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

def boundary(x):
    return x[0]<0.5

boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim-1, boundary)

v_cg2 = VectorElement('Lagrange', mesh.ufl_cell(), 2) # Order 2
V = FunctionSpace(mesh, v_cg2)
u = fem.Function(V)

dof_coords = V.tabulate_dof_coordinates()
bs = V.dofmap.bs
boundary_dofs = locate_dofs_topological(V,mesh.topology.dim-1, boundary_facets)


def velocity(x):
    return (np.sin(x[0]), x[1])
# Create initial condition
u_n = Function(V)
u_n.name = "u_n"
u_n.interpolate(velocity)

for dof in boundary_dofs:
    coord = dof_coords[dof]
    dofs = np.array([dof*bs+i for i in range(bs)], dtype=np.int32)
    exact = [np.sin(coord[0]), coord[1]]
    print("Coordinate", coord, "Function", u_n.x.array[dofs], "Exact", exact[0], exact[1])
    for i in range(bs):
        assert np.isclose(u_n.x.array[dofs[i]], exact[i])

Note that it is often useful to try to reduce your problem to something on a built-in mesh, and with a known function (usually something found through interpolation).