Continue Diamond mesh construction

- ItemValueUtils: min/max/sum functions require improvement or fixes
- Diamond meshes are built correctly in 2D and sequential
  - faces boundaries references lists are correctly translated
  - node boundaries (ref lists) are missing
  - 1D and 3D are to be done

src/mesh/DiamondDualMeshBuilder.cpp

@@ -3,16 +3,232 @@
 #include <mesh/Connectivity.hpp>
 #include <mesh/ConnectivityDescriptor.hpp>
 #include <mesh/ConnectivityDispatcher.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
 #include <mesh/RefId.hpp>
 #include <utils/Array.hpp>
 #include <utils/Messenger.hpp>
 
 template <size_t Dimension>
 void
-DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>& p_mesh)
+DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>&)
 {
-  static_assert(Dimension <= 3, "invalid dimension");
+  static_assert(Dimension <= 3, "invalid mesh dimension");
+}
+
+template <>
+void
+DiamondDualMeshBuilder::_buildDiamondMeshFrom<1>(const std::shared_ptr<const IMesh>& p_mesh)
+{
+  m_mesh = p_mesh;
+}
+
+template <>
+void
+DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMesh>& p_mesh)
+{
+  m_mesh = p_mesh;
+
+  using MeshType = Mesh<Connectivity2D>;
+
+  const MeshType& primal_mesh = dynamic_cast<const MeshType&>(*p_mesh);
+
+  const Connectivity2D& primal_connectivity = primal_mesh.connectivity();
+
+  ConnectivityDescriptor diamond_descriptor;
+
+  const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
+  const size_t primal_number_of_cells = primal_connectivity.numberOfCells();
+
+  const size_t diamond_number_of_nodes = primal_number_of_cells + primal_number_of_nodes;
+
+  diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes);
+
+  // const auto& primal_node_number = primal_connectivity.nodeNumber();
+  // const auto& primal_cell_number = primal_connectivity.cellNumber();
+
+  // for (size_t i_primal_node = 0; i_primal_node < primal_connectivity.numberOfNodes(); ++i_primal_node) {
+  //   diamond_descriptor.node_number_vector[i_primal_node] = primal_node_number[i_primal_node];
+  // }
+
+#warning fix max and compute diamond node numbers correctly in parallel
+  //  const size_t max_node_number =  max(primal_node_number);
+
+  for (size_t i = 0; i < diamond_number_of_nodes; ++i) {
+    diamond_descriptor.node_number_vector[i] = i;
+  }
+
+  const size_t diamond_number_of_cells = primal_connectivity.numberOfFaces();
+  diamond_descriptor.cell_number_vector.resize(diamond_number_of_cells);
+
+  const auto& primal_face_number = primal_connectivity.faceNumber();
+
+  for (FaceId i_primal_face = 0; i_primal_face < primal_connectivity.numberOfFaces(); ++i_primal_face) {
+    diamond_descriptor.cell_number_vector[i_primal_face] = primal_face_number[i_primal_face];
+  }
+
+  diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);
+
+  diamond_descriptor.cell_type_vector.resize(diamond_number_of_cells);
+
+  const auto& primal_face_to_cell_matrix = primal_connectivity.faceToCellMatrix();
+
+  for (FaceId i_face = 0; i_face < primal_connectivity.numberOfFaces(); ++i_face) {
+    const size_t i_cell               = i_face;
+    const auto& primal_face_cell_list = primal_face_to_cell_matrix[i_face];
+
+    if (primal_face_cell_list.size() == 1) {
+      diamond_descriptor.cell_type_vector[i_cell] = CellType::Triangle;
+    } else {
+      Assert(primal_face_cell_list.size() == 2);
+      diamond_descriptor.cell_type_vector[i_cell] = CellType::Quadrangle;
+    }
+  }
+
+  diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);
+
+  const auto& primal_face_to_node_matrix              = primal_connectivity.faceToNodeMatrix();
+  const auto& primal_face_local_number_in_their_cells = primal_connectivity.faceLocalNumbersInTheirCells();
+  const auto& cell_face_is_reversed                   = primal_connectivity.cellFaceIsReversed();
+  for (FaceId i_face = 0; i_face < primal_connectivity.numberOfFaces(); ++i_face) {
+    const size_t& i_diamond_cell      = i_face;
+    const auto& primal_face_cell_list = primal_face_to_cell_matrix[i_face];
+    const auto& primal_face_node_list = primal_face_to_node_matrix[i_face];
+    if (primal_face_cell_list.size() == 1) {
+      diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(3);
+
+      const CellId cell_id      = primal_face_cell_list[0];
+      const auto i_face_in_cell = primal_face_local_number_in_their_cells(i_face, 0);
+
+      diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = primal_face_node_list[0];
+      diamond_descriptor.cell_to_node_vector[i_diamond_cell][1] = primal_face_node_list[1];
+      diamond_descriptor.cell_to_node_vector[i_diamond_cell][2] = primal_number_of_nodes + cell_id;
+
+      if (cell_face_is_reversed(cell_id, i_face_in_cell)) {
+        std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][0],
+                  diamond_descriptor.cell_to_node_vector[i_diamond_cell][1]);
+      }
+    } else {
+      Assert(primal_face_cell_list.size() == 2);
+      diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(4);
+
+      const CellId cell0_id     = primal_face_cell_list[0];
+      const CellId cell1_id     = primal_face_cell_list[1];
+      const auto i_face_in_cell = primal_face_local_number_in_their_cells(i_face, 0);
+
+      diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = primal_number_of_nodes + cell0_id;
+      diamond_descriptor.cell_to_node_vector[i_diamond_cell][1] = primal_face_node_list[0];
+      diamond_descriptor.cell_to_node_vector[i_diamond_cell][2] = primal_number_of_nodes + cell1_id;
+      diamond_descriptor.cell_to_node_vector[i_diamond_cell][3] = primal_face_node_list[1];
+
+      if (cell_face_is_reversed(cell0_id, i_face_in_cell)) {
+        std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][1],
+                  diamond_descriptor.cell_to_node_vector[i_diamond_cell][3]);
+      }
+    }
+  }
+
+  MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(diamond_descriptor);
+
+  {
+    using Face = ConnectivityFace<2>;
+
+    const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map = [&] {
+      std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
+      for (FaceId l = 0; l < diamond_descriptor.face_to_node_vector.size(); ++l) {
+        const auto& node_vector = diamond_descriptor.face_to_node_vector[l];
+
+        face_to_id_map[Face(node_vector, diamond_descriptor.node_number_vector)] = l;
+      }
+      return face_to_id_map;
+    }();
+
+    for (size_t i_face_list = 0; i_face_list < primal_connectivity.numberOfRefItemList<ItemType::face>();
+         ++i_face_list) {
+      const auto& primal_ref_face_list = primal_connectivity.refItemList<ItemType::face>(i_face_list);
+      std::cout << "treating " << primal_ref_face_list.refId() << '\n';
+      const auto& primal_face_list = primal_ref_face_list.list();
+
+      const std::vector<FaceId> diamond_face_list = [&]() {
+        std::vector<FaceId> diamond_face_list;
+        diamond_face_list.reserve(primal_face_list.size());
+        for (size_t i = 0; i < primal_face_list.size(); ++i) {
+          FaceId primal_face_id = primal_face_list[i];
+
+          const auto& primal_face_node_list = primal_face_to_node_matrix[primal_face_id];
+
+          const auto i_diamond_face = [&]() {
+            std::vector<unsigned int> node_list(primal_face_node_list.size());
+            for (size_t i = 0; i < primal_face_node_list.size(); ++i) {
+              node_list[i] = primal_face_node_list[i];
+            }
+            return face_to_id_map.find(Face(node_list, diamond_descriptor.node_number_vector));
+          }();
+
+          if (i_diamond_face != face_to_id_map.end()) {
+            diamond_face_list.push_back(i_diamond_face->second);
+          }
+        }
+        return diamond_face_list;
+      }();
+
+      if (parallel::allReduceOr(diamond_face_list.size() > 0)) {
+        Array<FaceId> face_array(diamond_face_list.size());
+        for (size_t i = 0; i < diamond_face_list.size(); ++i) {
+          face_array[i] = diamond_face_list[i];
+        }
+        diamond_descriptor.addRefItemList(RefFaceList{primal_ref_face_list.refId(), face_array});
+      }
+    }
+  }
+
+#warning completly wrong
+  diamond_descriptor.node_owner_vector.resize(diamond_descriptor.node_number_vector.size());
+  std::fill(diamond_descriptor.node_owner_vector.begin(), diamond_descriptor.node_owner_vector.end(), parallel::rank());
+
+  diamond_descriptor.cell_owner_vector.resize(diamond_descriptor.cell_number_vector.size());
+  std::fill(diamond_descriptor.cell_owner_vector.begin(), diamond_descriptor.cell_owner_vector.end(), parallel::rank());
+
+  std::shared_ptr p_diamond_connectivity = Connectivity2D::build(diamond_descriptor);
+  Connectivity2D& diamond_connectivity   = *p_diamond_connectivity;
+
+  NodeValue<TinyVector<2>> diamond_xr{diamond_connectivity};
+
+  const auto primal_xr = primal_mesh.xr();
+  MeshData<Mesh<Connectivity2D>> primal_mesh_data{primal_mesh};
+  const auto primal_xj = primal_mesh_data.xj();
 
