Feature element boundary [3/3]

src/t8_forest/t8_forest.cxx

@@ -1993,6 +1993,7 @@ t8_forest_element_is_leaf (const t8_forest_t forest, const t8_element_t *element
 {
   T8_ASSERT (t8_forest_is_committed (forest));
   T8_ASSERT (t8_forest_tree_is_local (forest, local_tree));
+  T8_ASSERT (element != NULL);
 
   /* We get the array of the tree's elements and then search in the array of elements for our 
    * element candidate. */
@@ -2023,6 +2024,75 @@ t8_forest_element_is_leaf (const t8_forest_t forest, const t8_element_t *element
   return (scheme->element_is_equal (tree_class, element, check_element));
 }
 
+int
+t8_forest_leaf_is_boundary (const t8_forest_t forest, t8_locidx_t local_tree, const t8_element_t *leaf, int face)
+{
+  T8_ASSERT (t8_forest_is_committed (forest));
+  T8_ASSERT (t8_forest_element_is_leaf (forest, leaf, local_tree));
+  T8_ASSERT (leaf != NULL);
+
+  const t8_eclass_t tree_class = t8_forest_get_tree_class (forest, local_tree);
+  const t8_eclass_scheme_c *scheme = t8_forest_get_eclass_scheme (forest, tree_class);
+
+  /* Check whether this leaf is at the boundary of its tree. */
+  const int is_root_boundary = scheme->t8_element_is_root_boundary (leaf, face);
+
+  if (is_root_boundary) {
+    /* This leaf is at a tree's boundary.
+     * If the respective tree face is at the domain boundary,
+     * then the element is as well.
+     * If the tree face is not at the domain boundary, the element's face
+     * could still be at an inner boundary. */
+    const int cmesh_face = scheme->t8_element_tree_face (leaf, face);
+    const t8_cmesh_t cmesh = t8_forest_get_cmesh (forest);
+    const t8_locidx_t cmesh_local_tree = t8_forest_ltreeid_to_cmesh_ltreeid (forest, local_tree);
+    int tree_boundary = t8_cmesh_tree_face_is_boundary (cmesh, cmesh_local_tree, cmesh_face);
+    if (tree_boundary) {
+      return 1;
+    }
+  }
+
+  /* This leaf is not at the tree's boundary.
+    * If the forest has holes, we need to check whether this leaf is at an internal boundary.*/
+
+  if (!forest->incomplete_trees) {
+    /* The forest has no holes, thus the leaf cannot be a boundary leaf. */
+    return 0;
+  }
+
+  /*
+   * The remaining code handles the case that the forest has holes and the element may thus 
+   * be an inner boundary.
+   * This case is not yet support due to issue #825. Hence, we currently abort.
+   * Once the issue is resolved, the abort message can be removed and the code should work.
+   * */
+
+  SC_ABORT ("This forest has holes and a computation of boundary elements is not supported. Once "
+            "https://github.com/DLR-AMR/t8code/issues/825 is resolved, the function will be available.\n");
+
+  /* we need to compute the face neighbors to know whether the element is a boundary element. */
+  const int is_balanced = t8_forest_is_balanced (forest);
+  int num_neighbors;
+  t8_element_t **neighbor_leaves;
+  t8_locidx_t *pelement_indices;
+  t8_eclass_scheme_c *pneigh_scheme;
+  /* The forest has holes, the leaf could lie inside a tree but its neighbor was deleted. */
+  t8_forest_leaf_face_neighbors (forest, local_tree, leaf, &neighbor_leaves, face, NULL, &num_neighbors,
+                                 &pelement_indices, &pneigh_scheme, is_balanced);
+
+  if (num_neighbors == 0) {
+    /* The element has no neighbors, it is a boundary element. */
+    return 1;
+  }
+  else {
+    /* If neighbors were found, these arrays were allocated and need clean-up. */
+    T8_FREE (neighbor_leaves);
+    T8_FREE (pelement_indices);
+    /* This leaf is not a boundary leaf. */
+    return 0;
+  }
+}
+
 /* Check if an element is owned by a specific rank */
 int
 t8_forest_element_check_owner (t8_forest_t forest, t8_element_t *element, t8_gloidx_t gtreeid, t8_eclass_t eclass,

src/t8_forest/t8_forest_geometrical.h

@@ -170,6 +170,17 @@ void
 t8_forest_element_face_normal (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, int face,
                                double normal[3]);
 
+/** Compute whether a given face of a given leaf in a forest lies at the domain boundary.
+ * This includes inner boundaries if the forest has deleted elements.
+ * \param [in]      forest    The forest.
+ * \param [in]      local_tree A local tree of \a forest.
+ * \param [in]      leaf      A leaf element of \a local_tree in \a forest.
+ * \param [in]      face      The face number of a face of \a leaf to check.
+ * \return True (non-zero) if face \a face of \a leaf lies at the domain boundary.
+*/
+int
+t8_forest_leaf_is_boundary (const t8_forest_t forest, t8_locidx_t local_tree, const t8_element_t *leaf, int face);
+
 T8_EXTERN_C_END ();
 
 #endif /* !T8_FOREST_GEOMETRICAL_H */

test/CMakeLists.txt

@@ -102,6 +102,7 @@ add_t8_test( NAME t8_gtest_balance_parallel             SOURCES t8_gtest_main.cx
 add_t8_test( NAME t8_gtest_forest_commit_parallel       SOURCES t8_gtest_main.cxx t8_forest/t8_gtest_forest_commit.cxx )
 add_t8_test( NAME t8_gtest_forest_face_normal_serial    SOURCES t8_gtest_main.cxx t8_forest/t8_gtest_forest_face_normal.cxx )
 add_t8_test( NAME t8_gtest_element_is_leaf_serial       SOURCES t8_gtest_main.cxx t8_forest/t8_gtest_element_is_leaf.cxx )
+add_t8_test( NAME t8_gtest_element_is_boundary_parallel SOURCES t8_gtest_main.cxx t8_forest/t8_gtest_element_is_boundary.cxx )
 add_t8_test( NAME t8_gtest_partition_data_parallel      SOURCES t8_gtest_main.cxx t8_forest/t8_gtest_partition_data.cxx )
 
 add_t8_test( NAME t8_gtest_permute_hole_serial          SOURCES t8_gtest_main.cxx t8_forest_incomplete/t8_gtest_permute_hole.cxx )

test/Makefile.am

@@ -79,6 +79,7 @@ t8code_googletest_programs = \
   test/t8_forest/t8_gtest_forest_commit \
   test/t8_forest/t8_gtest_balance \
   test/t8_forest/t8_gtest_element_is_leaf \
+  test/t8_forest/t8_gtest_element_is_boundary \
   test/t8_IO/t8_gtest_vtk_reader \
   test/t8_IO/t8_gtest_vtk_writer \
   test/t8_forest_incomplete/t8_gtest_permute_hole \
@@ -306,6 +307,10 @@ test_t8_forest_t8_gtest_element_is_leaf_SOURCES = \
   test/t8_gtest_main.cxx \
   test/t8_forest/t8_gtest_element_is_leaf.cxx
 
+test_t8_forest_t8_gtest_element_is_boundary_SOURCES = \
+  test/t8_gtest_main.cxx \
+  test/t8_forest/t8_gtest_element_is_boundary.cxx
+
 test_t8_forest_incomplete_t8_gtest_permute_hole_SOURCES = \
   test/t8_gtest_main.cxx \
   test/t8_forest_incomplete/t8_gtest_permute_hole.cxx
@@ -574,6 +579,10 @@ test_t8_forest_t8_gtest_element_is_leaf_LDADD = $(t8_gtest_target_ld_add)
 test_t8_forest_t8_gtest_element_is_leaf_LDFLAGS = $(t8_gtest_target_ld_flags)
 test_t8_forest_t8_gtest_element_is_leaf_CPPFLAGS = $(t8_gtest_target_cpp_flags)
 
+test_t8_forest_t8_gtest_element_is_boundary_LDADD = $(t8_gtest_target_ld_add)
+test_t8_forest_t8_gtest_element_is_boundary_LDFLAGS = $(t8_gtest_target_ld_flags)
+test_t8_forest_t8_gtest_element_is_boundary_CPPFLAGS = $(t8_gtest_target_cpp_flags)
+
 test_t8_IO_t8_gtest_vtk_reader_LDADD = $(t8_gtest_target_ld_add)
 test_t8_IO_t8_gtest_vtk_reader_LDFLAGS = $(t8_gtest_target_ld_flags)
 test_t8_IO_t8_gtest_vtk_reader_CPPFLAGS = $(t8_gtest_target_cpp_flags)
@@ -683,6 +692,7 @@ test_t8_forest_t8_gtest_ghost_and_owner_CPPFLAGS += $(t8_gtest_target_mpi_cpp_fl
 test_t8_forest_t8_gtest_forest_commit_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_forest_t8_gtest_balance_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_forest_t8_gtest_element_is_leaf_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
+test_t8_forest_t8_gtest_element_is_boundary_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_IO_t8_gtest_vtk_reader_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_forest_incomplete_t8_gtest_permute_hole_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_forest_incomplete_t8_gtest_recursive_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)