Connectivity.cpp

#include <mesh/Connectivity.hpp>

#include <mesh/ConnectivityDescriptor.hpp>
#include <utils/Messenger.hpp>

#include <map>

template <size_t Dimension>
Connectivity<Dimension>::Connectivity()
{}

template <size_t Dimension>
void
Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor)
{
  Assert(descriptor.cell_to_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());
    Assert(descriptor.face_to_node_vector.size() == descriptor.face_number_vector.size());
    Assert(descriptor.face_owner_vector.size() == descriptor.face_number_vector.size());
  }

  auto& cell_to_node_matrix = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::node)];
  cell_to_node_matrix       = descriptor.cell_to_node_vector;

  {
    WeakCellValue<CellType> cell_type(*this);
    parallel_for(
      this->numberOfCells(), PUGS_LAMBDA(CellId j) { cell_type[j] = descriptor.cell_type_vector[j]; });
    m_cell_type = cell_type;
  }

  m_cell_number = WeakCellValue<int>(*this, convert_to_array(descriptor.cell_number_vector));

  Array node_number_array = convert_to_array(descriptor.node_number_vector);
  m_node_number           = WeakNodeValue<int>(*this, node_number_array);

  {
    WeakCellValue<int> cell_global_index(*this);
    int first_index = 0;
    parallel_for(
      this->numberOfCells(), PUGS_LAMBDA(CellId j) { cell_global_index[j] = first_index + j; });
    m_cell_global_index = cell_global_index;
  }

  m_cell_owner = WeakCellValue<int>(*this, convert_to_array(descriptor.cell_owner_vector));

  {
    const int rank = parallel::rank();
    WeakCellValue<bool> cell_is_owned(*this);
    parallel_for(
      this->numberOfCells(), PUGS_LAMBDA(CellId j) { cell_is_owned[j] = (m_cell_owner[j] == rank); });
    m_cell_is_owned = cell_is_owned;
  }

  Array node_owner_array = convert_to_array(descriptor.node_owner_vector);
  m_node_owner           = WeakNodeValue<int>{*this, node_owner_array};

  Array<bool> node_is_owned_array(this->numberOfNodes());
  {
    const int rank = parallel::rank();
    WeakNodeValue<bool> node_is_owned(*this, node_is_owned_array);
    parallel_for(
      this->numberOfNodes(), PUGS_LAMBDA(NodeId r) { node_is_owned[r] = (m_node_owner[r] == rank); });
    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) {
    // faces are similar to nodes
    m_face_number   = WeakFaceValue<int>(*this, node_number_array);
    m_face_owner    = WeakFaceValue<int>(*this, node_owner_array);
    m_face_is_owned = WeakFaceValue<bool>(*this, node_is_owned_array);

    // edges are similar to nodes
    m_edge_number   = WeakEdgeValue<int>(*this, node_number_array);
    m_edge_owner    = WeakEdgeValue<int>(*this, node_owner_array);
    m_edge_is_owned = WeakEdgeValue<bool>(*this, node_is_owned_array);
  } else {
    m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)] = descriptor.face_to_node_vector;

    m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)] = descriptor.cell_to_face_vector;

    {
      FaceValuePerCell<bool> cell_face_is_reversed(*this);
      for (CellId j = 0; j < descriptor.cell_face_is_reversed_vector.size(); ++j) {
        const auto& face_cells_vector = descriptor.cell_face_is_reversed_vector[j];
        for (unsigned short lj = 0; lj < face_cells_vector.size(); ++lj) {
          cell_face_is_reversed(j, lj) = face_cells_vector[lj];
        }
      }
      m_cell_face_is_reversed = cell_face_is_reversed;
    }

    Array face_number_array = convert_to_array(descriptor.face_number_vector);
    m_face_number           = WeakFaceValue<int>(*this, face_number_array);

    Array face_owner_array = convert_to_array(descriptor.face_owner_vector);
    m_face_owner           = WeakFaceValue<int>(*this, face_owner_array);

    Array<bool> face_is_owned_array(this->numberOfFaces());
    {
      const int rank = parallel::rank();
      WeakFaceValue<bool> face_is_owned(*this, face_is_owned_array);
      parallel_for(
        this->numberOfFaces(), PUGS_LAMBDA(FaceId l) { face_is_owned[l] = (m_face_owner[l] == rank); });
      m_face_is_owned = face_is_owned;
    }

    m_ref_face_list_vector = descriptor.template refItemListVector<ItemType::face>();

    if constexpr (Dimension == 2) {
      // edges are similar to faces
      m_edge_number   = WeakEdgeValue<int>(*this, face_number_array);
      m_edge_owner    = WeakEdgeValue<int>(*this, face_owner_array);
      m_edge_is_owned = WeakEdgeValue<bool>(*this, face_is_owned_array);

    } else {
      m_item_to_item_matrix[itemTId(ItemType::edge)][itemTId(ItemType::node)] = descriptor.edge_to_node_vector;

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

      m_edge_number = WeakEdgeValue<int>(*this, convert_to_array(descriptor.edge_number_vector));
      m_edge_owner  = WeakEdgeValue<int>(*this, convert_to_array(descriptor.edge_owner_vector));

      {
        const int rank = parallel::rank();
        WeakEdgeValue<bool> edge_is_owned(*this);
        parallel_for(
          this->numberOfEdges(), PUGS_LAMBDA(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>();
    }
  }
}

template Connectivity<1>::Connectivity();
template Connectivity<2>::Connectivity();
template Connectivity<3>::Connectivity();

template void Connectivity<1>::_buildFrom(const ConnectivityDescriptor&);
template void Connectivity<2>::_buildFrom(const ConnectivityDescriptor&);
template void Connectivity<3>::_buildFrom(const ConnectivityDescriptor&);