diff --git a/cpp/dolfinx/mesh/Topology.cpp b/cpp/dolfinx/mesh/Topology.cpp
index 37fa3ad3781..fa7f27fefb8 100644
--- a/cpp/dolfinx/mesh/Topology.cpp
+++ b/cpp/dolfinx/mesh/Topology.cpp
@@ -1,4 +1,4 @@
-// Copyright (C) 2006-2022 Anders Logg and Garth N. Wells
+// Copyright (C) 2006-2024 Anders Logg and Garth N. Wells
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -43,8 +43,6 @@ determine_sharing_ranks(MPI_Comm comm, std::span<const std::int64_t> indices)
 {
   common::Timer timer("Topology: determine shared index ownership");
 
-  const int size = dolfinx::MPI::size(comm);
-
   // FIXME: use sensible name
   std::int64_t global_range = 0;
   {
@@ -57,6 +55,7 @@ determine_sharing_ranks(MPI_Comm comm, std::span<const std::int64_t> indices)
   // Build {dest, pos} list, and sort
   std::vector<std::array<int, 2>> dest_to_index;
   {
+    const int size = dolfinx::MPI::size(comm);
     dest_to_index.reserve(indices.size());
     for (auto idx : indices)
     {
@@ -136,10 +135,8 @@ determine_sharing_ranks(MPI_Comm comm, std::span<const std::int64_t> indices)
   // Build {global index, pos, src} list
   std::vector<std::array<std::int64_t, 3>> indices_list;
   for (std::size_t p = 0; p < recv_disp0.size() - 1; ++p)
-  {
     for (std::int32_t i = recv_disp0[p]; i < recv_disp0[p + 1]; ++i)
       indices_list.push_back({recv_buffer0[i], i, int(p)});
-  }
   std::ranges::sort(indices_list);
 
   // Find which ranks have each index
@@ -173,7 +170,7 @@ determine_sharing_ranks(MPI_Comm comm, std::span<const std::int64_t> indices)
       // Update number of items to be sent to each rank and record owner
       for (auto itx = it; itx != it1; ++itx)
       {
-        auto& data = *itx;
+        const std::array<std::int64_t, 3>& data = *itx;
         num_items_per_pos1[data[1]] = num + 1;
         num_items_per_dest1[data[2]] += num + 1;
       }
@@ -380,8 +377,7 @@ exchange_indexing(MPI_Comm comm, std::span<const std::int64_t> indices,
   std::vector<int> src, dest;
   for (std::int32_t i = 0; i < index_to_ranks.num_nodes(); ++i)
   {
-    auto ranks = index_to_ranks.links(i);
-    if (ranks.front() == mpi_rank)
+    if (auto ranks = index_to_ranks.links(i); ranks.front() == mpi_rank)
       dest.insert(dest.end(), std::next(ranks.begin()), ranks.end());
     else
       src.push_back(ranks.front());
@@ -406,8 +402,7 @@ exchange_indexing(MPI_Comm comm, std::span<const std::int64_t> indices,
   {
     // Get (global) ranks that share this vertex. Note that first rank
     // is the owner.
-    auto ranks = index_to_ranks.links(i);
-    if (ranks.front() == mpi_rank)
+    if (auto ranks = index_to_ranks.links(i); ranks.front() == mpi_rank)
     {
       // Get local vertex index
       std::int64_t idx_old = indices[i];
@@ -452,7 +447,7 @@ exchange_indexing(MPI_Comm comm, std::span<const std::int64_t> indices,
 
     std::vector<std::int64_t> sbuffer;
     sbuffer.reserve(send_disp.back());
-    for (auto& data : send_buffer)
+    for (const std::vector<std::int64_t>& data : send_buffer)
       sbuffer.insert(sbuffer.end(), data.begin(), data.end());
 
     // Get receive sizes
@@ -529,9 +524,9 @@ std::vector<std::array<std::int64_t, 3>> exchange_ghost_indexing(
 
   // --
 
-  // For each rank, list of owned vertices that are ghosted by other ranks
+  // For each rank, list of owned vertices that are ghosted by other
+  // ranks
   std::vector<std::vector<std::int64_t>> shared_vertices_fwd(dest.size());
-
   {
     // -- Send cell ghost indices to owner
     MPI_Comm comm1;
@@ -557,14 +552,17 @@ std::vector<std::array<std::int64_t, 3>> exchange_ghost_indexing(
 
     // Compute send sizes and displacements
     std::vector<int> send_sizes, send_disp{0};
-    auto it = owner_to_ghost.begin();
-    while (it != owner_to_ghost.end())
     {
-      auto it1 = std::find_if(it, owner_to_ghost.end(),
-                              [r = it->first](auto x) { return x.first != r; });
-      send_sizes.push_back(std::distance(it, it1));
-      send_disp.push_back(send_disp.back() + send_sizes.back());
-      it = it1;
+      auto it = owner_to_ghost.begin();
+      while (it != owner_to_ghost.end())
+      {
+        auto it1
+            = std::find_if(it, owner_to_ghost.end(),
+                           [r = it->first](auto x) { return x.first != r; });
+        send_sizes.push_back(std::distance(it, it1));
+        send_disp.push_back(send_disp.back() + send_sizes.back());
+        it = it1;
+      }
     }
 
     // Exchange number of indices to send/receive from each rank
@@ -588,7 +586,7 @@ std::vector<std::array<std::int64_t, 3>> exchange_ghost_indexing(
     MPI_Comm_free(&comm1);
 
     // Iterate over ranks that ghost cells owned by this rank
-    auto local_range = map0.local_range();
+    std::array<std::int64_t, 2> local_range = map0.local_range();
     for (std::size_t r = 0; r < recv_disp.size() - 1; ++r)
     {
       assert(r < shared_vertices_fwd.size());
@@ -632,10 +630,10 @@ std::vector<std::array<std::int64_t, 3>> exchange_ghost_indexing(
     const int mpi_rank = dolfinx::MPI::rank(comm);
 
     // Iterate over each rank to send vertex data to
-    for (const auto& vertices_old : shared_vertices_fwd)
+    for (const std::vector<std::int64_t>& vertices_old : shared_vertices_fwd)
     {
       // Iterate over vertex indices (old) for current destination rank
-      for (auto vertex_old : vertices_old)
+      for (std::int64_t vertex_old : vertices_old)
       {
         // Find new vertex index and determine owning rank
         auto it = std::ranges::lower_bound(
@@ -645,7 +643,6 @@ std::vector<std::array<std::int64_t, 3>> exchange_ghost_indexing(
         assert(it != global_local_entities1.end());
         assert(it->first == vertex_old);
         assert(it->second != -1);
-
         std::int64_t global_idx = it->second < nlocal1
                                       ? it->second + offset1
                                       : ghost_entities1[it->second - nlocal1];
@@ -713,129 +710,92 @@ std::vector<std::int32_t> convert_to_local_indexing(
 
 //-----------------------------------------------------------------------------
 Topology::Topology(
-    MPI_Comm comm, CellType cell_type,
-    std::shared_ptr<const common::IndexMap> vertex_map,
-    std::shared_ptr<const common::IndexMap> cell_map,
-    std::shared_ptr<graph::AdjacencyList<std::int32_t>> cells,
-    const std::optional<std::vector<std::int64_t>>& original_index)
-    : Topology(comm, {cell_type}, vertex_map, {cell_map}, {cells},
-               original_index
-                   ? std::vector<std::vector<std::int64_t>>{*original_index}
-                   : std::optional<std::vector<std::vector<std::int64_t>>>(
-                         std::nullopt))
-{
-}
-//-----------------------------------------------------------------------------
-Topology::Topology(
-    MPI_Comm comm, std::vector<CellType> cell_types,
+    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_index)
     : original_cell_index(original_index
                               ? *original_index
-                              : std::vector<std::vector<std::int64_t>>()),
-      _comm(comm), _entity_types({mesh::CellType::point}),
-      _entity_type_offsets({0, 1}), _interprocess_facets(1)
+                              : std::vector<std::vector<std::int64_t>>())
 {
   assert(!cell_types.empty());
   int tdim = cell_dim(cell_types.front());
-
 #ifndef NDEBUG
   for (auto ct : cell_types)
     assert(cell_dim(ct) == tdim);
 #endif
 
-  // Create all the entity types in the mesh
-  if (tdim > 1)
-  {
-    // In 2D, the facet is an interval
-    _entity_types.push_back(CellType::interval);
-    _entity_type_offsets.push_back(_entity_types.size());
+  _entity_types.resize(tdim + 1);
+  _entity_types[0] = {mesh::CellType::point};
+  _entity_types[tdim] = cell_types;
 
-    if (tdim == 3)
+  // Set data
+  _index_maps.insert({{0, 0}, vertex_map});
+  _connectivity.insert(
+      {{{0, 0}, {0, 0}},
+       std::make_shared<graph::AdjacencyList<std::int32_t>>(
+           vertex_map->size_local() + vertex_map->num_ghosts())});
+  if (tdim > 0)
+  {
+    for (std::size_t i = 0; i < cell_types.size(); ++i)
     {
-      // Find all facet types
-      std::set<mesh::CellType> facet_types;
-      for (auto c : cell_types)
-      {
-        assert(cell_dim(c) == tdim);
-        for (int i = 0; i < cell_num_entities(c, tdim - 1); ++i)
-          facet_types.insert(cell_facet_type(c, i));
-      }
-      _entity_types.insert(_entity_types.end(), facet_types.begin(),
-                           facet_types.end());
-      _interprocess_facets.resize(facet_types.size());
-      _entity_type_offsets.push_back(_entity_types.size());
+      _index_maps.insert({{tdim, (int)i}, cell_maps[i]});
+      _connectivity.insert({{{tdim, int(i)}, {0, 0}}, cells[i]});
     }
   }
 
-  // Cell Types
-  _entity_types.insert(_entity_types.end(), cell_types.begin(),
-                       cell_types.end());
-  if (tdim > 0)
-    _entity_type_offsets.push_back(_entity_types.size());
-
-  int conn_size = _entity_type_offsets.back();
-  _index_map.resize(conn_size);
-
-  // Create square list of lists
-  _connectivity.resize(conn_size);
-  for (auto& c : _connectivity)
-    c.resize(conn_size);
-
-  // Set data
-  this->set_index_map(0, vertex_map);
-  this->set_connectivity(
-      std::make_shared<graph::AdjacencyList<std::int32_t>>(
-          vertex_map->size_local() + vertex_map->num_ghosts()),
-      0, 0);
-  for (std::size_t i = 0; i < cell_types.size(); ++i)
+  // FIXME: This is a hack for setting _interprocess_facets when
+  // tdim==1, i.e. the 'facets' are vertices
+  if (tdim == 1)
   {
-    this->set_index_map(tdim, i, cell_maps[i]);
-    this->set_connectivity(cells[i], {tdim, int(i)}, {0, 0});
+    auto [cell_entity, entity_vertex, index_map, interprocess_entities]
+        = compute_entities(*this, 0, CellType::point);
+    std::ranges::sort(interprocess_entities);
+    _interprocess_facets.push_back(std::move(interprocess_entities));
   }
 }
 //-----------------------------------------------------------------------------
-int Topology::dim() const noexcept { return _entity_type_offsets.size() - 2; }
+int Topology::dim() const noexcept
+{
+  return mesh::cell_dim(_entity_types.back().front());
+}
 //-----------------------------------------------------------------------------
-std::vector<CellType> Topology::entity_types(int dim) const
+const std::vector<CellType>& Topology::entity_types(int dim) const
 {
-  assert(dim < (int)_entity_type_offsets.size() - 1 and dim >= 0);
-  return std::vector<CellType>(
-      std::next(_entity_types.begin(), _entity_type_offsets[dim]),
-      std::next(_entity_types.begin(), _entity_type_offsets[dim + 1]));
+  return _entity_types.at(dim);
 }
 //-----------------------------------------------------------------------------
-mesh::CellType Topology::cell_type() const { return _entity_types.back(); }
+mesh::CellType Topology::cell_type() const
+{
+  return _entity_types.back().at(0);
+}
 //-----------------------------------------------------------------------------
 std::vector<std::shared_ptr<const common::IndexMap>>
 Topology::index_maps(int dim) const
 {
-  assert(dim < (int)_entity_type_offsets.size() - 1);
-  std::vector maps(
-      std::next(_index_map.begin(), _entity_type_offsets[dim]),
-      std::next(_index_map.begin(), _entity_type_offsets[dim + 1]));
+  std::vector<std::shared_ptr<const common::IndexMap>> maps;
+  for (std::size_t i = 0; i < _entity_types[dim].size(); ++i)
+  {
+    auto it = _index_maps.find({dim, int(i)});
+    assert(it != _index_maps.end());
+    maps.push_back(it->second);
+  }
   return maps;
 }
 //-----------------------------------------------------------------------------
 std::shared_ptr<const common::IndexMap> Topology::index_map(int dim) const
 {
-  assert(dim < (int)_entity_type_offsets.size() - 1);
-  return _index_map[_entity_type_offsets[dim]];
+  return this->index_maps(dim).at(0);
 }
 //-----------------------------------------------------------------------------
 std::shared_ptr<const graph::AdjacencyList<std::int32_t>>
 Topology::connectivity(std::array<int, 2> d0, std::array<int, 2> d1) const
 {
-  int dim0 = d0[0];
-  int dim1 = d1[0];
-  assert(dim0 < (int)_entity_type_offsets.size() - 1);
-  assert(d0[1] < (_entity_type_offsets[dim0 + 1] - _entity_type_offsets[dim0]));
-  assert(dim1 < (int)_entity_type_offsets.size() - 1);
-  assert(d1[1] < (_entity_type_offsets[dim1 + 1] - _entity_type_offsets[dim1]));
-  return _connectivity[_entity_type_offsets[dim0] + d0[1]]
-                      [_entity_type_offsets[dim1] + d1[1]];
+  if (auto it = _connectivity.find({d0, d1}); it == _connectivity.end())
+    return nullptr;
+  else
+    return it->second;
 }
 //-----------------------------------------------------------------------------
 std::shared_ptr<const graph::AdjacencyList<std::int32_t>>
@@ -875,6 +835,8 @@ const std::vector<std::uint8_t>& Topology::get_facet_permutations() const
 //-----------------------------------------------------------------------------
 const std::vector<std::int32_t>& Topology::interprocess_facets(int index) const
 {
+  if (_interprocess_facets.empty())
+    throw std::runtime_error("Interprocess facets have not been computed.");
   return _interprocess_facets.at(index);
 }
 //-----------------------------------------------------------------------------
@@ -883,84 +845,64 @@ const std::vector<std::int32_t>& Topology::interprocess_facets() const
   return this->interprocess_facets(0);
 }
 //-----------------------------------------------------------------------------
-void Topology::set_index_map(int dim, int i,
-                             std::shared_ptr<const common::IndexMap> map)
-{
-  assert(dim < (int)_entity_type_offsets.size() - 1);
-  assert(i < (_entity_type_offsets[dim + 1] - _entity_type_offsets[dim]));
-  _index_map[_entity_type_offsets[dim] + i] = map;
-}
-//-----------------------------------------------------------------------------
-void Topology::set_index_map(int dim,
-                             std::shared_ptr<const common::IndexMap> map)
-{
-  if (_entity_type_offsets[dim + 1] - _entity_type_offsets[dim] != 1)
-    throw std::runtime_error("Cannot set IndexMap on mixed topology mesh");
-  this->set_index_map(dim, 0, map);
-}
-//-----------------------------------------------------------------------------
-void Topology::set_connectivity(
-    std::shared_ptr<graph::AdjacencyList<std::int32_t>> c,
-    std::array<int, 2> d0, std::array<int, 2> d1)
-{
-  auto [dim0, i0] = d0;
-  auto [dim1, i1] = d1;
-  assert(dim0 < (int)_entity_type_offsets.size() - 1);
-  assert(i0 < (_entity_type_offsets[dim0 + 1] - _entity_type_offsets[dim0]));
-  assert(dim1 < (int)_entity_type_offsets.size() - 1);
-  assert(i1 < (_entity_type_offsets[dim1 + 1] - _entity_type_offsets[dim1]));
-  _connectivity[_entity_type_offsets[dim0] + i0]
-               [_entity_type_offsets[dim1] + i1]
-      = c;
-}
-//-----------------------------------------------------------------------------
-void Topology::set_connectivity(
-    std::shared_ptr<graph::AdjacencyList<std::int32_t>> c, int d0, int d1)
-{
-  this->set_connectivity(c, {d0, 0}, {d1, 0});
-}
-//-----------------------------------------------------------------------------
-std::int32_t Topology::create_entities(int dim)
+bool Topology::create_entities(int dim)
 {
   // TODO: is this check sufficient/correct? Does not catch the
   // cell_entity entity case. Should there also be a check for
   // connectivity(this->dim(), dim)?
+
   // Skip if already computed (vertices (dim=0) should always exist)
   if (connectivity(dim, 0))
-    return -1;
+    return false;
+
+  int tdim = this->dim();
+  if (dim == 1 and tdim > 1)
+    _entity_types[1] = {mesh::CellType::interval};
+  else if (dim == 2 and tdim > 2)
+  {
+    //  Find all facet types
+    std::set<mesh::CellType> e_types;
+    for (auto c : _entity_types[tdim])
+      for (int i = 0; i < cell_num_entities(c, dim); ++i)
+        e_types.insert(cell_facet_type(c, i));
+    _entity_types[dim] = std::vector(e_types.begin(), e_types.end());
+  }
 
-  for (std::size_t index = 0; index < this->entity_types(dim).size(); ++index)
+  // for (std::size_t index = 0; index < this->entity_types(dim).size();
+  // ++index)
+  for (auto entity = this->entity_types(dim).begin();
+       entity != this->entity_types(dim).end(); ++entity)
   {
+    int index = std::distance(this->entity_types(dim).begin(), entity);
+
     // Create local entities
     auto [cell_entity, entity_vertex, index_map, interprocess_entities]
-        = compute_entities(_comm.comm(), *this, dim, index);
+        = compute_entities(*this, dim, *entity);
     for (std::size_t k = 0; k < cell_entity.size(); ++k)
     {
       if (cell_entity[k])
       {
-        set_connectivity(cell_entity[k], {this->dim(), int(k)},
-                         {dim, int(index)});
+        _connectivity.insert(
+            {{{this->dim(), int(k)}, {dim, int(index)}}, cell_entity[k]});
       }
     }
 
     // TODO: is this check necessary? Seems redundant after the "skip
     // check"
     if (entity_vertex)
-      set_connectivity(entity_vertex, {dim, int(index)}, {0, 0});
+      _connectivity.insert({{{dim, int(index)}, {0, 0}}, entity_vertex});
 
-    assert(index_map);
-    this->set_index_map(dim, index, index_map);
+    _index_maps.insert({{dim, int(index)}, index_map});
 
-    // Store boundary facets
+    // Store interprocess facets
     if (dim == this->dim() - 1)
     {
       std::ranges::sort(interprocess_entities);
-      assert(index < _interprocess_facets.size());
-      _interprocess_facets[index] = std::move(interprocess_entities);
+      _interprocess_facets.push_back(std::move(interprocess_entities));
     }
   }
 
-  return this->index_maps(dim)[0]->size_local();
+  return true;
 }
 //-----------------------------------------------------------------------------
 void Topology::create_connectivity(int d0, int d1)
@@ -995,9 +937,10 @@ void Topology::create_connectivity(int d0, int d1)
 
       // Attach connectivities
       if (c_d0_d1)
-        set_connectivity(c_d0_d1, {d0, i0}, {d1, i1});
+        _connectivity.insert({{{d0, i0}, {d1, i1}}, c_d0_d1});
+
       if (c_d1_d0)
-        set_connectivity(c_d1_d0, {d1, i1}, {d0, i0});
+        _connectivity.insert({{{d1, i1}, {d0, i0}}, c_d1_d0});
     }
   }
 }
@@ -1007,11 +950,12 @@ void Topology::create_entity_permutations()
   if (!_cell_permutations.empty())
     return;
 
-  int tdim = this->dim();
+  // FIXME: Is creating all entities always required? Could it be made
+  // cheaper by doing a local version? This call does quite a lot of
+  // parallel work.
 
-  // FIXME: Is this always required? Could it be made cheaper by doing a
-  // local version? This call does quite a lot of parallel work
   // Create all mesh entities
+  int tdim = this->dim();
   for (int d = 0; d < tdim; ++d)
     create_entities(d);
 
@@ -1021,7 +965,12 @@ void Topology::create_entity_permutations()
   _cell_permutations = std::move(cell_permutations);
 }
 //-----------------------------------------------------------------------------
-MPI_Comm Topology::comm() const { return _comm.comm(); }
+MPI_Comm Topology::comm() const
+{
+  auto it = _index_maps.find({this->dim(), 0});
+  assert(it != _index_maps.end());
+  return it->second->comm();
+}
 //-----------------------------------------------------------------------------
 Topology mesh::create_topology(
     MPI_Comm comm, const std::vector<CellType>& cell_types,
@@ -1036,8 +985,9 @@ Topology mesh::create_topology(
   assert(ghost_owners.size() == cells.size());
   assert(original_cell_index.size() == cells.size());
 
+  // Check cell data consistency and compile spans of owned and ghost
+  // cells
   spdlog::info("Create topology (generalised)");
-  // Check cell data consistency and compile spans of owned and ghost cells
   std::vector<std::int32_t> num_local_cells(cell_types.size());
   std::vector<std::span<const std::int64_t>> owned_cells;
   std::vector<std::span<const std::int64_t>> ghost_cells;
@@ -1103,8 +1053,8 @@ Topology mesh::create_topology(
     {
       for (auto vtx : cells[i])
       {
-        auto it = std::ranges::lower_bound(owned_vertices, vtx);
-        if (it != owned_vertices.end() and *it == vtx)
+        if (auto it = std::ranges::lower_bound(owned_vertices, vtx);
+            it != owned_vertices.end() and *it == vtx)
         {
           std::size_t pos = std::distance(owned_vertices.begin(), it);
           if (local_vertex_indices[pos] < 0)
@@ -1185,8 +1135,7 @@ Topology mesh::create_topology(
       // ghost vertices that are not on the process boundary. Data is
       // communicated via ghost cells. Note that the ghost cell owner
       // (who we get the vertex index from) is not necessarily the
-      // vertex owner.
-      // Repeat for each cell type
+      // vertex owner. Repeat for each cell type,
       std::vector<std::array<std::int64_t, 3>> recv_data;
       for (std::size_t i = 0; i < cell_types.size(); ++i)
       {
@@ -1290,11 +1239,11 @@ Topology mesh::create_topology(
 
   // Save original cell index
   std::vector<std::vector<std::int64_t>> orig_index;
-  for (std::span<const std::int64_t> idx : original_cell_index)
-    orig_index.emplace_back(idx.begin(), idx.end());
+  std::transform(original_cell_index.begin(), original_cell_index.end(),
+                 std::back_inserter(orig_index), [](auto idx)
+                 { return std::vector<std::int64_t>(idx.begin(), idx.end()); });
 
-  return Topology(comm, cell_types, index_map_v, index_map_c, cells_c,
-                  orig_index);
+  return Topology(cell_types, index_map_v, index_map_c, cells_c, orig_index);
 }
 //-----------------------------------------------------------------------------
 Topology
@@ -1387,7 +1336,7 @@ mesh::create_subtopology(const Topology& topology, int dim,
   auto sub_e_to_v = std::make_shared<graph::AdjacencyList<std::int32_t>>(
       std::move(sub_e_to_v_vec), std::move(sub_e_to_v_offsets));
 
-  return {Topology(topology.comm(), entity_type, submap0, submap, sub_e_to_v),
+  return {Topology({entity_type}, submap0, {submap}, {sub_e_to_v}),
           std::move(subentities), std::move(subvertices0)};
 }
 //-----------------------------------------------------------------------------
diff --git a/cpp/dolfinx/mesh/Topology.h b/cpp/dolfinx/mesh/Topology.h
index 1b0517b22d1..7043bfa5e66 100644
--- a/cpp/dolfinx/mesh/Topology.h
+++ b/cpp/dolfinx/mesh/Topology.h
@@ -1,4 +1,4 @@
-// Copyright (C) 2006-2022 Anders Logg and Garth N. Wells
+// Copyright (C) 2006-2024 Anders Logg and Garth N. Wells
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -9,10 +9,12 @@
 #include <array>
 #include <cstdint>
 #include <dolfinx/common/MPI.h>
+#include <map>
 #include <memory>
 #include <optional>
 #include <span>
 #include <tuple>
+#include <utility>
 #include <vector>
 
 namespace dolfinx::common
@@ -38,33 +40,19 @@ enum class CellType;
 /// where dim is the topological dimension and i is the index of the
 /// entity within that topological dimension.
 ///
-/// @todo Rework memory management and associated API. Currently, there
-/// is no clear caching policy implemented and no way of discarding
-/// cached data.
+/// @todo Rework memory management and associated API. Currently, the
+/// caching policy is not clear.
 class Topology
 {
 public:
-  /// @brief Topology constructor.
-  /// @param[in] comm MPI communicator.
-  /// @param[in] cell_type Type of cell.
-  /// @param[in] vertex_map Index map describing the distribution of
-  /// mesh vertices.
-  /// @param[in] cell_map Index map describing the distribution of mesh
-  /// cells.
-  /// @param[in] cells Cell-to-vertex connectivity.
-  /// @param[in] original_index Original index for each cell in `cells`.
-  Topology(MPI_Comm comm, CellType cell_type,
-           std::shared_ptr<const common::IndexMap> vertex_map,
-           std::shared_ptr<const common::IndexMap> cell_map,
-           std::shared_ptr<graph::AdjacencyList<std::int32_t>> cells,
-           const std::optional<std::vector<std::int64_t>>& original_index
-           = std::nullopt);
-
-  /// @brief Create empty mesh topology with multiple cell types.
+  /// @brief Create a mesh topology.
   ///
-  /// @warning Experimental
+  /// A Topology represents the connectivity of a mesh. Mesh entities,
+  /// i.e. vertices, edges, faces and cells, are defined in terms of
+  /// their vertices. Connectivity represents the relationships between
+  /// entities, e.g. the cells that are connected to a given edge in the
+  /// mesh.
   ///
-  /// @param[in] comm MPI communicator.
   /// @param[in] cell_types Types of cells.
   /// @param[in] vertex_map Index map describing the distribution of
   /// mesh vertices.
@@ -75,7 +63,7 @@ class Topology
   /// @param[in] original_cell_index Original indices for each cell in
   /// `cells`.
   Topology(
-      MPI_Comm comm, std::vector<CellType> cell_types,
+      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,
@@ -98,13 +86,13 @@ class Topology
   /// Assignment
   Topology& operator=(Topology&& topology) = default;
 
-  /// @brief Return the topological dimension of the mesh.
+  /// @brief Topological dimension of the mesh.
   int dim() const noexcept;
 
   /// @brief Entity types in the topology for a given dimension.
   /// @param[in] dim Topological dimension.
   /// @return Entity types.
-  std::vector<CellType> entity_types(int dim) const;
+  const std::vector<CellType>& entity_types(int dim) const;
 
   /// @brief Cell type.
   ///
@@ -116,8 +104,6 @@ class Topology
   /// @brief Get the index maps that described the parallel distribution
   /// of the mesh entities of a given topological dimension.
   ///
-  /// @warning Experimental
-  ///
   /// @param[in] dim Topological dimension.
   /// @return Index maps, one for each cell type.
   std::vector<std::shared_ptr<const common::IndexMap>>
@@ -131,34 +117,35 @@ class Topology
   /// `nullptr` if index map has not been set.
   std::shared_ptr<const common::IndexMap> index_map(int dim) const;
 
-  /// @brief Get the connectivity from entities of topological
-  /// dimension d0 to dimension d1.
+  /// @brief Get the connectivity from entities of topological dimension
+  /// `d0` to dimension `d1`.
   ///
-  /// The entity type, and incident entity type are each described by a
-  /// pair (dim, index). The index within a topological dimension `dim`,
-  /// is that of the cell type given in `entity_types(dim)`.
+  /// The entity type and incident entity type are each described by a
+  /// pair `(dim, index)`. The index within a topological dimension
+  /// `dim`, is that of the cell type given in `entity_types(dim)`.
   ///
   /// @param[in] d0 Pair of (topological dimension of entities, index of
   /// "entity type" within topological dimension).
   /// @param[in] d1 Pair of (topological dimension of entities, index of
   /// incident "entity type" within topological dimension).
-  /// @return AdjacencyList of connectivity from entity type in d0 to
-  /// entity types in d1, or nullptr if not yet computed.
+  /// @return AdjacencyList of connectivity from entity type in `d0` to
+  /// entity types in `d1`, or `nullptr` if not yet computed.
   std::shared_ptr<const graph::AdjacencyList<std::int32_t>>
   connectivity(std::array<int, 2> d0, std::array<int, 2> d1) const;
 
-  /// @brief Return connectivity from entities of dimension d0 to
-  /// entities of dimension d1. Assumes only one entity type per dimension.
+  /// @brief Return connectivity from entities of dimension `d0` to
+  /// entities of dimension `d1`. Assumes only one entity type per
+  /// dimension.
   ///
-  /// @param[in] d0
-  /// @param[in] d1
-  /// @return The adjacency list that for each entity of dimension d0
-  /// gives the list of incident entities of dimension d1. Returns
+  /// @param[in] d0 Topological dimension.
+  /// @param[in] d1 Topological dimension.
+  /// @return The adjacency list that for each entity of dimension `d0`
+  /// gives the list of incident entities of dimension `d1`. Returns
   /// `nullptr` if connectivity has not been computed.
   std::shared_ptr<const graph::AdjacencyList<std::int32_t>>
   connectivity(int d0, int d1) const;
 
-  /// @brief Returns the permutation information
+  /// @brief Returns the permutation information.
   const std::vector<std::uint32_t>& get_cell_permutation_info() const;
 
   /// @brief Get the numbers that encode the number of permutations to
@@ -201,49 +188,11 @@ class Topology
   /// been computed.
   const std::vector<std::int32_t>& interprocess_facets() const;
 
-  /// @todo Merge with set_connectivity
-  ///
-  /// @brief Set the IndexMap for the `i`th `celltype` of dimension `dim`.
-  ///
-  /// @warning This is experimental and likely to change.
-  ///
-  /// @param[in] dim Topological dimension.
-  /// @param[in] i Index of cell type within dimension `dim`. The cell types
-  /// in the mesh for a given dimension are returned by ::entity_types.
-  /// @param[in] map Index map to set.
-  void set_index_map(int dim, int i,
-                     std::shared_ptr<const common::IndexMap> map);
-  /// @todo Merge with set_connectivity
-  ///
-  /// @brief Set the IndexMap for dimension dim
-  /// @warning This is experimental and likely to change
-  void set_index_map(int dim, std::shared_ptr<const common::IndexMap> map);
-
-  /// @brief Set connectivity for given pair of entity types, defined by
-  /// dimension and index, as listed in ::entity_types.
-  ///
-  /// General version for mixed topology. Connectivity from d0 to d1.
-  ///
-  /// @warning Experimental.
-  ///
-  /// @param[in] c Connectivity.
-  /// @param[in] d0 Pair of (topological dimension of entities, index of
-  /// "entity type" within topological dimension).
-  /// @param[in] d1 Pair of (topological dimension of incident entities,
-  /// index of incident "entity type" within topological dimension).
-  void set_connectivity(std::shared_ptr<graph::AdjacencyList<std::int32_t>> c,
-                        std::array<int, 2> d0, std::array<int, 2> d1);
-
-  /// @todo Merge with set_index_map
-  /// @brief Set connectivity for given pair of topological dimensions.
-  void set_connectivity(std::shared_ptr<graph::AdjacencyList<std::int32_t>> c,
-                        int d0, int d1);
-
   /// @brief Create entities of given topological dimension.
   /// @param[in] dim Topological dimension of entities to compute.
-  /// @return Number of newly created entities, returns -1 if entities
-  /// already existed
-  std::int32_t create_entities(int dim);
+  /// @return True if entities are created, false if entities already
+  /// existed.
+  bool create_entities(int dim);
 
   /// @brief Create connectivity between given pair of dimensions, `d0
   /// -> d1`.
@@ -257,36 +206,30 @@ class Topology
   /// Original cell index for each cell type
   std::vector<std::vector<std::int64_t>> original_cell_index;
 
-  /// Mesh MPI communicator
-  /// @return The communicator on which the topology is distributed
+  /// @brief Mesh MPI communicator.
+  /// @return Communicator on which the topology is distributed.
   MPI_Comm comm() const;
 
 private:
-  // MPI communicator
-  dolfinx::MPI::Comm _comm;
-
-  // Cell types for entities in Topology, as follows:
-  //
-  // [CellType::point, edge_types..., facet_types..., cell_types...]
-  //
-  // Only one type is expected for vertices, (and usually edges), but
-  // facets and cells can be a list of multiple types, e.g.
-  // [quadrilateral, triangle] for facets. Offsets are position in the
-  // list for each entity dimension, in AdjacencyList style.
-  std::vector<CellType> _entity_types;
-  std::vector<std::int8_t> _entity_type_offsets;
+  // 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;
-
-  // Connectivity between entity dimensions and cell types, arranged as
-  // a 2D array. The indexing follows the order in _entity_types, i.e.
-  // increasing in topological dimension. There may be multiple types in
-  // each dimension, e.g. triangle and quadrilateral facets.
-  // Connectivity between different entity types of same dimension will
-  // always be nullptr.
-  std::vector<std::vector<std::shared_ptr<graph::AdjacencyList<std::int32_t>>>>
+  // 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))
@@ -299,7 +242,8 @@ class Topology
   std::vector<std::uint32_t> _cell_permutations;
 
   // List of facets that are on the inter-process boundary for each
-  // facet type
+  // facet type. _interprocess_facets[i] is the inter-process facets of
+  // facet type i.
   std::vector<std::vector<std::int32_t>> _interprocess_facets;
 };
 
@@ -308,24 +252,28 @@ class Topology
 /// This function creates a Topology from cells that have been already
 /// distributed to the processes that own or ghost the cell.
 ///
-/// @param[in] comm Communicator across which the topology is
+/// @param[in] comm Communicator across which the topology will be
 /// distributed.
 /// @param[in] cell_types List of cell types in the topology.
 /// @param[in] cells Cell topology (list of vertices for each cell) for
 /// each cell type using global indices for the vertices. The cell type
-/// for `cells[i]` is `cell_types[i]`. Each `cells[i]` contains cells
-/// that have been distributed to this rank, e.g. via a graph
-/// partitioner. It must also contain all ghost cells via facet, i.e.
-/// cells that are on a neighboring process and which share a facet with
-/// a local cell. Ghost cells are the last `n` entries in `cells[i]`, where
-/// `n` is given by the length of `ghost_owners[i]`.
-/// @param[in] original_cell_index Input cell index for each cell type.
-/// @param[in] ghost_owners Owning rank for ghost cells (at end of each list of
-/// cells).
+/// for `cells[i]` is `cell_types[i]`, using row-major storage and where
+/// the row `cells[i][j]` is the vertices for cell `j` of cell type `i`.
+/// Each `cells[i]` contains cells that have been distributed to this
+/// rank, e.g. via a graph partitioner. It must also contain all ghost
+/// cells via facet, i.e. cells that are on a neighboring process and
+/// which share a facet with a local cell. Ghost cells are the last `n`
+/// entries in `cells[i]`, where `n` is given by the length of
+/// `ghost_owners[i]`.
+/// @param[in] original_cell_index Input cell index for each cell type,
+/// e.g. the cell index in an input file. This index remains associated
+/// with the cell after any re-ordering and parallel (re)distribution.
+/// @param[in] ghost_owners Owning rank for ghost cells (ghost cells are
+/// at end of each list of cells).
 /// @param[in] boundary_vertices Vertices on the 'exterior' (boundary)
 /// of the local topology. These vertices might appear on other
 /// processes.
-/// @return A distributed mesh topology
+/// @return A distributed mesh topology.
 Topology
 create_topology(MPI_Comm comm, const std::vector<CellType>& cell_types,
                 std::vector<std::span<const std::int64_t>> cells,
@@ -335,8 +283,10 @@ create_topology(MPI_Comm comm, const std::vector<CellType>& cell_types,
 
 /// @brief Create a mesh topology for a single cell type.
 ///
+/// This function provides a simplified interface to ::create_topology
+/// for the case that a mesh has one cell type only,
 ///
-/// @param[in] comm Communicator across which the topology is
+/// @param[in] comm Communicator across which the topology will be
 /// distributed.
 /// @param[in] cells Cell topology (list of vertices for each cell)
 /// using global indices for the vertices. It contains cells that have
@@ -353,7 +303,7 @@ create_topology(MPI_Comm comm, const std::vector<CellType>& cell_types,
 /// @param[in] boundary_vertices Vertices on the 'exterior' (boundary)
 /// of the local topology. These vertices might appear on other
 /// processes.
-/// @return A distributed mesh topology
+/// @return A distributed mesh topology.
 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,
@@ -379,12 +329,11 @@ create_subtopology(const Topology& topology, int dim,
 ///
 /// @warning This function may be removed in the future.
 ///
-/// @param[in] topology The mesh topology
-/// @param[in] dim Topological dimension of the entities
-/// @param[in] entities The mesh entities defined by their vertices
-/// @return The index of the ith entity in `entities`
-/// @note If an entity cannot be found on this rank, -1 is returned as
-/// the index.
+/// @param[in] topology Mesh topology.
+/// @param[in] dim Topological dimension of the entities.
+/// @param[in] entities Mesh entities defined by their vertices.
+/// @return Index of the ith entity in `entities`, If an entity
+/// cannot be found on this rank, -1 is returned as the index.
 std::vector<std::int32_t>
 entities_to_index(const Topology& topology, int dim,
                   std::span<const std::int32_t> entities);
diff --git a/cpp/dolfinx/mesh/topologycomputation.cpp b/cpp/dolfinx/mesh/topologycomputation.cpp
index 8b09085bae6..9a39836977a 100644
--- a/cpp/dolfinx/mesh/topologycomputation.cpp
+++ b/cpp/dolfinx/mesh/topologycomputation.cpp
@@ -1,4 +1,4 @@
-// Copyright (C) 2006-2020 Anders Logg, Garth N. Wells and Chris Richardson
+// Copyright (C) 2006-2024 Anders Logg, Garth N. Wells and Chris Richardson
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -17,6 +17,7 @@
 #include <dolfinx/common/sort.h>
 #include <dolfinx/graph/AdjacencyList.h>
 #include <memory>
+#include <mpi.h>
 #include <numeric>
 #include <random>
 #include <string>
@@ -753,6 +754,7 @@ compute_from_map(const graph::AdjacencyList<std::int32_t>& c_d0_0,
       connections.push_back(it->second);
     }
   }
+
   connections.shrink_to_fit();
   return graph::AdjacencyList(std::move(connections), std::move(offsets));
 }
@@ -763,27 +765,30 @@ compute_from_map(const graph::AdjacencyList<std::int32_t>& c_d0_0,
 std::tuple<std::vector<std::shared_ptr<graph::AdjacencyList<std::int32_t>>>,
            std::shared_ptr<graph::AdjacencyList<std::int32_t>>,
            std::shared_ptr<common::IndexMap>, std::vector<std::int32_t>>
-mesh::compute_entities(MPI_Comm comm, const Topology& topology, int dim,
-                       int index)
+mesh::compute_entities(const Topology& topology, int dim, CellType entity_type)
 {
   spdlog::info("Computing mesh entities of dimension {}", dim);
-  const int tdim = topology.dim();
 
   // Vertices must always exist
   if (dim == 0)
-    return {std::vector<std::shared_ptr<graph::AdjacencyList<std::int32_t>>>(),
-            nullptr, nullptr, std::vector<std::int32_t>()};
-
-  if (topology.connectivity({dim, index}, {0, 0}))
   {
     return {std::vector<std::shared_ptr<graph::AdjacencyList<std::int32_t>>>(),
             nullptr, nullptr, std::vector<std::int32_t>()};
   }
 
-  auto vertex_map = topology.index_map(0);
-  assert(vertex_map);
+  {
+    auto idx = std::ranges::find(topology.entity_types(dim), entity_type);
+    assert(idx != topology.entity_types(dim).end());
+    int index = std::distance(topology.entity_types(dim).begin(), idx);
+    if (topology.connectivity({dim, index}, {0, 0}))
+    {
+      return {
+          std::vector<std::shared_ptr<graph::AdjacencyList<std::int32_t>>>(),
+          nullptr, nullptr, std::vector<std::int32_t>()};
+    }
+  }
 
-  CellType entity_type = topology.entity_types(dim)[index];
+  const int tdim = topology.dim();
 
   // Lists of all cells by cell type
   std::vector<CellType> cell_types = topology.entity_types(tdim);
@@ -803,8 +808,12 @@ mesh::compute_entities(MPI_Comm comm, const Topology& topology, int dim,
     cell_lists[i] = {cell_types[i], cells, cell_map};
   }
 
+  auto vertex_map = topology.index_map(0);
+  assert(vertex_map);
+
+  // c->e, e->v
   auto [d0, d1, im, interprocess_entities] = compute_entities_by_key_matching(
-      comm, cell_lists, *vertex_map, entity_type, dim);
+      topology.comm(), cell_lists, *vertex_map, entity_type, dim);
 
   return {d0,
           std::make_shared<graph::AdjacencyList<std::int32_t>>(std::move(d1)),
diff --git a/cpp/dolfinx/mesh/topologycomputation.h b/cpp/dolfinx/mesh/topologycomputation.h
index 61f64716010..16a1365da3c 100644
--- a/cpp/dolfinx/mesh/topologycomputation.h
+++ b/cpp/dolfinx/mesh/topologycomputation.h
@@ -1,4 +1,4 @@
-// Copyright (C) 2006-2020 Anders Logg and Garth N. Wells
+// Copyright (C) 2006-2024 Anders Logg and Garth N. Wells
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -6,10 +6,10 @@
 
 #pragma once
 
+#include "cell_types.h"
 #include <array>
 #include <cstdint>
 #include <memory>
-#include <mpi.h>
 #include <tuple>
 #include <vector>
 
@@ -29,33 +29,35 @@ namespace dolfinx::mesh
 class Topology;
 
 /// @brief Compute mesh entities of given topological dimension by
-/// computing entity-to-vertex connectivity `(dim, 0)`, and
-/// cell-to-entity connectivity `(tdim, dim)`.
+/// computing cell-to-entity `(tdim, i) -> `(dim, entity_type)` and
+/// entity-to-vertex connectivity `(dim, entity_type) -> `(0, 0)`
+/// connectivity.
 ///
 /// Computed entities are oriented such that their local (to the
 /// process) orientation agrees with their global orientation
 ///
-/// @param[in] comm MPI Communicator
-/// @param[in] topology Mesh topology
-/// @param[in] dim The dimension of the entities to create
-/// @param[in] index Index of entity in dimension `dim` as listed in
-/// `Topology::entity_types(dim)`.
-/// @return Tuple of (cell-entity connectivity, entity-vertex
-/// connectivity, index map, list of interprocess entities).
-/// Interprocess entities lie on the "true" boundary between owned cells
-/// of each process. If the entities already exists, then {nullptr,
-/// nullptr, nullptr, std::vector()} is returned.
+/// @param[in] topology Mesh topology.
+/// @param[in] dim Dimension of the entities to create.
+/// @param[in] entity_type Entity type in dimension `dim` to create.
+/// Entity type must be in the list returned by Topology::entity_types.
+/// @return Tuple of (cell->entity connectivity, entity->vertex
+/// connectivity, index map for created entities, list of interprocess
+/// entities). Interprocess entities lie on the "true" boundary between
+/// owned cells of each process. If entities of type `entity_type`
+/// already exists, then {nullptr, nullptr, nullptr, std::vector()} is
+/// returned.
 std::tuple<std::vector<std::shared_ptr<graph::AdjacencyList<std::int32_t>>>,
            std::shared_ptr<graph::AdjacencyList<std::int32_t>>,
            std::shared_ptr<common::IndexMap>, std::vector<std::int32_t>>
-compute_entities(MPI_Comm comm, const Topology& topology, int dim, int index);
+compute_entities(const Topology& topology, int dim, CellType entity_type);
 
 /// @brief Compute connectivity (d0 -> d1) for given pair of entity
 /// types, given by topological dimension and index, as found in
 /// `Topology::entity_types()`
 /// @param[in] topology The topology
-/// @param[in] d0 The dimension and index of the entities
-/// @param[in] d1 The dimension and index of the incident entities
+/// @param[in] d0 Dimension and index of the entities, `(dim0, i)`.
+/// @param[in] d1 Dimension and index of the incident entities, `(dim1,
+/// j)`.
 /// @returns The connectivities [(d0 -> d1), (d1 -> d0)] if they are
 /// computed. If (d0, d1) already exists then a nullptr is returned. If
 /// (d0, d1) is computed and the computation of (d1, d0) was required as
diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index 48114afb546..6dace4d3981 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -109,22 +109,23 @@ def comm(self):
         return self._cpp_object.comm
 
     def create_connectivity(self, d0: int, d1: int):
-        """Create connectivity between given pair of dimensions, ``d0`` and  ``d1``.
+        """Build entity connectivity ``d0 -> d1``.
 
         Args:
-            d0: Dimension of entities one is mapping from
-            d1: Dimension of entities one is mapping to
+            d0: Dimension of entities connectivity is from.
+            d1: Dimension of entities connectivity is to.
         """
         self._cpp_object.create_connectivity(d0, d1)
 
-    def create_entities(self, dim: int) -> int:
+    def create_entities(self, dim: int) -> bool:
         """Create entities of given topological dimension.
 
         Args:
-            dim: Topological dimension
+            dim: Topological dimension of entities to create.
 
         Returns:
-            Number of newly created entities, returns -1 if entities already existed
+            ``True` is entities are created, ``False`` is if entities
+            already existed.
         """
         return self._cpp_object.create_entities(dim)
 
@@ -168,10 +169,10 @@ def get_facet_permutations(self) -> npt.NDArray[np.uint8]:
         return self._cpp_object.get_facet_permutations()
 
     def index_map(self, dim: int) -> _cpp.common.IndexMap:
-        """Get the IndexMap that described the parallel distribution of the mesh entities.
+        """Get the IndexMap that describes the parallel distribution of the mesh entities.
 
         Args:
-            dim: Topological dimension
+            dim: Topological dimension.
 
         Returns:
             Index map for the entities of dimension ``dim``.
@@ -193,28 +194,9 @@ def original_cell_index(self) -> npt.NDArray[np.int64]:
         """Get the original cell index"""
         return self._cpp_object.original_cell_index
 
-    def set_connectivity(self, graph: _cpp.graph.AdjacencyList_int32, d0: int, d1: int):
-        """Set connectivity for given pair of topological dimensions.
-
-        Args:
-            graph: Connectivity graph
-            d0: Topological dimension mapping from
-            d1: Topological dimension mapping to
-        """
-        self._cpp_object.set_connectivity(graph, d0, d1)
-
-    def set_index_map(self, dim: int, index_map: _cpp.common.IndexMap):
-        """Set the IndexMap for dimension ``dim``.
-
-        Args:
-            dim: Topological dimension of entity.
-            index_map: Index map to store.
-        """
-        return self._cpp_object.set_index_map(dim, index_map)
-
     @property
     def cell_type(self) -> CellType:
-        """Get the cell type of the topology"""
+        """Get the cell type of the topology."""
         return self._cpp_object.cell_type
 
 
diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp
index 487f75209b9..f0198352adb 100644
--- a/python/dolfinx/wrappers/mesh.cpp
+++ b/python/dolfinx/wrappers/mesh.cpp
@@ -559,13 +559,10 @@ void mesh(nb::module_& m)
   // dolfinx::mesh::TopologyComputation
   m.def(
       "compute_entities",
-      [](MPICommWrapper comm, const dolfinx::mesh::Topology& topology, int dim,
-         int index)
-      {
-        return dolfinx::mesh::compute_entities(comm.get(), topology, dim,
-                                               index);
-      },
-      nb::arg("comm"), nb::arg("topology"), nb::arg("dim"), nb::arg("index"));
+      [](const dolfinx::mesh::Topology& topology, int dim,
+         dolfinx::mesh::CellType entity_type)
+      { return dolfinx::mesh::compute_entities(topology, dim, entity_type); },
+      nb::arg("topology"), nb::arg("dim"), nb::arg("entity_type"));
   m.def("compute_connectivity", &dolfinx::mesh::compute_connectivity,
         nb::arg("topology"), nb::arg("d0"), nb::arg("d1"));
 
@@ -574,8 +571,7 @@ void mesh(nb::module_& m)
                                       "Topology object")
       .def(
           "__init__",
-          [](dolfinx::mesh::Topology* t, MPICommWrapper comm,
-             dolfinx::mesh::CellType cell_type,
+          [](dolfinx::mesh::Topology* t, dolfinx::mesh::CellType cell_type,
              std::shared_ptr<const dolfinx::common::IndexMap> vertex_map,
              std::shared_ptr<const dolfinx::common::IndexMap> cell_map,
              std::shared_ptr<dolfinx::graph::AdjacencyList<std::int32_t>> cells,
@@ -583,28 +579,18 @@ void mesh(nb::module_& m)
                  nb::ndarray<const std::int64_t, nb::ndim<1>, nb::c_contig>>
                  original_index)
           {
-            std::optional<std::vector<std::int64_t>> idx
-                = original_index
-                      ? std::vector<std::int64_t>(original_index->data(),
-                                                  original_index->data()
-                                                      + original_index->size())
-                      : std::optional<std::vector<std::int64_t>>(std::nullopt);
-            new (t) dolfinx::mesh::Topology(comm.get(), cell_type, vertex_map,
-                                            cell_map, cells, idx);
+            using U = std::vector<std::vector<std::int64_t>>;
+            using V = std::optional<U>;
+            V idx = original_index
+                        ? U(1, std::vector(original_index->data(),
+                                           original_index->data()
+                                               + original_index->size()))
+                        : V(std::nullopt);
+            new (t) dolfinx::mesh::Topology({cell_type}, vertex_map, {cell_map},
+                                            {cells}, idx);
           },
-          nb::arg("comm"), nb::arg("cell_type"), nb::arg("vertex_map"),
-          nb::arg("cell_map"), nb::arg("cells"),
-          nb::arg("original_index").none())
-      .def("set_connectivity",
-           nb::overload_cast<
-               std::shared_ptr<dolfinx::graph::AdjacencyList<std::int32_t>>,
-               int, int>(&dolfinx::mesh::Topology::set_connectivity),
-           nb::arg("c"), nb::arg("d0"), nb::arg("d1"))
-      .def("set_index_map",
-           nb::overload_cast<int,
-                             std::shared_ptr<const dolfinx::common::IndexMap>>(
-               &dolfinx::mesh::Topology::set_index_map),
-           nb::arg("dim"), nb::arg("map"))
+          nb::arg("cell_type"), nb::arg("vertex_map"), nb::arg("cell_map"),
+          nb::arg("cells"), nb::arg("original_index").none())
       .def("create_entities", &dolfinx::mesh::Topology::create_entities,
            nb::arg("dim"))
       .def("create_entity_permutations",
diff --git a/python/test/unit/mesh/test_mixed_topology.py b/python/test/unit/mesh/test_mixed_topology.py
index 842d28e467d..3af3e27744e 100644
--- a/python/test/unit/mesh/test_mixed_topology.py
+++ b/python/test/unit/mesh/test_mixed_topology.py
@@ -28,6 +28,7 @@ def test_mixed_topology_mesh():
 
     maps = topology.index_maps(topology.dim)
     assert len(maps) == 2
+
     # Two triangles and one quad
     assert maps[0].size_local == 2
     assert maps[1].size_local == 1
@@ -39,12 +40,19 @@ def test_mixed_topology_mesh():
 
     entity_types = topology.entity_types
     assert len(entity_types[0]) == 1
+
+    topology.create_entities(1)
+    entity_types = topology.entity_types
     assert len(entity_types[1]) == 1
-    assert len(entity_types[2]) == 2
     assert CellType.interval in entity_types[1]
+
+    entity_types = topology.entity_types
+    assert len(entity_types[2]) == 2
+
     # Two triangle cells
     assert entity_types[2][0] == CellType.triangle
     assert topology.connectivity((2, 0), (0, 0)).num_nodes == 2
+
     # One quadrlilateral cell
     assert entity_types[2][1] == CellType.quadrilateral
     assert topology.connectivity((2, 1), (0, 0)).num_nodes == 1
@@ -81,8 +89,15 @@ def test_mixed_topology_mesh_3d():
 
     entity_types = topology.entity_types
     assert len(entity_types[0]) == 1
+
+    topology.create_entities(1)
+    entity_types = topology.entity_types
     assert len(entity_types[1]) == 1
+
+    topology.create_entities(2)
+    entity_types = topology.entity_types
     assert len(entity_types[2]) == 2
+
     assert len(entity_types[3]) == 3
 
     # Create triangle and quadrilateral facets