diff --git a/src/mesh/Connectivity2D.hpp b/src/mesh/Connectivity2D.hpp
index b64d77a5ad25540c975c1584117a154a948152bb..e4d8c9c0c8bcf5c12dce93e497e4158d07be480b 100644
--- a/src/mesh/Connectivity2D.hpp
+++ b/src/mesh/Connectivity2D.hpp
@@ -16,10 +16,10 @@
 
 class Connectivity2D
 {
-public:
+ public:
   static constexpr size_t dimension = 2;
 
-private:
+ private:
   std::vector<RefFaceList> m_ref_face_list;
   std::vector<RefNodeList> m_ref_node_list;
 
@@ -317,9 +317,9 @@ private:
 
   Connectivity2D(const Kokkos::View<const unsigned short*> cell_nb_nodes,
                  const Kokkos::View<const unsigned int**> cell_nodes)
-    : m_number_of_cells (cell_nodes.extent(0)),
-      m_cell_nb_nodes(cell_nb_nodes),
-      m_cell_nodes (cell_nodes)
+      : m_number_of_cells (cell_nodes.extent(0)),
+        m_cell_nb_nodes(cell_nb_nodes),
+        m_cell_nodes (cell_nodes)
   {
     {
       Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", m_number_of_cells);
diff --git a/src/mesh/Connectivity3D.hpp b/src/mesh/Connectivity3D.hpp
index c677fdefad8043740123a23b13bdafd630bdea7a..5b65c12e81d52f3a8a573eb71997523648ff3dc8 100644
--- a/src/mesh/Connectivity3D.hpp
+++ b/src/mesh/Connectivity3D.hpp
@@ -14,6 +14,7 @@
 #include <RefNodeList.hpp>
 #include <RefFaceList.hpp>
 
+
 class Connectivity3D
 {
 public:
@@ -48,137 +49,206 @@ private:
 
   size_t  m_max_nb_node_per_cell;
 
-  // struct Face
-  // {
-  //   const unsigned int m_node0_id;
-  //   const unsigned int m_node1_id;
-
-  //   KOKKOS_INLINE_FUNCTION
-  //   bool operator<(const Face& f) const
-  //   {
-  //     return ((m_node0_id<f.m_node0_id) or
-  //             ((m_node0_id == f.m_node0_id) and
-  //              (m_node1_id<f.m_node1_id)));
-  //   }
-
-  //   KOKKOS_INLINE_FUNCTION
-  //   Face& operator=(const Face&) = default;
-
-  //   KOKKOS_INLINE_FUNCTION
-  //   Face& operator=(Face&&) = default;
-
-  //   KOKKOS_INLINE_FUNCTION
-  //   Face(const Face&) = default;
-
-  //   KOKKOS_INLINE_FUNCTION
-  //   Face(Face&&) = default;
-
-  //   KOKKOS_INLINE_FUNCTION
-  //   Face(unsigned int node0_id,
-  //        unsigned int node1_id)
-  //       : m_node0_id(node0_id),
-  //         m_node1_id(node1_id)
-  //   {
-  //     ;
-  //   }
-
-  //   KOKKOS_INLINE_FUNCTION
-  //   ~Face() = default;
-  // };
-
-  // void _computeFaceCellConnectivities()
-  // {
-  //   // In 2D faces are simply define
-  //   typedef std::pair<unsigned int, unsigned short> CellFaceId;
-  //   std::map<Face, std::vector<CellFaceId>> face_cells_map;
-  //   for (unsigned int j=0; j<m_number_of_cells; ++j) {
-  //     const unsigned short cell_nb_nodes = m_cell_nb_nodes[j];
-  //     for (unsigned short r=0; r<cell_nb_nodes; ++r) {
-  //       unsigned int node0_id = m_cell_nodes(j,r);
-  //       unsigned int node1_id = m_cell_nodes(j,(r+1)%cell_nb_nodes);
-  //       if (node1_id<node0_id) {
-  //         std::swap(node0_id, node1_id);
-  //       }
-  //       face_cells_map[Face(node0_id, node1_id)].push_back(std::make_pair(j, r));
-  //     }
-  //   }
-
-  //   m_number_of_faces = face_cells_map.size();
-  //   Kokkos::View<unsigned short*> face_nb_nodes("face_nb_nodes", m_number_of_faces);
-  //   Kokkos::parallel_for(m_number_of_faces, KOKKOS_LAMBDA(const unsigned int& l) {
-  //       face_nb_nodes[l] = 2;
-  //     });
-  //   m_face_nb_nodes = face_nb_nodes;
-
-  //   Kokkos::View<unsigned int*[2]> face_nodes("face_nodes", m_number_of_faces, 2);
-  //   {
-  //     int l=0;
-  //     for (const auto& face_cells_vector : face_cells_map) {
-  //       const Face& face = face_cells_vector.first;
-  //       face_nodes(l,0) = face.m_node0_id;
-  //       face_nodes(l,1) = face.m_node1_id;
-  //       ++l;
-  //     }
-  //   }
-  //   m_face_nodes = face_nodes;
-
-  //   Kokkos::View<unsigned short*> face_nb_cells("face_nb_cells", m_number_of_faces);
-  //   {
-  //     int l=0;
-  //     for (const auto& face_cells_vector : face_cells_map) {
-  //       const auto& cells_vector = face_cells_vector.second;
-  //       face_nb_cells[l] = cells_vector.size();
-  //       ++l;
-  //     }
-  //   }
-  //   m_face_nb_cells = face_nb_cells;
-
-  //   Kokkos::View<unsigned int**> face_cells("face_cells", face_cells_map.size(), 2);
-  //   {
-  //     int l=0;
-  //     for (const auto& face_cells_vector : face_cells_map) {
-  //       const auto& cells_vector = face_cells_vector.second;
-  //       for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
-  //         unsigned int cell_number = cells_vector[lj].first;
-  //         face_cells(l,lj) = cell_number;
-  //       }
-  //       ++l;
-  //     }
-  //   }
-  //   m_face_cells = face_cells;
-
-  //   // In 2d cell_nb_faces = cell_nb_node
-  //   m_cell_nb_faces = m_cell_nb_nodes;
-  //   Kokkos::View<unsigned int**> cell_faces("cell_faces", m_number_of_faces, m_max_nb_node_per_cell);
-  //   {
-  //     int l=0;
-  //     for (const auto& face_cells_vector : face_cells_map) {
-  //       const auto& cells_vector = face_cells_vector.second;
-  //       for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
-  //         unsigned int cell_number = cells_vector[lj].first;
-  //         unsigned short cell_local_face = cells_vector[lj].second;
-  //         cell_faces(cell_number,cell_local_face) = l;
-  //       }
-  //       ++l;
-  //     }
-  //   }
-  //   m_cell_faces = cell_faces;
-
-  //   Kokkos::View<unsigned short**> face_cell_local_face("face_cell_local_face",
-  //                                                       m_number_of_faces, m_max_nb_node_per_cell);
-  //   {
-  //     int l=0;
-  //     for (const auto& face_cells_vector : face_cells_map) {
-  //       const auto& cells_vector = face_cells_vector.second;
-  //       for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
-  //         unsigned short cell_local_face = cells_vector[lj].second;
-  //         face_cell_local_face(l,lj) = cell_local_face;
-  //       }
-  //       ++l;
-  //     }
-  //   }
-  //   m_face_cell_local_face = face_cell_local_face;
-  // }
+  struct Face
+  {
+    KOKKOS_INLINE_FUNCTION
+    std::vector<unsigned int> _sort(const std::vector<unsigned int>& node_list) const
+    {
+      auto min_id = std::min_element(node_list.begin(), node_list.end());
+      const size_t shift = std::distance(min_id, node_list.begin());
+
+      std::vector<unsigned int> rotated_node_list(node_list.size());
+      if (node_list[(shift+1)%node_list.size()] < node_list[(shift+node_list.size()-1)%node_list.size()]) {
+        for (size_t i=0; i<node_list.size(); ++i) {
+          const size_t reverse_shift = node_list.size()-shift;
+          rotated_node_list[(reverse_shift+node_list.size()-i)%node_list.size()] = node_list[i];
+        }
+      } else {
+        for (size_t i=0; i<node_list.size(); ++i) {
+          rotated_node_list[(shift+i)%node_list.size()] = node_list[i];
+        }
+      }
+
+      return rotated_node_list;
+    }
+
+    std::vector<unsigned int> m_node_id_list;
+
+    KOKKOS_INLINE_FUNCTION
+    bool operator<(const Face& f) const
+    {
+      const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
+      for (size_t i=0; i<min_nb_nodes; ++i) {
+        if (m_node_id_list[i] <  f.m_node_id_list[i]) return true;
+        if (m_node_id_list[i] != f.m_node_id_list[i]) return false;
+      }
+      return m_node_id_list.size() < f.m_node_id_list.size();
+    }
+
+    KOKKOS_INLINE_FUNCTION
+    Face& operator=(const Face&) = default;
+
+    KOKKOS_INLINE_FUNCTION
+    Face& operator=(Face&&) = default;
+
+    KOKKOS_INLINE_FUNCTION
+    Face(const Face&) = default;
+
+    KOKKOS_INLINE_FUNCTION
+    Face(Face&&) = default;
+
+    KOKKOS_INLINE_FUNCTION
+    Face(const std::vector<unsigned int>& given_node_id_list)
+        : m_node_id_list(_sort(given_node_id_list))
+    {
+      ;
+    }
+
+    KOKKOS_INLINE_FUNCTION
+    ~Face() = default;
+  };
+
+  void _computeFaceCellConnectivities()
+  {
+    // In 2D faces are simply define
+    typedef std::pair<unsigned int, unsigned short> CellFaceId;
+    std::map<Face, std::vector<CellFaceId>> face_cells_map;
+    for (unsigned int j=0; j<m_number_of_cells; ++j) {
+      const unsigned short cell_nb_nodes = m_cell_nb_nodes[j];
+      switch (cell_nb_nodes) {
+        case 4: { // tetrahedron
+          std::cerr << "tetrahedron cell type NIY!\n";
+          std::exit(0);
+          break;
+        }
+        case 8: { // hexahedron
+          // face 0
+          face_cells_map[Face({m_cell_nodes(j,0),
+                               m_cell_nodes(j,1),
+                               m_cell_nodes(j,2),
+                               m_cell_nodes(j,3)})].push_back(std::make_pair(j,0));
+          // face 1
+          face_cells_map[Face({m_cell_nodes(j,4),
+                               m_cell_nodes(j,5),
+                               m_cell_nodes(j,6),
+                               m_cell_nodes(j,7)})].push_back(std::make_pair(j,1));
+          // face 2
+          face_cells_map[Face({m_cell_nodes(j,0),
+                               m_cell_nodes(j,3),
+                               m_cell_nodes(j,7),
+                               m_cell_nodes(j,4)})].push_back(std::make_pair(j,2));
+          // face 3
+          face_cells_map[Face({m_cell_nodes(j,1),
+                               m_cell_nodes(j,2),
+                               m_cell_nodes(j,6),
+                               m_cell_nodes(j,5)})].push_back(std::make_pair(j,3));
+          // face 4
+          face_cells_map[Face({m_cell_nodes(j,0),
+                               m_cell_nodes(j,1),
+                               m_cell_nodes(j,5),
+                               m_cell_nodes(j,4)})].push_back(std::make_pair(j,4));
+          // face 1
+          face_cells_map[Face({m_cell_nodes(j,3),
+                               m_cell_nodes(j,2),
+                               m_cell_nodes(j,6),
+                               m_cell_nodes(j,7)})].push_back(std::make_pair(j,5));
+          break;
+        }
+        default: {
+          std::cerr << "unexpected cell type!\n";
+          std::exit(0);
+        }
+      }
+      for (unsigned short r=0; r<cell_nb_nodes; ++r) {
+        unsigned int node0_id = m_cell_nodes(j,r);
+        unsigned int node1_id = m_cell_nodes(j,(r+1)%cell_nb_nodes);
+        if (node1_id<node0_id) {
+          std::swap(node0_id, node1_id);
+        }
+        // face_cells_map[Face(node0_id, node1_id)].push_back(std::make_pair(j, r));
+      }
+    }
+
+    m_number_of_faces = face_cells_map.size();
+    Kokkos::View<unsigned short*> face_nb_nodes("face_nb_nodes", m_number_of_faces);
+    Kokkos::parallel_for(m_number_of_faces, KOKKOS_LAMBDA(const unsigned int& l) {
+        face_nb_nodes[l] = 4;
+      });
+    m_face_nb_nodes = face_nb_nodes;
+
+    std::cerr << __FILE__ << ':' << __LINE__ << ':'
+              << rang::fg::red << "m_max_nb_node_per_face forced to 8" << rang::fg::reset;
+    unsigned int m_max_nb_node_per_face = 8;
+    Kokkos::View<unsigned int**> face_nodes("face_nodes", m_number_of_faces, m_max_nb_node_per_face);
+    {
+      int l=0;
+      for (const auto& face_cells_vector : face_cells_map) {
+        const Face& face = face_cells_vector.first;
+        for(size_t r=0; r<face.m_node_id_list.size(); ++r) {
+          face_nodes(l,r) = face.m_node_id_list[r];
+        }
+        ++l;
+      }
+    }
+    m_face_nodes = face_nodes;
+
+    Kokkos::View<unsigned short*> face_nb_cells("face_nb_cells", m_number_of_faces);
+    {
+      int l=0;
+      for (const auto& face_cells_vector : face_cells_map) {
+        const auto& cells_vector = face_cells_vector.second;
+        face_nb_cells[l] = cells_vector.size();
+        ++l;
+      }
+    }
+    m_face_nb_cells = face_nb_cells;
+
+    Kokkos::View<unsigned int**> face_cells("face_cells", face_cells_map.size(), 2);
+    {
+      int l=0;
+      for (const auto& face_cells_vector : face_cells_map) {
+        const auto& cells_vector = face_cells_vector.second;
+        for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
+          unsigned int cell_number = cells_vector[lj].first;
+          face_cells(l,lj) = cell_number;
+        }
+        ++l;
+      }
+    }
+    m_face_cells = face_cells;
+
+    // In 2d cell_nb_faces = cell_nb_node
+    m_cell_nb_faces = m_cell_nb_nodes;
+    Kokkos::View<unsigned int**> cell_faces("cell_faces", m_number_of_faces, m_max_nb_node_per_cell);
+    {
+      int l=0;
+      for (const auto& face_cells_vector : face_cells_map) {
+        const auto& cells_vector = face_cells_vector.second;
+        for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
+          unsigned int cell_number = cells_vector[lj].first;
+          unsigned short cell_local_face = cells_vector[lj].second;
+          cell_faces(cell_number,cell_local_face) = l;
+        }
+        ++l;
+      }
+    }
+    m_cell_faces = cell_faces;
+
+    Kokkos::View<unsigned short**> face_cell_local_face("face_cell_local_face",
+                                                        m_number_of_faces, m_max_nb_node_per_cell);
+    {
+      int l=0;
+      for (const auto& face_cells_vector : face_cells_map) {
+        const auto& cells_vector = face_cells_vector.second;
+        for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
+          unsigned short cell_local_face = cells_vector[lj].second;
+          face_cell_local_face(l,lj) = cell_local_face;
+        }
+        ++l;
+      }
+    }
+    m_face_cell_local_face = face_cell_local_face;
+  }
 
  public:
   void addRefFaceList(const RefFaceList& ref_face_list)
@@ -313,7 +383,7 @@ private:
     return 0;
   }
 
