Introduce populate functionality according to cmesh borders

src/t8_forest/t8_forest.c

@@ -60,6 +60,7 @@ t8_forest_init (t8_forest_t *pforest)
   forest->maxlevel_existing = -1;
   forest->stats_computed = 0;
   forest->incomplete_trees = -1;
+  forest->set_initial_partition_according_to_cmesh = 0;
 }
 
 int
@@ -139,6 +140,13 @@ t8_forest_set_level (t8_forest_t forest, int level)
   forest->set_level = level;
 }
 
+void
+t8_forest_set_uniform_partition_from_cmesh (t8_forest_t forest, t8_cmesh_t cmesh, sc_MPI_Comm comm)
+{
+  t8_forest_set_cmesh (forest, cmesh, comm);
+  forest->set_initial_partition_according_to_cmesh = 1;
+}
+
 void
 t8_forest_set_copy (t8_forest_t forest, const t8_forest_t set_from)
 {
@@ -381,11 +389,10 @@ t8_forest_refines_irregular (t8_forest_t forest)
  * \param[in] forest  The forest to populate
 */
 static void
-t8_forest_populate_irregular (t8_forest_t forest)
+t8_forest_populate_irregular (t8_forest_t forest, int repartition)
 {
   t8_forest_t forest_zero;
   t8_forest_t forest_tmp;
-  t8_forest_t forest_tmp_partition;
   t8_cmesh_ref (forest->cmesh);
   t8_scheme_cxx_ref (forest->scheme_cxx);
   /* We start with a level 0 uniform refinement */
@@ -400,16 +407,15 @@ t8_forest_populate_irregular (t8_forest_t forest)
     t8_forest_init (&forest_tmp);
     t8_forest_set_level (forest_tmp, i);
     t8_forest_set_adapt (forest_tmp, forest_zero, t8_forest_refine_everything, 0);
+    if (!forest->set_initial_partition_according_to_cmesh) {
+      t8_forest_set_partition (forest_tmp, forest_zero, 0);
+    }
     t8_forest_commit (forest_tmp);
-    /* Partition the forest to even the load */
-    t8_forest_init (&forest_tmp_partition);
-    t8_forest_set_partition (forest_tmp_partition, forest_tmp, 0);
-    t8_forest_commit (forest_tmp_partition);
-    forest_zero = forest_tmp_partition;
+    forest_zero = forest_tmp;
   }
   /* Copy all elements over to the original forest. */
   t8_forest_copy_trees (forest, forest_zero, 1);
-  t8_forest_unref (&forest_tmp_partition);
+  t8_forest_unref (&forest_tmp); /* same as forest zero */
 }
 
 void
@@ -450,12 +456,15 @@ t8_forest_commit (t8_forest_t forest)
     mpiret = sc_MPI_Comm_rank (forest->mpicomm, &forest->mpirank);
     SC_CHECK_MPI (mpiret);
     /* Compute the maximum allowed refinement level */
+    t8_debugf ("compute maxlevel of forest:\n");
     t8_forest_compute_maxlevel (forest);
+    t8_debugf ("computed maxlevel: %i\n", forest->maxlevel);
     T8_ASSERT (forest->set_level <= forest->maxlevel);
     /* populate a new forest with tree and quadrant objects */
     if (t8_forest_refines_irregular (forest) && forest->set_level > 0) {
       /* On root level we will also use the normal algorithm */
-      t8_forest_populate_irregular (forest);
+      t8_forest_populate_irregular (forest, 1);
+      //      t8_forest_populate_irregular (forest, !forest->set_initial_partition_according_to_cmesh);
     }
     else {
       t8_forest_populate (forest);
@@ -1433,6 +1442,34 @@ t8_forest_new_uniform (t8_cmesh_t cmesh, t8_scheme_cxx_t *scheme, const int leve
   return forest;
 }
 
+t8_forest_t
+t8_forest_new_uniform_with_cmesh_partition (t8_cmesh_t cmesh, t8_scheme_cxx_t *scheme, const int level,
+                                            const int do_face_ghost, sc_MPI_Comm comm)
+{
+  t8_forest_t forest;
+
+  T8_ASSERT (t8_cmesh_is_committed (cmesh));
+  T8_ASSERT (scheme != NULL);
+  T8_ASSERT (0 <= level);
+
+  /* Initialize the forest */
+  t8_forest_init (&forest);
+  /* Set the cmesh, scheme and level */
+  t8_forest_set_uniform_partition_from_cmesh (forest, cmesh, comm);
+  t8_forest_set_scheme (forest, scheme);
+  t8_forest_set_level (forest, level);
+  if (do_face_ghost) {
+    t8_forest_set_ghost (forest, 1, T8_GHOST_FACES);
+  }
+
+  /* commit the forest */
+  t8_forest_commit (forest);
+  t8_global_productionf ("Constructed uniform forest with %lli global elements.\n",
+                         (long long) forest->global_num_elements);
+
+  return forest;
+}
+
 t8_forest_t
 t8_forest_new_adapt (t8_forest_t forest_from, t8_forest_adapt_t adapt_fn, int recursive, int do_face_ghost,
                      void *user_data)

src/t8_forest/t8_forest_cxx.cxx

@@ -192,6 +192,9 @@ t8_forest_compute_maxlevel (t8_forest_t forest)
       }
     }
   }
