Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
DLR-AMR
GitHub Repository: DLR-AMR/t8code
Path: blob/main/src/t8_cmesh/t8_cmesh_io/t8_cmesh_save.cxx
927 views
/*
  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) 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.
*/

/** \file t8_cmesh_save.cxx
 * We define routines to save and load a cmesh to/from the file system.
 */

#include <t8_version.h>
#include <t8_eclass.h>
#include <t8_cmesh/t8_cmesh_internal/t8_cmesh_types.h>
#include <t8_cmesh/t8_cmesh_internal/t8_cmesh_trees.h>
#include <t8_cmesh/t8_cmesh_io/t8_cmesh_save.h>
#include <t8_cmesh/t8_cmesh_internal/t8_cmesh_partition.h>
#include <t8_cmesh/t8_cmesh_internal/t8_cmesh_offset.h>
#include <t8_geometry/t8_geometry.h>
#include <t8_geometry/t8_geometry_base.h>
#include <t8_geometry/t8_geometry_handler.hxx>
#include <t8_geometry/t8_geometry_implementations/t8_geometry_linear.h>

/**
 * This macro is called to check a condition and if not fulfilled
 * close the file and exit the function
 */
#define T8_SAVE_CHECK_CLOSE(x, fp) \
  if (!(x)) { \
    t8_errorf ("file i/o error. Condition %s not fulfilled. " \
               "Line %i in file %s\n", \
               #x, __LINE__, __FILE__); \
    fclose (fp); \
    return 0; \
  }

/* Write the neighbor data of all ghosts */
static int
t8_cmesh_save_ghost_neighbors (const t8_cmesh_t cmesh, FILE *fp)
{
  t8_locidx_t ighost;
  t8_cghost_t ghost;
  t8_gloidx_t *face_neigh;
  int8_t *ttf;
  int iface, num_faces;
  int ret;
  char buffer[BUFSIZ];

  ret = fprintf (fp, "\n--- Ghost neighbor section ---\n\n");
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  for (ighost = 0; ighost < cmesh->num_ghosts; ighost++) {
    /* Reset the buffer */
    buffer[0] = '\0';
    ghost = t8_cmesh_trees_get_ghost_ext (cmesh->trees, ighost, &face_neigh, &ttf);
    /* Iterate over all faces to write the face neighbor information */
    num_faces = t8_eclass_num_faces[ghost->eclass];
    for (iface = 0; iface < num_faces; iface++) {
      snprintf (buffer + strlen (buffer), BUFSIZ - strlen (buffer), "%lli %i%s", (long long) face_neigh[iface],
                ttf[iface], iface == num_faces - 1 ? "" : ", ");
    }
    ret = fprintf (fp, "%s\n", buffer);
    T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  }
  return 1;
}

static int
t8_cmesh_load_ghost_attributes (const t8_cmesh_t cmesh, FILE *fp)
{

  t8_locidx_t ighost;
  t8_cghost_t ghost;
  t8_gloidx_t *face_neigh;
  int8_t *ttf;
  int iface, num_faces, ttf_entry;
  long long neigh;
  int ret;

  ret = fscanf (fp, "\n--- Ghost neighbor section ---\n\n");
  T8_SAVE_CHECK_CLOSE (ret == 0, fp);

  for (ighost = 0; ighost < cmesh->num_ghosts; ighost++) {
    /* Get a pointer to the ghost */
    ghost = t8_cmesh_trees_get_ghost_ext (cmesh->trees, ighost, &face_neigh, &ttf);
    /* Read the neighbor information */
    num_faces = t8_eclass_num_faces[ghost->eclass];
    for (iface = 0; iface < num_faces; iface++) {
      ret = fscanf (fp, "%lli %i%*c", &neigh, &ttf_entry);
      T8_SAVE_CHECK_CLOSE (ret == 2, fp);
      ttf[iface] = ttf_entry;
      face_neigh[iface] = neigh;
    }
  }
  return 1;
}

static int
t8_cmesh_save_ghosts (const t8_cmesh_t cmesh, FILE *fp)
{
  t8_locidx_t ighost;
  t8_cghost_t ghost;
  int ret;

  ret = fprintf (fp, "\n--- Ghost section ---");
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  ret = fprintf (fp, "\n\nGhosts %li\n\n", (long) cmesh->num_ghosts);
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  for (ighost = 0; ighost < cmesh->num_ghosts; ighost++) {
    ghost = t8_cmesh_trees_get_ghost (cmesh->trees, ighost);
    ret = fprintf (fp, "treeid %lli\neclass %i\n\n", (long long) ghost->treeid, (int) ghost->eclass);
    T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  }
  return 1;
}

static int
t8_cmesh_load_ghosts (const t8_cmesh_t cmesh, FILE *fp)
{
  t8_locidx_t ighost;
  t8_cghost_t ghost;
  int ret;
  int eclass;
  long num_ghosts;
  long long global_id;

  ret = fscanf (fp, "\n--- Ghost section ---\n");
  T8_SAVE_CHECK_CLOSE (ret == 0, fp);
  /* Check whether the number of ghosts is correct */
  ret = fscanf (fp, "Ghosts %li\n", &num_ghosts);
  T8_SAVE_CHECK_CLOSE (ret == 1, fp);
  T8_SAVE_CHECK_CLOSE (num_ghosts == cmesh->num_ghosts, fp);
  for (ighost = 0; ighost < cmesh->num_ghosts; ighost++) {
    /* Get a pointer to the ghost */
    ghost = t8_cmesh_trees_get_ghost (cmesh->trees, ighost);
    /* Read the ghost's global id */
    ret = fscanf (fp, "treeid %lli\n", &global_id);
    T8_SAVE_CHECK_CLOSE (ret == 1, fp);
    T8_SAVE_CHECK_CLOSE (0 <= global_id && global_id < cmesh->num_trees, fp);
    ghost->treeid = global_id;
    /* read the ghost's eclass */
    ret = fscanf (fp, "eclass %i\n", &eclass);
    T8_SAVE_CHECK_CLOSE (ret == 1, fp);
    ghost->eclass = (t8_eclass_t) eclass;
    /* Ignore the neighbor offset */
    ret = fscanf (fp, "neigh_offset %*i\n");
    T8_SAVE_CHECK_CLOSE (ret == 0, fp);
    ret = fscanf (fp, "\n");
    T8_SAVE_CHECK_CLOSE (ret == 0, fp);
  }
  return 1;
}

/* Load all attributes that were stored in a file */
static int
t8_cmesh_load_tree_attributes (const t8_cmesh_t cmesh, FILE *fp)
{
  double *vertices = NULL;
  t8_locidx_t itree;
  long treeid, neighbor;
  t8_ctree_t tree;
  int att, i, num_vertices, num_faces, iface;
  int ret, ttf_entry;
  t8_locidx_t *face_neighbors;
  int8_t *ttf;
  t8_stash_attribute_struct_t att_struct;

  ret = fscanf (fp, "\n--- Tree attribute section ---\n");
  T8_SAVE_CHECK_CLOSE (ret == 0, fp);
  /* loop over all trees */
  for (itree = 0; itree < cmesh->num_local_trees; itree++) {
    tree = t8_cmesh_trees_get_tree_ext (cmesh->trees, itree, &face_neighbors, &ttf);
    /* Allocate memory for the temporary vertices array.
     * This memory is reused for each tree */
    num_vertices = t8_eclass_num_vertices[tree->eclass];
    ret = fscanf (fp, "tree %li\n", &treeid);
    T8_SAVE_CHECK_CLOSE (ret == 1, fp);
    /* Load the tree neighbors */
    num_faces = t8_eclass_num_faces[tree->eclass];
    ret = fscanf (fp, "Neighbors:");
    T8_SAVE_CHECK_CLOSE (ret == 0, fp);
    for (iface = 0; iface < num_faces; iface++) {
      /* Read the face neighbors tree id and the tree_to_face entry */
      ret = fscanf (fp, "%li %i%*c", &neighbor, &ttf_entry);
      T8_SAVE_CHECK_CLOSE (ret == 2, fp);
      face_neighbors[iface] = neighbor;
      ttf[iface] = ttf_entry;
    }
    vertices = SC_REALLOC (vertices, double, 3 * num_vertices);
    for (att = 0; att < tree->num_attributes; att++) {
      /* Loop over all attributes of this tree that we need to read */
      /* Check whether this attribute really belongs to tree */
      T8_SAVE_CHECK_CLOSE (treeid == (long) itree, fp);
      ret = fscanf (fp, "id %i\nkey %i\n", &att_struct.package_id, &att_struct.key);
      T8_SAVE_CHECK_CLOSE (ret == 2, fp);
      /* We currently only support vertices as attributes.
       * Those have t8 package id and key 0 */
      T8_SAVE_CHECK_CLOSE (att_struct.package_id > 0 && att_struct.key == 0, fp);
      /* TODO: We set the package id to match the one of t8code manually.
       * See the comment above. */
      att_struct.package_id = t8_get_package_id ();

      /* read the size of the attribute */
      ret = fscanf (fp, "size %zu\n", &att_struct.attr_size);
      T8_SAVE_CHECK_CLOSE (ret == 1, fp);
      /* Read the vertices */
      for (i = 0; i < num_vertices; i++) {
        ret = fscanf (fp, "%lf %lf %lf\n", vertices + 3 * i, vertices + 3 * i + 1, vertices + 3 * i + 2);
        T8_SAVE_CHECK_CLOSE (ret == 3, fp);
      }
      att_struct.attr_data = vertices;
      att_struct.is_owned = 0;
      att_struct.id = itree + cmesh->first_tree;
      /* Now we read the attribute and can add it to the tree */
      t8_cmesh_trees_add_attribute (cmesh->trees, 0, &att_struct, itree, 0);
    }
  }
  SC_FREE (vertices);
  return 1;
}

static int
t8_cmesh_save_tree_attribute (const t8_cmesh_t cmesh, FILE *fp)
{
  double *vertices;
  t8_locidx_t itree;
  t8_ctree_t tree;
  int num_vertices;
  int ret, i;
  size_t att_size;
  t8_locidx_t *face_neigh;
  int8_t *ttf;
  int num_faces;
  char buffer[BUFSIZ] = "";

  ret = fprintf (fp, "\n--- Tree attribute section ---\n");
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  /* For each tree, write its attribute */
  for (itree = 0; itree < cmesh->num_local_trees; itree++) {
    /* TODO: Currently we only support storing of the tree vertices as only attribute.
     *        This attribute has package id t8_get_package_id() and key 0.
     */
    tree = t8_cmesh_trees_get_tree_ext (cmesh->trees, itree, &face_neigh, &ttf);
    ret = fprintf (fp, "\ntree %li\n", (long) itree);
    /* Iterate over all faces to write the face neighbor information */
    num_faces = t8_eclass_num_faces[tree->eclass];
    for (i = 0; i < num_faces; i++) {
      snprintf (buffer + strlen (buffer), BUFSIZ - strlen (buffer), "%li %i%s", (long) face_neigh[i], ttf[i],
                i == num_faces - 1 ? "" : ", ");
    }
    ret = fprintf (fp, "Neighbors: %s\n", buffer);
    /* Clear the buffer such that strlen returns 0 */
    buffer[0] = '\0';
    T8_SAVE_CHECK_CLOSE (ret > 0, fp);
    num_vertices = t8_eclass_num_vertices[tree->eclass];
    /* Write the attributes that are vertices */
    vertices = (double *) t8_cmesh_trees_get_attribute (cmesh->trees, itree, t8_get_package_id (), 0, &att_size, 0);
    if (vertices != NULL) {
      /* We have an attribute that is stored with key 0, we treat it as tree vertices */
      num_vertices = t8_eclass_num_vertices[tree->eclass];
      /* additional check that the attribute has as many bytes as the
       * coordinates for this tree */
      if (att_size != num_vertices * 3 * sizeof (double)) {
        /* TODO: We currently do not support saving of attributes different
         * to the tree vertices */
        fclose (fp);
        t8_errorf ("We do not support saving cmeshes with trees that "
                   "have attributes different to the tree vertices.\n");
        return 0;
      }
      ret = fprintf (fp, "id %i\nkey %i\n", t8_get_package_id (), 0);
      T8_SAVE_CHECK_CLOSE (ret > 0, fp);
      ret = fprintf (fp, "size %zd\n", num_vertices * 3 * sizeof (double));
      T8_ASSERT (strlen (buffer) == 0);
      for (i = 0; i < num_vertices; i++) {
        /* For each vertex, we write its three coordinates in a line */
        snprintf (buffer + strlen (buffer), BUFSIZ - strlen (buffer), "%e %e %e\n", vertices[3 * i],
                  vertices[3 * i + 1], vertices[3 * i + 2]);
      }
      fprintf (fp, "%s", buffer);
      /* Clear the buffer such that strlen returns 0 */
      buffer[0] = '\0';
    }
  }
  return 1;
}

static int
t8_cmesh_save_trees (const t8_cmesh_t cmesh, FILE *fp)
{
  t8_locidx_t itree;
  t8_ctree_t tree;
  int ret;

  T8_ASSERT (fp != NULL);
  ret = fprintf (fp, "Total bytes for trees and ghosts %zd\n", t8_cmesh_trees_size (cmesh->trees));
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);

  ret = fprintf (fp, "\n--- Tree section ---");
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  ret = fprintf (fp, "\n\nTrees %li\n\n", (long) cmesh->num_local_trees);
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);

  /* For each tree, write its metadata */
  for (itree = 0; itree < cmesh->num_local_trees; itree++) {
    tree = t8_cmesh_trees_get_tree (cmesh->trees, itree);
    ret = fprintf (fp, "eclass %i\n", (int) tree->eclass);
    T8_SAVE_CHECK_CLOSE (ret > 0, fp);
    if (tree->num_attributes != 1) {
      /* TODO: We currently do not support saving of attributes different
       * to the tree vertices */
      fclose (fp);
      t8_errorf ("We do not support saving cmeshes with trees that have attributes different to the tree vertices.\n");
      return 0;
    }
    ret = fprintf (fp, "num_attributes %i\nSize of attributes %zd\n\n", tree->num_attributes,
                   t8_cmesh_trees_attribute_size (tree));
    T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  }
  return 1;
}

