Merge branch 'main' into remove-fortran-dep-when-not-needed

CMakeLists.txt

@@ -55,9 +55,9 @@ set( CMAKE_C_STANDARD_REQUIRED ON )
 set( CMAKE_C_EXTENSIONS OFF )
 
 if( NOT DEFINED CMAKE_CXX_STANDARD ) 
-  set( CMAKE_CXX_STANDARD 17 )
-elseif( CMAKE_CXX_STANDARD LESS 17 )
-  message( FATAL_ERROR "The CMAKE_CXX_STANDARD variable has been set to ${CMAKE_CXX_STANDARD}, but t8code requires C++17" )
+  set( CMAKE_CXX_STANDARD 20 )
+elseif( CMAKE_CXX_STANDARD LESS 20 )
+  message( FATAL_ERROR "The CMAKE_CXX_STANDARD variable has been set to ${CMAKE_CXX_STANDARD}, but t8code requires C++20" )
 endif()
 set( CMAKE_CXX_STANDARD_REQUIRED ON )
 set( CMAKE_CXX_EXTENSIONS OFF )

src/t8.h

@@ -73,6 +73,13 @@ T8_EXTERN_C_BEGIN ();
 #define T8_ASSERT(c) SC_NOOP ()
 #endif
 
+/**Extended T8_ASSERT assertion with custom error message. Only active in debug-mode. */
+#ifdef T8_ENABLE_DEBUG
+#define T8_ASSERTF(c, msg) SC_CHECK_ABORT ((c), "Assertion '" #c "': " msg)
+#else
+#define T8_ASSERTF(c, msg) SC_NOOP ()
+#endif
+
 /** Allocate a \a t-array with \a n elements. */
 #define T8_ALLOC(t, n) (t *) sc_malloc (t8_get_package_id (), (n) * sizeof (t))
 

src/t8_cmesh/t8_cmesh_helpers.cxx

@@ -69,6 +69,9 @@ t8_cmesh_set_join_by_vertices (t8_cmesh_t cmesh, const t8_gloidx_t ntrees, const
     }
   }
 
+  /* Scaled tolerance to decide if two doubles are equal. */
+  const double tolerance = 10.0 * T8_PRECISION_EPS * std::abs (max_coord - min_coord);
+
   /* Setup hash table `faces` mapping a hash key to a pair containing `(itree, iface)`. */
   std::multimap<unsigned long, std::pair<int, int>> faces;
 
@@ -154,7 +157,7 @@ t8_cmesh_set_join_by_vertices (t8_cmesh_t cmesh, const t8_gloidx_t ntrees, const
                                                                    neigh_itree, neigh_ivert, icoord)];
 
               /* Compare the coordinates with some tolerance. */
-              if (fabs (face_vert - neigh_face_vert) < 10.0 * T8_PRECISION_EPS) {
+              if (std::abs (face_vert - neigh_face_vert) < tolerance) {
                 match_count_per_coord++;
               }
             }

src/t8_data/t8_containers.cxx

@@ -114,12 +114,15 @@ t8_element_array_init_size (t8_element_array_t *element_array, t8_eclass_scheme_
 }
 
 void
