Fix pyramid face corner

src/t8_schemes/t8_default/t8_default_pyramid/t8_dpyramid_bits.c

@@ -1401,7 +1401,7 @@ t8_dpyramid_get_face_corner (const t8_dpyramid_t *pyra, int face, int corner)
     return t8_dtet_face_corner[face][corner];
   }
   else {
-    int corner_number = t8_dpyramid_face_corner[face][corner];
+    const int corner_number = t8_dpyramid_face_corner[pyra->pyramid.type-T8_DPYRAMID_FIRST_TYPE][face][corner];
     T8_ASSERT (0 <= corner_number && corner_number < T8_DPYRAMID_FACES);
     return corner_number;
   }

src/t8_schemes/t8_default/t8_default_pyramid/t8_dpyramid_connectivity.c

@@ -128,10 +128,19 @@ const int t8_dpyramid_tritype_rootface_to_face[2][4] = {
   { 2, 0, 2, 0 }, 
   { 1, 0, 1, 0 } };
 
-const int t8_dpyramid_face_corner[5][4] = { 
+const int t8_dpyramid_face_corner[2][5][4] = { 
+{
   { 0, 2, 4, -1 }, 
   { 1, 3, 4, -1 }, 
   { 0, 1, 4, -1 }, 
   { 2, 3, 4, -1 }, 
-  { 0, 1, 2, 3 } };
+  { 0, 1, 2, 3 }
+},{
+  { 2, 3, 4, -1 }, 
+  { 0, 1, 4, -1 }, 
+  { 1, 3, 4, -1 }, 
+  { 0, 2, 4, -1 }, 
+  { 0, 1, 2, 3 }
+}
+};
 /* clang-format on */

src/t8_schemes/t8_default/t8_default_pyramid/t8_dpyramid_connectivity.h

@@ -111,8 +111,8 @@ extern const int t8_dpyramid_tritype_rootface_to_face[2][4];
 extern const int t8_dtet_type_cid_to_pyramid_parenttype[6][8];
 
 /** The corner numbers of the face of a pyramid
- * corner_number = A(face, corner_number_at_face)
+ * corner_number = A(pyramid_type, face, corner_number_at_face)
 */
-extern const int t8_dpyramid_face_corner[5][4];
+extern const int t8_dpyramid_face_corner[2][5][4];
 
 #endif  // T8_DPYRAMID_CONNECTIVITY_H

test/Makefile.am

@@ -40,6 +40,7 @@ t8code_googletest_programs = \
   test/t8_data/t8_gtest_shmem \
   test/t8_forest/t8_gtest_half_neighbors \
   test/t8_forest/t8_gtest_find_owner \
+  test/t8_forest/t8_gtest_forest_face_normal \
   test/t8_schemes/t8_gtest_face_descendant \
   test/t8_geometry/t8_gtest_point_inside \
   test/t8_forest/t8_gtest_user_data \
@@ -198,6 +199,10 @@ test_t8_forest_t8_gtest_find_owner_SOURCES = \
   test/t8_gtest_main.cxx \
   test/t8_forest/t8_gtest_find_owner.cxx
 
+test_t8_forest_t8_gtest_forest_face_normal_SOURCES = \
+  test/t8_gtest_main.cxx \
+  test/t8_forest/t8_gtest_forest_face_normal.cxx
+
 test_t8_schemes_t8_gtest_face_descendant_SOURCES = \
   test/t8_gtest_main.cxx \
   test/t8_schemes/t8_gtest_face_descendant.cxx
@@ -395,6 +400,10 @@ test_t8_forest_t8_gtest_find_owner_LDADD = $(t8_gtest_target_ld_add)
 test_t8_forest_t8_gtest_find_owner_LDFLAGS = $(t8_gtest_target_ld_flags)
 test_t8_forest_t8_gtest_find_owner_CPPFLAGS = $(t8_gtest_target_cpp_flags)
 
+test_t8_forest_t8_gtest_forest_face_normal_LDADD = $(t8_gtest_target_ld_add)
+test_t8_forest_t8_gtest_forest_face_normal_LDFLAGS = $(t8_gtest_target_ld_flags)
+test_t8_forest_t8_gtest_forest_face_normal_CPPFLAGS = $(t8_gtest_target_cpp_flags)
+
 test_t8_schemes_t8_gtest_face_descendant_LDADD = $(t8_gtest_target_ld_add)
 test_t8_schemes_t8_gtest_face_descendant_LDFLAGS = $(t8_gtest_target_ld_flags)
 test_t8_schemes_t8_gtest_face_descendant_CPPFLAGS = $(t8_gtest_target_cpp_flags)
@@ -486,6 +495,7 @@ test_t8_gtest_vtk_linkage_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_data_t8_gtest_shmem_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_forest_t8_gtest_half_neighbors_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_forest_t8_gtest_find_owner_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
+test_t8_forest_t8_gtest_forest_face_normal_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_schemes_t8_gtest_face_descendant_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_geometry_t8_gtest_point_inside_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)
 test_t8_forest_t8_gtest_user_data_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags)

