Compare revisions

Stéphane Del Pino · Stéphane Del Pino · Stéphane Del Pino · Stéphane Del Pino · Stéphane Del Pino · Stéphane Del Pino
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -271,8 +271,9 @@ add_executable(
 # Libraries
 target_link_libraries(
  pastis
+  PastisMesh
+  PastisUtils
  kokkos
  ${PARMETIS_LIBRARIES}
  ${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES}
-  PastisMesh
+)
-  PastisUtils)
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -296,7 +296,7 @@ template <size_t N, typename T>
 PASTIS_INLINE
 constexpr T det(const TinyMatrix<N,T>& A)
 {
-  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
  static_assert(std::is_floating_point<T>::value, "determinent for arbitrary dimension N is defined for floating point types only");
  TinyMatrix<N,T> M = A;
@@ -340,7 +340,7 @@ template <typename T>
 PASTIS_INLINE
 constexpr T det(const TinyMatrix<1,T>& A)
 {
-  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
  return A(0,0);
 }
@@ -348,7 +348,7 @@ template <typename T>
 PASTIS_INLINE
 constexpr T det(const TinyMatrix<2,T>& A)
 {
-  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
  return A(0,0)*A(1,1)-A(1,0)*A(0,1);
 }
@@ -356,7 +356,7 @@ template <typename T>
 PASTIS_INLINE
 constexpr T det(const TinyMatrix<3,T>& A)
 {
-  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
  return
    A(0,0)*(A(1,1)*A(2,2)-A(2,1)*A(1,2))
    -A(1,0)*(A(0,1)*A(2,2)-A(2,1)*A(0,2))
@@ -399,7 +399,7 @@ template <typename T>
 PASTIS_INLINE
 constexpr TinyMatrix<1,T> inverse(const TinyMatrix<1,T>& A)
 {
-  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
  static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
  TinyMatrix<1,T> A_1(1./A(0,0));
@@ -412,7 +412,7 @@ constexpr T cofactor(const TinyMatrix<N,T>& A,
                     const size_t& i,
                     const size_t& j)
 {
-  static_assert(std::is_arithmetic<T>::value, "cofactor is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "cofactor is not defined for non-arithmetic types");
  const T sign = ((i+j)%2) ? -1 : 1;
  return sign * det(getMinor(A, i, j));
@@ -422,7 +422,7 @@ template <typename T>
 PASTIS_INLINE
 constexpr TinyMatrix<2,T> inverse(const TinyMatrix<2,T>& A)
 {
-  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
  static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
  const T determinent = det(A);
@@ -438,7 +438,7 @@ template <typename T>
 PASTIS_INLINE
 constexpr TinyMatrix<3,T> inverse(const TinyMatrix<3,T>& A)
 {
-  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
  static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
  const T determinent = det(A);

--- a/src/main.cpp
+++ b/src/main.cpp
@@ -74,9 +74,9 @@ int main(int argc, char *argv[])
              case BoundaryConditionDescriptor::Type::symmetry: {
                const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
                    = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-                for (size_t i_ref_node_list=0; i_ref_node_list<mesh.connectivity().numberOfRefNodeList();
+                for (size_t i_ref_node_list=0; i_ref_node_list<mesh.connectivity().numberOfRefItemList<ItemType::node>();
                     ++i_ref_node_list) {
-                  const RefNodeList& ref_node_list = mesh.connectivity().refNodeList(i_ref_node_list);
+                  const RefNodeList& ref_node_list = mesh.connectivity().refItemList<ItemType::node>(i_ref_node_list);
                  const RefId& ref = ref_node_list.refId();
                  if (ref == sym_bc_descriptor.boundaryDescriptor()) {
                    SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc
@@ -191,9 +191,9 @@ int main(int argc, char *argv[])
              case BoundaryConditionDescriptor::Type::symmetry: {
                const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
                    = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-                for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
+                for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefItemList<ItemType::face>();
                     ++i_ref_face_list) {
-                  const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
+                  const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
                  const RefId& ref = ref_face_list.refId();
                  if (ref == sym_bc_descriptor.boundaryDescriptor()) {
                    SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc
@@ -295,9 +295,9 @@ int main(int argc, char *argv[])
              case BoundaryConditionDescriptor::Type::symmetry: {
                const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
                    = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-                for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
+                for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefItemList<ItemType::face>();
                     ++i_ref_face_list) {
-                  const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
+                  const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
                  const RefId& ref = ref_face_list.refId();
                  if (ref == sym_bc_descriptor.boundaryDescriptor()) {
                    SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc

--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -3,6 +3,8 @@
 #include <Messenger.hpp>
+#include <ConnectivityDescriptor.hpp>
 template<size_t Dimension>
 Connectivity<Dimension>::Connectivity() {}
@@ -10,8 +12,7 @@ template<size_t Dimension>
 void Connectivity<Dimension>::
 _buildFrom(const ConnectivityDescriptor& descriptor)
 {
-#warning All of these should be checked by ConnectivityDescriptor
+  Assert(descriptor.cell_to_node_vector.size() == descriptor.cell_type_vector.size());
-  Assert(descriptor.cell_by_node_vector.size() == descriptor.cell_type_vector.size());
  Assert(descriptor.cell_number_vector.size() == descriptor.cell_type_vector.size());
  if constexpr (Dimension>1) {
    Assert(descriptor.cell_to_face_vector.size() == descriptor.cell_type_vector.size());
@@ -21,7 +22,7 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
  auto& cell_to_node_matrix
      = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::node)];
-  cell_to_node_matrix = descriptor.cell_by_node_vector;
+  cell_to_node_matrix = descriptor.cell_to_node_vector;
  {
    WeakCellValue<CellType> cell_type(*this);
@@ -45,8 +46,6 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
  {
    WeakCellValue<int> cell_global_index(*this);
-#warning index must start accounting number of global indices of other procs
-#warning must take care of ghost cells
    int first_index = 0;
    parallel_for(this->numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
      cell_global_index[j] = first_index+j;
@@ -54,16 +53,6 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
    m_cell_global_index = cell_global_index;
  }
-  {
-    WeakCellValue<double> inv_cell_nb_nodes(*this);
-    parallel_for(this->numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
-        const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
-        inv_cell_nb_nodes[j] = 1./cell_nodes.length;
-      });
-    m_inv_cell_nb_nodes = inv_cell_nb_nodes;
-  }
  {
    WeakCellValue<int> cell_owner(*this);
    cell_owner = convert_to_array(descriptor.cell_owner_vector);
@@ -94,14 +83,13 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
    m_node_is_owned = node_is_owned;
  }
+  m_ref_node_list_vector = descriptor.template refItemListVector<ItemType::node>();
+  m_ref_cell_list_vector = descriptor.template refItemListVector<ItemType::cell>();
  if constexpr (Dimension>1) {
-    auto& face_to_node_matrix
+    m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)] = descriptor.face_to_node_vector;
-        = m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)];
-    face_to_node_matrix = descriptor.face_to_node_vector;
-    auto& cell_to_face_matrix
+    m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)] = descriptor.cell_to_face_vector;
-        = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)];
-    cell_to_face_matrix = descriptor.cell_to_face_vector;
    {
      FaceValuePerCell<bool> cell_face_is_reversed(*this);
@@ -111,14 +99,15 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
          cell_face_is_reversed(j,lj) = face_cells_vector[lj];
        }
      }
      m_cell_face_is_reversed = cell_face_is_reversed;
    }
    {
      WeakFaceValue<int> face_number(*this);
      face_number = convert_to_array(descriptor.face_number_vector);
      m_face_number = face_number;
    }
    {
      WeakFaceValue<int> face_owner(*this);
      face_owner = convert_to_array(descriptor.face_owner_vector);
@@ -134,8 +123,49 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
      m_face_is_owned = face_is_owned;
    }
+    m_ref_face_list_vector = descriptor.template refItemListVector<ItemType::face>();
+    if constexpr (Dimension>2) {
+      m_item_to_item_matrix[itemTId(ItemType::edge)][itemTId(ItemType::node)] = descriptor.edge_to_node_vector;
-    m_ref_face_list = descriptor.refFaceList();
+      m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::edge)] = descriptor.face_to_edge_vector;
+      m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::edge)] = descriptor.cell_to_edge_vector;
+      {
+        EdgeValuePerFace<bool> face_edge_is_reversed(*this);
+        for (FaceId l=0; l<descriptor.face_edge_is_reversed_vector.size(); ++l) {
+          const auto& edge_faces_vector = descriptor.face_edge_is_reversed_vector[l];
+          for (unsigned short el=0; el<edge_faces_vector.size(); ++el) {
+            face_edge_is_reversed(l,el) = edge_faces_vector[el];
+          }
+        }
+        m_face_edge_is_reversed = face_edge_is_reversed;
+      }
+      {
+        WeakEdgeValue<int> edge_number(*this);
+        edge_number = convert_to_array(descriptor.edge_number_vector);
+        m_edge_number = edge_number;
+      }
+      {
+        WeakEdgeValue<int> edge_owner(*this);
+        edge_owner = convert_to_array(descriptor.edge_owner_vector);
+        m_edge_owner = edge_owner;
+      }
+      {
+        const int rank = parallel::rank();
+        WeakEdgeValue<bool> edge_is_owned(*this);
+        parallel_for(this->numberOfEdges(), PASTIS_LAMBDA(const EdgeId& e) {
+            edge_is_owned[e] = (m_edge_owner[e] == rank);
+          });
+        m_edge_is_owned = edge_is_owned;
+      }
+      m_ref_edge_list_vector = descriptor.template refItemListVector<ItemType::edge>();
+    }
  }
 }

--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -7,6 +7,8 @@
 #include <PastisOStream.hpp>
 #include <PastisUtils.hpp>
+#include <PastisTraits.hpp>
 #include <TinyVector.hpp>
 #include <ItemValue.hpp>
@@ -28,58 +30,14 @@
 #include <RefId.hpp>
 #include <ItemType.hpp>
-#include <RefNodeList.hpp>
+#include <RefItemList.hpp>
-#include <RefFaceList.hpp>
 #include <SynchronizerManager.hpp>
-#include <tuple>
-#include <algorithm>
 #include <set>
-class ConnectivityDescriptor
+class ConnectivityDescriptor;
-{
- private:
-  std::vector<RefFaceList> m_ref_face_list;
- public:
-  std::vector<std::vector<unsigned int>> cell_by_node_vector;
-#warning should be set in 2d
-  std::vector<std::vector<unsigned int>> cell_to_face_vector;
-  std::vector<std::vector<bool>> cell_face_is_reversed_vector;
-  // std::vector<std::vector<unsigned int>> face_to_cell_vector;
-  std::vector<std::vector<unsigned int>> face_to_node_vector;
-  std::vector<CellType> cell_type_vector;
-  std::vector<int> cell_number_vector;
-  std::vector<int> face_number_vector;
-  std::vector<int> node_number_vector;
-  std::vector<int> cell_owner_vector;
-  std::vector<int> face_owner_vector;
-  std::vector<int> node_owner_vector;
-  const std::vector<RefFaceList>& refFaceList() const
-  {
-    return m_ref_face_list;
-  }
-  void addRefFaceList(const RefFaceList& ref_face_list)
-  {
-    m_ref_face_list.push_back(ref_face_list);
-  }
-  ConnectivityDescriptor& operator=(const ConnectivityDescriptor&) = delete;
-  ConnectivityDescriptor& operator=(ConnectivityDescriptor&&) = delete;
-  ConnectivityDescriptor() = default;
-  ConnectivityDescriptor(const ConnectivityDescriptor&) = default;
-  ConnectivityDescriptor(ConnectivityDescriptor&&) = delete;
-  ~ConnectivityDescriptor() = default;
-};
 template <size_t Dim>
 class Connectivity final
@@ -105,12 +63,10 @@ class Connectivity final
  ConnectivityMatrix m_item_to_item_matrix[Dimension+1][Dimension+1];
  WeakCellValue<const CellType> m_cell_type;
-#warning is m_cell_global_index really necessary? should it be computed on demand instead?
  WeakCellValue<const int> m_cell_global_index;
  WeakCellValue<const int> m_cell_number;
  WeakFaceValue<const int> m_face_number;
-#warning check that m_edge_number is filled correctly
  WeakEdgeValue<const int> m_edge_number;
  WeakNodeValue<const int> m_node_number;
@@ -127,6 +83,7 @@ class Connectivity final
  WeakNodeValue<const bool> m_node_is_owned;
  WeakFaceValuePerCell<const bool> m_cell_face_is_reversed;
+  WeakEdgeValuePerFace<const bool> m_face_edge_is_reversed;
  WeakNodeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_nodes;
  WeakEdgeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_edges;
@@ -146,8 +103,10 @@ class Connectivity final
  ConnectivityComputer m_connectivity_computer;
-  std::vector<RefFaceList> m_ref_face_list;
+  std::vector<RefCellList> m_ref_cell_list_vector;
-  std::vector<RefNodeList> m_ref_node_list;
+  std::vector<RefFaceList> m_ref_face_list_vector;
+  std::vector<RefEdgeList> m_ref_edge_list_vector;
+  std::vector<RefNodeList> m_ref_node_list_vector;
  WeakCellValue<const double> m_inv_cell_nb_nodes;
@@ -228,8 +187,7 @@ class Connectivity final
    } else if constexpr(item_type == ItemType::node) {
      return m_node_number;
    } else {
-      static_assert(item_type == ItemType::cell, "unknown ItemType");
+      static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
-      return m_cell_number;
    }
  }
