diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
index f775544dd9f..0491a8714d5 100644
--- a/cpp/demo/hyperelasticity/main.cpp
+++ b/cpp/demo/hyperelasticity/main.cpp
@@ -222,7 +222,7 @@ int main(int argc, char* argv[])
     // Construct appropriate Basix element for stress
     constexpr auto family = basix::element::family::P;
     const auto cell_type
-        = mesh::cell_type_to_basix_type(mesh->topology().cell_type());
+        = mesh::cell_type_to_basix_type(mesh->topology().cell_types()[0]);
     constexpr int k = 0;
     constexpr bool discontinuous = true;
 
diff --git a/cpp/demo/interpolation_different_meshes/main.cpp b/cpp/demo/interpolation_different_meshes/main.cpp
index 49e3e46375a..218674555c5 100644
--- a/cpp/demo/interpolation_different_meshes/main.cpp
+++ b/cpp/demo/interpolation_different_meshes/main.cpp
@@ -33,13 +33,13 @@ int main(int argc, char* argv[])
                          mesh::CellType::hexahedron));
 
     basix::FiniteElement element_tet = basix::element::create_lagrange(
-        mesh::cell_type_to_basix_type(mesh_tet->topology().cell_type()), 1,
+        mesh::cell_type_to_basix_type(mesh_tet->topology().cell_types()[0]), 1,
         basix::element::lagrange_variant::equispaced, false);
     auto V_tet = std::make_shared<fem::FunctionSpace<double>>(
         fem::create_functionspace(mesh_tet, element_tet, 3));
 
     basix::FiniteElement element_hex = basix::element::create_lagrange(
-        mesh::cell_type_to_basix_type(mesh_hex->topology().cell_type()), 2,
+        mesh::cell_type_to_basix_type(mesh_hex->topology().cell_types()[0]), 2,
         basix::element::lagrange_variant::equispaced, false);
     auto V_hex = std::make_shared<fem::FunctionSpace<double>>(
         fem::create_functionspace(mesh_hex, element_hex, 3));
