Refined Pyramid Face Normal Consistency Bug: inverted pyramid child (sub-element ID 8)

[ackirby88]

Refined pyramid sub-element number 8 has incorrect/inconsistent face normal values.
This program reads a gmsh file with a single pyramid.
Below shows the test code, gmsh file, and output of the program.

Execution of the test code (we modified the step3 program):
mpirun -np 1 ./t8_step3_adapt_forest box_pyr1_mesh

Program Output:
(Note the normals for, e.g., elem 8, neigh 6 are not the same (w/in sign) of elem 6, neigh 8):

[t8 0] elem:0 neigh:1  centroid_diff:  0.250000  0.000000  0.150000  normal:  0.894427 -0.000000  0.447214  dot:  0.997054
[t8 0] elem:0 neigh:3  centroid_diff:  0.000000  0.250000  0.150000  normal: -0.000000  0.894427  0.447214  dot:  0.997054
[t8 0] elem:1 neigh:8  centroid_diff:  0.000000  0.250000  0.150000  normal: -0.000000  0.894427  0.447214  dot:  0.997054
[t8 0] elem:1 neigh:2  centroid_diff:  0.250000  0.000000 -0.150000  normal:  0.894427 -0.000000 -0.447214  dot:  0.997054
[t8 0] elem:1 neigh:0  centroid_diff: -0.250000  0.000000 -0.150000  normal: -0.894427  0.000000 -0.447214  dot:  0.997054
[t8 0] elem:2 neigh:1  centroid_diff: -0.250000  0.000000  0.150000  normal: -0.894427 -0.000000  0.447214  dot:  0.997054
[t8 0] elem:2 neigh:5  centroid_diff:  0.000000  0.250000  0.150000  normal: -0.000000  0.894427  0.447214  dot:  0.997054
[t8 0] elem:3 neigh:8  centroid_diff:  0.250000  0.000000  0.150000  normal:  0.894427  0.000000  0.447214  dot:  0.997054
[t8 0] elem:3 neigh:4  centroid_diff:  0.000000  0.250000 -0.150000  normal:  0.000000  0.894427 -0.447214  dot:  0.997054
[t8 0] elem:3 neigh:0  centroid_diff:  0.000000 -0.250000 -0.150000  normal:  0.000000 -0.894427 -0.447214  dot:  0.997054
[t8 0] elem:4 neigh:6  centroid_diff:  0.250000  0.000000  0.150000  normal:  0.894427 -0.000000  0.447214  dot:  0.997054
[t8 0] elem:4 neigh:3  centroid_diff:  0.000000 -0.250000  0.150000  normal:  0.000000 -0.894427  0.447214  dot:  0.997054
[t8 0] elem:5 neigh:8  centroid_diff: -0.250000  0.000000  0.150000  normal: -0.894427  0.000000  0.447214  dot:  0.997054
[t8 0] elem:5 neigh:7  centroid_diff:  0.000000  0.250000 -0.150000  normal:  0.000000  0.894427 -0.447214  dot:  0.997054
[t8 0] elem:5 neigh:2  centroid_diff:  0.000000 -0.250000 -0.150000  normal:  0.000000 -0.894427 -0.447214  dot:  0.997054
[t8 0] elem:6 neigh:8  centroid_diff:  0.000000 -0.250000  0.150000  normal: -0.000000 -0.894427  0.447214  dot:  0.997054
[t8 0] elem:6 neigh:7  centroid_diff:  0.250000  0.000000 -0.150000  normal:  0.894427 -0.000000 -0.447214  dot:  0.997054
[t8 0] elem:6 neigh:4  centroid_diff: -0.250000  0.000000 -0.150000  normal: -0.894427  0.000000 -0.447214  dot:  0.997054
[t8 0] elem:7 neigh:6  centroid_diff: -0.250000  0.000000  0.150000  normal: -0.894427 -0.000000  0.447214  dot:  0.997054
[t8 0] elem:7 neigh:5  centroid_diff:  0.000000 -0.250000  0.150000  normal:  0.000000 -0.894427  0.447214  dot:  0.997054
[t8 0] elem:8 neigh:6  centroid_diff:  0.000000  0.250000 -0.150000  normal: -0.894427  0.000000 -0.447214  dot:  0.230089
[t8 0] elem:8 neigh:1  centroid_diff:  0.000000 -0.250000 -0.150000  normal:  0.894427 -0.000000 -0.447214  dot:  0.230089
[t8 0] elem:8 neigh:5  centroid_diff:  0.250000  0.000000 -0.150000  normal:  0.000000 -0.894427 -0.447214  dot:  0.230089
[t8 0] elem:8 neigh:3  centroid_diff: -0.250000  0.000000 -0.150000  normal:  0.000000  0.894427 -0.447214  dot:  0.230089
[t8 0] elem:8 neigh:9  centroid_diff:  0.000000  0.000000  0.200000  normal:  0.000000  0.000000  1.000000  dot:  1.000000
[t8 0] elem:9 neigh:8  centroid_diff:  0.000000  0.000000 -0.200000  normal: -0.000000 -0.000000 -1.000000  dot:  1.000000