/* Load all tree data (eclasses, neighbors, vertex coordinates,...)
 * from a cmesh file into a cmesh */
static int
t8_cmesh_load_trees (const t8_cmesh_t cmesh, FILE *fp)
{
  size_t bytes_for_trees, att_bytes;
  t8_locidx_t itree;
  int eclass, num_atts;
  int ret;
  long num_trees;

  ret = fscanf (fp, "Total bytes for trees and ghosts %zu\n", &bytes_for_trees);
  T8_SAVE_CHECK_CLOSE (ret == 1, fp); /* The bytes_for_trees data is currently not used */
  ret = fscanf (fp, "\n--- Tree section ---\n");
  T8_SAVE_CHECK_CLOSE (ret == 0, fp);
  t8_cmesh_trees_init (&cmesh->trees, 1, cmesh->num_local_trees, cmesh->num_ghosts);
  t8_cmesh_trees_start_part (cmesh->trees, 0, 0, cmesh->num_local_trees, 0, cmesh->num_ghosts, 1);
  /* Read the number of trees to come and check for consistency */
  ret = fscanf (fp, "Trees %li\n", &num_trees);
  T8_SAVE_CHECK_CLOSE (ret == 1, fp);
  T8_SAVE_CHECK_CLOSE (num_trees == cmesh->num_local_trees, fp);
  /* Read all the trees from the file and add them to the
   * trees structure */
  for (itree = 0; itree < cmesh->num_local_trees; itree++) {
    ret = fscanf (fp, "eclass %i\n", &eclass);
    T8_SAVE_CHECK_CLOSE (ret == 1, fp);
    t8_cmesh_trees_add_tree (cmesh->trees, itree, 0, (t8_eclass_t) eclass);

    /* After adding the tree, we set its face neighbors and face orientation */
    (void) t8_cmesh_trees_get_tree (cmesh->trees, itree);
    /* Check whether the number of attribute is really 1 */
    ret = fscanf (fp, "num_attributes %i\nSize of attributes %zu\n", &num_atts, &att_bytes);
    T8_SAVE_CHECK_CLOSE (ret == 2, fp);
    T8_SAVE_CHECK_CLOSE (num_atts == 1, fp);
    T8_SAVE_CHECK_CLOSE (att_bytes == 3 * sizeof (double) * t8_eclass_num_vertices[eclass], fp);
    /* If there really is one attribute it must be the trees vertices and thus
     * we initialize the tree's attribute accordingly */
    t8_cmesh_trees_init_attributes (cmesh->trees, itree, 1, att_bytes);
  }
  return 1;
}

