DiamondDualConnectivityBuilder.cpp

#include <mesh/DiamondDualConnectivityBuilder.hpp>

#include <mesh/Connectivity.hpp>
#include <mesh/ConnectivityDescriptor.hpp>
#include <mesh/ConnectivityDispatcher.hpp>
#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/RefId.hpp>
#include <utils/Array.hpp>
#include <utils/ArrayUtils.hpp>
#include <utils/Messenger.hpp>

#include <vector>

template <size_t Dimension>
void
DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<Dimension>& primal_connectivity,
                                                                    ConnectivityDescriptor& diamond_descriptor)
{
  const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
  const size_t primal_number_of_faces = primal_connectivity.numberOfFaces();
  const size_t primal_number_of_cells = primal_connectivity.numberOfCells();

  const auto& primal_node_number = primal_connectivity.nodeNumber();
  const auto& primal_cell_number = primal_connectivity.cellNumber();

  const size_t diamond_number_of_nodes = primal_number_of_cells + primal_number_of_nodes;
  const size_t diamond_number_of_cells = primal_number_of_faces;

  {
    m_primal_node_to_dual_node_map = NodeIdToNodeIdMap{primal_number_of_nodes};

    NodeId diamond_node_id = 0;
    for (NodeId primal_node_id = 0; primal_node_id < primal_number_of_nodes; ++primal_node_id) {
      m_primal_node_to_dual_node_map[primal_node_id] = std::make_pair(primal_node_id, diamond_node_id++);
    }

    m_primal_cell_to_dual_node_map = CellIdToNodeIdMap{primal_number_of_cells};
    for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
      m_primal_cell_to_dual_node_map[primal_cell_id] = std::make_pair(primal_cell_id, diamond_node_id++);
    }
  }

  diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes);

  parallel_for(m_primal_node_to_dual_node_map.size(), [&](size_t i) {
    const auto [primal_node_id, diamond_dual_node_id] = m_primal_node_to_dual_node_map[i];

    diamond_descriptor.node_number_vector[diamond_dual_node_id] = primal_node_number[primal_node_id];
  });

  const size_t cell_number_shift = max(primal_node_number) + 1;
  parallel_for(primal_number_of_cells, [&](size_t i) {
    const auto [primal_cell_id, diamond_dual_node_id] = m_primal_cell_to_dual_node_map[i];

    diamond_descriptor.node_number_vector[diamond_dual_node_id] =
      primal_cell_number[primal_cell_id] + cell_number_shift;
  });

  {
    m_primal_face_to_dual_cell_map = FaceIdToCellIdMap{primal_number_of_faces};
    CellId diamond_cell_id         = 0;
    for (FaceId primal_face_id = 0; primal_face_id < primal_number_of_faces; ++primal_face_id) {
      m_primal_face_to_dual_cell_map[primal_face_id] = std::make_pair(primal_face_id, diamond_cell_id++);
    }
  }

  diamond_descriptor.cell_number_vector.resize(diamond_number_of_cells);
  const auto& primal_face_number = primal_connectivity.faceNumber();
  parallel_for(diamond_number_of_cells, [&](size_t i) {
    const auto [primal_face_id, dual_cell_id]           = m_primal_face_to_dual_cell_map[i];
    diamond_descriptor.cell_number_vector[dual_cell_id] = primal_face_number[primal_face_id];
  });

  if constexpr (Dimension == 3) {
    const size_t number_of_edges = diamond_descriptor.edge_to_node_vector.size();
    diamond_descriptor.edge_number_vector.resize(number_of_edges);
    for (size_t i_edge = 0; i_edge < number_of_edges; ++i_edge) {
      diamond_descriptor.edge_number_vector[i_edge] = i_edge;
    }
    if (parallel::size() > 1) {
      throw NotImplementedError("parallel edge numbering is undefined");
    }
  }

  diamond_descriptor.cell_type_vector.resize(diamond_number_of_cells);

  const auto& primal_face_to_cell_matrix = primal_connectivity.faceToCellMatrix();

  parallel_for(primal_number_of_faces, [&](FaceId face_id) {
    const size_t i_cell               = face_id;
    const auto& primal_face_cell_list = primal_face_to_cell_matrix[face_id];

    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;
    }

    if constexpr (Dimension == 3) {
      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);

  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();
  parallel_for(primal_number_of_faces, [&](FaceId face_id) {
    const size_t& i_diamond_cell      = face_id;
    const auto& primal_face_cell_list = primal_face_to_cell_matrix[face_id];
    const auto& primal_face_node_list = primal_face_to_node_matrix[face_id];
    if (primal_face_cell_list.size() == 1) {
      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(face_id, 0);

      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(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_cell0 = primal_face_local_number_in_their_cells(face_id, 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;
        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_cell0)) {
          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_cell0)) {
          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 <>
void
DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor<1>(const Connectivity<1>& primal_connectivity,
                                                                       ConnectivityDescriptor& diamond_descriptor)
{
  const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
  const size_t primal_number_of_faces = primal_connectivity.numberOfFaces();
  const size_t primal_number_of_cells = primal_connectivity.numberOfCells();

  const size_t diamond_number_of_nodes = 2 * primal_number_of_nodes - primal_number_of_cells;

  const size_t diamond_number_of_cells = primal_number_of_faces;
  const size_t number_of_kept_nodes    = 2 * (diamond_number_of_nodes - diamond_number_of_cells);

  const auto& primal_node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
  size_t next_kept_node_id               = 0;

  diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes);

  const auto& primal_node_number = primal_connectivity.nodeNumber();

  for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
    const auto& primal_node_cell_list = primal_node_to_cell_matrix[primal_node_id];
    if (primal_node_cell_list.size() == 1) {
      diamond_descriptor.node_number_vector[next_kept_node_id++] = primal_node_number[primal_node_id];
    }
  }

  const size_t primal_number_of_kept_nodes = next_kept_node_id;

  const auto& primal_cell_number = primal_connectivity.cellNumber();

  const size_t cell_number_shift = max(primal_node_number) + 1;

  for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
    diamond_descriptor.node_number_vector[primal_number_of_kept_nodes + primal_cell_id] =
      primal_cell_number[primal_cell_id] + cell_number_shift;
  }

  if (number_of_kept_nodes != next_kept_node_id) {
    throw UnexpectedError("unexpected number of kept node" + std::to_string(next_kept_node_id) +
                          " != " + std::to_string(number_of_kept_nodes));
  }

  diamond_descriptor.cell_number_vector.resize(diamond_number_of_cells);

  for (NodeId primal_node_id = 0; primal_node_id < primal_number_of_nodes; ++primal_node_id) {
    diamond_descriptor.cell_number_vector[primal_node_id] = primal_node_number[primal_node_id];
  }

  diamond_descriptor.cell_type_vector.resize(diamond_number_of_cells);

  for (NodeId node_id = 0; node_id < primal_connectivity.numberOfNodes(); ++node_id) {
    const size_t i_cell = node_id;

    diamond_descriptor.cell_type_vector[i_cell] = CellType::Line;
  }

  diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);

  const auto& primal_node_local_number_in_their_cells = primal_connectivity.nodeLocalNumbersInTheirCells();

  {
    size_t next_kept_node_id = 0;
    for (NodeId i_node = 0; i_node < primal_connectivity.numberOfNodes(); ++i_node) {
      const size_t& i_diamond_cell      = i_node;
      const auto& primal_node_cell_list = primal_node_to_cell_matrix[i_node];
      diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(2);
      if (primal_node_cell_list.size() == 1) {
        const auto i_node_in_cell = primal_node_local_number_in_their_cells(i_node, 0);

        diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node_in_cell] = next_kept_node_id++;
        diamond_descriptor.cell_to_node_vector[i_diamond_cell][1 - i_node_in_cell] =
          number_of_kept_nodes + primal_node_cell_list[0];
      } else {
        diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = number_of_kept_nodes + primal_node_cell_list[0];
        diamond_descriptor.cell_to_node_vector[i_diamond_cell][1] = number_of_kept_nodes + primal_node_cell_list[1];
      }
    }
  }
}