+  int global_maxlevel;
+  sc_MPI_Allreduce (&forest->maxlevel, &global_maxlevel, 1, sc_MPI_INT, sc_MPI_MAX, forest->mpicomm);
+  forest->maxlevel = global_maxlevel;
   T8_ASSERT (forest->maxlevel >= 0);
   t8_debugf ("Computed maxlevel %i\n", forest->maxlevel);
 }
@@ -1153,86 +1156,155 @@ t8_forest_compute_desc (t8_forest_t forest)
   }
 }
 
+void
+t8_forest_populate_according_to_cmesh (t8_forest_t forest)
+{
+  SC_CHECK_ABORT (!forest->cmesh->first_tree_shared, "Cmesh needs to be partitioned without shared elements \n");
+  forest->first_local_tree = t8_cmesh_get_first_treeid (forest->cmesh);
+  int num_local_trees = t8_cmesh_get_num_local_trees (forest->cmesh);
+  forest->last_local_tree = forest->first_local_tree + num_local_trees - 1;
+
+  forest->global_num_elements = forest->local_num_elements = 0;
+  int is_empty = (num_local_trees == 0);
+  t8_linearidx_t count_elements;
+  /* create only the non-empty tree objects */
+  if (is_empty) {
+    t8_debugf ("is empty\n");
+    /* This processor is empty
+     * we still set the tree array to store 0 as the number of trees here */
+    forest->trees = sc_array_new (sizeof (t8_tree_struct_t));
+    count_elements = 0;
+    /* Set the first local tree larger than the last local tree to
+     * indicate empty forest */
+    forest->first_local_tree = forest->last_local_tree + 1;
+  }
+  else {
+    t8_debugf ("has %i trees\n", num_local_trees);
+    /* for each tree, allocate elements */
+    T8_ASSERT (num_local_trees == forest->last_local_tree - forest->first_local_tree + 1);
+
+    forest->trees = sc_array_new_count (sizeof (t8_tree_struct_t), num_local_trees);
+    t8_gloidx_t first_ctree = t8_cmesh_get_first_treeid (forest->cmesh);
+    count_elements = 0; /** apparently this cannot be in the for loop initialisation. */
+    for (t8_gloidx_t jt = forest->first_local_tree; jt <= forest->last_local_tree; jt++) {
+      t8_tree_t tree = (t8_tree_t) t8_sc_array_index_locidx (forest->trees, jt - forest->first_local_tree);
+      t8_eclass_t tree_class = tree->eclass = t8_cmesh_get_tree_class (forest->cmesh, jt - first_ctree);
+      tree->elements_offset = count_elements;
+      t8_eclass_scheme_c *eclass_scheme = forest->scheme_cxx->eclass_schemes[tree_class];
+      T8_ASSERT (eclass_scheme != NULL);
+      t8_element_array_t *telements = &tree->elements;
+      /* calculate first and last element on this tree */
+      t8_locidx_t start = 0;
+      t8_locidx_t end = eclass_scheme->t8_element_count_leaves_from_root (forest->set_level);
+      t8_debugf ("end: %i\n", end);
+      /* Allocate elements for this processor. */
+      t8_element_array_init_size (telements, eclass_scheme, end - start);
+      t8_element_t *element = t8_element_array_index_locidx_mutable (telements, 0);
+      eclass_scheme->t8_element_set_linear_id (element, forest->set_level, start);
+      count_elements++;
+      for (t8_locidx_t et = start + 1; et < end; et++, count_elements++) {
+        t8_element_t *element_succ = t8_element_array_index_locidx_mutable (telements, et - start);
+        T8_ASSERT (eclass_scheme->t8_element_level (element) == forest->set_level);
+        eclass_scheme->t8_element_successor (element, element_succ);
+        element = element_succ;
+      }
+      t8_debugf ("count_elements in loop:%li\n", count_elements);
+    }
+    t8_debugf ("count_elements:%li\n", count_elements);
+  }
+  forest->local_num_elements = count_elements;
+  /* TODO: if no tree has pyramid type we can optimize this to global_num_elements = global_num_trees * 2^(dim*level) */
+  t8_forest_comm_global_num_elements (forest);
+  /* TODO: figure out global_first_position, global_first_quadrant without comm */
+}
+
 /* Create the elements on this process given a uniform partition of the coarse mesh. */
 void
 t8_forest_populate (t8_forest_t forest)
 {
-  t8_gloidx_t child_in_tree_begin;
-  t8_gloidx_t child_in_tree_end;
-  t8_locidx_t count_elements;
-  t8_locidx_t num_tree_elements;
-  t8_locidx_t num_local_trees;
-  t8_gloidx_t jt, first_ctree;
-  t8_gloidx_t start, end, et;
-  t8_tree_t tree;
-  t8_element_t *element, *element_succ;
-  t8_element_array_t *telements;
-  t8_eclass_t tree_class;
-  t8_eclass_scheme_c *eclass_scheme;
-  t8_gloidx_t cmesh_first_tree, cmesh_last_tree;
-  int is_empty;
-
   SC_CHECK_ABORT (forest->set_level <= forest->maxlevel, "Given refinement level exceeds the maximum.\n");
-  /* TODO: create trees and quadrants according to uniform refinement */
-  t8_cmesh_uniform_bounds (forest->cmesh, forest->set_level, forest->scheme_cxx, &forest->first_local_tree,
-                           &child_in_tree_begin, &forest->last_local_tree, &child_in_tree_end, NULL);
 
-  /* True if the forest has no elements */
-  is_empty = forest->first_local_tree > forest->last_local_tree
-             || (forest->first_local_tree == forest->last_local_tree && child_in_tree_begin >= child_in_tree_end);
+  t8_gloidx_t cmesh_first_tree = t8_cmesh_get_first_treeid (forest->cmesh);
+  t8_gloidx_t cmesh_last_tree = cmesh_first_tree + t8_cmesh_get_num_local_trees (forest->cmesh) - 1;
+  t8_gloidx_t
+    child_in_tree_begin; /** does not get filled when the cmesh does not contain any shared elements between processes */
+  t8_gloidx_t child_in_tree_end; /** see above*/
 
-  cmesh_first_tree = t8_cmesh_get_first_treeid (forest->cmesh);
-  cmesh_last_tree = cmesh_first_tree + t8_cmesh_get_num_local_trees (forest->cmesh) - 1;
+  if (forest->set_initial_partition_according_to_cmesh) {
+    SC_CHECK_ABORT (!forest->cmesh->first_tree_shared, "Cmesh needs to be partitioned without shared elements for a "
+                                                       "forest that is partitioned in the same way as the cmesh \n");
+    forest->first_local_tree = cmesh_first_tree;
+    forest->last_local_tree = cmesh_last_tree;
+  }
+  else {
+    t8_cmesh_uniform_bounds (forest->cmesh, forest->set_level, forest->scheme_cxx, &forest->first_local_tree,
+                             &child_in_tree_begin, &forest->last_local_tree, &child_in_tree_end, NULL);
+  }
 
-  if (!is_empty) {
+  /* True if the forest has no elements */
+  int is_empty
+    = forest->first_local_tree > forest->last_local_tree
+      || (!forest->set_initial_partition_according_to_cmesh && forest->first_local_tree == forest->last_local_tree
+          && child_in_tree_begin >= child_in_tree_end);
+
+  if (!is_empty && !forest->set_initial_partition_according_to_cmesh) {
     SC_CHECK_ABORT (forest->first_local_tree >= cmesh_first_tree && forest->last_local_tree <= cmesh_last_tree,
-                    "cmesh partition does not match the planned forest partition");
+                    "cmesh partition does not match the planned uniform equally distributed forest partition");
   }
 
   forest->global_num_elements = forest->local_num_elements = 0;
   /* create only the non-empty tree objects */
   if (is_empty) {
+    forest->local_num_elements = 0;
     /* This processor is empty
      * we still set the tree array to store 0 as the number of trees here */
     forest->trees = sc_array_new (sizeof (t8_tree_struct_t));
-    count_elements = 0;
     /* Set the first local tree larger than the last local tree to
-     * indicate empty forest */
+     * indicate empty forest, this is sometimes used instead of the more accurate check that local_num_element = 0. TODO: check that line can be removed */
     forest->first_local_tree = forest->last_local_tree + 1;
   }
   else {
     /* for each tree, allocate elements */
-    num_local_trees = forest->last_local_tree - forest->first_local_tree + 1;
+    t8_locidx_t num_local_trees = forest->last_local_tree - forest->first_local_tree + 1;
     forest->trees = sc_array_new_count (sizeof (t8_tree_struct_t), num_local_trees);
-    first_ctree = t8_cmesh_get_first_treeid (forest->cmesh);
-    for (jt = forest->first_local_tree, count_elements = 0; jt <= forest->last_local_tree; jt++) {
-      tree = (t8_tree_t) t8_sc_array_index_locidx (forest->trees, jt - forest->first_local_tree);
-      tree_class = tree->eclass = t8_cmesh_get_tree_class (forest->cmesh, jt - first_ctree);
+    t8_gloidx_t first_ctree = t8_cmesh_get_first_treeid (forest->cmesh);
+    t8_locidx_t count_elements = 0;
+    for (t8_gloidx_t jt = forest->first_local_tree; jt <= forest->last_local_tree; jt++) {
+      t8_tree_t tree = (t8_tree_t) t8_sc_array_index_locidx (forest->trees, jt - forest->first_local_tree);
+      t8_eclass_t tree_class = tree->eclass = t8_cmesh_get_tree_class (forest->cmesh, jt - first_ctree);
       tree->elements_offset = count_elements;
-      eclass_scheme = forest->scheme_cxx->eclass_schemes[tree_class];
+      t8_eclass_scheme_c *eclass_scheme = forest->scheme_cxx->eclass_schemes[tree_class];
       T8_ASSERT (eclass_scheme != NULL);
-      telements = &tree->elements;
+      t8_element_array_t *telements = &tree->elements;
       /* calculate first and last element on this tree */
-      start = (jt == forest->first_local_tree) ? child_in_tree_begin : 0;
-      end = (jt == forest->last_local_tree) ? child_in_tree_end
-                                            : eclass_scheme->t8_element_count_leaves_from_root (forest->set_level);
+      t8_locidx_t start, end, num_tree_elements;
+      if (forest->set_initial_partition_according_to_cmesh) {
+        start = 0;
+        end = eclass_scheme->t8_element_count_leaves_from_root (forest->set_level);
+      }
+      else {
+        start = (jt == forest->first_local_tree) ? child_in_tree_begin : 0;
+        end = (jt == forest->last_local_tree) ? child_in_tree_end
+                                              : eclass_scheme->t8_element_count_leaves_from_root (forest->set_level);
+      }
       num_tree_elements = end - start;
+      t8_debugf ("num_tree_elements: %i \n", num_tree_elements);
       T8_ASSERT (num_tree_elements > 0);
       /* Allocate elements for this processor. */
       t8_element_array_init_size (telements, eclass_scheme, num_tree_elements);
-      element = t8_element_array_index_locidx_mutable (telements, 0);
+      t8_element_t *element = t8_element_array_index_locidx_mutable (telements, 0);
       eclass_scheme->t8_element_set_linear_id (element, forest->set_level, start);
       count_elements++;
-      for (et = start + 1; et < end; et++, count_elements++) {
-        element_succ = t8_element_array_index_locidx_mutable (telements, et - start);
+      for (t8_locidx_t et = start + 1; et < end; et++, count_elements++) {
+        t8_element_t *element_succ = t8_element_array_index_locidx_mutable (telements, et - start);
         T8_ASSERT (eclass_scheme->t8_element_level (element) == forest->set_level);
         eclass_scheme->t8_element_successor (element, element_succ);
         /* TODO: process elements here */
         element = element_succ;
       }
     }
+    forest->local_num_elements = count_elements;
   }
-  forest->local_num_elements = count_elements;
   /* TODO: if no tree has pyramid type we can optimize this to global_num_elements = global_num_trees * 2^(dim*level) */
   t8_forest_comm_global_num_elements (forest);
   /* TODO: figure out global_first_position, global_first_quadrant without comm */

src/t8_forest/t8_forest_general.h

@@ -213,6 +213,9 @@ t8_forest_set_scheme (t8_forest_t forest, t8_scheme_cxx_t *scheme);
 void
 t8_forest_set_level (t8_forest_t forest, int level);
 
+void
+t8_forest_set_uniform_partition_from_cmesh (t8_forest_t forest, t8_cmesh_t cmesh, sc_MPI_Comm comm);
+
 /** Set a forest as source for copying on committing.
  * By default, the forest takes ownership of the source \b from such that it will
  * be destroyed on calling \ref t8_forest_commit.  To keep ownership of \b
@@ -840,6 +843,10 @@ t8_forest_t
 t8_forest_new_uniform (t8_cmesh_t cmesh, t8_scheme_cxx_t *scheme, const int level, const int do_face_ghost,
                        sc_MPI_Comm comm);
 
+t8_forest_t
+t8_forest_new_uniform_with_cmesh_partition (t8_cmesh_t cmesh, t8_scheme_cxx_t *scheme, const int level,
+                                            const int do_face_ghost, sc_MPI_Comm comm);
+
 /** Build a adapted forest from another forest.
  * \param [in]    forest_from The forest to refine
  * \param [in]    adapt_fn    Adapt function to use

src/t8_forest/t8_forest_private.h

@@ -68,6 +68,9 @@ t8_forest_compute_desc (t8_forest_t forest);
 void
 t8_forest_populate (t8_forest_t forest);
 
+void
+t8_forest_populate_according_to_cmesh (t8_forest_t forest);
+
 /** Return the eclass scheme of a given element class associated to a forest.
  * This function does not check whether the given forest is committed, use with
  * caution and only if you are sure that the eclass_scheme was set.

src/t8_forest/t8_forest_types.h

@@ -69,7 +69,9 @@ typedef struct t8_forest
 {
   t8_refcount_t rc; /**< Reference counter. */
 
-  int set_level;          /**< Level to use in new construction. */
+  int set_level; /**< Level to use in new construction. */
+  int
+    set_initial_partition_according_to_cmesh; /** Decide whether to partition the elements according to a given cmesh without shared elements or equally */
   int set_for_coarsening; /**< Change partition to allow
                                                      for one round of coarsening */
 

test/Makefile.am

@@ -71,6 +71,7 @@ t8code_googletest_programs = \
   test/t8_schemes/t8_gtest_face_descendant \
   test/t8_geometry/t8_gtest_point_inside \
   test/t8_forest/t8_gtest_user_data \
+  test/t8_forest/t8_gtest_forest_new_from_part_cmesh \
   test/t8_forest/t8_gtest_transform \
   test/t8_forest/t8_gtest_ghost_exchange \
   test/t8_forest/t8_gtest_ghost_delete \
@@ -270,6 +271,10 @@ test_t8_forest_t8_gtest_user_data_SOURCES = \
   test/t8_gtest_main.cxx \
   test/t8_forest/t8_gtest_user_data.cxx
 
+test_t8_forest_t8_gtest_forest_new_from_part_cmesh_SOURCES = \
+  test/t8_gtest_main.cxx \
+  test/t8_forest/t8_gtest_forest_new_from_part_cmesh.cxx
+
 test_t8_forest_t8_gtest_transform_SOURCES = \
   test/t8_gtest_main.cxx \
   test/t8_forest/t8_gtest_transform.cxx
@@ -526,6 +531,10 @@ test_t8_forest_t8_gtest_user_data_LDADD = $(t8_gtest_target_ld_add)
 test_t8_forest_t8_gtest_user_data_LDFLAGS = $(t8_gtest_target_ld_flags)
 test_t8_forest_t8_gtest_user_data_CPPFLAGS = $(t8_gtest_target_cpp_flags)
 
+test_t8_forest_t8_gtest_forest_new_from_part_cmesh_LDADD = $(t8_gtest_target_ld_add)
+test_t8_forest_t8_gtest_forest_new_from_part_cmesh_LDFLAGS = $(t8_gtest_target_ld_flags)
+test_t8_forest_t8_gtest_forest_new_from_part_cmesh_CPPFLAGS = $(t8_gtest_target_cpp_flags)
+
 test_t8_forest_t8_gtest_transform_LDADD = $(t8_gtest_target_ld_add)
 test_t8_forest_t8_gtest_transform_LDFLAGS = $(t8_gtest_target_ld_flags)
 test_t8_forest_t8_gtest_transform_CPPFLAGS = $(t8_gtest_target_cpp_flags)
@@ -648,6 +657,7 @@ test_t8_forest_t8_gtest_forest_face_normal_CPPFLAGS += $(t8_gtest_target_mpi_cpp
 test_t8_schemes_t8_gtest_face_descendant_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_geometry_t8_gtest_point_inside_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_forest_t8_gtest_user_data_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
+test_t8_forest_t8_gtest_forest_new_from_part_cmesh_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_forest_t8_gtest_transform_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_forest_t8_gtest_ghost_exchange_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_forest_t8_gtest_ghost_delete_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)