-t8_element_array_init_view (t8_element_array_t *view, t8_element_array_t *array, size_t offset, size_t length)
+t8_element_array_init_view (t8_element_array_t *view, const t8_element_array_t *array, const size_t offset,
+                            const size_t length)
 {
   T8_ASSERT (t8_element_array_is_valid (array));
 
-  /* Initialize the element array */
-  sc_array_init_view (&view->array, &array->array, offset, length);
+  /* Initialize the element array.
+   * Unfortunately, we have to cast away the constness to pass to sc_array_init_view.
+   */
+  sc_array_init_view (&view->array, &((t8_element_array_t *) array)->array, offset, length);
   /* Set the scheme */
   view->scheme = array->scheme;
   T8_ASSERT (t8_element_array_is_valid (view));

src/t8_data/t8_containers.h

@@ -91,7 +91,8 @@ t8_element_array_init_size (t8_element_array_t *element_array, t8_eclass_scheme_
  *                       It is not necessary to call sc_array_reset later.
  */
 void
-t8_element_array_init_view (t8_element_array_t *view, t8_element_array_t *array, size_t offset, size_t length);
+t8_element_array_init_view (t8_element_array_t *view, const t8_element_array_t *array, const size_t offset,
+                            const size_t length);
 
 /** Initializes an already allocated (or static) view from given plain C data
  * (array of t8_element_t).

src/t8_forest/t8_forest.cxx

@@ -2508,10 +2508,15 @@ t8_forest_element_owners_at_face_recursion (t8_forest_t forest, t8_gloidx_t gtre
     else {
       last_owner_entry = -1;
     }
-    if (first_owner != last_owner_entry) {
+
+    if (first_owner > last_owner_entry) {
       /* We did not count this process as an owner, thus we add it */
       *(int *) sc_array_push (owners) = first_owner;
     }
+    if (last_owner > last_owner_entry) {
+      /* We did not count this process as an owner, thus we add it */
+      *(int *) sc_array_push (owners) = last_owner;
+    }
     T8_ASSERT (t8_forest_element_check_owner (forest, first_face_desc, gtreeid, eclass, first_owner, 1));
     T8_ASSERT (t8_forest_element_check_owner (forest, last_face_desc, gtreeid, eclass, first_owner, 1));
     /* free memory */
@@ -3380,6 +3385,7 @@ t8_forest_commit (t8_forest_t forest)
   }
 
   if (forest->mpisize > 1) {
+    sc_MPI_Barrier (forest->mpicomm);
     /* Construct a ghost layer, if desired */
     if (forest->do_ghost) {
       /* TODO: ghost type */

src/t8_forest/t8_forest_ghost.cxx

@@ -1459,19 +1459,6 @@ t8_forest_ghost_create_ext (t8_forest_t forest, int unbalanced_version)
     t8_forest_ghost_send_end (forest, ghost, send_info, requests);
   }
 
-  if (create_element_array) {
-    /* Free the offset memory, if created */
-    t8_shmem_array_destroy (&forest->element_offsets);
-  }
-  if (create_tree_array) {
-    /* Free the offset memory, if created */
-    t8_shmem_array_destroy (&forest->tree_offsets);
-  }
-  if (create_gfirst_desc_array) {
-    /* Free the offset memory, if created */
-    t8_shmem_array_destroy (&forest->global_first_desc);
-  }
-
   if (forest->profile != NULL) {
     /* If profiling is enabled, we measure the runtime of ghost_create */
     forest->profile->ghost_runtime += sc_MPI_Wtime ();
@@ -1491,6 +1478,19 @@ t8_forest_ghost_create_ext (t8_forest_t forest, int unbalanced_version)
     t8_global_productionf ("End ghost at %f  %f\n", sc_MPI_Wtime (), forest->profile->ghost_runtime);
   }
 
+  if (create_element_array) {
+    /* Free the offset memory, if created */
+    t8_shmem_array_destroy (&forest->element_offsets);
+  }
+  if (create_tree_array) {
+    /* Free the offset memory, if created */
+    t8_shmem_array_destroy (&forest->tree_offsets);
+  }
+  if (create_gfirst_desc_array) {
+    /* Free the offset memory, if created */
+    t8_shmem_array_destroy (&forest->global_first_desc);
+  }
+
   t8_global_productionf ("Done t8_forest_ghost with %i local elements and %i"
                          " ghost elements.\n",
                          t8_forest_get_local_num_elements (forest), t8_forest_get_num_ghosts (forest));

src/t8_forest/t8_forest_iterate.cxx

@@ -51,7 +51,7 @@ t8_forest_determine_child_type (sc_array_t *leaf_elements, size_t index, void *d
 }
 
 void
-t8_forest_split_array (const t8_element_t *element, t8_element_array_t *leaf_elements, size_t *offsets)
+t8_forest_split_array (const t8_element_t *element, const t8_element_array_t *leaf_elements, size_t *offsets)
 {
   sc_array_t offset_view;
   t8_forest_child_type_query_t query_data;
@@ -62,21 +62,22 @@ t8_forest_split_array (const t8_element_t *element, t8_element_array_t *leaf_ele
   query_data.level = ts->t8_element_level (element);
   query_data.ts = ts;
 
-  sc_array_t *element_array = t8_element_array_get_array_mutable (leaf_elements);
+  const sc_array_t *element_array = t8_element_array_get_array (leaf_elements);
   /* Split the elements array according to the elements' ancestor id at
    * the given level. In other words for each child C of element, find
    * the indices i, j such that all descendants of C are
    * elements[i], ..., elements[j-1]
    */
   sc_array_init_data (&offset_view, offsets, sizeof (size_t), query_data.num_children + 1);
-  sc_array_split (element_array, &offset_view, query_data.num_children, t8_forest_determine_child_type,
+  // Unfortunately, we have to cast away constness to pass to sc_array_split
+  sc_array_split ((sc_array_t *) element_array, &offset_view, query_data.num_children, t8_forest_determine_child_type,
                   (void *) &query_data);
 }
 
 void
 t8_forest_iterate_faces (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, int face,
-                         t8_element_array_t *leaf_elements, void *user_data, t8_locidx_t tree_lindex_of_first_leaf,
-                         t8_forest_iterate_face_fn callback)
+                         const t8_element_array_t *leaf_elements, void *user_data,
+                         t8_locidx_t tree_lindex_of_first_leaf, t8_forest_iterate_face_fn callback)
 {
   t8_eclass_scheme_c *ts;
   t8_eclass_t eclass;

src/t8_forest/t8_forest_iterate.h

@@ -82,9 +82,21 @@ typedef void (*t8_forest_query_fn) (t8_forest_t forest, const t8_locidx_t ltreei
 
 T8_EXTERN_C_BEGIN ();
 
-/* TODO: Document */
+/** Split an array of elements according to the children of a given element E.
+ *  In other words for each child C of E, find
+ *  the index i, j, such that all descendants of C are
+ *  elements[i], ..., elements[j-1]. 
+ * 
+ * \param [in] element       An element.
+ * \param [in] leaf_elements An array of leaf elements of \a element. Thus, all
+ *                           elements must be descendants. Sorted by linear index.
+ * \param [in,out] offsets   On input an allocated array of \a num_children_of_E + 1 entries.  
+ *                           On output entry i indicates the position in \a leaf_elements
+ *                           where the descandents of the i-th child of E start.
+ */
+
 void
-t8_forest_split_array (const t8_element_t *element, t8_element_array_t *leaf_elements, size_t *offsets);
+t8_forest_split_array (const t8_element_t *element, const t8_element_array_t *leaf_elements, size_t *offsets);
 
 /* TODO: comment */
 /* Iterate over all leaves of an element that touch a given face of the element */
@@ -96,8 +108,8 @@ t8_forest_split_array (const t8_element_t *element, t8_element_array_t *leaf_ele
  * If it returns false, the current element is not traversed further */
 void
 t8_forest_iterate_faces (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, int face,
-                         t8_element_array_t *leaf_elements, void *user_data, t8_locidx_t tree_lindex_of_first_leaf,
-                         t8_forest_iterate_face_fn callback);
+                         const t8_element_array_t *leaf_elements, void *user_data,
+                         t8_locidx_t tree_lindex_of_first_leaf, t8_forest_iterate_face_fn callback);
 
 /* Perform a top-down search of the forest, executing a callback on each
  * intermediate element. The search will enter each tree at least once.

src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.cxx

@@ -156,7 +156,7 @@ t8_geometry_tessellated_spherical_surface::t8_geom_evaluate (t8_cmesh_t cmesh, t
   const double distance = std::abs (t8_vec_dot (active_tree_vertices, normal));
 
   // Compute actual radius of the sphere.
-  const double radius = distance * std::cbrt (1.0);
+  const double radius = distance * std::sqrt (3.0);
 
   // Compute orthogonal coordinate system anchored on the cmesh element.
   t8_vec_orthogonal_tripod (normal, tangent1, tangent2);
@@ -212,10 +212,10 @@ t8_geometry_cubed_spherical_shell::t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx
   const double distance = std::abs (t8_vec_dot (active_tree_vertices, normal));
 
   // Compute actual radius of the sphere.
-  const double CBRT = std::cbrt (1.0);
-  const double inner_radius = distance * CBRT;
+  const double SQRT3 = std::sqrt (3.0);
+  const double inner_radius = distance * SQRT3;
   const double shell_thickness
-    = std::abs (t8_vec_dot (active_tree_vertices + t8_eclass_num_vertices[active_tree_class] * 3 / 2, normal)) * CBRT
+    = std::abs (t8_vec_dot (active_tree_vertices + t8_eclass_num_vertices[active_tree_class] * 3 / 2, normal)) * SQRT3
       - inner_radius;
 
   // Compute orthogonal coordinate system anchored on the cmesh element.