gmsh file (name box_pyr1_mesh.msh):

$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
5
1 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00
2 0.1000000000000000E+01 0.0000000000000000E+00 0.0000000000000000E+00
3 0.1000000000000000E+01 0.1000000000000000E+01 0.0000000000000000E+00
4 0.0000000000000000E+00 0.1000000000000000E+01 0.0000000000000000E+00
5 0.5000000000000000E+00 0.5000000000000000E+00 0.1000000000000000E+01
$EndNodes
$Elements
6
1 2 2   2  2 1 5 2
2 2 2   3  3 2 5 3
3 2 2   4  4 3 5 4
4 2 2   5  5 4 5 1
5 3 2   1  1 1 2 3 4
6 7 2 6 6 1 2 3 4 5
$EndElements
$PhysicalNames
6
2 1 "FACE-1"
2 2 "FACE-2"
2 3 "FACE-3"
2 4 "FACE-4"
2 5 "FACE-5"
2 6 "VOL1"
$EndPhysicalNames

Test code (following step3):

/*
  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.
*/

/* This is step3 of the t8code tutorials.
 * After generating a coarse mesh (step1) and building a uniform forest
 * on it (step2), we will now adapt (= refine and coarsen) the forest
 * according to our own criterion.
 *
 * The geometry (coarse mesh) is again a cube, this time modelled with
 * 6 tetrahedra, 6 prisms and 4 cubes.
 * We refine an element if its midpoint is within a sphere of given radius
 * around the point (0.5, 0.5, 1) and we coarsen outside of a given radius.
 * We will use non-recursive refinement, that means that the refinement level
 * of any element will change by at most +-1.
 *
 * How you can experiment here:
 *   - Look at the paraview output files of the unifomr and the adapted forest.
 *     For the adapted forest you can apply a slice filter to look into the cube.
 *   - Run the program with different process numbers. You should see that refining is
 *     independent of the number of processes, but coarsening is not.
 *     This is due to the face that a family can only be coarsened if it is completely
 *     local to a single process and the distribution among the process may break this property.
 *   - Change the midpoint coordinates and the radii.
 *   - Change the adaptation criterion such that elements inside the sphere are coarsened
 *     and elements outside are refined.
 *   - Use t8_productionf to print the local number of elements on each process.
 *     Notice, that the uniform forest is evenly distributed, but that the adapted forest
 *     is not. This is due to the fact that we do not repartition our forest here.
 *   - Add a maximum refinement level to the adapt_data struct and use non-recursive refinement.
 *     Do not refine an element if it has reached the maximum level. (Hint: ts->t8_element_level)
 */

#include <t8.h>                                     /* General t8code header, always include this. */
#include <t8_cmesh.h>                               /* cmesh definition and basic interface. */
#include <t8_cmesh/t8_cmesh_examples.h>             /* A collection of exemplary cmeshes */
#include <t8_forest/t8_forest_general.h>            /* forest definition and basic interface. */
#include <t8_forest/t8_forest_io.h>                 /* save forest */
#include <t8_forest/t8_forest_geometrical.h>        /* geometrical information of the forest */
#include <t8_schemes/t8_default/t8_default_cxx.hxx> /* default refinement scheme. */
#include <t8_vec.h>                                 /* Basic operations on 3D vectors. */
#include <tutorials/general/t8_step3.h>
#include <t8_cmesh_readmshfile.h>
#include <t8_cmesh/t8_cmesh_trees.h>

T8_EXTERN_C_BEGIN ();

/* This is our own defined data that we will pass on to the
 * adaptation callback. */

