Use SubItemValuePerItem to store local sub item number in item

src/mesh/Connectivity.cpp

@@ -150,6 +150,8 @@ void Connectivity<3>::_computeFaceCellConnectivities()
   m_connectivity_computer.computeInverseConnectivityMatrix(m_cell_to_face_matrix,
                                                            m_face_to_cell_matrix);
 
+  m_face_to_cell_local_face_matrix = CellValuePerFace<unsigned short>(*this);
+
   m_connectivity_computer.computeLocalChildItemNumberInItem(m_cell_to_face_matrix,
                                                             m_face_to_cell_matrix,
                                                             m_face_to_cell_local_face_matrix);

src/mesh/Connectivity.hpp

@@ -8,6 +8,7 @@
 
 #include <ConnectivityMatrix.hpp>
 #include <ConnectivityComputer.hpp>
+#include <SubItemValuePerItem.hpp>
 
 #include <vector>
 #include <unordered_map>
@@ -21,6 +22,8 @@
 #include <tuple>
 #include <algorithm>
 
+#include <IConnectivity.hpp>
+
 template <size_t Dimension>
 class Connectivity;
 
@@ -218,7 +221,8 @@ class ConnectivityFace<3>
 };
 
 template <size_t Dimension>
-class Connectivity
+class Connectivity final
+    : public IConnectivity
 {
 public:
   static constexpr size_t dimension = Dimension;
@@ -231,16 +235,19 @@ public:
   ConnectivityMatrixShort m_cell_to_face_is_reversed_matrix;
 
   ConnectivityMatrix m_face_to_cell_matrix;
-  ConnectivityMatrixShort m_face_to_cell_local_face_matrix;
+  CellValuePerFace<unsigned short> m_face_to_cell_local_face_matrix;
   ConnectivityMatrix m_face_to_node_matrix;
 
   ConnectivityMatrix m_node_to_cell_matrix;
-  ConnectivityMatrixShort m_node_to_cell_local_node_matrix;
+  CellValuePerNode<unsigned short> m_node_to_cell_local_node_matrix;
 
   template <TypeOfItem SubItemType,
             TypeOfItem ItemType>
   const ConnectivityMatrix& itemToItemMatrix() const = delete;
 
+  const ConnectivityMatrix& itemToItemMatrix(const TypeOfItem& item_type_0,
+                                             const TypeOfItem& item_type_1) const final;
+
 private:
   ConnectivityComputer m_connectivity_computer;
 
@@ -344,6 +351,8 @@ private:
     m_connectivity_computer.computeInverseConnectivityMatrix(m_cell_to_node_matrix,
                                                              m_node_to_cell_matrix);
 
+    m_node_to_cell_local_node_matrix = CellValuePerNode<unsigned short>(*this);
+
     m_connectivity_computer.computeLocalChildItemNumberInItem(m_cell_to_node_matrix,
                                                               m_node_to_cell_matrix,
                                                               m_node_to_cell_local_node_matrix);
@@ -361,6 +370,15 @@ private:
 
 using Connectivity3D = Connectivity<3>;
 
+template <>
+template <>
+inline const ConnectivityMatrix&
+Connectivity<3>::itemToItemMatrix<TypeOfItem::cell,
+                                  TypeOfItem::face>() const
+{
+  return m_cell_to_face_matrix;
+}
+
 template <>
 template <>
 inline const ConnectivityMatrix&
@@ -370,8 +388,36 @@ Connectivity<3>::itemToItemMatrix<TypeOfItem::cell,
   return m_cell_to_node_matrix;
 }
 
+template <>
+template <>
+inline const ConnectivityMatrix&
+Connectivity<3>::itemToItemMatrix<TypeOfItem::face,
+                                  TypeOfItem::cell>() const
+{
+  return m_face_to_cell_matrix;
+}
+
+template <>
+template <>
+inline const ConnectivityMatrix&
+Connectivity<3>::itemToItemMatrix<TypeOfItem::node,
+                                  TypeOfItem::cell>() const
+{
+  return m_node_to_cell_matrix;
+}
+
+
 using Connectivity2D = Connectivity<2>;
 