@@ -272,8 +230,7 @@ class Connectivity final
    } else if constexpr(item_type == ItemType::node) {
      return m_node_owner;
    } else {
-      static_assert(item_type == ItemType::cell, "unknown ItemType");
+      static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
-      return m_cell_owner;
    }
  }
@@ -316,8 +273,7 @@ class Connectivity final
    } else if constexpr(item_type == ItemType::node) {
      return m_node_is_owned;
    } else {
-      static_assert(item_type == ItemType::cell, "unknown ItemType");
+      static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
-      return m_cell_is_owned;
    }
  }
@@ -418,6 +374,13 @@ class Connectivity final
    return m_cell_face_is_reversed;
  }
+  PASTIS_INLINE
+  const auto& faceEdgeIsReversed() const
+  {
+    static_assert(Dimension>2, "reversed edges makes no sense in dimension 1 or 2");
+    return m_face_edge_is_reversed;
+  }
  PASTIS_INLINE
  const auto& cellLocalNumbersInTheirNodes() const
  {
@@ -512,29 +475,53 @@ class Connectivity final
    return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_faces);
  }
-  size_t numberOfRefFaceList() const
+  template <ItemType item_type>
+  size_t numberOfRefItemList() const
  {
-    return m_ref_face_list.size();
+    if constexpr (item_type == ItemType::cell) {
+      return m_ref_cell_list_vector.size();
+    } else if constexpr (item_type == ItemType::face) {
+      return m_ref_face_list_vector.size();
+    } else if constexpr (item_type == ItemType::edge) {
+      return m_ref_edge_list_vector.size();
+    } else if constexpr (item_type == ItemType::node) {
+      return m_ref_node_list_vector.size();
+    } else {
+      static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
    }
-  const RefFaceList& refFaceList(const size_t& i) const
-  {
-    return m_ref_face_list[i];
  }
-  void addRefNodeList(const RefNodeList& ref_node_list)
+  template <ItemType item_type>
+  const RefItemList<item_type>& refItemList(const size_t& i) const
  {
-    m_ref_node_list.push_back(ref_node_list);
+    if constexpr (item_type == ItemType::cell) {
+      return m_ref_cell_list_vector[i];
+    } else if constexpr (item_type == ItemType::face) {
+      return m_ref_face_list_vector[i];
+    } else if constexpr (item_type == ItemType::edge) {
+      return m_ref_edge_list_vector[i];
+    } else if constexpr (item_type == ItemType::node) {
+      return m_ref_node_list_vector[i];
+    } else {
+      static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
    }
-  size_t numberOfRefNodeList() const
-  {
-    return m_ref_node_list.size();
  }
-  const RefNodeList& refNodeList(const size_t& i) const
+  template <ItemType item_type>
+  void addRefItemList(const RefItemList<item_type>& ref_item_list)
  {
-    return m_ref_node_list[i];
+    if constexpr (item_type == ItemType::cell) {
+      m_ref_cell_list_vector.push_back(ref_item_list);
+    } else if constexpr (item_type == ItemType::face) {
+      m_ref_face_list_vector.push_back(ref_item_list);
+    } else if constexpr (item_type == ItemType::edge) {
+      m_ref_edge_list_vector.push_back(ref_item_list);
+    } else if constexpr (item_type == ItemType::node) {
+      m_ref_node_list_vector.push_back(ref_item_list);
+    } else {
+      static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
+    }
  }
  PASTIS_INLINE
@@ -624,7 +611,18 @@ class Connectivity final
  CellValue<const double> invCellNbNodes() const
  {
-#warning add calculation on demand when variables will be defined
+    if (not m_inv_cell_nb_nodes.isBuilt()) {
+      const auto& cell_to_node_matrix
+          = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::node)];
+      WeakCellValue<double> inv_cell_nb_nodes(*this);
+      parallel_for(this->numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+          const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+          inv_cell_nb_nodes[j] = 1./cell_nodes.length;
+        });
+      const_cast<WeakCellValue<const double>&>(m_inv_cell_nb_nodes) = inv_cell_nb_nodes;
+    }
    return m_inv_cell_nb_nodes;
  }

