Feature: extend the partition algorithm to handle weights

benchmarks/t8_time_prism_adapt.cxx

@@ -181,7 +181,7 @@ main (int argc, char **argv)
 
   if (sreturnA > BUFSIZ || sreturnB > BUFSIZ) {
     /* The usage string or help message was truncated */
-    /* Note: gcc >= 7.1 prints a warning if we 
+    /* Note: gcc >= 7.1 prints a warning if we
      * do not check the return value of snprintf. */
     t8_debugf ("Warning: Truncated usage string and help message to '%s' and '%s'\n", usage, help);
   }

example/forest/t8_test_face_iterate.cxx

@@ -212,7 +212,7 @@ main (int argc, char **argv)
 
   if (sreturnA > BUFSIZ || sreturnB > BUFSIZ) {
     /* The usage string or help message was truncated */
-    /* Note: gcc >= 7.1 prints a warning if we 
+    /* Note: gcc >= 7.1 prints a warning if we
      * do not check the return value of snprintf. */
     t8_debugf ("Warning: Truncated usage string and help message to '%s' and '%s'\n", usage, help);
   }

example/forest/t8_test_ghost.cxx

@@ -303,7 +303,7 @@ main (int argc, char **argv)
 
   if (sreturnA > BUFSIZ || sreturnB > BUFSIZ) {
     /* The usage string or help message was truncated */
-    /* Note: gcc >= 7.1 prints a warning if we 
+    /* Note: gcc >= 7.1 prints a warning if we
      * do not check the return value of snprintf. */
     t8_debugf ("Warning: Truncated usage string and help message to '%s' and '%s'\n", usage, help);
   }

example/forest/t8_test_ghost_large_level_diff.cxx

@@ -219,7 +219,7 @@ main (int argc, char *argv[])
 
   if (sreturn >= BUFSIZ) {
     /* The help message was truncated */
-    /* Note: gcc >= 7.1 prints a warning if we 
+    /* Note: gcc >= 7.1 prints a warning if we
      * do not check the return value of snprintf. */
     t8_debugf ("Warning: Truncated help message to '%s'\n", help);
   }

example/remove/t8_example_empty_trees.cxx

@@ -115,7 +115,7 @@ main (int argc, char **argv)
 
   if (sreturnA > BUFSIZ || sreturnB > BUFSIZ) {
     /* The usage string or help message was truncated */
-    /* Note: gcc >= 7.1 prints a warning if we 
+    /* Note: gcc >= 7.1 prints a warning if we
      * do not check the return value of snprintf. */
     t8_debugf ("Warning: Truncated usage string and help message to '%s' and '%s'\n", usage, help);
   }

src/t8_forest/t8_forest.cxx

