Calculate Normal to Surface with Complex Geometries

I’m trying to find a vector for each cell in my mesh, showing the direction normal to the surface. Previously I’ve been using cylinders, so have been able to calculate it using basic geometry with the cell location, but am now using less regular geometries. Below is an example mesh. Each layer is 1 cell thick and is set up as a separate volume within the mesh. I’ve also created surfaces for the outer and inner faces.

I’ve tried using the FacetNormal command, but when I try to export the data to verify it, I always get an error. Using this code:

U = VectorFunctionSpace(mesh, “DG”, 0)

n0=FacetNormal(mesh)
n0_p=project(n0, U)

File(‘Results/n0_p.vtk’) << n0_p

I get the error:

UFLException: Integral of type cell cannot contain a ReferenceNormal.

Any ideas on how I can either make FacetNormal work, or an alternative way to get the result I need? Thanks in advance.

See: Obtain velocity normal to boundary - #2 by dokken
or: How to plot normal unit vector of faces in a 2D mesh? - #2 by dokken

Thank you so much, I’ll look at them now!