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;

  static constexpr size_t dimension = MeshType::dimension;
  static_assert(dimension>0, "dimension must be strictly positive");

  static constexpr double inv_dimension = 1./dimension;

private:
  const MeshType& m_mesh;
  Kokkos::View<Rd**> m_Cjr;
  Kokkos::View<Rd*>  m_xj;
  Kokkos::View<double*>   m_Vj;

  KOKKOS_INLINE_FUNCTION
  void _updateCenter()
  { // Computes vertices isobarycenter
    if(dimension == 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)]);
	});
    } else {
      const Kokkos::View<const Rd*> xr = m_mesh.xr();
      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 double*>& inv_cell_nb_nodes
	= m_mesh.connectivity().invCellNbNodes();
      Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	  Rd X = zero;
	  for (int R=0; R<cell_nb_nodes[j]; ++R) {
	   X += xr[cell_nodes(j,R)];
	  }
	  m_xj[j] = inv_cell_nb_nodes[j]*X;
	});
    }
  }

  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){
	double sum_cjr_xr = 0;
	for (int R=0; R<cell_nb_nodes[j]; ++R) {
	  sum_cjr_xr += (xr[cell_nodes(j,R)], m_Cjr(j,R));
	}
	m_Vj[j] = inv_dimension * sum_cjr_xr;
      });
  }

  KOKKOS_INLINE_FUNCTION
  void _updateCjr() {
    if(dimension == 1) {
      // Cjr are constant overtime
    }
    static_assert(dimension==1, "only 1d is implemented");
  }

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 double*> Vj() const
  {
    return m_Vj;
  }

  void updateAllData()
  {
    this->_updateCenter();
    this->_updateVolume();
    this->_updateCjr();
  }

  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())
  {
    if (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;
	});
    }
    this->updateAllData();
  }

  ~MeshData()
  {
    ;
  }
};

#endif // MESH_DATA_HPP