--- a/src/mesh/ConnectivityDescriptor.hpp
+++ b/src/mesh/ConnectivityDescriptor.hpp
+#ifndef CONNECTIVITY_DESCRIPTOR_HPP
+#define CONNECTIVITY_DESCRIPTOR_HPP
+#include <RefItemList.hpp>
+#include <PastisTraits.hpp>
+#include <ItemOfItemType.hpp>
+#include <vector>
+class ConnectivityDescriptor
+{
+ private:
+  std::vector<RefCellList> m_ref_cell_list_vector;
+  std::vector<RefFaceList> m_ref_face_list_vector;
+  std::vector<RefEdgeList> m_ref_edge_list_vector;
+  std::vector<RefNodeList> m_ref_node_list_vector;
+ public:
+  std::vector<std::vector<unsigned int>> cell_to_node_vector;
+  std::vector<std::vector<unsigned int>> cell_to_face_vector;
+  std::vector<std::vector<unsigned int>> cell_to_edge_vector;
+  std::vector<std::vector<unsigned int>> face_to_node_vector;
+  std::vector<std::vector<unsigned int>> face_to_edge_vector;
+  std::vector<std::vector<unsigned int>> edge_to_node_vector;
+  template <typename ItemOfItemT>
+  auto& itemOfItemVector()
+  {
+    if constexpr (std::is_same_v<ItemOfItemT,NodeOfCell>) {
+      return cell_to_node_vector;
+    } else if constexpr (std::is_same_v<ItemOfItemT,FaceOfCell>) {
+      return cell_to_face_vector;
+    } else if constexpr (std::is_same_v<ItemOfItemT,EdgeOfCell>) {
+      return cell_to_edge_vector;
+    } else if constexpr (std::is_same_v<ItemOfItemT,EdgeOfFace>) {
+      return face_to_edge_vector;
+    } else if constexpr (std::is_same_v<ItemOfItemT,NodeOfFace>) {
+      return face_to_node_vector;
+    } else if constexpr (std::is_same_v<ItemOfItemT,NodeOfEdge>) {
+      return edge_to_node_vector;
+    } else {
+      static_assert(is_false_v<ItemOfItemT>, "Unexpected item of item type");
+    }
+  }
+  std::vector<Array<bool>> cell_face_is_reversed_vector;
+  std::vector<Array<bool>> face_edge_is_reversed_vector;
+  std::vector<CellType> cell_type_vector;
+  std::vector<int> cell_number_vector;
+  std::vector<int> face_number_vector;
+  std::vector<int> edge_number_vector;
+  std::vector<int> node_number_vector;
+  template <ItemType item_type>
+  const std::vector<int>& itemNumberVector() const
+  {
+    if constexpr (item_type == ItemType::cell) {
+      return cell_number_vector;
+    } else if constexpr (item_type == ItemType::face) {
+      return face_number_vector;
+    } else if constexpr (item_type == ItemType::edge) {
+      return edge_number_vector;
+    } else if constexpr (item_type == ItemType::node) {
+      return node_number_vector;
+    } else {
+      static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
+    }
+  }
+  std::vector<int> cell_owner_vector;
+  std::vector<int> face_owner_vector;
+  std::vector<int> edge_owner_vector;
+  std::vector<int> node_owner_vector;
+  template <ItemType item_type>
+  const std::vector<RefItemList<item_type>>& refItemListVector() const
+  {
+    if constexpr (item_type == ItemType::cell) {
+      return m_ref_cell_list_vector;
+    } else if constexpr (item_type == ItemType::face) {
+      return m_ref_face_list_vector;
+    } else if constexpr (item_type == ItemType::edge) {
+      return m_ref_edge_list_vector;
+    } else if constexpr (item_type == ItemType::node) {
+      return m_ref_node_list_vector;
+    } else {
+      static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
+    }
+  }
+  template <ItemType item_type>
+  void addRefItemList(const RefItemList<item_type>& ref_item_list)
+  {
+    if constexpr (item_type == ItemType::cell) {
+      m_ref_cell_list_vector.push_back(ref_item_list);
+    } else if constexpr (item_type == ItemType::face) {
+      m_ref_face_list_vector.push_back(ref_item_list);
+    } else if constexpr (item_type == ItemType::edge) {
+      m_ref_edge_list_vector.push_back(ref_item_list);
+    } else if constexpr (item_type == ItemType::node) {
+      m_ref_node_list_vector.push_back(ref_item_list);
+    } else {
+      static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
+    }
+  }
+  ConnectivityDescriptor& operator=(const ConnectivityDescriptor&) = delete;
+  ConnectivityDescriptor& operator=(ConnectivityDescriptor&&) = delete;
+  ConnectivityDescriptor() = default;
+  ConnectivityDescriptor(const ConnectivityDescriptor&) = default;
+  ConnectivityDescriptor(ConnectivityDescriptor&&) = delete;
+  ~ConnectivityDescriptor() = default;
+};
+#endif // CONNECTIVITY_DESCRIPTOR_HPP
--- a/src/mesh/ConnectivityDispatcher.cpp
+++ b/src/mesh/ConnectivityDispatcher.cpp
--- a/src/mesh/ConnectivityDispatcher.hpp
+++ b/src/mesh/ConnectivityDispatcher.hpp
@@ -6,6 +6,7 @@
 #include <ItemValueUtils.hpp>
 #include <unordered_map>