+template <>
+template <>
+inline const ConnectivityMatrix&
+Connectivity<2>::itemToItemMatrix<TypeOfItem::cell,
+                                  TypeOfItem::face>() const
+{
+  return m_cell_to_face_matrix;
+}
+
 template <>
 template <>
 inline const ConnectivityMatrix&
@@ -381,6 +427,24 @@ Connectivity<2>::itemToItemMatrix<TypeOfItem::cell,
   return m_cell_to_node_matrix;
 }
 
+template <>
+template <>
+inline const ConnectivityMatrix&
+Connectivity<2>::itemToItemMatrix<TypeOfItem::face,
+                                  TypeOfItem::cell>() const
+{
+  return m_face_to_cell_matrix;
+}
+
+template <>
+template <>
+inline const ConnectivityMatrix&
+Connectivity<2>::itemToItemMatrix<TypeOfItem::node,
+                                  TypeOfItem::cell>() const
+{
+  return m_node_to_cell_matrix;
+}
+
 using Connectivity1D = Connectivity<1>;
 
 template <>
@@ -392,4 +456,82 @@ Connectivity<1>::itemToItemMatrix<TypeOfItem::cell,
   return m_cell_to_node_matrix;
 }
 
+template <>
+template <>
+inline const ConnectivityMatrix&
+Connectivity<1>::itemToItemMatrix<TypeOfItem::cell,
+                                  TypeOfItem::face>() const
+{
+#warning in 1d, faces and node are the same
+  return m_cell_to_face_matrix;
+}
+
+template <>
+template <>
+inline const ConnectivityMatrix&
+Connectivity<1>::itemToItemMatrix<TypeOfItem::face,
+                                  TypeOfItem::cell>() const
+{
+#warning in 1d, faces and node are the same
+  return m_face_to_cell_matrix;
+}
+
+template <>
+template <>
+inline const ConnectivityMatrix&
+Connectivity<1>::itemToItemMatrix<TypeOfItem::node,
+                                  TypeOfItem::cell>() const
+{
+#warning in 1d, faces and node are the same
+  return m_node_to_cell_matrix;
+}
+
+template <size_t Dimension>
+const ConnectivityMatrix&
+Connectivity<Dimension>::
+itemToItemMatrix(const TypeOfItem& item_type_0,
+                 const TypeOfItem& item_type_1) const
+{
+  switch (item_type_0) {
+    case TypeOfItem::cell: {
+      switch (item_type_1) {
+        case TypeOfItem::node: {
+          return itemToItemMatrix<TypeOfItem::cell, TypeOfItem::node>();
+        }
+        default: {
+          std::cerr << __FILE__ << ":" << __LINE__ << ": NIY " << int(item_type_1) << "\n";
+          std::exit(1);
+        }
+      }
+    }
+    case TypeOfItem::face: {
+      switch (item_type_1) {
+        case TypeOfItem::cell: {
+          return itemToItemMatrix<TypeOfItem::face, TypeOfItem::cell>();
+        }
+        default: {
+          std::cerr << __FILE__ << ":" << __LINE__ << ": NIY " << int(item_type_1) << "\n";
+          std::exit(1);
+        }
+      }
+    }
+    case TypeOfItem::node: {
+      switch (item_type_1) {
+        case TypeOfItem::cell: {
+          return itemToItemMatrix<TypeOfItem::node, TypeOfItem::cell>();
+        }
+        default: {
+          std::cerr << __FILE__ << ":" << __LINE__ << ": NIY " << int(item_type_1) << "\n";
+          std::exit(1);
+        }
+      }
+    }
+    default: {
+      std::cerr << __FILE__ << ":" << __LINE__ << ": NIY " << int(item_type_0) << "\n";
+      std::exit(1);
+    }
+  }
+}
+
+
 #endif // CONNECTIVITY_HPP

