ConnectivityComputer.cpp

#include <Connectivity.hpp>

#include <ConnectivityComputer.hpp>
#include <map>


template <typename ConnectivityType>
KOKKOS_INLINE_FUNCTION
ConnectivityMatrix ConnectivityComputer::
computeConnectivityMatrix(const ConnectivityType& connectivity,
                          const ItemType& item_type,
                          const ItemType& child_item_type) const
{
  ConnectivityMatrix item_to_child_item_matrix;
  if (connectivity.isConnectivityMatrixBuilt(child_item_type, item_type)) {
    const ConnectivityMatrix& child_to_item_matrix
        = connectivity.getMatrix(child_item_type, item_type);

    std::cout << "computing connectivity "
              << itemName(item_type) << " -> " << itemName(child_item_type) << '\n';

    item_to_child_item_matrix
        = this->_computeInverse(child_to_item_matrix);
  } else {
    std::cerr << "unable to compute connectivity "
              << itemName(item_type) << " -> " << itemName(child_item_type) << '\n';
    std::exit(0);
  }

  return item_to_child_item_matrix;
}

template
ConnectivityMatrix ConnectivityComputer::
computeConnectivityMatrix(const Connectivity1D&, const ItemType&, const ItemType&) const;

template
ConnectivityMatrix ConnectivityComputer::
computeConnectivityMatrix(const Connectivity2D&, const ItemType&, const ItemType&) const;


template
ConnectivityMatrix ConnectivityComputer::
computeConnectivityMatrix(const Connectivity3D&, const ItemType&, const ItemType&) const;

ConnectivityMatrix
ConnectivityComputer::
_computeInverse(const ConnectivityMatrix& item_to_child_matrix) const
{
  std::map<unsigned int, std::vector<unsigned int>> child_to_item_vector_map;
  const size_t& number_of_items = item_to_child_matrix.numRows();

  for (unsigned int j=0; j<number_of_items; ++j) {
    const auto& item_to_child = item_to_child_matrix.rowConst(j);
    for (unsigned int r=0; r<item_to_child.length; ++r) {
      child_to_item_vector_map[item_to_child(r)].push_back(j);
    }
  }

  {
    size_t i=0;
    for (const auto& [child_item_id, item_vector] : child_to_item_vector_map) {
      if (child_item_id != i) {
        std::cerr << "sparse item numerotation NIY\n";
        std::exit(0);
      }
      ++i;
    }
  }

  std::vector<std::vector<unsigned int>> child_to_items_vector(child_to_item_vector_map.size());
  for (const auto& [child_item_id, item_vector] : child_to_item_vector_map) {
    child_to_items_vector[child_item_id] = item_vector;
  }

  return ConnectivityMatrix(child_to_items_vector);
}

template <ItemType item_type,
          ItemType child_item_type,
          typename ConnectivityType>
SubItemValuePerItem<const unsigned short, child_item_type, item_type>
ConnectivityComputer::computeLocalItemNumberInChildItem(const ConnectivityType& connectivity) const
{
  const ConnectivityMatrix& child_item_to_items_matrix
      = connectivity.getMatrix(child_item_type, item_type);
  const ConnectivityMatrix& item_to_child_items_matrix
      = connectivity.getMatrix(item_type, child_item_type);

  SubItemValuePerItem<unsigned short, child_item_type, item_type> item_number_in_child_item(connectivity);
  for (unsigned int r=0; r<item_to_child_items_matrix.numRows(); ++r) {
    const auto& item_to_child_items = item_to_child_items_matrix.rowConst(r);
    for (unsigned short J=0; J<item_to_child_items.length; ++J) {
      const unsigned int j = item_to_child_items(J);
      const auto& child_item_to_items = child_item_to_items_matrix.rowConst(j);

      for (unsigned int R=0; R<child_item_to_items.length; ++R) {
        if (child_item_to_items(R) == r) {
          item_number_in_child_item(r,J) = R;
          break;
        }
      }
    }
  }

  return item_number_in_child_item;
}

// 1D

template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::node>
ConnectivityComputer::
computeLocalItemNumberInChildItem<ItemType::node,
                                  ItemType::cell>(const Connectivity1D&) const;

template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::face>
ConnectivityComputer::
computeLocalItemNumberInChildItem<ItemType::face,
                                  ItemType::cell>(const Connectivity1D&) const;

template SubItemValuePerItem<const unsigned short, ItemType::node, ItemType::cell>
ConnectivityComputer::
computeLocalItemNumberInChildItem<ItemType::cell,
                                  ItemType::node>(const Connectivity1D&) const;

// 2D

template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::node>
ConnectivityComputer::
computeLocalItemNumberInChildItem<ItemType::node,
                                  ItemType::cell>(const Connectivity2D&) const;

template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::face>
ConnectivityComputer::
computeLocalItemNumberInChildItem<ItemType::face,
                                  ItemType::cell>(const Connectivity2D&) const;

template SubItemValuePerItem<const unsigned short, ItemType::node, ItemType::cell>
ConnectivityComputer::
computeLocalItemNumberInChildItem<ItemType::cell,
                                  ItemType::node>(const Connectivity2D&) const;

// 3D

template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::node>
ConnectivityComputer::
computeLocalItemNumberInChildItem<ItemType::node,
                                  ItemType::cell>(const Connectivity3D&) const;

template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::face>
ConnectivityComputer::
computeLocalItemNumberInChildItem<ItemType::face,
                                  ItemType::cell>(const Connectivity3D&) const;

template SubItemValuePerItem<const unsigned short, ItemType::node, ItemType::cell>
ConnectivityComputer::
computeLocalItemNumberInChildItem<ItemType::cell,
                                  ItemType::node>(const Connectivity3D&) const;