/* The adaptation callback function. This function will be called once for each element
 * and the return value decides whether this element should be refined or not.
 *   return > 0 -> This element should get refined.
 *   return = 0 -> This element should not get refined.
 * If the current element is the first element of a family (= all level l elements that arise from refining
 * the same level l-1 element) then this function is called with the whole family of elements
 * as input and the return value additionally decides whether the whole family should get coarsened.
 *   return > 0 -> The first element should get refined.
 *   return = 0 -> The first element should not get refined.
 *   return < 0 -> The whole family should get coarsened.
 *
 * \param [in] forest       The current forest that is in construction.
 * \param [in] forest_from  The forest from which we adapt the current forest (in our case, the uniform forest)
 * \param [in] which_tree   The process local id of the current tree.
 * \param [in] lelement_id  The tree local index of the current element (or the first of the family).
 * \param [in] ts           The refinement scheme for this tree's element class.
 * \param [in] is_family    if 1, the first \a num_elements entries in \a elements form a family. If 0, they do not.
 * \param [in] num_elements The number of entries in \a elements elements that are defined.
 * \param [in] elements     The element or family of elements to consider for refinement/coarsening.
 */
int
t8_step3_adapt_callback (t8_forest_t forest, t8_forest_t forest_from, t8_locidx_t which_tree, t8_locidx_t lelement_id,
                         t8_eclass_scheme_c *ts, const int is_family, const int num_elements, t8_element_t *elements[])
{
  /* Our adaptation criterion is to look at the midpoint coordinates of the current element and if
   * they are inside a sphere around a given midpoint we refine, if they are outside, we coarsen. */
  double centroid[3]; /* Will hold the element midpoint. */
  /* In t8_step3_adapt_forest we pass a t8_step3_adapt_data pointer as user data to the
   * t8_forest_new_adapt function. This pointer is stored as the used data of the new forest
   * and we can now access it with t8_forest_get_user_data (forest). */
  double refine_r = 1.0;

  /* Compute the element's centroid coordinates. */
  t8_forest_element_centroid (forest_from, which_tree, elements[0], centroid);

  double dl = sqrt(centroid[0]*centroid[0] + centroid[2]*centroid[2]);

  if(dl < refine_r) return 1;

  double  zi = (centroid[0]-0)/(5-0)*(3-0.5)+0.5;
  if (centroid[0] > 0.0 && centroid[0] < 5 &&
    centroid[2] > -zi && centroid[2] < zi) {
    return 1;
  }

  /* Do not change this element. */
  return 0;
}

t8_forest_t
t8_step3_adapt_forest (t8_forest_t forest)
{
  t8_forest_t balanced_forest;

  /* Check that forest is a committed, that is valid and usable, forest. */
  T8_ASSERT (t8_forest_is_committed (forest));

  /* Create a new forest that is adapted from \a forest with our adaptation callback.
   * We provide the adapt_data as user data that is stored as the used_data pointer of the
   * new forest (see also t8_forest_set_user_data).
   * The 0, 0 arguments are flags that control
   *   recursive  -    If non-zero adaptation is recursive, thus if an element is adapted the children
   *                   or parents are plugged into the callback again recursively until the forest does not
   *                   change any more. If you use this you should ensure that refinement will stop eventually.
   *                   One way is to check the element's level against a given maximum level.
   *   do_face_ghost - If non-zero additionally a layer of ghost elements is created for the forest.
   *                   We will discuss ghost in later steps of the tutorial.
   */

  t8_forest_ref(forest);
  t8_forest_init(&balanced_forest);
  t8_forest_set_adapt(balanced_forest, forest, t8_step3_adapt_callback, 0);
  // t8_forest_set_partition (balanced_forest, forest, 0);
  t8_forest_set_balance(balanced_forest, forest, 1);
  t8_forest_set_ghost(balanced_forest, 1, T8_GHOST_FACES);
  t8_forest_commit(balanced_forest);

  t8_forest_unref(&forest);
  return balanced_forest;
}

/* Print the local and global number of elements of a forest. */
void
t8_step3_print_forest_information (t8_forest_t forest)
{
  t8_locidx_t local_num_elements;
  t8_gloidx_t global_num_elements;

  /* Check that forest is a committed, that is valid and usable, forest. */
  T8_ASSERT (t8_forest_is_committed (forest));

  /* Get the local number of elements. */
  local_num_elements = t8_forest_get_local_num_elements (forest);
  /* Get the global number of elements. */
  global_num_elements = t8_forest_get_global_num_elements (forest);
  t8_global_productionf (" [step3] Local number of elements:\t\t%i\n", local_num_elements);
  t8_global_productionf (" [step3] Global number of elements:\t%li\n", global_num_elements);
}

