diff --git a/src/t8_forest/t8_forest.cxx b/src/t8_forest/t8_forest.cxx index 0a8fb6b29e..34592f67d8 100644 --- a/src/t8_forest/t8_forest.cxx +++ b/src/t8_forest/t8_forest.cxx @@ -2047,6 +2047,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. */ @@ -2076,6 +2077,75 @@ t8_forest_element_is_leaf (const t8_forest_t forest, const t8_element_t *element return (scheme->t8_element_equal (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, diff --git a/src/t8_forest/t8_forest_geometrical.h b/src/t8_forest/t8_forest_geometrical.h index 1accee5a5c..8353749248 100644 --- a/src/t8_forest/t8_forest_geometrical.h +++ b/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 */ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index ef9f9ba71a..b4fb049c97 100644 --- a/test/CMakeLists.txt +++ b/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 ) diff --git a/test/Makefile.am b/test/Makefile.am index 6ee6f4d17a..4fd16b846e 100644 --- a/test/Makefile.am +++ b/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) diff --git a/test/t8_forest/t8_gtest_element_is_boundary.cxx b/test/t8_forest/t8_gtest_element_is_boundary.cxx new file mode 100644 index 0000000000..51fbba8310 --- /dev/null +++ b/test/t8_forest/t8_gtest_element_is_boundary.cxx @@ -0,0 +1,255 @@ +/* + This file is part of t8code. + t8code is a C library to manage a collection (a forest) of multiple + connected adaptive space-trees of general element classes in parallel. + + Copyright (C) 2024 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +#include +#include +#include +#include +#include +#include +#include "test/t8_cmesh_generator/t8_cmesh_example_sets.hxx" +#include + +/* In this test we check the t8_forest_element_is_boundary function. + * Iterating over all cmesh test cases, we create a uniform and an adaptive forest. + * For each forest, we check that for each element t8_forest_element_is_boundary returns true + * if and only if the element is at the domain boundary. + */ + +/* Maximum uniform level for forest. */ +#define T8_IS_BOUNDARY_MAX_LVL 3 + +/* Adapt a forest such that always the first child of a + * family is refined and no other elements. This results in a highly + * imbalanced forest. */ +static int +t8_test_adapt_first_child (t8_forest_t forest, t8_forest_t forest_from, t8_locidx_t which_tree, t8_locidx_t lelement_id, + t8_eclass_scheme_c *ts, const int is_family, const int num_elements, + t8_element_t *elements[]) +{ + int level = ts->t8_element_level (elements[0]); + + /* we set a maximum refinement level as forest user data */ + int maxlevel = *(int *) t8_forest_get_user_data (forest); + if (level >= maxlevel) { + /* Do not refine after the maxlevel */ + return 0; + } + int child_id = ts->t8_element_child_id (elements[0]); + if (child_id == 1) { + return 1; + } + return 0; +} + +/* In a quad forest remove the first and fourth child of each element. + * This will reside in a forest where each element and each face is a boundary element. + * */ +static int +t8_test_adapt_quad_remove_first_and_fourth_child (t8_forest_t forest, t8_forest_t forest_from, t8_locidx_t which_tree, + t8_locidx_t lelement_id, t8_eclass_scheme_c *ts, const int is_family, + const int num_elements, t8_element_t *elements[]) +{ + /* For the purpose of this test, this function should only get called with quad elements. */ + SC_CHECK_ABORT (ts->t8_element_shape (elements[0]) == T8_ECLASS_QUAD, + "Special test adapt function must only be used on quad elements.\n"); + int child_id = ts->t8_element_child_id (elements[0]); + /* Remove child_id 0 and 3, do not change any other element. */ + if (child_id == 0 || child_id == 3) { + return -2; + } + return 0; +} + +class element_is_boundary: public testing::TestWithParam> { + protected: + void + SetUp () override + { + /* Construct a cmesh */ + const int level = std::get<0> (GetParam ()); + cmesh = std::get<1> (GetParam ())->cmesh_create (); + if (t8_cmesh_is_empty (cmesh)) { + /* forest_commit does not support empty cmeshes, we skip this case */ + GTEST_SKIP (); + } + /* Build the default scheme (TODO: Test this with all schemes) */ + scheme = t8_scheme_new_default_cxx (); + forest = t8_forest_new_uniform (cmesh, scheme, level, 0, sc_MPI_COMM_WORLD); + t8_forest_ref (forest); + int maxlevel = 7; + const int recursive_adapt = 1; + forest_adapt = t8_forest_new_adapt (forest, t8_test_adapt_first_child, recursive_adapt, 0, &maxlevel); + /* Since the cmesh may be partitioned when building forest and forest_adapt, the old + * cmesh might get deleted and replaced by a new partitioned cmesh. + * We need to update our internal cmesh pointer to accommodate. */ + cmesh = t8_forest_get_cmesh (forest_adapt); + } + + void + TearDown () override + { + if (t8_cmesh_is_empty (cmesh)) { + t8_cmesh_destroy (&cmesh); + } + else { + t8_forest_unref (&forest); + t8_forest_unref (&forest_adapt); + } + } + + t8_cmesh_t cmesh; + t8_forest_t forest, forest_adapt; + t8_scheme_cxx_t *scheme; +}; + +void +t8_test_element_is_boundary_for_forest (t8_forest_t forest, t8_cmesh_t cmesh, + const int each_element_face_is_expected_boundary) +{ + const t8_locidx_t num_local_trees = t8_forest_get_num_local_trees (forest); + + for (t8_locidx_t itree = 0; itree < num_local_trees; ++itree) { + const t8_locidx_t num_elements_in_tree = t8_forest_get_tree_num_elements (forest, itree); + const t8_eclass_t tree_class = t8_forest_get_tree_class (forest, itree); + const t8_eclass_scheme_c *scheme = t8_forest_get_eclass_scheme (forest, tree_class); + /* Iterate over all the tree's leaf elements, check whether the leaf + * is correctly identified by t8_forest_element_is_boundary, + * build its parent and its first child (if they exist), and verify + * that t8_forest_element_is_boundary returns false. */ + for (t8_locidx_t ielement = 0; ielement < num_elements_in_tree; ++ielement) { + const t8_element_t *leaf_element = t8_forest_get_element_in_tree (forest, itree, ielement); + const int num_element_faces = scheme->t8_element_num_faces (leaf_element); + for (int iface = 0; iface < num_element_faces; ++iface) { + /* Iterate over all faces */ + int face_is_at_boundary = 0; + if (!each_element_face_is_expected_boundary) { + /* Manually check whether this element is at the boundary */ + if (scheme->t8_element_is_root_boundary (leaf_element, iface)) { + /* This face is at a tree boundary, so it might be at the domain boundary. */ + /* We check whether the tree is at the domain boundary. */ + const int tree_face = scheme->t8_element_tree_face (leaf_element, iface); + const t8_locidx_t cmesh_local_tree = t8_forest_ltreeid_to_cmesh_ltreeid (forest, itree); + if (t8_cmesh_tree_face_is_boundary (cmesh, cmesh_local_tree, tree_face)) { + /* This element is at the domain boundary. */ + face_is_at_boundary = 1; + } + } + } + else { + /* each element and face is an expected boundary */ + face_is_at_boundary = 1; + } + EXPECT_EQ (t8_forest_leaf_is_boundary (forest, itree, leaf_element, iface), face_is_at_boundary); + } + } + } +} + +TEST_P (element_is_boundary, element_is_boundary) +{ + t8_test_element_is_boundary_for_forest (forest, cmesh, 0); +} + +TEST_P (element_is_boundary, element_is_boundary_adapt) +{ + t8_test_element_is_boundary_for_forest (forest_adapt, cmesh, 0); +} + +TEST (element_is_boundary, quad_forest_with_holes) +{ + /* This test is deactivated as long as Ghost does not work in combination + * with deleted elements. + * Once this is fixed, reactivate this test. + * See https://github.com/DLR-AMR/t8code/issues/825 + */ + GTEST_SKIP (); + /* Create a 10 x 5 2D brick cmesh, periodic in x direction. */ + t8_cmesh_t cmesh = t8_cmesh_new_brick_2d (10, 5, 1, 0, sc_MPI_COMM_WORLD); + + t8_scheme_cxx_t *scheme = t8_scheme_new_default_cxx (); + t8_forest_t forest = t8_forest_new_uniform (cmesh, scheme, T8_IS_BOUNDARY_MAX_LVL, 0, sc_MPI_COMM_WORLD); + t8_forest_t forest_adapt = t8_forest_new_adapt (forest, t8_test_adapt_quad_remove_first_and_fourth_child, 0, 1, NULL); + + t8_forest_write_vtk (forest_adapt, "test_quad_w_holes"); + + t8_test_element_is_boundary_for_forest (forest_adapt, cmesh, 1); + t8_forest_unref (&forest_adapt); +} + +/* This test class creates a single tree cmesh for each eclass and + * builds a uniform level 0 forest on it. + * For the single element in that forest all faces lie on the boundary. */ +class element_is_boundary_known_boundary: public testing::TestWithParam { + protected: + void + SetUp () override + { + eclass = GetParam (); + cmesh = t8_cmesh_new_from_class (eclass, sc_MPI_COMM_WORLD); + t8_scheme_cxx_t *scheme = t8_scheme_new_default_cxx (); + forest = t8_forest_new_uniform (cmesh, scheme, 0, 0, sc_MPI_COMM_WORLD); + } + + void + TearDown () override + { + t8_forest_unref (&forest); + } + + t8_cmesh_t cmesh; + t8_forest_t forest; + t8_eclass_t eclass; +}; + +/* For a level 0 forest of a single class cmesh all faces should + * be at the boundary. */ +TEST_P (element_is_boundary_known_boundary, level_0) +{ + const t8_locidx_t num_elements = t8_forest_get_local_num_elements (forest); + if (num_elements > 0) { + T8_ASSERT (num_elements == 1); + const t8_element_t *element = t8_forest_get_element_in_tree (forest, 0, 0); + const t8_eclass_scheme_c *scheme = t8_forest_get_eclass_scheme (forest, eclass); + const int num_faces = scheme->t8_element_num_faces (element); + for (int iface = 0; iface < num_faces; ++iface) { + EXPECT_TRUE (t8_forest_leaf_is_boundary (forest, 0, element, iface)); + } + } +} + +/* Define a lambda to beautify gtest output for tuples . + * This will set the correct level and cmesh name as part of the test case name. */ +auto pretty_print_level_and_cmesh_params + = [] (const testing::TestParamInfo> &info) { + std::string name = std::string ("Level_") + std::to_string (std::get<0> (info.param)); + std::string cmesh_name; + std::get<1> (info.param)->param_to_string (cmesh_name); + name += std::string ("_") + cmesh_name; + return name; + }; + +INSTANTIATE_TEST_SUITE_P (t8_gtest_element_is_boundary, element_is_boundary, + testing::Combine (testing::Range (0, T8_IS_BOUNDARY_MAX_LVL), AllCmeshsParam), + pretty_print_level_and_cmesh_params); + +INSTANTIATE_TEST_SUITE_P (t8_gtest_element_is_boundary, element_is_boundary_known_boundary, AllEclasses, print_eclass); diff --git a/test/t8_forest/t8_gtest_element_is_leaf.cxx b/test/t8_forest/t8_gtest_element_is_leaf.cxx index 51d32755ae..2422792767 100644 --- a/test/t8_forest/t8_gtest_element_is_leaf.cxx +++ b/test/t8_forest/t8_gtest_element_is_leaf.cxx @@ -145,7 +145,7 @@ TEST_P (element_is_leaf, element_is_leaf_adapt) t8_test_element_is_leaf_for_forest (forest_adapt); } -/* Define a lambda to beatify gtest output for tuples . +/* Define a lambda to beautify gtest output for tuples . * This will set the correct level and cmesh name as part of the test case name. */ auto pretty_print_level_and_cmesh_params = [] (const testing::TestParamInfo> &info) {