static int
t8_cmesh_save_header (const t8_cmesh_t cmesh, FILE *fp)
{
  int ret;
  int eclass;

  T8_ASSERT (fp != NULL);
  ret = fprintf (fp, "This is %s, file format version %u.\n\n", t8_get_package_string (), T8_CMESH_FORMAT);
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);

  /* Write 0 for replicated and 1 for partitioned cmesh */
  ret = fprintf (fp, "Partitioned %i\n", cmesh->set_partition != 0);
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);

  /* Write the rank of this process and the total number pf processes */
  ret = fprintf (fp, "Rank %i of %i\n", cmesh->mpirank, cmesh->mpisize);
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);

  /* Write the dimension of the cmesh */
  ret = fprintf (fp, "dim %i\n", cmesh->dimension);
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);

  /* Write the number of global and local trees and the number of ghosts */
  ret = fprintf (fp, "num_trees %lli\n", (long long) cmesh->num_trees);
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  ret = fprintf (fp, "num_local_trees %li\n", (long) cmesh->num_local_trees);
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  ret = fprintf (fp, "num_ghosts %li\n", (long) cmesh->num_ghosts);
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);

  /* Write the number of local trees for each eclass */
  ret = fprintf (fp, "num_local_trees_per_eclass ");
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  for (eclass = T8_ECLASS_ZERO; eclass < T8_ECLASS_COUNT; eclass++) {
    ret = fprintf (fp, "%li%s", (long) cmesh->num_local_trees_per_eclass[eclass],
                   eclass == T8_ECLASS_COUNT - 1 ? "\n" : ", ");
    T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  }
  /* Write the number of trees for each eclass */
  ret = fprintf (fp, "num_trees_per_eclass ");
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  for (eclass = T8_ECLASS_ZERO; eclass < T8_ECLASS_COUNT; eclass++) {
    ret = fprintf (fp, "%lli%s", (long long) cmesh->num_trees_per_eclass[eclass],
                   eclass == T8_ECLASS_COUNT - 1 ? "\n" : ", ");
    T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  }

  /* Write the global id of the first local tree and if the first tree is shared */
  ret = fprintf (fp, "first_tree %lli\n", (long long) cmesh->first_tree);
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  ret = fprintf (fp, "first_tree_shared %i\n", cmesh->first_tree_shared);
  T8_SAVE_CHECK_CLOSE (ret > 0, fp);
  return 1;
}

