Feature new for hybrid

example/cmesh/Makefile.am

@@ -6,6 +6,7 @@ bin_PROGRAMS += \
   example/cmesh/t8_cmesh_partition \
   example/cmesh/t8_cmesh_create_partitioned \
   example/cmesh/t8_cmesh_refine \
+  example/cmesh/t8_cmesh_hybrid_new \
   example/cmesh/t8_cmesh_set_join_by_vertices \
   example/cmesh/t8_cmesh_geometry_examples \
   example/cmesh/t8_cmesh_hypercube_pad
@@ -15,4 +16,6 @@ example_cmesh_t8_cmesh_refine_SOURCES = example/cmesh/t8_cmesh_refine.cxx
 example_cmesh_t8_cmesh_set_join_by_vertices_SOURCES = example/cmesh/t8_cmesh_set_join_by_vertices.cxx
 example_cmesh_t8_cmesh_geometry_examples_SOURCES = example/cmesh/t8_cmesh_geometry_examples.cxx
 example_cmesh_t8_cmesh_create_partitioned_SOURCES = example/cmesh/t8_cmesh_create_partitioned.cxx
+example_cmesh_t8_cmesh_hybrid_new_SOURCES = example/cmesh/t8_cmesh_hybrid_new.cxx
 example_cmesh_t8_cmesh_hypercube_pad_SOURCES = example/cmesh/t8_cmesh_hypercube_pad.cxx
+

example/cmesh/t8_cmesh_hybrid_new.cxx