src/mesh/ConnectivityComputer.cpp

@@ -36,23 +36,43 @@ computeInverseConnectivityMatrix(const ConnectivityMatrix& item_to_child_item_ma
 void ConnectivityComputer::
 computeLocalChildItemNumberInItem(const ConnectivityMatrix& item_to_child_items_matrix,
                                   const ConnectivityMatrix& child_item_to_items_matrix,
-                                  ConnectivityMatrixShort&  child_item_number_in_item_matrix) const
+                                  CellValuePerNode<unsigned short>& child_item_number_in_item_matrix) const
 {
-  std::vector<std::vector<unsigned int>> child_item_number_in_item_vector(child_item_to_items_matrix.numRows());
   for (unsigned int r=0; r<child_item_to_items_matrix.numRows(); ++r) {
     const auto& child_item_to_items = child_item_to_items_matrix.rowConst(r);
-    child_item_number_in_item_vector[r].resize(child_item_to_items.length);
     for (unsigned short J=0; J<child_item_to_items.length; ++J) {
       const unsigned int j = child_item_to_items(J);
       const auto& item_to_child_items = item_to_child_items_matrix.rowConst(j);
 
       for (unsigned int R=0; R<item_to_child_items.length; ++R) {
         if (item_to_child_items(R) == r) {
-          child_item_number_in_item_vector[r][J] = R;
+          child_item_number_in_item_matrix(r,J) = R;
+          break;
+        }
+      }
+    }
+  }
+}
+
+
+
+void ConnectivityComputer::
+computeLocalChildItemNumberInItem(const ConnectivityMatrix& item_to_child_items_matrix,
+                                  const ConnectivityMatrix& child_item_to_items_matrix,
+                                  CellValuePerFace<unsigned short>& child_item_number_in_item_matrix) const
+{
+  for (unsigned int r=0; r<child_item_to_items_matrix.numRows(); ++r) {
+    const auto& child_item_to_items = child_item_to_items_matrix.rowConst(r);
+    for (unsigned short J=0; J<child_item_to_items.length; ++J) {
+      const unsigned int j = child_item_to_items(J);
+      const auto& item_to_child_items = item_to_child_items_matrix.rowConst(j);
+
+      for (unsigned int R=0; R<item_to_child_items.length; ++R) {
+        if (item_to_child_items(R) == r) {
+          child_item_number_in_item_matrix(r,J) = R;
           break;
         }
       }
     }
   }
-  child_item_number_in_item_matrix = child_item_number_in_item_vector;
 }

src/mesh/ConnectivityComputer.hpp

@@ -2,6 +2,7 @@
 #define CONNECTIVITY_COMPUTER_HPP
 
 #include <ConnectivityMatrix.hpp>
+#include <SubItemValuePerItem.hpp>
 
 struct ConnectivityComputer
 {
@@ -10,7 +11,11 @@ struct ConnectivityComputer
 
   void computeLocalChildItemNumberInItem(const ConnectivityMatrix& item_to_child_items_matrix,
                                          const ConnectivityMatrix& child_item_to_items_matrix,
-                                         ConnectivityMatrixShort&  child_item_number_in_item_matrix) const;
+                                         CellValuePerNode<unsigned short>&  child_item_number_in_item_matrix) const;
+
+  void computeLocalChildItemNumberInItem(const ConnectivityMatrix& item_to_child_items_matrix,
+                                         const ConnectivityMatrix& child_item_to_items_matrix,
+                                         CellValuePerFace<unsigned short>&  child_item_number_in_item_matrix) const;
 };
 
 #endif // CONNECTIVITY_COMPUTER_HPP

src/mesh/IConnectivity.hpp

