Topology.h

// Copyright (C) 2006-2024 Anders Logg and Garth N. Wells
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier:    LGPL-3.0-or-later
 
#pragma once
 
#include <array>
#include <concepts>
#include <cstdint>
#include <dolfinx/common/MPI.h>
#include <dolfinx/graph/AdjacencyList.h>
#include <map>
#include <memory>
#include <optional>
#include <span>
#include <thread>
#include <tuple>
#include <utility>
#include <vector>
 
namespace dolfinx::common
{
class IndexMap;
}
 
namespace dolfinx::mesh
{
template <typename R>
concept CellRange = std::ranges::input_range<R> and std::ranges::sized_range<R>
                    and std::is_integral_v<
                        std::remove_const_t<std::ranges::range_value_t<R>>>;
 
enum class CellType : std::int8_t;

class Topology
{
public:
  Topology(
      std::vector<CellType> cell_types,
      std::shared_ptr<const common::IndexMap> vertex_map,
      std::vector<std::shared_ptr<const common::IndexMap>> cell_maps,
      std::vector<std::shared_ptr<graph::AdjacencyList<std::int32_t>>> cells,
      const std::optional<std::vector<std::vector<std::int64_t>>>&
          original_cell_index = std::nullopt,
      int num_threads = 1);

  Topology(const Topology& topology) = default;

  Topology(Topology&& topology) = default;

  ~Topology() = default;

  Topology& operator=(const Topology& topology) = delete;

  Topology& operator=(Topology&& topology) = default;

  int dim() const noexcept;

  const std::vector<CellType>& entity_types(int dim) const;

  CellType cell_type() const;

  std::vector<std::shared_ptr<const common::IndexMap>>
  index_maps(int dim) const;

  std::shared_ptr<const common::IndexMap> index_map(int dim) const;

  std::shared_ptr<const graph::AdjacencyList<std::int32_t>>
  connectivity(std::array<int, 2> d0, std::array<int, 2> d1) const;

  std::shared_ptr<const graph::AdjacencyList<std::int32_t>>
  connectivity(int d0, int d1) const;

  const std::vector<std::uint32_t>& get_cell_permutation_info() const;

  const std::vector<std::uint8_t>& get_facet_permutations() const;

  std::vector<CellType> cell_types() const;

  const std::vector<std::int32_t>& interprocess_facets(int index) const;

  const std::vector<std::int32_t>& interprocess_facets() const;

  bool create_entities(int dim, int num_threads = 1);

  void create_connectivity(int d0, int d1);

  void create_entity_permutations();

  std::vector<std::vector<std::int64_t>> original_cell_index;

  MPI_Comm comm() const;
 
private:
  // Cell types for entities in Topology, where _entity_types_new[d][i]
  // is the ith entity type of dimension d
  std::vector<std::vector<CellType>> _entity_types;
 
  // Parallel layout of entities for each dimension and cell type
  // flattened in the same layout as _entity_types above.
  // std::vector<std::shared_ptr<const common::IndexMap>> _index_map;
 
  // _index_maps[(d, i) is the index map for the ith entity type of
  // dimension d
  std::map<std::array<int, 2>, std::shared_ptr<const common::IndexMap>>
      _index_maps;
 
  // Connectivity between cell types _connectivity_new[(dim0, i0),
  // (dim1, i1)] is the connection from (dim0, i0) -> (dim1, i1),
  // where dim0 and dim1 are topological dimensions and i0 and i1
  // are the indices of cell types (following the order in _entity_types).
  std::map<std::pair<std::array<int, 2>, std::array<int, 2>>,
           std::shared_ptr<graph::AdjacencyList<std::int32_t>>>
      _connectivity;
 
  // The facet permutations (local facet, cell))
  // [cell0_0, cell0_1, ,cell0_2, cell1_0, cell1_1, ,cell1_2, ...,
  // celln_0, celln_1, ,celln_2,]
  std::vector<std::uint8_t> _facet_permutations;
 
  // Cell permutation info. See the documentation for
  // get_cell_permutation_info for documentation of how this is encoded.
  std::vector<std::uint32_t> _cell_permutations;
 
  // List of facets that are on the inter-process boundary for each
  // facet type. _interprocess_facets[i] is the inter-process facets of
  // facet type i.
  std::vector<std::vector<std::int32_t>> _interprocess_facets;
};

Topology
create_topology(MPI_Comm comm, const std::vector<CellType>& cell_types,
                std::vector<std::span<const std::int64_t>> cells,
                std::vector<std::span<const std::int64_t>> original_cell_index,
                std::vector<std::span<const int>> ghost_owners,
                std::span<const std::int64_t> boundary_vertices,
                int num_threads);

Topology create_topology(MPI_Comm comm, std::span<const std::int64_t> cells,
                         std::span<const std::int64_t> original_cell_index,
                         std::span<const int> ghost_owners, CellType cell_type,
                         std::span<const std::int64_t> boundary_vertices,
                         int num_threads);

std::tuple<Topology, std::vector<int32_t>, std::vector<int32_t>>
create_subtopology(const Topology& topology, int dim,
                   std::span<const std::int32_t> entities);

std::vector<std::int32_t>
entities_to_index(const Topology& topology, int dim,
                  std::span<const std::int32_t> entities);

std::vector<std::vector<std::int32_t>>
compute_mixed_cell_pairs(const Topology& topology, mesh::CellType facet_type);
 
} // namespace dolfinx::mesh