diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 0dd42769dc5fcd4b4864ee552f2c005716e5bc22..2093edd338aff76d9e55893d133964591b619a70 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -55,6 +55,15 @@ public:
     return std::move(tA);
   }
 
+  KOKKOS_INLINE_FUNCTION
+  TinyMatrix& operator*=(const T& t)
+  {
+    for (size_t i=0; i<N*N; ++i) {
+      m_values[i] *= t;
+    }
+    return *this;
+  }
+
   KOKKOS_INLINE_FUNCTION
   TinyMatrix operator*(const TinyMatrix& B) const
   {
@@ -358,4 +367,47 @@ T det(const TinyMatrix<3,T>& A)
     +A(2,0)*(A(0,1)*A(1,2)-A(1,1)*A(0,2));
 }
 
+template <size_t N, typename T>
+KOKKOS_INLINE_FUNCTION
+TinyMatrix<N,T> inverse(const TinyMatrix<N,T>& A);
+
+template <typename T>
+KOKKOS_INLINE_FUNCTION
+TinyMatrix<1,T> inverse(const TinyMatrix<1,T>& A)
+{
+  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non arithmetic types");
+  static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
+
+  TinyMatrix<1,T> A_1(1./A(0,0));
+  return std::move(A_1);
+}
+
+template <typename T>
+KOKKOS_INLINE_FUNCTION
+TinyMatrix<2,T> inverse(const TinyMatrix<2,T>& A)
+{
+  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non arithmetic types");
+  static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
+
+  const T determinent = det(A);
+  const T inv_determinent = 1./determinent;
+
+  TinyMatrix<2,T> cofactor_T(A(1,1), -A(1,0),
+                             -A(0,1), A(0,0));
+  return std::move(cofactor_T *= inv_determinent);
+}
+
+template <typename T>
+KOKKOS_INLINE_FUNCTION
+TinyMatrix<3,T> inverse(const TinyMatrix<3,T>& A)
+{
+  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non arithmetic types");
+  static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
+
+  std::cerr << "inverse not implemented in 3D\n";
+  std::exit(1);
+
+  return std::move(TinyMatrix<3,T>(identity));
+}
+
 #endif // TINYMATRIX_HPP
diff --git a/src/main.cpp b/src/main.cpp
index 6d1016fb58ac752c94cff745958c2ae8b0108c0a..4d4e70be7564dd4624562bf46eb9e3052d33dedd 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -12,6 +12,9 @@
 // #include <AcousticSolverTest.hpp>
 
 #include <Connectivity1D.hpp>
+#include <Connectivity2D.hpp>
+#include <Connectivity3D.hpp>
+
 #include <Mesh.hpp>
 #include <BoundaryCondition.hpp>
 #include <AcousticSolver.hpp>
@@ -340,6 +343,65 @@ int main(int argc, char *argv[])
           break;
         }
         case 3: {
+          typedef Connectivity3D ConnectivityType;
+          typedef Mesh<ConnectivityType> MeshType;
+          typedef MeshData<MeshType> MeshDataType;
+          typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType;
+
+          const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+
+          Kokkos::Timer timer;
+          timer.reset();
+          MeshDataType mesh_data(mesh);
+
+          std::vector<BoundaryConditionHandler> bc_list;
+
+          UnknownsType unknowns(mesh_data);
+
+          unknowns.initializeSod();
+
+          AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
+
+          typedef TinyVector<MeshType::dimension> Rd;
+
+          const Kokkos::View<const double*> Vj = mesh_data.Vj();
+          const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
+
+          const double tmax=0.2;
+          double t=0;
+
+          int itermax=std::numeric_limits<int>::max();
+          int iteration=0;
+
+          Kokkos::View<double*> rhoj = unknowns.rhoj();
+          Kokkos::View<double*> ej = unknowns.ej();
+          Kokkos::View<double*> pj = unknowns.pj();
+          Kokkos::View<double*> gammaj = unknowns.gammaj();
+          Kokkos::View<double*> cj = unknowns.cj();
+
+          BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+
+          VTKWriter vtk_writer("mesh", 0.01);
+
+          while((t<tmax) and (iteration<itermax)) {
+            vtk_writer.write(mesh, t);
+            double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
+            if (t+dt>tmax) {
+              dt=tmax-t;
+            }
+            acoustic_solver.computeNextStep(t,dt, unknowns);
+
+            block_eos.updatePandCFromRhoE();
+
+            t += dt;
+            ++iteration;
+          }
+          vtk_writer.write(mesh, t, true); // forces last output
+
+          std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
+                    << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+
+          method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
           break;
         }
       }