@@ -0,0 +1,244 @@
+/*
+  This file is part of t8code.
+  t8code is a C library to manage a collection (a forest) of multiple
+  connected adaptive space-trees of general element types in parallel.
+
+  Copyright (C) 2015 the developers
+
+  t8code is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2 of the License, or
+  (at your option) any later version.
+
+  t8code is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with t8code; if not, write to the Free Software Foundation, Inc.,
+  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+*/
+
+#include <sc_options.h>
+#include <sc_flops.h>
+#include <sc_refcount.h>
+#include <t8_schemes/t8_default/t8_default_cxx.hxx>
+#include <t8_schemes/t8_default/t8_default_tet/t8_dtet.h>
+#include <t8_schemes/t8_default/t8_default_prism/t8_dprism.h>
+#include <t8_schemes/t8_default/t8_default_pyramid/t8_dpyramid.h>
+#include <t8_forest/t8_forest.h>
+#include <t8_cmesh/t8_cmesh_examples.h>
+#include <t8_cmesh/t8_cmesh_partition.h>
+#include <t8_cmesh_readmshfile.h>
+#include <t8_cmesh_vtk_writer.h>
+#include <sc_statistics.h>
+
+/**
+ * Adaptation criterion for a forest containing hexahedra, tetrahedra, prisms and/or pryramids. 
+ * Refines every second element, except for pyramids. 
+ * Refines every pyramid
+ * 
+ * \param [in] forest the forest
+ * \param [in] forest_from 
+ * \param [in] which_tree The local id of the tree
+ * \param [in] lelement_id the local id of the element
+ * \param [in] ts the scheme to use
+ * \param [in] is_family flag, if the \a elements form a family
+ * \param [in] num_elements number of elements
+ * \param [in] elements A single element or a collection of \a num_elements
+ * \return int 
+ */
+static int
+t8_basic_hybrid_refine (t8_forest_t forest, t8_forest_t forest_from, t8_locidx_t which_tree, t8_locidx_t lelement_id,
+                        t8_eclass_scheme_c *ts, int is_family, int num_elements, t8_element_t *elements[])
+{
+  int level, id;
+  /*If the level is equal or higher than given by the user do not refine*/
+  level = ts->t8_element_level (elements[0]);
+  if (level >= *(int *) t8_forest_get_user_data (forest)) {
+    return 0;
+  }
+  else {
+    /*Refine every second element */
+    switch (ts->t8_element_shape (elements[0])) {
+    case T8_ECLASS_HEX:
+      id = ts->t8_element_child_id (elements[0]);
+      return id % 2 == 0 ? 1 : 0;
+    case T8_ECLASS_TET:
+      id = ts->t8_element_child_id (elements[0]);
+      return id % 2 == 0 ? 1 : 0;
+    case T8_ECLASS_PRISM:
+      id = ts->t8_element_child_id (elements[0]);
+      return id % 2 == 0 ? 1 : 0;
+    case T8_ECLASS_PYRAMID:
+      /* Refine every element. */
+      return 1;
+    default:
+      return 1;
+    }
+  }
+}
+
+/**
+ * Creates a cmesh and forest using this cmesh and times all high-level algorithms (new, adapt, ghost, partition, ...)
+ * 
+ * \param[in] level The initial level of forest
+ * \param[in] endlvl The maximal level of the forest
+ * \param[in] do_vtk Flag, if vtk-output should be produced
+ */
+static void
+t8_basic_hybrid (const int level, int endlvl, const int do_vtk)
+{
+  /* Create and initialize timing stats */
+  sc_statinfo_t times[6];
+  sc_stats_init (&times[0], "new");
+  sc_stats_init (&times[1], "adapt");
+  sc_stats_init (&times[2], "ghost");
+  sc_stats_init (&times[3], "partition");
+  sc_stats_init (&times[4], "total");
+  double total_time = -sc_MPI_Wtime ();
+
+  t8_cmesh_t cmesh;
+  /* Create the cmesh */
+  t8_global_productionf ("Constructing full hybrid mesh.\n");
+  cmesh = t8_cmesh_new_full_hybrid (sc_MPI_COMM_WORLD);
+
+  /* Partition the cmesh */
+  t8_cmesh_t cmesh_partition;
+  t8_cmesh_init (&cmesh_partition);
+  t8_cmesh_set_derive (cmesh_partition, cmesh);
+  t8_cmesh_set_partition_uniform (cmesh_partition, level, t8_scheme_new_default_cxx ());
+  t8_cmesh_commit (cmesh_partition, sc_MPI_COMM_WORLD);
+  cmesh = cmesh_partition;
+
+  /* Save the cmesh and create vtk output if the flag is set. */
+  char vtuname[BUFSIZ];
+  char cmesh_file[BUFSIZ];
+  snprintf (cmesh_file, BUFSIZ, "cmesh_hybrid");
+  snprintf (vtuname, BUFSIZ, "cmesh_hybrid");
+  t8_cmesh_save (cmesh, cmesh_file);
+  if (t8_cmesh_vtk_write_file (cmesh, vtuname, 1.0) == 0) {
+    t8_debugf ("Output to %s\n", vtuname);
+  }
+  else {
+    t8_debugf ("Error in output\n");
+  }
+
+  /* Create the forest and measure the time needed for the initial refinement*/
+  t8_forest_t forest;
+  t8_forest_init (&forest);
+  t8_forest_set_profiling (forest, 1);
+  t8_forest_set_cmesh (forest, cmesh, sc_MPI_COMM_WORLD);
+  t8_forest_set_scheme (forest, t8_scheme_new_default_cxx ());
+  t8_forest_set_level (forest, level);
+  double new_time = -sc_MPI_Wtime ();
+  t8_forest_commit (forest);
+  new_time += sc_MPI_Wtime ();
+
+  if (do_vtk) {
+    /* write the initially refined forest into a vtu-file*/
+    snprintf (vtuname, BUFSIZ, "forest_hybrid");
+    t8_forest_write_vtk (forest, vtuname);
+  }
+
+  /* Initialize the adapted forest*/
+  t8_forest_t forest_adapt;
+  t8_forest_init (&forest_adapt);
+  t8_forest_set_user_data (forest_adapt, &endlvl);
+  t8_forest_set_profiling (forest_adapt, 1);
+
+  /* Set the adaptation criterion to use. */
+  t8_forest_set_adapt (forest_adapt, forest, t8_basic_hybrid_refine, 1);
+
+  /* Measure the time needed to adapt */
+  double adapt_time = -sc_MPI_Wtime ();
+  t8_forest_commit (forest_adapt);
+  adapt_time += sc_MPI_Wtime ();
+
+  /* Initialize the partitioned forest */
+  t8_forest_t forest_partition;
+  t8_forest_init (&forest_partition);
+  t8_forest_set_partition (forest_partition, forest_adapt, 0);
+
+  /* We want to compute the ghosts, too */
+  t8_forest_set_ghost (forest_partition, 1, T8_GHOST_FACES);
+  /* Use the profiling algorithm to measure the time needed and also how many ghost elements are sent, 
+   *  how many procs we sent to and how many rounds we need to balance the forest. */
+  t8_forest_set_profiling (forest_partition, 1);
+  t8_forest_commit (forest_partition);
+
+  /* Write vtk-output*/
+  if (do_vtk) {
+    snprintf (vtuname, BUFSIZ, "forest_hybrid_partition");
+    t8_forest_write_vtk (forest_partition, vtuname);
+  }
+  total_time += sc_MPI_Wtime ();
+  int procs_sent;         /* Procs we sent to during forest_partition */
+  t8_locidx_t ghost_sent; /* Number of ghost-elements we sent */
+  t8_forest_print_profile (forest_partition);
+  const double partition_time = t8_forest_profile_get_partition_time (forest_partition, &procs_sent);
+  const double ghost_time = t8_forest_profile_get_ghost_time (forest_partition, &ghost_sent);
+
+  /* Accumulate the timings over all procs and print the stats */
+  sc_stats_accumulate (&times[0], new_time);
+  sc_stats_accumulate (&times[1], adapt_time);
+  sc_stats_accumulate (&times[2], ghost_time);
+  sc_stats_accumulate (&times[3], partition_time);
+  sc_stats_accumulate (&times[4], total_time);
+  sc_stats_compute (sc_MPI_COMM_WORLD, 5, times);
+  sc_stats_print (t8_get_package_id (), SC_LP_ESSENTIAL, 5, times, 1, 1);
+
+  t8_forest_unref (&forest_partition);
+}
+
+int
+main (int argc, char **argv)
+{
+  int level;
+  int endlvl;
+  int helpme;
+  int do_vtk;
+  sc_options_t *opt;
+  char help[BUFSIZ];
+
+  /* long help message */
+  snprintf (help, BUFSIZ,
+            "ADD EXAMPLE DESCRIPTION.\n\n"
+            "Usage:\t%s <OPTIONS>\n\t%s -h\t"
+            "for a brief overview of all options.\n",
+            basename (argv[0]), basename (argv[0]));
+  int mpiret = sc_MPI_Init (&argc, &argv);
+  SC_CHECK_MPI (mpiret);
+
+  sc_init (sc_MPI_COMM_WORLD, 1, 1, NULL, SC_LP_ESSENTIAL);
+  t8_init (SC_LP_DEFAULT);
+
+  opt = sc_options_new (argv[0]);
+  sc_options_add_switch (opt, 'h', "help", &helpme, "Display a short help message.");
+  sc_options_add_int (opt, 'l', "level", &level, 0, "The refinement level of the mesh.");
+  sc_options_add_int (opt, 'f', "final-level", &endlvl, 1, "The final refinement level of the mesh.");
+  sc_options_add_switch (opt, 'v', "vtk", &do_vtk, "Enable vtk-output.");
+
+  const int parsed = sc_options_parse (t8_get_package_id (), SC_LP_DEFAULT, opt, argc, argv);
+  if (helpme) {
+    /* display help message and usage */
+    t8_global_productionf ("%s\n", help);
+    sc_options_print_usage (t8_get_package_id (), SC_LP_ERROR, opt, NULL);
+  }
+  else if (parsed >= 0 && 0 <= level) {
+    t8_basic_hybrid (level, endlvl, do_vtk);
+  }
+  else {
+    /* wrong usage */
+    t8_global_productionf ("\n\t ERROR: Wrong usage.\n\n");
+    sc_options_print_usage (t8_get_package_id (), SC_LP_ERROR, opt, NULL);
+  }
+  sc_options_destroy (opt);
+  sc_finalize ();
+
+  mpiret = sc_MPI_Finalize ();
+  SC_CHECK_MPI (mpiret);
+
+  return 0;
+}

