Thanks for your reply

Sure thing. My dolfinx version → 0.7.1 from docker image. As for the example: lets assume that the problem is a simple mass transfer with Dbc = 0 on top and = 1 on the bottom edge :

```
from mpi4py import MPI
from dolfinx import mesh
import numpy as np
from dolfinx import default_scalar_type
from dolfinx.fem import (Constant, Function, FunctionSpace, Expression,
assemble_scalar, dirichletbc, form, locate_dofs_geometrical)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square
from dolfinx.plot import vtk_mesh
from mpi4py import MPI
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, FacetNormal, inner
import numpy as np
import pyvista
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
V = FunctionSpace(domain, ("Lagrange", 1))
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
def boundary_Dtop(x):
return np.isclose(x[1], 1)
dofs_Dzero = locate_dofs_geometrical(V, boundary_Dtop)
u_bc_zero = Constant(domain, 0.0)
bc_top = dirichletbc(u_bc_zero, dofs_Dzero, V)
def boundary_Dbottom(x):
return np.isclose(x[1], 0)
dofs_Dzero = locate_dofs_geometrical(V, boundary_Dbottom)
u_bc_zero = Constant(domain, 1.0)
bc_bottom = dirichletbc(u_bc_zero, dofs_Dzero, V)
bc = [bc_top, bc_bottom]
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v)) * dx
x = SpatialCoordinate(domain)
f = Constant(domain, default_scalar_type(0))
L = f * v * dx
problem = LinearProblem(a, L, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
```

I would like to calculate spatial distribution of the flux across the border.

The issue ive encountered is with this simple approach to calculate the grad(u) * n and interpolate it on some function space:

```
n = FacetNormal(domain)
flux_function = Function(V)
flux = Expression(inner(grad(uh),n), V.element.interpolation_points())
flux_function.interpolate(flux)
```

as i did in the past f.eg. when i wanted to assess the vector norm. But here i get error as stated in the previous post.

The error probably comes from the fact that the normal facet is not a vector field as i thought, but is some abstract representation calculated during the integration.

Never the less maybe someone have an idea how to solve this error or calculate the spatial distribution of flux accross the boundary.