Add better error handling for locate_dofs_topological

cpp/dolfinx/fem/DirichletBC.cpp

@@ -23,13 +23,13 @@ using namespace dolfinx::fem;
 
 namespace
 {
-
 /// @brief Find the cell (local to process) and index of an entity
 /// (local to cell) for a list of entities.
+///
 /// @param[in] mesh The mesh
 /// @param[in] entities The list of entities
 /// @param[in] dim The dimension of the entities
-/// @returns A list of (cell_index, entity_index) pairs for each input
+/// @returns A list of `(cell_index, entity_index)` pairs for each input
 /// entity.
 std::vector<std::pair<std::int32_t, int>>
 find_local_entity_index(const mesh::Topology& topology,
@@ -55,15 +55,23 @@ find_local_entity_index(const mesh::Topology& topology,
 
   std::vector<std::pair<std::int32_t, int>> entity_indices;
   entity_indices.reserve(entities.size());
+  std::int32_t num_entities_local = e_to_c->num_nodes();
   for (std::int32_t e : entities)
   {
     // Get first attached cell
+    if (e >= num_entities_local)
+    {
+      throw std::out_of_range(
+          "Input entity " + std::to_string(e)
+          + " is larger than the number of entities on this process ("
+          + std::to_string(num_entities_local) + ").");
+    }
     assert(e_to_c->num_links(e) > 0);
-    const int cell = e_to_c->links(e).front();
 
     // Get local index of facet with respect to the cell
+    std::int32_t cell = e_to_c->links(e).front();
     auto entities_d = c_to_e->links(cell);
-    auto it = std::find(entities_d.begin(), entities_d.end(), e);
+    auto it = std::ranges::find(entities_d, e);
     assert(it != entities_d.end());
     std::size_t entity_local_index = std::distance(entities_d.begin(), it);
     entity_indices.emplace_back(cell, entity_local_index);
@@ -74,7 +82,7 @@ find_local_entity_index(const mesh::Topology& topology,
 //-----------------------------------------------------------------------------
 
 /// Find all DOFs on this process that have been detected on another
-/// process
+/// process.
 /// @param[in] comm A symmetric communicator
 /// @param[in] map The index map with the dof layout
 /// @param[in] bs_map The block size of the index map, i.e. the dof

python/test/unit/fem/test_bcs.py

@@ -406,3 +406,16 @@ def test_blocked_dof_ownership(shape):
         owned_sub_dofs = boundary_dofs_V[boundary_dofs_V < num_owned_blocked * bs]
         assert len(owned_sub_dofs) == num_owned_sub
         assert len(unrolled_dofs_sub) == len(boundary_dofs_V)
+
+
+def test_bc_index_out_of_range():
+    mesh = create_unit_square(MPI.COMM_WORLD, 4, 4)
+    mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
+    facet_map = mesh.topology.index_map(mesh.topology.dim - 1)
+    num_facets_on_proc = facet_map.size_local + facet_map.num_ghosts
+    facets = np.arange(num_facets_on_proc + 1, dtype=np.int32)
+    V = functionspace(mesh, ("Lagrange", 1))
+    with pytest.raises(
+        IndexError, match=r".*is larger than the number of entities on this process.*"
+    ):
+        locate_dofs_topological(V, mesh.topology.dim - 1, facets)