/* Read the number of trees, dimension, etc. from a saved cmesh file.
 * If anything goes wrong, the file is closed and 0 is returned */
static int
t8_cmesh_load_header (const t8_cmesh_t cmesh, FILE *fp)
{
  int file_format, save_rank, save_mpisize;
  int ieclass;
  int ret;
  long long global_num_trees, tree_per_class, first_local_tree;
  long local_num_trees, local_num_ghosts;
  int first_shared;

  /* Check whether the file was saved with the current file format */
  ret = fscanf (fp, "This is %*[^ ] %*[^ ] file format version %i.\n", &file_format);
  T8_SAVE_CHECK_CLOSE (ret == 1, fp);
  if (file_format != T8_CMESH_FORMAT) {
    /* The file was saved with an old format and we cannot read it any more */
    t8_errorf ("Input file is in an old format that we cannot read anymore.\n");
    fclose (fp);
    return 0;
  }

  /* Read whether the cmesh is partitioned, its dimension and the number of
   * trees. Also read which MPI rank wrote the file */
  ret = fscanf (fp, "Partitioned %i\n", &cmesh->set_partition);
  T8_SAVE_CHECK_CLOSE (ret == 1, fp);
  ret = fscanf (fp, "Rank %i of %i\n", &save_rank, &save_mpisize);
  T8_SAVE_CHECK_CLOSE (ret == 2, fp);
  /* Check if the rank and mpisize stored were valid */
  T8_SAVE_CHECK_CLOSE (0 <= save_rank && save_rank < save_mpisize, fp);
  /* It does not make sense to load a cmesh on a rank smaller than the one that
   * saved it. */
  T8_SAVE_CHECK_CLOSE (cmesh->mpirank <= save_rank && cmesh->mpisize <= save_mpisize, fp);
  ret = fscanf (fp, "dim %i\n", &cmesh->dimension);
  T8_SAVE_CHECK_CLOSE (ret == 1, fp);
  /* Check if the read dimension is in the correct range */
  T8_SAVE_CHECK_CLOSE (cmesh->dimension >= 0 && cmesh->dimension <= 3, fp);
  /* Since t8_gloidx_t and t8_locidx_t are integer datatypes that are not fixed,
   * we first read the tree numbers into fixed datatypes (long long for gloidx and
   * long for locidx) */
  ret = fscanf (fp, "num_trees %lli\nnum_local_trees %li\n", &global_num_trees, &local_num_trees);
  T8_SAVE_CHECK_CLOSE (ret == 2, fp);
  T8_SAVE_CHECK_CLOSE (local_num_trees <= global_num_trees, fp);
  cmesh->num_trees = (t8_gloidx_t) global_num_trees;
  cmesh->num_local_trees = (t8_locidx_t) local_num_trees;
  /* Read the number of ghost trees */
  ret = fscanf (fp, "num_ghosts %li\n", &local_num_ghosts);
  T8_SAVE_CHECK_CLOSE (ret == 1, fp);
  T8_SAVE_CHECK_CLOSE (0 <= local_num_ghosts && local_num_ghosts < cmesh->num_trees, fp);
  cmesh->num_ghosts = local_num_ghosts;

  /* Read the number of local trees per eclass */
  ret = fscanf (fp, "num_local_trees_per_eclass");
  T8_SAVE_CHECK_CLOSE (ret == 0, fp);
  for (ieclass = T8_ECLASS_ZERO; ieclass < T8_ECLASS_COUNT; ieclass++) {
    ret = fscanf (fp, "%lli,", &tree_per_class);
    T8_SAVE_CHECK_CLOSE (ret == 1, fp);
    cmesh->num_local_trees_per_eclass[ieclass] = (t8_locidx_t) tree_per_class;
  }
  /* Read the number of trees per eclass */
  ret = fscanf (fp, "\nnum_trees_per_eclass");
  T8_SAVE_CHECK_CLOSE (ret == 0, fp);
  for (ieclass = T8_ECLASS_ZERO; ieclass < T8_ECLASS_COUNT; ieclass++) {
    ret = fscanf (fp, "%lli,", &tree_per_class);
    T8_SAVE_CHECK_CLOSE (ret == 1, fp);
    cmesh->num_trees_per_eclass[ieclass] = (t8_gloidx_t) tree_per_class;
  }
  /* Read the first local tree id and whether the first tree is shared */
  ret = fscanf (fp, "\nfirst_tree %lli\n", &first_local_tree);
  T8_SAVE_CHECK_CLOSE (ret == 1, fp);
  T8_SAVE_CHECK_CLOSE (0 <= first_local_tree && first_local_tree < global_num_trees, fp);
  cmesh->first_tree = first_local_tree;
  ret = fscanf (fp, "first_tree_shared %i\n", &first_shared);
  T8_SAVE_CHECK_CLOSE (ret == 1, fp);
  cmesh->first_tree_shared = first_shared;
  T8_SAVE_CHECK_CLOSE (!cmesh->set_partition || cmesh->first_tree_shared == 0 || cmesh->first_tree_shared == 1, fp);
  return 1;
}

