diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl
index 9011ac91215..4843f211089 100644
--- a/src/semidiscretization/semidiscretization_coupled.jl
+++ b/src/semidiscretization/semidiscretization_coupled.jl
@@ -603,10 +603,14 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{
     i_cell_start, i_cell_step = index_to_start_step_2d(indices[1], axes(mesh_other, 1))
     j_cell_start, j_cell_step = index_to_start_step_2d(indices[2], axes(mesh_other, 2))
 
-    i_cell = i_cell_start
-    j_cell = j_cell_start
+    # We need indices starting at 1 for the handling of `i_cell` etc.
+    Base.require_one_based_indexing(cells)
+
+    @threaded for i in eachindex(cells)
+        cell = cells[i]
+        i_cell = i_cell_start + (i - 1) * i_cell_step
+        j_cell = j_cell_start + (i - 1) * j_cell_step
 
-    for cell in cells
         i_node = i_node_start
         j_node = j_node_start
         element_id = linear_indices[i_cell, j_cell]
@@ -628,9 +632,6 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{
             i_node += i_node_step
             j_node += j_node_step
         end
-
-        i_cell += i_cell_step
-        j_cell += j_cell_step
     end
 end