+#include <ConnectivityDescriptor.hpp>
 template <int Dimension>
 class ConnectivityDispatcher
@@ -23,8 +24,8 @@ class ConnectivityDispatcher
  {
    using ItemId = ItemIdT<item_type>;
    ItemValue<const int, item_type> m_new_owner;
+    Array<const unsigned int> m_list_to_send_size_by_proc;
    std::vector<Array<const ItemId>> m_list_to_send_by_proc;
-#warning is m_list_to_recv_size_by_proc really necessary?
    Array<const unsigned int> m_list_to_recv_size_by_proc;
    std::unordered_map<int, int> m_number_to_id_map;
    std::vector<Array<const ItemId>> m_recv_id_correspondance_by_proc;
@@ -65,36 +66,138 @@ class ConnectivityDispatcher
    }
  }
+  template <typename ItemToItem>
+  struct DispatchedItemOfItemInfo
+  {
+    std::vector<Array<const int>> m_number_of_sub_item_per_item_to_recv_by_proc;
+    std::vector<Array<const int>> m_sub_item_numbers_to_recv_by_proc;
+  };
+  DispatchedItemOfItemInfo<NodeOfCell> m_dispatched_node_of_cell_info;
+  DispatchedItemOfItemInfo<EdgeOfCell> m_dispatched_edge_of_cell_info;
+  DispatchedItemOfItemInfo<FaceOfCell> m_dispatched_face_of_cell_info;
+  DispatchedItemOfItemInfo<NodeOfEdge> m_dispatched_node_of_edge_info;
+  DispatchedItemOfItemInfo<FaceOfEdge> m_dispatched_face_of_edge_info;
+  DispatchedItemOfItemInfo<CellOfEdge> m_dispatched_cell_of_edge_info;
+  DispatchedItemOfItemInfo<NodeOfFace> m_dispatched_node_of_face_info;
+  DispatchedItemOfItemInfo<EdgeOfFace> m_dispatched_edge_of_face_info;
+  DispatchedItemOfItemInfo<CellOfFace> m_dispatched_cell_of_face_info;
+  DispatchedItemOfItemInfo<EdgeOfNode> m_dispatched_edge_of_node_info;
+  DispatchedItemOfItemInfo<FaceOfNode> m_dispatched_face_of_node_info;
+  DispatchedItemOfItemInfo<CellOfNode> m_dispatched_cell_of_node_info;
+  template <typename ItemOfItem>
+  PASTIS_INLINE
+  DispatchedItemOfItemInfo<ItemOfItem>& _dispatchedInfo()
+  {
+    if constexpr (std::is_same_v<NodeOfCell, ItemOfItem>) {
+      return m_dispatched_node_of_cell_info;
+    } else if constexpr (std::is_same_v<EdgeOfCell, ItemOfItem>) {
+      return m_dispatched_edge_of_cell_info;
+    } else if constexpr (std::is_same_v<FaceOfCell, ItemOfItem>) {
+      return m_dispatched_face_of_cell_info;
+    } else if constexpr (std::is_same_v<NodeOfEdge, ItemOfItem>) {
+      return m_dispatched_node_of_edge_info;
+    } else if constexpr (std::is_same_v<FaceOfEdge, ItemOfItem>) {
+      return m_dispatched_face_of_edge_info;
+    } else if constexpr (std::is_same_v<CellOfEdge, ItemOfItem>) {
+      return m_dispatched_cell_of_edge_info;
+    } else if constexpr (std::is_same_v<NodeOfFace, ItemOfItem>) {
+      return m_dispatched_node_of_face_info;
+    } else if constexpr (std::is_same_v<EdgeOfFace, ItemOfItem>) {
+      return m_dispatched_edge_of_face_info;
+    } else if constexpr (std::is_same_v<CellOfFace, ItemOfItem>) {
+      return m_dispatched_cell_of_face_info;
+    } else if constexpr (std::is_same_v<EdgeOfNode, ItemOfItem>) {
+      return m_dispatched_edge_of_node_info;
+    } else if constexpr (std::is_same_v<FaceOfNode, ItemOfItem>) {
+      return m_dispatched_face_of_node_info;
+    } else if constexpr (std::is_same_v<CellOfNode, ItemOfItem>) {
+      return m_dispatched_cell_of_node_info;
+    } else {
+      static_assert(is_false_v<ItemOfItem>, "Unexpected ItemOfItem type");
+    }
+  }
+  template <typename ItemOfItem>
+  PASTIS_INLINE
+  const DispatchedItemOfItemInfo<ItemOfItem>& _dispatchedInfo() const
+  {
+    if constexpr (std::is_same_v<NodeOfCell, ItemOfItem>) {
+      return m_dispatched_node_of_cell_info;
+    } else if constexpr (std::is_same_v<EdgeOfCell, ItemOfItem>) {
+      return m_dispatched_edge_of_cell_info;
+    } else if constexpr (std::is_same_v<FaceOfCell, ItemOfItem>) {
+      return m_dispatched_face_of_cell_info;
+    } else if constexpr (std::is_same_v<NodeOfEdge, ItemOfItem>) {
+      return m_dispatched_node_of_edge_info;
+    } else if constexpr (std::is_same_v<FaceOfEdge, ItemOfItem>) {
+      return m_dispatched_face_of_edge_info;
+    } else if constexpr (std::is_same_v<CellOfEdge, ItemOfItem>) {
+      return m_dispatched_cell_of_edge_info;
+    } else if constexpr (std::is_same_v<NodeOfFace, ItemOfItem>) {
+      return m_dispatched_node_of_face_info;
+    } else if constexpr (std::is_same_v<EdgeOfFace, ItemOfItem>) {
+      return m_dispatched_edge_of_face_info;
+    } else if constexpr (std::is_same_v<CellOfFace, ItemOfItem>) {
+      return m_dispatched_cell_of_face_info;
+    } else if constexpr (std::is_same_v<EdgeOfNode, ItemOfItem>) {
+      return m_dispatched_edge_of_node_info;
+    } else if constexpr (std::is_same_v<FaceOfNode, ItemOfItem>) {
+      return m_dispatched_face_of_node_info;
+    } else if constexpr (std::is_same_v<CellOfNode, ItemOfItem>) {
+      return m_dispatched_cell_of_node_info;
+    } else {
+      static_assert(is_false_v<ItemOfItem>, "Unexpected ItemOfItem type");
+    }
+  }
  template <ItemType item_type>
  void _buildNewOwner();
  template <ItemType item_type>
  void _buildItemListToSend();
