### Normal component of a vector function on interior facets

`u`∈[

`V`]

^{2}, I want to find the projection of its trace on and perpendicular to the interior facets, $u\cdot n_F$

`u`·

`n`

_{F}, $\forall F\in\Gamma_{int}$∀

`F`∈Γ

_{int}. Rather naively I have constructed the following minimal example:

```
from dolfin import *
mesh = UnitSquareMesh(2, 2)
n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, "CG", 1)
u = interpolate(Constant((1, 1)), V)
Vscl = FunctionSpace(mesh, "CG", 1)
uscl, vscl = TrialFunction(Vscl), TestFunction(Vscl)
projection_eqn = (uscl*vscl*dx == (dot(u, n)*vscl)("+")*dS)
normal_pressure = Function(Vscl)
solve(projection_eqn, normal_pressure)
```

Would anyone have any advice or more rigorous methods of computing such a function?

### 1 Answer

`Discontinuous Lagrange Trace`

is the finite element family you are looking for. It is a function space living just on interior and exterior facets and being discontinuous piecewise polynomial.```
from dolfin import *
mesh = UnitSquareMesh(1, 1)
n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, "CG", 1)
u = interpolate(Constant((1, 1)), V)
Vscl = FunctionSpace(mesh, "Discontinuous Lagrange Trace", 0)
u_n, v = TrialFunction(Vscl), TestFunction(Vscl)
normal_trace = Function(Vscl)
projection_eqn = ((u_n * v)("+") * dS + (u_n * v) * ds ==
(dot(u, n) * v)("+") * dS + (dot(u, n) * v) * ds)
solve(projection_eqn, normal_trace)
print(normal_trace.vector()[:])
```

Mathematical note: In infinite dimensional setting, the most general space with well defined normal trace is `H(\text{div}, \Omega)` having trace in `H^{-1/2}(\partial \Omega)`. These `H^{-1/2}` are ugly, not even functions, and defined via Green's theorem. Fortunately, in finite dimension we work with `H(\text{div}, \Omega)` conforming subspaces and they have normal traces piecewise polynomial discontinuous.

For more see e.g.: "A Simple Introduction to the Mixed Finite Element Method" from Gatica.

Another tricky way would be to interpolate into Raviart-Thomas (degree=1) space and scale its degrees of freedom by facets area. You would get just normal component evaluations.

When looking at the output vector, do you know if there is a simple way to associate which value belongs to which facet in the mesh.

You cannot use eval on this FunctionSpace, so I was thinking of just iterating over cells:

```
facet_to_dofmap = {}
for c in cells(mesh):
i = c.index()
fac = c.entities(1)
dofs = Vscl.dofmap().cell_dofs(i)
for f,d in zip(fac, dofs):
facet_to_dofmap[f] = d
print(facet_to_dofmap)
```

But it won't work for higher order spaces, there are more dofs per facet.