In the mixed-poisson example of Dolfinx demo, the way normal flux BCs (which is essential in this context) are imposed is by interpolating a position dependent function h(x) into the Hdiv-like space (e.g RT, BDM) and then passing it to the fem.dirichletbc with the corresponding dofs and Hdiv part of the mixed space.
Being specific, in the aforementioned example, only y direction value are set because the \vec{n} = [0, \pm 1]^T. Mathematically we generalise this by taking \vec{g}(\vec{x}, \vec{n}) = g_n(\vec{x}) \vec{n}(\vec{x}), with \vec{n}(\vec{x}) spatially dependent (think in curved boundary). We could think in programming a vectorial function like this to imporse the normal flux g_n given.
In legacy Fenics, I used to program that using a UserExpression, implementing the method eval_cell with x and ufc_cell as arguments (the latter contains information about the normal vector) full code
def eval_cell(self, values, x, ufc_cell):
if(ufc_cell.local_facet>-1):
cell = df.Cell(self.mesh, ufc_cell.index)
n = cell.normal(ufc_cell.local_facet)
values[:] = self.g*n.array()[:self.mesh.geometric_dimension()]
Is there a way to do the same dolfinx?