diff --git a/src/main.cpp b/src/main.cpp
index 4ffdd9b0ab49363011b0a963941fdc9b8bff231b..ba6f6c504ae02f7859d722de46295e192f92fffe 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -236,7 +236,7 @@ int main(int argc, char *argv[])
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) {
+ for (CellId j=0; j<mesh.numberOfCells(); ++j) {
fout << xj[j][0] << ' ' << rhoj[j] << '\n';
}
}
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 6e19dcc6a99bf426247d4f40ed3fe0fa9ca96ed0..e4e003d79531daf4420e9a18662157a017f01d89 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -11,7 +11,7 @@ void Connectivity<3>::_computeCellFaceAndFaceNodeConnectivities()
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) {
+ for (CellId j=0; j<this->numberOfCells(); ++j) {
const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
switch (m_cell_type[j]) {
@@ -99,7 +99,7 @@ void Connectivity<3>::_computeCellFaceAndFaceNodeConnectivities()
auto& cell_to_face_matrix
= m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)];
std::vector<std::vector<unsigned int>> cell_to_face_vector(this->numberOfCells());
- for (size_t j=0; j<cell_to_face_vector.size(); ++j) {
+ for (CellId j=0; j<cell_to_face_vector.size(); ++j) {
cell_to_face_vector[j].resize(cell_nb_faces[j]);
}
int l=0;
@@ -120,7 +120,7 @@ void Connectivity<3>::_computeCellFaceAndFaceNodeConnectivities()
const auto& cells_vector = face_cells_vector.second;
for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
- cell_face_is_reversed(cell_number, cell_local_face) = reversed;
+ cell_face_is_reversed(CellId{cell_number}, cell_local_face) = reversed;
}
}
@@ -228,14 +228,14 @@ Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
{
CellValue<CellType> cell_type(*this);
- Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
cell_type[j] = cell_type_vector[j];
});
m_cell_type = cell_type;
}
{
CellValue<double> inv_cell_nb_nodes(*this);
- Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const CellId& 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/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp
index 0c87ccb6d64f840cd450eb5ded67c6a337a5cf25..672502070399313bdcaf17df068ab42675bb3f01 100644
--- a/src/mesh/ConnectivityComputer.cpp
+++ b/src/mesh/ConnectivityComputer.cpp
@@ -88,7 +88,7 @@ ConnectivityComputer::computeLocalItemNumberInChildItem(const ConnectivityType&
= connectivity.getMatrix(item_type, child_item_type);
SubItemValuePerItem<unsigned short, child_item_type, item_type> item_number_in_child_item(connectivity);
- for (unsigned int r=0; r<item_to_child_items_matrix.numRows(); ++r) {
+ for (ItemIdT<item_type> r=0; r<item_to_child_items_matrix.numRows(); ++r) {
const auto& item_to_child_items = item_to_child_items_matrix.rowConst(r);
for (unsigned short J=0; J<item_to_child_items.length; ++J) {
const unsigned int j = item_to_child_items(J);
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index b14e9c493ff6b55d95402acb482a1a9a7f6fd665..be77c6ac16eefb092c858ac80231747d4e3ae0a0 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -841,7 +841,7 @@ GmshReader::__proceedData()
using Rd = TinyVector<3, double>;
NodeValue<Rd> xr(connectivity);
- for (size_t i=0; i<__vertices.size(); ++i) {
+ for (NodeId i=0; i<__vertices.size(); ++i) {
xr[i][0] = __vertices[i][0];
xr[i][1] = __vertices[i][1];
xr[i][2] = __vertices[i][2];
@@ -913,7 +913,7 @@ GmshReader::__proceedData()
using Rd = TinyVector<2, double>;
NodeValue<Rd> xr(connectivity);
- for (size_t i=0; i<__vertices.size(); ++i) {
+ for (NodeId i=0; i<__vertices.size(); ++i) {
xr[i][0] = __vertices[i][0];
xr[i][1] = __vertices[i][1];
}
@@ -957,7 +957,7 @@ GmshReader::__proceedData()
using Rd = TinyVector<1, double>;
NodeValue<Rd> xr(connectivity);
- for (size_t i=0; i<__vertices.size(); ++i) {
+ for (NodeId i=0; i<__vertices.size(); ++i) {
xr[i][0] = __vertices[i][0];
}
diff --git a/src/mesh/ItemId.hpp b/src/mesh/ItemId.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..20f538a0f2e1b83f8e3c4426458fe847c82341f8
--- /dev/null
+++ b/src/mesh/ItemId.hpp
@@ -0,0 +1,78 @@
+#ifndef ITEM_ID_HPP
+#define ITEM_ID_HPP
+
+#include <ItemType.hpp>
+
+template <ItemType item_type>
+class ItemIdT
+{
+ private:
+ unsigned int m_id;
+
+ public:
+ operator const unsigned int&() const
+ {
+ return m_id;
+ }
+
+ ItemIdT operator++(int)
+ {
+ ItemIdT item_id(m_id);
+ ++m_id;
+ return std::move(item_id);
+ }
+
+ ItemIdT& operator++()
+ {
+ ++m_id;
+ return *this;
+ }
+
+ ItemIdT& operator=(const unsigned int& id)
+ {
+ m_id = id;
+ return *this;
+ }
+
+ ItemIdT& operator=(const ItemIdT&) = default;
+ ItemIdT& operator=(ItemIdT&&) = default;
+
+ ItemIdT(const unsigned int& id) : m_id{id}{}
+
+ ItemIdT(const ItemIdT&) = default;
+
+ ItemIdT(ItemIdT&&) = default;
+ ItemIdT() = default;
+ ~ItemIdT() = default;
+
+ // forbidden rules
+ template <ItemType other_item_type>
+ ItemIdT(const ItemIdT<other_item_type>&)
+ {
+ static_assert(other_item_type == item_type,
+ "wrong type of item");
+ }
+
+ template <ItemType other_item_type>
+ ItemIdT& operator=(const ItemIdT<other_item_type>&)
+ {
+ static_assert(other_item_type == item_type,
+ "wrong type of item");
+ }
+
+
+ template <ItemType other_item_type>
+ ItemIdT& operator=(ItemIdT<other_item_type>&&)
+ {
+ static_assert(other_item_type == item_type,
+ "wrong type of item");
+ }
+
+};
+
+using NodeId = ItemIdT<ItemType::node>;
+using EdgeId = ItemIdT<ItemType::edge>;
+using FaceId = ItemIdT<ItemType::face>;
+using CellId = ItemIdT<ItemType::cell>;
+
+#endif // ITEM_ID_HPP
diff --git a/src/mesh/ItemValue.hpp b/src/mesh/ItemValue.hpp
index 6f951654308070ee5d47f603f5ae1b0bcc8004c7..77eb1cb37d9998fb1af6e4a27d05ff7fcdb60077 100644
--- a/src/mesh/ItemValue.hpp
+++ b/src/mesh/ItemValue.hpp
@@ -1,11 +1,12 @@
#ifndef ITEM_VALUE_HPP
#define ITEM_VALUE_HPP
-//#include <Kokkos_Core.hpp>
+#include <PastisAssert.hpp>
+
#include <Array.hpp>
#include <ItemType.hpp>
-#include <PastisAssert.hpp>
+#include <ItemId.hpp>
#include <IConnectivity.hpp>
@@ -17,6 +18,9 @@ class ItemValue
static const ItemType item_t{item_type};
using data_type = DataType;
+ using ItemId = ItemIdT<item_type>;
+ using index_type = ItemId;
+
private:
bool m_is_built{false};
@@ -43,12 +47,20 @@ class ItemValue
// 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
+ DataType& operator[](const ItemId& i) const
{
Assert(m_is_built);
return m_values[i];
}
+ template <typename IndexType>
+ DataType& operator[](const IndexType& i) const
+ {
+ static_assert(std::is_same<IndexType,ItemId>(),
+ "ItemValue must be indexed by ItemId");
+ return m_values[i];
+ }
+
KOKKOS_INLINE_FUNCTION
size_t numberOfItems() const
{
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 70c73a004383732f618f5a750781cc7e4977a291..2c86899269dfbbcb1d01c679b136547341bc1078 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -41,9 +41,9 @@ class MeshData
ItemType::node);
CellValue<Rd> xj(m_mesh.connectivity());
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
- xj[j] = 0.5*(xr[cell_nodes(0)]+xr[cell_nodes(1)]);
+ xj[j] = 0.5*(xr[NodeId{cell_nodes(0)}]+xr[NodeId{cell_nodes(1)}]);
});
m_xj = xj;
} else {
@@ -56,11 +56,11 @@ class MeshData
= 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){
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& 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)];
+ X += xr[NodeId{cell_nodes(R)}];
}
xj[j] = inv_cell_nb_nodes[j]*X;
});
@@ -76,12 +76,12 @@ class MeshData
= m_mesh.connectivity().getMatrix(ItemType::cell,
ItemType::node);
CellValue<double> Vj(m_mesh.connectivity());
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
double sum_cjr_xr = 0;
const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
for (size_t R=0; R<cell_nodes.length; ++R) {
- sum_cjr_xr += (xr[cell_nodes(R)], m_Cjr(j,R));
+ sum_cjr_xr += (xr[NodeId{cell_nodes(R)}], m_Cjr(j,R));
}
Vj[j] = inv_dimension * sum_cjr_xr;
});
@@ -101,12 +101,12 @@ class MeshData
{
NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& 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)]);
+ Rd half_xrp_xrm = 0.5*(xr[NodeId{cell_nodes(Rp1)}]-xr[NodeId{cell_nodes(Rm1)}]);
Cjr(j,R) = Rd{-half_xrp_xrm[1], half_xrp_xrm[0]};
}
});
@@ -115,16 +115,16 @@ class MeshData
{
NodeValuePerCell<double> ljr(m_mesh.connectivity());
- Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& j){
- ljr[j] = l2Norm(m_Cjr[j]);
+ Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const size_t& 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& j){
- njr[j] = (1./m_ljr[j])*m_Cjr[j];
+ Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const size_t& jr){
+ njr[jr] = (1./m_ljr[jr])*m_Cjr[jr];
});
m_njr = njr;
}
@@ -136,22 +136,24 @@ class MeshData
= m_mesh.connectivity().getMatrix(ItemType::face,
ItemType::node);
- Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
+ Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const FaceId& l) {
const auto& face_nodes = face_to_node_matrix.rowConst(l);
const size_t nb_nodes = face_nodes.length;
std::vector<Rd> dxr(nb_nodes);
for (size_t r=0; r<nb_nodes; ++r) {
- dxr[r] = xr[face_nodes((r+1)%nb_nodes)] - xr[face_nodes((r+nb_nodes-1)%nb_nodes)];
+ dxr[r]
+ = xr[NodeId{face_nodes((r+1)%nb_nodes)}]
+ - xr[NodeId{face_nodes((r+nb_nodes-1)%nb_nodes)}];
}
const double inv_12_nb_nodes = 1./(12.*nb_nodes);
for (size_t r=0; r<nb_nodes; ++r) {
Rd Nr = zero;
const Rd two_dxr = 2*dxr[r];
for (size_t s=0; s<nb_nodes; ++s) {
- Nr += crossProduct((two_dxr - dxr[s]), xr[face_nodes(s)]);
+ Nr += crossProduct((two_dxr - dxr[s]), xr[NodeId{face_nodes(s)}]);
}
Nr *= inv_12_nb_nodes;
- Nr -= (1./6.)*crossProduct(dxr[r], xr[face_nodes(r)]);
+ Nr -= (1./6.)*crossProduct(dxr[r], xr[NodeId{face_nodes(r)}]);
Nlr(l,r) = Nr;
}
});
@@ -168,11 +170,11 @@ class MeshData
{
NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
- Kokkos::parallel_for(Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
+ Kokkos::parallel_for(Cjr.numberOfValues(), KOKKOS_LAMBDA(const size_t& jr){
Cjr[jr] = zero;
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
const auto& cell_faces = cell_to_face_matrix.rowConst(j);
@@ -213,7 +215,7 @@ class MeshData
{
NodeValuePerCell<double> ljr(m_mesh.connectivity());
- Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
+ Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const size_t& jr){
ljr[jr] = l2Norm(m_Cjr[jr]);
});
m_ljr = ljr;
@@ -221,7 +223,7 @@ class MeshData
{
NodeValuePerCell<Rd> njr(m_mesh.connectivity());
- Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
+ Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const size_t& jr){
njr[jr] = (1./m_ljr[jr])*m_Cjr[jr];
});
m_njr = njr;
@@ -275,7 +277,7 @@ class MeshData
// in 1d Cjr are computed once for all
{
NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
Cjr(j,0)=-1;
Cjr(j,1)= 1;
});
@@ -285,7 +287,7 @@ class MeshData
m_njr = m_Cjr;
{
NodeValuePerCell<double> ljr(m_mesh.connectivity());
- Kokkos::parallel_for(ljr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
+ Kokkos::parallel_for(ljr.numberOfValues(), KOKKOS_LAMBDA(const size_t& jr){
ljr[jr] = 1;
});
m_ljr = ljr;
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 16c016166eb441540701243bb04a99338a0dc072..768671115bfb3def2e05add89413e425434ff322 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -157,7 +157,7 @@ _checkBoundaryIsFlat(const TinyVector<2,double>& normal,
const NodeValue<const R2>& xr = mesh.xr();
Kokkos::parallel_for(m_node_list.size(), KOKKOS_LAMBDA(const size_t& r) {
- const R2& x = xr[m_node_list[r]];
+ const R2& x = xr[NodeId{m_node_list[r]}];
if ((x-origin,normal)>1E-13*length) {
std::cerr << "this FlatBoundary is not flat!\n";
std::exit(1);
@@ -200,7 +200,7 @@ _getNormal(const MeshType& mesh)
-std::numeric_limits<double>::max());
for (size_t r=0; r<m_node_list.size(); ++r) {
- const R2& x = xr[m_node_list[r]];
+ const R2& x = xr[NodeId{m_node_list[r]}];
if ((x[0]<xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
xmin = x;
}
@@ -247,7 +247,7 @@ _getNormal(const MeshType& mesh)
const NodeValue<const R3>& xr = mesh.xr();
for (size_t r=0; r<m_node_list.size(); ++r) {
- const R3& x = xr[m_node_list[r]];
+ const R3& x = xr[NodeId{m_node_list[r]}];
if (x[0]<xmin[0]) {
xmin = x;
}
@@ -333,7 +333,7 @@ _getOutgoingNormal(const MeshType& mesh)
const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
double max_height = 0;
for (size_t r=0; r<j0_nodes.length; ++r) {
- const double height = (xr[j0_nodes(r)]-xr[r0], normal);
+ const double height = (xr[NodeId{j0_nodes(r)}]-xr[NodeId{r0}], normal);
if (std::abs(height) > std::abs(max_height)) {
max_height = height;
}
@@ -370,7 +370,7 @@ _getOutgoingNormal(const MeshType& mesh)
const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
double max_height = 0;
for (size_t r=0; r<j0_nodes.length; ++r) {
- const double height = (xr[j0_nodes(r)]-xr[r0], normal);
+ const double height = (xr[NodeId{j0_nodes(r)}]-xr[NodeId{r0}], normal);
if (std::abs(height) > std::abs(max_height)) {
max_height = height;
}
@@ -407,7 +407,7 @@ _getOutgoingNormal(const MeshType& mesh)
const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
double max_height = 0;
for (size_t r=0; r<j0_nodes.length; ++r) {
- const double height = (xr[j0_nodes(r)]-xr[r0], normal);
+ const double height = (xr[NodeId{j0_nodes(r)}]-xr[NodeId{r0}], normal);
if (std::abs(height) > std::abs(max_height)) {
max_height = height;
}
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index b578fc5faae3d94659be87e29a2d13e54e4347ef..1c9a46ac0ef7b5a55723049dc177915dfee7e886 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -3,9 +3,12 @@
#include <Kokkos_StaticCrsGraph.hpp>
+#include <PastisAssert.hpp>
+
#include <Array.hpp>
#include <ItemType.hpp>
-#include <PastisAssert.hpp>
+
+#include <ItemId.hpp>
#include <ConnectivityMatrix.hpp>
#include <IConnectivity.hpp>
@@ -29,11 +32,13 @@ class SubItemValuePerItem<DataType,
static const ItemType sub_item_t{sub_item_type};
using data_type = DataType;
+ using ItemId = ItemIdT<item_type>;
+ using index_type = ItemId;
private:
bool m_is_built{false};
- ConnectivityMatrix::HostRowType m_host_row_map;
+ typename ConnectivityMatrix::HostRowType m_host_row_map;
Array<DataType> m_values;
// Allows const version to access our data
@@ -98,10 +103,19 @@ class SubItemValuePerItem<DataType,
// 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 size_t& r) const
+ DataType& operator()(const ItemId& j, const size_t& r) const
{
Assert(m_is_built);
- return m_values[m_host_row_map(j)+r];
+ return m_values[m_host_row_map(size_t{j})+r];
+ }
+
+ template <typename IndexType>
+ KOKKOS_FORCEINLINE_FUNCTION
+ DataType& operator()(const IndexType& j, const size_t& r) const
+ {
+ static_assert(std::is_same<IndexType, size_t>(),
+ "SubItemValuePerItem indexed by ItemId");
+ return m_values[m_host_row_map(size_t{j})+r];
}
KOKKOS_INLINE_FUNCTION
@@ -114,12 +128,20 @@ class SubItemValuePerItem<DataType,
// 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
+ DataType& operator[](const size_t& i) const
{
Assert(m_is_built);
return m_values[i];
}
+ template <typename IndexType>
+ DataType& operator[](const IndexType& i) const
+ {
+ static_assert(std::is_same<IndexType, size_t>(),
+ "Access to SubItemValuePerItem's array must be indexed by size_t");
+ return m_values[i];
+ }
+
KOKKOS_INLINE_FUNCTION
size_t numberOfItems() const
{
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index 24865f43decc9ccce1514d619fa01a871df13176..f39fd3415a00d3d92e8e72437a600485cefd275d 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -44,21 +44,21 @@ class VTKWriter
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 (NodeId 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) {
- for (unsigned int r=0; r<mesh.numberOfNodes(); ++r) {
+ for (NodeId r=0; r<mesh.numberOfNodes(); ++r) {
for (unsigned short i=0; i<2; ++i) {
fout << xr[r][i] << ' ';
}
fout << "0 "; // VTK requires 3 components
}
} else {
- for (unsigned int r=0; r<mesh.numberOfNodes(); ++r) {
+ for (NodeId 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 fe9801155b587510d13a00e479df9a76ab74f660..6139b022b3fc5b15f47418cf4882f189e212d4b0 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -39,7 +39,7 @@ class AcousticSolver
computeRhoCj(const CellValue<const double>& rhoj,
const CellValue<const double>& cj)
{
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
m_rhocj[j] = rhoj[j]*cj[j];
});
return m_rhocj;
@@ -51,7 +51,7 @@ class AcousticSolver
const NodeValuePerCell<const double>& ljr,
const NodeValuePerCell<const Rd>& njr)
{
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
const size_t& nb_nodes =m_Ajr.numberOfSubValues(j);
const double& rho_c = rhocj[j];
for (size_t r=0; r<nb_nodes; ++r) {
@@ -69,14 +69,14 @@ class AcousticSolver
const auto& node_local_numbers_in_their_cells
= m_connectivity.nodeLocalNumbersInTheirCells();
- Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
+ Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const NodeId& r) {
Rdd sum = zero;
const auto& node_to_cell = node_to_cell_matrix.rowConst(r);
const auto& node_local_number_in_its_cells
= node_local_numbers_in_their_cells.itemValues(r);
for (size_t j=0; j<node_to_cell.length; ++j) {
- const unsigned int J = node_to_cell(j);
+ const CellId J = node_to_cell(j);
const unsigned int R = node_local_number_in_its_cells[j];
sum += Ajr(J,R);
}
@@ -98,14 +98,14 @@ class AcousticSolver
const auto& node_local_numbers_in_their_cells
= m_connectivity.nodeLocalNumbersInTheirCells();
- Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
+ Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const NodeId& r) {
Rd& br = m_br[r];
br = zero;
const auto& node_to_cell = node_to_cell_matrix.rowConst(r);
const auto& node_local_number_in_its_cells
= node_local_numbers_in_their_cells.itemValues(r);
for (size_t j=0; j<node_to_cell.length; ++j) {
- const unsigned int J = node_to_cell(j);
+ const CellId J = node_to_cell(j);
const unsigned int R = node_local_number_in_its_cells[j];
br += Ajr(J,R)*uj[J] + pj[J]*Cjr(J,R);
}
@@ -147,7 +147,7 @@ class AcousticSolver
const Array<const unsigned int>& node_list
= symmetry_bc.nodeList();
Kokkos::parallel_for(symmetry_bc.numberOfNodes(), KOKKOS_LAMBDA(const int& r_number) {
- const int r = node_list[r_number];
+ const NodeId r = node_list[r_number];
m_Ar[r] = P*m_Ar[r]*P + nxn;
m_br[r] = P*m_br[r];
@@ -164,7 +164,7 @@ class AcousticSolver
{
inverse(Ar, m_inv_Ar);
const NodeValue<const Rdd> invAr = m_inv_Ar;
- Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
+ Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const NodeId& r) {
m_ur[r]=invAr[r]*br[r];
});
@@ -182,11 +182,11 @@ class AcousticSolver
= m_mesh.connectivity().getMatrix(ItemType::cell,
ItemType::node);
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
for (size_t r=0; r<cell_nodes.length; ++r) {
- m_Fjr(j,r) = Ajr(j,r)*(uj[j]-ur[cell_nodes(r)])+pj[j]*Cjr(j,r);
+ m_Fjr(j,r) = Ajr(j,r)*(uj[j]-ur[NodeId{cell_nodes(r)}])+pj[j]*Cjr(j,r);
}
});
}
@@ -194,7 +194,7 @@ class AcousticSolver
void inverse(const NodeValue<const Rdd>& A,
NodeValue<Rdd>& inv_A) const
{
- Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
+ Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const NodeId& r) {
inv_A[r] = ::inverse(A[r]);
});
}
@@ -262,7 +262,7 @@ class AcousticSolver
= m_mesh.connectivity().getMatrix(ItemType::cell,
ItemType::node);
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
double S = 0;
@@ -275,7 +275,6 @@ class AcousticSolver
return ReduceMin(m_Vj_over_cj);
}
-
void computeNextStep(const double& t, const double& dt,
UnknownsType& unknowns)
{
@@ -303,13 +302,13 @@ class AcousticSolver
ItemType::node);
const CellValue<const double> inv_mj = unknowns.invMj();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
Rd momentum_fluxes = zero;
double energy_fluxes = 0;
for (size_t R=0; R<cell_nodes.length; ++R) {
- const unsigned int r=cell_nodes(R);
+ const NodeId r=cell_nodes(R);
momentum_fluxes += Fjr(j,R);
energy_fluxes += (Fjr(j,R), ur[r]);
}
@@ -317,18 +316,18 @@ class AcousticSolver
Ej[j] -= (dt*inv_mj[j]) * energy_fluxes;
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
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){
+ Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const NodeId& r){
mutable_xr[r] += dt*ur[r];
});
m_mesh_data.updateAllData();
const CellValue<const double> mj = unknowns.mj();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
rhoj[j] = mj[j]/Vj[j];
});
}
diff --git a/src/scheme/BlockPerfectGas.hpp b/src/scheme/BlockPerfectGas.hpp
index 0ec05d266f46fe80e38ea9bd4153f19996f4ab12..704787756bc97ac60b48699c8215c9444cb87bf2 100644
--- a/src/scheme/BlockPerfectGas.hpp
+++ b/src/scheme/BlockPerfectGas.hpp
@@ -30,11 +30,11 @@ public:
void updatePandCFromRhoE()
{
const size_t nj = m_ej.size();
- const CellValue<const double> rho = m_rhoj;
- const CellValue<const double> e = m_ej;
- const CellValue<const double> gamma = m_gammaj;
+ 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){
+ Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const CellId& 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]);
@@ -44,16 +44,16 @@ public:
void updateEandCFromRhoP()
{
const size_t nj = m_ej.size();
- const CellValue<const double> rho = m_rhoj;
- const CellValue<const double> p = m_pj;
- const CellValue<const double> gamma = m_gammaj;
+ 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 size_t& j){
+ Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const CellId& j){
m_ej[j] = p[j]/(rho[j]*(gamma[j]-1));
});
- const CellValue<const double> e = m_ej;
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const size_t& j){
+ const CellValue<const double>& e = m_ej;
+ Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const CellId& 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 e20317d19811b9fa3a454719091d31b371efa6b6..2eaaf474dbbe520da30cc2c5b0ad5b4737e12538 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -124,7 +124,7 @@ public:
{
const CellValue<const Rd>& xj = m_mesh_data.xj();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
if (xj[j][0]<0.5) {
m_rhoj[j]=1;
} else {
@@ -132,7 +132,7 @@ public:
}
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
if (xj[j][0]<0.5) {
m_pj[j]=1;
} else {
@@ -140,27 +140,27 @@ public:
}
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
m_uj[j] = zero;
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
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){
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& 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){
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
m_mj[j] = m_rhoj[j] * Vj[j];
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
m_inv_mj[j] = 1./m_mj[j];
});
}
diff --git a/src/utils/Array.hpp b/src/utils/Array.hpp
index 410fce99af45005386d9cd5504992bd0815428fc..e72efb042c5c9006960a15dfd20c87ed29c29d82 100644
--- a/src/utils/Array.hpp
+++ b/src/utils/Array.hpp
@@ -9,6 +9,7 @@ class Array
{
public:
using data_type = DataType;
+ using index_type = size_t;
private:
Kokkos::View<DataType*> m_values;
@@ -24,7 +25,7 @@ class Array
}
KOKKOS_INLINE_FUNCTION
- DataType& operator[](const size_t& i) const
+ DataType& operator[](const index_type& i) const
{
Assert(i<m_values.extent(0));
return m_values[i];
diff --git a/src/utils/ArrayUtils.hpp b/src/utils/ArrayUtils.hpp
index 7cb9915369d08b3e89243d48a8370a957a2c3541..9a55abfebde94b4ab53adff473c52907ec5385f8 100644
--- a/src/utils/ArrayUtils.hpp
+++ b/src/utils/ArrayUtils.hpp
@@ -9,6 +9,7 @@ class ReduceMin
private:
const ArrayType& m_array;
using data_type = std::remove_const_t<typename ArrayType::data_type>;
+ using index_type = typename ArrayType::index_type;
public:
KOKKOS_INLINE_FUNCTION
@@ -20,7 +21,7 @@ class ReduceMin
}
KOKKOS_INLINE_FUNCTION
- void operator()(const size_t& i, data_type& data) const
+ void operator()(const index_type& i, data_type& data) const
{
if (m_array[i] < data) {
data = m_array[i];
@@ -60,6 +61,7 @@ class ReduceMax
private:
const ArrayType& m_array;
using data_type = std::remove_const_t<typename ArrayType::data_type>;
+ using index_type = typename ArrayType::index_type;
public:
KOKKOS_INLINE_FUNCTION
@@ -71,7 +73,7 @@ class ReduceMax
}
KOKKOS_INLINE_FUNCTION
- void operator()(const size_t& i, data_type& data) const
+ void operator()(const index_type& i, data_type& data) const
{
if (m_array[i] > data) {
data = m_array[i];