-  Array<const unsigned int> _buildNbCellToSend();
  void _buildCellNumberIdMap();
  template <typename ItemOfItemT>
-  void _buildSubItemNumberToIdMap(const std::vector<Array<const int>>& recv_cell_node_number_by_proc);
+  void _buildSubItemNumberToIdMap();
+  template <ItemType item_type>
+  void _buildItemToExchangeLists();
+  template <ItemType item_type>
+  void _buildNumberOfItemToExchange();
+  template <typename ItemOfItemT>
+  void _buildItemToSubItemDescriptor();
+  void _dispatchEdges();
  void _dispatchFaces();
  template<typename DataType, ItemType item_type, typename ConnectivityPtr>
  void _gatherFrom(const ItemValue<DataType, item_type, ConnectivityPtr>& data_to_gather,
                   std::vector<std::remove_const_t<DataType>>& gathered_vector);
+  template<typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+  void _gatherFrom(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& data_to_gather,
+                   std::vector<Array<std::remove_const_t<DataType>>>& gathered_vector);
  template <typename SubItemOfItemT>
-  std::vector<Array<const int>>
+  void _buildNumberOfSubItemPerItemToRecvByProc();
-  _getRecvNumberOfSubItemPerItemByProc();
  template <typename SubItemOfItemT>