int
t8_cmesh_save (const t8_cmesh_t cmesh, const char *fileprefix)
{
  FILE *fp;
  char filename[BUFSIZ];
  int has_linear_geom = 0;

  T8_ASSERT (t8_cmesh_is_committed (cmesh));
  if (!cmesh->set_partition && cmesh->mpirank != 0) {
    /* If the cmesh is replicated, only rank 0 writes it */
    return 1;
  }

  /* Check that the only registered geometry is the linear geometry and
   * that this geometry is used for all trees. */
  if (cmesh->geometry_handler->get_num_geometries () == 1) {
    /* Get the stored geometry and the linear geometry and compare their names. */
    const t8_geometry *geom = cmesh->geometry_handler->get_unique_geometry ();
    has_linear_geom = geom->t8_geom_get_type () == T8_GEOMETRY_TYPE_LINEAR;
  }
  if (!has_linear_geom) {
    /* This cmesh does not have the linear geometry for all trees. */
    t8_errorf ("Error when saving cmesh. Cmesh has more than one geometry or the geometry is not linear.\n");
    return 0;
  }

  /* Create the output filename as fileprefix_RANK.cmesh,
   * where we write RANK with 4 significant figures. */
  snprintf (filename, BUFSIZ, "%s_%04i.cmesh", fileprefix, cmesh->mpirank);

  /* Open the file in write mode */
  fp = fopen (filename, "w");
  if (fp == NULL) {
    /* Could not open file */
    t8_errorf ("Error when opening file %s.\n", filename);
    return 0;
  }
  /* Write all metadata of the cmesh */
  if (!t8_cmesh_save_header (cmesh, fp)) {
    t8_errorf ("Error when opening file %s.\n", filename);
    return 0;
  }
  /* Write all metadata of the trees */
  if (!t8_cmesh_save_trees (cmesh, fp)) {
    t8_errorf ("Error when opening file %s.\n", filename);
    return 0;
  }
  if (cmesh->set_partition) {
    /* Write all ghost metadata */
    if (!t8_cmesh_save_ghosts (cmesh, fp)) {
      t8_errorf ("Error when opening file %s.\n", filename);
      return 0;
    }
  }
  if (!t8_cmesh_save_tree_attribute (cmesh, fp)) {
    /* Write all tree attributes */
    t8_errorf ("Error when opening file %s.\n", filename);
    return 0;
  }
  if (cmesh->set_partition) {
    /* Write all ghost neighbor data */
    if (!t8_cmesh_save_ghost_neighbors (cmesh, fp)) {
      t8_errorf ("Error when opening file %s.\n", filename);
      return 0;
    }
  }
  /* Close the file */
  fclose (fp);
  return 1;
}

