How to get average facet normals for each node?

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).