-  std::vector<Array<const int>>
+  void _buildSubItemNumbersToRecvByProc();
-  _getRecvItemSubItemNumberingByProc(const std::vector<Array<const int>>& recv_number_of_sub_item_per_item_by_proc);
  template <ItemType item_type>
  void _buildRecvItemIdCorrespondanceByProc();
+  template <ItemType item_type>
+  void _buildItemReferenceList();
 public:
  std::shared_ptr<const ConnectivityType>
  dispatchedConnectivity() const

--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -262,6 +262,9 @@ private:
  template <size_t Dimension>
  void __computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor);
+  template <size_t Dimension>
+  void __computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor);
  template <int Dimension>
  void _dispatch();

--- a/src/mesh/ItemToItemMatrix.hpp
+++ b/src/mesh/ItemToItemMatrix.hpp
@@ -83,7 +83,7 @@ class ItemToItemMatrix
  PASTIS_INLINE
  const auto& operator[](const IndexType& source_id) const
  {
-    static_assert(std::is_same<IndexType, SourceItemId>(),
+    static_assert(std::is_same_v<IndexType, SourceItemId>,
                  "ItemToItemMatrix must be indexed using correct ItemId");
    using RowType = decltype(m_connectivity_matrix.rowConst(source_id));
    return SubItemList<RowType>(m_connectivity_matrix.rowConst(source_id));

--- a/src/mesh/ItemType.hpp
+++ b/src/mesh/ItemType.hpp
@@ -118,4 +118,8 @@ struct ItemTypeId<3>
  }
 };
