added MeshData class (contains updated geometric info)

experimental/AcousticSolverWithMesh.hpp

@@ -9,11 +9,14 @@
 #include <TinyVector.hpp>
 #include <TinyMatrix.hpp>
 #include <Mesh.hpp>
+#include <MeshData.hpp>
 
-
-template<typename MeshType>
+template<typename MeshData>
 class AcousticSolverWithMesh 
 {
+  typedef typename MeshData::MeshType MeshType;
+
+  MeshData& m_mesh_data;
   const MeshType& m_mesh;
   const typename MeshType::Connectivity& m_connectivity;
 
@@ -225,11 +228,12 @@ private:
 
 
 public:
-  AcousticSolverWithMesh(MeshType& mesh)
-    : m_mesh(mesh),
+  AcousticSolverWithMesh(MeshData& mesh_data)
+    : m_mesh_data(mesh_data),
+      m_mesh(mesh_data.mesh()),
       m_connectivity(m_mesh.connectivity()),
-      m_nj(mesh.numberOfCells()),
-      m_nr(mesh.numberOfNodes()),
+      m_nj(m_mesh.numberOfCells()),
+      m_nr(m_mesh.numberOfNodes()),
       m_br("br", m_nr),
       m_Ajr("Ajr", m_nj, m_connectivity.maxNbNodePerCell()),
       m_Ar("Ar", m_nr),
@@ -239,7 +243,6 @@ public:
       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);
@@ -247,28 +250,17 @@ public:
     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 = mesh.xr();
-
-    const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes();
-    Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){
-	xj[j] = 0.5*(xr[cell_nodes(j,0)]+xr[cell_nodes(j,1)]);
-      });
+    Kokkos::View<Rd*> xr = m_mesh.xr();
 
-    Kokkos::View<Rd*[2]> Cjr("Cjr",m_nj);
-    Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) {
-	Cjr(j,0)=-1;
-	Cjr(j,1)= 1;
-      });
+    const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
+    const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
 
-    Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){
-	Vj[j] = (xr[cell_nodes(j,1)], Cjr(j,1))
-	  + (xr[cell_nodes(j,0)], Cjr(j,0));
-      });
+    const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes();
+    const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
 
     Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){
 	if (xj[j][0]<0.5) {
@@ -351,15 +343,7 @@ public:
 	  xr[r] += dt*ur[r];
 	});
 
-      Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){
-	  xj[j] = 0.5*(xr[j]+xr[j+1]);
-	});
-
-      Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){
-	  Vj[j]
-	    = (xr[cell_nodes(j,1)], Cjr(j,1))
-	    + (xr[cell_nodes(j,0)], Cjr(j,0));
-	});
+      m_mesh_data.updateAllData();
 
       Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){
 	  rhoj[j] = mj[j]/Vj[j];

experimental/Mesh.hpp

@@ -39,7 +39,7 @@ public:
   }
 
   #warning PASSER CES NOUVEAUX VECTEURS EN CONST
-  Kokkos::View<Rd*> xr()
+  Kokkos::View<Rd*> xr() const
   {
     return m_xr;
   }

experimental/MeshData.hpp

+#ifndef MESH_DATA_HPP
+#define MESH_DATA_HPP
+
+#include <Kokkos_Core.hpp>
+#include <TinyVector.hpp>
+
+template <typename M>
+class MeshData
+{
+public:
+  typedef M MeshType;
+
+  typedef TinyVector<dimension> Rd;
+  typedef double R;
+
+  static constexpr size_t dimension = MeshType::dimension;
+  static_assert(dimension>0, "dimension must be strictly positive");
+
+  static constexpr R inv_dimension = 1./dimension;
+
+private:
+  const MeshType& m_mesh;
+  Kokkos::View<Rd**> m_Cjr;
+  Kokkos::View<Rd*>  m_xj;
+  Kokkos::View<R*>   m_Vj;
+
+  template<size_t Dim>
+  KOKKOS_INLINE_FUNCTION
+  void _updateCenter();
+
+  template <>
+  KOKKOS_INLINE_FUNCTION
+  void _updateCenter<1>()
+  {
+    const Kokkos::View<const Rd*> xr = m_mesh.xr();
+
+    const Kokkos::View<const unsigned int**>& cell_nodes
+      = m_mesh.connectivity().cellNodes();
+    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+	m_xj[j] = 0.5*(xr[cell_nodes(j,0)]+xr[cell_nodes(j,1)]);
+      });
+  }
+
+  template<size_t Dim>
+  KOKKOS_INLINE_FUNCTION
+  void _updateVolume()
+  {
+    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();
+
+    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+	R 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;
+      });
+  }
+
+  template<>
+  KOKKOS_INLINE_FUNCTION
+  void _updateVolume<1>()
+  {
+    const Kokkos::View<const Rd*> xr = m_mesh.xr();
+
+    const Kokkos::View<const unsigned int**>& cell_nodes
+      = m_mesh.connectivity().cellNodes();
+
+    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+	m_Vj[j]
+	  = (xr[cell_nodes(j,1)], m_Cjr(j,1))
+	  + (xr[cell_nodes(j,0)], m_Cjr(j,0));
+      });
+  }
+
+  template<size_t Dim>
+  KOKKOS_INLINE_FUNCTION
+  void _updateCjr();
+
+  template<>
+  KOKKOS_INLINE_FUNCTION
+  void _updateCjr<1>() {}
+
+public:
+
+  const MeshType& mesh() const
+  {
+    return m_mesh;
+  }
+
+
+  const Kokkos::View<const Rd**> Cjr() const
+  {
+    return m_Cjr;
+  }
+
+  const Kokkos::View<const Rd*> xj() const
+  {
+    return m_xj;
+  }
+
+  const Kokkos::View<const R*> Vj() const
+  {
+    return m_Vj;
+  }
+
+  void updateAllData()
+  {
+    this->_updateCenter<dimension>();
+    this->_updateVolume<dimension>();
+    this->_updateCjr<dimension>();
+  }
+
+  MeshData(const MeshType& mesh)
+    : m_mesh(mesh),
+      m_Cjr("Cjr", mesh.numberOfCells(), mesh.connectivity().maxNbNodePerCell()),
+      m_xj("xj", mesh.numberOfCells()),
+      m_Vj("Vj", mesh.numberOfCells())
+  {
+    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+	m_Cjr(j,0)=-1;
+	m_Cjr(j,1)= 1;
+      });
+    this->updateAllData();
+  }
+
+  ~MeshData()
+  {
+    ;
+  }
+};
+
+#endif // MESH_DATA_HPP

main.cpp

@@ -118,8 +118,13 @@ int main(int argc, char *argv[])
     Kokkos::Timer timer;
     timer.reset();
     Connectivity1D connectivity(number);
-    Mesh<Connectivity1D> mesh(connectivity);
-    AcousticSolverWithMesh<Mesh<Connectivity1D>> acoustic_solver(mesh);
+    typedef Mesh<Connectivity1D> MeshType;
+    typedef MeshData<MeshType> MeshDataType;
+
+    MeshType mesh(connectivity);
+    MeshDataType mesh_data(mesh);
+    AcousticSolverWithMesh<MeshDataType> acoustic_solver(mesh_data);
+
     method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
   }