Generalize diamond dual mesh to 3D

src/mesh/DiamondDualMeshBuilder.cpp

@@ -12,32 +12,9 @@
 
 template <size_t Dimension>
 void
-DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>&)
+DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<Dimension>& primal_connectivity,
+                                                            ConnectivityDescriptor& diamond_descriptor)
 {
-  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();
 
@@ -45,18 +22,19 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMe
 
   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();
+  const auto& primal_node_number = primal_connectivity.nodeNumber();
 
-  // 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];
-  // }
+  for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
+    diamond_descriptor.node_number_vector[primal_node_id] = primal_node_number[primal_node_id];
+  }
+
+  const auto& primal_cell_number = primal_connectivity.cellNumber();
 
-#warning fix max and compute diamond node numbers correctly in parallel
-  //  const size_t max_node_number =  max(primal_node_number);
+  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;
+  for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
+    diamond_descriptor.node_number_vector[primal_number_of_nodes + primal_cell_id] =
+      primal_cell_number[primal_cell_id] + max_node_number;
   }
 
   const size_t diamond_number_of_cells = primal_connectivity.numberOfFaces();
@@ -78,12 +56,25 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMe
     const size_t i_cell               = i_face;
     const auto& primal_face_cell_list = primal_face_to_cell_matrix[i_face];
 
+    if constexpr (Dimension == 1) {
+      throw NotImplementedError("dimension 1");
+    } else if constexpr (Dimension == 2) {
       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;
       }
+    } else {
+      static_assert(Dimension == 3, "unexpected dimension");
+
+      if (primal_face_cell_list.size() == 1) {
+        diamond_descriptor.cell_type_vector[i_cell] = CellType::Pyramid;
+      } else {
+        Assert(primal_face_cell_list.size() == 2);
+        diamond_descriptor.cell_type_vector[i_cell] = CellType::Diamond;
+      }
+    }
   }
 
   diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);
@@ -96,27 +87,40 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMe
     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);
+      diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(primal_face_node_list.size() + 1);
 
       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;
+      for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
+        diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node] = primal_face_node_list[i_node];
+      }
+      diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size()] =
+        primal_number_of_nodes + cell_id;
 
       if (cell_face_is_reversed(cell_id, i_face_in_cell)) {
+        if constexpr (Dimension == 2) {
           std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][0],
                     diamond_descriptor.cell_to_node_vector[i_diamond_cell][1]);
+
+        } else {
+          for (size_t i_node = 0; i_node < primal_face_node_list.size() / 2; ++i_node) {
+            std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node],
+                      diamond_descriptor
+                        .cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() - 1 - i_node]);
+          }
+        }
       }
     } else {
       Assert(primal_face_cell_list.size() == 2);
-      diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(4);
+      diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(primal_face_node_list.size() + 2);
 
       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);
 
+      if constexpr (Dimension == 2) {
+        Assert(primal_face_node_list.size() == 2);
         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;
@@ -126,13 +130,45 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMe
           std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][1],
                     diamond_descriptor.cell_to_node_vector[i_diamond_cell][3]);
         }
+      } else {
+        diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = primal_number_of_nodes + cell0_id;
+        for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
+          diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node + 1] = primal_face_node_list[i_node];
+        }
+        diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() + 1] =
+          primal_number_of_nodes + cell1_id;
+
+        if (cell_face_is_reversed(cell0_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][primal_face_node_list.size() + 1]);
+        }
+      }
     }
   }
+}
+
+template <size_t Dimension>
+void
+DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>& p_mesh)
+{
+  static_assert(Dimension <= 3, "invalid mesh dimension");
+  using ConnectivityType = Connectivity<Dimension>;
+  using MeshType         = Mesh<ConnectivityType>;
+
+  const MeshType& primal_mesh = dynamic_cast<const MeshType&>(*p_mesh);
+
+  const ConnectivityType& primal_connectivity = primal_mesh.connectivity();
+
+  ConnectivityDescriptor diamond_descriptor;
+
+  this->_buildDiamondConnectivityDescriptor(primal_connectivity, diamond_descriptor);
 
-  MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(diamond_descriptor);
+  MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(diamond_descriptor);
 
   {
-    using Face = ConnectivityFace<2>;
+    const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix();
+
+    using Face = ConnectivityFace<Dimension>;
 
     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;
@@ -144,9 +180,9 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMe
       return face_to_id_map;
     }();
 
-    for (size_t i_face_list = 0; i_face_list < primal_connectivity.numberOfRefItemList<ItemType::face>();
+    for (size_t i_face_list = 0; i_face_list < primal_connectivity.template numberOfRefItemList<ItemType::face>();
          ++i_face_list) {
-      const auto& primal_ref_face_list = primal_connectivity.refItemList<ItemType::face>(i_face_list);
+      const auto& primal_ref_face_list = primal_connectivity.template 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();
 
@@ -183,24 +219,55 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMe
     }
   }
 
-#warning completly wrong
+  const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
+  const size_t primal_number_of_cells = primal_connectivity.numberOfCells();
+
   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());
