From 2949bd657d3ce3ab539da316abcfe5fe14291da7 Mon Sep 17 00:00:00 2001 From: "Dreyer, Lukas" Date: Fri, 24 Jan 2025 13:40:53 +0100 Subject: [PATCH 1/5] Add vertex functionality to scheme interface --- .../t8_default_common/t8_default_common.hxx | 19 +++++++++++++ src/t8_schemes/t8_scheme.hxx | 27 +++++++++++++++++++ 2 files changed, 46 insertions(+) diff --git a/src/t8_schemes/t8_default/t8_default_common/t8_default_common.hxx b/src/t8_schemes/t8_default/t8_default_common/t8_default_common.hxx index 9d9bdb7d74..3f0d387789 100644 --- a/src/t8_schemes/t8_default/t8_default_common/t8_default_common.hxx +++ b/src/t8_schemes/t8_default/t8_default_common/t8_default_common.hxx @@ -273,6 +273,25 @@ class t8_default_scheme_common: public t8_crtp { return count_leaves_from_level (0, level, dim); } + inline int + element_max_num_vertex_neighbors () const + { + SC_ABORT ("Not implemented for this eclass\n"); + } + + inline void + element_vertex_neighbors (const t8_element_t *element, const int vertex, int *num_neighbors, t8_element_t **neighbors, + int *neigh_ivertices) const + { + SC_ABORT ("Not implemented for this eclass\n"); + } + + inline void + element_corner_descendant (const t8_element_t *element, int vertex, int level, t8_element_t *descendant) const + { + SC_ABORT ("Not implemented for this eclass\n"); + } + #if T8_ENABLE_DEBUG inline void element_debug_print (const t8_element_t *elem) const diff --git a/src/t8_schemes/t8_scheme.hxx b/src/t8_schemes/t8_scheme.hxx index dac78e4835..9afa900d09 100644 --- a/src/t8_schemes/t8_scheme.hxx +++ b/src/t8_schemes/t8_scheme.hxx @@ -837,6 +837,33 @@ class t8_scheme { eclass_schemes[tree_class]); }; + inline int + element_max_num_vertex_neighbors (const t8_eclass_t tree_class) const + { + return std::visit ([&] (auto &&scheme) { return scheme.element_max_num_vertex_neighbors (); }, + eclass_schemes[tree_class]); + } + + inline void + element_vertex_neighbors (const t8_eclass_t tree_class, const t8_element_t *element, const int vertex, + int *num_neighbors, t8_element_t **neighbors, int *neigh_ivertices) const + { + return std::visit ( + [&] (auto &&scheme) { + return scheme.element_vertex_neighbors (element, vertex, num_neighbors, neighbors, neigh_ivertices); + }, + eclass_schemes[tree_class]); + } + + inline void + element_corner_descendant (const t8_eclass_t tree_class, const t8_element_t *element, int vertex, int level, + t8_element_t *descendant) const + { + return std::visit ( + [&] (auto &&scheme) { return scheme.element_corner_descendant (element, vertex, level, descendant); }, + eclass_schemes[tree_class]); + } + /** Compute the coordinates of a given element vertex inside a reference tree * that is embedded into [0,1]^d (d = dimension). * \param [in] tree_class The eclass of the current tree. From 9a1c1d1d4a1c6d6e42114f10c43ed90fda8ab960 Mon Sep 17 00:00:00 2001 From: "Dreyer, Lukas" Date: Fri, 24 Jan 2025 13:44:15 +0100 Subject: [PATCH 2/5] Implement functionality for default quads --- .../t8_default_quad/t8_default_quad.hxx | 54 +++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/src/t8_schemes/t8_default/t8_default_quad/t8_default_quad.hxx b/src/t8_schemes/t8_default/t8_default_quad/t8_default_quad.hxx index 4908500395..7374c68cd9 100644 --- a/src/t8_schemes/t8_default/t8_default_quad/t8_default_quad.hxx +++ b/src/t8_schemes/t8_default/t8_default_quad/t8_default_quad.hxx @@ -546,6 +546,60 @@ class t8_default_scheme_quad: public t8_default_scheme_commonlevel); + *num_neighbors = 0; + const int dim = 2; + for (int icube = 0; icube < (1 << dim); icube++) { + int idim = 0; + p4est_qcoord_t shift = (vertex & 1 << idim) + (icube & 1 << idim); + shift >>= idim; + shift -= 1; + shift *= len; + ((p4est_quadrant_t *) neighbors[*num_neighbors])->x = elem->x + shift; + idim = 1; + shift = (vertex & 1 << idim) + (icube & 1 << idim); + shift >>= idim; + shift -= 1; + shift *= len; + ((p4est_quadrant_t *) neighbors[*num_neighbors])->y = elem->y + shift; + + ((p4est_quadrant_t *) neighbors[*num_neighbors])->level = elem->level; + + const int neigh_cube_vertex = (1 << dim) - 1 - icube; + neigh_ivertices[*num_neighbors] = neigh_cube_vertex; + if (element_is_equal (element, neighbors[*num_neighbors]) + || ((p4est_quadrant_t *) neighbors[*num_neighbors])->x < 0 + || ((p4est_quadrant_t *) neighbors[*num_neighbors])->x >= P4EST_ROOT_LEN + || ((p4est_quadrant_t *) neighbors[*num_neighbors])->y < 0 + || ((p4est_quadrant_t *) neighbors[*num_neighbors])->y >= P4EST_ROOT_LEN) { + continue; + } + ++(*num_neighbors); + } + } + + inline void + element_corner_descendant (const t8_element_t *element, int vertex, int level, t8_element_t *descendant) const + { + p4est_quadrant_t *elem = (p4est_quadrant_t *) element; + p4est_quadrant_t *desc = (p4est_quadrant_t *) descendant; + p4est_qcoord_t coord_offset = P4EST_QUADRANT_LEN (elem->level) - P4EST_QUADRANT_LEN (level); + desc->x = elem->x + coord_offset * ((vertex & 1 << 0) >> 0); + desc->y = elem->y + coord_offset * ((vertex & 1 << 1) >> 1); + desc->level = level; + } + /** Compute the integer coordinates of a given element vertex. * The default scheme implements the Morton type SFCs. In these SFCs the * elements are positioned in a cube [0,1]^(dL) with dimension d (=0,1,2,3) and From d6d69dff9bb8de49c0fc3d388eae91b8baa01c46 Mon Sep 17 00:00:00 2001 From: "Dreyer, Lukas" Date: Fri, 24 Jan 2025 13:45:12 +0100 Subject: [PATCH 3/5] enable vertex ghosts --- src/t8_forest/t8_forest.cxx | 2 +- src/t8_forest/t8_forest_ghost.cxx | 80 +++++++++++++++++++++++++++++-- 2 files changed, 76 insertions(+), 6 deletions(-) diff --git a/src/t8_forest/t8_forest.cxx b/src/t8_forest/t8_forest.cxx index 3b1f8bac36..8a167b854e 100644 --- a/src/t8_forest/t8_forest.cxx +++ b/src/t8_forest/t8_forest.cxx @@ -2899,7 +2899,7 @@ t8_forest_set_ghost_ext (t8_forest_t forest, int do_ghost, t8_ghost_type_t ghost { T8_ASSERT (t8_forest_is_initialized (forest)); /* We currently only support face ghosts */ - SC_CHECK_ABORT (do_ghost == 0 || ghost_type == T8_GHOST_FACES, + SC_CHECK_ABORT (do_ghost == 0 || ghost_type == T8_GHOST_FACES || ghost_type == T8_GHOST_VERTICES, "Ghost neighbors other than face-neighbors are not supported.\n"); SC_CHECK_ABORT (1 <= ghost_version && ghost_version <= 3, "Invalid choice for ghost version. Choose 1, 2, or 3.\n"); diff --git a/src/t8_forest/t8_forest_ghost.cxx b/src/t8_forest/t8_forest_ghost.cxx index 7cbbde1edc..c5078d9f2c 100644 --- a/src/t8_forest/t8_forest_ghost.cxx +++ b/src/t8_forest/t8_forest_ghost.cxx @@ -173,7 +173,7 @@ t8_forest_ghost_init (t8_forest_ghost_t *pghost, t8_ghost_type_t ghost_type) t8_forest_ghost_t ghost; /* We currently only support face-neighbor ghosts */ - T8_ASSERT (ghost_type == T8_GHOST_FACES); + T8_ASSERT (ghost_type == T8_GHOST_FACES || ghost_type == T8_GHOST_VERTICES); /* Allocate memory for ghost */ ghost = *pghost = T8_ALLOC_ZERO (t8_forest_ghost_struct_t, 1); @@ -502,8 +502,66 @@ typedef struct } t8_forest_ghost_boundary_data_t; static int -t8_forest_ghost_search_boundary (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, - const int is_leaf, const t8_element_array_t *leaves, const t8_locidx_t tree_leaf_index) +t8_forest_ghost_search_vertex_boundary (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, + const int is_leaf, const t8_element_array_t *leaves, + const t8_locidx_t tree_leaf_index) +{ + T8_ASSERT (t8_forest_get_num_global_trees (forest) == 1); + T8_ASSERT (ltreeid == 0); + + if (!is_leaf) { + /** TODO: Cut off search if all neighbors are owned by ourselves. + * Use bounds to efficiently check */ + return 1; + } + + /** + * Determine all vertex neighbors and determine their owner. + */ + + t8_eclass_t eclass = t8_forest_get_eclass (forest, ltreeid); + const t8_scheme *scheme = t8_forest_get_scheme (forest); + + t8_debugf ("enter t8_forest_ghost_search_vertex_boundary for :\n"); + scheme->element_debug_print (eclass, element); + + const int max_neighs = scheme->element_max_num_vertex_neighbors (eclass); + t8_element_t **vertex_neighbors = T8_ALLOC (t8_element_t *, max_neighs); + scheme->element_new (eclass, max_neighs, vertex_neighbors); + int *neigh_ivertices = T8_ALLOC (int, max_neighs); + t8_element_t *vertex_descendant; + scheme->element_new (eclass, 1, &vertex_descendant); + for (int ivertex = 0; ivertex < scheme->element_get_num_corners (eclass, element); ivertex++) { + int num_neighbors; + t8_debugf ("determine vertex neighbors around vertex %i\n", ivertex); + scheme->element_vertex_neighbors (eclass, element, ivertex, &num_neighbors, vertex_neighbors, neigh_ivertices); + t8_debugf ("Found %i possible neighbors\n", num_neighbors); + for (int ineigh = 0; ineigh < num_neighbors; ineigh++) { + t8_debugf ("neigh %i, neigh_ivertex%i\n", ineigh, neigh_ivertices[ineigh]); + scheme->element_corner_descendant (eclass, vertex_neighbors[ineigh], neigh_ivertices[ineigh], + scheme->get_maxlevel (eclass), vertex_neighbors[ineigh]); + t8_debugf ("corner_descendant of neigh %i:\n", ineigh); + scheme->element_debug_print (eclass, vertex_neighbors[ineigh]); + const int remote_rank + = t8_forest_element_find_owner_ext (forest, t8_forest_global_tree_id (forest, ltreeid), + vertex_neighbors[ineigh], eclass, 0, forest->mpisize - 1, 0, 1); + t8_debugf ("add remote for rank %i\n", remote_rank); + if (remote_rank != forest->mpirank) { + t8_ghost_add_remote (forest, forest->ghosts, remote_rank, ltreeid, element, tree_leaf_index); + } + } + } + scheme->element_destroy (eclass, 1, &vertex_descendant); + scheme->element_destroy (eclass, max_neighs, vertex_neighbors); + T8_FREE (neigh_ivertices); + T8_FREE (vertex_neighbors); + return 0; +} + +static int +t8_forest_ghost_search_face_boundary (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, + const int is_leaf, const t8_element_array_t *leaves, + const t8_locidx_t tree_leaf_index) { t8_forest_ghost_boundary_data_t *data = (t8_forest_ghost_boundary_data_t *) t8_forest_get_user_data (forest); int num_faces, iface, faces_totally_owned, level; @@ -654,7 +712,17 @@ t8_forest_ghost_fill_remote_v3 (t8_forest_t forest) /* Set the user data for the search routine */ t8_forest_set_user_data (forest, &data); /* Loop over the trees of the forest */ - t8_forest_search (forest, t8_forest_ghost_search_boundary, NULL, NULL); + + switch (forest->ghost_type) { + case T8_GHOST_FACES: + t8_forest_search (forest, t8_forest_ghost_search_face_boundary, NULL, NULL); + break; + case T8_GHOST_VERTICES: + t8_forest_search (forest, t8_forest_ghost_search_vertex_boundary, NULL, NULL); + break; + default: + SC_ABORT ("Edge ghosts not yet implemented!"); + } /* Reset the user data from before search */ t8_forest_set_user_data (forest, store_user_data); @@ -1428,17 +1496,19 @@ t8_forest_ghost_create_ext (t8_forest_t forest, int unbalanced_version) return; } /* Currently we only support face ghosts */ - T8_ASSERT (forest->ghost_type == T8_GHOST_FACES); + T8_ASSERT (forest->ghost_type == T8_GHOST_FACES || forest->ghost_type == T8_GHOST_VERTICES); /* Initialize the ghost structure */ t8_forest_ghost_init (&forest->ghosts, forest->ghost_type); ghost = forest->ghosts; if (unbalanced_version == -1) { + t8_debugf ("remote v3\n"); t8_forest_ghost_fill_remote_v3 (forest); } else { /* Construct the remote elements and processes. */ + t8_debugf ("remote old\n"); t8_forest_ghost_fill_remote (forest, ghost, unbalanced_version != 0); } From 5d85bb5bf54cc83548513f3e2f06b951b457ec4e Mon Sep 17 00:00:00 2001 From: "Dreyer, Lukas" Date: Fri, 24 Jan 2025 13:46:11 +0100 Subject: [PATCH 4/5] adapt tutorial to vertex ghost --- .../t8_step4_partition_balance_ghost.cxx | 62 ++++--------------- 1 file changed, 12 insertions(+), 50 deletions(-) diff --git a/tutorials/general/t8_step4_partition_balance_ghost.cxx b/tutorials/general/t8_step4_partition_balance_ghost.cxx index 80d08ac7f4..509799dc02 100644 --- a/tutorials/general/t8_step4_partition_balance_ghost.cxx +++ b/tutorials/general/t8_step4_partition_balance_ghost.cxx @@ -149,40 +149,19 @@ t8_step4_partition_ghost (t8_forest_t forest) * We currently support ghost mode T8_GHOST_FACES that creates face neighbor ghost elements * and will in future also support other modes for edge/vertex neighbor ghost elements. */ - t8_forest_set_ghost (new_forest, 1, T8_GHOST_FACES); + t8_forest_set_ghost (new_forest, 1, T8_GHOST_VERTICES); /* Commit the forest, this step will perform the partitioning and ghost layer creation. */ t8_forest_commit (new_forest); return new_forest; } -/* In this function we adapt a forest as in step3 and balance it. - * In our main program the input forest is already adapted and then the resulting twice adapted forest will be unbalanced. - */ -static t8_forest_t -t8_step4_balance (t8_forest_t forest) +int +t8_refine_first_child (t8_forest_t forest, t8_forest_t forest_from, t8_locidx_t which_tree, t8_eclass_t tree_class, + t8_locidx_t lelement_id, const t8_scheme *scheme, const int is_family, const int num_elements, + t8_element_t *elements[]) { - t8_forest_t balanced_forest; - /* Adapt the input forest. */ - t8_forest_t unbalanced_forest = t8_step3_adapt_forest (forest); - - /* Output to vtk. */ - t8_forest_write_vtk (unbalanced_forest, "t8_step4_unbalanced_forest"); - t8_global_productionf (" [step4] Wrote unbalanced forest to vtu files: %s*\n", "t8_step4_unbalanced_forest"); - - /* Initialize new forest. */ - t8_forest_init (&balanced_forest); - /* Specify that this forest should result from balancing unbalanced_forest. - * The last argument is the flag 'no_repartition'. - * Since balancing will refine elements, the load-balance will be broken afterwards. - * Setting this flag to false (no_repartition = false -> yes repartition) will repartition - * the forest after balance, such that every process has the same number of elements afterwards. - */ - t8_forest_set_balance (balanced_forest, unbalanced_forest, 0); - /* Commit the forest. */ - t8_forest_commit (balanced_forest); - - return balanced_forest; + return !scheme->element_get_child_id (tree_class, elements[0]); } int @@ -196,7 +175,6 @@ t8_step4_main (int argc, char **argv) const char *prefix_uniform = "t8_step4_uniform_forest"; const char *prefix_adapt = "t8_step4_adapted_forest"; const char *prefix_partition_ghost = "t8_step4_partitioned_ghost_forest"; - const char *prefix_balance = "t8_step4_balanced_forest"; /* The uniform refinement level of the forest. */ const int level = 3; @@ -208,7 +186,7 @@ t8_step4_main (int argc, char **argv) /* Initialize the sc library, has to happen before we initialize t8code. */ sc_init (sc_MPI_COMM_WORLD, 1, 1, NULL, SC_LP_ESSENTIAL); /* Initialize t8code with log level SC_LP_PRODUCTION. See sc.h for more info on the log levels. */ - t8_init (SC_LP_PRODUCTION); + t8_init (SC_LP_DEBUG); /* Print a message on the root process. */ t8_global_productionf (" [step4] \n"); @@ -230,7 +208,7 @@ t8_step4_main (int argc, char **argv) t8_global_productionf (" [step4] Creating an adapted forest as in step3.\n"); t8_global_productionf (" [step4] \n"); /* Build a cube cmesh with tet, hex, and prism trees. */ - cmesh = t8_cmesh_new_hypercube_hybrid (comm, 0, 0); + cmesh = t8_cmesh_new_hypercube (T8_ECLASS_QUAD, comm, 0, 0, 0); t8_global_productionf (" [step4] Created coarse mesh.\n"); forest = t8_forest_new_uniform (cmesh, t8_scheme_new_default (), level, 0, comm); @@ -242,15 +220,13 @@ t8_step4_main (int argc, char **argv) t8_global_productionf (" [step4] Wrote uniform level %i forest to vtu files: %s*\n", level, prefix_uniform); /* - * Adapt the forest. + * Adapt. */ - /* Adapt the forest. We can reuse the forest variable, since the new adapted - * forest will take ownership of the old forest and destroy it. - * Note that the adapted forest is a new forest, though. */ - forest = t8_step3_adapt_forest (forest); + forest = t8_forest_new_adapt (forest, t8_refine_first_child, 0, 0, NULL); + forest = t8_forest_new_adapt (forest, t8_refine_first_child, 0, 0, NULL); - /* Print information of our new forest. */ + /* Print information of the forest. */ t8_step3_print_forest_information (forest); /* Write forest to vtu files. */ @@ -270,20 +246,6 @@ t8_step4_main (int argc, char **argv) /* Write forest to vtu files. */ t8_forest_write_vtk_ext (forest, prefix_partition_ghost, 1, 1, 1, 1, 1, 0, 1, 0, NULL); - /* - * Balance - */ - - t8_global_productionf (" [step4] \n"); - t8_global_productionf (" [step4] Creating an unbalanced forest and balancing it.\n"); - t8_global_productionf (" [step4] \n"); - forest = t8_step4_balance (forest); - t8_global_productionf (" [step4] Balanced forest.\n"); - t8_step3_print_forest_information (forest); - /* Write forest to vtu files. */ - t8_forest_write_vtk (forest, prefix_balance); - t8_global_productionf (" [step4] Wrote balanced forest to vtu files: %s*\n", prefix_balance); - /* * clean-up */ From 4218bec56650e3bacf1a4f61955ced63f74b4d90 Mon Sep 17 00:00:00 2001 From: "Dreyer, Lukas" Date: Fri, 24 Jan 2025 13:54:21 +0100 Subject: [PATCH 5/5] move debug print into T8_ENABLE_DEBUG --- src/t8_forest/t8_forest_ghost.cxx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/t8_forest/t8_forest_ghost.cxx b/src/t8_forest/t8_forest_ghost.cxx index c5078d9f2c..818ff28591 100644 --- a/src/t8_forest/t8_forest_ghost.cxx +++ b/src/t8_forest/t8_forest_ghost.cxx @@ -523,8 +523,9 @@ t8_forest_ghost_search_vertex_boundary (t8_forest_t forest, t8_locidx_t ltreeid, const t8_scheme *scheme = t8_forest_get_scheme (forest); t8_debugf ("enter t8_forest_ghost_search_vertex_boundary for :\n"); +#ifdef T8_ENABLE_DEBUG scheme->element_debug_print (eclass, element); - +#endif const int max_neighs = scheme->element_max_num_vertex_neighbors (eclass); t8_element_t **vertex_neighbors = T8_ALLOC (t8_element_t *, max_neighs); scheme->element_new (eclass, max_neighs, vertex_neighbors); @@ -541,7 +542,9 @@ t8_forest_ghost_search_vertex_boundary (t8_forest_t forest, t8_locidx_t ltreeid, scheme->element_corner_descendant (eclass, vertex_neighbors[ineigh], neigh_ivertices[ineigh], scheme->get_maxlevel (eclass), vertex_neighbors[ineigh]); t8_debugf ("corner_descendant of neigh %i:\n", ineigh); +#ifdef T8_ENABLE_DEBUG scheme->element_debug_print (eclass, vertex_neighbors[ineigh]); +#endif const int remote_rank = t8_forest_element_find_owner_ext (forest, t8_forest_global_tree_id (forest, ltreeid), vertex_neighbors[ineigh], eclass, 0, forest->mpisize - 1, 0, 1);