diff --git a/cpp/demo/mixed_topology/CMakeLists.txt b/cpp/demo/mixed_topology/CMakeLists.txt
new file mode 100644
index 00000000000..3b4687b3d66
--- /dev/null
+++ b/cpp/demo/mixed_topology/CMakeLists.txt
@@ -0,0 +1,24 @@
+cmake_minimum_required(VERSION 3.16)
+project(mixed_topology_testing)
+
+# Find DOLFINx config file
+find_package(DOLFINX REQUIRED)
+
+# Set C++20 standard
+set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+
+include(CheckSymbolExists)
+set(CMAKE_REQUIRED_INCLUDES ${PETSC_INCLUDE_DIRS})
+check_symbol_exists(PETSC_USE_COMPLEX petscsystypes.h PETSC_SCALAR_COMPLEX)
+
+# Add target to compile UFL files
+if(PETSC_SCALAR_COMPLEX EQUAL 1)
+  set(SCALAR_TYPE "--scalar_type=double _Complex")
+endif()
+
+set(CMAKE_INCLUDE_CURRENT_DIR ON)
+
+add_executable(${PROJECT_NAME} main.cpp)
+target_link_libraries(${PROJECT_NAME} dolfinx)
diff --git a/cpp/demo/mixed_topology/main.cpp b/cpp/demo/mixed_topology/main.cpp
new file mode 100644
index 00000000000..2b500d56ddc
--- /dev/null
+++ b/cpp/demo/mixed_topology/main.cpp
@@ -0,0 +1,149 @@
+#include <dolfinx/common/IndexMap.h>
+#include <dolfinx/common/MPI.h>
+#include <dolfinx/graph/AdjacencyList.h>
+#include <dolfinx/mesh/Geometry.h>
+#include <dolfinx/mesh/Mesh.h>
+#include <dolfinx/mesh/Topology.h>
+#include <dolfinx/mesh/cell_types.h>
+#include <iostream>
+#include <vector>
+
+// Note: this demo is not currently intended to provide a fully functional
+// example of using a mixed-topology mesh, but shows only the
+// basic constrution. Experimental.
+
+int main(int argc, char* argv[])
+{
+  dolfinx::init_logging(argc, argv);
+  MPI_Init(&argc, &argv);
+
+  const int nx_s = 2;
+  const int nx_t = 2;
+  const int ny = 4;
+
+  const int num_s = nx_s * ny;
+  const int num_t = 2 * nx_t * ny;
+
+  std::vector<double> x;
+  for (int i = 0; i < nx_s + nx_t + 1; ++i)
+  {
+    for (int j = 0; j < ny + 1; ++j)
+    {
+      x.push_back(i);
+      x.push_back(j);
+      x.push_back(0);
+    }
+  }
+
+  std::vector<std::int64_t> cells;
+  std::vector<std::int32_t> offsets{0};
+  for (int i = 0; i < nx_s + nx_t; ++i)
+  {
+    for (int j = 0; j < ny; ++j)
+    {
+      const int v_0 = j + i * (ny + 1);
+      const int v_1 = v_0 + 1;
+      const int v_2 = v_0 + ny + 1;
+      const int v_3 = v_0 + ny + 2;
+
+      if (i < nx_s)
+      {
+        cells.push_back(v_0);
+        cells.push_back(v_1);
+        cells.push_back(v_3);
+        cells.push_back(v_2);
+
+        offsets.push_back(cells.size());
+      }
+      else
+      {
+        cells.push_back(v_0);
+        cells.push_back(v_1);
+        cells.push_back(v_2);
+
+        offsets.push_back(cells.size());
+
+        cells.push_back(v_1);
+        cells.push_back(v_2);
+        cells.push_back(v_3);
+
+        offsets.push_back(cells.size());
+      }
+    }
+  }
+
+  dolfinx::graph::AdjacencyList<std::int64_t> cells_list(cells, offsets);
+
+  std::vector<std::int64_t> original_global_index(num_s + num_t);
+  std::iota(original_global_index.begin(), original_global_index.end(), 0);
+
+  std::vector<int> ghost_owners;
+
+  std::vector<std::int32_t> cell_group_offsets{0, num_s, num_s + num_t,
+                                               num_s + num_t, num_s + num_t};
+
+  std::vector<std::int64_t> boundary_vertices;
+  for (int j = 0; j < ny + 1; ++j)
+  {
+    boundary_vertices.push_back(j);
+    boundary_vertices.push_back(j + (nx_s + nx_t + 1) * ny);
+  }
+  for (int i = 0; i < nx_s + nx_t + 1; ++i)
+  {
+    boundary_vertices.push_back((ny + 1) * i);
+    boundary_vertices.push_back(ny + (ny + 1) * i);
+  }
+
+  std::sort(boundary_vertices.begin(), boundary_vertices.end());
+  boundary_vertices.erase(
+      std::unique(boundary_vertices.begin(), boundary_vertices.end()),
+      boundary_vertices.end());
+
+  std::vector<dolfinx::mesh::CellType> cell_types{
+      dolfinx::mesh::CellType::quadrilateral,
+      dolfinx::mesh::CellType::triangle};
+  std::vector<dolfinx::fem::CoordinateElement> elements;
+  for (auto ct : cell_types)
+    elements.push_back(dolfinx::fem::CoordinateElement(ct, 1));
+
+  {
+    auto topo = dolfinx::mesh::create_topology(
+        MPI_COMM_WORLD, cells_list, original_global_index, ghost_owners,
+        cell_types, cell_group_offsets, boundary_vertices);
+
+    auto topo_cells = topo.connectivity(2, 0);
+
+    for (int i = 0; i < topo_cells->num_nodes(); ++i)
+    {
+      std::cout << i << " [";
+      for (auto q : topo_cells->links(i))
+        std::cout << q << " ";
+      std::cout << "]\n";
+    }
+
+    topo.create_connectivity(1, 0);
+
+    auto topo_facets = topo.connectivity(1, 0);
+    for (int i = 0; i < topo_facets->num_nodes(); ++i)
+    {
+      std::cout << i << " [";
+      for (auto q : topo_facets->links(i))
+        std::cout << q << " ";
+      std::cout << "]\n";
+    }
+
+    auto geom = dolfinx::mesh::create_geometry(MPI_COMM_WORLD, topo, elements,
+                                               cells_list, x, 2);
+
+    dolfinx::mesh::Mesh mesh(MPI_COMM_WORLD, topo, geom);
+    std::cout << "num cells = " << mesh.topology().index_map(2)->size_local()
+              << "\n";
+    for (auto q : mesh.topology().entity_group_offsets(2))
+      std::cout << q << " ";
+    std::cout << "\n";
+  }
+
+  MPI_Finalize();
+
+  return 0;
+}
diff --git a/cpp/dolfinx/fem/DirichletBC.cpp b/cpp/dolfinx/fem/DirichletBC.cpp
index b5ebfd952dd..0bfed840eca 100644
--- a/cpp/dolfinx/fem/DirichletBC.cpp
+++ b/cpp/dolfinx/fem/DirichletBC.cpp
@@ -182,10 +182,15 @@ fem::locate_dofs_topological(mesh::Topology& topology, const DofMap& dofmap,
                              int dim, std::span<const std::int32_t> entities,
                              bool remote)
 {
+  auto cell_types = topology.cell_types();
+  if (cell_types.size() > 1)
+  {
+    throw std::runtime_error("Multiple cell types in DirichletBC");
+  }
+
   // Prepare an element - local dof layout for dofs on entities of the
   // entity_dim
-  const int num_cell_entities
-      = mesh::cell_num_entities(topology.cell_type(), dim);
+  const int num_cell_entities = mesh::cell_num_entities(cell_types.back(), dim);
   std::vector<std::vector<int>> entity_dofs;
   for (int i = 0; i < num_cell_entities; ++i)
   {
@@ -298,9 +303,13 @@ std::array<std::vector<std::int32_t>, 2> fem::locate_dofs_topological(
   // Check that dof layouts are the same
   assert(dofmap0.element_dof_layout() == dofmap1.element_dof_layout());
 
+  auto cell_types = topology.cell_types();
+  if (cell_types.size() > 1)
+  {
+    throw std::runtime_error("Multiple cell types in DirichletBC");
+  }
   // Build vector of local dofs for each cell entity
-  const int num_cell_entities
-      = mesh::cell_num_entities(topology.cell_type(), dim);
+  const int num_cell_entities = mesh::cell_num_entities(cell_types.back(), dim);
   std::vector<std::vector<int>> entity_dofs;
   for (int i = 0; i < num_cell_entities; ++i)
   {
diff --git a/cpp/dolfinx/fem/DofMap.cpp b/cpp/dolfinx/fem/DofMap.cpp
index bd86f2887aa..12825a6b2f4 100644
--- a/cpp/dolfinx/fem/DofMap.cpp
+++ b/cpp/dolfinx/fem/DofMap.cpp
@@ -233,8 +233,8 @@ std::pair<DofMap, std::vector<std::int32_t>> DofMap::collapse(
       // Create new element dof layout and reset parent
       ElementDofLayout collapsed_dof_layout = layout.copy();
 
-      auto [_index_map, bs, dofmap]
-          = build_dofmap_data(comm, topology, collapsed_dof_layout, reorder_fn);
+      auto [_index_map, bs, dofmap] = build_dofmap_data(
+          comm, topology, {collapsed_dof_layout}, reorder_fn);
       auto index_map
           = std::make_shared<common::IndexMap>(std::move(_index_map));
       return DofMap(layout, index_map, bs, std::move(dofmap), bs);
diff --git a/cpp/dolfinx/fem/Expression.h b/cpp/dolfinx/fem/Expression.h
index 161a48357d6..6f229a37c4a 100644
--- a/cpp/dolfinx/fem/Expression.h
+++ b/cpp/dolfinx/fem/Expression.h
@@ -141,7 +141,12 @@ class Expression
     // Prepare cell geometry
     const graph::AdjacencyList<std::int32_t>& x_dofmap
         = _mesh->geometry().dofmap();
-    const std::size_t num_dofs_g = _mesh->geometry().cmap().dim();
+
+    // Get geometry data
+    auto cmaps = _mesh->geometry().cmaps();
+    assert(cmaps.size() == 1);
+
+    const std::size_t num_dofs_g = cmaps.back().dim();
     std::span<const U> x_g = _mesh->geometry().x();
 
     // Create data structures used in evaluation
diff --git a/cpp/dolfinx/fem/FiniteElement.h b/cpp/dolfinx/fem/FiniteElement.h
index 07a089db44a..3d66ce1350f 100644
--- a/cpp/dolfinx/fem/FiniteElement.h
+++ b/cpp/dolfinx/fem/FiniteElement.h
@@ -168,7 +168,7 @@ class FiniteElement
   /// @return True if the map is the identity
   bool map_ident() const noexcept;
 
-  /// @brief Points on the reference cell at which an expression need to
+  /// @brief Points on the reference cell at which an expression needs to
   /// be evaluated in order to interpolate the expression in the finite
   /// element space.
   ///
diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
index 2562062135c..c06fb1adba4 100644
--- a/cpp/dolfinx/fem/Function.h
+++ b/cpp/dolfinx/fem/Function.h
@@ -393,15 +393,20 @@ class Function
     const std::size_t tdim = mesh->topology().dim();
     auto map = mesh->topology().index_map(tdim);
 
+    // Get coordinate map
+    if (mesh->geometry().cmaps().size() > 1)
+    {
+      throw std::runtime_error(
+          "Function with multiple geometry maps not implemented.");
+    }
+    const CoordinateElement& cmap = mesh->geometry().cmaps()[0];
+
     // Get geometry data
     const graph::AdjacencyList<std::int32_t>& x_dofmap
         = mesh->geometry().dofmap();
-    const std::size_t num_dofs_g = mesh->geometry().cmap().dim();
+    const std::size_t num_dofs_g = cmap.dim();
     std::span<const U> x_g = mesh->geometry().x();
 
-    // Get coordinate map
-    const CoordinateElement& cmap = mesh->geometry().cmap();
-
     // Get element
     std::shared_ptr<const FiniteElement> element = _function_space->element();
     assert(element);
diff --git a/cpp/dolfinx/fem/FunctionSpace.h b/cpp/dolfinx/fem/FunctionSpace.h
index 3790bcb378e..96d224b34f9 100644
--- a/cpp/dolfinx/fem/FunctionSpace.h
+++ b/cpp/dolfinx/fem/FunctionSpace.h
@@ -216,12 +216,16 @@ class FunctionSpace
     const auto [X, Xshape] = _element->interpolation_points();
 
     // Get coordinate map
-    const CoordinateElement& cmap = _mesh->geometry().cmap();
+    if (_mesh->geometry().cmaps().size() > 1)
+    {
+      throw std::runtime_error("Mixed topology not supported");
+    }
+    const CoordinateElement& cmap = _mesh->geometry().cmaps()[0];
 
     // Prepare cell geometry
     const graph::AdjacencyList<std::int32_t>& x_dofmap
         = _mesh->geometry().dofmap();
-    const std::size_t num_dofs_g = _mesh->geometry().cmap().dim();
+    const std::size_t num_dofs_g = cmap.dim();
     std::span<const T> x_g = _mesh->geometry().x();
 
     // Array to hold coordinates to return
diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h
index 5b19f98b15b..144459aaa7b 100644
--- a/cpp/dolfinx/fem/assemble_matrix_impl.h
+++ b/cpp/dolfinx/fem/assemble_matrix_impl.h
@@ -48,7 +48,7 @@ void assemble_cells(
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Iterate over active cells
@@ -144,7 +144,7 @@ void assemble_exterior_facets(
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Data structures used in assembly
@@ -240,7 +240,7 @@ void assemble_interior_facets(
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Data structures used in assembly
@@ -445,7 +445,11 @@ void assemble_matrix(
     else
       get_perm = [](std::size_t) { return 0; };
 
-    int num_cell_facets = mesh::cell_num_entities(mesh->topology().cell_type(),
+    auto cell_types = mesh->topology().cell_types();
+    if (cell_types.size() > 1)
+      throw std::runtime_error("Multiple cell types in the assembler.");
+
+    int num_cell_facets = mesh::cell_num_entities(cell_types.back(),
                                                   mesh->topology().dim() - 1);
     const std::vector<int> c_offsets = a.coefficient_offsets();
     for (int i : a.integral_ids(IntegralType::interior_facet))
diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h
index 9111c62444a..dcf531dc24c 100644
--- a/cpp/dolfinx/fem/assemble_scalar_impl.h
+++ b/cpp/dolfinx/fem/assemble_scalar_impl.h
@@ -34,7 +34,7 @@ T assemble_cells(const mesh::Geometry<scalar_value_type_t<T>>& geometry,
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Create data structures used in assembly
@@ -74,7 +74,7 @@ T assemble_exterior_facets(
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Create data structures used in assembly
@@ -117,7 +117,7 @@ T assemble_interior_facets(
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Create data structures used in assembly
@@ -199,7 +199,12 @@ T assemble_scalar(
     const std::vector<std::uint8_t>& perms
         = mesh->topology().get_facet_permutations();
 
-    int num_cell_facets = mesh::cell_num_entities(mesh->topology().cell_type(),
+    auto cell_types = mesh->topology().cell_types();
+    if (cell_types.size() > 1)
+    {
+      throw std::runtime_error("MUltiple cell types in the assembler");
+    }
+    int num_cell_facets = mesh::cell_num_entities(cell_types.back(),
                                                   mesh->topology().dim() - 1);
     const std::vector<int> c_offsets = M.coefficient_offsets();
     for (int i : M.integral_ids(IntegralType::interior_facet))
diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h
index 3ebfa96c450..f61e918ce7a 100644
--- a/cpp/dolfinx/fem/assemble_vector_impl.h
+++ b/cpp/dolfinx/fem/assemble_vector_impl.h
@@ -60,7 +60,7 @@ void _lift_bc_cells(
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Data structures used in bc application
@@ -211,7 +211,7 @@ void _lift_bc_exterior_facets(
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Data structures used in bc application
@@ -319,7 +319,7 @@ void _lift_bc_interior_facets(
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Data structures used in assembly
@@ -498,7 +498,7 @@ void assemble_cells(
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // FIXME: Add proper interface for num_dofs
@@ -569,7 +569,7 @@ void assemble_exterior_facets(
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // FIXME: Add proper interface for num_dofs
@@ -636,7 +636,7 @@ void assemble_interior_facets(
 {
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Create data structures used in assembly
@@ -828,7 +828,12 @@ void lift_bc(std::span<T> b, const Form<T, U>& a,
     else
       get_perm = [](std::size_t) { return 0; };
 
-    int num_cell_facets = mesh::cell_num_entities(mesh->topology().cell_type(),
+    auto cell_types = mesh->topology().cell_types();
+    if (cell_types.size() > 1)
+    {
+      throw std::runtime_error("MUltiple cell types in the assembler");
+    }
+    int num_cell_facets = mesh::cell_num_entities(cell_types.back(),
                                                   mesh->topology().dim() - 1);
     for (int i : a.integral_ids(IntegralType::interior_facet))
     {
@@ -1027,7 +1032,13 @@ void assemble_vector(
     else
       get_perm = [](std::size_t) { return 0; };
 
-    int num_cell_facets = mesh::cell_num_entities(mesh->topology().cell_type(),
+    auto cell_types = mesh->topology().cell_types();
+    if (cell_types.size() > 1)
+    {
+      throw std::runtime_error("MUltiple cell types in the assembler");
+    }
+
+    int num_cell_facets = mesh::cell_num_entities(cell_types.back(),
                                                   mesh->topology().dim() - 1);
     for (int i : L.integral_ids(IntegralType::interior_facet))
     {
diff --git a/cpp/dolfinx/fem/discreteoperators.h b/cpp/dolfinx/fem/discreteoperators.h
index ce65a3119c9..78b3668940f 100644
--- a/cpp/dolfinx/fem/discreteoperators.h
+++ b/cpp/dolfinx/fem/discreteoperators.h
@@ -196,7 +196,10 @@ void interpolation_matrix(const FunctionSpace<U>& V0,
   const std::size_t value_size1 = e1->value_size() / bs1;
 
   // Get geometry data
-  const CoordinateElement& cmap = mesh->geometry().cmap();
+  auto cmaps = mesh->geometry().cmaps();
+  assert(cmaps.size() == 1);
+
+  const CoordinateElement& cmap = cmaps.back();
   const graph::AdjacencyList<std::int32_t>& x_dofmap
       = mesh->geometry().dofmap();
   const std::size_t num_dofs_g = cmap.dim();
diff --git a/cpp/dolfinx/fem/dofmapbuilder.cpp b/cpp/dolfinx/fem/dofmapbuilder.cpp
index be4fa7b54d8..f3e0791da19 100644
--- a/cpp/dolfinx/fem/dofmapbuilder.cpp
+++ b/cpp/dolfinx/fem/dofmapbuilder.cpp
@@ -128,8 +128,9 @@ reorder_owned(const graph::AdjacencyList<std::int32_t>& dofmap,
 /// {dimension, mesh entity index} for each local dof i}
 std::tuple<graph::AdjacencyList<std::int32_t>, std::vector<std::int64_t>,
            std::vector<std::pair<std::int8_t, std::int32_t>>>
-build_basic_dofmap(const mesh::Topology& topology,
-                   const fem::ElementDofLayout& element_dof_layout)
+build_basic_dofmap(
+    const mesh::Topology& topology,
+    const std::vector<fem::ElementDofLayout>& element_dof_layouts)
 {
   // Start timer for dofmap initialization
   common::Timer t0("Init dofmap from element dofmap");
@@ -137,13 +138,55 @@ build_basic_dofmap(const mesh::Topology& topology,
   // Topological dimension
   const int D = topology.dim();
 
+  // Checks for mixed topology
+  if (element_dof_layouts.size() != topology.cell_types().size())
+    throw std::runtime_error("Mixed topology mismatch: cell types");
+
+  // Should be 2*n+1 cell group offsets in Topology for n elements
+  if (element_dof_layouts.size() * 2 + 1
+      != topology.entity_group_offsets(D).size())
+  {
+    throw std::runtime_error("Mixed topology mismatch: groups");
+  }
+
+  // Mixed topology can only manage one dof (e.g. on vertex, or on edge, i.e. P1
+  // or P2) for now.
+  if (element_dof_layouts.size() > 1)
+  {
+    for (auto e : element_dof_layouts)
+    {
+      for (int d = 0; d <= D; ++d)
+      {
+
+        if (e.num_entity_dofs(d) > 1)
+          throw std::runtime_error("Mixed topology: dofmapbuilder does not yet "
+                                   "support elements with more than "
+                                   "one dof per entity dimension (P1/P2)");
+      }
+    }
+
+    std::array<int, 3> nd;
+    for (int d = 0; d < D; ++d)
+      nd[d] = element_dof_layouts[0].num_entity_dofs(d);
+    for (int d = 0; d < D; ++d)
+      for (std::size_t i = 1; i < element_dof_layouts.size(); ++i)
+      {
+        if (element_dof_layouts[i].num_entity_dofs(d) != nd[d])
+          throw std::runtime_error(
+              "Mixed topology: dofmapbuilder - incompatible elements "
+              "(different number of dofs on shared entity)");
+      }
+  }
+
   // Generate and number required mesh entities
   std::vector<std::int8_t> needs_entities(D + 1, false);
   std::vector<std::int32_t> num_mesh_entities_local(D + 1, 0),
       num_mesh_entities_global(D + 1, 0);
   for (int d = 0; d <= D; ++d)
   {
-    if (element_dof_layout.num_entity_dofs(d) > 0)
+    // FIXME: Mixed-topology - P2/Q2 - Q has dof on interior, P does not - need
+    // to check all elements here.
+    if (element_dof_layouts[0].num_entity_dofs(d) > 0)
     {
       if (!topology.connectivity(d, 0))
       {
@@ -177,34 +220,50 @@ build_basic_dofmap(const mesh::Topology& topology,
   }
 
   // Number of dofs on this process
+  // FIXME: This assumes that the number of dofs on the entities of each element
+  // must be the same. This must be true for entities which are shared
+  // (vertices, edges, facets) but may differ for cell.
   std::int32_t local_size(0), d(0);
   for (std::int32_t n : num_mesh_entities_local)
-    local_size += n * element_dof_layout.num_entity_dofs(d++);
-
-  // Number of dofs per cell
-  const int local_dim = element_dof_layout.num_dofs();
+    local_size += n * element_dof_layouts[0].num_entity_dofs(d++);
 
   // Allocate dofmap memory
   const int num_cells = topology.connectivity(D, 0)->num_nodes();
-  std::vector<std::int32_t> dofs(num_cells * local_dim);
-  std::vector<std::int32_t> cell_ptr(num_cells + 1, local_dim);
-  cell_ptr[0] = 0;
-  std::partial_sum(std::next(cell_ptr.begin(), 1), cell_ptr.end(),
-                   std::next(cell_ptr.begin(), 1));
+  const std::vector<int>& group_offsets = topology.entity_group_offsets(D);
+
+  std::vector<std::int32_t> cell_ptr;
+  cell_ptr.reserve(num_cells + 1);
+  cell_ptr.push_back(0);
+
+  const std::size_t nelem = element_dof_layouts.size();
+  // Go through elements twice: regular cells followed by ghost cells.
+  for (std::size_t i = 0; i < 2 * nelem; ++i)
+  {
+    // Number of dofs per cell for this element layout
+    const int local_dim = element_dof_layouts[i % nelem].num_dofs();
+    for (int j = group_offsets[i]; j < group_offsets[i + 1]; ++j)
+      cell_ptr.push_back(cell_ptr.back() + local_dim);
+  }
+  std::vector<std::int32_t> dofs(cell_ptr.back());
 
   // Allocate entity indices array
   std::vector<std::vector<int32_t>> entity_indices_local(D + 1);
   std::vector<std::vector<int64_t>> entity_indices_global(D + 1);
   for (int d = 0; d <= D; ++d)
   {
-    const int num_entities = mesh::cell_num_entities(topology.cell_type(), d);
-    entity_indices_local[d].resize(num_entities);
-    entity_indices_global[d].resize(num_entities);
+    for (auto cell_type : topology.cell_types())
+    {
+      const std::size_t num_entities = mesh::cell_num_entities(cell_type, d);
+      entity_indices_local[d].resize(
+          std::max(num_entities, entity_indices_local[d].size()));
+      entity_indices_global[d].resize(entity_indices_local[d].size());
+    }
   }
 
-  // Entity dofs on cell (dof = entity_dofs[dim][entity][index])
-  const std::vector<std::vector<std::vector<int>>>& entity_dofs
-      = element_dof_layout.entity_dofs_all();
+  // Entity dofs on cell (dof = entity_dofs[element][dim][entity][index])
+  std::vector<std::vector<std::vector<std::vector<int>>>> entity_dofs(nelem);
+  for (std::size_t i = 0; i < element_dof_layouts.size(); ++i)
+    entity_dofs[i] = element_dof_layouts[i].entity_dofs_all();
 
   // Storage for local-to-global map
   std::vector<std::int64_t> local_to_global(local_size);
@@ -212,60 +271,72 @@ build_basic_dofmap(const mesh::Topology& topology,
   // Dof -> (dim, entity index) marker
   std::vector<std::pair<std::int8_t, std::int32_t>> dof_entity(local_size);
 
-  // Loops over cells and build dofmaps from ElementDofmap
-  for (int c = 0; c < connectivity[0]->num_nodes(); ++c)
+  assert(group_offsets.back() == connectivity[0]->num_nodes());
+
+  // Loop over cells, group by group, and build dofmaps from respective
+  // ElementDofmap
+  for (std::size_t i = 0; i < 2 * nelem; ++i)
   {
-    // Get local (process) and global indices for each cell entity
-    for (int d = 0; d < D; ++d)
+    for (int c = group_offsets[i]; c < group_offsets[i + 1]; ++c)
     {
-      if (needs_entities[d])
+      // Get local (process) and global indices for each cell entity
+      for (int d = 0; d < D; ++d)
       {
-        auto entities = connectivity[d]->links(c);
-        for (std::size_t i = 0; i < entities.size(); ++i)
+        if (needs_entities[d])
         {
-          entity_indices_local[d][i] = entities[i];
-          entity_indices_global[d][i] = global_indices[d][entities[i]];
+          auto entities = connectivity[d]->links(c);
+          for (std::size_t i = 0; i < entities.size(); ++i)
+          {
+            entity_indices_local[d][i] = entities[i];
+            entity_indices_global[d][i] = global_indices[d][entities[i]];
+          }
         }
       }
-    }
 
-    // Handle cell index separately because cell.entities(D) doesn't work.
-    if (needs_entities[D])
-    {
-      entity_indices_global[D][0] = global_indices[D][c];
-      entity_indices_local[D][0] = c;
-    }
-
-    // Iterate over each topological dimension
-    std::int32_t offset_local = 0;
-    for (auto e_dofs_d = entity_dofs.begin(); e_dofs_d != entity_dofs.end();
-         ++e_dofs_d)
-    {
-      const std::size_t d = std::distance(entity_dofs.begin(), e_dofs_d);
+      // Handle cell index separately because cell.entities(D) doesn't work.
+      if (needs_entities[D])
+      {
+        entity_indices_global[D][0] = global_indices[D][c];
+        entity_indices_local[D][0] = c;
+      }
 
-      // Iterate over each entity of current dimension d
-      for (auto e_dofs = e_dofs_d->begin(); e_dofs != e_dofs_d->end(); ++e_dofs)
+      // Iterate over each topological dimension for this element (twice, once
+      // for regular, and later for ghosts).
+      const int elem = i % nelem;
+      std::int32_t offset_local = 0;
+      for (auto e_dofs_d = entity_dofs[elem].begin();
+           e_dofs_d != entity_dofs[elem].end(); ++e_dofs_d)
       {
-        // Get entity indices (local to cell, local to process, and
-        // global)
-        const std::size_t e = std::distance(e_dofs_d->begin(), e_dofs);
-        const std::int32_t e_index_local = entity_indices_local[d][e];
-
-        // Loop over dofs belong to entity e of dimension d (d, e)
-        // d: topological dimension
-        // e: local entity index
-        // dof_local: local index of dof at (d, e)
-        const std::size_t num_entity_dofs = e_dofs->size();
-        for (auto dof_local = e_dofs->begin(); dof_local != e_dofs->end();
-             ++dof_local)
+        const std::size_t d
+            = std::distance(entity_dofs[elem].begin(), e_dofs_d);
+
+        // Iterate over each entity of current dimension d
+        for (auto e_dofs = e_dofs_d->begin(); e_dofs != e_dofs_d->end();
+             ++e_dofs)
         {
-          const std::int32_t count = std::distance(e_dofs->begin(), dof_local);
-          const std::int32_t dof
-              = offset_local + num_entity_dofs * e_index_local + count;
-          dofs[cell_ptr[c] + *dof_local] = dof;
+          // Get entity indices (local to cell, local to process, and
+          // global)
+          const std::size_t e = std::distance(e_dofs_d->begin(), e_dofs);
+          const std::int32_t e_index_local = entity_indices_local[d][e];
+
+          // Loop over dofs belong to entity e of dimension d (d, e)
+          // d: topological dimension
+          // e: local entity index
+          // dof_local: local index of dof at (d, e)
+          const std::size_t num_entity_dofs = e_dofs->size();
+          for (auto dof_local = e_dofs->begin(); dof_local != e_dofs->end();
+               ++dof_local)
+          {
+            const std::int32_t count
+                = std::distance(e_dofs->begin(), dof_local);
+            const std::int32_t dof
+                = offset_local + num_entity_dofs * e_index_local + count;
+            dofs[cell_ptr[c] + *dof_local] = dof;
+          }
         }
+        offset_local
+            += entity_dofs[elem][d][0].size() * num_mesh_entities_local[d];
       }
-      offset_local += entity_dofs[d][0].size() * num_mesh_entities_local[d];
     }
   }
 
@@ -279,7 +350,8 @@ build_basic_dofmap(const mesh::Topology& topology,
     if (needs_entities[d])
     {
       // NOTE This assumes all entities have the same number of dofs
-      auto num_entity_dofs = entity_dofs[d][0].size();
+      // Mixed topology: probably OK for P1/Q1 but not higher...
+      auto num_entity_dofs = entity_dofs[0][d][0].size();
       for (std::int32_t e_index_local = 0;
            e_index_local < num_mesh_entities_local[d]; ++e_index_local)
       {
@@ -568,7 +640,7 @@ std::pair<std::vector<std::int64_t>, std::vector<int>> get_global_indices(
 std::tuple<common::IndexMap, int, graph::AdjacencyList<std::int32_t>>
 fem::build_dofmap_data(
     MPI_Comm comm, const mesh::Topology& topology,
-    const ElementDofLayout& element_dof_layout,
+    const std::vector<ElementDofLayout>& element_dof_layouts,
     const std::function<std::vector<int>(
         const graph::AdjacencyList<std::int32_t>&)>& reorder_fn)
 {
@@ -581,17 +653,17 @@ fem::build_dofmap_data(
   // pair {dimension, mesh entity index} giving the mesh entity that dof
   // i is associated with.
   const auto [node_graph0, local_to_global0, dof_entity0]
-      = build_basic_dofmap(topology, element_dof_layout);
+      = build_basic_dofmap(topology, element_dof_layouts);
 
   // Compute global dofmap offset
   std::int64_t offset = 0;
   for (int d = 0; d <= D; ++d)
   {
-    if (element_dof_layout.num_entity_dofs(d) > 0)
+    if (element_dof_layouts[0].num_entity_dofs(d) > 0)
     {
       assert(topology.index_map(d));
       offset += topology.index_map(d)->local_range()[0]
-                * element_dof_layout.num_entity_dofs(d);
+                * element_dof_layouts[0].num_entity_dofs(d);
     }
   }
 
@@ -624,13 +696,9 @@ fem::build_dofmap_data(
       dofmap[local_dim0 * cell + j] = new_node;
     }
   }
+  graph::AdjacencyList<std::int32_t> dmap(dofmap, node_graph0.offsets());
 
-  int dofs_per_cell
-      = dofmap.empty() ? 0 : dofmap.size() / node_graph0.num_nodes();
-  graph::AdjacencyList<std::int32_t> dmap
-      = graph::regular_adjacency_list(std::move(dofmap), dofs_per_cell);
-
-  return {std::move(index_map), element_dof_layout.block_size(),
+  return {std::move(index_map), element_dof_layouts[0].block_size(),
           std::move(dmap)};
 }
 //-----------------------------------------------------------------------------
diff --git a/cpp/dolfinx/fem/dofmapbuilder.h b/cpp/dolfinx/fem/dofmapbuilder.h
index 62f13ab9a19..db045f6b3ec 100644
--- a/cpp/dolfinx/fem/dofmapbuilder.h
+++ b/cpp/dolfinx/fem/dofmapbuilder.h
@@ -34,14 +34,14 @@ class ElementDofLayout;
 /// Build dofmap data for an element on a mesh topology
 /// @param[in] comm MPI communicator
 /// @param[in] topology The mesh topology
-/// @param[in] element_dof_layout The element dof layout for the
+/// @param[in] element_dof_layouts The element dof layouts for the
 /// function space
 /// @param[in] reorder_fn Graph reordering function that is applied to
 /// the dofmap
 /// @return The index map and local to global DOF data for the DOF map
 std::tuple<common::IndexMap, int, graph::AdjacencyList<std::int32_t>>
 build_dofmap_data(MPI_Comm comm, const mesh::Topology& topology,
-                  const ElementDofLayout& element_dof_layout,
+                  const std::vector<ElementDofLayout>& element_dof_layouts,
                   const std::function<std::vector<int>(
                       const graph::AdjacencyList<std::int32_t>&)>& reorder_fn);
 
diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h
index 85a190928b0..401117af6ad 100644
--- a/cpp/dolfinx/fem/interpolate.h
+++ b/cpp/dolfinx/fem/interpolate.h
@@ -46,7 +46,13 @@ std::vector<T> interpolation_coords(const fem::FiniteElement& element,
   const std::size_t gdim = geometry.dim();
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
   std::span<const T> x_g = geometry.x();
-  const CoordinateElement& cmap = geometry.cmap();
+
+  if (geometry.cmaps().size() > 1)
+  {
+    throw std::runtime_error("Mixed topology not supported");
+  }
+
+  const CoordinateElement& cmap = geometry.cmaps()[0];
   const std::size_t num_dofs_g = cmap.dim();
 
   // Get the interpolation points on the reference cells
@@ -444,7 +450,10 @@ void interpolate_nonmatching_maps(Function<T, U>& u1, const Function<T, U>& u0,
   const std::size_t value_size0 = element0->value_size() / bs0;
 
   // Get geometry data
-  const CoordinateElement& cmap = mesh->geometry().cmap();
+  if (mesh->geometry().cmaps().size() > 1)
+    throw std::runtime_error("Multiple cmaps");
+
+  const CoordinateElement& cmap = mesh->geometry().cmaps()[0];
   const graph::AdjacencyList<std::int32_t>& x_dofmap
       = mesh->geometry().dofmap();
   const std::size_t num_dofs_g = cmap.dim();
@@ -856,7 +865,9 @@ void interpolate(Function<T, U>& u, std::span<const T> f,
       throw std::runtime_error("Interpolation data has the wrong shape.");
 
     // Get coordinate map
-    const CoordinateElement& cmap = mesh->geometry().cmap();
+    if (mesh->geometry().cmaps().size() > 1)
+      throw std::runtime_error("Multiple cmaps");
+    const CoordinateElement& cmap = mesh->geometry().cmaps()[0];
 
     // Get geometry data
     const graph::AdjacencyList<std::int32_t>& x_dofmap
diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp
index 4a6574d0fa4..af8c858ada9 100644
--- a/cpp/dolfinx/fem/utils.cpp
+++ b/cpp/dolfinx/fem/utils.cpp
@@ -153,7 +153,7 @@ fem::create_dofmap(MPI_Comm comm, const ElementDofLayout& layout,
   }
 
   auto [_index_map, bs, dofmap]
-      = build_dofmap_data(comm, topology, layout, reorder_fn);
+      = build_dofmap_data(comm, topology, {layout}, reorder_fn);
   auto index_map = std::make_shared<common::IndexMap>(std::move(_index_map));
 
   // If the element's DOF transformations are permutations, permute the
diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h
index 5e331be0a1e..8c5b2694a23 100644
--- a/cpp/dolfinx/fem/utils.h
+++ b/cpp/dolfinx/fem/utils.h
@@ -716,13 +716,18 @@ create_functionspace(ufcx_function_space* (*fptr)(const char*),
   ufcx_finite_element* ufcx_element = space->finite_element;
   assert(ufcx_element);
 
-  if (space->geometry_degree != mesh->geometry().cmap().degree()
+  if (mesh->geometry().cmaps().size() > 1)
+  {
+    throw std::runtime_error("Not supported for Mixed Topology");
+  }
+
+  if (space->geometry_degree != mesh->geometry().cmaps()[0].degree()
       or static_cast<basix::cell::type>(space->geometry_basix_cell)
              != mesh::cell_type_to_basix_type(
-                 mesh->geometry().cmap().cell_shape())
+                 mesh->geometry().cmaps()[0].cell_shape())
       or static_cast<basix::element::lagrange_variant>(
              space->geometry_basix_variant)
-             != mesh->geometry().cmap().variant())
+             != mesh->geometry().cmaps()[0].variant())
   {
     throw std::runtime_error("UFL mesh and CoordinateElement do not match.");
   }
@@ -732,7 +737,7 @@ create_functionspace(ufcx_function_space* (*fptr)(const char*),
   ufcx_dofmap* ufcx_map = space->dofmap;
   assert(ufcx_map);
   ElementDofLayout layout
-      = create_element_dof_layout(*ufcx_map, mesh->topology().cell_type());
+      = create_element_dof_layout(*ufcx_map, mesh->topology().cell_types()[0]);
   return FunctionSpace(
       mesh, element,
       std::make_shared<DofMap>(create_dofmap(
@@ -801,12 +806,13 @@ void pack(std::span<T> coeffs, std::int32_t cell, int bs, std::span<const T> v,
 /// @private
 /// @brief  Concepts for function that returns cell index
 template <typename F>
-concept FetchCells = requires(F&& f, std::span<const std::int32_t> v) {
-  requires std::invocable<F, std::span<const std::int32_t>>;
-  {
-    f(v)
-  } -> std::convertible_to<std::int32_t>;
-};
+concept FetchCells
+    = requires(F&& f, std::span<const std::int32_t> v) {
+        requires std::invocable<F, std::span<const std::int32_t>>;
+        {
+          f(v)
+          } -> std::convertible_to<std::int32_t>;
+      };
 
 /// @brief Pack a single coefficient for a set of active entities.
 ///
diff --git a/cpp/dolfinx/geometry/gjk.h b/cpp/dolfinx/geometry/gjk.h
index 7043dc202ca..bee951ccaa4 100644
--- a/cpp/dolfinx/geometry/gjk.h
+++ b/cpp/dolfinx/geometry/gjk.h
@@ -183,7 +183,7 @@ nearest_simplex(std::span<const T> s)
     {
       if (f_inside[0]) // The origin is inside the tetrahedron
         return {std::vector<T>(s.begin(), s.end()), {0, 0, 0}};
-      else             // The origin projection P faces BCD
+      else // The origin projection P faces BCD
         return nearest_simplex<T>(s.template subspan<0, 3 * 3>());
     }
 
diff --git a/cpp/dolfinx/geometry/utils.h b/cpp/dolfinx/geometry/utils.h
index 780cfc9a9fb..c5062044737 100644
--- a/cpp/dolfinx/geometry/utils.h
+++ b/cpp/dolfinx/geometry/utils.h
@@ -38,6 +38,12 @@ std::vector<T> shortest_vector(const mesh::Mesh<T>& mesh, int dim,
 {
   const int tdim = mesh.topology().dim();
   const mesh::Geometry<T>& geometry = mesh.geometry();
+
+  if (geometry.cmaps().size() > 1)
+  {
+    throw std::runtime_error("Mixed topology not supported");
+  }
+
   std::span<const T> geom_dofs = geometry.x();
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
   std::vector<T> shortest_vectors(3 * entities.size());
@@ -84,7 +90,7 @@ std::vector<T> shortest_vector(const mesh::Mesh<T>& mesh, int dim,
       // Tabulate geometry dofs for the entity
       auto dofs = x_dofmap.links(c);
       const std::vector<int> entity_dofs
-          = geometry.cmap().create_dof_layout().entity_closure_dofs(
+          = geometry.cmaps()[0].create_dof_layout().entity_closure_dofs(
               dim, local_cell_entity);
       std::vector<T> nodes(3 * entity_dofs.size());
       for (std::size_t i = 0; i < entity_dofs.size(); i++)
@@ -492,6 +498,11 @@ std::int32_t compute_first_colliding_cell(const mesh::Mesh<T>& mesh,
   std::vector<std::int32_t> cell_candidates;
   impl::_compute_collisions_point<T>(tree, point, cell_candidates);
 
+  if (mesh.geometry().cmaps().size() > 1)
+  {
+    throw std::runtime_error("Mixed topology not supported");
+  }
+
   if (cell_candidates.empty())
     return -1;
   else
@@ -500,7 +511,7 @@ std::int32_t compute_first_colliding_cell(const mesh::Mesh<T>& mesh,
     const mesh::Geometry<T>& geometry = mesh.geometry();
     std::span<const T> geom_dofs = geometry.x();
     const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-    const std::size_t num_nodes = geometry.cmap().dim();
+    const std::size_t num_nodes = geometry.cmaps()[0].dim();
     std::vector<T> coordinate_dofs(num_nodes * 3);
     for (auto cell : cell_candidates)
     {
diff --git a/cpp/dolfinx/graph/AdjacencyList.h b/cpp/dolfinx/graph/AdjacencyList.h
index 5666ad3eca1..241e34e9847 100644
--- a/cpp/dolfinx/graph/AdjacencyList.h
+++ b/cpp/dolfinx/graph/AdjacencyList.h
@@ -172,10 +172,10 @@ AdjacencyList(T, U) -> AdjacencyList<typename T::value_type>;
 /// @return An adjacency list
 template <typename U>
   requires requires {
-    typename std::decay_t<U>::value_type;
-    requires std::convertible_to<
-        U, std::vector<typename std::decay_t<U>::value_type>>;
-  }
+             typename std::decay_t<U>::value_type;
+             requires std::convertible_to<
+                 U, std::vector<typename std::decay_t<U>::value_type>>;
+           }
 AdjacencyList<typename std::decay_t<U>::value_type>
 regular_adjacency_list(U&& data, int degree)
 {
diff --git a/cpp/dolfinx/io/ADIOS2Writers.cpp b/cpp/dolfinx/io/ADIOS2Writers.cpp
index 3c3546f6f35..4b9937f3fe9 100644
--- a/cpp/dolfinx/io/ADIOS2Writers.cpp
+++ b/cpp/dolfinx/io/ADIOS2Writers.cpp
@@ -100,7 +100,6 @@ adios2::Variable<T> define_variable(adios2::IO& io, std::string name,
     return io.DefineVariable<T>(name, shape, start, count);
 }
 //-----------------------------------------------------------------------------
-
 /// Extract name of functions and split into real and imaginary component
 template <typename T>
 std::vector<std::string>
@@ -178,7 +177,8 @@ std::vector<T> pack_function_data(const fem::Function<T, U>& u)
 
   // The Function and the mesh must have identical element_dof_layouts
   // (up to the block size)
-  assert(dofmap->element_dof_layout() == geometry.cmap().create_dof_layout());
+  assert(dofmap->element_dof_layout()
+         == geometry.cmaps()[0].create_dof_layout());
 
   const int tdim = topology.dim();
   auto cell_map = topology.index_map(tdim);
@@ -314,9 +314,9 @@ void fides_write_mesh(adios2::IO& io, adios2::Engine& engine,
   // cell, and compute 'VTK' connectivity
   const int tdim = topology.dim();
   const std::int32_t num_cells = topology.index_map(tdim)->size_local();
-  const int num_nodes = geometry.cmap().dim();
+  const int num_nodes = geometry.cmaps()[0].dim();
   const auto [cells, shape] = io::extract_vtk_connectivity(
-      mesh.geometry().dofmap(), mesh.topology().cell_type());
+      mesh.geometry().dofmap(), mesh.topology().cell_types()[0]);
 
   // "Put" topology data in the result in the ADIOS2 file
   adios2::Variable<std::int64_t> local_topology = define_variable<std::int64_t>(
@@ -337,9 +337,9 @@ void fides_initialize_mesh_attributes(adios2::IO& io, const mesh::Mesh<T>& mesh)
   const mesh::Topology& topology = mesh.topology();
 
   // Check that mesh is first order mesh
-  const int num_dofs_g = geometry.cmap().dim();
+  const int num_dofs_g = geometry.cmaps()[0].dim();
   const int num_vertices_per_cell
-      = mesh::cell_num_entities(topology.cell_type(), 0);
+      = mesh::cell_num_entities(topology.cell_types()[0], 0);
   if (num_dofs_g != num_vertices_per_cell)
     throw std::runtime_error("Fides only supports lowest-order meshes.");
 
@@ -353,7 +353,7 @@ void fides_initialize_mesh_attributes(adios2::IO& io, const mesh::Mesh<T>& mesh)
   define_attribute<std::string>(io, "Fides_Connecticity_Variable",
                                 "connectivity");
 
-  std::string cell_type = to_fides_cell(topology.cell_type());
+  std::string cell_type = to_fides_cell(topology.cell_types()[0]);
   define_attribute<std::string>(io, "Fides_Cell_Type", cell_type);
 
   define_attribute<std::string>(io, "Fides_Time_Variable", "step");
@@ -579,8 +579,11 @@ FidesWriter::FidesWriter(MPI_Comm comm, const std::filesystem::path& filename,
   }
 
   // Check that all functions are first order Lagrange
-  int num_vertices_per_cell
-      = mesh::cell_num_entities(mesh->topology().cell_type(), 0);
+  auto cell_types = mesh->topology().cell_types();
+  if (cell_types.size() > 1)
+    throw std::runtime_error("Multiple cell types in IO.");
+
+  int num_vertices_per_cell = mesh::cell_num_entities(cell_types.back(), 0);
   for (auto& v : _u)
   {
     std::visit(
@@ -630,4 +633,4 @@ void FidesWriter::write(double t)
   _engine->EndStep();
 }
 
-#endif
\ No newline at end of file
+#endif
diff --git a/cpp/dolfinx/io/ADIOS2Writers.h b/cpp/dolfinx/io/ADIOS2Writers.h
index eab1a59fe70..4380a7b5ed0 100644
--- a/cpp/dolfinx/io/ADIOS2Writers.h
+++ b/cpp/dolfinx/io/ADIOS2Writers.h
@@ -374,7 +374,7 @@ void vtx_write_mesh(adios2::IO& io, adios2::Engine& engine,
   engine.Put<std::uint32_t>(vertices, num_vertices);
 
   const auto [vtkcells, shape] = io::extract_vtk_connectivity(
-      mesh.geometry().dofmap(), mesh.topology().cell_type());
+      mesh.geometry().dofmap(), mesh.topology().cell_types()[0]);
 
   // Add cell metadata
   const int tdim = topology.dim();
@@ -385,7 +385,8 @@ void vtx_write_mesh(adios2::IO& io, adios2::Engine& engine,
   adios2::Variable<std::uint32_t> celltype_variable
       = define_variable<std::uint32_t>(io, "types");
   engine.Put<std::uint32_t>(
-      celltype_variable, cells::get_vtk_cell_type(topology.cell_type(), tdim));
+      celltype_variable,
+      cells::get_vtk_cell_type(topology.cell_types()[0], tdim));
 
   // Pack mesh 'nodes'. Output is written as [N0, v0_0,...., v0_N0, N1,
   // v1_0,...., v1_N1,....], where N is the number of cell nodes and v0,
@@ -470,7 +471,8 @@ void vtx_write_mesh_from_space(adios2::IO& io, adios2::Engine& engine,
   engine.Put<std::uint32_t>(vertices, num_dofs);
   engine.Put<std::uint32_t>(elements, vtkshape[0]);
   engine.Put<std::uint32_t>(
-      cell_type, cells::get_vtk_cell_type(mesh->topology().cell_type(), tdim));
+      cell_type,
+      cells::get_vtk_cell_type(mesh->topology().cell_types()[0], tdim));
   engine.Put<T>(local_geometry, x.data());
   engine.Put<std::int64_t>(local_topology, cells.data());
 
diff --git a/cpp/dolfinx/io/VTKFile.cpp b/cpp/dolfinx/io/VTKFile.cpp
index 8e4c6bd7477..1bfd840d6b1 100644
--- a/cpp/dolfinx/io/VTKFile.cpp
+++ b/cpp/dolfinx/io/VTKFile.cpp
@@ -422,7 +422,7 @@ void write_function(
   {
     std::vector<std::int64_t> tmp;
     std::tie(tmp, cshape) = io::extract_vtk_connectivity(
-        mesh0->geometry().dofmap(), mesh0->topology().cell_type());
+        mesh0->geometry().dofmap(), mesh0->topology().cell_types()[0]);
     cells.assign(tmp.begin(), tmp.end());
     const mesh::Geometry<double>& geometry = mesh0->geometry();
     x.assign(geometry.x().begin(), geometry.x().end());
@@ -442,10 +442,17 @@ void write_function(
   piece_node.append_attribute("NumberOfPoints") = xshape[0];
   piece_node.append_attribute("NumberOfCells") = cshape[0];
 
+  // FIXME
+  auto cell_types = mesh0->topology().cell_types();
+  if (cell_types.size() > 1)
+  {
+    throw std::runtime_error("Multiple cell types in IO");
+  }
+
   // Add mesh data to "Piece" node
   int tdim = mesh0->topology().dim();
   add_mesh(x, xshape, x_id, x_ghost, cells, cshape,
-           *mesh0->topology().index_map(tdim), mesh0->topology().cell_type(),
+           *mesh0->topology().index_map(tdim), cell_types.back(),
            mesh0->topology().dim(), piece_node);
 
   // FIXME: is this actually setting the first?
@@ -772,14 +779,20 @@ void io::VTKFile::write(const mesh::Mesh<double>& mesh, double time)
   piece_node.append_attribute("NumberOfPoints") = num_points;
   piece_node.append_attribute("NumberOfCells") = num_cells;
 
+  auto cell_types = topology.cell_types();
+  if (cell_types.size() > 1)
+  {
+    throw std::runtime_error("Multiple cell types in IO");
+  }
+
   // Add mesh data to "Piece" node
-  const auto [cells, cshape] = extract_vtk_connectivity(
-      mesh.geometry().dofmap(), mesh.topology().cell_type());
+  const auto [cells, cshape]
+      = extract_vtk_connectivity(mesh.geometry().dofmap(), cell_types[0]);
   std::array<std::size_t, 2> xshape = {geometry.x().size() / 3, 3};
   std::vector<std::uint8_t> x_ghost(xshape[0], 0);
   std::fill(std::next(x_ghost.begin(), xmap->size_local()), x_ghost.end(), 1);
   add_mesh(geometry.x(), xshape, geometry.input_global_indices(), x_ghost,
-           cells, cshape, *topology.index_map(tdim), topology.cell_type(),
+           cells, cshape, *topology.index_map(tdim), cell_types[0],
            topology.dim(), piece_node);
 
   // Create filepath for a .vtu file
diff --git a/cpp/dolfinx/io/XDMFFile.cpp b/cpp/dolfinx/io/XDMFFile.cpp
index 9e8738f24b8..78c8480e32d 100644
--- a/cpp/dolfinx/io/XDMFFile.cpp
+++ b/cpp/dolfinx/io/XDMFFile.cpp
@@ -240,9 +240,8 @@ mesh::Mesh<double> XDMFFile::read_mesh(const fem::CoordinateElement& element,
 
   graph::AdjacencyList<std::int64_t> cells_adj(std::move(cells),
                                                std::move(offset));
-
   mesh::Mesh<double> mesh
-      = mesh::create_mesh(_comm.comm(), cells_adj, element, x, xshape, mode);
+      = mesh::create_mesh(_comm.comm(), cells_adj, {element}, x, xshape, mode);
   mesh.name = name;
   return mesh;
 }
@@ -349,10 +348,13 @@ XDMFFile::read_meshtags(std::shared_ptr<const mesh::Mesh<double>> mesh,
       entities_values = xdmf_utils::distribute_entity_data(
           *mesh, mesh::cell_dim(cell_type), entities1, values);
 
+  auto cell_types = mesh->topology().cell_types();
+  if (cell_types.size() > 1)
+    throw std::runtime_error("cell type IO");
+
   LOG(INFO) << "XDMF create meshtags";
   const std::size_t num_vertices_per_entity = mesh::cell_num_entities(
-      mesh::cell_entity_type(mesh->topology().cell_type(),
-                             mesh::cell_dim(cell_type), 0),
+      mesh::cell_entity_type(cell_types.back(), mesh::cell_dim(cell_type), 0),
       0);
   const graph::AdjacencyList<std::int32_t> entities_adj
       = graph::regular_adjacency_list(std::move(entities_values.first),
diff --git a/cpp/dolfinx/io/vtk_utils.h b/cpp/dolfinx/io/vtk_utils.h
index 113ca7d6dfb..aa5c5712852 100644
--- a/cpp/dolfinx/io/vtk_utils.h
+++ b/cpp/dolfinx/io/vtk_utils.h
@@ -82,7 +82,7 @@ tabulate_lagrange_dof_coordinates(const fem::FunctionSpace<T>& V)
   // Get the dof coordinates on the reference element and the  mesh
   // coordinate map
   const auto [X, Xshape] = element->interpolation_points();
-  const fem::CoordinateElement& cmap = mesh->geometry().cmap();
+  const fem::CoordinateElement& cmap = mesh->geometry().cmaps()[0];
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& dofmap_x
@@ -204,7 +204,7 @@ vtk_mesh_from_space(const fem::FunctionSpace<T>& V)
   const std::uint32_t num_nodes
       = V.element()->space_dimension() / element_block_size;
   const std::vector<std::uint8_t> vtkmap = io::cells::transpose(
-      io::cells::perm_vtk(mesh->topology().cell_type(), num_nodes));
+      io::cells::perm_vtk(mesh->topology().cell_types()[0], num_nodes));
 
   // Extract topology for all local cells as
   // [v0_0, ...., v0_N0, v1_0, ...., v1_N1, ....]
diff --git a/cpp/dolfinx/io/xdmf_mesh.cpp b/cpp/dolfinx/io/xdmf_mesh.cpp
index 5307dee150a..6e8b829f9d4 100644
--- a/cpp/dolfinx/io/xdmf_mesh.cpp
+++ b/cpp/dolfinx/io/xdmf_mesh.cpp
@@ -29,15 +29,23 @@ void xdmf_mesh::add_topology_data(MPI_Comm comm, pugi::xml_node& xml_node,
 
   const int tdim = topology.dim();
 
-  if (tdim == 2 and topology.cell_type() == mesh::CellType::prism)
+  // FIXME
+  auto cell_types = topology.cell_types();
+  if (cell_types.size() > 1)
+    throw std::runtime_error("cell type IO");
+
+  if (tdim == 2 and cell_types.back() == mesh::CellType::prism)
     throw std::runtime_error("More work needed for prism cell");
 
   // Get entity 'cell' type
   const mesh::CellType entity_cell_type
-      = mesh::cell_entity_type(topology.cell_type(), dim, 0);
+      = mesh::cell_entity_type(cell_types.back(), dim, 0);
 
+  if (geometry.cmaps().size() > 1)
+    throw std::runtime_error(
+        "XDMF I/O with multiple geometry maps not implemented.");
   const fem::ElementDofLayout cmap_dof_layout
-      = geometry.cmap().create_dof_layout();
+      = geometry.cmaps()[0].create_dof_layout();
 
   // Get number of nodes per entity
   const int num_nodes_per_entity = cmap_dof_layout.num_entity_closure_dofs(dim);
@@ -90,9 +98,13 @@ void xdmf_mesh::add_topology_data(MPI_Comm comm, pugi::xml_node& xml_node,
     if (!c_to_e)
       throw std::runtime_error("Mesh is missing cell-entity connectivity.");
 
+    auto cell_types = topology.cell_types();
+    if (cell_types.size() > 1)
+      throw std::runtime_error("cell type IO");
+
     // Tabulate geometry dofs for local entities
     std::vector<std::vector<int>> entity_dofs;
-    for (int e = 0; e < mesh::cell_num_entities(topology.cell_type(), dim); ++e)
+    for (int e = 0; e < mesh::cell_num_entities(cell_types.back(), dim); ++e)
       entity_dofs.push_back(cmap_dof_layout.entity_closure_dofs(dim, e));
 
     for (std::int32_t e : entities)
diff --git a/cpp/dolfinx/io/xdmf_utils.cpp b/cpp/dolfinx/io/xdmf_utils.cpp
index c0357088665..7b6e2cfa982 100644
--- a/cpp/dolfinx/io/xdmf_utils.cpp
+++ b/cpp/dolfinx/io/xdmf_utils.cpp
@@ -94,7 +94,13 @@ compute_point_values(const fem::Function<T, double>& u)
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap
       = mesh->geometry().dofmap();
-  const int num_dofs_g = mesh->geometry().cmap().dim();
+
+  if (mesh->geometry().cmaps().size() > 1)
+  {
+    throw std::runtime_error(
+        "XDMF I/O with multiple geometry maps not implemented.");
+  }
+  const int num_dofs_g = mesh->geometry().cmaps()[0].dim();
 
   // Interpolate point values on each cell (using last computed value if
   // not continuous, e.g. discontinuous Galerkin methods)
@@ -451,16 +457,25 @@ xdmf_utils::distribute_entity_data(const mesh::Mesh<double>& mesh,
 {
   LOG(INFO) << "XDMF distribute entity data";
 
+  auto cell_types = mesh.topology().cell_types();
+  if (cell_types.size() > 1)
+    throw std::runtime_error("cell type IO");
+
   // Use ElementDofLayout of the cell to get vertex dof indices (local
   // to a cell), i.e. build a map from local vertex index to associated
   // local dof index
   std::vector<int> cell_vertex_dofs;
   {
     // Get layout of dofs on 0th cell entity of dimension entity_dim
+
+    if (mesh.geometry().cmaps().size() > 1)
+    {
+      throw std::runtime_error(
+          "XDMF I/O with multiple geometry maps not implemented.");
+    }
     const fem::ElementDofLayout cmap_dof_layout
-        = mesh.geometry().cmap().create_dof_layout();
-    for (int i = 0; i < mesh::cell_num_entities(mesh.topology().cell_type(), 0);
-         ++i)
+        = mesh.geometry().cmaps()[0].create_dof_layout();
+    for (int i = 0; i < mesh::cell_num_entities(cell_types.back(), 0); ++i)
     {
       const std::vector<int>& local_index = cmap_dof_layout.entity_dofs(0, i);
       assert(local_index.size() == 1);
@@ -526,14 +541,18 @@ xdmf_utils::distribute_entity_data(const mesh::Mesh<double>& mesh,
     const int comm_size = dolfinx::MPI::size(comm);
     const std::int64_t num_nodes_g = mesh.geometry().index_map()->size_global();
 
+    auto cell_types = mesh.topology().cell_types();
+    if (cell_types.size() > 1)
+      throw std::runtime_error("cell type IO");
+
     const std::size_t num_vert_per_entity = mesh::cell_num_entities(
-        mesh::cell_entity_type(mesh.topology().cell_type(), entity_dim, 0), 0);
+        mesh::cell_entity_type(cell_types.back(), entity_dim, 0), 0);
     auto c_to_v = mesh.topology().connectivity(mesh.topology().dim(), 0);
     if (!c_to_v)
       throw std::runtime_error("Missing cell-vertex connectivity.");
 
     const fem::ElementDofLayout cmap_dof_layout
-        = mesh.geometry().cmap().create_dof_layout();
+        = mesh.geometry().cmaps()[0].create_dof_layout();
     const std::vector<int> entity_layout
         = cmap_dof_layout.entity_closure_dofs(entity_dim, 0);
 
@@ -612,8 +631,12 @@ xdmf_utils::distribute_entity_data(const mesh::Mesh<double>& mesh,
     const MPI_Comm comm = mesh.comm();
     const int comm_size = dolfinx::MPI::size(comm);
 
+    auto cell_types = mesh.topology().cell_types();
+    if (cell_types.size() > 1)
+      throw std::runtime_error("cell type IO");
+
     const std::size_t num_vert_per_entity = mesh::cell_num_entities(
-        mesh::cell_entity_type(mesh.topology().cell_type(), entity_dim, 0), 0);
+        mesh::cell_entity_type(cell_types.back(), entity_dim, 0), 0);
 
     // Build map from global node index to ranks that have the node
     std::multimap<std::int64_t, int> node_to_rank;
@@ -686,8 +709,12 @@ xdmf_utils::distribute_entity_data(const mesh::Mesh<double>& mesh,
     // Build map from input global indices to local vertex numbers
     LOG(INFO) << "XDMF build map";
 
+    auto cell_types = mesh.topology().cell_types();
+    if (cell_types.size() > 1)
+      throw std::runtime_error("cell type IO");
+
     const std::size_t num_vert_per_entity = mesh::cell_num_entities(
-        mesh::cell_entity_type(mesh.topology().cell_type(), entity_dim, 0), 0);
+        mesh::cell_entity_type(cell_types.back(), entity_dim, 0), 0);
     auto c_to_v = mesh.topology().connectivity(mesh.topology().dim(), 0);
     if (!c_to_v)
       throw std::runtime_error("Missing cell-vertex connectivity.");
diff --git a/cpp/dolfinx/mesh/Geometry.h b/cpp/dolfinx/mesh/Geometry.h
index 67da4737876..4ce46983d21 100644
--- a/cpp/dolfinx/mesh/Geometry.h
+++ b/cpp/dolfinx/mesh/Geometry.h
@@ -38,8 +38,8 @@ class Geometry
   ///
   /// @param[in] index_map Index map associated with the geometry dofmap
   /// @param[in] dofmap The geometry (point) dofmap. For a cell, it
-  /// gives the position in the point array of each local geometry node.
-  /// @param[in] element The element that describes the cell geometry map.
+  /// gives the position in the point array of each local geometry node
+  /// @param[in] elements The elements that describes the cell geometry maps
   /// @param[in] x The point coordinates. The shape is `(num_points, 3)`
   /// and the storage is row-major.
   /// @param[in] dim The geometric dimension (`0 < dim <= 3`).
@@ -50,10 +50,10 @@ class Geometry
             std::convertible_to<std::vector<T>> V,
             std::convertible_to<std::vector<std::int64_t>> W>
   Geometry(std::shared_ptr<const common::IndexMap> index_map, U&& dofmap,
-           const fem::CoordinateElement& element, V&& x, int dim,
+           const std::vector<fem::CoordinateElement>& elements, V&& x, int dim,
            W&& input_global_indices)
       : _dim(dim), _dofmap(std::forward<U>(dofmap)), _index_map(index_map),
-        _cmap(element), _x(std::forward<V>(x)),
+        _cmaps(elements), _x(std::forward<V>(x)),
         _input_global_indices(std::forward<W>(input_global_indices))
   {
     assert(_x.size() % 3 == 0);
@@ -80,7 +80,7 @@ class Geometry
   template <typename U>
   Geometry<U> astype() const
   {
-    return Geometry<U>(_index_map, _dofmap, _cmap,
+    return Geometry<U>(_index_map, _dofmap, _cmaps,
                        std::vector<U>(_x.begin(), _x.end()), _dim,
                        _input_global_indices);
   }
@@ -110,10 +110,10 @@ class Geometry
   /// (num_points, 3)
   std::span<T> x() { return _x; }
 
-  /// @brief The element that describes the geometry map.
+  /// @brief The elements that describes the geometry maps.
   ///
-  /// @return The coordinate/geometry element
-  const fem::CoordinateElement& cmap() const { return _cmap; }
+  /// @return The coordinate/geometry elements
+  const std::vector<fem::CoordinateElement>& cmaps() const { return _cmaps; }
 
   /// Global user indices
   const std::vector<std::int64_t>& input_global_indices() const
@@ -131,8 +131,8 @@ class Geometry
   // IndexMap for geometry 'dofmap'
   std::shared_ptr<const common::IndexMap> _index_map;
 
-  // The coordinate element
-  fem::CoordinateElement _cmap;
+  // The coordinate elements
+  std::vector<fem::CoordinateElement> _cmaps;
 
   // Coordinates for all points stored as a contiguous array (row-major,
   // column size = 3)
@@ -150,7 +150,7 @@ class Geometry
 ///
 /// @param[in] comm The MPI communicator to build the Geometry on
 /// @param[in] topology The mesh topology
-/// @param[in] element The element that defines the geometry map for
+/// @param[in] elements The elements that defines the geometry map for
 /// each cell
 /// @param[in] cell_nodes The mesh cells, including higher-order
 /// geometry 'nodes'
@@ -165,7 +165,7 @@ class Geometry
 template <typename U>
 mesh::Geometry<typename std::remove_reference_t<typename U::value_type>>
 create_geometry(MPI_Comm comm, const Topology& topology,
-                const fem::CoordinateElement& element,
+                const std::vector<fem::CoordinateElement>& elements,
                 const graph::AdjacencyList<std::int64_t>& cell_nodes,
                 const U& x, int dim,
                 const std::function<std::vector<int>(
@@ -175,22 +175,28 @@ create_geometry(MPI_Comm comm, const Topology& topology,
   // TODO: make sure required entities are initialised, or extend
   // fem::build_dofmap_data
 
+  std::vector<fem::ElementDofLayout> dof_layouts;
+  for (auto e : elements)
+    dof_layouts.push_back(e.create_dof_layout());
+
   //  Build 'geometry' dofmap on the topology
-  auto [_dof_index_map, bs, dofmap] = fem::build_dofmap_data(
-      comm, topology, element.create_dof_layout(), reorder_fn);
+  auto [_dof_index_map, bs, dofmap]
+      = fem::build_dofmap_data(comm, topology, dof_layouts, reorder_fn);
   auto dof_index_map
       = std::make_shared<common::IndexMap>(std::move(_dof_index_map));
 
   // If the mesh has higher order geometry, permute the dofmap
-  if (element.needs_dof_permutations())
+  if (elements[0].needs_dof_permutations())
   {
+    if (elements.size() > 1)
+      throw std::runtime_error("Unsupported for Mixed Topology");
     const int D = topology.dim();
     const int num_cells = topology.connectivity(D, 0)->num_nodes();
     const std::vector<std::uint32_t>& cell_info
         = topology.get_cell_permutation_info();
 
     for (std::int32_t cell = 0; cell < num_cells; ++cell)
-      element.unpermute_dofs(dofmap.links(cell), cell_info[cell]);
+      elements[0].unpermute_dofs(dofmap.links(cell), cell_info[cell]);
   }
 
   auto remap_data
@@ -240,7 +246,7 @@ create_geometry(MPI_Comm comm, const Topology& topology,
   }
 
   return Geometry<typename std::remove_reference_t<typename U::value_type>>(
-      dof_index_map, std::move(dofmap), element, std::move(xg), dim,
+      dof_index_map, std::move(dofmap), elements, std::move(xg), dim,
       std::move(igi));
 }
 
@@ -257,9 +263,14 @@ std::pair<mesh::Geometry<T>, std::vector<int32_t>>
 create_subgeometry(const Topology& topology, const Geometry<T>& geometry,
                    int dim, std::span<const std::int32_t> subentity_to_entity)
 {
+  if (geometry.cmaps().size() > 1)
+  {
+    throw std::runtime_error("Mixed topology not supported");
+  }
+
   // Get the geometry dofs in the sub-geometry based on the entities in
   // sub-geometry
-  const fem::ElementDofLayout layout = geometry.cmap().create_dof_layout();
+  const fem::ElementDofLayout layout = geometry.cmaps()[0].create_dof_layout();
   // NOTE: Unclear what this return for prisms
   const std::size_t num_entity_dofs = layout.num_entity_closure_dofs(dim);
 
@@ -355,9 +366,10 @@ create_subgeometry(const Topology& topology, const Geometry<T>& geometry,
 
   // Create sub-geometry coordinate element
   CellType sub_coord_cell
-      = cell_entity_type(geometry.cmap().cell_shape(), dim, 0);
-  fem::CoordinateElement sub_coord_ele(sub_coord_cell, geometry.cmap().degree(),
-                                       geometry.cmap().variant());
+      = cell_entity_type(geometry.cmaps()[0].cell_shape(), dim, 0);
+  fem::CoordinateElement sub_coord_ele(sub_coord_cell,
+                                       geometry.cmaps()[0].degree(),
+                                       geometry.cmaps()[0].variant());
 
   // Sub-geometry input_global_indices
   // TODO: Check this
@@ -370,7 +382,7 @@ create_subgeometry(const Topology& topology, const Geometry<T>& geometry,
 
   // Create geometry
   return {Geometry<T>(sub_x_dof_index_map, std::move(sub_x_dofmap),
-                      sub_coord_ele, std::move(sub_x), geometry.dim(),
+                      {sub_coord_ele}, std::move(sub_x), geometry.dim(),
                       std::move(sub_igi)),
           std::move(subx_to_x_dofmap)};
 }
diff --git a/cpp/dolfinx/mesh/Topology.cpp b/cpp/dolfinx/mesh/Topology.cpp
index 583ab8838d5..8eb3b64d336 100644
--- a/cpp/dolfinx/mesh/Topology.cpp
+++ b/cpp/dolfinx/mesh/Topology.cpp
@@ -717,18 +717,29 @@ graph::AdjacencyList<std::int32_t> convert_to_local_indexing(
 } // namespace
 
 //-----------------------------------------------------------------------------
-Topology::Topology(MPI_Comm comm, CellType type)
-    : _comm(comm), _cell_type(type),
+Topology::Topology(MPI_Comm comm, std::vector<CellType> types)
+    : _comm(comm), _cell_types(types),
       _connectivity(
-          cell_dim(type) + 1,
+          cell_dim(types[0]) + 1,
           std::vector<std::shared_ptr<graph::AdjacencyList<std::int32_t>>>(
-              cell_dim(type) + 1))
+              cell_dim(types[0]) + 1))
 {
   // Do nothing
 }
 //-----------------------------------------------------------------------------
 int Topology::dim() const noexcept { return _connectivity.size() - 1; }
 //-----------------------------------------------------------------------------
+void Topology::set_entity_group_offsets(
+    int dim, const std::vector<std::int32_t>& offsets)
+{
+  _entity_group_offsets[dim] = offsets;
+}
+//-----------------------------------------------------------------------------
+const std::vector<std::int32_t>& Topology::entity_group_offsets(int dim) const
+{
+  return _entity_group_offsets[dim];
+}
+//-----------------------------------------------------------------------------
 void Topology::set_index_map(int dim,
                              std::shared_ptr<const common::IndexMap> map)
 {
@@ -873,30 +884,39 @@ const std::vector<std::int32_t>& Topology::interprocess_facets() const
   return _interprocess_facets;
 }
 //-----------------------------------------------------------------------------
-mesh::CellType Topology::cell_type() const noexcept { return _cell_type; }
+std::vector<mesh::CellType> Topology::cell_types() const noexcept
+{
+  return _cell_types;
+}
 //-----------------------------------------------------------------------------
 MPI_Comm Topology::comm() const { return _comm.comm(); }
 //-----------------------------------------------------------------------------
 Topology mesh::create_topology(
     MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& cells,
     std::span<const std::int64_t> original_cell_index,
-    std::span<const int> ghost_owners, const CellType& cell_type,
+    std::span<const int> ghost_owners, const std::vector<CellType>& cell_type,
+    const std::vector<std::int32_t>& cell_group_offsets,
     std::span<const std::int64_t> boundary_vertices)
 {
   common::Timer timer("Topology: create");
 
   LOG(INFO) << "Create topology";
-  if (cells.num_nodes() > 0
-      and cells.num_links(0) != num_cell_vertices(cell_type))
+  for (std::size_t i = 0; i < cell_type.size(); i++)
   {
-    throw std::runtime_error(
-        "Inconsistent number of cell vertices. Got "
-        + std::to_string(cells.num_links(0)) + ", expected "
-        + std::to_string(num_cell_vertices(cell_type)) + ".");
+    std::int32_t offset = cell_group_offsets[i];
+    int num_vertices = num_cell_vertices(cell_type[i]);
+    if (cells.num_nodes() > 0 and cells.num_links(offset) != num_vertices)
+      throw std::runtime_error("Inconsistent number of cell vertices. Got "
+                               + std::to_string(cells.num_links(offset))
+                               + ", expected " + std::to_string(num_vertices)
+                               + ".");
   }
 
   const std::int32_t num_local_cells = cells.num_nodes() - ghost_owners.size();
 
+  if (num_local_cells != cell_group_offsets[cell_type.size()])
+    throw std::runtime_error("Inconsistent offset or ghost number.");
+
   // Create sets of owned and unowned vertices from the cell ownership
   // and the list of boundary vertices
   auto [owned_vertices, unowned_vertices]
@@ -1130,6 +1150,8 @@ Topology mesh::create_topology(
       original_cell_index.begin(),
       std::next(original_cell_index.begin(), num_local_nodes));
 
+  topology.set_entity_group_offsets(tdim, cell_group_offsets);
+
   return topology;
 }
 //-----------------------------------------------------------------------------
@@ -1199,7 +1221,8 @@ mesh::create_subtopology(const Topology& topology, int dim,
       submap0->size_local() + submap0->num_ghosts());
 
   // Sub-topology entity to vertex connectivity
-  const CellType entity_type = cell_entity_type(topology.cell_type(), dim, 0);
+  const CellType entity_type
+      = cell_entity_type(topology.cell_types()[0], dim, 0);
   int num_vertices_per_entity = cell_num_entities(entity_type, 0);
   auto e_to_v = topology.connectivity(dim, 0);
   assert(e_to_v);
@@ -1232,11 +1255,15 @@ mesh::create_subtopology(const Topology& topology, int dim,
       std::move(sub_e_to_v_vec), std::move(sub_e_to_v_offsets));
 
   // Create sub-topology
-  Topology subtopology(topology.comm(), entity_type);
+  Topology subtopology(topology.comm(), {entity_type});
   subtopology.set_index_map(0, submap0);
   subtopology.set_index_map(dim, submap);
   subtopology.set_connectivity(sub_v_to_v, 0, 0);
   subtopology.set_connectivity(sub_e_to_v, dim, 0);
+  // Set groups (only one cell type)
+  subtopology.set_entity_group_offsets(
+      dim,
+      {0, submap->size_local(), submap->size_local() + submap->num_ghosts()});
 
   return {std::move(subtopology), std::move(subentities),
           std::move(subvertices0)};
@@ -1259,8 +1286,12 @@ mesh::entities_to_index(const Topology& topology, int dim,
   auto e_to_v = topology.connectivity(dim, 0);
   assert(e_to_v);
 
+  auto cell_types = topology.cell_types();
+  if (cell_types.size() > 1)
+    throw std::runtime_error("multiple cell types entities_to_index");
+
   const int num_vertices_per_entity
-      = cell_num_entities(cell_entity_type(topology.cell_type(), dim, 0), 0);
+      = cell_num_entities(cell_entity_type(cell_types.back(), dim, 0), 0);
 
   // Build map from ordered local vertex indices (key) to entity index
   // (value)
@@ -1299,4 +1330,4 @@ mesh::entities_to_index(const Topology& topology, int dim,
 
   return indices;
 }
-//-----------------------------------------------------------------------------
\ No newline at end of file
+//-----------------------------------------------------------------------------
diff --git a/cpp/dolfinx/mesh/Topology.h b/cpp/dolfinx/mesh/Topology.h
index dfc7b09bb27..27bcf025b65 100644
--- a/cpp/dolfinx/mesh/Topology.h
+++ b/cpp/dolfinx/mesh/Topology.h
@@ -44,7 +44,7 @@ class Topology
 {
 public:
   /// Create empty mesh topology
-  Topology(MPI_Comm comm, CellType type);
+  Topology(MPI_Comm comm, std::vector<CellType> type);
 
   /// Copy constructor
   Topology(const Topology& topology) = default;
@@ -64,6 +64,26 @@ class Topology
   /// @brief Return the topological dimension of the mesh.
   int dim() const noexcept;
 
+  /// @brief Set the offsets for each group of entities of a particular
+  /// dimension. See `entity_group_offsets`.
+  ///
+  /// @param dim Dimension of the entities
+  /// @param offsets The offsets
+  void set_entity_group_offsets(int dim,
+                                const std::vector<std::int32_t>& offsets);
+
+  /// @brief Get the offsets for each group of entities of a particular
+  /// dimension. The topology may consist of more than one cell type
+  /// or facet type. In that case, the cells of the same types are
+  /// grouped together in blocks, firstly for regular cells, then
+  /// repeated for ghost cells. For example, a mesh with two
+  /// triangles, three quads and no ghosts would have offsets: 0, 2,
+  /// 5, 5, 5. A mesh with twenty tetrahedra and two ghost tetrahedra
+  /// has offsets: 0, 20, 22.
+  /// @param dim Dimension of the entities
+  /// @return The offsets
+  const std::vector<std::int32_t>& entity_group_offsets(int dim) const;
+
   /// @todo Merge with set_connectivity
   ///
   /// Set the IndexMap for dimension dim
@@ -112,8 +132,8 @@ class Topology
   const std::vector<std::uint8_t>& get_facet_permutations() const;
 
   /// Cell type
-  /// @return Cell type that the topology is for
-  CellType cell_type() const noexcept;
+  /// @return Cell types that the topology is for
+  std::vector<CellType> cell_types() const noexcept;
 
   /// @brief Create entities of given topological dimension.
   /// @param[in] dim Topological dimension
@@ -144,8 +164,11 @@ class Topology
   // MPI communicator
   dolfinx::MPI::Comm _comm;
 
-  // Cell type
-  CellType _cell_type;
+  // Cell types
+  std::vector<CellType> _cell_types;
+
+  // Entity group offsets
+  std::array<std::vector<std::int32_t>, 4> _entity_group_offsets;
 
   // Parallel layout of entities for each dimension
   std::array<std::shared_ptr<const common::IndexMap>, 4> _index_map;
@@ -180,7 +203,9 @@ class Topology
 /// with each cell
 /// @param[in] ghost_owners The owning rank of each ghost cell (ghost
 /// cells are always at the end of the list of `cells`)
-/// @param[in] cell_type The cell shape
+/// @param[in] cell_type A vector with cell shapes
+/// @param[in] cell_group_offsets vector with each group offset, including
+/// ghosts.
 /// @param[in] boundary_vertices List of vertices on the exterior of the
 /// local mesh which may be shared with other processes.
 /// @return A distributed mesh topology
@@ -188,7 +213,8 @@ Topology create_topology(MPI_Comm comm,
                          const graph::AdjacencyList<std::int64_t>& cells,
                          std::span<const std::int64_t> original_cell_index,
                          std::span<const int> ghost_owners,
-                         const CellType& cell_type,
+                         const std::vector<CellType>& cell_type,
+                         const std::vector<std::int32_t>& cell_group_offsets,
                          std::span<const std::int64_t> boundary_vertices);
 
 /// @brief Create a topology for a subset of entities of a given
diff --git a/cpp/dolfinx/mesh/generation.h b/cpp/dolfinx/mesh/generation.h
index 446ab85de61..6090b10b37e 100644
--- a/cpp/dolfinx/mesh/generation.h
+++ b/cpp/dolfinx/mesh/generation.h
@@ -120,7 +120,7 @@ Mesh<T> create_interval(MPI_Comm comm, std::size_t nx, std::array<double, 2> x,
   {
     return create_mesh(
         comm, graph::regular_adjacency_list(std::vector<std::int64_t>(), 2),
-        element, std::vector<T>(), {0, 1}, partitioner);
+        {element}, std::vector<T>(), {0, 1}, partitioner);
   }
 
   const T a = x[0];
@@ -154,7 +154,7 @@ Mesh<T> create_interval(MPI_Comm comm, std::size_t nx, std::array<double, 2> x,
       cells[2 * ix + j] = ix + j;
 
   return create_mesh(comm, graph::regular_adjacency_list(std::move(cells), 2),
-                     element, geom, {geom.size(), 1}, partitioner);
+                     {element}, geom, {geom.size(), 1}, partitioner);
 }
 
 /// Create a uniform mesh::Mesh over the rectangle spanned by the two
@@ -328,7 +328,7 @@ Mesh<T> build_tet(MPI_Comm comm, const std::array<std::array<double, 3>, 2>& p,
 
   fem::CoordinateElement element(CellType::tetrahedron, 1);
   return create_mesh(comm, graph::regular_adjacency_list(std::move(cells), 4),
-                     element, geom, {geom.size() / 3, 3}, partitioner);
+                     {element}, geom, {geom.size() / 3, 3}, partitioner);
 }
 
 template <typename T>
@@ -372,7 +372,7 @@ mesh::Mesh<T> build_hex(MPI_Comm comm,
 
   fem::CoordinateElement element(CellType::hexahedron, 1);
   return create_mesh(comm, graph::regular_adjacency_list(std::move(cells), 8),
-                     element, geom, {geom.size() / 3, 3}, partitioner);
+                     {element}, geom, {geom.size() / 3, 3}, partitioner);
 }
 
 template <typename T>
@@ -420,7 +420,7 @@ Mesh<T> build_prism(MPI_Comm comm,
 
   fem::CoordinateElement element(CellType::prism, 1);
   return create_mesh(comm, graph::regular_adjacency_list(std::move(cells), 6),
-                     element, geom, {geom.size() / 3, 3}, partitioner);
+                     {element}, geom, {geom.size() / 3, 3}, partitioner);
 }
 
 template <typename T>
@@ -436,7 +436,7 @@ Mesh<T> build_tri(MPI_Comm comm, const std::array<std::array<double, 2>, 2>& p,
   {
     return create_mesh(
         comm, graph::regular_adjacency_list(std::vector<std::int64_t>(), 3),
-        element, std::vector<T>(), {0, 2}, partitioner);
+        {element}, std::vector<T>(), {0, 2}, partitioner);
   }
 
   const std::array<double, 2> p0 = p[0];
@@ -606,7 +606,7 @@ Mesh<T> build_tri(MPI_Comm comm, const std::array<std::array<double, 2>, 2>& p,
   }
 
   return create_mesh(comm, graph::regular_adjacency_list(std::move(cells), 3),
-                     element, geom, {geom.size() / 2, 2}, partitioner);
+                     {element}, geom, {geom.size() / 2, 2}, partitioner);
 }
 
 template <typename T>
@@ -621,7 +621,7 @@ Mesh<T> build_quad(MPI_Comm comm, const std::array<std::array<double, 2>, 2> p,
   {
     return create_mesh(
         comm, graph::regular_adjacency_list(std::vector<std::int64_t>(), 4),
-        element, std::vector<T>(), {0, 2}, partitioner);
+        {element}, std::vector<T>(), {0, 2}, partitioner);
   }
 
   const std::size_t nx = n[0];
@@ -664,7 +664,7 @@ Mesh<T> build_quad(MPI_Comm comm, const std::array<std::array<double, 2>, 2> p,
   }
 
   return create_mesh(comm, graph::regular_adjacency_list(std::move(cells), 4),
-                     element, geom, {geom.size() / 2, 2}, partitioner);
+                     {element}, geom, {geom.size() / 2, 2}, partitioner);
 }
 
 } // namespace impl
diff --git a/cpp/dolfinx/mesh/permutationcomputation.cpp b/cpp/dolfinx/mesh/permutationcomputation.cpp
index bc021f8290c..817a6f6584f 100644
--- a/cpp/dolfinx/mesh/permutationcomputation.cpp
+++ b/cpp/dolfinx/mesh/permutationcomputation.cpp
@@ -230,8 +230,13 @@ template <int BITSETSIZE>
 std::vector<std::bitset<BITSETSIZE>>
 compute_edge_reflections(const mesh::Topology& topology)
 {
+
+  auto cell_types = topology.cell_types();
+  if (cell_types.size() > 1)
+    throw std::runtime_error("cell type permutation");
+
   const int tdim = topology.dim();
-  const mesh::CellType cell_type = topology.cell_type();
+  const mesh::CellType cell_type = cell_types.back();
   const int edges_per_cell = cell_num_entities(cell_type, 1);
 
   const std::int32_t num_cells = topology.connectivity(tdim, 0)->num_nodes();
@@ -295,7 +300,12 @@ compute_face_permutations(const mesh::Topology& topology)
 
   auto im = topology.index_map(0);
   assert(im);
-  const mesh::CellType cell_type = topology.cell_type();
+
+  auto cell_types = topology.cell_types();
+  if (cell_types.size() > 1)
+    throw std::runtime_error("cell type permutation");
+
+  const mesh::CellType cell_type = cell_types.back();
   const int faces_per_cell = cell_num_entities(cell_type, 2);
   if (cell_type == mesh::CellType::triangle
       or cell_type == mesh::CellType::tetrahedron)
@@ -317,7 +327,11 @@ std::pair<std::vector<std::uint8_t>, std::vector<std::uint32_t>>
 mesh::compute_entity_permutations(const mesh::Topology& topology)
 {
   const int tdim = topology.dim();
-  const CellType cell_type = topology.cell_type();
+  auto cell_types = topology.cell_types();
+  if (cell_types.size() > 1)
+    throw std::runtime_error("cell type permutation");
+  const CellType cell_type = cell_types.back();
+
   const std::int32_t num_cells = topology.connectivity(tdim, 0)->num_nodes();
   const int facets_per_cell = cell_num_entities(cell_type, tdim - 1);
 
diff --git a/cpp/dolfinx/mesh/topologycomputation.cpp b/cpp/dolfinx/mesh/topologycomputation.cpp
index 9f553f70512..56dcbf2b0bd 100644
--- a/cpp/dolfinx/mesh/topologycomputation.cpp
+++ b/cpp/dolfinx/mesh/topologycomputation.cpp
@@ -541,7 +541,7 @@ compute_entities_by_key_matching(
     for (std::size_t i = 0; i < entity_list_shape0; ++i)
     {
       auto it = std::next(entity_list_sorted.begin(), i * entity_list_shape1);
-      std::sort(it, std::next(it, entity_list_shape1), std::greater<>());
+      std::sort(it, std::next(it, entity_list_shape1), std::less<>());
     }
 
     // Sort the list and label uniquely
@@ -591,7 +591,6 @@ compute_entities_by_key_matching(
   std::vector<int> size_ev(entity_count);
   for (std::size_t i = 0; i < entity_list_shape0; ++i)
   {
-    // if (entity_list(i, max_vertices_per_entity - 1) == -1)
     if (entity_list[i * entity_list_shape1 + entity_list_shape1 - 1] == -1)
       size_ev[local_index[i]] = max_vertices_per_entity - 1;
     else
@@ -762,8 +761,9 @@ mesh::compute_entities(MPI_Comm comm, const Topology& topology, int dim)
   assert(vertex_map);
   auto cell_map = topology.index_map(tdim);
   assert(cell_map);
+
   auto [d0, d1, im, interprocess_facets] = compute_entities_by_key_matching(
-      comm, *cells, *vertex_map, *cell_map, topology.cell_type(), dim);
+      comm, *cells, *vertex_map, *cell_map, topology.cell_types().back(), dim);
 
   return {std::make_shared<graph::AdjacencyList<std::int32_t>>(std::move(d0)),
           std::make_shared<graph::AdjacencyList<std::int32_t>>(std::move(d1)),
diff --git a/cpp/dolfinx/mesh/utils.h b/cpp/dolfinx/mesh/utils.h
index 6931d48d5b6..987b4f9bc01 100644
--- a/cpp/dolfinx/mesh/utils.h
+++ b/cpp/dolfinx/mesh/utils.h
@@ -35,6 +35,38 @@ enum class GhostMode : int
   shared_vertex
 };
 
+namespace
+{
+/// Re-order an adjacency list
+template <typename T>
+graph::AdjacencyList<T> reorder_list(const graph::AdjacencyList<T>& list,
+                                     std::span<const std::int32_t> nodemap)
+{
+  // Copy existing data to keep ghost values (not reordered)
+  std::vector<T> data(list.array());
+  std::vector<std::int32_t> offsets(list.offsets().size());
+
+  // Compute new offsets (owned and ghost)
+  offsets[0] = 0;
+  for (std::size_t n = 0; n < nodemap.size(); ++n)
+    offsets[nodemap[n] + 1] = list.num_links(n);
+  for (std::size_t n = nodemap.size(); n < (std::size_t)list.num_nodes(); ++n)
+    offsets[n + 1] = list.num_links(n);
+  std::partial_sum(offsets.begin(), offsets.end(), offsets.begin());
+  graph::AdjacencyList<T> list_new(std::move(data), std::move(offsets));
+
+  for (std::size_t n = 0; n < nodemap.size(); ++n)
+  {
+    auto links_old = list.links(n);
+    auto links_new = list_new.links(nodemap[n]);
+    assert(links_old.size() == links_new.size());
+    std::copy(links_old.begin(), links_old.end(), links_new.begin());
+  }
+
+  return list_new;
+}
+} // namespace
+
 // See https://github.com/doxygen/doxygen/issues/9552
 /// Signature for the cell partitioning function. The function should
 /// compute the destination rank for cells currently on this rank.
@@ -248,18 +280,24 @@ std::vector<T> cell_normals(const Mesh<T>& mesh, int dim,
   if (entities.empty())
     return std::vector<T>();
 
-  if (mesh.topology().cell_type() == CellType::prism and dim == 2)
+  if (mesh.topology().cell_types().size() > 1)
+  {
+    throw std::runtime_error("Mixed topology not supported");
+  }
+
+  if (mesh.topology().cell_types()[0] == CellType::prism and dim == 2)
     throw std::runtime_error("More work needed for prism cell");
 
   const int gdim = mesh.geometry().dim();
-  const CellType type = cell_entity_type(mesh.topology().cell_type(), dim, 0);
+  const CellType type
+      = cell_entity_type(mesh.topology().cell_types()[0], dim, 0);
 
   // Find geometry nodes for topology entities
   std::span<const T> x = mesh.geometry().x();
 
   // Orient cells if they are tetrahedron
   bool orient = false;
-  if (mesh.topology().cell_type() == CellType::tetrahedron)
+  if (mesh.topology().cell_types()[0] == CellType::tetrahedron)
     orient = true;
 
   std::vector<std::int32_t> geometry_entities
@@ -596,7 +634,11 @@ std::vector<std::int32_t>
 entities_to_geometry(const Mesh<T>& mesh, int dim,
                      std::span<const std::int32_t> entities, bool orient)
 {
-  CellType cell_type = mesh.topology().cell_type();
+  if (mesh.topology().cell_types().size() > 1)
+  {
+    throw std::runtime_error("Mixed topology not supported");
+  }
+  CellType cell_type = mesh.topology().cell_types()[0];
   if (cell_type == CellType::prism and dim == 2)
     throw std::runtime_error("More work needed for prism cells");
   if (orient and (cell_type != CellType::tetrahedron or dim != 2))
@@ -704,14 +746,14 @@ compute_incident_entities(const Topology& topology,
 template <typename U>
 Mesh<typename std::remove_reference_t<typename U::value_type>>
 create_mesh(MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& cells,
-            const fem::CoordinateElement& element, const U& x,
+            const std::vector<fem::CoordinateElement>& elements, const U& x,
             std::array<std::size_t, 2> xshape,
             const CellPartitionFunction& cell_partitioner)
 {
-  const fem::ElementDofLayout dof_layout = element.create_dof_layout();
+  const fem::ElementDofLayout dof_layout = elements[0].create_dof_layout();
 
   // Function top build geometry. Used to scope memory operations.
-  auto build_topology = [](auto comm, auto& element, auto& dof_layout,
+  auto build_topology = [](auto comm, auto& elements, auto& dof_layout,
                            auto& cells, auto& cell_partitioner)
   {
     // -- Partition topology
@@ -730,10 +772,10 @@ create_mesh(MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& cells,
     // Compute the destination rank for cells on this process via graph
     // partitioning.
     const int size = dolfinx::MPI::size(comm);
-    const int tdim = cell_dim(element.cell_shape());
-    const graph::AdjacencyList dest = cell_partitioner(
+    const int tdim = cell_dim(elements[0].cell_shape());
+    const graph::AdjacencyList<std::int32_t> dest = cell_partitioner(
         comm, size, tdim,
-        extract_topology(element.cell_shape(), dof_layout, cells));
+        extract_topology(elements[0].cell_shape(), dof_layout, cells));
 
     // -- Distribute cells (topology, includes higher-order 'nodes')
 
@@ -748,8 +790,9 @@ create_mesh(MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& cells,
 
     // Extract cell 'topology', i.e. extract the vertices for each cell
     // and discard any 'higher-order' nodes
-    graph::AdjacencyList cells_extracted
-        = extract_topology(element.cell_shape(), dof_layout, cell_nodes);
+
+    graph::AdjacencyList<std::int64_t> cells_extracted
+        = extract_topology(elements[0].cell_shape(), dof_layout, cell_nodes);
 
     // -- Re-order cells
 
@@ -768,38 +811,6 @@ create_mesh(MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& cells,
 
     const std::vector<int> remap = graph::reorder_gps(graph);
 
-    /// Re-order an adjacency list
-    auto reorder_list
-        = [](const auto& list, std::span<const std::int32_t> nodemap)
-    {
-      using X =
-          typename std::remove_reference_t<decltype(list.array())>::value_type;
-
-      // Copy existing data to keep ghost values (not reordered)
-      std::vector<X> data(list.array());
-      std::vector<std::int32_t> offsets(list.offsets().size());
-
-      // Compute new offsets (owned and ghost)
-      offsets[0] = 0;
-      for (std::size_t n = 0; n < nodemap.size(); ++n)
-        offsets[nodemap[n] + 1] = list.num_links(n);
-      for (std::size_t n = nodemap.size(); n < (std::size_t)list.num_nodes();
-           ++n)
-        offsets[n + 1] = list.num_links(n);
-      std::partial_sum(offsets.begin(), offsets.end(), offsets.begin());
-      graph::AdjacencyList<X> list_new(std::move(data), std::move(offsets));
-
-      for (std::size_t n = 0; n < nodemap.size(); ++n)
-      {
-        auto links_old = list.links(n);
-        auto links_new = list_new.links(nodemap[n]);
-        assert(links_old.size() == links_new.size());
-        std::copy(links_old.begin(), links_old.end(), links_new.begin());
-      }
-
-      return list_new;
-    };
-
     // Create re-ordered cell lists (leaves ghosts unchanged)
     std::vector<std::int64_t> original_cell_index(original_cell_index0.size());
     for (std::size_t i = 0; i < remap.size(); ++i)
@@ -826,14 +837,20 @@ create_mesh(MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& cells,
     // Create cells and vertices with the ghosting requested. Input
     // topology includes cells shared via facet, but ghosts will be
     // removed later if not required by ghost_mode.
+
+    std::vector<std::int32_t> cell_group_offsets
+        = {0, std::int32_t(cells_extracted.num_nodes() - ghost_owners.size()),
+           cells_extracted.num_nodes()};
+
+    std::vector<mesh::CellType> cell_type = {elements[0].cell_shape()};
     return std::pair{create_topology(comm, cells_extracted, original_cell_index,
-                                     ghost_owners, element.cell_shape(),
-                                     boundary_vertices),
+                                     ghost_owners, cell_type,
+                                     cell_group_offsets, boundary_vertices),
                      std::move(cell_nodes)};
   };
 
   auto [topology, cell_nodes]
-      = build_topology(comm, element, dof_layout, cells, cell_partitioner);
+      = build_topology(comm, elements, dof_layout, cells, cell_partitioner);
 
   // Create connectivity required to compute the Geometry (extra
   // connectivities for higher-order geometries)
@@ -844,11 +861,12 @@ create_mesh(MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& cells,
       topology.create_entities(e);
   }
 
-  if (element.needs_dof_permutations())
+  if (elements[0].needs_dof_permutations())
     topology.create_entity_permutations();
 
   Geometry geometry
-      = create_geometry(comm, topology, element, cell_nodes, x, xshape[1]);
+      = create_geometry(comm, topology, elements, cell_nodes, x, xshape[1]);
+
   return Mesh<typename U::value_type>(comm, std::move(topology),
                                       std::move(geometry));
 }
@@ -866,7 +884,7 @@ create_mesh(MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& cells,
 /// indices). For lowest order cells this will be just the cell
 /// vertices. For higher-order cells, other cells 'nodes' will be
 /// included.
-/// @param[in] element The coordinate element that describes the
+/// @param[in] elements The coordinate elements that describe the
 /// geometric mapping for cells
 /// @param[in] x The coordinates of mesh nodes
 /// @param[in] xshape The shape of `x`. It should be `(num_points, gdim)`.
@@ -875,10 +893,10 @@ create_mesh(MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& cells,
 template <typename U>
 Mesh<typename std::remove_reference_t<typename U::value_type>>
 create_mesh(MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& cells,
-            const fem::CoordinateElement& element, const U& x,
+            const std::vector<fem::CoordinateElement>& elements, const U& x,
             std::array<std::size_t, 2> xshape, GhostMode ghost_mode)
 {
-  return create_mesh(comm, cells, element, x, xshape,
+  return create_mesh(comm, cells, elements, x, xshape,
                      create_cell_partitioner(ghost_mode));
 }
 
diff --git a/cpp/dolfinx/refinement/plaza.h b/cpp/dolfinx/refinement/plaza.h
index 6762c088c58..09b3b5f363f 100644
--- a/cpp/dolfinx/refinement/plaza.h
+++ b/cpp/dolfinx/refinement/plaza.h
@@ -32,7 +32,7 @@ enum class Option : int
 {
   none = 0, /*!< No extra data */
   parent_cell
-  = 1,      /*!< Compute list with the parent cell index for each new cell  */
+  = 1, /*!< Compute list with the parent cell index for each new cell  */
   parent_facet
   = 2, /*!< Compute list of the cell-local facet indices in the parent cell of
           each facet in each new cell (or -1 if no match) */
@@ -423,12 +423,17 @@ template <typename T>
 std::tuple<mesh::Mesh<T>, std::vector<std::int32_t>, std::vector<std::int8_t>>
 refine(const mesh::Mesh<T>& mesh, bool redistribute, Option option)
 {
+  if (mesh.geometry().cmaps().size() > 1)
+  {
+    throw std::runtime_error("Mixed topology not supported");
+  }
+
   auto [cell_adj, new_coords, xshape, parent_cell, parent_facet]
       = compute_refinement_data(mesh, option);
 
   if (dolfinx::MPI::size(mesh.comm()) == 1)
   {
-    return {mesh::create_mesh(mesh.comm(), cell_adj, mesh.geometry().cmap(),
+    return {mesh::create_mesh(mesh.comm(), cell_adj, mesh.geometry().cmaps(),
                               new_coords, xshape, mesh::GhostMode::none),
             std::move(parent_cell), std::move(parent_facet)};
   }
@@ -470,12 +475,17 @@ refine(const mesh::Mesh<T>& mesh, std::span<const std::int32_t> edges,
        bool redistribute, Option option)
 {
 
+  if (mesh.geometry().cmaps().size() > 1)
+  {
+    throw std::runtime_error("Mixed topology not supported");
+  }
+
   auto [cell_adj, new_vertex_coords, xshape, parent_cell, parent_facet]
       = compute_refinement_data(mesh, edges, option);
 
   if (dolfinx::MPI::size(mesh.comm()) == 1)
   {
-    return {mesh::create_mesh(mesh.comm(), cell_adj, mesh.geometry().cmap(),
+    return {mesh::create_mesh(mesh.comm(), cell_adj, mesh.geometry().cmaps(),
                               new_vertex_coords, xshape, mesh::GhostMode::none),
             std::move(parent_cell), std::move(parent_facet)};
   }
@@ -517,8 +527,13 @@ compute_refinement_data(const mesh::Mesh<T>& mesh, Option option)
 {
   common::Timer t0("PLAZA: refine");
 
-  if (mesh.topology().cell_type() != mesh::CellType::triangle
-      and mesh.topology().cell_type() != mesh::CellType::tetrahedron)
+  if (mesh.topology().cell_types().size() > 1)
+  {
+    throw std::runtime_error("Mixed topology not supported");
+  }
+
+  if (mesh.topology().cell_types()[0] != mesh::CellType::triangle
+      and mesh.topology().cell_types()[0] != mesh::CellType::tetrahedron)
   {
     throw std::runtime_error("Cell type not supported");
   }
@@ -582,8 +597,13 @@ compute_refinement_data(const mesh::Mesh<T>& mesh,
 {
   common::Timer t0("PLAZA: refine");
 
-  if (mesh.topology().cell_type() != mesh::CellType::triangle
-      and mesh.topology().cell_type() != mesh::CellType::tetrahedron)
+  if (mesh.topology().cell_types().size() > 1)
+  {
+    throw std::runtime_error("Mixed topology not supported");
+  }
+
+  if (mesh.topology().cell_types()[0] != mesh::CellType::triangle
+      and mesh.topology().cell_types()[0] != mesh::CellType::tetrahedron)
   {
     throw std::runtime_error("Cell type not supported");
   }
diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index 900e0d5e2b1..53bd80021cb 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -24,8 +24,13 @@ namespace dolfinx::refinement
 template <typename T>
 mesh::Mesh<T> refine(const mesh::Mesh<T>& mesh, bool redistribute = true)
 {
-  if (mesh.topology().cell_type() != mesh::CellType::triangle
-      and mesh.topology().cell_type() != mesh::CellType::tetrahedron)
+  if (mesh.topology().cell_types().size() > 1)
+  {
+    throw std::runtime_error("Mixed topology not supported");
+  }
+
+  if (mesh.topology().cell_types()[0] != mesh::CellType::triangle
+      and mesh.topology().cell_types()[0] != mesh::CellType::tetrahedron)
   {
     throw std::runtime_error("Refinement only defined for simplices");
   }
@@ -59,8 +64,13 @@ mesh::Mesh<T> refine(const mesh::Mesh<T>& mesh,
                      std::span<const std::int32_t> edges,
                      bool redistribute = true)
 {
-  if (mesh.topology().cell_type() != mesh::CellType::triangle
-      and mesh.topology().cell_type() != mesh::CellType::tetrahedron)
+  if (mesh.topology().cell_types().size() > 1)
+  {
+    throw std::runtime_error("Mixed topology not supported");
+  }
+
+  if (mesh.topology().cell_types()[0] != mesh::CellType::triangle
+      and mesh.topology().cell_types()[0] != mesh::CellType::tetrahedron)
   {
     throw std::runtime_error("Refinement only defined for simplices");
   }
diff --git a/cpp/dolfinx/refinement/utils.h b/cpp/dolfinx/refinement/utils.h
index ed4b1bc82f9..a36a7db0c2e 100644
--- a/cpp/dolfinx/refinement/utils.h
+++ b/cpp/dolfinx/refinement/utils.h
@@ -67,7 +67,7 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> create_new_geometry(
   auto map_c = mesh.topology().index_map(tdim);
 
   assert(map_c);
-  auto dof_layout = mesh.geometry().cmap().create_dof_layout();
+  auto dof_layout = mesh.geometry().cmaps()[0].create_dof_layout();
   auto entity_dofs_all = dof_layout.entity_dofs_all();
   for (int c = 0; c < map_c->size_local() + map_c->num_ghosts(); ++c)
   {
@@ -279,7 +279,7 @@ mesh::Mesh<T> partition(const mesh::Mesh<T>& old_mesh,
   if (redistribute)
   {
     return mesh::create_mesh(old_mesh.comm(), cell_topology,
-                             old_mesh.geometry().cmap(), new_coords, xshape,
+                             old_mesh.geometry().cmaps(), new_coords, xshape,
                              ghost_mode);
   }
   else
@@ -298,7 +298,7 @@ mesh::Mesh<T> partition(const mesh::Mesh<T>& old_mesh,
     };
 
     return mesh::create_mesh(old_mesh.comm(), cell_topology,
-                             old_mesh.geometry().cmap(), new_coords, xshape,
+                             old_mesh.geometry().cmaps(), new_coords, xshape,
                              partitioner);
   }
 }
diff --git a/cpp/test/mesh/distributed_mesh.cpp b/cpp/test/mesh/distributed_mesh.cpp
index 58e41232b56..d3e3ef84dbf 100644
--- a/cpp/test/mesh/distributed_mesh.cpp
+++ b/cpp/test/mesh/distributed_mesh.cpp
@@ -113,12 +113,17 @@ void test_distributed_mesh(mesh::CellPartitionFunction partitioner)
                     [](std::int64_t i) { return (i != -1); });
   external_vertices.erase(external_vertices.begin(), it);
 
+  std::vector<int> cell_group_offsets
+      = {0, std::int32_t(cell_nodes.num_nodes() - ghost_owners.size()),
+         cell_nodes.num_nodes()};
+
+  std::vector<mesh::CellType> cell_types = {cmap.cell_shape()};
   mesh::Topology topology = mesh::create_topology(
-      mpi_comm, cell_nodes, original_cell_index, ghost_owners,
-      cmap.cell_shape(), external_vertices);
+      mpi_comm, cell_nodes, original_cell_index, ghost_owners, cell_types,
+      cell_group_offsets, external_vertices);
   int tdim = topology.dim();
 
-  mesh::Geometry geometry = mesh::create_geometry(mpi_comm, topology, cmap,
+  mesh::Geometry geometry = mesh::create_geometry(mpi_comm, topology, {cmap},
                                                   cell_nodes, x, xshape[1]);
 
   auto mesh = std::make_shared<mesh::Mesh<double>>(
diff --git a/python/dolfinx/io/gmshio.py b/python/dolfinx/io/gmshio.py
index fb084ef0d04..adcbb52f3a7 100644
--- a/python/dolfinx/io/gmshio.py
+++ b/python/dolfinx/io/gmshio.py
@@ -261,8 +261,8 @@ def model_to_mesh(model: gmsh.model, comm: _MPI.Comm, rank: int, gdim: int = 3,
         if has_facet_data:
             # Permute facets from MSH to DOLFINx ordering
             # FIXME: This does not work for prism meshes
-            if topology.cell_type == CellType.prism or topology.cell_type == CellType.pyramid:
-                raise RuntimeError(f"Unsupported cell type {topology.cell_type}")
+            if topology.cell_types[0] == CellType.prism or topology.cell_types[0] == CellType.pyramid:
+                raise RuntimeError(f"Unsupported cell type {topology.cell_types[0]}")
 
             facet_type = _cpp.mesh.cell_entity_type(_cpp.mesh.to_type(str(ufl_domain.ufl_cell())), topology.dim - 1, 0)
             gmsh_facet_perm = cell_perm_array(facet_type, num_facet_nodes)
diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index 314e6cd7848..77a7405e203 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -265,9 +265,10 @@ def create_mesh(comm: _MPI.Comm, cells: typing.Union[np.ndarray, _cpp.graph.Adja
 
 def create_submesh(msh, dim, entities):
     submsh, entity_map, vertex_map, geom_map = _cpp.mesh.create_submesh(msh._cpp_object, dim, entities)
+    assert len(submsh.geometry.cmaps) == 1
     submsh_ufl_cell = ufl.Cell(submsh.topology.cell_name(), geometric_dimension=submsh.geometry.dim)
     submsh_domain = ufl.Mesh(basix.ufl.element(
-        "Lagrange", submsh_ufl_cell.cellname(), submsh.geometry.cmap.degree, submsh.geometry.cmap.variant,
+        "Lagrange", submsh_ufl_cell.cellname(), submsh.geometry.cmaps[0].degree, submsh.geometry.cmaps[0].variant,
         shape=(submsh.geometry.dim, ), gdim=submsh.geometry.dim))
     return (Mesh(submsh, submsh_domain), entity_map, vertex_map, geom_map)
 
diff --git a/python/dolfinx/plot.py b/python/dolfinx/plot.py
index 7df396a9065..357cd804b84 100644
--- a/python/dolfinx/plot.py
+++ b/python/dolfinx/plot.py
@@ -41,7 +41,9 @@ def create_vtk_mesh(msh: mesh.Mesh, dim: typing.Optional[int] = None, entities=N
 
     tdim = msh.topology.dim
 
-    cell_type = _cpp.mesh.cell_entity_type(msh.topology.cell_type, dim, 0)
+    if len(msh.topology.cell_types) != 1:
+        raise RuntimeError("Multiple cell types")
+    cell_type = _cpp.mesh.cell_entity_type(msh.topology.cell_types[0], dim, 0)
     degree = msh.geometry.cmap.degree
     if cell_type == mesh.CellType.prism:
         raise RuntimeError("Plotting of prism meshes not supported")
@@ -99,7 +101,9 @@ def _(V: fem.FunctionSpace, entities=None):
 
     dofmap = V.dofmap
     num_dofs_per_cell = V.dofmap.dof_layout.num_dofs
-    cell_type = msh.topology.cell_type
+    if len(msh.topology.cell_types) != 1:
+        raise RuntimeError("Multiple cell types")
+    cell_type = msh.topology.cell_types[0]
     perm = np.argsort(_cpp.io.perm_vtk(cell_type, num_dofs_per_cell))
 
     vtk_type = _first_order_vtk[cell_type] if degree == 1 else _cpp.io.get_vtk_cell_type(cell_type, tdim)
diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp
index f4a9a5d1bdf..2997f5da886 100644
--- a/python/dolfinx/wrappers/fem.cpp
+++ b/python/dolfinx/wrappers/fem.cpp
@@ -609,8 +609,10 @@ void fem(py::module& m)
       {
         ufcx_dofmap* p = reinterpret_cast<ufcx_dofmap*>(dofmap);
         assert(p);
+
         dolfinx::fem::ElementDofLayout layout
-            = dolfinx::fem::create_element_dof_layout(*p, topology.cell_type());
+            = dolfinx::fem::create_element_dof_layout(
+                *p, topology.cell_types().back());
         return dolfinx::fem::create_dofmap(comm.get(), layout, topology,
                                            nullptr, element);
       },
@@ -623,7 +625,7 @@ void fem(py::module& m)
          const dolfinx::fem::ElementDofLayout& layout)
       {
         auto [map, bs, dofmap] = dolfinx::fem::build_dofmap_data(
-            comm.get(), topology, layout,
+            comm.get(), topology, {layout},
             [](const dolfinx::graph::AdjacencyList<std::int32_t>& g)
             { return dolfinx::graph::reorder_gps(g); });
         return std::tuple(std::move(map), bs, std::move(dofmap));
diff --git a/python/dolfinx/wrappers/io.cpp b/python/dolfinx/wrappers/io.cpp
index c2c293f9a28..d294729fac6 100644
--- a/python/dolfinx/wrappers/io.cpp
+++ b/python/dolfinx/wrappers/io.cpp
@@ -74,7 +74,7 @@ void io(py::module& m)
                 std::span(values.data(), values.size()));
 
         std::size_t num_vert_per_entity = dolfinx::mesh::cell_num_entities(
-            dolfinx::mesh::cell_entity_type(mesh.topology().cell_type(),
+            dolfinx::mesh::cell_entity_type(mesh.topology().cell_types().back(),
                                             entity_dim, 0),
             0);
         std::array shape_e
diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp
index 33c3eb0f8a9..315879a5076 100644
--- a/python/dolfinx/wrappers/mesh.cpp
+++ b/python/dolfinx/wrappers/mesh.cpp
@@ -226,7 +226,7 @@ void mesh(py::module& m)
         std::size_t shape1 = x.ndim() == 1 ? 1 : x.shape()[1];
         std::vector shape{std::size_t(x.shape(0)), shape1};
         return dolfinx::mesh::create_mesh(
-            comm.get(), cells, element, std::span(x.data(), x.size()),
+            comm.get(), cells, {element}, std::span(x.data(), x.size()),
             {static_cast<std::size_t>(x.shape(0)),
              static_cast<std::size_t>(x.shape(1))},
             partitioner_wrapper);
@@ -267,8 +267,8 @@ void mesh(py::module& m)
           },
           "Return coordinates of all geometry points. Each row is the "
           "coordinate of a point.")
-      .def_property_readonly("cmap", &dolfinx::mesh::Geometry<double>::cmap,
-                             "The coordinate map")
+      .def_property_readonly("cmaps", &dolfinx::mesh::Geometry<double>::cmaps,
+                             "The coordinate maps")
       .def_property_readonly(
           "input_global_indices",
           &dolfinx::mesh::Geometry<double>::input_global_indices);
@@ -287,9 +287,11 @@ void mesh(py::module& m)
   py::class_<dolfinx::mesh::Topology, std::shared_ptr<dolfinx::mesh::Topology>>(
       m, "Topology", py::dynamic_attr(), "Topology object")
       .def(py::init([](const MPICommWrapper comm,
-                       const dolfinx::mesh::CellType cell_type)
+                       const std::vector<dolfinx::mesh::CellType> cell_type)
                     { return dolfinx::mesh::Topology(comm.get(), cell_type); }),
            py::arg("comm"), py::arg("cell_type"))
+      .def("entity_group_offsets",
+           &dolfinx::mesh::Topology::entity_group_offsets)
       .def("set_connectivity", &dolfinx::mesh::Topology::set_connectivity,
            py::arg("c"), py::arg("d0"), py::arg("d1"))
       .def("set_index_map", &dolfinx::mesh::Topology::set_index_map,
@@ -330,9 +332,14 @@ void mesh(py::module& m)
                                        py::const_),
            py::arg("d0"), py::arg("d1"))
       .def("index_map", &dolfinx::mesh::Topology::index_map, py::arg("dim"))
-      .def_property_readonly("cell_type", &dolfinx::mesh::Topology::cell_type)
-      .def("cell_name", [](const dolfinx::mesh::Topology& self)
-           { return dolfinx::mesh::to_string(self.cell_type()); })
+      .def_property_readonly("cell_types", &dolfinx::mesh::Topology::cell_types)
+      .def("cell_name",
+           [](const dolfinx::mesh::Topology& self)
+           {
+             if (self.cell_types().size() > 1)
+               throw std::runtime_error("Multiple cell types not supported");
+             return dolfinx::mesh::to_string(self.cell_types()[0]);
+           })
       .def("interprocess_facets", &dolfinx::mesh::Topology::interprocess_facets)
       .def_property_readonly("comm", [](dolfinx::mesh::Mesh<double>& self)
                              { return MPICommWrapper(self.comm()); });
@@ -438,7 +445,11 @@ void mesh(py::module& m)
       {
         std::vector<std::int32_t> idx = dolfinx::mesh::entities_to_geometry(
             mesh, dim, std::span(entities.data(), entities.size()), orient);
-        dolfinx::mesh::CellType cell_type = mesh.topology().cell_type();
+
+        if (mesh.topology().cell_types().size() > 1)
+          throw std::runtime_error("Multiple cell type not supported.");
+
+        dolfinx::mesh::CellType cell_type = mesh.topology().cell_types()[0];
         std::size_t num_vertices = dolfinx::mesh::num_cell_vertices(
             cell_entity_type(cell_type, dim, 0));
         std::array<std::size_t, 2> shape
diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py
index b1fead5ee38..6b2176c2cb5 100644
--- a/python/test/unit/fem/test_assembler.py
+++ b/python/test/unit/fem/test_assembler.py
@@ -933,7 +933,7 @@ def partitioner(comm, nparts, local_graph, num_ghost_nodes):
         cells = graph.create_adjacencylist(cells)
         x = np.array([[0., 0.], [1., 0.], [1., 1.], [0., 1.]])
     else:
-        # No cells onm other ranks
+        # No cells on other ranks
         cells = graph.create_adjacencylist(np.empty((0, 3), dtype=np.int64))
         x = np.empty((0, 2), dtype=np.float64)
 
diff --git a/python/test/unit/fem/test_dof_permuting.py b/python/test/unit/fem/test_dof_permuting.py
index c2265b3ced8..16b47f55505 100644
--- a/python/test/unit/fem/test_dof_permuting.py
+++ b/python/test/unit/fem/test_dof_permuting.py
@@ -120,7 +120,8 @@ def test_dof_positions(cell_type, space_type):
     # for each global dof number
     coord_dofs = mesh.geometry.dofmap
     x_g = mesh.geometry.x
-    cmap = mesh.geometry.cmap
+    assert len(mesh.geometry.cmaps) == 1
+    cmap = mesh.geometry.cmaps[0]
     tdim = mesh.topology.dim
 
     V = FunctionSpace(mesh, space_type)
diff --git a/python/test/unit/fem/test_dofmap.py b/python/test/unit/fem/test_dofmap.py
index cf99e52d898..1b8f805b59f 100644
--- a/python/test/unit/fem/test_dofmap.py
+++ b/python/test/unit/fem/test_dofmap.py
@@ -253,7 +253,8 @@ def test_higher_order_coordinate_map(points, celltype, order):
     X = V.element.interpolation_points()
     coord_dofs = mesh.geometry.dofmap
     x_g = mesh.geometry.x
-    cmap = mesh.geometry.cmap
+    assert len(mesh.geometry.cmaps) == 1
+    cmap = mesh.geometry.cmaps[0]
 
     x_coord_new = np.zeros([len(points), mesh.geometry.dim])
 
@@ -300,7 +301,8 @@ def test_higher_order_tetra_coordinate_map(order):
     coord_dofs = mesh.geometry.dofmap
     x_g = mesh.geometry.x
 
-    cmap = mesh.geometry.cmap
+    assert len(mesh.geometry.cmaps) == 1
+    cmap = mesh.geometry.cmaps[0]
     x_coord_new = np.zeros([len(points), mesh.geometry.dim])
 
     i = 0
diff --git a/python/test/unit/fem/test_expression.py b/python/test/unit/fem/test_expression.py
index 680c5e182d6..fd49ecedfd5 100644
--- a/python/test/unit/fem/test_expression.py
+++ b/python/test/unit/fem/test_expression.py
@@ -339,11 +339,13 @@ def e_exact(x):
 
     Q_dofs_unrolled = bs * np.repeat(Q_dofs, bs).reshape(-1, bs) + np.arange(bs)
     Q_dofs_unrolled = Q_dofs_unrolled.reshape(-1, bs * quadrature_points.shape[0]).astype(Q_dofs.dtype)
+    assert len(mesh.geometry.cmaps) == 1
+
     with e_Q.vector.localForm() as local:
         e_exact_eval = np.zeros_like(local.array)
         for cell in range(num_cells):
             xg = x_g[coord_dofs.links(cell), :tdim]
-            x = mesh.geometry.cmap.push_forward(quadrature_points, xg)
+            x = mesh.geometry.cmaps[0].push_forward(quadrature_points, xg)
             e_exact_eval[Q_dofs_unrolled[cell]] = e_exact(x.T).T.flatten()
         assert np.allclose(local.array, e_exact_eval)
 
diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py
index 7f469b268e8..519b2a9dd4e 100644
--- a/python/test/unit/fem/test_interpolation.py
+++ b/python/test/unit/fem/test_interpolation.py
@@ -57,7 +57,8 @@ def random_point_in_reference(cell_type):
 
 
 def random_point_in_cell(mesh):
-    cell_type = mesh.topology.cell_type
+    assert len(mesh.topology.cell_types) == 1
+    cell_type = mesh.topology.cell_types[0]
     point = random_point_in_reference(cell_type)
 
     if cell_type == CellType.interval:
diff --git a/python/test/unit/io/test_xdmf_meshdata.py b/python/test/unit/io/test_xdmf_meshdata.py
index d6c95f38de2..227d5847fd8 100644
--- a/python/test/unit/io/test_xdmf_meshdata.py
+++ b/python/test/unit/io/test_xdmf_meshdata.py
@@ -47,7 +47,8 @@ def test_read_mesh_data(tempdir, tdim, n):
         cells = file.read_topology_data()
         x = file.read_geometry_data()
 
-    assert cell_shape == mesh.topology.cell_type
+    assert len(mesh.topology.cell_types) == 1
+    assert cell_shape == mesh.topology.cell_types[0]
     assert cell_degree == 1
     assert mesh.topology.index_map(tdim).size_global == mesh.comm.allreduce(cells.shape[0], op=MPI.SUM)
     assert mesh.geometry.index_map().size_global == mesh.comm.allreduce(x.shape[0], op=MPI.SUM)