static t8_forest_t
t8_step3_partition_ghost (t8_forest_t forest)
{
  t8_forest_t new_forest;

  /* Check that forest is a committed, that is valid and usable, forest. */
  T8_ASSERT (t8_forest_is_committed (forest));

  /* Initialize */
  t8_forest_init (&new_forest);

  /* Tell the new_forest that is should partition the existing forest.
   * This will change the distribution of the forest elements among the processes
   * in such a way that afterwards each process has the same number of elements
   * (+- 1 if the number of elements is not divisible by the number of processes).
   *
   * The third 0 argument is the flag 'partition_for_coarsening' which is currently not
   * implemented. Once it is, this will ensure that a family of elements will not be split
   * across multiple processes and thus one level coarsening is always possible (see also the
   * comments on coarsening in t8_step3).
   */
  t8_forest_set_partition (new_forest, forest, 0);
  /* Tell the new_forest to create a ghost layer.
   * This will gather those face neighbor elements of process local element that reside
   * on a different process.
   *
   * We currently support ghost mode T8_GHOST_FACES that creates face neighbor ghost elements
   * and will in future also support other modes for edge/vertex neighbor ghost elements.
   */
  t8_forest_set_ghost (new_forest, 1, T8_GHOST_FACES);
  /* Commit the forest, this step will perform the partitioning and ghost layer creation. */
  t8_forest_commit (new_forest);

  return new_forest;
}

static void
t8_step3_exchange_ghost_data (t8_forest_t forest, double *data, int nvar)
{
  sc_array *sc_array_wrapper;
  t8_locidx_t num_elements = t8_forest_get_local_num_elements (forest);
  t8_locidx_t num_ghosts = t8_forest_get_num_ghosts (forest);

  /* t8_forest_ghost_exchange_data expects an sc_array (of length num_local_elements + num_ghosts).
   * We wrap our data array to an sc_array. */
  sc_array_wrapper = sc_array_new_data (data, nvar*sizeof (double), num_elements + num_ghosts);

  /* Carry out the data exchange. The entries with indices > num_local_elements will get overwritten. */
  t8_forest_ghost_exchange_data (forest, sc_array_wrapper);

  /* Destroy the wrapper array. This will not free the data memory since we used sc_array_new_data. */
  sc_array_destroy (sc_array_wrapper);
}

double *
t8_step3_create_element_data (t8_forest_t forest)
{
  t8_locidx_t num_local_elements;
  t8_locidx_t num_ghost_elements;
  double *element_data;

  /* Check that forest is a committed, that is valid and usable, forest. */
  T8_ASSERT (t8_forest_is_committed (forest));

  /* Get the number of local elements of forest. */
  num_local_elements = t8_forest_get_local_num_elements (forest);
  /* Get the number of ghost elements of forest. */
  num_ghost_elements = t8_forest_get_num_ghosts (forest);

  /* Now we need to build an array of our data that is as long as the number
   * of elements plus the number of ghosts. You can use any allocator such as
   * new, malloc or the t8code provide allocation macro T8_ALLOC.
   * Note that in the latter case you need
   * to use T8_FREE in order to free the memory.
   */
  element_data = T8_ALLOC (double, 3*(num_local_elements + num_ghost_elements));
  {
    t8_locidx_t itree, num_local_trees;
    t8_locidx_t current_index;
    t8_locidx_t ielement, num_elements_in_tree;
    t8_eclass_t tree_class;
    t8_eclass_scheme_c *eclass_scheme;
    const t8_element_t *element;

    /* Get the number of trees that have elements of this process. */
    num_local_trees = t8_forest_get_num_local_trees (forest);
    for (itree = 0, current_index = 0; itree < num_local_trees; ++itree) {
      /* This loop iterates through all local trees in the forest. */
      /* Each tree may have a different element class (quad/tri/hex/tet etc.) and therefore
       * also a different way to interpret its elements. In order to be able to handle elements
       * of a tree, we need to get its eclass_scheme, and in order to so we first get its eclass. */
      tree_class = t8_forest_get_tree_class (forest, itree);
      eclass_scheme = t8_forest_get_eclass_scheme (forest, tree_class);
      /* Get the number of elements of this tree. */
      num_elements_in_tree = t8_forest_get_tree_num_elements (forest, itree);
      for (ielement = 0; ielement < num_elements_in_tree; ++ielement, ++current_index) {
        element = t8_forest_get_element_in_tree (forest, itree, ielement);

        t8_forest_element_centroid(forest,itree,element,&element_data[current_index*3]);
      }
    }
  }
  return element_data;
}