-  Connectivity3D(const Connectivity2D&) = delete;
+  Connectivity3D(const Connectivity3D&) = delete;
 
   Connectivity3D(const Kokkos::View<const unsigned short*> cell_nb_nodes,
                  const Kokkos::View<const unsigned int**> cell_nodes)
@@ -340,7 +410,7 @@ private:
                                       m_node_cells,
                                       m_node_cell_local_node);
 
-    // this->_computeFaceCellConnectivities();
+    this->_computeFaceCellConnectivities();
   }
 
   ~Connectivity3D()
diff --git a/src/mesh/ConnectivityUtils.hpp b/src/mesh/ConnectivityUtils.hpp
index 8b74b23a4f036777dd972c5b4a4ccc1d45003506..99d78ddea7a61964781b0ed1230ee44b1f92f66d 100644
--- a/src/mesh/ConnectivityUtils.hpp
+++ b/src/mesh/ConnectivityUtils.hpp
@@ -2,30 +2,31 @@
 #define CONNECTIVITY_UTILS_HPP
 
 #include <map>
+#include <Kokkos_Core.hpp>
 
 class ConnectivityUtils
 {
-public:
+ public:
   void computeNodeCellConnectivity(const Kokkos::View<const unsigned int**> cell_nodes,
-				   const Kokkos::View<const unsigned short*> cell_nb_nodes,
-				   const size_t& number_of_cells,
-				   size_t& max_nb_node_per_cell,
-				   size_t& number_of_nodes,
-				   Kokkos::View<unsigned short*>& node_nb_cells,
-				   Kokkos::View<unsigned int**>& node_cells,
-				   Kokkos::View<unsigned short**>& node_cell_local_node)
+                                   const Kokkos::View<const unsigned short*> cell_nb_nodes,
+                                   const size_t& number_of_cells,
+                                   size_t& max_nb_node_per_cell,
+                                   size_t& number_of_nodes,
+                                   Kokkos::View<unsigned short*>& node_nb_cells,
+                                   Kokkos::View<unsigned int**>& node_cells,
+                                   Kokkos::View<unsigned short**>& node_cell_local_node)
   {
     // Computes inefficiently node->cells connectivity [Version 0]
     std::multimap<unsigned int, unsigned int> node_cells_map;
     using namespace Kokkos::Experimental;
     Kokkos::parallel_reduce(number_of_cells, KOKKOS_LAMBDA(const int& j, size_t& nb_max) {
-	const size_t n = cell_nb_nodes[j];
-	if (n > nb_max) nb_max = n;
+        const size_t n = cell_nb_nodes[j];
+        if (n > nb_max) nb_max = n;
       }, Kokkos::Max<size_t>(max_nb_node_per_cell));
 
     for (unsigned int j=0; j<number_of_cells; ++j) {
       for (unsigned int r=0; r<cell_nb_nodes[j]; ++r) {
-	node_cells_map.insert(std::make_pair(cell_nodes(j,r),j));
+        node_cells_map.insert(std::make_pair(cell_nodes(j,r),j));
       }
     }
 
@@ -42,7 +43,7 @@ public:
     if ((node_ids[0] != 0) or (node_ids[node_ids.size()-1] != node_ids.size()-1)) {
       std::cerr << "sparse node numerotation NIY\n";
       for (size_t i=0; i<node_ids.size(); ++i) {
-	std::cout << "node_ids[" << i << "] = " << node_ids[i] << '\n';
+        std::cout << "node_ids[" << i << "] = " << node_ids[i] << '\n';
       }
       std::exit(0);
     }
@@ -59,7 +60,7 @@ public:
       const size_t nb_cells = cells_vector.size();
       node_nb_cells[i] = nb_cells;
       if (nb_cells > max_node_cells) {
-	max_node_cells = nb_cells;
+        max_node_cells = nb_cells;
       }
     }
 
