How to get average facet normals for each node?

I/O
1

Hello,

I use FEniCS to solve a simple elasticity problem and for a certain surface of my geometry, I want to calculate the projection (i.e. dot product) of the displacement vectors to the local surface normal vector. What is the best way to achieve that?

My approach currently works only for vertices, but ideally I would like to do this for all nodes:

  1. I collect all facet indices of the surface I am interested in (with tag tag) via
   all_facets = np.array(list(facets(mesh)))
   surfacefacets = all_facets[np.where(boundaries.array() == tag)[0]]
   surfacevertexindices = set.union(*[set(facet.entities(0)) for facet in surfacefacets])
  1. Then for each vertex that is on that surface, I calculate the average normal vector (each vertex is connected to multiple facets, each with its own facet normal, so the average facet normal is OK for me) and project the displacement vector on the average normal vector:
   for vertexindex, vertexcoords in enumerate(mesh.coordinates()):
       if vertexindex in surfacevertexindices:
           vertexfacets = list(filter(lambda x: vertexindex in x.entities(0), surfacefacets))
           average_normal_vector = np.array([[facet.normal().x(),
                                              facet.normal().y(),
                                              facet.normal().z()]
                                              for facet in vertexfacets]
                                              ).mean(axis=0)
           average_normal_vector /= np.linalg.norm(average_normal_vector)
           displacement = u(vertexcoords)
           normal_displacement = np.dot(displacement, average_normal_vector)

This implementation seems clumsy to me and I like to have something that works not only for vertices, but for all nodes on the surface.

What is the FEniCS way to get the normal deformations for all vertices? Or even better, for all surfaces nodes?

Thanks for your help!

Edit: Changed misleading variable name.

2

After digging through the forum for any snippets that might apply to my problem, I found how to get node coordinates on a specific boundary.

With this, I came up with the following procedure:

  1. Get the facet indices for all facets on the boundary with tag tag:
    surfacefacetindices = np.where(boundaries.array() == tag)[0]
  2. Get the node (=dof) indices of the facet indices (since boundary is a 2d surface, dimension=2):
    dofs = V.dofmap().entity_closure_dofs(mesh, dimension, surfacefacetindices)
  3. Build a map dof index → dof coordinate triple (x, y, z):
    dof2coords = get_dofcoords_from_dof_index(V, list(set(dofs)))
    here, list(set(dofs)) yields the unique dofs (a dof can be on multiple facets, e.g. if a node is at a vertex where multiple facets touch).
    uniquedofs = dof2coords.keys()
  4. We need to know, what facets a dof belongs to in order to calculate the normal vectors at a dof. So, build a map dof index → facet index:
    dof2facet = get_facet_index_from_dof_index(V, mesh, surfacefacetindices)
  5. Finally, build a map dof index → average normal vector:
    dof2normal = get_average_normal_vector_from_dof(uniquedofs, dof2facet, mesh)

And here are the functions I used:

def get_dofcoords_from_dof_index(V, dofs):
    doftab = V.tabulate_dof_coordinates()
    dof2coords = dict()
    for dof in dofs:
        coords = list(doftab[dof])
        if not coords in dof2coords.values():
            dof2coords[dof] = coords
    for dof, coords in dof2coords.items():
        dof2coords[dof] = np.array(coords)
    return dof2coords

and

def get_facet_index_from_dof_index(V, mesh, facetindices, dimension=2):
    dof2facet = defaultdict(set)
    for facetindex in facetindices:
        dofs = V.dofmap().entity_closure_dofs(mesh, dimension, [facetindex])
        for dof in dofs:
            dof2facet[dof].add(facetindex)
    return dof2facet

and finally

def get_average_normal_vector_from_dof(dofs, dof2facet, mesh):
    dof2normal = dict()
    all_facets = np.array(list(facets(mesh)))
    for dof in dofs:
        doffacets = np.array(list(dof2facet[dof]))
        average_normal_vector = \ 
            np.array([[facet.normal().x(), facet.normal().y(), facet.normal().z()]
                      for facet in all_facets[doffacets]]).mean(axis=0)
        average_normal_vector /= np.linalg.norm(average_normal_vector)
        dof2normal[dof] = average_normal_vector
    return dof2normal

So now I am able to calculate the node normal vectors. But I am still a bit unclear on a few of the details:

  • When writing get_dofcoords_from_dof_index(), I did not expect that different dof indices would have identical coordinates. Why is that? This is especially strange to me, since I fed the function of unique dofs. A dof can be on multiple facets, I get that, but unique dofs should have unique coordinates. What am I missing here?
  • About entity_closure_dofs():
    • What does this method do, exactly? Why does entity_dofs() not work, for example? The manual did not tell me much.
    • Did I correctly interpret the 2nd argument (dimension) as the dimension of the boundary? That is, for a 2D mesh the boundary would be 1D and hence dimension=1?
  • In the above link to a thread, where dokken posted how to get node coordinates, he goes from cells to facets in order to build up the list of dofs: He uses mesh.topology and dofmap, but
    dofs = V.dofmap().entity_closure_dofs(mesh, dimension, surfacefacetindices)
    does the same thing without all the loops. Did I miss something?

Anyway, I’m glad I was able to figure this out, though it took me quite a while to get to it.

As a sidenote: I went through lots of try-and-error because I found little documentation in the tutorials, docstrings, or in the forum. Don’t get me wrong, the manual and tutorials are great if you can apply them directly to your problem, and I learned a lot by reading them. FEniCS could still be far more accessible, in my opinion, by docstrings that consistently told you what the methods do (in addition to its signatures and outputs, which of course helps, too).

3

Hi @steve,

Sorry this is a little late, since it seems you’ve already found a solution. But for your awareness, I wanted to point you to the following thread, which seems like it might accomplish what you’re looking for in a way that uses more FEniCS and less numpy: Obtain velocity normal to boundary - #2 by dokken

Regarding the first of your questions below, it does make sense to have multiple DOFs sharing the same coordinate. Since you mentioned elasticity, I assume your dependent variable is displacement (a vector). Thus for a given vertex (i), you will have three displacement components u_x^{(i)}, u_y^{(i)}, u_z^{(i)} at each mesh vertex, and thus three DOFs which all share the same coordinate.

Regarding your second and third questions, unfortunately, I don’t have an answer.

4

Hi @conpierce8,

Thanks for your reply! The post you linked to indeed solves my problem much more nicely and I was able to implement it in my code. I actually ran across this thread while searching, but only now that I’ve been able to solve my problem as described above, I understand that projecting the FacetNormals manually does indeed achieve the same goal.

Your remark about dofs sharing coordinate values helped me, too. It really cleared up how tabulate_dof_coordinates works and I was able to confirm my intuition with a toy example using a simple BoxMesh together with a FunctionSpace and a VectorFunctionSpace of different degrees. Plus, now I know what block size is useful for.

So thanks again for your reply!

1 Like