Feature vertex ghosts

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");
 

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,69 @@ 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");
+#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);
+  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);
+#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);
+      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 +715,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 +1499,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);
     }
 

src/t8_schemes/t8_default/t8_default_common/t8_default_common.hxx

@@ -273,6 +273,25 @@ class t8_default_scheme_common: public t8_crtp<TUnderlyingEclassScheme> {
     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

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_common<t8_default_scheme_
   void
   element_get_anchor (const t8_element_t *elem, int anchor[3]) const;
 
+  inline int
+  element_max_num_vertex_neighbors () const
+  {
+    return 4;
+  }
+
+  inline void
+  element_vertex_neighbors (const t8_element_t *element, const int vertex, int *num_neighbors, t8_element_t **neighbors,
+                            int *neigh_ivertices) const
+  {
+    p4est_quadrant_t *elem = (p4est_quadrant_t *) element;
+    p4est_qcoord_t len = P4EST_QUADRANT_LEN (elem->level);
+    *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 

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.

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
    */