#undef T8_SAVE_CHECK_CLOSE

t8_cmesh_t
t8_cmesh_load (const char *filename, sc_MPI_Comm comm)
{
  FILE *fp;
  t8_cmesh_t cmesh;
  int mpiret;

  /* Open the file in read mode */
  fp = fopen (filename, "r");
  if (fp == NULL) {
    /* Could not open file */
    t8_errorf ("Error when opening file %s.\n", filename);
    return NULL;
  }
  t8_cmesh_init (&cmesh);
  /* Read all metadata of the cmesh */
  if (!t8_cmesh_load_header (cmesh, fp)) {
    t8_errorf ("Error when opening file %s.\n", filename);
    t8_cmesh_destroy (&cmesh);
    return NULL;
  }
  /* Read all metadata of the trees */
  if (!t8_cmesh_load_trees (cmesh, fp)) {
    t8_errorf ("Error when opening file %s.\n", filename);
    t8_cmesh_destroy (&cmesh);
    return NULL;
  }
  if (cmesh->set_partition) {
    /* Write all ghost metadata */
    if (!t8_cmesh_load_ghosts (cmesh, fp)) {
      t8_errorf ("Error when opening file %s.\n", filename);
      t8_cmesh_destroy (&cmesh);
      return NULL;
    }
  }
  t8_cmesh_trees_finish_part (cmesh->trees, 0);
  if (!t8_cmesh_load_tree_attributes (cmesh, fp)) {
    t8_errorf ("Error when opening file %s.\n", filename);
    t8_cmesh_destroy (&cmesh);
    return NULL;
  }
  if (cmesh->set_partition) {
    /* Write all ghost metadata */
    if (!t8_cmesh_load_ghost_attributes (cmesh, fp)) {
      t8_errorf ("Error when opening file %s.\n", filename);
      t8_cmesh_destroy (&cmesh);
      return NULL;
    }
  }
  /* Close the file */
  fclose (fp);

  cmesh->committed = 1;
  mpiret = sc_MPI_Comm_rank (comm, &cmesh->mpirank);
  SC_CHECK_MPI (mpiret);
  mpiret = sc_MPI_Comm_size (comm, &cmesh->mpisize);
  SC_CHECK_MPI (mpiret);
  t8_stash_destroy (&cmesh->stash);
  return cmesh;
}

/* Query whether a given process will open a cmesh saved file.
 * This depends on the number of processes, the number of files and
 * the load mode.
 * In load mode simple, if rank i is smaller than the number of files,
 * rank i opens file number i.
 * In load mode bgq, there must be at leas num_files compute nodes
 * and on the i-th node exactly one process opens file i.
 * This function returns true, if the given rank will load a file and
 * it stores the number of the file in file_to_load.
 */