src/t8.c

@@ -165,3 +165,10 @@ t8_sc_array_index_locidx (const sc_array_t *array, const t8_locidx_t it)
 
   return array->array + array->elem_size * (size_t) it;
 }
+
+void *
+t8_sc_array_index_gloidx (const sc_array_t *array, const t8_gloidx_t it)
+{
+  T8_ASSERT (it >= 0 && (size_t) it < array->elem_count);
+  return array->array + array->elem_size * (size_t) it;
+}

src/t8.h

@@ -100,6 +100,8 @@ typedef int64_t t8_gloidx_t;
 #define T8_MPI_GLOIDX sc_MPI_LONG_LONG_INT
 /** Macro to get the absolute value of a t8_gloidx_t */
 #define T8_GLOIDX_ABS(x) ((t8_gloidx_t) llabs ((long long) (x)))
+/** Maximum possible value of a t8_gloidx_t*/
+#define T8_GLOIDX_MAX INT64_MAX
 /** Comparison function for t8_gloidx_t */
 #define t8_compare_gloidx(v, w) sc_int64_compare (v, w)
 
@@ -129,6 +131,8 @@ typedef enum {
   T8_MPI_PARTITION_FOREST,              /**< Used for forest partitioning */
   T8_MPI_GHOST_FOREST,                  /**< Used for for ghost layer creation */
   T8_MPI_GHOST_EXC_FOREST,              /**< Used for ghost data exchange */
+  T8_MPI_CMESH_UNIFORM_BOUNDS_START,    /**< Used for cmesh uniform bounds computation. */
+  T8_MPI_CMESH_UNIFORM_BOUNDS_END,      /**< Used for cmesh uniform bounds computation. */
   T8_MPI_TAG_LAST
 } t8_MPI_tag_t;
 