+  {
+#warning to recode
+    NodeId i_node = 0;
+    for (; i_node < primal_number_of_nodes; ++i_node) {
+      diamond_xr[i_node] = primal_xr[i_node];
+    }
+
+    for (CellId i_cell = 0; i_cell < primal_number_of_cells; ++i_cell) {
+      diamond_xr[i_node++] = primal_xj[i_cell];
+    }
+  }
+
+  std::shared_ptr p_diamond_mesh = std::make_shared<Mesh<Connectivity2D>>(p_diamond_connectivity, diamond_xr);
+
+  MeshData<Mesh<Connectivity2D>> dual_mesh_data{*p_diamond_mesh};
+
+  double sum = 0;
+  for (CellId cell_id = 0; cell_id < p_diamond_mesh->numberOfCells(); ++cell_id) {
+    sum += dual_mesh_data.Vj()[cell_id];
+  }
+
+  std::cout << "volume = " << sum << '\n';
+
+  m_mesh = std::make_shared<Mesh<Connectivity2D>>(p_diamond_connectivity, diamond_xr);
+}
+
+template <>
+void
+DiamondDualMeshBuilder::_buildDiamondMeshFrom<3>(const std::shared_ptr<const IMesh>& p_mesh)
+{
   m_mesh = p_mesh;
 }
 