test/t8_forest/t8_gtest_forest_face_normal.cxx

@@ -0,0 +1,120 @@
+/*
+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) 2023 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/src/sc_functions.h>
+#include <gtest/gtest.h>
+#include <t8_eclass.h>
+#include <t8_schemes/t8_default/t8_default_cxx.hxx>
+#include <t8_schemes/t8_default/t8_default_pyramid/t8_dpyramid_bits.h>
+#include <t8_cmesh/t8_cmesh_examples.h>
+#include <t8_forest/t8_forest_general.h>
+#include <t8_forest/t8_forest_geometrical.h>
+#include <test/t8_gtest_macros.hxx>
+
+/**
+ * This file tests the face normal computation of elements.
+ */
+
+class class_forest_face_normal: public testing::TestWithParam<std::tuple<t8_eclass_t, int>> {
+ protected:
+  void
+  SetUp () override
+  {
+    eclass = std::get<0> (GetParam ());
+    level = std::get<1> (GetParam ());
+    scheme = t8_scheme_new_default_cxx ();
+    t8_cmesh_t cmesh = t8_cmesh_new_hypercube (eclass, sc_MPI_COMM_WORLD, 0, 0, 0);
+    const int do_face_ghost = 1;
+    forest = t8_forest_new_uniform (cmesh, scheme, level, do_face_ghost, sc_MPI_COMM_WORLD);
+  }
+  void
+  TearDown () override
+  {
+    t8_forest_unref (&forest);
+  }
+  t8_forest_t forest;
+  t8_scheme_cxx *scheme;
+  t8_eclass_t eclass;
+  int level;
+};
+
+TEST_P (class_forest_face_normal, back_and_forth)
+{
+  /** Iterate over all elements of a uniformly refined forest. For all faceneighbors elements, if they are local,
+    * check if their facenormal is the negative of the corresponding facenormal of the neighbor elements.
+    */
+  const t8_locidx_t local_num_trees = t8_forest_get_num_local_trees (forest);
+  /* Iterate over all elements. */
+  for (t8_locidx_t itree = 0; itree < local_num_trees; itree++) {
+    const t8_locidx_t tree_elements = t8_forest_get_tree_num_elements (forest, itree);
+    const t8_eclass_t tree_eclass = t8_forest_get_tree_class (forest, itree);
+    ASSERT_EQ(eclass, tree_eclass);
+    const t8_eclass_scheme_c* escheme = t8_forest_get_eclass_scheme(forest, tree_eclass);
+    for (t8_locidx_t ielement = 0; ielement < tree_elements; ielement++) {
+      const t8_element_t *element = t8_forest_get_element_in_tree (forest, itree, ielement);
+      const int num_faces = escheme->t8_element_num_faces(element);
+      for (int iface = 0; iface < num_faces; iface++){
+        /* Compute facenormal */
+        double face_normal[3] = {0,0,0};
+        t8_forest_element_face_normal(forest,itree,element,iface,face_normal);
+
+        /* Get all face neighbors */
+
+        t8_element_t **neighbors;
+        int num_neighbors;
+        const int forest_is_balanced = 1;
+        t8_eclass_scheme_c *neigh_scheme;
+        int *dual_faces;
+        t8_locidx_t *neigh_ids;
+
+
+        t8_forest_leaf_face_neighbors (forest, itree, element, &neighbors, iface, &dual_faces, &num_neighbors,
+                                       &neigh_ids, &neigh_scheme, forest_is_balanced);
+
+        /* Iterate and compute their facenormal */
+
+        for(int ineigh = 0; ineigh < num_neighbors; ineigh++){
+          double neigh_face_normal[3];
+          t8_locidx_t ineightree;
+          /* Skip ghost elements, because those can't be passed to t8_forest_element_face_normal */
+          if (neigh_ids[ineigh] >= t8_forest_get_local_num_elements(forest))continue;
+
+          t8_element_t *neigh_elem = t8_forest_get_element(forest, neigh_ids[ineigh], &ineightree);
+          t8_forest_element_face_normal(forest, ineightree, neigh_elem, dual_faces[ineigh], neigh_face_normal);
+          EXPECT_NEAR(face_normal[0], -neigh_face_normal[0], 10*T8_PRECISION_EPS);
+          EXPECT_NEAR(face_normal[1], -neigh_face_normal[1], 10*T8_PRECISION_EPS);
+          EXPECT_NEAR(face_normal[2], -neigh_face_normal[2], 10*T8_PRECISION_EPS);
+        }
+
+        if(num_neighbors > 0) {
+          neigh_scheme->t8_element_destroy (num_neighbors, neighbors);
+          T8_FREE (neigh_ids);
+          T8_FREE (neighbors);
+          T8_FREE (dual_faces);
+        }
+      }
+    }
+  }
+}
+
+INSTANTIATE_TEST_SUITE_P (t8_gtest_forest_face_normal, class_forest_face_normal,
+                          testing::Combine (AllEclasses, testing::Range (0, 2)));