+template <ItemType item_type>
+PASTIS_INLINE
+constexpr bool is_false_item_type_v = false;
 #endif // ITEM_TYPE_HPP
--- a/src/mesh/ItemValue.hpp
+++ b/src/mesh/ItemValue.hpp
@@ -95,13 +95,10 @@ class ItemValue
  }
  template <typename IndexType>
-  DataType& operator[](const IndexType& i) const noexcept(NO_ASSERT)
+  DataType& operator[](const IndexType&) const noexcept(NO_ASSERT)
  {
    static_assert(std::is_same_v<IndexType,ItemId>,
                  "ItemValue must be indexed by ItemId");
-    static_assert(not std::is_const_v<DataType>,
-                  "Cannot modify ItemValue of const");
-    return m_values[i];
  }
  PASTIS_INLINE

--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -11,7 +11,6 @@
 struct IMesh
 {
  virtual size_t dimension() const = 0;
-  // virtual CSRGraph cellToCellGraph() const = 0;
  ~IMesh() = default;
 };
@@ -66,7 +65,6 @@ public:
    return m_xr;
  }
-  [[deprecated("should rework this class: quite ugly")]]
  PASTIS_INLINE
  NodeValue<Rd> mutableXr() const
  {

--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -178,13 +178,11 @@ class MeshData
              const FaceId& l = cell_faces[L];
              const auto& face_nodes = face_to_node_matrix[l];
-#warning should this lambda be replaced by a precomputed correspondance?
              auto local_node_number_in_cell
                  = [&](const NodeId& node_number) {
                      for (size_t i_node=0; i_node<cell_nodes.size(); ++i_node) {
                        if (node_number == cell_nodes[i_node]) {
                          return i_node;
-                          break;
                        }
                      }
                      return std::numeric_limits<size_t>::max();
@@ -227,31 +225,37 @@ class MeshData
  }
 public:
+  PASTIS_INLINE
  const MeshType& mesh() const
  {
    return m_mesh;
  }
+  PASTIS_INLINE
  const NodeValuePerCell<const Rd>& Cjr() const
  {
    return m_Cjr;
  }
+  PASTIS_INLINE
  const NodeValuePerCell<const double>& ljr() const
  {
    return m_ljr;
  }
+  PASTIS_INLINE
  const NodeValuePerCell<const Rd>& njr() const
  {
    return m_njr;
  }
+  PASTIS_INLINE
  const CellValue<const Rd>& xj() const
  {
    return m_xj;
  }
+  PASTIS_INLINE
  const CellValue<const double>& Vj() const
  {
    return m_Vj;

--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -7,8 +7,7 @@
 #include <Kokkos_Vector.hpp>
 #include <TinyVector.hpp>
-#include <RefNodeList.hpp>
+#include <RefItemList.hpp>
-#include <RefFaceList.hpp>
 #include <ConnectivityMatrix.hpp>
 #include <IConnectivity.hpp>
@@ -38,7 +37,7 @@ class MeshNodeBoundary
    const auto& face_to_cell_matrix
        = mesh.connectivity().faceToCellMatrix();
-    const Array<const FaceId>& face_list = ref_face_list.faceList();
+    const Array<const FaceId>& face_list = ref_face_list.list();
    parallel_for(face_list.size(), PASTIS_LAMBDA(const int& l){
        const auto& face_cells = face_to_cell_matrix[face_list[l]];
        if (face_cells.size()>1) {
@@ -74,7 +73,7 @@ class MeshNodeBoundary
  template <typename MeshType>
  MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list)
-      : m_node_list(ref_node_list.nodeList())
+      : m_node_list(ref_node_list.list())
  {
    static_assert(Dimension == MeshType::Dimension);
  }
@@ -288,7 +287,6 @@ _getNormal(const MeshType& mesh)
      zmax = x;
    }
  }
-#warning re work this part to avoir parallelism dependance
  Array<R3> xmin_array = parallel::allGather(xmin);
  Array<R3> xmax_array = parallel::allGather(xmax);
  Array<R3> ymin_array = parallel::allGather(ymin);
@@ -321,7 +319,6 @@ _getNormal(const MeshType& mesh)
    if (x[2] > zmax[2]) { zmax = x; }
  }
  const R3 u = xmax-xmin;
  const R3 v = ymax-ymin;
  const R3 w = zmax-zmin;
@@ -356,7 +353,6 @@ _getNormal(const MeshType& mesh)
  normal *= 1./sqrt(normal_l2);
-#warning Add flatness test
  // this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh);
  return normal;

--- a/src/mesh/RefFaceList.hpp
+++ b/src/mesh/RefFaceList.hpp
-#ifndef REF_FACE_LIST_HPP
-#define REF_FACE_LIST_HPP
-#include <Array.hpp>
-#include <RefId.hpp>
-class RefFaceList
-{
- private:
-  RefId m_ref_id;
-  Array<const FaceId> m_face_id_list;
- public:
-  const RefId& refId() const
-  {
-    return m_ref_id;
-  }
-  const Array<const FaceId>& faceList() const
-  {
-    return m_face_id_list;
-  }
-  RefFaceList(const RefId& ref_id,
-              const Array<const FaceId>& face_id_list)
-      : m_ref_id(ref_id),
-        m_face_id_list(face_id_list)
-  {
-    ;
-  }
-  RefFaceList& operator=(const RefFaceList&) = default;
-  RefFaceList& operator=(RefFaceList&&) = default;
-  RefFaceList() = default;
-  RefFaceList(const RefFaceList&) = default;
-  ~RefFaceList() = default;
-};
-#endif // REF_FACE_LIST_HPP
--- a/src/mesh/RefItemList.hpp
+++ b/src/mesh/RefItemList.hpp
+#ifndef REF_ITEM_LIST_HPP
+#define REF_ITEM_LIST_HPP
+#include <Array.hpp>
+#include <RefId.hpp>
+#include <ItemId.hpp>
+template <ItemType item_type>
+class RefItemList
+{
+ public:
+  using ItemId = ItemIdT<item_type>;
+ private:
+  RefId m_ref_id;
+  Array<const ItemId> m_item_id_list;
+ public:
+  const RefId& refId() const
+  {
+    return m_ref_id;
+  }
+  const Array<const ItemId>& list() const
+  {
+    return m_item_id_list;
+  }
+  RefItemList(const RefId& ref_id,
+              const Array<const ItemId>& item_id_list)
+      : m_ref_id(ref_id),
+        m_item_id_list(item_id_list)
+  {
+    ;
+  }
+  RefItemList& operator=(const RefItemList&) = default;
+  RefItemList& operator=(RefItemList&&) = default;
+  RefItemList() = default;
+  RefItemList(const RefItemList&) = default;
+  ~RefItemList() = default;
+};
+using RefNodeList = RefItemList<ItemType::node>;
+using RefEdgeList = RefItemList<ItemType::edge>;
+using RefFaceList = RefItemList<ItemType::face>;
+using RefCellList = RefItemList<ItemType::cell>;
+#endif // REF_ITEM_LIST_HPP
--- a/src/mesh/RefNodeList.hpp
+++ b/src/mesh/RefNodeList.hpp
-#ifndef REF_NODE_LIST_HPP
-#define REF_NODE_LIST_HPP
-#include <Array.hpp>
-#include <RefId.hpp>
-class RefNodeList
-{
- private:
-  RefId m_ref_id;
-  Array<const NodeId> m_node_id_list;
- public:
-  const RefId& refId() const
-  {
-    return m_ref_id;
-  }
-  const Array<const NodeId>& nodeList() const
-  {
-    return m_node_id_list;
-  }
-  RefNodeList(const RefId& ref_id,
-              const Array<const NodeId>& node_id_list)
-      : m_ref_id(ref_id),
-        m_node_id_list(node_id_list)
-  {
-    ;
-  }
-  RefNodeList& operator=(const RefNodeList&) = default;
-  RefNodeList& operator=(RefNodeList&&) = default;
-  RefNodeList() = default;
-  RefNodeList(const RefNodeList&) = default;
-  ~RefNodeList() = default;
-};
-#endif // REF_NODE_LIST_HPP
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -17,20 +17,10 @@
 #include <memory>
-// template <typename DataType,
-//           typename ItemOfItem,
-//           typename ConnectivityPtr = std::shared_ptr<const IConnectivity>,
-//           typename Allowed=void>
-// class SubItemValuePerItem;
 template <typename DataType,
          typename ItemOfItem,
          typename ConnectivityPtr = std::shared_ptr<const IConnectivity>>
-class SubItemValuePerItem// <DataType,
+class SubItemValuePerItem
-                         //  ItemOfItem,
-                         //  ConnectivityPtr// ,
-                         //  // std::enable_if_t<sub_item_type != item_type>
-                         //  >
 {
 public:
  static constexpr ItemType item_type{ItemOfItem::item_type};
No results found