@@ -266,6 +270,13 @@ t8_init (int log_threshold);
 void *
 t8_sc_array_index_locidx (const sc_array_t *array, const t8_locidx_t it);
 
+/** Return a pointer to an array element indexed by a t8_gloidx_t.
+ * \param [in] index needs to be in [0]..[elem_count-1].
+ * \return           A void * pointing to entry \a it in \a array.
+ */
+void *
+t8_sc_array_index_gloidx (const sc_array_t *array, const t8_gloidx_t it);
+
 /* call this at the end of a header file to match T8_EXTERN_C_BEGIN (). */
 T8_EXTERN_C_END ();
 

src/t8_cmesh.h

@@ -668,13 +668,37 @@ t8_cmesh_get_partition_table (t8_cmesh_t cmesh);
  *                                  the calling processor anymore. Not computed if NULL.
  * \param [out]   first_tree_shared If not NULL, 1 or 0 is stored here depending on whether \a first_local_tree is the
  *                                 same as \a last_local_tree on the next process.
- * \a cmesh must be committed before calling this function. *
+ * \a cmesh must be committed before calling this function. 
  */
 void
 t8_cmesh_uniform_bounds (t8_cmesh_t cmesh, int level, t8_scheme_cxx_t *ts, t8_gloidx_t *first_local_tree,
                          t8_gloidx_t *child_in_tree_begin, t8_gloidx_t *last_local_tree, t8_gloidx_t *child_in_tree_end,
                          int8_t *first_tree_shared);
 
+/**
+ * Calculate the section of a uniform hybrid forest for the current rank. Need for hybrid meshes, especially 
+ * meshes where not all elements refine into 1:2^dim manner. The section is calculated without assuming such refinement
+ * and each process computes its number of elements on the given \var level communicates the number to other processes
+ * and the correct section is computed based on this information. 
+ * 
+ * \param [in] cmesh        The cmesh to be considered.
+ * \param [in] level        The uniform refinement level to be created.
+ * \param [in] scheme       The element scheme for which to compute the bounds.
+ * \param [out]   first_local_tree  The global index of the first tree that contains elements belonging to the calling processor.
+ * \param [out]   child_in_tree_begin The global index of the first element belonging to the calling processor. Not computed if NULL.
+ * \param [out]   last_local_tree  The global index of the last tree that contains elements belonging to the calling processor.
+ * \param [out]   child_in_tree_end The global index of the first element that does not belonging to
+ *                                  the calling processor anymore. Not computed if NULL.
+ * \param [out]   first_tree_shared If not NULL, 1 or 0 is stored here depending on whether \a first_local_tree is the
+ *                                 same as \a last_local_tree on the next process.
+ * \param [in] comm         The communicator 
+ */
+void
+t8_cmesh_uniform_bounds_hybrid (const t8_cmesh_t cmesh, const int level, const t8_scheme_cxx_t *scheme,
+                                t8_gloidx_t *first_local_tree, t8_gloidx_t *child_in_tree_begin,
+                                t8_gloidx_t *last_local_tree, t8_gloidx_t *child_in_tree_end, int8_t *first_tree_shared,
+                                sc_MPI_Comm comm);
+
 /** Increase the reference counter of a cmesh.
  * \param [in,out] cmesh        On input, this cmesh must exist with positive
  *                              reference count.  It may be in any state.