+#ifndef ICONNECTIVITY_HPP
+#define ICONNECTIVITY_HPP
+
+#include <TypeOfItem.hpp>
+#include <ConnectivityMatrix.hpp>
+
+struct IConnectivity
+{
+  virtual const ConnectivityMatrix&
+  itemToItemMatrix(const TypeOfItem& item_type_0,
+                   const TypeOfItem& item_type_1) const = 0;
+
+  IConnectivity() = default;
+  IConnectivity(const IConnectivity&) = delete;
+  virtual ~IConnectivity() = default;
+};
+
+#endif // ICONNECTIVITY_HPP

src/mesh/SubItemValuePerItem.hpp

@@ -3,8 +3,10 @@
 
 #include <Kokkos_StaticCrsGraph.hpp>
 #include <TypeOfItem.hpp>
+#include <PastisAssert.hpp>
 
 #include <ConnectivityMatrix.hpp>
+#include <IConnectivity.hpp>
 
 template <typename DataType,
           TypeOfItem SubItemType,
@@ -158,11 +160,10 @@ class SubItemValuePerItem
 
   SubItemValuePerItem() = default;
 
-  template <typename ConnectivityType>
-  SubItemValuePerItem(const ConnectivityType& connectivity)
+  SubItemValuePerItem(const IConnectivity& connectivity)
   {
     ConnectivityMatrix connectivity_matrix
-        = connectivity.template itemToItemMatrix<ItemType,SubItemType>();
+        = connectivity.itemToItemMatrix(ItemType,SubItemType);
     m_host_row_map = connectivity_matrix.rowsMap();
     m_values = Kokkos::View<DataType*>("values", connectivity_matrix.numEntries());
   }
@@ -173,4 +174,10 @@ class SubItemValuePerItem
 template <typename DataType>
 using NodeValuePerCell = SubItemValuePerItem<DataType, TypeOfItem::node, TypeOfItem::cell>;
 
+template <typename DataType>
+using CellValuePerNode = SubItemValuePerItem<DataType, TypeOfItem::cell, TypeOfItem::node>;
+
+template <typename DataType>
+using CellValuePerFace = SubItemValuePerItem<DataType, TypeOfItem::cell, TypeOfItem::face>;
+
 #endif // SUBITEM_VALUE_PER_ITEM_HPP

src/scheme/AcousticSolver.hpp

@@ -101,10 +101,10 @@ class AcousticSolver
     Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
         Rdd sum = zero;
         const auto& node_to_cell = m_connectivity.m_node_to_cell_matrix.rowConst(r);
-        const auto& node_to_cell_local_node = m_connectivity.m_node_to_cell_local_node_matrix.rowConst(r);
+        const auto& node_to_cell_local_node = m_connectivity.m_node_to_cell_local_node_matrix.itemValues(r);
         for (size_t j=0; j<node_to_cell.length; ++j) {
           const unsigned int J = node_to_cell(j);
-          const unsigned int R = node_to_cell_local_node(j);
+          const unsigned int R = node_to_cell_local_node[j];
           sum += Ajr(J,R);
         }
         m_Ar(r) = sum;
@@ -123,10 +123,10 @@ class AcousticSolver
         Rd& br = m_br(r);
         br = zero;
         const auto& node_to_cell = m_connectivity.m_node_to_cell_matrix.rowConst(r);
-        const auto& node_to_cell_local_node = m_connectivity.m_node_to_cell_local_node_matrix.rowConst(r);
+        const auto& node_to_cell_local_node = m_connectivity.m_node_to_cell_local_node_matrix.itemValues(r);
         for (size_t j=0; j<node_to_cell.length; ++j) {
           const unsigned int J = node_to_cell(j);
-          const unsigned int R = node_to_cell_local_node(j);
+          const unsigned int R = node_to_cell_local_node[j];
           br += Ajr(J,R)*uj(J) + pj(J)*Cjr(J,R);
         }
       });