Feature: add neighbor subface helper

src/t8_forest/t8_forest.cxx

@@ -1885,6 +1885,38 @@ t8_forest_leaf_face_neighbors (t8_forest_t forest, t8_locidx_t ltreeid, const t8
                                      pelement_indices, pneigh_eclass, forest_is_balanced, NULL, NULL);
 }
 
+int
+t8_forest_leaf_neighbor_subface (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *leaf, int face,
+                                 t8_eclass_t neighbor_tree_class, const t8_element_t *neighbor_leaf, int neighbor_face)
+{
+  t8_scheme const *scheme = t8_forest_get_scheme (forest);
+
+  t8_element_t *target = nullptr;
+  scheme->element_new (neighbor_tree_class, 1, &target);
+
+  int dummy;  // Can't pass a nullptr to t8_forest_element_face_neighbor below (see #2214)
+  t8_forest_element_face_neighbor (forest, ltreeid, leaf, target, neighbor_tree_class, face, &dummy);
+
+  int const num_children = scheme->element_get_num_face_children (neighbor_tree_class, neighbor_leaf, neighbor_face);
+
+  std::array<t8_element_t *, 4> children;  // assumes a 2:1 balanced forest
+  scheme->element_new (neighbor_tree_class, 4, children.begin ());
+
+  scheme->element_get_children_at_face (neighbor_tree_class, neighbor_leaf, neighbor_face, children.begin (),
+                                        num_children, nullptr);
+
+  int result = -1;
+  for (int i_child = 0; i_child < num_children; ++i_child) {
+    if (scheme->element_compare (neighbor_tree_class, target, children[i_child]) == 0) {
+      result = i_child;
+    }
+  }
+
+  scheme->element_destroy (neighbor_tree_class, 4, children.begin ());
+  scheme->element_destroy (neighbor_tree_class, 1, &target);
+  return result;
+}
+
 void
 t8_forest_print_all_leaf_neighbors (t8_forest_t forest)
 {

src/t8_forest/t8_forest_general.h

@@ -645,6 +645,21 @@ t8_forest_leaf_face_neighbors_ext (t8_forest_t forest, t8_locidx_t ltreeid, cons
                                    t8_locidx_t **pelement_indices, t8_eclass_t *pneigh_eclass, int forest_is_balanced,
                                    t8_gloidx_t *gneigh_tree, int *orientation);
 
+/** Compute the subface index for a coarser neighbor
+ * \param [in] forest       The forest. Must be committed and balanced.
+ * \param [in] ltreeid      A local tree id.
+ * \param [in] leaf         A leaf in \a ltreeid.
+ * \param [in] face         The face index of \a leaf to consider.
+ * \param [in] neighbor_tree_eclass The eclass of the neighbor element.
+ * \param [in] neighbor_leaf The leaf of \a forest on the other side of the face of index \a face of element \a leaf. Must be one level coarser than \a leaf.
+ * \param [in] neighbor_face The face index of \a neighbor_leaf (i.e. the dual face of \a face).
+ * \returns The index of the subface of \a neighbor_face which corresponds to \a face.
+ * \note This function is designed to be called after \ref t8_forest_leaf_face_neighbors_ext to complement its output.
+ */
+int
+t8_forest_leaf_neighbor_subface (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *leaf, int face,
+                                 t8_eclass_t neighbor_tree_class, const t8_element_t *neighbor_leaf, int neighbor_face);
+
 /** Exchange ghost information of user defined element data.
  * \param [in] forest       The forest. Must be committed.
  * \param [in] element_data An array of length num_local_elements + num_ghosts