Feature: Tutorial 5 mesh handle version

mesh_handle/README.md

@@ -12,5 +12,6 @@ The folder's most important files are:
 - The [element.hxx](element.hxx) defines the elements (mesh or ghost elements) of the mesh handle.
 - The [competences.hxx](competences.hxx) defines additional competences/functionality of an element to access additional data.
 - The [competence_pack.hxx](competence_pack.hxx) is needed to pack element or mesh competences to pass it as template parameter to the mesh.
+- The [mesh_io.hxx](mesh_io.hxx) allows to output results in vtk format.
 
 Headers in the [internal/](internal/) folder are only intended to implement details of the mesh handle, so they should not need to be included for any other purpose.

mesh_handle/mesh_io.hxx

@@ -0,0 +1,92 @@
+/*
+  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 classes in parallel.
+
+  Copyright (C) 2026 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.
+*/
+
+/** \file mesh_io.hxx
+ * In- and output of meshes. For example, writing a mesh to vtk format.
+ */
+
+#pragma once
+#include <t8_forest/t8_forest_io.h>
+
+namespace t8_mesh_handle
+{
+/** 
+ * Write the mesh in a parallel vtu format. Extended version.
+ * See \see write_mesh_to_vtk for the standard version of this function.
+ * Writes one master .pvtu file and each process writes in its own .vtu file.
+ * If linked and not otherwise specified, the VTK API is used.
+ * If the VTK library is not linked, an ASCII file is written.
+ * This may change in accordance with \a write_ghosts, \a write_curved and 
+ * \a do_not_use_API, because the export of ghosts is not yet available with 
+ * the VTK API and the export of curved elements is not available with the
+ * inbuilt function to write ASCII files. The function will for example
+ * still use the VTK API to satisfy \a write_curved, even if \a do_not_use_API 
+ * is set to true.
+ * This function is collective and must be called on each process.
+ * \param [in] mesh             The mesh to write.
+ * \param [in] fileprefix       The prefix of the files where the vtk will be stored.
+ *             The master file is then fileprefix.pvtu and the process with rank r writes in the file fileprefix_r.vtu.
+ * \param [in] num_data         Number of user defined double valued data fields to write.
+ * \param [in] data             Array of t8_vtk_data_field_t of length \a num_data providing the user defined 
+ *              per element data. If scalar and vector fields are used, all scalar fields must come first in the array.
+ * \param [in] write_treeid     If true, the global tree id of the underlying forest is written for each element.
+ * \param [in] write_mpirank    If true, the mpirank is written for each element.
+ * \param [in] write_level      If true, the refinement level is written for each element.
+ * \param [in] write_element_id If true, the global element id is written for each element.
+ * \param [in] write_ghosts     If true, each process additionally writes its ghost elements.
+ *              For ghost element the treeid of the underlying forest is -1.
+ * \param [in] write_curved     If true, write the elements as curved element types from vtk.
+ * \param [in] do_not_use_API   Do not use the VTK API, even if linked and available.
+ * \return  True if successful, false if not (process local).
+ */
+template <typename TMeshClass>
+int
+write_mesh_to_vtk_ext (TMeshClass &mesh, const char *fileprefix, const int num_data, t8_vtk_data_field_t *data,
+                       bool write_treeid = false, bool write_mpirank = true, bool write_level = true,
+                       bool write_element_id = true, bool write_ghosts = false, bool write_curved = false,
+                       bool do_not_use_API = false)
+{
+  return t8_forest_write_vtk_ext (mesh->get_forest (), fileprefix, write_treeid, write_mpirank, write_level,
+                                  write_element_id, write_ghosts, write_curved, do_not_use_API, num_data, data);
+}
+
+/** 
+ * Write the mesh in a parallel vtu format. Writes one master
+ * .pvtu file and each process writes in its own .vtu file.
+ * If linked, the VTK API is used.
+ * If the VTK library is not linked, an ASCII file is written.
+ * This function writes the elements, level, mpirank, element id and the tree id of the underlying forest as data.
+ * This function is collective and must be called on each process.
+ * For more options use \see write_mesh_to_vtk_ext.
+ * \param [in] mesh       The mesh to write.
+ * \param [in] fileprefix The prefix of the files where the vtk will be stored.
+ *             The master file is then fileprefix.pvtu and the process with rank r writes in the file fileprefix_r.vtu.
+ * \return  True if successful, false if not (process local).
+ */
+template <typename TMeshClass>
+int
+write_mesh_to_vtk (TMeshClass &mesh, const char *fileprefix)
+{
+  return t8_forest_write_vtk (mesh->get_forest (), fileprefix);
+}
+
+}  // namespace t8_mesh_handle