void check_consistency(t8_forest_t forest, double *centroid)
{
  t8_locidx_t num_local_trees;
  t8_locidx_t num_local_elements;
  t8_locidx_t num_ghost_elements;
  t8_eclass_scheme_c *neigh_scheme;

  t8_element_t **neighbors;
  int num_neighbors;
  int *dual_faces;
  t8_locidx_t *neighs;
  double diff[6][3];
  double normal[6][3];
  double diff_norm,normal_norm,dot;
  double this_center[3];

  int ibad = 0;
  double eps = 1.e-12;

  t8_cmesh_t cmesh = t8_forest_get_cmesh(forest);
  if(t8_cmesh_trees_is_face_consistent (cmesh, cmesh->trees)){
    t8_productionf(" passed cmesh consistency test\n");
  }
  else{
    t8_productionf(" failed cmesh consistency test************************\n");
  }

  /* Get the number of local elements of forest. */
  num_local_elements = t8_forest_get_local_num_elements (forest);
  /* Get the number of ghost elements of forest. */
  num_ghost_elements = t8_forest_get_num_ghosts (forest);
  num_local_trees = t8_forest_get_num_local_trees (forest);
  t8_productionf(" num_local:%i\n",num_local_elements);

  for(t8_locidx_t itree = 0, elem_index = 0; itree < num_local_trees; ++itree) {

    t8_eclass_t tree_class = t8_forest_get_tree_class (forest, itree);
    t8_eclass_scheme_c *eclass_scheme = t8_forest_get_eclass_scheme (forest, tree_class);
    /* Get the number of elements of this tree. */
    int num_elements_in_tree = t8_forest_get_tree_num_elements (forest, itree);

    for (t8_locidx_t ielement = 0; ielement < num_elements_in_tree; ++ielement, ++elem_index) {
      const t8_element_t *element = t8_forest_get_element_in_tree (forest, itree, ielement);
      const int num_faces = eclass_scheme->t8_element_num_faces(element);
      const t8_element_shape_t element_shape = eclass_scheme->t8_element_shape(element);
      const int num_corners = t8_eclass_num_vertices[element_shape];

      this_center[0]  = centroid[elem_index*3+0];
      this_center[1]  = centroid[elem_index*3+1];
      this_center[2]  = centroid[elem_index*3+2];

      double center[3] = {0.0, 0.0, 0.0};
      for (int icorner = 0; icorner < num_corners; icorner++) {
        double coords[3];
        t8_forest_element_coordinate (forest, itree, element, icorner, coords);
      //t8_productionf("elem:%i icorner:%i coord:%f %f %f\n",elem_index,icorner,coords[0],coords[1],coords[2]);
        t8_vec_axpy (coords, center, 1);
      }
      t8_vec_ax (center, 1. / num_corners);
      t8_vec_axpy (this_center, center, -1);
      double norm = t8_vec_norm (center);
      if(norm > eps) {
        ibad++;
        t8_productionf(" bad centroid:%i  %f %f %f\n",elem_index,center[0],center[1],center[2]);
      }

      /* Set the faces */
      for (int iface = 0; iface < num_faces; iface++) {
        /* Compute the indices of the face neighbors */
        //t8_element_shape_t face_shape = eclass_scheme->t8_element_face_shape(element,iface);

        t8_forest_leaf_face_neighbors (forest, itree, element,
                                       &neighbors, iface,
                                       &dual_faces,
                                       &num_neighbors,
                                       &neighs, &neigh_scheme, 1);

        if(num_neighbors > 0) {
          t8_locidx_t neigh_index = neighs[0];

          double neig_center[3];
          neig_center[0] = centroid[neigh_index*3+0];
          neig_center[1] = centroid[neigh_index*3+1];
          neig_center[2] = centroid[neigh_index*3+2];

          t8_vec_diff(neig_center, this_center, diff[iface]);
          diff_norm = t8_vec_norm(diff[iface]);

          double fnorm[3];
          t8_forest_element_face_normal(forest,itree,element,iface,fnorm);
          dot = t8_vec_dot(diff[iface], fnorm)/(diff_norm);

          t8_productionf("elem:%i neigh:%i  "
                         "centroid_diff: % 8.6f % 8.6f % 8.6f  "
                         "normal: % 8.6f % 8.6f % 8.6f  "
                         "dot: % 8.6f\n",
                          elem_index,neigh_index,
                          diff[iface][0],diff[iface][1],diff[iface][2],
                          fnorm[0],fnorm[1],fnorm[2],
                          dot);

          neigh_scheme->t8_element_destroy (num_neighbors, neighbors);
          T8_FREE (neighs);
          T8_FREE (neighbors);
          T8_FREE (dual_faces);
        }
      }
    }
  }
  t8_productionf("# bad centroids: %i\n",ibad);
}