@@ -83,7 +83,7 @@ t8_forest_is_incomplete_family (const t8_forest_t forest, const t8_locidx_t ltre
   scheme->element_new (tree_class, 1, &element_parent_current);
   scheme->element_new (tree_class, 1, &element_compare);
 
-  /* We first assume that we have an (in)complete family with the size of array elements. 
+  /* We first assume that we have an (in)complete family with the size of array elements.
    * In the following we try to disprove this. */
   int family_size = elements_size;
 
@@ -92,9 +92,9 @@ t8_forest_is_incomplete_family (const t8_forest_t forest, const t8_locidx_t ltre
   const int child_id_current = scheme->element_get_child_id (tree_class, elements[0]);
   scheme->element_get_parent (tree_class, elements[0], element_parent_current);
 
-  /* Elements of the current family could already be passed, so that 
+  /* Elements of the current family could already be passed, so that
    * the element/family currently under consideration can no longer be coarsened.
-   * Also, there may be successors of a hypothetical previous family member 
+   * Also, there may be successors of a hypothetical previous family member
    * that would be overlapped after coarsening.
    * */
   if (child_id_current > 0 && el_considered > 0) {
@@ -136,12 +136,12 @@ t8_forest_is_incomplete_family (const t8_forest_t forest, const t8_locidx_t ltre
   T8_ASSERT (family_size > 0);
   T8_ASSERT (family_size >= 0 && family_size <= elements_size);
 
-  /* There may be successors of a hypothetical later family member (with index 
+  /* There may be successors of a hypothetical later family member (with index
    * family_size in this family) that would be overlapped after coarsening. */
   if (family_size < elements_size) {
     /* Get level of element after last element of current possible family */
     const int level = scheme->element_get_level (tree_class, elements[family_size]);
-    /* Only elements with higher level then level of current element, can get 
+    /* Only elements with higher level then level of current element, can get
      * potentially be overlapped. */
     if (level > level_current) {
       /* Compare ancestors */
@@ -164,7 +164,7 @@ t8_forest_is_incomplete_family (const t8_forest_t forest, const t8_locidx_t ltre
   const int num_siblings = scheme->element_get_num_siblings (tree_class, elements[0]);
   T8_ASSERT (family_size <= num_siblings);
   /* If the first/last element at a process boundary is not the first/last
-   * element of a possible family, we are not guaranteed to consider all 
+   * element of a possible family, we are not guaranteed to consider all
    * family members.*/
   if (el_considered == 0 && child_id_current > 0 && ltree_id == 0 && forest->mpirank > 0) {
     return 0;
@@ -300,9 +300,9 @@ t8_forest_no_overlap ([[maybe_unused]] t8_forest_t forest)
        * More detailed:
        * Let e_a and e_b be two elements.
        * If the level of e_a is equal to the level of the nca of e_a and e_b,
-       * then e_b is a descendant of e_a. 
+       * then e_b is a descendant of e_a.
        * If the level of e_b is equal to the level of the nca of e_a and e_b,
-       * then e_a is a descendant of e_b. 
+       * then e_a is a descendant of e_b.
        * Thus e_a and e_b overlap in both cases.
        * Note: If e_a equals e_b, e_a is the descendant of e_b and vice versa.
        * */
@@ -1305,9 +1305,9 @@ t8_forest_tree_shared ([[maybe_unused]] t8_forest_t forest, [[maybe_unused]] int
     else {
       SC_ABORT ("For incomplete trees the method t8_forest_last_tree_shared aka "
                 "t8_forest_tree_shared(forest, 1) is not implemented.\n");
-      /* TODO: If last_local_tree is 0 of the current process and it gets 0 as the 
-       * first_local_tree of the bigger process, then it cannot be said whether 
-       * the tree with id 0 is shared or not, since the bigger process could also 
+      /* TODO: If last_local_tree is 0 of the current process and it gets 0 as the
+       * first_local_tree of the bigger process, then it cannot be said whether
+       * the tree with id 0 is shared or not, since the bigger process could also
        * carry an empty forest. */
     }
     /* If global_neighbour_tree_idx == forest->first_local_tree tree is shared */
@@ -1960,7 +1960,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));
 
-  /* We get the array of the tree's elements and then search in the array of elements for our 
+  /* We get the array of the tree's elements and then search in the array of elements for our
    * element candidate. */
   /* Get the array */
   const t8_element_array_t *elements = t8_forest_get_tree_leaf_element_array (forest, local_tree);
@@ -1982,7 +1982,7 @@ t8_forest_element_is_leaf (const t8_forest_t forest, const t8_element_t *element
     /* The element was not found. */
     return 0;
   }
-  /* An element was found but it may not be the candidate element. 
+  /* An element was found but it may not be the candidate element.
    * To identify whether the element was found, we compare these two. */
   const t8_element_t *check_element = t8_element_array_index_locidx (elements, search_result);
   T8_ASSERT (check_element != NULL);
@@ -2323,7 +2323,7 @@ t8_forest_element_find_owner_old (t8_forest_t forest, t8_gloidx_t gtreeid, t8_el
     return proc;
   }
   else {
-    /* Get the next owning process. Its first descendant is in fact an element of the tree. 
+    /* Get the next owning process. Its first descendant is in fact an element of the tree.
      * If it is bigger than the descendant we look for, then proc is the owning process of element. */
     proc_next = *(int *) sc_array_index (owners_of_tree, 1);
     if (*(t8_linearidx_t *) t8_shmem_array_index (forest->global_first_desc, (size_t) proc_next)
@@ -2818,6 +2818,7 @@ t8_forest_set_partition (t8_forest_t forest, const t8_forest_t set_from, int set
   T8_ASSERT (forest->scheme == NULL);
 
   forest->set_for_coarsening = set_for_coarsening;
+  forest->weight_function = nullptr;
 
   if (set_from != NULL) {
     /* If set_from = NULL, we assume a previous forest_from was set */
@@ -2833,6 +2834,12 @@ t8_forest_set_partition (t8_forest_t forest, const t8_forest_t set_from, int set
   }
 }
 
+void
+t8_forest_set_partition_weights (t8_forest_t forest, t8_weight_fcn_t *weight_callback)
+{
+  forest->weight_function = weight_callback;
+}
+
 void
 t8_forest_set_balance (t8_forest_t forest, const t8_forest_t set_from, int no_repartition)
 {
@@ -2966,9 +2973,9 @@ t8_forest_comm_global_num_leaf_elements (t8_forest_t forest)
  * Check if any tree in a forest refines irregularly.
  * An irregular refining tree is a tree with an element that does not
  * refine into 2^dim children. For example the default implementation
- * of pyramids. 
+ * of pyramids.
  * \note This function is MPI collective
- * 
+ *
  * \param[in] forest    The forest to check
  * \return          Non-zero if any tree refines irregular
  */
@@ -3159,6 +3166,7 @@ t8_forest_commit (t8_forest_t forest)
           t8_forest_ref (forest->set_from);
         }
         t8_forest_set_partition (forest_partition, forest->set_from, forest->set_for_coarsening);
+        t8_forest_set_partition_weights (forest_partition, forest->weight_function);
         /* activate profiling, if this forest has profiling */
         t8_forest_set_profiling (forest_partition, forest->profile != NULL);
         /* Commit the partitioned forest */

src/t8_forest/t8_forest_general.h

@@ -48,11 +48,18 @@ typedef enum {
   T8_GHOST_VERTICES  /**< Consider all vertex (codimension 3) and edge and face neighbors. */
 } t8_ghost_type_t;
 
-/** This typedef is needed as a helper construct to 
+/** This typedef is needed as a helper construct to
  * properly be able to define a function that returns
  * a pointer to a void fun(void) function. \see t8_forest_get_user_function.
  */
 typedef void (*t8_generic_function_pointer) (void);
+
+/**
+ * The prototype a weight function for the partition algorithm.
+ * The function should be pure, and return a positive weight given a forest, a local tree index and an element index within the local tree
+ */
+typedef double (t8_weight_fcn_t) (t8_forest_t, t8_locidx_t, t8_locidx_t);
+
 T8_EXTERN_C_BEGIN ();
 
 /** Callback function prototype to replace one set of elements with another.
@@ -77,12 +84,12 @@ T8_EXTERN_C_BEGIN ();
  * \param [in] first_incoming  The tree local index of the first incoming element.
  *                             0 <= first_incom < new_which_tree->num_elements
  *
- * If an element is being refined, \a refine and \a num_outgoing will be 1 and 
+ * If an element is being refined, \a refine and \a num_outgoing will be 1 and
  * \a num_incoming will be the number of children.
- * If a family is being coarsened, \a refine will be -1, \a num_outgoing will be 
- * the number of family members and \a num_incoming will be 1. 
- * If an element is being removed, \a refine and \a num_outgoing will be 1 and 
- * \a num_incoming will be 0. 
+ * If a family is being coarsened, \a refine will be -1, \a num_outgoing will be
+ * the number of family members and \a num_incoming will be 1.
+ * If an element is being removed, \a refine and \a num_outgoing will be 1 and
+ * \a num_incoming will be 0.
  * Else \a refine will be 0 and \a num_outgoing and \a num_incoming will both be 1.
  * \see t8_forest_iterate_replace
  */
@@ -96,7 +103,7 @@ typedef void (*t8_forest_replace_t) (t8_forest_t forest_old, t8_forest_t forest_
  * form a family and we decide whether this family should be coarsened
  * or only the first element should be refined.
  * Otherwise \a is_family must equal zero and we consider the first entry
- * of the element array for refinement. 
+ * of the element array for refinement.
  * Entries of the element array beyond the first \a num_elements are undefined.
  * \param [in] forest       The forest to which the new elements belong.
  * \param [in] forest_from  The forest that is adapted.
@@ -121,7 +128,7 @@ typedef int (*t8_forest_adapt_t) (t8_forest_t forest, t8_forest_t forest_from, t
 
 /** Create a new forest with reference count one.
  * This forest needs to be specialized with the t8_forest_set_* calls.
- * Currently it is mandatory to either call the functions \see t8_forest_set_mpicomm, 
+ * Currently it is mandatory to either call the functions \see t8_forest_set_mpicomm,
  * \ref t8_forest_set_cmesh, and \ref t8_forest_set_scheme,
  * or to call one of \ref t8_forest_set_copy, \ref t8_forest_set_adapt, or
  * \ref t8_forest_set_partition.  It is illegal to mix these calls, or to
@@ -159,7 +166,7 @@ t8_forest_is_committed (t8_forest_t forest);
  * \param [in] forest   The forest to consider.
  * \return              True if \a forest has no elements which are inside each other.
  * \note This function is collective, but only checks local overlapping on each process.
- * \see t8_forest_partition_test_boundary_element if you also want to test for 
+ * \see t8_forest_partition_test_boundary_element if you also want to test for
  * global overlap across the process boundaries.
  */
 int
@@ -320,10 +327,21 @@ t8_forest_get_user_function (const t8_forest_t forest);
  * false, it is not necessary to call \ref t8_forest_set_partition additionally.
  * \note This setting may not be combined with \ref t8_forest_set_copy and overwrites
  * this setting.
+ * \note Calling this function sets the forest for unweighted partitioning. If weighted
+ * partitioning is desired, \ref t8_forest_set_partition_weights must be called after.
  */
 void
 t8_forest_set_partition (t8_forest_t forest, const t8_forest_t set_from, int set_for_coarsening);
 
+/** Set a user-defined weight function to guide the partitioning.
+ * \param [in, out] forest  The forest.
+ * \param [in]      weight_callback A callback function defining element weights for the partitioning.
+ * \pre \a weight_callback must be free of side-effects, the behavior is undefined otherwise
+ * \note If \a weight_callback is null, then all the elements are assumed to have the same weight
+ */
+void
+t8_forest_set_partition_weights (t8_forest_t forest, t8_weight_fcn_t *weight_callback);
+
 /** Set a source forest to be balanced during commit.
  * A forest is said to be balanced if each element has face neighbors of level
  * at most +1 or -1 of the element's level.
@@ -379,8 +397,8 @@ t8_forest_set_ghost_ext (t8_forest_t forest, int do_ghost, t8_ghost_type_t ghost
 
 /**
  *  Use assertions and document that the forest_set (..., from) and
- *  set_load are mutually exclusive. 
- * 
+ *  set_load are mutually exclusive.
+ *
  *  TODO: Unused function -> remove?
  */
 void
@@ -452,7 +470,7 @@ t8_forest_get_eclass (const t8_forest_t forest, const t8_locidx_t ltreeid);
 
 /**
  * Check whether a given tree id belongs to a local tree in a forest.
- * 
+ *
  * \param [in]    forest The forest.
  * \param [in]    local_tree A tree id.
  * \return True if and only if the id \a local_tree belongs to a local tree of \a forest.
@@ -467,8 +485,8 @@ t8_forest_tree_is_local (const t8_forest_t forest, const t8_locidx_t local_tree)
  * \param [in]      forest The forest.
  * \param [in]      gtreeid The global id of a tree.
  * \return                 The tree's local id in \a forest, if it is a local tree.
- *                         A negative number if not. Ghosts trees are not considered 
- *                         as local. 
+ *                         A negative number if not. Ghosts trees are not considered
+ *                         as local.
  * \see t8_forest_get_local_or_ghost_id for ghost trees.
  * \see https://github.com/DLR-AMR/t8code/wiki/Tree-indexing for more details about tree indexing.
  */
@@ -522,7 +540,7 @@ t8_forest_get_coarse_tree (t8_forest_t forest, t8_locidx_t ltreeid);
 
 /**
  * Query whether a given element is a leaf in a forest.
- * 
+ *
  * \param [in]  forest    The forest.
  * \param [in]  element   An element of a local tree in \a forest.
  * \param [in]  local_tree A local tree id of \a forest.
@@ -585,7 +603,7 @@ t8_forest_leaf_face_neighbors (t8_forest_t forest, t8_locidx_t ltreeid, const t8
                                t8_element_t **pneighbor_leaves[], int face, int *dual_faces[], int *num_neighbors,
                                t8_locidx_t **pelement_indices, t8_eclass_t *pneigh_eclass, int forest_is_balanced);
 
-/** Like \ref t8_forest_leaf_face_neighbors but also provides information about the global neighbors and the orientation. 
+/** Like \ref t8_forest_leaf_face_neighbors but also provides information about the global neighbors and the orientation.
  * \param [in]    forest  The forest. Must have a valid ghost layer.
  * \param [in]    ltreeid A local tree id.
  * \param [in]    leaf    A leaf in tree \a ltreeid of \a forest.
@@ -603,8 +621,8 @@ t8_forest_leaf_face_neighbors (t8_forest_t forest, t8_locidx_t ltreeid, const t8
  * \param [in]    forest_is_balanced True if we know that \a forest is balanced, false
  *                        otherwise.
  * \param [out]   gneigh_tree  The global tree IDs of the neighbor trees.
- * \param [out]   orientation  If not NULL on input, the face orientation is computed and stored here. 
- *                                         Thus, if the face connection is an inter-tree connection the orientation of the tree-to-tree connection is stored. 
+ * \param [out]   orientation  If not NULL on input, the face orientation is computed and stored here.
+ *                                         Thus, if the face connection is an inter-tree connection the orientation of the tree-to-tree connection is stored.
  *                                         Otherwise, the value 0 is stored.
  * All other parameters and behavior are identical to \ref t8_forest_leaf_face_neighbors.
  * \note If there are no face neighbors, then *neighbor_leaves = NULL, num_neighbors = 0,
@@ -854,24 +872,24 @@ t8_forest_element_face_neighbor (t8_forest_t forest, t8_locidx_t ltreeid, const
 
 /**
  * TODO: Can be removed since it is unused.
- * 
+ *
  * \param[in] forest The forest.
  */
 void
 t8_forest_iterate (t8_forest_t forest);
 
 /** Query whether a batch of points lies inside an element. For bilinearly interpolated elements.
  * \note For 2D quadrilateral elements this function is only an approximation. It is correct
- *  if the four vertices lie in the same plane, but it may produce only approximate results if 
+ *  if the four vertices lie in the same plane, but it may produce only approximate results if
  *  the vertices do not lie in the same plane.
  * \param [in]      forest      The forest.
  * \param [in]      ltreeid     The forest local id of the tree in which the element is.
  * \param [in]      element     The element.
  * \param [in]      points      3-dimensional coordinates of the points to check
  * \param [in]      num_points  The number of points to check
- * \param [in, out] is_inside   An array of length \a num_points, filled with 0/1 on output. True (non-zero) if a \a point 
- *                              lies within an \a element, false otherwise. The return value is also true if the point 
- *                              lies on the element boundary. Thus, this function may return true for different leaf 
+ * \param [in, out] is_inside   An array of length \a num_points, filled with 0/1 on output. True (non-zero) if a \a point
+ *                              lies within an \a element, false otherwise. The return value is also true if the point
+ *                              lies on the element boundary. Thus, this function may return true for different leaf
  *                              elements, if they are neighbors and the point lies on the common boundary.
  * \param [in]      tolerance   Tolerance that we allow the point to not exactly match the element.
  *                              If this value is larger we detect more points.