+
+  const auto& primal_node_owner = primal_connectivity.nodeOwner();
+  for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
+    diamond_descriptor.node_owner_vector[primal_node_id] = primal_node_owner[primal_node_id];
+  }
+  const auto& primal_cell_owner = primal_connectivity.cellOwner();
+  for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
+    diamond_descriptor.node_owner_vector[primal_number_of_nodes + primal_cell_id] = primal_cell_owner[primal_cell_id];
+  }
 
   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());
+  const auto& primal_face_owner = primal_connectivity.faceOwner();
+  for (FaceId primal_face_id = 0; primal_face_id < primal_number_of_cells; ++primal_face_id) {
+    diamond_descriptor.cell_owner_vector[primal_face_id] = primal_face_owner[primal_face_id];
+  }
 
-  std::shared_ptr p_diamond_connectivity = Connectivity2D::build(diamond_descriptor);
-  Connectivity2D& diamond_connectivity   = *p_diamond_connectivity;
+  {
+    std::vector<int> face_cell_owner(diamond_descriptor.face_number_vector.size());
+    std::fill(std::begin(face_cell_owner), std::end(face_cell_owner), parallel::size());
+
+    for (size_t i_cell = 0; i_cell < diamond_descriptor.cell_to_face_vector.size(); ++i_cell) {
+      const auto& cell_face_list = diamond_descriptor.cell_to_face_vector[i_cell];
+      for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+        const size_t face_id     = cell_face_list[i_face];
+        face_cell_owner[face_id] = std::min(face_cell_owner[face_id], diamond_descriptor.cell_number_vector[i_cell]);
+      }
+    }
+
+    diamond_descriptor.face_owner_vector.resize(face_cell_owner.size());
+    for (size_t i_face = 0; i_face < face_cell_owner.size(); ++i_face) {
+      diamond_descriptor.face_owner_vector[i_face] = diamond_descriptor.cell_owner_vector[face_cell_owner[i_face]];
+    }
+  }
 
-  NodeValue<TinyVector<2>> diamond_xr{diamond_connectivity};
+  std::shared_ptr p_diamond_connectivity = ConnectivityType::build(diamond_descriptor);
+  ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
+
+  NodeValue<TinyVector<Dimension>> diamond_xr{diamond_connectivity};
 
   const auto primal_xr = primal_mesh.xr();
-  MeshData<Mesh<Connectivity2D>> primal_mesh_data{primal_mesh};
+  MeshData<MeshType> primal_mesh_data{primal_mesh};
   const auto primal_xj = primal_mesh_data.xj();
 
   {
-#warning to recode
+#warning define transfer functions
     NodeId i_node = 0;
     for (; i_node < primal_number_of_nodes; ++i_node) {
       diamond_xr[i_node] = primal_xr[i_node];
@@ -211,9 +278,11 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMe
     }
   }
 
-  std::shared_ptr p_diamond_mesh = std::make_shared<Mesh<Connectivity2D>>(p_diamond_connectivity, diamond_xr);
+  std::shared_ptr p_diamond_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
 
-  MeshData<Mesh<Connectivity2D>> dual_mesh_data{*p_diamond_mesh};
+#warning USELESS TEST
+  // -->>
+  MeshData<MeshType> dual_mesh_data{*p_diamond_mesh};
 
   double sum = 0;
   for (CellId cell_id = 0; cell_id < p_diamond_mesh->numberOfCells(); ++cell_id) {
@@ -221,15 +290,16 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMe
   }
 
   std::cout << "volume = " << sum << '\n';
+  // <<--
 
-  m_mesh = std::make_shared<Mesh<Connectivity2D>>(p_diamond_connectivity, diamond_xr);
+  m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
 }
 
 template <>
-void
-DiamondDualMeshBuilder::_buildDiamondMeshFrom<3>(const std::shared_ptr<const IMesh>& p_mesh)
+[[deprecated]] void
+DiamondDualMeshBuilder::_buildDiamondMeshFrom<1>(const std::shared_ptr<const IMesh>&)
 {
-  m_mesh = p_mesh;
+  m_mesh = 0;
 }
 
 DiamondDualMeshBuilder::DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>& p_mesh)

src/mesh/DiamondDualMeshBuilder.hpp

@@ -6,9 +6,17 @@
 
 #include <memory>
 
+template <size_t>
+class Connectivity;
+
+class ConnectivityDescriptor;
+
 class DiamondDualMeshBuilder : public MeshBuilderBase
 {
  private:
+  template <size_t Dimension>
+  void _buildDiamondConnectivityDescriptor(const Connectivity<Dimension>&, ConnectivityDescriptor&);
+
   template <size_t Dimension>
   void _buildDiamondMeshFrom(const std::shared_ptr<const IMesh>&);