int
t8_step3_main (int argc, char **argv)
{
  int mpiret;
  sc_MPI_Comm comm;
  t8_cmesh_t cmesh;
  t8_forest_t forest;
  /* The prefix for our output files. */
  const char *prefix_uniform = "t8_step3_uniform_forest";

  /* The uniform refinement level of the forest. */
  const int level = 3;

  /* Initialize MPI. This has to happen before we initialize sc or t8code. */
  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);

  /* Print a message on the root process. */
  t8_global_productionf (" [step3] \n");
  t8_global_productionf (" [step3] Hello, this is the step3 example of t8code.\n");
  t8_global_productionf (" [step3] In this example we will refine and coarsen a forest.\n");
  t8_global_productionf (" [step3] \n");

  /* We will use MPI_COMM_WORLD as a communicator. */
  comm = sc_MPI_COMM_WORLD;

  /*
   * Setup.
   * Build cmesh and uniform forest.
   */
  const char *mshfile = (argc == 2) ? argv[1] : "box_pyr1_mesh";

  t8_productionf(" argc:%i %s\n",argc,mshfile);
  cmesh = t8_cmesh_from_msh_file (mshfile, 1, sc_MPI_COMM_WORLD, 3, 0, 0);
  if(t8_cmesh_trees_is_face_consistent (cmesh, cmesh->trees)){
    t8_productionf(" pass consistent test\n");
  } else{
    t8_productionf(" not pass consistent test************************\n");
  }

  t8_cmesh_t cmesh_partition;
  t8_cmesh_init(&cmesh_partition);
  t8_cmesh_set_partition_uniform(cmesh_partition,0,t8_scheme_new_default_cxx());
  t8_cmesh_set_derive(cmesh_partition,cmesh);
  t8_cmesh_commit(cmesh_partition,comm);
  forest = t8_forest_new_uniform (cmesh_partition, t8_scheme_new_default_cxx (), 0, 1, comm);

  /* Print information of the forest. */
  t8_step3_print_forest_information (forest);

  /* Write forest to vtu files. */
  t8_forest_write_vtk (forest, prefix_uniform);
  t8_global_productionf (" [step3] Wrote uniform forest to vtu files: %s*\n", prefix_uniform);

  /*
   *  Adapt the forest.
   */

  /* Adapt the forest. We can reuse the forest variable, since the new adapted
   * forest will take ownership of the old forest and destroy it.
   * Note that the adapted forest is a new forest, though. */
  forest = t8_step3_adapt_forest (forest);

  /*
   *  Partition the forest.
   */
  forest = t8_step3_partition_ghost (forest);

  /* Print information of the forest. */
  t8_step3_print_forest_information (forest);

  double *centroid;
  int nvar = 3;
  centroid = t8_step3_create_element_data (forest);
  t8_step3_exchange_ghost_data (forest, centroid, nvar);

  /*
   * perform face normal consistency check
   */
  check_consistency(forest, centroid);

  /*
   * clean-up
   */
  T8_FREE (centroid);

  /* Destroy the forest. */
  t8_forest_unref (&forest);
  t8_global_productionf (" [step3] Destroyed forest.\n");

  sc_finalize ();

  mpiret = sc_MPI_Finalize ();
  SC_CHECK_MPI (mpiret);

  return 0;
}

T8_EXTERN_C_END ();

[lukasdreyer]

Thank you for bringing this bug to our attention. I have created a branch with a possible fix: #863. Does this resolve your problem?

[lukasdreyer]

The bugfix #863 should solve the inconsistencies with t8_forest_element_face_normal, but it brought to light another problem that you might stumble upon: The node and face enumeration are different for the two pyramid types that appear, as you have already observed for other element types. Therefore, also the mapping from face to node differs for the element different pyramid types.
As we discussed yesterday, this is something that we want to standardize for t8code elements. If you have any questions regarding the current state of vertex and face numbering for different element types, also feel free to discuss in mattermost.

[ackirby88]

Thanks, Lukas, for the quick fix! It appears it is working now and we are doing more tests on complex problems to be sure. For the face to node orderings, we perform a negative volume check to classify the types and do an internal reordering.