diff --git a/src/main.cpp b/src/main.cpp
index b1a3decd13b36a4b738c49c9815a0c4bfb327d5f..4bdc52d00d5714866d00bc4965215c74d6e3ec21 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -194,7 +194,7 @@ int main(int argc, char *argv[])
typedef TinyVector<MeshType::dimension> Rd;
- const Kokkos::View<const double*> Vj = mesh_data.Vj();
+ const ValuePerCell<const double>& Vj = mesh_data.Vj();
const double tmax=0.2;
double t=0;
@@ -301,7 +301,7 @@ int main(int argc, char *argv[])
AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
- const Kokkos::View<const double*> Vj = mesh_data.Vj();
+ const ValuePerCell<const double>& Vj = mesh_data.Vj();
const double tmax=0.2;
double t=0;
@@ -397,7 +397,7 @@ int main(int argc, char *argv[])
AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
- const Kokkos::View<const double*> Vj = mesh_data.Vj();
+ const ValuePerCell<const double>& Vj = mesh_data.Vj();
const double tmax=0.2;
double t=0;
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 0fe98256bb87d8d02b29899906e475737f4e835d..b03f8cda8c323a569115802f306d914ff0921929 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -4,6 +4,7 @@
#include <Kokkos_Core.hpp>
#include <TinyVector.hpp>
+#include <ValuePerItem.hpp>
#include <SubItemValuePerItem.hpp>
#include <map>
@@ -27,7 +28,7 @@ class MeshData
NodeValuePerCell<double> m_ljr;
NodeValuePerCell<Rd> m_njr;
Kokkos::View<Rd*> m_xj;
- Kokkos::View<double*> m_Vj;
+ ValuePerCell<const double> m_Vj;
KOKKOS_INLINE_FUNCTION
void _updateCenter()
@@ -70,7 +71,7 @@ class MeshData
const auto& cell_to_node_matrix
= m_mesh.connectivity().getMatrix(ItemType::cell,
ItemType::node);
-
+ ValuePerCell<double> Vj(m_mesh.connectivity());
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
double sum_cjr_xr = 0;
const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
@@ -78,8 +79,9 @@ class MeshData
for (size_t R=0; R<cell_nodes.length; ++R) {
sum_cjr_xr += (xr[cell_nodes(R)], m_Cjr(j,R));
}
- m_Vj[j] = inv_dimension * sum_cjr_xr;
+ Vj[j] = inv_dimension * sum_cjr_xr;
});
+ m_Vj = Vj;
}
KOKKOS_INLINE_FUNCTION
@@ -229,7 +231,7 @@ class MeshData
return m_xj;
}
- const Kokkos::View<const double*> Vj() const
+ const ValuePerCell<const double>& Vj() const
{
return m_Vj;
}
@@ -247,7 +249,7 @@ class MeshData
m_ljr(mesh.connectivity()),
m_njr(mesh.connectivity()),
m_xj("xj", mesh.numberOfCells()),
- m_Vj("Vj", mesh.numberOfCells())
+ m_Vj(mesh.connectivity())
{
if constexpr (dimension==1) {
// in 1d Cjr are computed once for all
diff --git a/src/mesh/ValuePerItem.hpp b/src/mesh/ValuePerItem.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..93efce25420dd9e3bb4d6f2202563feac5e6c701
--- /dev/null
+++ b/src/mesh/ValuePerItem.hpp
@@ -0,0 +1,128 @@
+#ifndef VALUE_PER_ITEM_HPP
+#define VALUE_PER_ITEM_HPP
+
+#include <ItemType.hpp>
+#include <PastisAssert.hpp>
+
+#include <IConnectivity.hpp>
+
+template <typename DataType,
+ ItemType item_type>
+class ValuePerItem
+{
+ public:
+ static const ItemType item_t{item_type};
+
+ private:
+ bool m_is_built{false};
+
+ Kokkos::View<DataType*> m_values;
+
+ // Allows const version to access our data
+ friend ValuePerItem<std::add_const_t<DataType>,
+ item_type>;
+
+ public:
+ KOKKOS_FORCEINLINE_FUNCTION
+ const bool& isBuilt() const
+ {
+ return m_is_built;
+ }
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ DataType& operator()(const size_t& j)
+ {
+ Assert(m_is_built);
+ return m_values(j);
+ }
+
+ // Following Kokkos logic, these classes are view and const view does allow
+ // changes in data
+ KOKKOS_FORCEINLINE_FUNCTION
+ DataType& operator()(const size_t& j) const
+ {
+ Assert(m_is_built);
+ return m_values(j);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ size_t numberOfValues() const
+ {
+ Assert(m_is_built);
+ return m_values.extent(0);
+ }
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ DataType& operator[](const size_t& i)
+ {
+ Assert(m_is_built);
+ return m_values[i];
+ }
+
+ // Following Kokkos logic, these classes are view and const view does allow
+ // changes in data
+ KOKKOS_FORCEINLINE_FUNCTION
+ DataType& operator[](const size_t & i) const
+ {
+ Assert(m_is_built);
+ return m_values[i];
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ size_t numberOfItems() const
+ {
+ Assert(m_is_built);
+ Assert(m_values.extent(0) != 0);
+ return m_values.extent(0);
+ }
+
+ template <typename DataType2>
+ KOKKOS_INLINE_FUNCTION
+ ValuePerItem&
+ operator=(const ValuePerItem<DataType2, item_type>& value_per_item)
+ {
+ // ensures that DataType is the same as source DataType2
+ static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+ "Cannot assign ValuePerItem of different type");
+ // ensures that const is not lost through copy
+ static_assert(((std::is_const<DataType2>() and std::is_const<DataType>())
+ or not std::is_const<DataType2>()),
+ "Cannot assign const ValuePerItem to a non const ValuePerItem");
+
+ m_values = value_per_item.m_values;
+ m_is_built = value_per_item.m_is_built;
+
+ return *this;
+ }
+
+ template <typename DataType2>
+ KOKKOS_INLINE_FUNCTION
+ ValuePerItem(const ValuePerItem<DataType2, item_type>& value_per_item)
+ {
+ this->operator=(value_per_item);
+ }
+
+ ValuePerItem() = default;
+
+ ValuePerItem(const IConnectivity& connectivity)
+ : m_is_built{true}
+ {
+ m_values = Kokkos::View<std::remove_const_t<DataType>*>("values", connectivity.numberOf<item_type>());
+ }
+
+ ~ValuePerItem() = default;
+};
+
+template <typename DataType>
+using ValuePerNode = ValuePerItem<DataType, ItemType::node>;
+
+template <typename DataType>
+using ValuePerEdge = ValuePerItem<DataType, ItemType::edge>;
+
+template <typename DataType>
+using ValuePerFace = ValuePerItem<DataType, ItemType::face>;
+
+template <typename DataType>
+using ValuePerCell = ValuePerItem<DataType, ItemType::cell>;
+
+#endif // VALUE_PER_ITEM_HPP
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 0cc6ec0df7a97b91999376873ed7b35a0c5acf6c..2732799763ea7dc1510e5f4c71078ffad5564532 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -239,7 +239,7 @@ class AcousticSolver
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 ValuePerCell<const double>& Vj,
const NodeValuePerCell<Rd>& Cjr,
const NodeValuePerCell<double>& ljr,
const NodeValuePerCell<Rd>& njr) {
@@ -286,7 +286,7 @@ class AcousticSolver
}
KOKKOS_INLINE_FUNCTION
- double acoustic_dt(const Kokkos::View<const double*>& Vj,
+ double acoustic_dt(const ValuePerCell<const double>& Vj,
const Kokkos::View<const double*>& cj) const
{
const NodeValuePerCell<double>& ljr = m_mesh_data.ljr();
@@ -324,7 +324,7 @@ class AcousticSolver
Kokkos::View<double*> cj = unknowns.cj();
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
- const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
+ const ValuePerCell<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();
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 617cb8b39759dcf4ae4ff6f4a41f09a7bc056584..aafb44aaaf47f3d275169e14ab447e2f5b90035c 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -152,7 +152,7 @@ public:
m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]);
});
- const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
+ const ValuePerCell<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];
});