src/mesh/DiamondDualMeshBuilder.hpp

@@ -10,10 +10,10 @@ class DiamondDualMeshBuilder : public MeshBuilderBase
 {
  private:
   template <size_t Dimension>
-  void _buildDiamondMeshFrom(const std::shared_ptr<const IMesh>& mesh);
+  void _buildDiamondMeshFrom(const std::shared_ptr<const IMesh>&);
 
  public:
-  DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>& mesh);
+  DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>&);
   ~DiamondDualMeshBuilder() = default;
 };
 

src/mesh/ItemValueUtils.hpp

@@ -10,11 +10,11 @@
 
 #include <iostream>
 
-template <typename DataType, ItemType item_type>
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
 std::remove_const_t<DataType>
-min(const ItemValue<DataType, item_type>& item_value)
+min(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
 {
-  using ItemValueType = ItemValue<DataType, item_type>;
+  using ItemValueType = ItemValue<DataType, item_type, ConnectivityPtr>;
   using data_type     = std::remove_const_t<typename ItemValueType::data_type>;
   using index_type    = typename ItemValueType::index_type;
 
@@ -100,11 +100,11 @@ min(const ItemValue<DataType, item_type>& item_value)
   return parallel::allReduceMin(local_min);
 }
 
-template <typename DataType, ItemType item_type>
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
 std::remove_const_t<DataType>
-max(const ItemValue<DataType, item_type>& item_value)
+max(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
 {
-  using ItemValueType = ItemValue<DataType, item_type>;
+  using ItemValueType = ItemValue<DataType, item_type, ConnectivityPtr>;
   using data_type     = std::remove_const_t<typename ItemValueType::data_type>;
   using index_type    = typename ItemValueType::index_type;
 
@@ -268,8 +268,7 @@ sum(const ItemValue<DataType, item_type>& item_value)
     }
 
     PUGS_INLINE
-    ItemValueSum(const ItemValueType& item_value, const IsOwnedType& is_owned)
-      : m_item_value(item_value), m_is_owned(is_owned)
+    ItemValueSum(const ItemValueType& item_value) : m_item_value(item_value), m_is_owned{is_owned(item_value)}
     {
       ;
     }