diff --git a/src/mesh/Connectivity3D.hpp b/src/mesh/Connectivity3D.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c677fdefad8043740123a23b13bdafd630bdea7a
--- /dev/null
+++ b/src/mesh/Connectivity3D.hpp
@@ -0,0 +1,352 @@
+#ifndef CONNECTIVITY_3D_HPP
+#define CONNECTIVITY_3D_HPP
+
+#include <Kokkos_Core.hpp>
+#include <PastisAssert.hpp>
+#include <TinyVector.hpp>
+
+#include <ConnectivityUtils.hpp>
+#include <vector>
+#include <map>
+#include <algorithm>
+
+#include <RefId.hpp>
+#include <RefNodeList.hpp>
+#include <RefFaceList.hpp>
+
+class Connectivity3D
+{
+public:
+  static constexpr size_t dimension = 3;
+
+private:
+  std::vector<RefFaceList> m_ref_face_list;
+  std::vector<RefNodeList> m_ref_node_list;
+
+  const size_t  m_number_of_cells;
+  size_t m_number_of_faces;
+  size_t m_number_of_nodes;
+
+  const Kokkos::View<const unsigned short*> m_cell_nb_nodes;
+  const Kokkos::View<const unsigned int**> m_cell_nodes;
+  Kokkos::View<double*> m_inv_cell_nb_nodes;
+
+  Kokkos::View<const unsigned short*> m_cell_nb_faces;
+  Kokkos::View<unsigned int**> m_cell_faces;
+
+  Kokkos::View<unsigned short*> m_node_nb_cells;
+  Kokkos::View<unsigned int**> m_node_cells;
+  Kokkos::View<unsigned short**> m_node_cell_local_node;
+
+  Kokkos::View<unsigned short*> m_face_nb_cells;
+  Kokkos::View<unsigned int**> m_face_cells;
+  Kokkos::View<unsigned short**> m_face_cell_local_face;
+
+  Kokkos::View<unsigned short*> m_face_nb_nodes;
+  Kokkos::View<unsigned int**> m_face_nodes;
+  Kokkos::View<unsigned short**> m_face_node_local_face;
+
+  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;
+  // }
+
+ public:
+  void addRefFaceList(const RefFaceList& ref_face_list)
+  {
+    m_ref_face_list.push_back(ref_face_list);
+  }
+
+  size_t numberOfRefFaceList() const
+  {
+    return m_ref_face_list.size();
+  }
+
+  const RefFaceList& refFaceList(const size_t& i) const
+  {
+    return m_ref_face_list[i];
+  }
+
+  void addRefNodeList(const RefNodeList& ref_node_list)
+  {
+    m_ref_node_list.push_back(ref_node_list);
+  }
+
+  size_t numberOfRefNodeList() const
+  {
+    return m_ref_node_list.size();
+  }
+
+  const RefNodeList& refNodeList(const size_t& i) const
+  {
+    return m_ref_node_list[i];
+  }
+
+  const size_t& numberOfNodes() const
+  {
+    return m_number_of_nodes;
+  }
+
+  const size_t& numberOfFaces() const
+  {
+    return m_number_of_faces;
+  }
+
+  const size_t& numberOfCells() const
+  {
+    return m_number_of_cells;
+  }
+
+  const size_t& maxNbNodePerCell() const
+  {
+    return m_max_nb_node_per_cell;
+  }
+
+  const Kokkos::View<const unsigned int**> cellNodes() const
+  {
+    return m_cell_nodes;
+  }
+
+  const Kokkos::View<const unsigned int**> cellFaces() const
+  {
+    return m_cell_faces;
+  }
+
+  const Kokkos::View<const unsigned short*> nodeNbCells() const
+  {
+    return m_node_nb_cells;
+  }
+
+  const Kokkos::View<const unsigned short*> cellNbNodes() const
+  {
+    return m_cell_nb_nodes;
+  }
+
+  const Kokkos::View<const double*> invCellNbNodes() const
+  {
+    return m_inv_cell_nb_nodes;
+  }
+
+  const Kokkos::View<const unsigned short*> cellNbFaces() const
+  {
+    return m_cell_nb_faces;
+  }
+
+  const Kokkos::View<const unsigned short*> faceNbCells() const
+  {
+    return m_face_nb_cells;
+  }
+
+  const Kokkos::View<const unsigned short*> faceNbNodes() const
+  {
+    return m_face_nb_nodes;
+  }
+
+  const Kokkos::View<const unsigned int**> nodeCells() const
+  {
+    return m_node_cells;
+  }
+
+  const Kokkos::View<const unsigned int**> faceCells() const
+  {
+    return m_face_cells;
+  }
+
+  const Kokkos::View<const unsigned int**> faceNodes() const
+  {
+    return m_face_nodes;
+  }
+
+  const Kokkos::View<const unsigned short**> nodeCellLocalNode() const
+  {
+    return m_node_cell_local_node;
+  }
+
+  const Kokkos::View<const unsigned short**> faceCellLocalFace() const
+  {
+    return m_face_cell_local_face;
+  }
+
+  unsigned int getFaceNumber(const unsigned int node0_id,
+                             const unsigned int node1_id) const
+  {
+#warning rewrite
+    const unsigned int n0_id = std::min(node0_id, node1_id);
+    const unsigned int n1_id = std::max(node0_id, node1_id);
+    // worst code ever
+    for (unsigned int l=0; l<m_number_of_faces; ++l) {
+      if ((m_face_nodes(l,0) == n0_id) and (m_face_nodes(l,1) == n1_id)) {
+        return l;
+      }
+    }
+    std::cerr << "Face not found!\n";
+    std::exit(0);
+    return 0;
+  }
+
+  Connectivity3D(const Connectivity2D&) = delete;
+
+  Connectivity3D(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)
+  {
+    {
+      Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", m_number_of_cells);
+      Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const int&j){
+          inv_cell_nb_nodes[j] = 1./m_cell_nb_nodes[j];
+        });
+      m_inv_cell_nb_nodes = inv_cell_nb_nodes;
+    }
+    Assert(m_number_of_cells>0);
+
+    ConnectivityUtils utils;
+    utils.computeNodeCellConnectivity(m_cell_nodes,
+                                      m_cell_nb_nodes,
+                                      m_number_of_cells,
+                                      m_max_nb_node_per_cell,
+                                      m_number_of_nodes,
+                                      m_node_nb_cells,
+                                      m_node_cells,
+                                      m_node_cell_local_node);
+
+    // this->_computeFaceCellConnectivities();
+  }
+
+  ~Connectivity3D()
+  {
+    ;
+  }
+};
+
+#endif // CONNECTIVITY_3D_HPP
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index fe528e180bf59e5770f49c97d239357decb6ea4f..2791afe32a3da5bf394272ee8a0592e57742388c 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -7,12 +7,15 @@
 
 #include <Connectivity1D.hpp>
 #include <Connectivity2D.hpp>
+#include <Connectivity3D.hpp>
+
 #include <Mesh.hpp>
 
 #include <RefFaceList.hpp>
 
 #include <map>
 #include <regex>
+#include <iomanip>
 
 template<typename T>
 inline std::string stringify(const T & t)
@@ -791,8 +794,57 @@ GmshReader::__proceedData()
   std::cout << "- dimension 2 entities : " << (dimension2_mask, elementNumber) << '\n';
   std::cout << "- dimension 3 entities : " << (dimension3_mask, elementNumber) << '\n';
   if ((dimension3_mask, elementNumber)>0) {
-    std::cerr << "*** using a 3d mesh (NIY)\n";
-    std::exit(0);
+    const size_t nb_cells = (dimension3_mask, elementNumber);
+    size_t max_nb_node_per_cell=4;
+    if (elementNumber[4] > 0) {
+      max_nb_node_per_cell = 8;
+    }
+    const Kokkos::View<unsigned int**> cell_nodes("cell_nodes", nb_cells, max_nb_node_per_cell);
+    const size_t nb_tetrahedra = __tetrahedra.extent(0);
+    for (size_t j=0; j<nb_tetrahedra; ++j) {
+      for (int r=0; r<4; ++r) {
+        cell_nodes(j,r) = __tetrahedra[j][r];
+      }
+    }
+    const size_t nb_hexahedra = __hexahedra.extent(0);
+    for (size_t j=0; j<nb_hexahedra; ++j) {
+      const size_t jh = j+nb_tetrahedra;
+      for (int r=0; r<8; ++r) {
+        cell_nodes(jh,r) = __hexahedra[j][r];
+      }
+    }
+    const Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", nb_cells);
+    for (size_t j=0; j<nb_tetrahedra; ++j) {
+      cell_nb_nodes[j] = 4;
+    }
+    for (size_t j=nb_tetrahedra; j<nb_tetrahedra+nb_hexahedra; ++j) {
+      cell_nb_nodes[j] = 8;
+    }
+
+    std::shared_ptr<Connectivity3D> p_connectivity(new Connectivity3D(cell_nb_nodes, cell_nodes));
+    Connectivity3D& connectivity = *p_connectivity;
+
+    // for (size_t j=0; j<nb_cells; ++j) {
+    //   std::cout << std::setw(5) << j << ": ";
+    //   for (size_t r=0; r<cell_nb_nodes[j]; ++r) {
+    //     std::cout << ' ' << cell_nodes(j,r);
+    //   }
+    //   std::cout << '\n';
+    // }
+
+    // std::cerr << "*** using a 3d mesh (NIY)\n";
+    // std::exit(0);
+    typedef Mesh<Connectivity3D> MeshType;
+    typedef TinyVector<3, double> Rd;
+
+    Kokkos::View<Rd*> xr("xr", __vertices.extent(0));
+    for (size_t i=0; i<__vertices.extent(0); ++i) {
+      xr[i][0] = __vertices[i][0];
+      xr[i][1] = __vertices[i][1];
+      xr[i][2] = __vertices[i][2];
+    }
+    m_mesh = std::shared_ptr<IMesh>(new MeshType(p_connectivity, xr));
+
   } else if ((dimension2_mask, elementNumber)>0) {
     const size_t nb_cells = (dimension2_mask, elementNumber);
     size_t max_nb_node_per_cell=3;
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 02b5c0a39c3883cdca519cfcb7c771fcaff37172..f65f3aa0999fd605f15836d135a80f345ba8a9c6 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -64,11 +64,11 @@ private:
     const Kokkos::View<const Rd*> xr = m_mesh.xr();
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	double sum_cjr_xr = 0;
-	for (int R=0; R<cell_nb_nodes[j]; ++R) {
-	  sum_cjr_xr += (xr[cell_nodes(j,R)], m_Cjr(j,R));
-	}
-	m_Vj[j] = inv_dimension * sum_cjr_xr;
+        double sum_cjr_xr = 0;
+        for (int R=0; R<cell_nb_nodes[j]; ++R) {
+          sum_cjr_xr += (xr[cell_nodes(j,R)], m_Cjr(j,R));
+        }
+        m_Vj[j] = inv_dimension * sum_cjr_xr;
       });
   }
 
@@ -77,40 +77,68 @@ private:
     if constexpr (dimension == 1) {
       // Cjr/njr/ljr are constant overtime
       }
-    else if (dimension == 2) {
+    else if constexpr (dimension == 2) {
       const Kokkos::View<const unsigned int**>& cell_nodes
-	= m_mesh.connectivity().cellNodes();
+          = m_mesh.connectivity().cellNodes();
       const Kokkos::View<const unsigned short*> cell_nb_nodes
-	= m_mesh.connectivity().cellNbNodes();
+          = m_mesh.connectivity().cellNbNodes();
 
       const Kokkos::View<const Rd*> xr = m_mesh.xr();
 
       Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	  for (int R=0; R<cell_nb_nodes[j]; ++R) {
-	    int Rp1 = (R+1)%cell_nb_nodes[j];
-	    int Rm1 = (R+cell_nb_nodes[j]-1)%cell_nb_nodes[j];
-	    Rd half_xrp_xrm = 0.5*(xr(cell_nodes(j,Rp1))-xr(cell_nodes(j,Rm1)));
-	    m_Cjr(j,R) = Rd{-half_xrp_xrm[1], half_xrp_xrm[0]};
-	  }
-	});
+          for (int R=0; R<cell_nb_nodes[j]; ++R) {
+            int Rp1 = (R+1)%cell_nb_nodes[j];
+            int Rm1 = (R+cell_nb_nodes[j]-1)%cell_nb_nodes[j];
+            Rd half_xrp_xrm = 0.5*(xr(cell_nodes(j,Rp1))-xr(cell_nodes(j,Rm1)));
+            m_Cjr(j,R) = Rd{-half_xrp_xrm[1], half_xrp_xrm[0]};
+          }
+        });
 
       Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	  for (int R=0; R<cell_nb_nodes[j]; ++R) {
-	    const Rd& Cjr = m_Cjr(j,R);
-	    m_ljr(j,R) = sqrt((Cjr,Cjr));
-	  }
-	});
+          for (int R=0; R<cell_nb_nodes[j]; ++R) {
+            const Rd& Cjr = m_Cjr(j,R);
+            m_ljr(j,R) = l2Norm(Cjr);
+          }
+        });
 
       Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	  for (int R=0; R<cell_nb_nodes[j]; ++R) {
-	    const Rd& Cjr = m_Cjr(j,R);
-	    const double inv_ljr = 1./m_ljr(j,R);
-	    m_njr(j,R) = inv_ljr*Cjr;
-	  }
-	});
+          for (int R=0; R<cell_nb_nodes[j]; ++R) {
+            const Rd& Cjr = m_Cjr(j,R);
+            const double inv_ljr = 1./m_ljr(j,R);
+            m_njr(j,R) = inv_ljr*Cjr;
+          }
+        });
+    } else if (dimension ==3) {
+      const Kokkos::View<const unsigned int**>& cell_nodes
+          = m_mesh.connectivity().cellNodes();
+      const Kokkos::View<const unsigned short*> cell_nb_nodes
+          = m_mesh.connectivity().cellNbNodes();
+
+      const Kokkos::View<const Rd*> xr = m_mesh.xr();
 
+      std::cerr << "Cjr are not computed in 3D!\n";
+      Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+          for (int R=0; R<cell_nb_nodes[j]; ++R) {
+            m_Cjr(j,R) = zero;
+          }
+        });
+
+      Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+          for (int R=0; R<cell_nb_nodes[j]; ++R) {
+            const Rd& Cjr = m_Cjr(j,R);
+            m_ljr(j,R) = l2Norm(Cjr);
+          }
+        });
+
+      Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+          for (int R=0; R<cell_nb_nodes[j]; ++R) {
+            const Rd& Cjr = m_Cjr(j,R);
+            const double inv_ljr = 1./m_ljr(j,R);
+            m_njr(j,R) = inv_ljr*Cjr;
+          }
+        });
     }
-    static_assert((dimension==1) or (dimension==2), "only 1d and 2d are implemented");
+    static_assert((dimension<=3), "only 1d, 2d and 3d are implemented");
   }
 
 public:
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index 070f24f22c6a5cb0d75c380b0724705938966f02..da7fc8bf2a6b1dce729dcb20c2cf153dafc1fdbb 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -46,7 +46,7 @@ public:
 	  }
 	  fout << "0 0 "; // VTK requires 3 components
 	}
-      } else if (MeshType::dimension ==2) {
+      } else if constexpr (MeshType::dimension ==2) {
 	const Kokkos::View<const TinyVector<2>*> xr = mesh.xr();
 	for (unsigned int r=0; r<mesh.numberOfNodes(); ++r) {
 	  for (unsigned short i=0; i<2; ++i) {
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index ebd715ab55ecd42904320aca22ea4f08763f3a6c..3e2d108d7eef1386aa4ec1c2679ecd26c80ec243 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -228,21 +228,10 @@ private:
   }
 
   void inverse(const Kokkos::View<const Rdd*>& A,
-	       Kokkos::View<Rdd*>& inv_A) const {
-    if constexpr (dimension==1) {
-	Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
-	    inv_A(r) = Rdd{1./(A(r)(0,0))};
-	  });
-      }
-    else {
-      Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
-	  const Rdd& M = A(r);
-	  double det = M(0,0)*M(1,1)-M(1,0)*M(0,1);
-	  double inv_det = 1./det;
-	  inv_A(r) = inv_det*Rdd(M(1,1),-M(0,1),
-				-M(1,0), M(0,0));
-	});
-    }
+               Kokkos::View<Rdd*>& inv_A) const {
+    Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
+        inv_A(r) = ::inverse(A(r));
+      });
   }
 
   void inverse(const Kokkos::View<const double*>& x,