static int
t8_cmesh_load_proc_loads (const int mpirank, const int mpisize, const int num_files, sc_MPI_Comm comm,
                          const t8_load_mode_t mode, int *file_to_load, int num_procs_per_node)
{
  sc_MPI_Comm inter = sc_MPI_COMM_NULL, intra = sc_MPI_COMM_NULL;
  int mpiret, interrank, intrarank, intersize;

  if (num_procs_per_node <= 0 && mode == T8_LOAD_STRIDE) {
    t8_global_infof ("number of processes per node set to 16\n");
    num_procs_per_node = 16;
  }
  /* Fill with invalid value */
  *file_to_load = -1;
  switch (mode) {
  case T8_LOAD_SIMPLE:
    /* In simple mode, process i opens file i, if i < num_files */
    if (mpirank < num_files) {
      /* Open file mpirank */
      *file_to_load = mpirank;
      return 1;
    }
    else {
      /* Do not open any files */
      return 0;
    }
    break;
  case T8_LOAD_BGQ:
    /* In bgq mode on each shared memory node is one process that
     * opens a file. */
    /* Compute internode and intranode communicator */
    sc_mpi_comm_attach_node_comms (comm, 0);
    /* Store internode and intranode communicator */
    sc_mpi_comm_get_node_comms (comm, &intra, &inter);
    /* Abort if we could not compute these communicators */
    SC_CHECK_ABORT (intra != sc_MPI_COMM_NULL && inter != sc_MPI_COMM_NULL, "Could not get proper internode "
                                                                            "and intranode communicators.\n");
    mpiret = sc_MPI_Comm_size (inter, &intersize);
    SC_CHECK_MPI (mpiret);
    /* Check if the number of nodes is at least as big as the number of files */
    SC_CHECK_ABORTF (intersize <= num_files,
                     "Must have more compute nodes "
                     "than files. %i nodes and %i fields are given.\n",
                     intersize, num_files);
    mpiret = sc_MPI_Comm_rank (inter, &interrank);
    SC_CHECK_MPI (mpiret);
    mpiret = sc_MPI_Comm_rank (intra, &intrarank);
    SC_CHECK_MPI (mpiret);
    /* If the current node number i is smaller than the number of files
     * than the first process on this node loads the file i */
    if (interrank < num_files && intrarank == 0) {
      *file_to_load = interrank;
      return 1;
    }
    else {
      return 0;
    }
    break;
  case T8_LOAD_STRIDE:
    /* In Juqueen mode, every 16-th process loads a file. The user should
     * control, that these processes reside on different compute nodes to
     * gain maximal efficiency. */
    SC_CHECK_ABORT (ceil (mpisize / (double) num_procs_per_node) >= num_files,
                    "Too many files for too few processes.\n");
    if (mpirank % num_procs_per_node == 0 && mpirank / num_procs_per_node < num_files) {
      *file_to_load = mpirank / num_procs_per_node;
      return 1;
    }
    else {
      return 0;
    }
  default:
    SC_ABORT_NOT_REACHED ();
  }
  SC_ABORT_NOT_REACHED ();
}

/* After we loaded a cmesh on a part of the processes, we have to correctly
 * set the first tree on the other, empty, cmeshes.
 * In load mode simple, the first tree is just the global number of trees.
 * In Load mode BGQ and JUQUEEN, we have to figure out the first tree of
 * the next bigger nonloading process.
 */
static int
t8_cmesh_load_bigger_nonloading (const int mpirank, const int mpisize, const int num_files, const t8_load_mode_t mode,
                                 sc_MPI_Comm comm, const int num_procs_per_node)
{
  int next_bigger_nonloading;
  sc_MPI_Comm inter = sc_MPI_COMM_NULL, intra = sc_MPI_COMM_NULL;
  sc_MPI_Group intragroup, commgroup;
  int mpiret, interrank, intrarank, intrasize, commrank;
  int rankzero;

  switch (mode) {
  case T8_LOAD_SIMPLE:
    /* In simple mode, the first num_files processes load the cmesh and
     * the rest is empty, this the next bigger nonloading rank is always
     * rank mpisize. */
    next_bigger_nonloading = mpisize;
    break;
  case T8_LOAD_BGQ:
    /* In BGQ mode, on the compute nodes with node-id bigger num_files,
     * the first rank on a node opens the file. */
    /* Get the inter and intra comms and our rank within these comms. */
    sc_mpi_comm_get_node_comms (comm, &intra, &inter);
    mpiret = sc_MPI_Comm_rank (inter, &interrank);
    SC_CHECK_MPI (mpiret);
    mpiret = sc_MPI_Comm_rank (intra, &intrarank);
    SC_CHECK_MPI (mpiret);

    if (interrank >= num_files - 1) {
      /* We are beyond all processes that loaded */
      next_bigger_nonloading = mpisize;
    }
    else {
      /* We need to get the rank in comm of the first process
       * on the next compute node. This is the process with rank 0 in our
       * intra communicator plus the mpisize in the intracommunicator. */
      /* First we get the groups from the communicators */
      mpiret = sc_MPI_Comm_group (intra, &intragroup);
      SC_CHECK_MPI (mpiret);
      mpiret = sc_MPI_Comm_group (comm, &commgroup);
      SC_CHECK_MPI (mpiret);
      /* We can now transform the rank 0 of the intragroup to get
       * the correct rank in the commgroup */
      rankzero = 0;
      mpiret = sc_MPI_Group_translate_ranks (intragroup, 1, &rankzero, commgroup, &commrank);
      SC_CHECK_MPI (mpiret);
      /* We get the size of the intracommunicator */
      mpiret = sc_MPI_Group_size (intragroup, &intrasize);
      SC_CHECK_MPI (mpiret);
      next_bigger_nonloading = commrank + intrasize;
    }
    break;
  case T8_LOAD_STRIDE:
    /* In Juqueen mode, every 16-th process has opened the file. */
    if (mpirank / num_procs_per_node < num_files - 1) {
      /* If we are in a multiple of 16, where a file was loaded,
       * the first process in the next group did load it. */
      next_bigger_nonloading
        = mpirank - mpirank % num_procs_per_node + num_procs_per_node; /* This is the next number divisible by
                                                                                                   num_procs_per_node */
    }
    else {
      /* Every rank that loaded a file is smaller than ours. */
      next_bigger_nonloading = mpisize;
    }
    break;
  default:
    SC_ABORT_NOT_REACHED ();
  }
  return next_bigger_nonloading;
}