template <size_t Dimension>
void
DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivity& i_primal_connectivity)
{
  static_assert(Dimension <= 3, "invalid connectivity dimension");
  using ConnectivityType = Connectivity<Dimension>;

  const ConnectivityType& primal_connectivity = dynamic_cast<const ConnectivityType&>(i_primal_connectivity);

  ConnectivityDescriptor diamond_descriptor;

  this->_buildDiamondConnectivityDescriptor(primal_connectivity, diamond_descriptor);

  if constexpr (Dimension > 1) {
    ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(diamond_descriptor);
    if constexpr (Dimension > 2) {
      ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(diamond_descriptor);
    }
  }

  {
    const std::unordered_map<unsigned int, NodeId> node_to_id_map = [&] {
      std::unordered_map<unsigned int, NodeId> node_to_id_map;
      for (size_t i_node = 0; i_node < diamond_descriptor.node_number_vector.size(); ++i_node) {
        node_to_id_map[diamond_descriptor.node_number_vector[i_node]] = i_node;
      }
      return node_to_id_map;
    }();

    for (size_t i_node_list = 0; i_node_list < primal_connectivity.template numberOfRefItemList<ItemType::node>();
         ++i_node_list) {
      const auto& primal_ref_node_list = primal_connectivity.template refItemList<ItemType::node>(i_node_list);
      const auto& primal_node_list     = primal_ref_node_list.list();

      const std::vector<NodeId> diamond_node_list = [&]() {
        std::vector<NodeId> diamond_node_list;

        for (size_t i_primal_node = 0; i_primal_node < primal_node_list.size(); ++i_primal_node) {
          auto primal_node_id       = primal_node_list[i_primal_node];
          const auto i_diamond_node = node_to_id_map.find(primal_connectivity.nodeNumber()[primal_node_id]);
          if (i_diamond_node != node_to_id_map.end()) {
            diamond_node_list.push_back(i_diamond_node->second);
          }
        }

        return diamond_node_list;
      }();

      if (parallel::allReduceOr(diamond_node_list.size() > 0)) {
        Array<NodeId> node_array(diamond_node_list.size());
        for (size_t i = 0; i < diamond_node_list.size(); ++i) {
          node_array[i] = diamond_node_list[i];
        }
        diamond_descriptor.addRefItemList(RefNodeList{primal_ref_node_list.refId(), node_array});
      }
    }
  }

  if constexpr (Dimension > 1) {
    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;
      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.template numberOfRefItemList<ItemType::face>();
         ++i_face_list) {
      const auto& primal_ref_face_list = primal_connectivity.template refItemList<ItemType::face>(i_face_list);
      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});
      }
    }
  }

  if constexpr (Dimension > 2) {
    const auto& primal_edge_to_node_matrix = primal_connectivity.edgeToNodeMatrix();
    using Edge                             = ConnectivityFace<2>;

    const std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map = [&] {
      std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map;
      for (EdgeId l = 0; l < diamond_descriptor.edge_to_node_vector.size(); ++l) {
        const auto& node_vector = diamond_descriptor.edge_to_node_vector[l];
        edge_to_id_map[Edge(node_vector, diamond_descriptor.node_number_vector)] = l;
      }
      return edge_to_id_map;
    }();

    for (size_t i_edge_list = 0; i_edge_list < primal_connectivity.template numberOfRefItemList<ItemType::edge>();
         ++i_edge_list) {
      const auto& primal_ref_edge_list = primal_connectivity.template refItemList<ItemType::edge>(i_edge_list);
      const auto& primal_edge_list     = primal_ref_edge_list.list();

      const std::vector<EdgeId> diamond_edge_list = [&]() {
        std::vector<EdgeId> diamond_edge_list;
        diamond_edge_list.reserve(primal_edge_list.size());
        for (size_t i = 0; i < primal_edge_list.size(); ++i) {
          EdgeId primal_edge_id = primal_edge_list[i];

          const auto& primal_edge_node_list = primal_edge_to_node_matrix[primal_edge_id];

          const auto i_diamond_edge = [&]() {
            std::vector<unsigned int> node_list(primal_edge_node_list.size());
            for (size_t i = 0; i < primal_edge_node_list.size(); ++i) {
              node_list[i] = primal_edge_node_list[i];
            }
            return edge_to_id_map.find(Edge(node_list, diamond_descriptor.node_number_vector));
          }();

          if (i_diamond_edge != edge_to_id_map.end()) {
            diamond_edge_list.push_back(i_diamond_edge->second);
          }
        }
        return diamond_edge_list;
      }();

      if (parallel::allReduceOr(diamond_edge_list.size() > 0)) {
        Array<EdgeId> edge_array(diamond_edge_list.size());
        for (size_t i = 0; i < diamond_edge_list.size(); ++i) {
          edge_array[i] = diamond_edge_list[i];
        }
        diamond_descriptor.addRefItemList(RefEdgeList{primal_ref_edge_list.refId(), edge_array});
      }
    }
  }

  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());

  if constexpr (Dimension == 1) {
    const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
    const auto& primal_node_owner   = primal_connectivity.nodeOwner();
    size_t next_kept_node_id        = 0;
    for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
      if (node_to_cell_matrix[primal_node_id].size() == 1) {
        diamond_descriptor.node_owner_vector[next_kept_node_id++] = primal_node_owner[primal_node_id];
      }
    }
    const size_t number_of_kept_nodes = next_kept_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[number_of_kept_nodes + primal_cell_id] = primal_cell_owner[primal_cell_id];
    }
  } else {
    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];
    }
  }

  if constexpr (Dimension == 1) {
    diamond_descriptor.cell_owner_vector.resize(diamond_descriptor.cell_number_vector.size());
    const auto& primal_node_owner = primal_connectivity.nodeOwner();
    for (NodeId primal_node_id = 0; primal_node_id < primal_number_of_nodes; ++primal_node_id) {
      diamond_descriptor.cell_owner_vector[primal_node_id] = primal_node_owner[primal_node_id];
    }
  } else {
    diamond_descriptor.cell_owner_vector.resize(diamond_descriptor.cell_number_vector.size());
    const size_t primal_number_of_faces = primal_connectivity.numberOfFaces();
    const auto& primal_face_owner       = primal_connectivity.faceOwner();
    for (FaceId primal_face_id = 0; primal_face_id < primal_number_of_faces; ++primal_face_id) {
      diamond_descriptor.cell_owner_vector[primal_face_id] = primal_face_owner[primal_face_id];
    }
  }

  {
    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]];
    }
  }

  if constexpr (Dimension == 3) {
    std::vector<int> edge_cell_owner(diamond_descriptor.edge_number_vector.size());
    std::fill(std::begin(edge_cell_owner), std::end(edge_cell_owner), parallel::size());

    for (size_t i_cell = 0; i_cell < diamond_descriptor.cell_to_face_vector.size(); ++i_cell) {
      const auto& cell_edge_list = diamond_descriptor.cell_to_edge_vector[i_cell];
      for (size_t i_edge = 0; i_edge < cell_edge_list.size(); ++i_edge) {
        const size_t edge_id     = cell_edge_list[i_edge];
        edge_cell_owner[edge_id] = std::min(edge_cell_owner[edge_id], diamond_descriptor.cell_number_vector[i_cell]);
      }
    }

    diamond_descriptor.edge_owner_vector.resize(edge_cell_owner.size());
    for (size_t i_edge = 0; i_edge < edge_cell_owner.size(); ++i_edge) {
      diamond_descriptor.edge_owner_vector[i_edge] = diamond_descriptor.cell_owner_vector[edge_cell_owner[i_edge]];
    }
  }

  m_connectivity = ConnectivityType::build(diamond_descriptor);

  if constexpr (Dimension == 1) {
    const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();

    NodeId dual_node_id            = 0;
    m_primal_node_to_dual_node_map = [&]() {
      std::vector<std::pair<NodeId, NodeId>> primal_node_to_dual_node_vector;
      for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
        if (node_to_cell_matrix[primal_node_id].size() == 1) {
          primal_node_to_dual_node_vector.push_back(std::make_pair(primal_node_id, dual_node_id++));
        }
      }
      return convert_to_array(primal_node_to_dual_node_vector);
    }();

    m_primal_cell_to_dual_node_map = [&]() {
      CellIdToNodeIdMap primal_cell_to_dual_node_map{primal_number_of_cells};
      for (CellId primal_cell_id = 0; primal_cell_id < primal_cell_to_dual_node_map.size(); ++primal_cell_id) {
        primal_cell_to_dual_node_map[primal_cell_id] = std::make_pair(primal_cell_id, dual_node_id++);
      }
      return primal_cell_to_dual_node_map;
    }();

    m_primal_face_to_dual_cell_map = [&]() {
      FaceIdToCellIdMap primal_face_to_dual_cell_map{primal_connectivity.numberOfFaces()};
      for (size_t id = 0; id < primal_face_to_dual_cell_map.size(); ++id) {
        const CellId dual_cell_id   = id;
        const FaceId primal_face_id = id;

        primal_face_to_dual_cell_map[id] = std::make_pair(primal_face_id, dual_cell_id);
      }
      return primal_face_to_dual_cell_map;
    }();
  }

  m_mapper =
    std::make_shared<ConnectivityToDiamondDualConnectivityDataMapper<Dimension>>(primal_connectivity,
                                                                                 dynamic_cast<const ConnectivityType&>(
                                                                                   *m_connectivity),
                                                                                 m_primal_node_to_dual_node_map,
                                                                                 m_primal_cell_to_dual_node_map,
                                                                                 m_primal_face_to_dual_cell_map);
}

DiamondDualConnectivityBuilder::DiamondDualConnectivityBuilder(const IConnectivity& connectivity)
{
  if (parallel::size() > 1) {
    throw NotImplementedError("Construction of diamond dual mesh is not implemented in parallel");
  }
  switch (connectivity.dimension()) {
  case 1: {
    this->_buildDiamondConnectivityFrom<1>(connectivity);
    break;
  }
  case 2: {
    this->_buildDiamondConnectivityFrom<2>(connectivity);
    break;
  }
  case 3: {
    this->_buildDiamondConnectivityFrom<3>(connectivity);
    break;
  }
  default: {
    throw UnexpectedError("invalid connectivity dimension: " + std::to_string(connectivity.dimension()));
  }
  }
}