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:
- Get the facet indices for all facets on the boundary with tag
tag
:
surfacefacetindices = np.where(boundaries.array() == tag)[0]
- 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)
- 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()
- 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)
- 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 hencedimension=1
?
- What does this method do, exactly? Why does
- 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).