AdjacencyList.h

// Copyright (C) 2019-2021 Garth N. Wells
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier:    LGPL-3.0-or-later
 
#pragma once
 
#include <cassert>
#include <numeric>
#include <sstream>
#include <utility>
#include <vector>
#include <xtensor/xtensor.hpp>
#include <xtl/xspan.hpp>
 
namespace dolfinx::graph
{
 
template <typename T>
auto create_adjacency_data(const xt::xtensor<T, 2>& array)
{
  std::vector<T> data(array.shape(0) * array.shape(1));
  std::vector<std::int32_t> offset(array.shape(0) + 1, 0);
  for (std::size_t i = 0; i < array.shape(0); ++i)
  {
    for (std::size_t j = 0; j < array.shape(1); ++j)
      data[i * array.shape(1) + j] = array(i, j);
    offset[i + 1] = offset[i] + array.shape(1);
  }
  return std::pair(std::move(data), std::move(offset));
}
 
template <typename T>
class AdjacencyList
{
public:
  explicit AdjacencyList(const std::int32_t n) : _array(n), _offsets(n + 1)
  {
    std::iota(_array.begin(), _array.end(), 0);
    std::iota(_offsets.begin(), _offsets.end(), 0);
  }
 
  template <
      typename U, typename V,
      typename = std::enable_if_t<
          std::is_same<std::vector<T>, std::decay_t<U>>::value
          && std::is_same<std::vector<std::int32_t>, std::decay_t<V>>::value>>
  AdjacencyList(U&& data, V&& offsets)
      : _array(std::forward<U>(data)), _offsets(std::forward<V>(offsets))
  {
    _array.reserve(_offsets.back());
    assert(_offsets.back() == (std::int32_t)_array.size());
  }
 
  template <typename X>
  explicit AdjacencyList(const std::vector<X>& data)
  {
    // Initialize offsets and compute total size
    _offsets.reserve(data.size() + 1);
    _offsets.push_back(0);
    for (auto row = data.begin(); row != data.end(); ++row)
      _offsets.push_back(_offsets.back() + row->size());
 
    _array.reserve(_offsets.back());
    for (auto e = data.begin(); e != data.end(); ++e)
      _array.insert(_array.end(), e->begin(), e->end());
  }
 
  AdjacencyList(const AdjacencyList& list) = default;
 
  AdjacencyList(AdjacencyList&& list) = default;
 
  ~AdjacencyList() = default;
 
  AdjacencyList& operator=(const AdjacencyList& list) = default;
 
  AdjacencyList& operator=(AdjacencyList&& list) = default;
 
  bool operator==(const AdjacencyList& list) const
  {
    return this->_array == list._array and this->_offsets == list._offsets;
  }
 
  std::int32_t num_nodes() const { return _offsets.size() - 1; }
 
  int num_links(int node) const
  {
    assert((node + 1) < (int)_offsets.size());
    return _offsets[node + 1] - _offsets[node];
  }
 
  xtl::span<T> links(int node)
  {
    return xtl::span<T>(_array.data() + _offsets[node],
                        _offsets[node + 1] - _offsets[node]);
  }
 
  xtl::span<const T> links(int node) const
  {
    return xtl::span<const T>(_array.data() + _offsets[node],
                              _offsets[node + 1] - _offsets[node]);
  }
 
  const std::vector<T>& array() const { return _array; }
 
  std::vector<T>& array() { return _array; }
 
  const std::vector<std::int32_t>& offsets() const { return _offsets; }
 
  std::vector<std::int32_t>& offsets() { return _offsets; }
 
  std::string str() const
  {
    std::stringstream s;
    s << "<AdjacencyList> with " + std::to_string(this->num_nodes()) + " nodes"
      << std::endl;
    for (std::size_t e = 0; e < _offsets.size() - 1; ++e)
    {
      s << "  " << e << ": [";
      for (auto link : this->links(e))
        s << link << " ";
      s << "]" << std::endl;
    }
    return s.str();
  }
 
private:
  // Connections for all entities stored as a contiguous array
  std::vector<T> _array;
 
  // Position of first connection for each entity (using local index)
  std::vector<std::int32_t> _offsets;
};
 
template <typename U>
AdjacencyList<typename U::value_type> regular_adjacency_list(U&& data,
                                                             int degree)
{
  if (degree == 0 and !data.empty())
  {
    throw std::runtime_error("Degree is zero but data is not empty for "
                             "constant degree AdjacencyList");
  }
 
  if (degree > 0 and data.size() % degree != 0)
  {
    throw std::runtime_error(
        "Incompatible data size and degree for constant degree AdjacencyList");
  }
 
  std::int32_t num_nodes = degree == 0 ? data.size() : data.size() / degree;
  std::vector<std::int32_t> offsets(num_nodes + 1, 0);
  for (std::size_t i = 1; i < offsets.size(); ++i)
    offsets[i] = offsets[i - 1] + degree;
  return AdjacencyList<typename U::value_type>(std::forward<U>(data),
                                               std::move(offsets));
}
 
} // namespace dolfinx::graph