/* Load the files fileprefix_0000.cmesh, ... , fileprefix_N.cmesh
 * (N = num_files - 1)
 * on N processes and repartition the cmesh to all calling processes.
 * If N = 1, the cmesh is broadcasted and not partitioned.
 */
t8_cmesh_t
t8_cmesh_load_and_distribute (const char *fileprefix, const int num_files, sc_MPI_Comm comm, const t8_load_mode_t mode,
                              const int procs_per_node)
{
  t8_cmesh_t cmesh;
  char buffer[BUFSIZ];
  int mpiret, mpirank, mpisize;
  int file_to_load;
  int next_bigger_nonloading;
  int did_load;

  mpiret = sc_MPI_Comm_rank (comm, &mpirank);
  SC_CHECK_MPI (mpiret);
  mpiret = sc_MPI_Comm_size (comm, &mpisize);
  SC_CHECK_MPI (mpiret);

  T8_ASSERT (mpisize >= num_files);

  /* Try to set the comm type */
  SC_CHECK_ABORT (t8_shmem_init (comm) > 0, "Error in shared memory setup. Could not load cmesh.");
  t8_shmem_set_type (comm, T8_SHMEM_BEST_TYPE);

  /* Use cmesh_bcast, if only one process loads the cmesh: */
  if (num_files == 1) {
    cmesh = NULL;
    if (mpirank == 0) {
      snprintf (buffer, BUFSIZ, "%s_%04d.cmesh", fileprefix, 0);
      cmesh = t8_cmesh_load (buffer, comm);
    }
    cmesh = t8_cmesh_bcast (cmesh, 0, comm);
  }
  else {
    /* More than one process loads a file */

    if (t8_cmesh_load_proc_loads (mpirank, mpisize, num_files, comm, mode, &file_to_load, procs_per_node)) {
      T8_ASSERT (fileprefix != NULL);
      T8_ASSERT (0 <= file_to_load && file_to_load < num_files);
      snprintf (buffer, BUFSIZ, "%s_%04d.cmesh", fileprefix, file_to_load);
      t8_infof ("Opening file %s\n", buffer);
      cmesh = t8_cmesh_load (buffer, comm);
      if (num_files == mpisize) {
        /* Each process has loaded the cmesh and we can return */
        return cmesh;
      }
      did_load = 1;
    }
    else {
      did_load = 0;
      /* On the processes that do not load the cmesh, initialize it
       * with zero local trees and ghosts */
      t8_cmesh_init (&cmesh);
      t8_cmesh_trees_init (&cmesh->trees, 0, 0, 0);
      cmesh->mpirank = mpirank;
      cmesh->mpisize = mpisize;
      t8_stash_destroy (&cmesh->stash);
      /* There are no faces, so we know all about them */
      cmesh->face_knowledge = 3;
      /* There are no tree, thus the first tree is not shared */
      cmesh->first_tree_shared = 0;
      cmesh->committed = 1;
      cmesh->set_partition = 1;
      cmesh->num_local_trees = 0;
    }
    /* The cmeshes on the processes that did not load have to
     * know the global number of trees */
    sc_MPI_Bcast (&cmesh->num_trees, 1, T8_MPI_GLOIDX, 0, comm);
    /* And the dimension */
    sc_MPI_Bcast (&cmesh->dimension, 1, sc_MPI_INT, 0, comm);
    t8_cmesh_gather_trees_per_eclass (cmesh, comm);
    T8_ASSERT (t8_cmesh_is_committed (cmesh));
    /* We now create the cmeshs offset in order to properly
     * set the first tree for the empty processes */
    t8_cmesh_gather_treecount (cmesh, comm);
    if (!did_load) {
      /* Calculate the next bigger nonloading rank. */
      next_bigger_nonloading
        = t8_cmesh_load_bigger_nonloading (mpirank, mpisize, num_files, mode, comm, procs_per_node);
      /* Set the first tree of this process to the first tree of the next nonloading
       * rank */
      cmesh->first_tree
        = t8_offset_first (next_bigger_nonloading, t8_shmem_array_get_gloidx_array (cmesh->tree_offsets));
    }
    /* Since we changed the first tree on some processes, we have to
     * regather the first trees on each process */
    t8_cmesh_gather_treecount (cmesh, comm);
  }
  T8_ASSERT (t8_cmesh_is_committed (cmesh));
  return cmesh;
}