This is quite easy to debug if you add:
traction_facets = mesh.locate_entities_boundary(domain, 1, lambda x : np.sqrt((x[0]-grasp_node[0])**2 + (x[1]-grasp_node[1])**2 + (x[2]-grasp_node[2])**2) < force_radius)
mt = mesh.meshtags(domain, domain.geometry.dim-1, traction_facets, 1)
import dolfinx
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "facet_marker.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_meshtags(mt, domain.geometry)
and look at the marked facets in Paraview. Then you get:
Then, looking at your code, you want to find the facets, the dimension in
is wrong. It shouldn’t be 1, but domain.topology.dim-1, which is 2 in the case of a 3D geometry with tetrahedrons.
Changing this, you get:
