From f686a9a260c9856e96b60b0fcedd5619fdd07ab2 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Wed, 5 Sep 2018 19:24:03 +0200
Subject: [PATCH] Use ItemValue instead of Kokkos::View where possible

---
 CMakeLists.txt                               |   6 +-
 src/CMakeLists.txt                           |   4 -
 src/experimental/AcousticSolverClass.cpp     | 381 ------------------
 src/experimental/AcousticSolverClass.hpp     |  84 ----
 src/experimental/AcousticSolverTest.cpp      | 402 -------------------
 src/experimental/AcousticSolverTest.hpp      |  96 -----
 src/experimental/CMakeLists.txt              |  17 -
 src/experimental/MeshLessAcousticSolver.cpp  | 244 -----------
 src/experimental/MeshLessAcousticSolver.hpp  |  31 --
 src/experimental/RawKokkosAcousticSolver.cpp | 248 ------------
 src/experimental/RawKokkosAcousticSolver.hpp |   9 -
 src/main.cpp                                 |  34 +-
 src/mesh/Connectivity.cpp                    |   6 +-
 src/mesh/Connectivity.hpp                    |   7 +-
 src/mesh/GmshReader.cpp                      |   6 +-
 src/mesh/ItemValue.hpp                       |   3 +
 src/mesh/Mesh.hpp                            |  29 +-
 src/mesh/MeshData.hpp                        | 213 +++++-----
 src/mesh/MeshNodeBoundary.hpp                |  14 +-
 src/mesh/SubItemValuePerItem.hpp             |   2 +-
 src/output/VTKWriter.hpp                     |  11 +-
 src/scheme/AcousticSolver.hpp                | 146 +++----
 src/scheme/BlockPerfectGas.hpp               |  68 ++--
 src/scheme/FiniteVolumesEulerUnknowns.hpp    | 107 ++---
 24 files changed, 344 insertions(+), 1824 deletions(-)
 delete mode 100644 src/experimental/AcousticSolverClass.cpp
 delete mode 100644 src/experimental/AcousticSolverClass.hpp
 delete mode 100644 src/experimental/AcousticSolverTest.cpp
 delete mode 100644 src/experimental/AcousticSolverTest.hpp
 delete mode 100644 src/experimental/CMakeLists.txt
 delete mode 100644 src/experimental/MeshLessAcousticSolver.cpp
 delete mode 100644 src/experimental/MeshLessAcousticSolver.hpp
 delete mode 100644 src/experimental/RawKokkosAcousticSolver.cpp
 delete mode 100644 src/experimental/RawKokkosAcousticSolver.hpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 341828ab7..e369d6694 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -136,8 +136,6 @@ include_directories(src/output)
 include_directories(src/utils)
 include_directories(src/scheme)
 
-include_directories(src/experimental)
-
 # Pastis tests
 
 set(CATCH_MODULE_PATH "${PASTIS_SOURCE_DIR}/packages/Catch2")
@@ -167,7 +165,6 @@ if("${CMAKE_BUILD_TYPE}" STREQUAL "Coverage")
   endif()
 
   set (GCOVR_EXCLUDE
-    -e "${PASTIS_SOURCE_DIR}/src/experimental"
     -e "${PASTIS_SOURCE_DIR}/src/main.cpp"
     -e "${PASTIS_SOURCE_DIR}/src/utils/BacktraceManager.cpp"
     -e "${PASTIS_SOURCE_DIR}/src/utils/BacktraceManager.hpp"
@@ -216,5 +213,4 @@ target_link_libraries(
   pastis
   kokkos
   PastisUtils
-  PastisMesh
-  PastisExperimental)
+  PastisMesh)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index fcc3e886e..ebf2a69d2 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -21,7 +21,3 @@ include_directories(scheme)
 
 # Pastis output
 include_directories(output)
-
-# Pastis experimental
-add_subdirectory(experimental)
-include_directories(experimental)
diff --git a/src/experimental/AcousticSolverClass.cpp b/src/experimental/AcousticSolverClass.cpp
deleted file mode 100644
index 051e148bb..000000000
--- a/src/experimental/AcousticSolverClass.cpp
+++ /dev/null
@@ -1,381 +0,0 @@
-#include <AcousticSolverClass.hpp>
-#include <rang.hpp>
-
-#include <BlockPerfectGas.hpp>
-
-typedef const double my_double;
-
-struct AcousticSolverClass::ReduceMin
-{
-private:
-  const Kokkos::View<my_double*> x_;
-
-public:
-  typedef Kokkos::View<my_double*>::non_const_value_type value_type;
-
-  ReduceMin(const Kokkos::View<my_double*>& x) : x_ (x) {}
-
-  typedef Kokkos::View<my_double*>::size_type size_type;
-    
-  KOKKOS_INLINE_FUNCTION void
-  operator() (const size_type i, value_type& update) const
-  {
-    if (x_(i) < update) {
-      update = x_(i);
-    }
-  }
-
-  KOKKOS_INLINE_FUNCTION void
-  join (volatile value_type& dst,
-	const volatile value_type& src) const
-  {
-    if (src < dst) {
-      dst = src;
-    }
-  }
-
-  KOKKOS_INLINE_FUNCTION void
-  init (value_type& dst) const
-  { // The identity under max is -Inf.
-    dst= Kokkos::reduction_identity<value_type>::min();
-  }
-};
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const double*>
-AcousticSolverClass::computeRhoCj(const Kokkos::View<const double*>& rhoj,
-				  const Kokkos::View<const double*>& cj)
-{
-  Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) {
-      m_rhocj[j] = rhoj[j]*cj[j];
-    });
-  return m_rhocj;
-}
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const double*[2]>
-AcousticSolverClass::computeAjr(const Kokkos::View<const double*>& rhocj,
-				const Kokkos::View<const double*[2]>& Cjr)
-{
-  Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) {
-      for (int r=0; r<2; ++r) {
-	m_Ajr(j,r) = rhocj(j)*Cjr(j,r)*Cjr(j,r);
-      }
-    });
-
-  return m_Ajr;
-}
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const double*>
-AcousticSolverClass::computeAr(const Kokkos::View<const double*[2]>& Ajr,
-			       const Kokkos::View<const int*[2]>& node_cells,
-			       const Kokkos::View<const int*[2]>& node_cell_local_node,
-			       const Kokkos::View<const int*>& node_nb_cells)
-{
-  Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) {
-      double sum = 0;
-      for (int j=0; j<node_nb_cells(r); ++j) {
-	const int J = node_cells(r,j);
-	const int R = node_cell_local_node(r,j);
-  	sum += Ajr(J,R);
-      }
-      m_Ar(r) = sum;
-    });
-
-  return m_Ar;
-}
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const double*>
-AcousticSolverClass::computeBr(const Kokkos::View<const double*[2]>& Ajr,
-			       const Kokkos::View<const double*[2]>& Cjr,
-			       const Kokkos::View<const double*>& uj,
-			       const Kokkos::View<const double*>& pj,
-			       const Kokkos::View<const int*[2]>& node_cells,
-			       const Kokkos::View<const int*[2]>& node_cell_local_node,
-			       const Kokkos::View<const int*>& node_nb_cells)
-{
-  Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) {
-      double sum = 0;
-      for (int j=0; j<node_nb_cells(r); ++j) {
-  	const int J = node_cells(r,j);
-  	const int R = node_cell_local_node(r,j);
-  	sum += Ajr(J,R)*uj(J) + Cjr(J,R)*pj(J);
-      }
-      m_br(r) = sum;
-    });
-
-  return m_br;
-}
-
-KOKKOS_INLINE_FUNCTION
-Kokkos::View<double*>
-AcousticSolverClass::computeUr(const Kokkos::View<const double*>& Ar,
-			       const Kokkos::View<const double*>& br)
-{
-  inverse(Ar, m_inv_Ar);
-  const Kokkos::View<const double*> invAr = m_inv_Ar;
-  Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) {
-      m_ur[r]=br(r)*invAr(r);
-    });
-  m_ur[0]=0;
-  m_ur[m_nr-1]=0;
-
-  return m_ur;
-}
-
-KOKKOS_INLINE_FUNCTION
-Kokkos::View<double*[2]>
-AcousticSolverClass::computeFjr(const Kokkos::View<const double*[2]>& Ajr,
-				const Kokkos::View<const double*>& ur,
-				const Kokkos::View<const double*[2]>& Cjr,
-				const Kokkos::View<const double*>& uj,
-				const Kokkos::View<const double*>& pj,
-				const Kokkos::View<const int*[2]>& cell_nodes)
-{
-  Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) {
-      for (int r=0; r<2; ++r) {
-	m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(j,r)))+Cjr(j,r)*pj(j);
-      }
-    });
-
-  return m_Fjr;
-}
-
-KOKKOS_INLINE_FUNCTION
-double AcousticSolverClass::
-acoustic_dt(const Kokkos::View<const double*>& Vj,
-	    const Kokkos::View<const double*>& cj) const
-{
-  Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){
-      m_Vj_over_cj[j] = Vj[j]/cj[j];
-    });
-
-  double dt = std::numeric_limits<double>::max();
-  Kokkos::parallel_reduce(m_nj, ReduceMin(m_Vj_over_cj), dt);
-
-  return dt;
-}
-
-KOKKOS_INLINE_FUNCTION
-void
-AcousticSolverClass::inverse(const Kokkos::View<const double*>& x,
-			     Kokkos::View<double*>& inv_x) const
-{
-  Kokkos::parallel_for(x.size(), KOKKOS_LAMBDA(const int& r) {
-      inv_x(r) = 1./x(r);
-    });
-}
-
-KOKKOS_INLINE_FUNCTION
-void AcousticSolverClass::computeExplicitFluxes(const Kokkos::View<const double*>& xr,
-						const Kokkos::View<const double*>& xj,
-						const Kokkos::View<const double*>& rhoj,
-						const Kokkos::View<const double*>& uj,
-						const Kokkos::View<const double*>& pj,
-						const Kokkos::View<const double*>& cj,
-						const Kokkos::View<const double*>& Vj,
-						const Kokkos::View<const double*[2]>& Cjr,
-						const Kokkos::View<const int*[2]>& cell_nodes,
-						const Kokkos::View<const int*[2]>& node_cells,
-						const Kokkos::View<const int*>& node_nb_cells,
-						const Kokkos::View<const int*[2]>& node_cell_local_node,
-						Kokkos::View<double*>& ur,
-						Kokkos::View<double*[2]>& Fjr)
-{
-  const Kokkos::View<const double*> rhocj  = computeRhoCj(rhoj, cj);
-  const Kokkos::View<const double*[2]> Ajr = computeAjr(rhocj, Cjr);
-
-  const Kokkos::View<const double*> Ar = computeAr(Ajr, node_cells, node_cell_local_node, node_nb_cells);
-  const Kokkos::View<const double*> br = computeBr(Ajr, Cjr, uj, pj,
-						   node_cells, node_cell_local_node, node_nb_cells);
-
-  ur  = computeUr(Ar, br);
-  Fjr = computeFjr(Ajr, ur, Cjr, uj, pj, cell_nodes);
-}
-
-AcousticSolverClass::AcousticSolverClass(const long int& nj)
-  : m_nj(nj),
-    m_nr(nj+1),
-    m_br("br", m_nr),
-    m_Ajr("Ajr", m_nj),
-    m_Ar("Ar", m_nr),
-    m_inv_Ar("inv_Ar", m_nr),
-    m_Fjr("Fjr", m_nj),
-    m_ur("ur", m_nr),
-    m_rhocj("rho_c", m_nj),
-    m_Vj_over_cj("Vj_over_cj", m_nj)
-{
-  Kokkos::View<double*> xj("xj",m_nj);
-  Kokkos::View<double*> rhoj("rhoj",m_nj);
-
-  Kokkos::View<double*> uj("uj",m_nj);
-
-  Kokkos::View<double*> Ej("Ej",m_nj);
-  Kokkos::View<double*> ej("ej",m_nj);
-  Kokkos::View<double*> pj("pj",m_nj);
-  Kokkos::View<double*> Vj("Vj",m_nj);
-  Kokkos::View<double*> gammaj("gammaj",m_nj);
-  Kokkos::View<double*> cj("cj",m_nj);
-  Kokkos::View<double*> mj("mj",m_nj);
-
-  Kokkos::View<double*>  xr("xr", m_nr);
-
-  const double delta_x = 1./m_nj;
-
-  Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
-      xr[r] = r*delta_x;
-    });
-
-  Kokkos::View<int*[2]> cell_nodes("cell_nodes",m_nj,2);
-  Kokkos::View<int*[2]> node_cells("node_cells",m_nr,2);
-  Kokkos::View<int*[2]> node_cell_local_node("node_cells",m_nr,2);
-  Kokkos::View<int*> node_nb_cells("node_cells",m_nr);
-
-  Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
-      node_nb_cells(r) = 2;
-    });
-  node_nb_cells(0) = 1;
-  node_nb_cells(m_nr-1) = 1;
-
-  node_cells(0,0) = 0;
-  Kokkos::parallel_for(m_nr-2, KOKKOS_LAMBDA(const int& r){
-      node_cells(r+1,0) = r;
-      node_cells(r+1,1) = r+1;
-    });
-  node_cells(m_nr-1,0) =m_nj-1;
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      cell_nodes(j,0) = j;
-      cell_nodes(j,1) = j+1;
-    });
-
-  Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
-      for (int J=0; J<node_nb_cells(r); ++J) {
-	int j = node_cells(r,J);
-	for (int R=0; R<2; ++R) {
-	  if (cell_nodes(j,R) == r) {
-	    node_cell_local_node(r,J) = R;
-	  }
-	}
-      }
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      xj[j] = 0.5*(xr[cell_nodes(j,0)]+xr[cell_nodes(j,1)]);
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      Vj[j] = xr[cell_nodes(j,1)]-xr[cell_nodes(j,0)];
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      if (xj[j]<0.5) {
-	rhoj[j]=1;
-      } else {
-	rhoj[j]=0.125;
-      }
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      if (xj[j]<0.5) {
-	pj[j]=1;
-      } else {
-	pj[j]=0.1;
-      }
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      gammaj[j] = 1.4;
-    });
-
-  BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
-
-  block_eos.updateEandCFromRhoP();
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      Ej[j] = ej[j]+0.5*uj[j]*uj[j];
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      mj[j] = rhoj[j] * Vj[j];
-    });
-
-  Kokkos::View<double*> inv_mj("inv_mj",m_nj);
-  inverse(mj, inv_mj);
-
-  const double tmax=0.2;
-  double t=0;
-
-  int itermax=std::numeric_limits<int>::max();
-  int iteration=0;
-
-  Kokkos::View<double*[2]> Cjr("Cjr",m_nj);
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
-      Cjr(j,0)=-1;
-      Cjr(j,1)= 1;
-    });
-
-  while((t<tmax) and (iteration<itermax)) {
-    double dt = 0.4*acoustic_dt(Vj, cj);
-    if (t+dt<tmax) {
-      t+=dt;
-    } else {
-      dt=tmax-t;
-      t=tmax;
-    }
-
-    computeExplicitFluxes(xr, xj,
-			  rhoj, uj, pj, cj, Vj, Cjr,
-			  cell_nodes, node_cells, node_nb_cells, node_cell_local_node,
-			  m_ur, m_Fjr);
-
-    const Kokkos::View<const double*[2]> Fjr = m_Fjr;
-    const Kokkos::View<const double*> ur = m_ur;
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-	double momentum_fluxes = 0;
-	double energy_fluxes = 0;
-	for (int R=0; R<2; ++R) {
-	  const int r=cell_nodes(j,R);
-	  momentum_fluxes += Fjr(j,R);
-	  energy_fluxes += Fjr(j,R)*ur[r];
-	}
-	uj[j] -= dt*inv_mj[j]*(momentum_fluxes);
-	Ej[j] -= dt*inv_mj[j]*(energy_fluxes);
-      });
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
-	ej[j] = Ej[j] - 0.5 * uj[j]*uj[j];
-      });
-
-    Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
-	xr[r] += dt*ur[r];
-      });
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-	xj[j] = 0.5*(xr[j]+xr[j+1]);
-	Vj[j] = xr[j+1]-xr[j];
-      });
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-	rhoj[j] = mj[j]/Vj[j];
-      });
-
-    block_eos.updatePandCFromRhoE();    
-    
-    ++iteration;
-  }
-
-  // {
-  //   std::ofstream fout("rho");
-  //   for (int j=0; j<nj; ++j) {
-  //     fout << xj[j] << ' ' << rhoj[j] << '\n';
-  //   }
-  // }
-
-  std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
-	    << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-
-}
diff --git a/src/experimental/AcousticSolverClass.hpp b/src/experimental/AcousticSolverClass.hpp
deleted file mode 100644
index 2c2f8c0ec..000000000
--- a/src/experimental/AcousticSolverClass.hpp
+++ /dev/null
@@ -1,84 +0,0 @@
-#ifndef ACOUSTIC_SOLVER_CLASS_HPP
-#define ACOUSTIC_SOLVER_CLASS_HPP
-
-#include <Kokkos_Core.hpp>
-
-class AcousticSolverClass
-{
-private:
-  inline const Kokkos::View<const double*>
-  computeRhoCj(const Kokkos::View<const double*>& rhoj,
-	       const Kokkos::View<const double*>& cj);
-
-  inline const Kokkos::View<const double*[2]>
-  computeAjr(const Kokkos::View<const double*>& rhocj,
-	     const Kokkos::View<const double*[2]>& Cjr);
-
-  inline const Kokkos::View<const double*>
-  computeAr(const Kokkos::View<const double*[2]>& Ajr,
-	    const Kokkos::View<const int*[2]>& node_cells,
-	    const Kokkos::View<const int*[2]>& node_cell_local_node,
-	    const Kokkos::View<const int*>& node_nb_cells);
-
-  inline const Kokkos::View<const double*>
-  computeBr(const Kokkos::View<const double*[2]>& Ajr,
-	    const Kokkos::View<const double*[2]>& Cjr,
-	    const Kokkos::View<const double*>& uj,
-	    const Kokkos::View<const double*>& pj,
-	    const Kokkos::View<const int*[2]>& node_cells,
-	    const Kokkos::View<const int*[2]>& node_cell_local_node,
-	    const Kokkos::View<const int*>& node_nb_cells);
-
-  Kokkos::View<double*>
-  computeUr(const Kokkos::View<const double*>& Ar,
-	    const Kokkos::View<const double*>& br);
-
-  Kokkos::View<double*[2]>
-  computeFjr(const Kokkos::View<const double*[2]>& Ajr,
-	     const Kokkos::View<const double*>& ur,
-	     const Kokkos::View<const double*[2]>& Cjr,
-	     const Kokkos::View<const double*>& uj,
-	     const Kokkos::View<const double*>& pj,
-	     const Kokkos::View<const int*[2]>& cell_nodes);
-
-  void inverse(const Kokkos::View<const double*>& x,
-	       Kokkos::View<double*>& inv_x) const;
-
-  inline double acoustic_dt(const Kokkos::View<const double*>& Vj,
-			    const Kokkos::View<const double*>& cj) const;
-
-  inline void computeExplicitFluxes(const Kokkos::View<const double*>& xr,
-				    const Kokkos::View<const double*>& xj,
-				    const Kokkos::View<const double*>& rhoj,
-				    const Kokkos::View<const double*>& uj,
-				    const Kokkos::View<const double*>& pj,
-				    const Kokkos::View<const double*>& cj,
-				    const Kokkos::View<const double*>& Vj,
-				    const Kokkos::View<const double*[2]>& Cjr,
-				    const Kokkos::View<const int*[2]>& cell_nodes,
-				    const Kokkos::View<const int*[2]>& node_cells,
-				    const Kokkos::View<const int*>& node_nb_cells,
-				    const Kokkos::View<const int*[2]>& node_cell_local_node,
-				    Kokkos::View<double*>& ur,
-				    Kokkos::View<double*[2]>& Fjr);
-
-  struct ReduceMin;
-
-  const long int m_nj;
-  const long int m_nr;
-
-  Kokkos::View<double*> m_br;
-  Kokkos::View<double*[2]> m_Ajr;
-  Kokkos::View<double*> m_Ar;
-  Kokkos::View<double*> m_inv_Ar;
-  Kokkos::View<double*[2]> m_Fjr;
-  Kokkos::View<double*> m_ur;
-  Kokkos::View<double*> m_rhocj;
-  Kokkos::View<double*> m_Vj_over_cj;
-
-
-public:
-  AcousticSolverClass(const long int& nj);
-};
-
-#endif // ACOUSTIC_SOLVER_CLASS_HPP
diff --git a/src/experimental/AcousticSolverTest.cpp b/src/experimental/AcousticSolverTest.cpp
deleted file mode 100644
index aaa9c0dbf..000000000
--- a/src/experimental/AcousticSolverTest.cpp
+++ /dev/null
@@ -1,402 +0,0 @@
-#include <AcousticSolverTest.hpp>
-#include <rang.hpp>
-
-#include <BlockPerfectGas.hpp>
-#include <fstream>
-
-typedef const double my_double;
-
-struct AcousticSolverTest::ReduceMin
-{
-private:
-  const Kokkos::View<my_double*> x_;
-
-public:
-  typedef Kokkos::View<my_double*>::non_const_value_type value_type;
-
-  ReduceMin(const Kokkos::View<my_double*>& x) : x_ (x) {}
-
-  typedef Kokkos::View<my_double*>::size_type size_type;
-
-  KOKKOS_INLINE_FUNCTION void
-  operator() (const size_type i, value_type& update) const
-  {
-    if (x_(i) < update) {
-      update = x_(i);
-    }
-  }
-
-  KOKKOS_INLINE_FUNCTION void
-  join (volatile value_type& dst,
-	const volatile value_type& src) const
-  {
-    if (src < dst) {
-      dst = src;
-    }
-  }
-
-  KOKKOS_INLINE_FUNCTION void
-  init (value_type& dst) const
-  { // The identity under max is -Inf.
-    dst= Kokkos::reduction_identity<value_type>::min();
-  }
-};
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const double*>
-AcousticSolverTest::computeRhoCj(const Kokkos::View<const double*>& rhoj,
-				 const Kokkos::View<const double*>& cj)
-{
-  Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) {
-      m_rhocj[j] = rhoj[j]*cj[j];
-    });
-  return m_rhocj;
-}
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const Rdd*[2]>
-AcousticSolverTest::computeAjr(const Kokkos::View<const double*>& rhocj,
-			       const Kokkos::View<const Rd*[2]>& Cjr)
-{
-  Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) {
-      for (int r=0; r<2; ++r) {
-	m_Ajr(j,r) = tensorProduct(rhocj(j)*Cjr(j,r), Cjr(j,r));
-      }
-    });
-
-  return m_Ajr;
-}
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const Rdd*>
-AcousticSolverTest::computeAr(const Kokkos::View<const Rdd*[2]>& Ajr,
-			      const Kokkos::View<const int*[2]>& node_cells,
-			      const Kokkos::View<const int*[2]>& node_cell_local_node,
-			      const Kokkos::View<const int*>& node_nb_cells)
-{
-  Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) {
-      Rdd sum = zero;
-      for (int j=0; j<node_nb_cells(r); ++j) {
-	const int J = node_cells(r,j);
-	const int R = node_cell_local_node(r,j);
-  	sum += Ajr(J,R);
-      }
-      m_Ar(r) = sum;
-    });
-
-  return m_Ar;
-}
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const Rd*>
-AcousticSolverTest::computeBr(const Kokkos::View<const Rdd*[2]>& Ajr,
-			      const Kokkos::View<const Rd*[2]>& Cjr,
-			      const Kokkos::View<const Rd*>& uj,
-			      const Kokkos::View<const double*>& pj,
-			      const Kokkos::View<const int*[2]>& node_cells,
-			      const Kokkos::View<const int*[2]>& node_cell_local_node,
-			      const Kokkos::View<const int*>& node_nb_cells)
-{
-  Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) {
-      Rd& br = m_br(r);
-      br = zero;
-      for (int j=0; j<node_nb_cells(r); ++j) {
-  	const int J = node_cells(r,j);
-  	const int R = node_cell_local_node(r,j);
-  	br += Ajr(J,R)*uj(J) + pj(J)*Cjr(J,R);
-      }
-    });
-
-  return m_br;
-}
-
-KOKKOS_INLINE_FUNCTION
-Kokkos::View<Rd*>
-AcousticSolverTest::computeUr(const Kokkos::View<const Rdd*>& Ar,
-			      const Kokkos::View<const Rd*>& br)
-{
-  inverse(Ar, m_inv_Ar);
-  const Kokkos::View<const Rdd*> invAr = m_inv_Ar;
-  Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) {
-      m_ur[r]=invAr(r)*br(r);
-    });
-  m_ur[0]=zero;
-  m_ur[m_nr-1]=zero;
-
-  return m_ur;
-}
-
-KOKKOS_INLINE_FUNCTION
-Kokkos::View<Rd*[2]>
-AcousticSolverTest::computeFjr(const Kokkos::View<const Rdd*[2]>& Ajr,
-			       const Kokkos::View<const Rd*>& ur,
-			       const Kokkos::View<const Rd*[2]>& Cjr,
-			       const Kokkos::View<const Rd*>& uj,
-			       const Kokkos::View<const double*>& pj,
-			       const Kokkos::View<const int*[2]>& cell_nodes)
-{
-  Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) {
-      for (int r=0; r<2; ++r) {
-	m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(j,r)))+pj(j)*Cjr(j,r);
-      }
-    });
-
-  return m_Fjr;
-}
-
-KOKKOS_INLINE_FUNCTION
-double AcousticSolverTest::
-acoustic_dt(const Kokkos::View<const double*>& Vj,
-	    const Kokkos::View<const double*>& cj) const
-{
-  Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){
-      m_Vj_over_cj[j] = Vj[j]/cj[j];
-    });
-
-  double dt = std::numeric_limits<double>::max();
-  Kokkos::parallel_reduce(m_nj, ReduceMin(m_Vj_over_cj), dt);
-
-  return dt;
-}
-
-KOKKOS_INLINE_FUNCTION
-void
-AcousticSolverTest::inverse(const Kokkos::View<const double*>& x,
-			    Kokkos::View<double*>& inv_x) const
-{
-  Kokkos::parallel_for(x.size(), KOKKOS_LAMBDA(const int& r) {
-      inv_x(r) = 1./x(r);
-    });
-}
-
-KOKKOS_INLINE_FUNCTION
-void
-AcousticSolverTest::inverse(const Kokkos::View<const Rdd*>& A,
-			    Kokkos::View<Rdd*>& inv_A) const
-{
-  Kokkos::parallel_for(A.size(), KOKKOS_LAMBDA(const int& r) {
-      inv_A(r) = Rdd{1./(A(r)(0,0))};
-    });
-}
-
-
-KOKKOS_INLINE_FUNCTION
-void AcousticSolverTest::computeExplicitFluxes(const Kokkos::View<const Rd*>& xr,
-					       const Kokkos::View<const Rd*>& xj,
-					       const Kokkos::View<const double*>& rhoj,
-					       const Kokkos::View<const Rd*>& uj,
-					       const Kokkos::View<const double*>& pj,
-					       const Kokkos::View<const double*>& cj,
-					       const Kokkos::View<const double*>& Vj,
-					       const Kokkos::View<const Rd*[2]>& Cjr,
-					       const Kokkos::View<const int*[2]>& cell_nodes,
-					       const Kokkos::View<const int*[2]>& node_cells,
-					       const Kokkos::View<const int*>& node_nb_cells,
-					       const Kokkos::View<const int*[2]>& node_cell_local_node,
-					       Kokkos::View<Rd*>& ur,
-					       Kokkos::View<Rd*[2]>& Fjr)
-{
-  const Kokkos::View<const double*> rhocj  = computeRhoCj(rhoj, cj);
-  const Kokkos::View<const Rdd*[2]> Ajr = computeAjr(rhocj, Cjr);
-
-  const Kokkos::View<const Rdd*> Ar = computeAr(Ajr, node_cells, node_cell_local_node, node_nb_cells);
-  const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, pj,
-						   node_cells, node_cell_local_node, node_nb_cells);
-
-  ur  = computeUr(Ar, br);
-  Fjr = computeFjr(Ajr, ur, Cjr, uj, pj, cell_nodes);
-}
-
-AcousticSolverTest::AcousticSolverTest(const long int& nj)
-  : m_nj(nj),
-    m_nr(nj+1),
-    m_br("br", m_nr),
-    m_Ajr("Ajr", m_nj),
-    m_Ar("Ar", m_nr),
-    m_inv_Ar("inv_Ar", m_nr),
-    m_Fjr("Fjr", m_nj),
-    m_ur("ur", m_nr),
-    m_rhocj("rho_c", m_nj),
-    m_Vj_over_cj("Vj_over_cj", m_nj)
-{
-  Kokkos::View<Rd*> xj("xj",m_nj);
-  Kokkos::View<double*> rhoj("rhoj",m_nj);
-
-  Kokkos::View<Rd*> uj("uj",m_nj);
-
-  Kokkos::View<double*> Ej("Ej",m_nj);
-  Kokkos::View<double*> ej("ej",m_nj);
-  Kokkos::View<double*> pj("pj",m_nj);
-  Kokkos::View<double*> Vj("Vj",m_nj);
-  Kokkos::View<double*> gammaj("gammaj",m_nj);
-  Kokkos::View<double*> cj("cj",m_nj);
-  Kokkos::View<double*> mj("mj",m_nj);
-
-  Kokkos::View<Rd*>  xr("xr", m_nr);
-
-  const double delta_x = 1./m_nj;
-
-  Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
-      xr[r][0] = r*delta_x;
-    });
-
-  Kokkos::View<int*[2]> cell_nodes("cell_nodes",m_nj,2);
-  Kokkos::View<int*[2]> node_cells("node_cells",m_nr,2);
-  Kokkos::View<int*[2]> node_cell_local_node("node_cells",m_nr,2);
-  Kokkos::View<int*> node_nb_cells("node_cells",m_nr);
-
-  Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
-      node_nb_cells(r) = 2;
-    });
-  node_nb_cells(0) = 1;
-  node_nb_cells(m_nr-1) = 1;
-
-  node_cells(0,0) = 0;
-  Kokkos::parallel_for(m_nr-2, KOKKOS_LAMBDA(const int& r){
-      node_cells(r+1,0) = r;
-      node_cells(r+1,1) = r+1;
-    });
-  node_cells(m_nr-1,0) =m_nj-1;
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      cell_nodes(j,0) = j;
-      cell_nodes(j,1) = j+1;
-    });
-
-  Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
-      for (int J=0; J<node_nb_cells(r); ++J) {
-	int j = node_cells(r,J);
-	for (int R=0; R<2; ++R) {
-	  if (cell_nodes(j,R) == r) {
-	    node_cell_local_node(r,J) = R;
-	  }
-	}
-      }
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      xj[j] = 0.5*(xr[cell_nodes(j,0)]+xr[cell_nodes(j,1)]);
-    });
-
-  Kokkos::View<Rd*[2]> Cjr("Cjr",m_nj);
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
-      Cjr(j,0)=-1;
-      Cjr(j,1)= 1;
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      Vj[j] = (xr[cell_nodes(j,1)], Cjr(j,1))
-	    + (xr[cell_nodes(j,0)], Cjr(j,0));
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      if (xj[j][0]<0.5) {
-	rhoj[j]=1;
-      } else {
-	rhoj[j]=0.125;
-      }
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      if (xj[j][0]<0.5) {
-	pj[j]=1;
-      } else {
-	pj[j]=0.1;
-      }
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      uj[j] = zero;
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      gammaj[j] = 1.4;
-    });
-
-  BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
-
-  block_eos.updateEandCFromRhoP();
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      Ej[j] = ej[j]+0.5*(uj[j],uj[j]);
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      mj[j] = rhoj[j] * Vj[j];
-    });
-
-  Kokkos::View<double*> inv_mj("inv_mj",m_nj);
-  inverse(mj, inv_mj);
-
-  const double tmax=0.2;
-  double t=0;
-
-  int itermax=std::numeric_limits<int>::max();
-  int iteration=0;
-
-  while((t<tmax) and (iteration<itermax)) {
-    double dt = 0.4*acoustic_dt(Vj, cj);
-    if (t+dt<tmax) {
-      t+=dt;
-    } else {
-      dt=tmax-t;
-      t=tmax;
-    }
-
-    computeExplicitFluxes(xr, xj,
-			  rhoj, uj, pj, cj, Vj, Cjr,
-			  cell_nodes, node_cells, node_nb_cells, node_cell_local_node,
-			  m_ur, m_Fjr);
-
-    const Kokkos::View<const Rd*[2]> Fjr = m_Fjr;
-    const Kokkos::View<const Rd*> ur = m_ur;
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
-	Rd momentum_fluxes = zero;
-	double energy_fluxes = 0;
-	for (int R=0; R<2; ++R) {
-	  const int r=cell_nodes(j,R);
-	  momentum_fluxes +=  Fjr(j,R);
-	  energy_fluxes   += ((Fjr(j,R), ur[r]));
-	}
-	uj[j] -= (dt*inv_mj[j]) * momentum_fluxes;
-	Ej[j] -= (dt*inv_mj[j]) * energy_fluxes;
-      });
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
-	ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
-      });
-
-    Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
-	xr[r] += dt*ur[r];
-      });
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-	xj[j] = 0.5*(xr[j]+xr[j+1]);
-      });
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      Vj[j] = (xr[cell_nodes(j,1)], Cjr(j,1))
-	    + (xr[cell_nodes(j,0)], Cjr(j,0));
-    });
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-	rhoj[j] = mj[j]/Vj[j];
-      });
-
-    block_eos.updatePandCFromRhoE();
-
-    ++iteration;
-  }
-
-  {
-    std::ofstream fout("rho");
-    for (int j=0; j<nj; ++j) {
-      fout << xj[j][0] << ' ' << rhoj[j] << '\n';
-    }
-  }
-
-  std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
-	    << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-
-}
diff --git a/src/experimental/AcousticSolverTest.hpp b/src/experimental/AcousticSolverTest.hpp
deleted file mode 100644
index 209a18747..000000000
--- a/src/experimental/AcousticSolverTest.hpp
+++ /dev/null
@@ -1,96 +0,0 @@
-#ifndef ACOUSTIC_SOLVER_TEST_HPP
-#define ACOUSTIC_SOLVER_TEST_HPP
-
-#include <Kokkos_Core.hpp>
-
-#include <TinyVector.hpp>
-#include <TinyMatrix.hpp>
-
-constexpr static size_t dimension = 1;
-
-typedef TinyVector<dimension> Rd;
-typedef TinyMatrix<dimension> Rdd;
-
-class AcousticSolverTest 
-{
-private:
-
-  inline const Kokkos::View<const double*>
-  computeRhoCj(const Kokkos::View<const double*>& rhoj,
-	       const Kokkos::View<const double*>& cj);
-
-  inline const Kokkos::View<const Rdd*[2]>
-  computeAjr(const Kokkos::View<const double*>& rhocj,
-	     const Kokkos::View<const Rd*[2]>& Cjr);
-
-  inline const Kokkos::View<const Rdd*>
-  computeAr(const Kokkos::View<const Rdd*[2]>& Ajr,
-	    const Kokkos::View<const int*[2]>& node_cells,
-	    const Kokkos::View<const int*[2]>& node_cell_local_node,
-	    const Kokkos::View<const int*>& node_nb_cells);
-
-  inline const Kokkos::View<const Rd*>
-  computeBr(const Kokkos::View<const Rdd*[2]>& Ajr,
-	    const Kokkos::View<const Rd*[2]>& Cjr,
-	    const Kokkos::View<const Rd*>& uj,
-	    const Kokkos::View<const double*>& pj,
-	    const Kokkos::View<const int*[2]>& node_cells,
-	    const Kokkos::View<const int*[2]>& node_cell_local_node,
-	    const Kokkos::View<const int*>& node_nb_cells);
-
-  Kokkos::View<Rd*>
-  computeUr(const Kokkos::View<const Rdd*>& Ar,
-	    const Kokkos::View<const Rd*>& br);
-
-  Kokkos::View<Rd*[2]>
-  computeFjr(const Kokkos::View<const Rdd*[2]>& Ajr,
-	     const Kokkos::View<const Rd*>& ur,
-	     const Kokkos::View<const Rd*[2]>& Cjr,
-	     const Kokkos::View<const Rd*>& uj,
-	     const Kokkos::View<const double*>& pj,
-	     const Kokkos::View<const int*[2]>& cell_nodes);
-
-  void inverse(const Kokkos::View<const Rdd*>& A,
-	       Kokkos::View<Rdd*>& inv_A) const;
-
-  void inverse(const Kokkos::View<const double*>& x,
-	       Kokkos::View<double*>& inv_x) const;
-
-  inline double acoustic_dt(const Kokkos::View<const double*>& Vj,
-			    const Kokkos::View<const double*>& cj) const;
-
-  inline void computeExplicitFluxes(const Kokkos::View<const Rd*>& xr,
-				    const Kokkos::View<const Rd*>& xj,
-				    const Kokkos::View<const double*>& rhoj,
-				    const Kokkos::View<const Rd*>& uj,
-				    const Kokkos::View<const double*>& pj,
-				    const Kokkos::View<const double*>& cj,
-				    const Kokkos::View<const double*>& Vj,
-				    const Kokkos::View<const Rd*[2]>& Cjr,
-				    const Kokkos::View<const int*[2]>& cell_nodes,
-				    const Kokkos::View<const int*[2]>& node_cells,
-				    const Kokkos::View<const int*>& node_nb_cells,
-				    const Kokkos::View<const int*[2]>& node_cell_local_node,
-				    Kokkos::View<Rd*>& ur,
-				    Kokkos::View<Rd*[2]>& Fjr);
-
-  struct ReduceMin;
-
-  const size_t m_nj;
-  const size_t m_nr;
-
-  Kokkos::View<Rd*> m_br;
-  Kokkos::View<Rdd*[2]> m_Ajr;
-  Kokkos::View<Rdd*> m_Ar;
-  Kokkos::View<Rdd*> m_inv_Ar;
-  Kokkos::View<Rd*[2]> m_Fjr;
-  Kokkos::View<Rd*> m_ur;
-  Kokkos::View<double*> m_rhocj;
-  Kokkos::View<double*> m_Vj_over_cj;
-
-
-public:
-  AcousticSolverTest(const long int& nj);
-};
-
-#endif // ACOUSTIC_SOLVER_TEST_HPP
diff --git a/src/experimental/CMakeLists.txt b/src/experimental/CMakeLists.txt
deleted file mode 100644
index 1c6043787..000000000
--- a/src/experimental/CMakeLists.txt
+++ /dev/null
@@ -1,17 +0,0 @@
-include_directories(${CMAKE_CURRENT_SOURCE_DIR})
-include_directories(${CMAKE_CURRENT_BINARY_DIR})
-
-include_directories(${Kokkos_INCLUDE_DIRS_RET})
-include_directories(${PASTIS_SOURCE_DIR}/packages/rang/include)
-
-# ------------------- Source files --------------------
-
-add_library(
-  PastisExperimental
-  AcousticSolverClass.cpp
-  AcousticSolverTest.cpp
-  RawKokkosAcousticSolver.cpp
-  MeshLessAcousticSolver.cpp
-  )
-
-
diff --git a/src/experimental/MeshLessAcousticSolver.cpp b/src/experimental/MeshLessAcousticSolver.cpp
deleted file mode 100644
index bf5014481..000000000
--- a/src/experimental/MeshLessAcousticSolver.cpp
+++ /dev/null
@@ -1,244 +0,0 @@
-#include <MeshLessAcousticSolver.hpp>
-#include <rang.hpp>
-
-#include <memory>
-
-#include <BlockPerfectGas.hpp>
-
-typedef const double my_double;
-
-struct MeshLessAcousticSolver::ReduceMin {
-private:
-  const Kokkos::View<my_double*> x_;
-
-public:
-  typedef Kokkos::View<my_double*>::non_const_value_type value_type;
-
-  ReduceMin(const Kokkos::View<my_double*>& x) : x_ (x) {}
-
-  typedef Kokkos::View<my_double*>::size_type size_type;
-    
-  KOKKOS_INLINE_FUNCTION void
-  operator() (const size_type i, value_type& update) const
-  {
-    if (x_(i) < update) {
-      update = x_(i);
-    }
-  }
-
-  KOKKOS_INLINE_FUNCTION void
-  join (volatile value_type& dst,
-	const volatile value_type& src) const
-  {
-    if (src < dst) {
-      dst = src;
-    }
-  }
-
-  KOKKOS_INLINE_FUNCTION void
-  init (value_type& dst) const
-  { // The identity under max is -Inf.
-    dst= Kokkos::reduction_identity<value_type>::min();
-  }
-};
-    
-
-KOKKOS_INLINE_FUNCTION
-double MeshLessAcousticSolver::
-acoustic_dt(const Kokkos::View<const double*>& Vj,
-	    const Kokkos::View<const double*>& cj) const
-{
-  const size_t nj = Vj.size();
-  double dt = std::numeric_limits<double>::max();
-
-  Kokkos::View<double*> Vj_cj("Vj_cj", nj);
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      Vj_cj[j] = Vj[j]/cj[j];
-    });
-
-  Kokkos::parallel_reduce(nj, ReduceMin(Vj_cj), dt);
-
-  return dt;
-}
-
-
-KOKKOS_INLINE_FUNCTION
-void MeshLessAcousticSolver::computeExplicitFluxes(const Kokkos::View<const double*>& xr,
-					   const Kokkos::View<const double*>& xj,
-					   const Kokkos::View<const double*>& rhoj,
-					   const Kokkos::View<const double*>& uj,
-					   const Kokkos::View<const double*>& pj,
-					   const Kokkos::View<const double*>& cj,
-					   const Kokkos::View<const double*>& Vj,
-					   Kokkos::View<double*>& ur,
-					   Kokkos::View<double*>& pr) const
-{
-  // calcul de ur
-  ur[0]=0;
-  const size_t nr = ur.size();
-  const size_t nj = uj.size();
-
-  Kokkos::parallel_for(nj-1, KOKKOS_LAMBDA(const int& j) {
-    const int r = j+1;
-    const int k = r;
-
-    const double ujr = uj[j];
-    const double ukr = uj[k];
-    const double pjr = pj[j];
-    const double pkr = pj[k];
-
-    ur[r]=(rhoj[j]*cj[j]*ujr + rhoj[k]*cj[k]*ukr + pjr-pkr)/(rhoj[j]*cj[j]+rhoj[k]*cj[k]);
-    });
-  ur[nr-1]=0;
-
-  // calcul de pr
-  pr[0] = pj[0] + rhoj[0]*cj[0]*(ur[0] - uj[0]);
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
-    const int r = j+1;
-
-    const double ujr = uj[j];
-    const double pjr = pj[j];
-
-    pr[r]=pjr+rhoj[j]*cj[j]*(ujr-ur[r]);
-    });
-}
-
-
-MeshLessAcousticSolver::MeshLessAcousticSolver(const long int& nj)
-{
-  Kokkos::View<double*> xj("xj", nj);
-  Kokkos::View<double*> rhoj("rhoj", nj);
-
-  Kokkos::View<double*> uj("uj", nj);
-
-  Kokkos::View<double*> Ej("Ej", nj);
-  Kokkos::View<double*> ej("ej", nj);
-  Kokkos::View<double*> pj("pj", nj);
-  Kokkos::View<double*> Vj("Vj", nj);
-  Kokkos::View<double*> gammaj("gammaj", nj);
-  Kokkos::View<double*> cj("cj", nj);
-  Kokkos::View<double*> mj("mj", nj);
-  Kokkos::View<double*> inv_mj("inv_mj", nj);
-
-  const int nr=nj+1;
-
-  Kokkos::View<double*>  xr("xr", nr);
-
-  const double delta_x = 1./nj;
-
-  Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){
-      xr[r] = r*delta_x;
-    });
-
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      xj[j] = 0.5*(xr[j]+xr[j+1]);
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      Vj[j] = xr[j+1]-xr[j];
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-    if (xj[j]<0.5) {
-      rhoj[j]=1;
-    } else {
-      rhoj[j]=0.125;
-    }
-  });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-    if (xj[j]<0.5) {
-      pj[j]=1;
-    } else {
-      pj[j]=0.1;
-    }
-  });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      gammaj[j] = 1.4;
-    });
-
-  BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
-
-  block_eos.updateEandCFromRhoP();
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      Ej[j] = ej[j]+0.5*uj[j]*uj[j];
-  });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-    mj[j] = rhoj[j] * Vj[j];
-  });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      inv_mj[j] = 1./mj[j];
-    });
-
-  const double tmax=0.2;
-  double t=0;
-
-  int itermax=std::numeric_limits<int>::max();
-  int iteration=0;
-
-  while((t<tmax) and (iteration<itermax)) {
-    double dt = 0.4*acoustic_dt(Vj, cj);
-    if (t+dt<tmax) {
-      t+=dt;
-    } else {
-      dt=tmax-t;
-      t=tmax;
-    }
-    
-    Kokkos::View<double*> ur("ur", nr);
-    Kokkos::View<double*> pr("pr", nr);
-
-    computeExplicitFluxes(xr, xj,
-			  rhoj, uj, pj, cj, Vj,
-			  ur, pr);
-
-    Kokkos::View<double*> new_uj("new_uj", nj);
-    Kokkos::View<double*> new_Ej("new_Ej", nj);
-    
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      int rm=j;
-      int rp=j+1;
-
-      new_uj[j] = uj[j] + dt*inv_mj[j]*(pr[rm]-pr[rp]);
-      new_Ej[j] = Ej[j] + dt*inv_mj[j]*(pr[rm]*ur[rm]-pr[rp]*ur[rp]);
-      });
-
-    uj=new_uj;
-    Ej=new_Ej;
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-	ej[j] = Ej[j] - 0.5 * uj[j]*uj[j];
-      });
-
-    Kokkos::View<double*> xr_new("new_xr", nr);
-
-    Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){
-      xr_new[r] = xr[r] + dt*ur[r];
-      });
-
-    xr = xr_new;
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      xj[j] = 0.5*(xr[j]+xr[j+1]);
-      Vj[j] = xr[j+1]-xr[j];
-      });
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-	rhoj[j] = mj[j]/Vj[j];
-      });
-
-    block_eos.updatePandCFromRhoE();    
-    
-    ++iteration;
-  }
-
-  std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
-	    << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-
-}
diff --git a/src/experimental/MeshLessAcousticSolver.hpp b/src/experimental/MeshLessAcousticSolver.hpp
deleted file mode 100644
index ccdba5838..000000000
--- a/src/experimental/MeshLessAcousticSolver.hpp
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef MESH_LESS_ACOUSTIC_SOLVER_HPP
-#define MESH_LESS_ACOUSTIC_SOLVER_HPP
-
-#include <Kokkos_Core.hpp>
-
-class MeshLessAcousticSolver 
-{
-private:
-
-  inline double e(double rho, double p, double gamma) const;
-  inline double p(double rho, double e, double gamma) const;
-
-  inline double acoustic_dt(const Kokkos::View<const double*>& Vj,
-			    const Kokkos::View<const double*>& cj) const;
-
-  inline void computeExplicitFluxes(const Kokkos::View<const double*>& xr,
-				    const Kokkos::View<const double*>& xj,
-				    const Kokkos::View<const double*>& rhoj,
-				    const Kokkos::View<const double*>& uj,
-				    const Kokkos::View<const double*>& pj,
-				    const Kokkos::View<const double*>& cj,
-				    const Kokkos::View<const double*>& Vj,
-				    Kokkos::View<double*>& ur,
-				    Kokkos::View<double*>& pr) const;
-
-  struct ReduceMin;
-public:
-  MeshLessAcousticSolver(const long int& nj);
-};
-
-#endif // MESH_LESS_ACOUSTIC_SOLVER_HPP
diff --git a/src/experimental/RawKokkosAcousticSolver.cpp b/src/experimental/RawKokkosAcousticSolver.cpp
deleted file mode 100644
index d27549d34..000000000
--- a/src/experimental/RawKokkosAcousticSolver.cpp
+++ /dev/null
@@ -1,248 +0,0 @@
-#include <RawKokkosAcousticSolver.hpp>
-#include <Kokkos_Core.hpp>
-#include <rang.hpp>
-
-namespace RawKokkos
-{
-
-inline double e(double rho, double p, double gamma)
-{
-  return p/(rho*(gamma-1));
-}
-
-inline double p(double rho, double e, double gamma)
-{
-  return (gamma-1)*rho*e;
-}
-
-typedef const double my_double;
-
-struct ReduceMin {
-private:
-  const Kokkos::View<my_double*> x_;
-
-public:
-  typedef Kokkos::View<my_double*>::non_const_value_type value_type;
-
-  ReduceMin(const Kokkos::View<my_double*>& x) : x_ (x) {}
-
-  typedef Kokkos::View<my_double*>::size_type size_type;
-    
-  KOKKOS_INLINE_FUNCTION void
-  operator() (const size_type i, value_type& update) const
-  {
-    if (x_(i) < update) {
-      update = x_(i);
-    }
-  }
-
-  KOKKOS_INLINE_FUNCTION void
-  join (volatile value_type& dst,
-	const volatile value_type& src) const
-  {
-    if (src < dst) {
-      dst = src;
-    }
-  }
-
-  KOKKOS_INLINE_FUNCTION void
-  init (value_type& dst) const
-  { // The identity under max is -Inf.
-    dst= Kokkos::reduction_identity<value_type>::min();
-  }
-};
-    
-
-double acoustic_dt(const Kokkos::View<double*>& Vj,
-		   const Kokkos::View<double*>& cj)
-{
-  const size_t nj = Vj.size();
-  double dt = std::numeric_limits<double>::max();
-
-  Kokkos::View<double*> Vj_cj("Vj_cj", nj);
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      Vj_cj[j] = Vj[j]/cj[j];
-    });
-
-  Kokkos::parallel_reduce(nj, ReduceMin(Vj_cj), dt);
-
-  // Kokkos::parallel_reduce(n, KOKKOS_LAMBDA(const long i, long& lcount) {
-  //   lcount += (i % 2) == 0;
-  // }, count);
-  return dt;
-}
-
-
-void computeExplicitFluxes(const Kokkos::View<double*>& xr,
-			   const Kokkos::View<double*>& xj,
-			   const Kokkos::View<double*>& rhoj,
-			   const Kokkos::View<double*>& uj,
-			   const Kokkos::View<double*>& pj,
-			   const Kokkos::View<double*>& cj,
-			   const Kokkos::View<double*>& Vj,
-			   Kokkos::View<double*>& ur,
-			   Kokkos::View<double*>& pr)
-{
-  // calcul de ur
-  ur[0]=0;
-  const size_t nr = ur.size();
-  const size_t nj = uj.size();
-
-  Kokkos::parallel_for(nj-1, KOKKOS_LAMBDA(const int& j) {
-    const int r = j+1;
-    const int k = r;
-    const double ujr = uj[j];
-    const double ukr = uj[k];
-    const double pjr = pj[j];
-    const double pkr = pj[k];
-
-    ur[r]=(rhoj[j]*cj[j]*ujr + rhoj[k]*cj[k]*ukr + pjr-pkr)/(rhoj[j]*cj[j]+rhoj[k]*cj[k]);
-    });
-  ur[nr-1]=0;
-
-  // calcul de pr
-  pr[0] = pj[0] + rhoj[0]*cj[0]*(ur[0] - uj[0]);
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
-    const int r = j+1;
-
-    const double ujr = uj[j];
-    const double pjr = pj[j];
-
-    pr[r]=pjr+rhoj[j]*cj[j]*(ujr-ur[r]);
-    });
-}
-
-
-void AcousticSolver(const long int& nj)
-{
-  Kokkos::View<double*> xj("xj", nj);
-  Kokkos::View<double*> rhoj("rhoj", nj);
-
-  Kokkos::View<double*> uj("uj", nj);
-
-  Kokkos::View<double*> Ej("Ej", nj);
-  Kokkos::View<double*> ej("ej", nj);
-  Kokkos::View<double*> pj("pj", nj);
-  Kokkos::View<double*> Vj("Vj", nj);
-  Kokkos::View<double*> gammaj("gammaj", nj);
-  Kokkos::View<double*> cj("cj", nj);
-  Kokkos::View<double*> mj("mj", nj);
-  Kokkos::View<double*> inv_mj("inv_mj", nj);
-
-  const int nr=nj+1;
-
-  Kokkos::View<double*>  xr("xr", nr);
-
-  const double delta_x = 1./nj;
-
-  Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){
-      xr[r] = r*delta_x;
-    });
-
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      xj[j] = 0.5*(xr[j]+xr[j+1]);
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      Vj[j] = xr[j+1]-xr[j];
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-    if (xj[j]<0.5) {
-      rhoj[j]=1;
-    } else {
-      rhoj[j]=0.125;
-    }
-  });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-    if (xj[j]<0.5) {
-      pj[j]=1;
-    } else {
-      pj[j]=0.1;
-    }
-  });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      gammaj[j] = 1.4;
-    });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-    ej[j] = e(rhoj[j],pj[j],gammaj[j]);
-  });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-    Ej[j] = ej[j]+0.5*uj[j]*uj[j];
-  });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-    cj[j] = std::sqrt(gammaj[j]*pj[j]/rhoj[j]);
-  });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-    mj[j] = rhoj[j] * Vj[j];
-  });
-
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      inv_mj[j] = 1./mj[j];
-  });
-
-  const double tmax=0.2;
-  double t=0;
-
-  int itermax=std::numeric_limits<int>::max();
-  int iteration=0;
-
-  while((t<tmax) and (iteration<itermax)) {
-    double dt = 0.4*acoustic_dt(Vj, cj);
-    if (t+dt<tmax) {
-      t+=dt;
-    } else {
-      dt=tmax-t;
-      t=tmax;
-    }
-    
-    Kokkos::View<double*> ur("ur", nr);
-    Kokkos::View<double*> pr("pr", nr);
-
-    computeExplicitFluxes(xr, xj,
-			  rhoj, uj, pj, cj, Vj,
-			  ur, pr);
-    
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      int rm=j;
-      int rp=j+1;
-
-      uj[j] += dt*inv_mj[j]*(pr[rm]-pr[rp]);
-      Ej[j] += dt*inv_mj[j]*(pr[rm]*ur[rm]-pr[rp]*ur[rp]);
-      });
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-	ej[j] = Ej[j] - 0.5 * uj[j]*uj[j];
-      });
-
-    Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){
-      xr[r] += dt*ur[r];
-      });
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-      xj[j] = 0.5*(xr[j]+xr[j+1]);
-      Vj[j] = xr[j+1]-xr[j];
-      });
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-	rhoj[j] = mj[j]/Vj[j];
-	pj[j] = p(rhoj[j], ej[j], gammaj[j]);
-	cj[j] = std::sqrt(gammaj[j]*pj[j]/rhoj[j]); // inv_mj*vj
-      });
-    
-    ++iteration;
-  }
-
-  std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
-	    << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-}
-
-}
diff --git a/src/experimental/RawKokkosAcousticSolver.hpp b/src/experimental/RawKokkosAcousticSolver.hpp
deleted file mode 100644
index e9c36ff43..000000000
--- a/src/experimental/RawKokkosAcousticSolver.hpp
+++ /dev/null
@@ -1,9 +0,0 @@
-#ifndef RAW_KOKKOS_ACOUSTIC_SOLVER_HPP
-#define RAW_KOKKOS_ACOUSTIC_SOLVER_HPP
-
-namespace RawKokkos
-{
-  void AcousticSolver(const long int& nj);
-}
-
-#endif // RAW_KOKKOS_ACOUSTIC_SOLVER_HPP
diff --git a/src/main.cpp b/src/main.cpp
index 8e640d841..06aec1142 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -202,11 +202,11 @@ int main(int argc, char *argv[])
           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();
+          CellValue<double>& rhoj = unknowns.rhoj();
+          CellValue<double>& ej = unknowns.ej();
+          CellValue<double>& pj = unknowns.pj();
+          CellValue<double>& gammaj = unknowns.gammaj();
+          CellValue<double>& cj = unknowns.cj();
 
           BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
 
@@ -233,8 +233,8 @@ int main(int argc, char *argv[])
           method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
 
           { // gnuplot output for density
-            const Kokkos::View<const Rd*> xj   = mesh_data.xj();
-            const Kokkos::View<const double*> rhoj = unknowns.rhoj();
+            const CellValue<const Rd>& xj = mesh_data.xj();
+            const CellValue<const double>& rhoj = unknowns.rhoj();
             std::ofstream fout("rho");
             for (size_t j=0; j<mesh.numberOfCells(); ++j) {
               fout << xj[j][0] << ' ' << rhoj[j] << '\n';
@@ -309,11 +309,11 @@ int main(int argc, char *argv[])
           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();
+          CellValue<double>& rhoj = unknowns.rhoj();
+          CellValue<double>& ej = unknowns.ej();
+          CellValue<double>& pj = unknowns.pj();
+          CellValue<double>& gammaj = unknowns.gammaj();
+          CellValue<double>& cj = unknowns.cj();
 
           BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
 
@@ -405,11 +405,11 @@ int main(int argc, char *argv[])
           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();
+          CellValue<double>& rhoj = unknowns.rhoj();
+          CellValue<double>& ej = unknowns.ej();
+          CellValue<double>& pj = unknowns.pj();
+          CellValue<double>& gammaj = unknowns.gammaj();
+          CellValue<double>& cj = unknowns.cj();
 
           BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
 
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index eb5fb87b3..c3a2c01a9 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -9,7 +9,7 @@ void Connectivity<3>::_computeCellFaceAndFaceNodeConnectivities()
   const auto& cell_to_node_matrix
       = this->getMatrix(ItemType::cell, ItemType::node);
 
-  Kokkos::View<unsigned short*> cell_nb_faces("cell_nb_faces", this->numberOfCells());
+  CellValue<unsigned short> cell_nb_faces(*this);
   std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
   for (unsigned int j=0; j<this->numberOfCells(); ++j) {
     const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
@@ -227,14 +227,14 @@ Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
   Assert(this->numberOfCells()>0);
 
   {
-    Kokkos::View<CellType*> cell_type("cell_type", this->numberOfCells());
+    CellValue<CellType> cell_type(*this);
     Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& j){
         cell_type[j] = cell_type_vector[j];
       });
     m_cell_type = cell_type;
   }
   {
-    Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", this->numberOfCells());
+    CellValue<double> inv_cell_nb_nodes(*this);
     Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& j){
         const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
         inv_cell_nb_nodes[j] = 1./cell_nodes.length;
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index b72a38abd..964daeb26 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -5,6 +5,7 @@
 #include <TinyVector.hpp>
 
 #include <Kokkos_Core.hpp>
+#include <ItemValue.hpp>
 
 #include <IConnectivity.hpp>
 
@@ -234,7 +235,7 @@ class Connectivity final
 
  private:
   ConnectivityMatrix m_item_to_item_matrix[Dimension+1][Dimension+1];
-  Kokkos::View<const CellType*> m_cell_type;
+  CellValue<const CellType> m_cell_type;
 
   FaceValuePerCell<const bool> m_cell_face_is_reversed;
 
@@ -259,7 +260,7 @@ class Connectivity final
   std::vector<RefFaceList> m_ref_face_list;
   std::vector<RefNodeList> m_ref_node_list;
 
-  Kokkos::View<double*> m_inv_cell_nb_nodes;
+  CellValue<const double> m_inv_cell_nb_nodes;
 
   using Face = ConnectivityFace<Dimension>;
 
@@ -470,7 +471,7 @@ class Connectivity final
     return cell_to_node_matrix.numRows();
   }
 
-  const Kokkos::View<const double*> invCellNbNodes() const
+  const CellValue<const double>& invCellNbNodes() const
   {
 #warning add calculation on demand when variables will be defined
     return m_inv_cell_nb_nodes;
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index f16ab2fa0..d483ee924 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -846,7 +846,7 @@ GmshReader::__proceedData()
     typedef Mesh<Connectivity3D> MeshType;
     typedef TinyVector<3, double> Rd;
 
-    Kokkos::View<Rd*> xr("xr", __vertices.extent(0));
+    NodeValue<Rd> xr(connectivity);
     for (size_t i=0; i<__vertices.extent(0); ++i) {
       xr[i][0] = __vertices[i][0];
       xr[i][1] = __vertices[i][1];
@@ -918,7 +918,7 @@ GmshReader::__proceedData()
     typedef Mesh<Connectivity2D> MeshType;
     typedef TinyVector<2, double> Rd;
 
-    Kokkos::View<Rd*> xr("xr", __vertices.extent(0));
+    NodeValue<Rd> xr(connectivity);
     for (size_t i=0; i<__vertices.extent(0); ++i) {
       xr[i][0] = __vertices[i][0];
       xr[i][1] = __vertices[i][1];
@@ -962,7 +962,7 @@ GmshReader::__proceedData()
     typedef Mesh<Connectivity1D> MeshType;
     typedef TinyVector<1, double> Rd;
 
-    Kokkos::View<Rd*> xr("xr", __vertices.extent(0));
+    NodeValue<Rd> xr(connectivity);
     for (size_t i=0; i<__vertices.extent(0); ++i) {
       xr[i][0] = __vertices[i][0];
     }
diff --git a/src/mesh/ItemValue.hpp b/src/mesh/ItemValue.hpp
index a00cb3bfd..5534df489 100644
--- a/src/mesh/ItemValue.hpp
+++ b/src/mesh/ItemValue.hpp
@@ -1,6 +1,8 @@
 #ifndef ITEM_VALUE_HPP
 #define ITEM_VALUE_HPP
 
+#include <Kokkos_Core.hpp>
+
 #include <ItemType.hpp>
 #include <PastisAssert.hpp>
 
@@ -12,6 +14,7 @@ class ItemValue
 {
  public:
   static const ItemType item_t{item_type};
+  typedef DataType data_type;
 
  private:
   bool m_is_built{false};
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 3aee2c287..383e92294 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -1,7 +1,7 @@
 #ifndef MESH_HPP
 #define MESH_HPP
 
-#include <Kokkos_Core.hpp>
+#include <ItemValue.hpp>
 #include <TinyVector.hpp>
 
 #include <memory>
@@ -22,15 +22,17 @@ public:
 
 private:
   const std::shared_ptr<Connectivity> m_connectivity;
-  Kokkos::View<Rd*> m_xr;
+  NodeValue<const Rd> m_xr;
+  NodeValue<Rd> m_mutable_xr;
 
 public:
-
+  KOKKOS_INLINE_FUNCTION
   const size_t meshDimension() const
   {
     return dimension;
   }
 
+  KOKKOS_INLINE_FUNCTION
   const Connectivity& connectivity() const
   {
     return *m_connectivity;
@@ -54,28 +56,25 @@ public:
     return m_connectivity->numberOfCells();
   }
 
-  Kokkos::View<Rd*> xr() const
+  KOKKOS_INLINE_FUNCTION
+  const NodeValue<const Rd>& xr() const
   {
     return m_xr;
   }
 
+  [[deprecated("should rework this class: quite ugly")]]
   KOKKOS_INLINE_FUNCTION
-  Mesh(const std::shared_ptr<Connectivity>& connectivity)
-    : m_connectivity(connectivity),
-      m_xr("xr", connectivity->numberOfNodes())
+  NodeValue<Rd> mutableXr() const
   {
-    static_assert (Connectivity::dimension ==1,"no automatic calculation of vertices in dim>1");
-    const double delta_x = 1./connectivity->numberOfCells();
-    Kokkos::parallel_for(connectivity->numberOfNodes(), KOKKOS_LAMBDA(const int& r){
-	m_xr[r][0] = r*delta_x;
-      });
+    return m_mutable_xr;
   }
 
   KOKKOS_INLINE_FUNCTION
   Mesh(const std::shared_ptr<Connectivity>& connectivity,
-       Kokkos::View<Rd*>& xr)
-    : m_connectivity(connectivity),
-      m_xr(xr)
+       NodeValue<Rd>& xr)
+      : m_connectivity{connectivity},
+        m_xr{xr},
+        m_mutable_xr{xr}
   {
     ;
   }
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 254ff42b1..aa24493e8 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -24,50 +24,54 @@ class MeshData
 
  private:
   const MeshType& m_mesh;
-  NodeValuePerCell<Rd> m_Cjr;
-  NodeValuePerCell<double> m_ljr;
-  NodeValuePerCell<Rd> m_njr;
-  Kokkos::View<Rd*>  m_xj;
+  NodeValuePerCell<const Rd> m_Cjr;
+  NodeValuePerCell<const double> m_ljr;
+  NodeValuePerCell<const Rd> m_njr;
+  CellValue<const Rd> m_xj;
   CellValue<const double> m_Vj;
 
   KOKKOS_INLINE_FUNCTION
   void _updateCenter()
   { // Computes vertices isobarycenter
     if constexpr (dimension == 1) {
-      const Kokkos::View<const Rd*> xr = m_mesh.xr();
+      const NodeValue<const Rd>& xr = m_mesh.xr();
 
       const auto& cell_to_node_matrix
           = m_mesh.connectivity().getMatrix(ItemType::cell,
                                             ItemType::node);
 
+      CellValue<Rd> xj(m_mesh.connectivity());
       Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
           const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
-          m_xj[j] = 0.5*(xr[cell_nodes(0)]+xr[cell_nodes(1)]);
+          xj[j] = 0.5*(xr[cell_nodes(0)]+xr[cell_nodes(1)]);
         });
+      m_xj = xj;
     } else {
-      const Kokkos::View<const Rd*> xr = m_mesh.xr();
-      const Kokkos::View<const double*>& inv_cell_nb_nodes
+      const NodeValue<const Rd>& xr = m_mesh.xr();
+
+      const CellValue<const double>& inv_cell_nb_nodes
           = m_mesh.connectivity().invCellNbNodes();
 
       const auto& cell_to_node_matrix
           = m_mesh.connectivity().getMatrix(ItemType::cell,
                                             ItemType::node);
-
+      CellValue<Rd> xj(m_mesh.connectivity());
       Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
           Rd X = zero;
           const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
           for (size_t R=0; R<cell_nodes.length; ++R) {
             X += xr[cell_nodes(R)];
           }
-          m_xj[j] = inv_cell_nb_nodes[j]*X;
+          xj[j] = inv_cell_nb_nodes[j]*X;
         });
+      m_xj = xj;
     }
   }
 
   KOKKOS_INLINE_FUNCTION
   void _updateVolume()
   {
-    const Kokkos::View<const Rd*> xr = m_mesh.xr();
+    const NodeValue<const Rd>& xr = m_mesh.xr();
     const auto& cell_to_node_matrix
         = m_mesh.connectivity().getMatrix(ItemType::cell,
                                           ItemType::node);
@@ -90,33 +94,42 @@ class MeshData
       // Cjr/njr/ljr are constant overtime
     }
     else if constexpr (dimension == 2) {
-      const Kokkos::View<const Rd*> xr = m_mesh.xr();
+      const NodeValue<const Rd>& xr = m_mesh.xr();
       const auto& cell_to_node_matrix
           = m_mesh.connectivity().getMatrix(ItemType::cell,
                                             ItemType::node);
 
-      Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-          const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
-          for (size_t R=0; R<cell_nodes.length; ++R) {
-            int Rp1 = (R+1)%cell_nodes.length;
-            int Rm1 = (R+cell_nodes.length-1)%cell_nodes.length;
-            Rd half_xrp_xrm = 0.5*(xr(cell_nodes(Rp1))-xr(cell_nodes(Rm1)));
-            m_Cjr(j,R) = Rd{-half_xrp_xrm[1], half_xrp_xrm[0]};
-          }
-        });
-
-      const NodeValuePerCell<Rd>& Cjr = m_Cjr;
-      Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& j){
-          m_ljr[j] = l2Norm(Cjr[j]);
-        });
-
-      const NodeValuePerCell<double>& ljr = m_ljr;
-      Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& j){
-          m_njr[j] = (1./ljr[j])*Cjr[j];
-        });
-
+      {
+        NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
+        Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+            const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+            for (size_t R=0; R<cell_nodes.length; ++R) {
+              int Rp1 = (R+1)%cell_nodes.length;
+              int Rm1 = (R+cell_nodes.length-1)%cell_nodes.length;
+              Rd half_xrp_xrm = 0.5*(xr(cell_nodes(Rp1))-xr(cell_nodes(Rm1)));
+              Cjr(j,R) = Rd{-half_xrp_xrm[1], half_xrp_xrm[0]};
+            }
+          });
+        m_Cjr = Cjr;
+      }
+
+      {
+        NodeValuePerCell<double> ljr(m_mesh.connectivity());
+        Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& j){
+            ljr[j] = l2Norm(m_Cjr[j]);
+          });
+        m_ljr = ljr;
+      }
+
+      {
+        NodeValuePerCell<Rd> njr(m_mesh.connectivity());
+        Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& j){
+            njr[j] = (1./m_ljr[j])*m_Cjr[j];
+          });
+        m_njr = njr;
+      }
     } else if (dimension ==3) {
-      const Kokkos::View<const Rd*> xr = m_mesh.xr();
+      const NodeValue<const Rd>& xr = m_mesh.xr();
 
       NodeValuePerFace<Rd> Nlr(m_mesh.connectivity());
       const auto& face_to_node_matrix
@@ -143,10 +156,6 @@ class MeshData
           }
         });
 
-      Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
-          m_Cjr[jr] = zero;
-        });
-
       const auto& cell_to_node_matrix
           = m_mesh.connectivity().getMatrix(ItemType::cell,
                                             ItemType::node);
@@ -156,51 +165,67 @@ class MeshData
                                             ItemType::face);
 
       const auto& cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
-      Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
-          const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
 
-          const auto& cell_faces = cell_to_face_matrix.rowConst(j);
-          const auto& face_is_reversed = cell_face_is_reversed.itemValues(j);
+      {
+        NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
+        Kokkos::parallel_for(Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
+            Cjr[jr] = zero;
+          });
+
+        Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+            const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
 
-          for (size_t L=0; L<cell_faces.length; ++L) {
-            const size_t l = cell_faces(L);
-            const auto& face_nodes = face_to_node_matrix.rowConst(l);
+            const auto& cell_faces = cell_to_face_matrix.rowConst(j);
+            const auto& face_is_reversed = cell_face_is_reversed.itemValues(j);
+
+            for (size_t L=0; L<cell_faces.length; ++L) {
+              const size_t l = cell_faces(L);
+              const auto& face_nodes = face_to_node_matrix.rowConst(l);
 
 #warning should this lambda be replaced by a precomputed correspondance?
-            std::function local_node_number_in_cell
-                = [&](const size_t& node_number) {
-                    for (size_t i_node=0; i_node<cell_nodes.length; ++i_node) {
-                      if (node_number == cell_nodes(i_node)) {
-                        return i_node;
-                        break;
+              std::function local_node_number_in_cell
+                  = [&](const size_t& node_number) {
+                      for (size_t i_node=0; i_node<cell_nodes.length; ++i_node) {
+                        if (node_number == cell_nodes(i_node)) {
+                          return i_node;
+                          break;
+                        }
                       }
-                    }
-                    return std::numeric_limits<size_t>::max();
-                  };
-
-            if (face_is_reversed[L]) {
-              for (size_t rl = 0; rl<face_nodes.length; ++rl) {
-                const size_t R = local_node_number_in_cell(face_nodes(rl));
-                m_Cjr(j, R) -= Nlr(l,rl);
-              }
-            } else {
-              for (size_t rl = 0; rl<face_nodes.length; ++rl) {
-                const size_t R = local_node_number_in_cell(face_nodes(rl));
-                m_Cjr(j, R) += Nlr(l,rl);
+                      return std::numeric_limits<size_t>::max();
+                    };
+
+              if (face_is_reversed[L]) {
+                for (size_t rl = 0; rl<face_nodes.length; ++rl) {
+                  const size_t R = local_node_number_in_cell(face_nodes(rl));
+                  Cjr(j, R) -= Nlr(l,rl);
+                }
+              } else {
+                for (size_t rl = 0; rl<face_nodes.length; ++rl) {
+                  const size_t R = local_node_number_in_cell(face_nodes(rl));
+                  Cjr(j, R) += Nlr(l,rl);
+                }
               }
             }
-          }
-        });
-
-      const NodeValuePerCell<Rd>& Cjr = m_Cjr;
-      Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
-          m_ljr[jr] = l2Norm(Cjr[jr]);
-        });
-
-      const NodeValuePerCell<double>& ljr = m_ljr;
-      Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
-          m_njr[jr] = (1./ljr[jr])*Cjr[jr];
-        });
+          });
+
+        m_Cjr = Cjr;
+      }
+
+      {
+        NodeValuePerCell<double> ljr(m_mesh.connectivity());
+        Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
+            ljr[jr] = l2Norm(m_Cjr[jr]);
+          });
+        m_ljr = ljr;
+      }
+
+      {
+        NodeValuePerCell<Rd> njr(m_mesh.connectivity());
+        Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
+            njr[jr] = (1./m_ljr[jr])*m_Cjr[jr];
+          });
+        m_njr = njr;
+      }
     }
     static_assert((dimension<=3), "only 1d, 2d and 3d are implemented");
   }
@@ -211,22 +236,22 @@ class MeshData
     return m_mesh;
   }
 
-  const NodeValuePerCell<Rd>& Cjr() const
+  const NodeValuePerCell<const Rd>& Cjr() const
   {
     return m_Cjr;
   }
 
-  const NodeValuePerCell<double>& ljr() const
+  const NodeValuePerCell<const double>& ljr() const
   {
     return m_ljr;
   }
 
-  const NodeValuePerCell<Rd>& njr() const
+  const NodeValuePerCell<const Rd>& njr() const
   {
     return m_njr;
   }
 
-  const Kokkos::View<const Rd*> xj() const
+  const CellValue<const Rd>& xj() const
   {
     return m_xj;
   }
@@ -244,25 +269,27 @@ class MeshData
   }
 
   MeshData(const MeshType& mesh)
-      : m_mesh(mesh),
-        m_Cjr(mesh.connectivity()),
-        m_ljr(mesh.connectivity()),
-        m_njr(mesh.connectivity()),
-        m_xj("xj", mesh.numberOfCells()),
-        m_Vj(mesh.connectivity())
+      : m_mesh(mesh)
   {
     if constexpr (dimension==1) {
       // in 1d Cjr are computed once for all
-      Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
-          m_Cjr(j,0)=-1;
-          m_Cjr(j,1)= 1;
-        });
+      {
+        NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
+        Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+            Cjr(j,0)=-1;
+            Cjr(j,1)= 1;
+          });
+        m_Cjr = Cjr;
+      }
       // in 1d njr=Cjr (here performs a shallow copy)
-      m_njr=m_Cjr;
-      Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
-          m_ljr(j,0)= 1;
-          m_ljr(j,1)= 1;
+      m_njr = m_Cjr;
+      {
+        NodeValuePerCell<double> ljr(m_mesh.connectivity());
+        Kokkos::parallel_for(ljr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
+            ljr[jr] = 1;
         });
+        m_ljr = ljr;
+      }
     }
     this->updateAllData();
   }
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 87c537439..91a0943de 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -2,6 +2,8 @@
 #define MESH_NODE_BOUNDARY_HPP
 
 #include <Kokkos_Core.hpp>
+#include <ItemValue.hpp>
+
 #include <Kokkos_Vector.hpp>
 #include <TinyVector.hpp>
 
@@ -149,7 +151,7 @@ _checkBoundaryIsFlat(const TinyVector<2,double>& normal,
   const R2 origin = 0.5*(xmin+xmax);
   const double length = l2Norm(xmax-xmin);
 
-  const Kokkos::View<const R2*> xr = mesh.xr();
+  const NodeValue<const R2>& xr = mesh.xr();
 
   Kokkos::parallel_for(m_node_list.extent(0), KOKKOS_LAMBDA(const size_t& r) {
       const R2& x = xr[m_node_list[r]];
@@ -186,7 +188,7 @@ _getNormal(const MeshType& mesh)
   static_assert(MeshType::dimension == 2);
   typedef TinyVector<2,double> R2;
 
-  const Kokkos::View<const R2*> xr = mesh.xr();
+  const NodeValue<const R2>& xr = mesh.xr();
 
   R2 xmin(std::numeric_limits<double>::max(),
           std::numeric_limits<double>::max());
@@ -239,7 +241,7 @@ _getNormal(const MeshType& mesh)
   R3 ymax = xmax;
   R3 zmax = xmax;
 
-  const Kokkos::View<const R3*> xr = mesh.xr();
+  const NodeValue<const R3>& xr = mesh.xr();
 
   for (size_t r=0; r<m_node_list.extent(0); ++r) {
     const R3& x = xr[m_node_list[r]];
@@ -314,7 +316,7 @@ _getOutgoingNormal(const MeshType& mesh)
 
   const R normal = this->_getNormal(mesh);
 
-  const Kokkos::View<const R*>& xr = mesh.xr();
+  const NodeValue<const R>& xr = mesh.xr();
   const auto& cell_to_node_matrix
       = mesh.connectivity().getMatrix(ItemType::cell,
                                       ItemType::node);
@@ -351,7 +353,7 @@ _getOutgoingNormal(const MeshType& mesh)
 
   const R2 normal = this->_getNormal(mesh);
 
-  const Kokkos::View<const R2*>& xr = mesh.xr();
+  const NodeValue<const R2>& xr = mesh.xr();
   const auto& cell_to_node_matrix
       = mesh.connectivity().getMatrix(ItemType::cell,
                                       ItemType::node);
@@ -388,7 +390,7 @@ _getOutgoingNormal(const MeshType& mesh)
 
   const R3 normal = this->_getNormal(mesh);
 
-  const Kokkos::View<const R3*>& xr = mesh.xr();
+  const NodeValue<const R3>& xr = mesh.xr();
   const auto& cell_to_node_matrix
       = mesh.connectivity().getMatrix(ItemType::cell,
                                       ItemType::node);
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index 52d91666b..be60e2634 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -182,7 +182,7 @@ class SubItemValuePerItem<DataType,
         = connectivity.getMatrix(item_type, sub_item_type);
 
     m_host_row_map = connectivity_matrix.rowsMap();
-    m_values = Kokkos::View<DataType*>("values", connectivity_matrix.numEntries());
+    m_values = Kokkos::View<std::remove_const_t<DataType>*>("values", connectivity_matrix.numEntries());
   }
 
   ~SubItemValuePerItem() = default;
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index a33fef0bf..24865f43d 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -8,6 +8,8 @@
 #include <TinyVector.hpp>
 #include <IConnectivity.hpp>
 
+#include <ItemValue.hpp>
+
 class VTKWriter
 {
  private:
@@ -39,16 +41,16 @@ class VTKWriter
 
     fout << "<Points>\n";
     fout << "<DataArray Name=\"Positions\" NumberOfComponents=\"3\" type=\"Float64\" format=\"ascii\">\n";
-    if constexpr(MeshType::dimension ==1) {
-      const Kokkos::View<const TinyVector<1>*> xr = mesh.xr();
+    using Rd = TinyVector<MeshType::dimension>;
+    const NodeValue<const Rd>& xr = mesh.xr();
+    if constexpr(MeshType::dimension == 1) {
       for (unsigned int r=0; r<mesh.numberOfNodes(); ++r) {
         for (unsigned short i=0; i<1; ++i) {
           fout << xr[r][i] << ' ';
         }
         fout << "0 0 "; // VTK requires 3 components
       }
-    } else if constexpr (MeshType::dimension ==2) {
-      const Kokkos::View<const TinyVector<2>*> xr = mesh.xr();
+    } else if constexpr (MeshType::dimension == 2) {
       for (unsigned int r=0; r<mesh.numberOfNodes(); ++r) {
         for (unsigned short i=0; i<2; ++i) {
           fout << xr[r][i] << ' ';
@@ -56,7 +58,6 @@ class VTKWriter
         fout << "0 "; // VTK requires 3 components
       }
     } else {
-      const Kokkos::View<const TinyVector<3>*> xr = mesh.xr();
       for (unsigned int r=0; r<mesh.numberOfNodes(); ++r) {
         for (unsigned short i=0; i<3; ++i) {
           fout << xr[r][i] << ' ';
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 6f83bdcaa..3bc6269ee 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -36,14 +36,15 @@ class AcousticSolver
   struct ReduceMin
   {
    private:
-    const Kokkos::View<const double*> x_;
+    const CellValue<const double> x_;
 
    public:
-    typedef Kokkos::View<const double*>::non_const_value_type value_type;
+    typedef std::remove_const_t<CellValue<const double>::data_type> value_type;
 
-    ReduceMin(const Kokkos::View<const double*>& x) : x_ (x) {}
+    ReduceMin(const CellValue<const double>& x) : x_ (x) {}
 
-    typedef Kokkos::View<const double*>::size_type size_type;
+#warning was: "typedef Kokkos::View<const double*>::size_type size_type;""
+    typedef size_t size_type;
 
     KOKKOS_INLINE_FUNCTION
     void operator() (const size_type i, value_type& update) const
@@ -70,9 +71,9 @@ class AcousticSolver
   };
 
   KOKKOS_INLINE_FUNCTION
-  const Kokkos::View<const double*>
-  computeRhoCj(const Kokkos::View<const double*>& rhoj,
-               const Kokkos::View<const double*>& cj)
+  const CellValue<const double>
+  computeRhoCj(const CellValue<const double>& rhoj,
+               const CellValue<const double>& cj)
   {
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
         m_rhocj[j] = rhoj[j]*cj[j];
@@ -81,10 +82,10 @@ class AcousticSolver
   }
 
   KOKKOS_INLINE_FUNCTION
-  void computeAjr(const Kokkos::View<const double*>& rhocj,
-                  const NodeValuePerCell<Rd>& Cjr,
-                  const NodeValuePerCell<double>& ljr,
-                  const NodeValuePerCell<Rd>& njr)
+  void computeAjr(const CellValue<const double>& rhocj,
+                  const NodeValuePerCell<const Rd>& Cjr,
+                  const NodeValuePerCell<const double>& ljr,
+                  const NodeValuePerCell<const Rd>& njr)
   {
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
         const size_t& nb_nodes =m_Ajr.numberOfSubValues(j);
@@ -96,8 +97,8 @@ class AcousticSolver
   }
 
   KOKKOS_INLINE_FUNCTION
-  const Kokkos::View<const Rdd*>
-  computeAr(const NodeValuePerCell<Rdd>& Ajr) {
+  const NodeValue<const Rdd>
+  computeAr(const NodeValuePerCell<const Rdd>& Ajr) {
     const auto& node_to_cell_matrix
         = m_connectivity.getMatrix(ItemType::node,
                                    ItemType::cell);
@@ -122,11 +123,11 @@ class AcousticSolver
   }
 
   KOKKOS_INLINE_FUNCTION
-  const Kokkos::View<const Rd*>
+  const NodeValue<const Rd>
   computeBr(const NodeValuePerCell<Rdd>& Ajr,
-            const NodeValuePerCell<Rd>& Cjr,
-            const Kokkos::View<const Rd*>& uj,
-            const Kokkos::View<const double*>& pj) {
+            const NodeValuePerCell<const Rd>& Cjr,
+            const CellValue<const Rd>& uj,
+            const CellValue<const double>& pj) {
     const auto& node_to_cell_matrix
         = m_connectivity.getMatrix(ItemType::node,
                                    ItemType::cell);
@@ -191,12 +192,12 @@ class AcousticSolver
     }
   }
 
-  Kokkos::View<Rd*>
-  computeUr(const Kokkos::View<const Rdd*>& Ar,
-            const Kokkos::View<const Rd*>& br)
+  NodeValue<Rd>
+  computeUr(const NodeValue<const Rdd>& Ar,
+            const NodeValue<const Rd>& br)
   {
     inverse(Ar, m_inv_Ar);
-    const Kokkos::View<const Rdd*> invAr = m_inv_Ar;
+    const NodeValue<const Rdd> invAr = m_inv_Ar;
     Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
         m_ur[r]=invAr(r)*br(r);
       });
@@ -206,10 +207,10 @@ class AcousticSolver
 
   void
   computeFjr(const NodeValuePerCell<Rdd>& Ajr,
-             const Kokkos::View<const Rd*>& ur,
-             const NodeValuePerCell<Rd>& Cjr,
-             const Kokkos::View<const Rd*>& uj,
-             const Kokkos::View<const double*>& pj)
+             const NodeValue<const Rd>& ur,
+             const NodeValuePerCell<const Rd>& Cjr,
+             const CellValue<const Rd>& uj,
+             const CellValue<const double>& pj)
   {
     const auto& cell_to_node_matrix
         = m_mesh.connectivity().getMatrix(ItemType::cell,
@@ -224,8 +225,8 @@ class AcousticSolver
       });
   }
 
-  void inverse(const Kokkos::View<const Rdd*>& A,
-               Kokkos::View<Rdd*>& inv_A) const
+  void inverse(const NodeValue<const Rdd>& A,
+               NodeValue<Rdd>& inv_A) const
   {
     Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
         inv_A(r) = ::inverse(A(r));
@@ -233,37 +234,38 @@ class AcousticSolver
   }
 
   KOKKOS_INLINE_FUNCTION
-  void computeExplicitFluxes(const Kokkos::View<const Rd*>& xr,
-                             const Kokkos::View<const Rd*>& xj,
-                             const Kokkos::View<const double*>& rhoj,
-                             const Kokkos::View<const Rd*>& uj,
-                             const Kokkos::View<const double*>& pj,
-                             const Kokkos::View<const double*>& cj,
+  void computeExplicitFluxes(const NodeValue<const Rd>& xr,
+                             const CellValue<const Rd>& xj,
+                             const CellValue<const double>& rhoj,
+                             const CellValue<const Rd>& uj,
+                             const CellValue<const double>& pj,
+                             const CellValue<const double>& cj,
                              const CellValue<const double>& Vj,
-                             const NodeValuePerCell<Rd>& Cjr,
-                             const NodeValuePerCell<double>& ljr,
-                             const NodeValuePerCell<Rd>& njr) {
-    const Kokkos::View<const double*> rhocj  = computeRhoCj(rhoj, cj);
+                             const NodeValuePerCell<const Rd>& Cjr,
+                             const NodeValuePerCell<const double>& ljr,
+                             const NodeValuePerCell<const Rd>& njr) {
+    const CellValue<const double> rhocj = computeRhoCj(rhoj, cj);
     computeAjr(rhocj, Cjr, ljr, njr);
 
-    const Kokkos::View<const Rdd*> Ar = computeAr(m_Ajr);
-    const Kokkos::View<const Rd*> br = computeBr(m_Ajr, Cjr, uj, pj);
+    NodeValuePerCell<const Rdd> Ajr = m_Ajr;;
+    const NodeValue<const Rdd> Ar = computeAr(Ajr);
+    const NodeValue<const Rd> br = computeBr(m_Ajr, Cjr, uj, pj);
 
     this->applyBoundaryConditions();
 
-    Kokkos::View<Rd*> ur = m_ur;
-    ur  = computeUr(Ar, br);
+    NodeValue<Rd>& ur = m_ur;
+    ur = computeUr(Ar, br);
     computeFjr(m_Ajr, ur, Cjr, uj, pj);
   }
 
-  Kokkos::View<Rd*> m_br;
+  NodeValue<Rd> m_br;
   NodeValuePerCell<Rdd> m_Ajr;
-  Kokkos::View<Rdd*> m_Ar;
-  Kokkos::View<Rdd*> m_inv_Ar;
+  NodeValue<Rdd> m_Ar;
+  NodeValue<Rdd> m_inv_Ar;
   NodeValuePerCell<Rd> m_Fjr;
-  Kokkos::View<Rd*> m_ur;
-  Kokkos::View<double*> m_rhocj;
-  Kokkos::View<double*> m_Vj_over_cj;
+  NodeValue<Rd> m_ur;
+  CellValue<double> m_rhocj;
+  CellValue<double> m_Vj_over_cj;
 
  public:
   AcousticSolver(MeshData& mesh_data,
@@ -273,23 +275,23 @@ class AcousticSolver
         m_mesh(mesh_data.mesh()),
         m_connectivity(m_mesh.connectivity()),
         m_boundary_condition_list(bc_list),
-        m_br("br", m_mesh.numberOfNodes()),
+        m_br(m_connectivity),
         m_Ajr(m_connectivity),
-        m_Ar("Ar", m_mesh.numberOfNodes()),
-        m_inv_Ar("inv_Ar", m_mesh.numberOfNodes()),
+        m_Ar(m_connectivity),
+        m_inv_Ar(m_connectivity),
         m_Fjr(m_connectivity),
-        m_ur("ur", m_mesh.numberOfNodes()),
-        m_rhocj("rho_c", m_mesh.numberOfCells()),
-        m_Vj_over_cj("Vj_over_cj", m_mesh.numberOfCells())
+        m_ur(m_connectivity),
+        m_rhocj(m_connectivity),
+        m_Vj_over_cj(m_connectivity)
   {
     ;
   }
 
   KOKKOS_INLINE_FUNCTION
   double acoustic_dt(const CellValue<const double>& Vj,
-                     const Kokkos::View<const double*>& cj) const
+                     const CellValue<const double>& cj) const
   {
-    const NodeValuePerCell<double>& ljr = m_mesh_data.ljr();
+    const NodeValuePerCell<const double>& ljr = m_mesh_data.ljr();
     const auto& cell_to_node_matrix
         = m_mesh.connectivity().getMatrix(ItemType::cell,
                                           ItemType::node);
@@ -314,31 +316,31 @@ class AcousticSolver
   void computeNextStep(const double& t, const double& dt,
                        UnknownsType& unknowns)
   {
-    Kokkos::View<double*> rhoj = unknowns.rhoj();
-    Kokkos::View<Rd*> uj = unknowns.uj();
-    Kokkos::View<double*> Ej = unknowns.Ej();
+    CellValue<double>& rhoj = unknowns.rhoj();
+    CellValue<Rd>& uj = unknowns.uj();
+    CellValue<double>& Ej = unknowns.Ej();
 
-    Kokkos::View<double*> ej = unknowns.ej();
-    Kokkos::View<double*> pj = unknowns.pj();
-    Kokkos::View<double*> gammaj = unknowns.gammaj();
-    Kokkos::View<double*> cj = unknowns.cj();
+    CellValue<double>& ej = unknowns.ej();
+    CellValue<double>& pj = unknowns.pj();
+    CellValue<double>& gammaj = unknowns.gammaj();
+    CellValue<double>& cj = unknowns.cj();
 
-    const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
+    const CellValue<const Rd>& xj = m_mesh_data.xj();
     const CellValue<const double>& Vj = m_mesh_data.Vj();
-    const NodeValuePerCell<Rd>& Cjr = m_mesh_data.Cjr();
-    const NodeValuePerCell<double>& ljr = m_mesh_data.ljr();
-    const NodeValuePerCell<Rd>& njr = m_mesh_data.njr();
-    Kokkos::View<Rd*> xr = m_mesh.xr();
+    const NodeValuePerCell<const Rd>& Cjr = m_mesh_data.Cjr();
+    const NodeValuePerCell<const double>& ljr = m_mesh_data.ljr();
+    const NodeValuePerCell<const Rd>& njr = m_mesh_data.njr();
+    const NodeValue<const Rd>& xr = m_mesh.xr();
 
     computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr, ljr, njr);
 
     const NodeValuePerCell<Rd>& Fjr = m_Fjr;
-    const Kokkos::View<const Rd*> ur = m_ur;
+    const NodeValue<const Rd> ur = m_ur;
     const auto& cell_to_node_matrix
         = m_mesh.connectivity().getMatrix(ItemType::cell,
                                           ItemType::node);
 
-    const Kokkos::View<const double*> inv_mj = unknowns.invMj();
+    const CellValue<const double> inv_mj = unknowns.invMj();
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
         const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
 
@@ -357,13 +359,13 @@ class AcousticSolver
         ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
       });
 
+    NodeValue<Rd> mutable_xr = m_mesh.mutableXr();
     Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r){
-        xr[r] += dt*ur[r];
+        mutable_xr[r] += dt*ur[r];
       });
-
     m_mesh_data.updateAllData();
 
-    const Kokkos::View<const double*> mj = unknowns.mj();
+    const CellValue<const double> mj = unknowns.mj();
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
         rhoj[j] = mj[j]/Vj[j];
       });
diff --git a/src/scheme/BlockPerfectGas.hpp b/src/scheme/BlockPerfectGas.hpp
index c018399e5..97a0256ef 100644
--- a/src/scheme/BlockPerfectGas.hpp
+++ b/src/scheme/BlockPerfectGas.hpp
@@ -1,58 +1,60 @@
 #ifndef BLOCK_PERFECTGAS_HPP
 #define BLOCK_PERFECTGAS_HPP
 
+#include <ItemValue.hpp>
+
 struct BlockPerfectGas
 {
 private:
-  Kokkos::View<double*> m_rhoj;
-  Kokkos::View<double*> m_ej;
-  Kokkos::View<double*> m_pj;
-  Kokkos::View<double*> m_gammaj;
-  Kokkos::View<double*> m_cj;
+  CellValue<double>& m_rhoj;
+  CellValue<double>& m_ej;
+  CellValue<double>& m_pj;
+  CellValue<double>& m_gammaj;
+  CellValue<double>& m_cj;
 
 public:
-  BlockPerfectGas(Kokkos::View<double*> rhoj,
-		 Kokkos::View<double*> ej,
-		 Kokkos::View<double*> pj,
-		 Kokkos::View<double*> gammaj,
-		 Kokkos::View<double*> cj)
-    : m_rhoj(rhoj),
-      m_ej(ej),
-      m_pj(pj),
-      m_gammaj(gammaj),
-      m_cj(cj)
+  BlockPerfectGas(CellValue<double>& rhoj,
+                  CellValue<double>& ej,
+                  CellValue<double>& pj,
+                  CellValue<double>& gammaj,
+                  CellValue<double>& cj)
+      : m_rhoj(rhoj),
+        m_ej(ej),
+        m_pj(pj),
+        m_gammaj(gammaj),
+        m_cj(cj)
   {
     ;
   }
 
   void updatePandCFromRhoE()
   {
-    const int nj = m_ej.size();
-    const Kokkos::View<const double*> rho = m_rhoj;
-    const Kokkos::View<const double*> e = m_ej;
-    const Kokkos::View<const double*> gamma = m_gammaj;
-
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-	const double gamma_minus_one = gamma[j]-1;
-	m_pj[j] = gamma_minus_one*rho[j]*e[j];
-	m_cj[j] = std::sqrt(gamma[j]*gamma_minus_one*e[j]);
+    const size_t nj = m_ej.numberOfValues();
+    const CellValue<const double> rho = m_rhoj;
+    const CellValue<const double> e = m_ej;
+    const CellValue<const double> gamma = m_gammaj;
+
+    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const size_t& j){
+        const double gamma_minus_one = gamma[j]-1;
+        m_pj[j] = gamma_minus_one*rho[j]*e[j];
+        m_cj[j] = std::sqrt(gamma[j]*gamma_minus_one*e[j]);
       });
   }
 
   void updateEandCFromRhoP()
   {
-    const int nj = m_ej.size();
-    const Kokkos::View<const double*> rho = m_rhoj;
-    const Kokkos::View<const double*> p = m_pj;
-    const Kokkos::View<const double*> gamma = m_gammaj;
+    const size_t nj = m_ej.numberOfValues();
+    const CellValue<const double> rho = m_rhoj;
+    const CellValue<const double> p = m_pj;
+    const CellValue<const double> gamma = m_gammaj;
 
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-	m_ej[j] = p[j]/(rho[j]*(gamma[j]-1));
+    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const size_t& j){
+        m_ej[j] = p[j]/(rho[j]*(gamma[j]-1));
       });
 
-    const Kokkos::View<const double*> e = m_ej;
-    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-	m_cj[j] = std::sqrt(gamma[j]*(gamma[j]-1)*e[j]);
+    const CellValue<const double> e = m_ej;
+    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const size_t& j){
+        m_cj[j] = std::sqrt(gamma[j]*(gamma[j]-1)*e[j]);
       });
   }
 };
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index fd58774a6..ba0684ad2 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -1,6 +1,9 @@
 #ifndef FINITE_VOLUMES_EULER_UNKNOWNS_HPP
 #define FINITE_VOLUMES_EULER_UNKNOWNS_HPP
 
+#include <TinyVector.hpp>
+#include <ItemValue.hpp>
+
 template <typename TMeshData>
 class FiniteVolumesEulerUnknowns
 {
@@ -15,165 +18,165 @@ private:
   const MeshDataType& m_mesh_data;
   const MeshType& m_mesh;
 
-  Kokkos::View<double*> m_rhoj;
-  Kokkos::View<Rd*> m_uj;
-  Kokkos::View<double*> m_Ej;
+  CellValue<double> m_rhoj;
+  CellValue<Rd> m_uj;
+  CellValue<double> m_Ej;
 
-  Kokkos::View<double*> m_ej;
-  Kokkos::View<double*> m_pj;
-  Kokkos::View<double*> m_gammaj;
-  Kokkos::View<double*> m_cj;
-  Kokkos::View<double*> m_mj;
-  Kokkos::View<double*> m_inv_mj;
+  CellValue<double> m_ej;
+  CellValue<double> m_pj;
+  CellValue<double> m_gammaj;
+  CellValue<double> m_cj;
+  CellValue<double> m_mj;
+  CellValue<double> m_inv_mj;
 
 public:
-  Kokkos::View<double*> rhoj()
+  CellValue<double>& rhoj()
   {
     return m_rhoj;
   }
 
-  const Kokkos::View<const double*> rhoj() const
+  const CellValue<const double> rhoj() const
   {
     return m_rhoj;
   }
 
-  Kokkos::View<Rd*> uj()
+  CellValue<Rd>& uj()
   {
     return m_uj;
   }
 
-  const Kokkos::View<const Rd*> uj() const
+  const CellValue<const Rd> uj() const
   {
     return m_uj;
   }
 
-  Kokkos::View<double*> Ej()
+  CellValue<double>& Ej()
   {
     return m_Ej;
   }
 
-  const Kokkos::View<const double*> Ej() const
+  const CellValue<const double> Ej() const
   {
     return m_Ej;
   }
 
-  Kokkos::View<double*> ej()
+  CellValue<double>& ej()
   {
     return m_ej;
   }
 
-  const Kokkos::View<const double*> ej() const
+  const CellValue<const double> ej() const
   {
     return m_ej;
   }
 
-  Kokkos::View<double*> pj()
+  CellValue<double>& pj()
   {
     return m_pj;
   }
 
-  const Kokkos::View<const double*> pj() const
+  const CellValue<const double> pj() const
   {
     return m_pj;
   }
 
-  Kokkos::View<double*> gammaj()
+  CellValue<double>& gammaj()
   {
     return m_gammaj;
   }
 
-  const Kokkos::View<const double*> gammaj() const
+  const CellValue<const double> gammaj() const
   {
     return m_gammaj;
   }
 
-  Kokkos::View<double*> cj()
+  CellValue<double>& cj()
   {
     return m_cj;
   }
 
-  const Kokkos::View<const double*> cj() const
+  const CellValue<const double> cj() const
   {
     return m_cj;
   }
 
-  Kokkos::View<double*> mj()
+  CellValue<double>& mj()
   {
     return m_mj;
   }
 
-  const Kokkos::View<const double*> mj() const
+  const CellValue<const double> mj() const
   {
     return m_mj;
   }
 
-  Kokkos::View<double*> invMj()
+  CellValue<double>& invMj()
   {
     return m_inv_mj;
   }
 
-  const Kokkos::View<const double*> invMj() const
+  const CellValue<const double> invMj() const
   {
     return m_inv_mj;
   }
 
   void initializeSod()
   {
-    const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
+    const CellValue<const Rd>& xj = m_mesh_data.xj();
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	if (xj[j][0]<0.5) {
-	  m_rhoj[j]=1;
-	} else {
-	  m_rhoj[j]=0.125;
-	}
+        if (xj[j][0]<0.5) {
+          m_rhoj[j]=1;
+        } else {
+          m_rhoj[j]=0.125;
+        }
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	if (xj[j][0]<0.5) {
-	  m_pj[j]=1;
-	} else {
-	  m_pj[j]=0.1;
-	}
+        if (xj[j][0]<0.5) {
+          m_pj[j]=1;
+        } else {
+          m_pj[j]=0.1;
+        }
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	m_uj[j] = zero;
+        m_uj[j] = zero;
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	m_gammaj[j] = 1.4;
+        m_gammaj[j] = 1.4;
       });
 
     BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
     block_eos.updateEandCFromRhoP();
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]);
+        m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]);
       });
 
     const CellValue<const double>& Vj = m_mesh_data.Vj();
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	m_mj[j] = m_rhoj[j] * Vj[j];
+        m_mj[j] = m_rhoj[j] * Vj[j];
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	m_inv_mj[j] = 1./m_mj[j];
+        m_inv_mj[j] = 1./m_mj[j];
       });
   }
 
   FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data)
     : m_mesh_data(mesh_data),
       m_mesh(m_mesh_data.mesh()),
-      m_rhoj("rhoj",m_mesh.numberOfCells()),
-      m_uj("uj",m_mesh.numberOfCells()),
-      m_Ej("Ej",m_mesh.numberOfCells()),
-      m_ej("ej",m_mesh.numberOfCells()),
-      m_pj("pj",m_mesh.numberOfCells()),
-      m_gammaj("gammaj",m_mesh.numberOfCells()),
-      m_cj("cj",m_mesh.numberOfCells()),
-      m_mj("mj",m_mesh.numberOfCells()),
-      m_inv_mj("inv_mj",m_mesh.numberOfCells())
+      m_rhoj(m_mesh.connectivity()),
+      m_uj(m_mesh.connectivity()),
+      m_Ej(m_mesh.connectivity()),
+      m_ej(m_mesh.connectivity()),
+      m_pj(m_mesh.connectivity()),
+      m_gammaj(m_mesh.connectivity()),
+      m_cj(m_mesh.connectivity()),
+      m_mj(m_mesh.connectivity()),
+      m_inv_mj(m_mesh.connectivity())
   {
     ;
   }
-- 
GitLab