Compute normal vector to cells

Question

I have a surface mesh in 3D and I want to define a vertex_normals matrix that represents the average of the normal vectors (properly aligned) of all the cells that intersect with each vertex. My approach is not very fast and I wonder if there's a better (or preexisting) way to do it.

My approach has been the following:

• make a list (vertex_neighbours) of all the the cells that contain each vertex.
• make a list (normal_faces) of the normal vector to each cell (cross product of vertex coordinates, not sure if fenics has this functionality already).
• combine both lists to generate the desired vertex_normals:

vertex_normals = [np.sum(normal_faces[vertex_neighbours[vertex]])/
len(vertex_neighbours[vertex]) for vertex in xrange(num_vertices)]

This is unfortunately not very fast for large meshes, and so far I haven't been able to figure out a way to vectorize it.

Comments

Which of the three parts of the algorithm performs poorly?

---

The third part, specifically the code for vertex_normals. I found a way to vectorize the other two parts, making them about 40x faster than with a for-loop, but for the third part the best I could do was list comprehension.

---

This is rather pythonic than FEniCS question. Nevertheless, does this help vertex_normals = np.average(normal_faces[vertex_neighbours], axis=1)?

---

I'll give it a try, thanks. I was wondering if there was a preexisting function, there's something called FacetNormal, but I couldn't figure out how to use it.

---

As I know

n = FacetNormal(mesh)

can be used only in forms as it does not have eval method - FFC generates code evaluating it directly. Maybe it can be projected or interpolated to some suitable space - is there a vector space having one degree of freedom corresponding to each facet?

---

FunctionSpace(mesh, "RT", 1)

---

From the picture that looks about right. Could you give me an example?, say, suppose you wanted to find the normal vector to cell 7.

---

V = FunctionSpace(mesh, "RT", 1)
n = FacetNormal(mesh)
n = project(n, V)

Then you need to access correct dofs of n. You can use V.dofmap().dof_to_vertex_map() or parameters['reorder_dofs_serial'] = False. See this question.

---

I'm sorry, this is of course wrong. Here DOFs are associated with facets so obviously you can't access them by dof_to_vertex_map(). But some other helper function within V.dofmap() could help - check DofMap documentation and dolfin/fem/DofMap.cpp in case of emergency. Alternatively you can do (in serial only) parameters['reorder_dofs_serial'] = False which should ensure that dofs are numbered by facet indeces.