tutorials/CMakeLists.txt

@@ -33,6 +33,10 @@ function( add_t8_tutorial )
 
     add_executable( ${ADD_T8_TUTORIAL_NAME} ${ADD_T8_TUTORIAL_SOURCES} )
     target_link_libraries( ${ADD_T8_TUTORIAL_NAME} PRIVATE T8 SC::SC )
+    string( FIND ${ADD_T8_TUTORIAL_NAME} "mesh_handle" is_mesh_handle_file )
+    if ( (${is_mesh_handle_file} GREATER_EQUAL 0) )
+        target_link_libraries(${ADD_T8_TUTORIAL_NAME} PRIVATE T8_MESH_HANDLE )
+    endif ()
     target_include_directories( ${ADD_T8_TUTORIAL_NAME} PRIVATE ${CMAKE_CURRENT_LIST_DIR}/.. )
 
     set_target_properties(${ADD_T8_TUTORIAL_NAME} PROPERTIES
@@ -69,3 +73,7 @@ copy_tutorial_file (features/t8_features_curved_meshes_generate_cmesh_hex.geo)
 copy_tutorial_file (features/t8_features_curved_meshes_generate_cmesh_quad.geo)
 copy_tutorial_file (features/t8_features_curved_meshes_generate_cmesh_tet.geo)
 copy_tutorial_file (features/t8_features_curved_meshes_generate_cmesh_tri.geo)
+
+if( T8CODE_BUILD_MESH_HANDLE )
+    add_t8_tutorial( NAME t8_step5_mesh_handle SOURCES general/t8_step5_mesh_handle.cxx )
+endif()

tutorials/general/t8_step5_mesh_handle.cxx

@@ -0,0 +1,239 @@
+/*
+  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) 2026 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.
+*/
+
+/** \file t8_step5_mesh_handle.cxx
+ *
+ * This is the same as t8_step5_element_data.cxx but using the mesh handle interface instead of the forest interface.
+ */
+
+#include <t8.h>
+
+#include <mesh_handle/mesh.hxx>
+#include <mesh_handle/competence_pack.hxx>
+#include <mesh_handle/constructor_wrappers.hxx>
+#include <mesh_handle/mesh_io.hxx>
+#include <t8_types/t8_vec.hxx>
+#include <memory>
+#include <span>
+
+/* The data that we want to store for each element.
+ * In this example we want to store the element's level and volume. */
+struct t8_step5_data_per_element
+{
+  int level;
+  double volume;
+};
+
+/** User data type we will pass to the adapt callback. */
+struct user_data
+{
+  t8_3D_vec midpoint;               /**< The midpoint of our sphere. */
+  double refine_if_inside_radius;   /**< If an element's center is smaller than this value, we refine the element. */
+  double coarsen_if_outside_radius; /**< If an element's center is larger this value, we coarsen its family. */
+};
+
+/** The adaptation callback function. This will refine elements inside of a given sphere and coarsen the elements 
+ * outside of a given sphere.
+ * \tparam TMeshClass    The mesh handle class.
+ * \param [in] mesh      The mesh that should be adapted.
+ * \param [in] elements  One element or a family of elements to consider for adaptation.
+ * \param [in] user_data The user data to be used during the adaptation process.
+ * \return 1 if the first entry in \a elements should be refined,
+ *        -1 if the family \a elements shall be coarsened,
+ *         0 else.
+ */
+template <typename TMeshClass>
+int
+adapt_callback ([[maybe_unused]] const TMeshClass &mesh, std::span<const typename TMeshClass::element_class> elements,
+                const user_data &user_data)
+{
+  auto element_centroid = elements[0].get_centroid ();
+  double dist = t8_dist<t8_3D_vec, t8_3D_vec> (element_centroid, user_data.midpoint);
+  if (dist < user_data.refine_if_inside_radius) {
+    return 1;
+  }
+  // Check if we got a family and if yes, if we should coarsen.
+  if ((elements.size () > 1) && (dist > user_data.coarsen_if_outside_radius)) {
+    return -1;
+  }
+  return 0;
+}
+
+/** Build a mesh with initial uniform refinement level \a level which is adapted according to \ref adapt_callback, 
+ * partitioned and balanced afterwards  and ghost elements are set.
+ * \tparam TMeshClass    The mesh handle class.
+ * \param [in] comm     MPI communicator to use.
+ * \param [in] level Initial refinement level.
+ * \return Unique pointer to the mesh created.
+ */
+template <typename TMesh>
+std::unique_ptr<TMesh>
+t8_step5_build_mesh (sc_MPI_Comm comm, int level)
+{
+  auto mesh_handle = t8_mesh_handle::handle_hypercube_hybrid_uniform_default<TMesh> (level, comm);
+  struct user_data adapt_data = {
+    { 0.5, 0.5, 1 }, /* Midpoint of the sphere. */
+    0.2,             /* Refine if inside this radius. */
+    0.4              /* Coarsen if outside this radius. */
+  };
+  /* Adapt, partition, balance and create ghost elements. */
+  mesh_handle->set_balance ();
+  mesh_handle->set_partition ();
+  mesh_handle->set_adapt (TMesh::template mesh_adapt_callback_wrapper<user_data> (adapt_callback<TMesh>, adapt_data),
+                          false);
+  mesh_handle->set_ghost ();
+  mesh_handle->commit ();
+  return mesh_handle;
+}
+
+/** Set element data to the mesh handle. */
+template <typename TMeshClass>
+void
+t8_step5_set_element_data (TMeshClass &mesh)
+{
+  for (auto &elem : *mesh) {
+    elem.set_element_data ({ elem.get_level (), elem.get_volume () });
+  }
+}
+
+/** Exchange element data set in \ref t8_step5_set_element_data for ghost elements. */
+template <typename TMeshClass>
+void
+t8_step5_exchange_ghost_data (TMeshClass &mesh)
+{
+  mesh->exchange_ghost_data ();
+}
+
+/** Write the mesh as vtu and also write the element's volumes in the file.
+ * t8code supports writing element based data to vtu as long as its stored
+ * as doubles. Each of the data fields to write has to be provided in its own
+ * array of length num_local_elements.
+ * We support two types: T8_VTK_SCALAR - One double per element
+ *                  and  T8_VTK_VECTOR - 3 doubles per element
+ */
+template <typename TMeshClass>
+static void
+t8_step5_output_data_to_vtu (TMeshClass &mesh, const char *prefix)
+{
+  t8_locidx_t num_elements = mesh->get_num_local_elements ();
+  /* We need to allocate a new array to store the volumes on their own.
+   * This array has one entry per local element. */
+  double *element_volumes = T8_ALLOC (double, num_elements);
+  /* The number of user defined data fields to write. */
+  int num_data = 1;
+  /* For each user defined data field we need one t8_vtk_data_field_t variable */
+  t8_vtk_data_field_t vtk_data;
+  /* Set the type of this variable. Since we have one value per element, we pick T8_VTK_SCALAR. */
+  vtk_data.type = T8_VTK_SCALAR;
+  /* The name of the field as should be written to the file. */
+  strcpy (vtk_data.description, "Element volume");
+  vtk_data.data = element_volumes;
+  /* Copy the element's volumes from our data array to the output array. */
+  for (t8_locidx_t ielem = 0; ielem < num_elements; ++ielem) {
+    element_volumes[ielem] = ((*mesh)[ielem]).get_element_data ().volume;
+  }
+  /* To write user defined data, we need the extended output function write_mesh_to_vtk_ext. 
+   * Despite writing user data, it also offers more control over which properties to write. */
+  t8_mesh_handle::write_mesh_to_vtk_ext (mesh, prefix, num_data, &vtk_data);
+  T8_FREE (element_volumes);
+}
+
+/** Main function. */
+int
+t8_step5_main (int argc, char **argv)
+{
+  /* The prefix for our output files. */
+  const char *prefix_mesh = "t8_step5_handle";
+  const char *prefix_mesh_with_data = "t8_step5_handle_with_volume_data";
+  /* The initial uniform refinement level. */
+  const int level = 3;
+
+  /* Initialize MPI. This has to happen before we initialize sc or t8code. */
+  int mpiret = sc_MPI_Init (&argc, &argv);
+  /* Error check the MPI return value. */
+  SC_CHECK_MPI (mpiret);
+  /* 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);
+  /* We will use MPI_COMM_WORLD as a communicator. */
+  sc_MPI_Comm comm = sc_MPI_COMM_WORLD;
+
+  /* Print a message on the root process. */
+  t8_global_productionf (" [step5] \n");
+  t8_global_productionf (" [step5] Hello, this is the step5 example of t8code using the mesh handle.\n");
+  t8_global_productionf (
+    " [step5] In this example we will store data on our elements and exchange the data of ghost elements.\n");
+  t8_global_productionf (" [step5] \n");
+
+  /* Setup: Build cmesh and adapt uniformly. Adapt similar to step 3 & 4. */
+  t8_global_productionf (" [step5] \n");
+  t8_global_productionf (" [step5] Creating an adapted mesh as in step3.\n");
+  t8_global_productionf (" [step5] \n");
+  { /* We put the mesh in its own scope so that it is automatically destroyed at the end of the scope. 
+     * This is only necessary because sc_finalize checks if there are leftover references. 
+     * This unique pointer would have been destroyed automatically at the end of the programme. */
+    using mesh_class = t8_mesh_handle::mesh<t8_mesh_handle::competence_pack<>, t8_step5_data_per_element>;
+    auto mesh = t8_step5_build_mesh<mesh_class> (comm, level);
+
+    t8_mesh_handle::write_mesh_to_vtk (mesh, prefix_mesh);
+    t8_global_productionf (" [step5] Wrote mesh to vtu files: %s*\n", prefix_mesh);
+
+    t8_step5_set_element_data (mesh);
+    t8_global_productionf (" [step5] Computed level and volume data for local elements.\n");
+    if (mesh->get_num_local_elements () > 0) {
+      /* Output the stored data of the first local element (if it exists). */
+      t8_global_productionf (" [step5] Element 0 has level %i and volume %e.\n", ((*mesh)[0]).get_element_data ().level,
+                             ((*mesh)[0]).get_element_data ().volume);
+    }
+
+    /* Exchange the data values of the ghost elements. */
+    t8_step5_exchange_ghost_data (mesh);
+    t8_global_productionf (" [step5] Exchanged ghost data.\n");
+    if (mesh->get_num_ghosts () > 0) {
+      /* Output the data of the first ghost element (if it exists). */
+      t8_locidx_t first_ghost_index = mesh->get_num_local_elements ();
+      t8_global_productionf (" [step5] Ghost 0 has level %i and volume %e.\n",
+                             ((*mesh)[first_ghost_index]).get_element_data ().level,
+                             ((*mesh)[first_ghost_index]).get_element_data ().volume);
+    }
+
+    /* Output the volume data to vtu. */
+    t8_step5_output_data_to_vtu (mesh, prefix_mesh_with_data);
+    t8_global_productionf (" [step5] Wrote mesh and volume data to %s*.\n", prefix_mesh_with_data);
+
+    /* Cleanup. */
+  }  // End scope of mesh
+  sc_finalize ();
+
+  mpiret = sc_MPI_Finalize ();
+  SC_CHECK_MPI (mpiret);
+
+  return 0;
+}
+
+/** Entry point of the program. */
+int
+main (int argc, char **argv)
+{
+  return t8_step5_main (argc, argv);
+}