At the risk of reading too much into the question, I will first assume that you specifically need \nabla u\cdot\mathbf{n} on the boundary, rather than the full vector \nabla u (which avoids the extra complication of how to get the normal \mathbf{n} on `boundarymesh`

). If you have \nabla\cdot\mathbf{n} as a nodal `Function`

on `mesh`

, then you can interpolate it in the space `W`

on `boundarymesh`

relatively easily. The question is then “What is the best way to obtain a nodal function that represents \nabla u\cdot\mathbf{n}?” If I assume further that \nabla u\cdot\mathbf{n} is the flux of a PDE that was solved on `mesh`

, then the following is a good way to extract a nodal flux that is consistent with the interior solution (cf. Exercise 8 in Section 2.12 of Tom Hughes’s FEM book):

```
from fenics import *
mesh = UnitCubeMesh(15, 15, 15)
boundarymesh = BoundaryMesh(mesh, 'exterior')
n = FacetNormal(mesh)
V = FunctionSpace(mesh, "CG", 1)
u = Function(V)
W = FunctionSpace(boundarymesh, "CG", 1)
# Replace with the actual residual of the PDE used to compute u; in this case
# the flux for this problem will be the normal gradient, as in the MWE.
x = SpatialCoordinate(mesh)
f = sin(pi*x[0])
def PDE_res(u,v):
return inner(grad(u),grad(v))*dx - f*v*dx
# Dirichlet BC on boundary where you want to extract flux:
bcStr = ""
for i in range(0,3):
bcStr += "near(x["+str(i)+"],0.0)||near(x["+str(i)+"],1.0)"
if(i<2):
bcStr += "||"
bcStr = "("+bcStr+")"
uBC = DirichletBC(V,Constant(0.0),bcStr)
# Calculate something for u:
v = TestFunction(V)
res = PDE_res(TrialFunction(V),v)
solve(lhs(res)==rhs(res),u,uBC)
# Calculate a nodal flux consistent with the conservation law, i.e., find a
# nodal field `h` such that, if it were applied as Neumann data, it would have
# produced the solution `u`, which is considered given in this problem.
resFlux = PDE_res(u,v) - TrialFunction(V)*v*ds
h = Function(V)
# (Need to pass special `keep_diagonal` argument to LHS assembly, because
# the TrialFunction is only assembled on the boundary, and there will
# otherwise be a zero block in the sparse structure when we go to apply BCs
# on the interior, leading to errors.)
ABdry = assemble(lhs(resFlux),keep_diagonal=True)
bBdry = assemble(rhs(resFlux))
# The flux is defined only on boundary nodes, so set everything else to zero
# to obtain a well-posed linear algebra problem.
hBC = DirichletBC(V,Constant(0.0),"!"+bcStr)
hBC.apply(ABdry,bBdry)
solve(ABdry,h.vector(),bBdry)
# Interpolate the flux `h` defined on boundary nodes of `mesh` onto the space
# `W` defined on `boundarymesh`. (This step is exact, because `W` is the
# trace space of `V`.)
h.set_allow_extrapolation(True)
hbdry = interpolate(h, W)
# Test: Net flux should balance net source term from corresponding PDE
# (a.k.a. the "compatibility condition"):
print("Flux: "+str(assemble(hbdry*dx)))
print("Source term: "+str(assemble(f*dx(domain=mesh))))
# (This also verifies that flux `hbdry` interpolated on boundary mesh can be
# used in variational forms.)
```