@@ -67,22 +68,22 @@ public:
     for (size_t i=0; i<node_cells_vector.size(); ++i) {
       const auto& cells_vector = node_cells_vector[i];
       for (size_t j=0; j<cells_vector.size(); ++j) {
-	node_cells(i,j) = cells_vector[j];
+        node_cells(i,j) = cells_vector[j];
       }
     }
 
     node_cell_local_node = Kokkos::View<unsigned short**>("node_cell_local_node",
-							  node_ids.size(), max_node_cells);
+                                                          node_ids.size(), max_node_cells);
     Kokkos::parallel_for(number_of_nodes, KOKKOS_LAMBDA(const unsigned int& r){
-	for (unsigned short J=0; J<node_nb_cells[r]; ++J) {
-	  const unsigned int j = node_cells(r,J);
-	  for (unsigned int R=0; R<cell_nb_nodes[j]; ++R) {
-	    if (cell_nodes(j,R) == r) {
-	      node_cell_local_node(r,J)=R;
-	      break;
-	    }
-	  }
-	}
+        for (unsigned short J=0; J<node_nb_cells[r]; ++J) {
+          const unsigned int j = node_cells(r,J);
+          for (unsigned int R=0; R<cell_nb_nodes[j]; ++R) {
+            if (cell_nodes(j,R) == r) {
+              node_cell_local_node(r,J)=R;
+              break;
+            }
+          }
+        }
       });
   }
 };