From d2d5d2c9e7297eba4fb80742015043cac516f190 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Wed, 4 Jun 2025 11:50:22 +0100
Subject: [PATCH 01/28] wip

---
 cpp/dolfinx/refinement/CMakeLists.txt | 18 ++++++++++--------
 python/CMakeLists.txt                 |  2 +-
 python/dolfinx/fem/function.py        | 11 ++++++++++-
 python/dolfinx/mesh.py                | 18 ++++++++++++++++++
 python/dolfinx/wrappers/fem.cpp       |  1 +
 5 files changed, 40 insertions(+), 10 deletions(-)

diff --git a/cpp/dolfinx/refinement/CMakeLists.txt b/cpp/dolfinx/refinement/CMakeLists.txt
index 9560e45c823..0d9d997c490 100644
--- a/cpp/dolfinx/refinement/CMakeLists.txt
+++ b/cpp/dolfinx/refinement/CMakeLists.txt
@@ -1,14 +1,16 @@
 set(HEADERS_refinement
-    ${CMAKE_CURRENT_SOURCE_DIR}/dolfinx_refinement.h
-    ${CMAKE_CURRENT_SOURCE_DIR}/interval.h
-    ${CMAKE_CURRENT_SOURCE_DIR}/plaza.h
-    ${CMAKE_CURRENT_SOURCE_DIR}/refine.h
-    ${CMAKE_CURRENT_SOURCE_DIR}/utils.h
-    ${CMAKE_CURRENT_SOURCE_DIR}/option.h
-    PARENT_SCOPE
+  ${CMAKE_CURRENT_SOURCE_DIR}/dolfinx_refinement.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/interval.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/plaza.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/refine.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/utils.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/option.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/uniform.h
+  PARENT_SCOPE
 )
 
 target_sources(
   dolfinx PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/plaza.cpp
-                  ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/uniform.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp
 )
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index 42a5b6ef045..68e45482cab 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.21)
 
-project(dolfinx_nanobind LANGUAGES CXX)
+project(dolfinx_nanobind LANGUAGES CXX C)
 
 if(WIN32)
   # Windows requires all symbols to be manually exported. This flag exports all
diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index 0374a6fa24d..13325068825 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -36,7 +36,10 @@ class Constant(ufl.Constant):
     ]
 
     def __init__(
-        self, domain, c: typing.Union[np.ndarray, typing.Sequence, np.floating, np.complexfloating]
+        self,
+        domain,
+        c: typing.Union[np.ndarray, typing.Sequence, np.floating, np.complexfloating],
+        name: typing.Optional[str] = None,
     ):
         """A constant with respect to a domain.
 
@@ -60,6 +63,12 @@ def __init__(
         except AttributeError:
             raise AttributeError("Constant value must have a dtype attribute.")
 
+        # Set name
+        if name is None:
+            self.name = "c"
+        else:
+            self.name = name
+
     @property
     def value(self):
         """The value of the constant"""
diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index 46385f0d345..8062e637076 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -185,6 +185,24 @@ def index_map(self, dim: int) -> _cpp.common.IndexMap:
                 f"Call `dolfinx.mesh.Topology.create_entities({dim}) first."
             )
 
+    def index_maps(self, dim: int) -> list[_cpp.common.IndexMap]:
+        """Get the IndexMaps that describes the parallel distribution of the mesh entities,
+           for each entity type of the dimension.
+
+        Args:
+            dim: Topological dimension.
+
+        Returns:
+            List of IndexMaps for the entities of dimension ``dim``.
+        """
+        if (imaps := self._cpp_object.index_maps(dim)) is not None:
+            return imaps
+        else:
+            raise RuntimeError(
+                f"Entities of dimension {dim} have not been computed."
+                f"Call `dolfinx.mesh.Topology.create_entities({dim}) first."
+            )
+
     def interprocess_facets(self) -> npt.NDArray[np.int32]:
         """List of inter-process facets, if facet topology has been computed."""
         return self._cpp_object.interprocess_facets()
diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp
index 594624f86b5..0f71e0d1010 100644
--- a/python/dolfinx/wrappers/fem.cpp
+++ b/python/dolfinx/wrappers/fem.cpp
@@ -734,6 +734,7 @@ void declare_form(nb::module_& m, std::string type)
       .def_prop_ro("dtype", [](const dolfinx::fem::Form<T, U>&)
                    { return dolfinx_wrappers::numpy_dtype<T>(); })
       .def_prop_ro("coefficients", &dolfinx::fem::Form<T, U>::coefficients)
+      .def_prop_ro("constants", &dolfinx::fem::Form<T, U>::constants)
       .def_prop_ro("rank", &dolfinx::fem::Form<T, U>::rank)
       .def_prop_ro("mesh", &dolfinx::fem::Form<T, U>::mesh)
       .def_prop_ro("function_spaces",

From 1159a851664431d2593026d78c22ab815b974d0c Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Sat, 14 Jun 2025 21:10:25 +0100
Subject: [PATCH 02/28] Add files

---
 cpp/dolfinx/refinement/uniform.cpp | 186 +++++++++++++++++++++++++++++
 cpp/dolfinx/refinement/uniform.h   |   7 ++
 2 files changed, 193 insertions(+)
 create mode 100644 cpp/dolfinx/refinement/uniform.cpp
 create mode 100644 cpp/dolfinx/refinement/uniform.h

diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
new file mode 100644
index 00000000000..ff4d66ddc22
--- /dev/null
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -0,0 +1,186 @@
+#include <dolfinx/common/IndexMap.h>
+#include <dolfinx/common/Scatterer.h>
+#include <dolfinx/mesh/Mesh.h>
+#include <dolfinx/mesh/utils.h>
+#include <iterator>
+#include <vector>
+
+using namespace dolfinx;
+
+std::vector<std::int64_t>
+hex_subdivision(const std::vector<std::int64_t>& entities)
+{
+  int facet_list[8][3] = {{0, 1, 2}, {0, 1, 3}, {0, 2, 4}, {0, 3, 4},
+                          {1, 2, 5}, {1, 3, 5}, {2, 4, 5}, {3, 4, 5}};
+
+  int edge_list[8][3] = {{0, 1, 2}, {0, 3, 4},  {1, 5, 6},  {3, 5, 7},
+                         {2, 8, 9}, {4, 8, 10}, {6, 9, 11}, {7, 10, 11}};
+
+  std::vector<std::int64_t> topology;
+  for (int vi = 0; vi < 8; ++vi)
+  {
+    int edge_offset = 8;
+    int facet_offset = edge_offset + 12;
+    auto ee = edge_list[vi];
+    auto ff = facet_list[vi];
+    std::array<std::int64_t, 8> new_cell = {entities[vi],
+                                            entities[ee[1] + edge_offset],
+                                            entities[ee[0] + edge_offset],
+                                            entities[ff[0] + facet_offset],
+                                            entities[ee[2] + edge_offset],
+                                            entities[ff[2] + facet_offset],
+                                            entities[ff[1] + facet_offset],
+                                            entities.back()};
+    topology.insert(topology.end(), new_cell.begin(), new_cell.end());
+  }
+  return topology;
+}
+
+mesh::Mesh<double> uniform_refine(const mesh::Mesh<double>& mesh)
+{
+  // Requires edges and facets to be built already
+  auto topology = mesh.topology();
+
+  // Collect up entity types on each dimension
+  std::vector<std::vector<mesh::CellType>> entity_types;
+  for (int i = 0; i < 4; ++i)
+    entity_types.push_back(topology->entity_types(i));
+
+  // Get index maps for vertices and edges and create maps for new vertices
+  std::vector<std::vector<std::int64_t>> new_v;
+  std::vector<std::shared_ptr<const common::IndexMap>> index_maps;
+
+  // Indices for vertices and edges on dim 0 and dim 1.
+  std::vector<int> e_index = {0, 0};
+  // Check for quadrilateral faces and get index, if any.
+  if (auto it = std::find(entity_types[2].begin(), entity_types[2].end(),
+                          mesh::CellType::quadrilateral);
+      it != entity_types[2].end())
+    e_index.push_back(std::distance(entity_types[2].begin(), it));
+
+  // Check for hexahedral cells, and get index, if any.
+  if (auto it = std::find(entity_types[3].begin(), entity_types[3].end(),
+                          mesh::CellType::hexahedron);
+      it != entity_types[3].end())
+    e_index.push_back(std::distance(entity_types[3].begin(), it));
+
+  // Add up all local vertices, edges, quad facets and hex cells.
+  std::int64_t nlocal = 0;
+  for (std::size_t dim = 0; dim < e_index.size(); ++dim)
+  {
+    index_maps.push_back(topology->index_maps(dim)[e_index[dim]]);
+    new_v.push_back(std::vector<std::int64_t>(
+        index_maps.back()->size_local() + index_maps.back()->num_ghosts()));
+    nlocal += index_maps.back()->size_local();
+  }
+
+  // std::int64_t nv_global = index_maps[0]->size_global();
+
+  // Get current geometry and put into new array for vertices
+  std::vector<std::int32_t> vertex_to_x(index_maps[0]->size_local()
+                                        + index_maps[0]->num_ghosts());
+
+  // Iterate over cells
+  for (int j = 0; j < static_cast<int>(entity_types[3].size()); ++j)
+  {
+    // Get geometry for each cell type
+    auto x_dofmap = mesh.geometry().dofmap(j);
+    auto c_to_v = topology->connectivity({3, j}, {0, 0});
+    auto dof_layout = mesh.geometry().cmaps().at(j).create_dof_layout();
+    std::vector<int> entity_dofs(dof_layout.num_dofs());
+    for (int k = 0; k < dof_layout.num_dofs(); ++k)
+      entity_dofs[k] = dof_layout.entity_dofs(0, k)[0];
+
+    // Iterate over cells of this type
+    for (std::int32_t c = 0; c < topology->index_maps(3)[j]->size_local()
+                                     + topology->index_maps(3)[j]->num_ghosts();
+         ++c)
+    {
+      auto vertices = c_to_v->links(c);
+      for (std::size_t i = 0; i < vertices.size(); ++i)
+        vertex_to_x[vertices[i]] = x_dofmap(c, entity_dofs[i]);
+    }
+  }
+
+  // Copy existing vertices
+  std::vector<double> new_x(nlocal * 3);
+  auto x_g = mesh.geometry().x();
+  for (int i = 0; i < index_maps[0]->size_local(); ++i)
+    for (int j = 0; j < 3; ++j)
+      new_x[i * 3 + j] = x_g[vertex_to_x[i] * 3 + j];
+
+  // Get this process global offset and range
+  std::int64_t nscan;
+  MPI_Scan(&nlocal, &nscan, 1, MPI_INT64_T, MPI_SUM, mesh.comm());
+  std::array<std::int64_t, 2> local_range = {nscan - nlocal, nscan};
+
+  std::int32_t w_off = 0;
+  for (int j = 0; j < static_cast<int>(index_maps.size()); ++j)
+  {
+    std::int32_t num_entities = index_maps[j]->size_local();
+    std::iota(new_v[j].begin(), std::next(new_v[j].begin(), num_entities),
+              local_range[0] + w_off);
+
+    if (j > 0)
+    {
+      auto e_to_v = topology->connectivity({j, e_index[j]}, {0, 0});
+      for (std::int32_t w = 0; w < num_entities; ++w)
+      {
+        auto vt = e_to_v->links(w);
+        std::size_t nv_ent = vt.size();
+        std::array<double, 3> v = {0, 0, 0};
+        for (std::size_t i = 0; i < nv_ent; ++i)
+        {
+          for (int k = 0; k < 3; ++k)
+            v[k] += x_g[3 * vertex_to_x[vt[i]] + k];
+        }
+        for (int k = 0; k < 3; ++k)
+          new_x[(w + w_off) * 3 + k] = v[k] / static_cast<double>(nv_ent);
+      }
+    }
+    common::Scatterer sc(*index_maps[j], 1);
+    sc.scatter_fwd(std::span<const std::int64_t>(new_v[j]),
+                   std::span(std::next(new_v[j].begin(), num_entities),
+                             index_maps[j]->num_ghosts()));
+    w_off += num_entities;
+  }
+
+  // Create new topology...
+
+  std::vector<std::vector<std::int64_t>> mixed_topology(entity_types[3].size());
+  for (int k = 0; k < static_cast<int>(entity_types[3].size()); ++k)
+  {
+    auto c_to_v = topology->connectivity({3, k}, {0, 0});
+    auto c_to_e = topology->connectivity({3, k}, {1, 0});
+
+    for (int c = 0; c < topology->index_maps(3)[k]->size_local(); ++c)
+    {
+      std::vector<std::int64_t> entities;
+      for (std::int32_t i : c_to_v->links(c))
+        entities.push_back(new_v[0][i]);
+      for (std::int32_t i : c_to_e->links(c))
+        entities.push_back(new_v[1][i]);
+      if (e_index.size() > 2)
+      {
+        for (std::int32_t i :
+             topology->connectivity({3, k}, {2, e_index[2]})->links(c))
+          entities.push_back(new_v[2][i]);
+      }
+      if (e_index.size() > 3)
+        entities.push_back(new_v[3][c]);
+
+      if (entity_types[3][k] == mesh::CellType::hexahedron)
+      {
+        auto new_cells = hex_subdivision(entities);
+        mixed_topology[k].insert(mixed_topology[k].end(), new_cells.begin(),
+                                 new_cells.end());
+      }
+    }
+  }
+
+  auto partitioner = mesh::create_cell_partitioner(mesh::GhostMode::none);
+  mesh::Mesh new_mesh = mesh::create_mesh(mesh.comm(),
+      mesh.comm(), mixed_topology, mesh.geometry().cmaps(), mesh.comm(), new_x, {new_x.size()/3, 3}, partitioner);
+
+  return new_mesh;
+}
diff --git a/cpp/dolfinx/refinement/uniform.h b/cpp/dolfinx/refinement/uniform.h
new file mode 100644
index 00000000000..422c4cdfc90
--- /dev/null
+++ b/cpp/dolfinx/refinement/uniform.h
@@ -0,0 +1,7 @@
+#include <dolfinx/mesh/Mesh.h>
+
+#pragma once
+
+/// @brief Uniform refinement of a 3D mesh containing hex, tet, prism or pyramids.
+/// @param mesh Input mesh
+mesh::Mesh<double> uniform_refine(const mesh::Mesh<double>& mesh);

From 6b0502c6daefd383822ec6839da68d1fd7133ad2 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Tue, 17 Jun 2025 13:16:13 +0100
Subject: [PATCH 03/28] working?

---
 cpp/dolfinx/refinement/uniform.cpp | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
index ff4d66ddc22..8e24abb2b8d 100644
--- a/cpp/dolfinx/refinement/uniform.cpp
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -179,8 +179,14 @@ mesh::Mesh<double> uniform_refine(const mesh::Mesh<double>& mesh)
   }
 
   auto partitioner = mesh::create_cell_partitioner(mesh::GhostMode::none);
-  mesh::Mesh new_mesh = mesh::create_mesh(mesh.comm(),
-      mesh.comm(), mixed_topology, mesh.geometry().cmaps(), mesh.comm(), new_x, {new_x.size()/3, 3}, partitioner);
+
+  std::vector<std::span<const std::int64_t>> topo_span;
+  for (auto v : mixed_topology)
+    topo_span.push_back(std::span(v));
+
+  mesh::Mesh new_mesh = mesh::create_mesh(
+      mesh.comm(), mesh.comm(), topo_span, mesh.geometry().cmaps(), mesh.comm(),
+      new_x, {new_x.size() / 3, 3}, partitioner);
 
   return new_mesh;
 }

From 08a0287b20b303da8b493eaf703daf3de88e7039 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Tue, 17 Jun 2025 14:26:11 +0100
Subject: [PATCH 04/28] debug

---
 cpp/dolfinx/refinement/dolfinx_refinement.h |  1 +
 cpp/dolfinx/refinement/uniform.cpp          | 11 ++++++++++-
 cpp/dolfinx/refinement/uniform.h            |  5 +++++
 python/dolfinx/wrappers/refinement.cpp      |  3 +++
 4 files changed, 19 insertions(+), 1 deletion(-)

diff --git a/cpp/dolfinx/refinement/dolfinx_refinement.h b/cpp/dolfinx/refinement/dolfinx_refinement.h
index 46588b267b4..802cdc2eff6 100644
--- a/cpp/dolfinx/refinement/dolfinx_refinement.h
+++ b/cpp/dolfinx/refinement/dolfinx_refinement.h
@@ -12,3 +12,4 @@ namespace dolfinx::refinement
 
 #include <dolfinx/refinement/interval.h>
 #include <dolfinx/refinement/refine.h>
+#include <dolfinx/refinement/uniform.h>
diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
index 8e24abb2b8d..42e9c2f5616 100644
--- a/cpp/dolfinx/refinement/uniform.cpp
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -1,3 +1,4 @@
+
 #include <dolfinx/common/IndexMap.h>
 #include <dolfinx/common/Scatterer.h>
 #include <dolfinx/mesh/Mesh.h>
@@ -5,6 +6,8 @@
 #include <iterator>
 #include <vector>
 
+#include "uniform.h"
+
 using namespace dolfinx;
 
 std::vector<std::int64_t>
@@ -36,7 +39,7 @@ hex_subdivision(const std::vector<std::int64_t>& entities)
   return topology;
 }
 
-mesh::Mesh<double> uniform_refine(const mesh::Mesh<double>& mesh)
+mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
 {
   // Requires edges and facets to be built already
   auto topology = mesh.topology();
@@ -182,7 +185,13 @@ mesh::Mesh<double> uniform_refine(const mesh::Mesh<double>& mesh)
 
   std::vector<std::span<const std::int64_t>> topo_span;
   for (auto v : mixed_topology)
+  {
+    std::cout << "[";
     topo_span.push_back(std::span(v));
+    for (auto q : topo_span.back())
+      std::cout << q << " ";
+    std::cout << "]\n";
+  }
 
   mesh::Mesh new_mesh = mesh::create_mesh(
       mesh.comm(), mesh.comm(), topo_span, mesh.geometry().cmaps(), mesh.comm(),
diff --git a/cpp/dolfinx/refinement/uniform.h b/cpp/dolfinx/refinement/uniform.h
index 422c4cdfc90..7ddafbbfa79 100644
--- a/cpp/dolfinx/refinement/uniform.h
+++ b/cpp/dolfinx/refinement/uniform.h
@@ -2,6 +2,11 @@
 
 #pragma once
 
+namespace dolfinx::refinement
+{
+
 /// @brief Uniform refinement of a 3D mesh containing hex, tet, prism or pyramids.
 /// @param mesh Input mesh
 mesh::Mesh<double> uniform_refine(const mesh::Mesh<double>& mesh);
+
+} // namespace dolfinx::refinement
diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp
index 1e7c73a8e96..10b8860eb83 100644
--- a/python/dolfinx/wrappers/refinement.cpp
+++ b/python/dolfinx/wrappers/refinement.cpp
@@ -14,6 +14,7 @@
 #include <dolfinx/refinement/interval.h>
 #include <dolfinx/refinement/option.h>
 #include <dolfinx/refinement/refine.h>
+#include <dolfinx/refinement/uniform.h>
 #include <dolfinx/refinement/utils.h>
 #include <functional>
 #include <nanobind/nanobind.h>
@@ -33,6 +34,8 @@ namespace
 template <std::floating_point T>
 void export_refinement(nb::module_& m)
 {
+  m.def("uniform_refine", dolfinx::refinement::uniform_refine);
+
   m.def(
       "refine",
       [](const dolfinx::mesh::Mesh<T>& mesh,

From 0b3ee95d5a64dccdbcad4f8facb7e756e8390152 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Tue, 17 Jun 2025 14:29:58 +0100
Subject: [PATCH 05/28] try to fix formatitng

---
 python/dolfinx/mesh.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index d3533efef85..8109794ec8e 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -190,8 +190,8 @@ def index_map(self, dim: int) -> _cpp.common.IndexMap:
             )
 
     def index_maps(self, dim: int) -> list[_cpp.common.IndexMap]:
-        """Get the IndexMaps that describes the parallel distribution of the mesh entities,
-           for each entity type of the dimension.
+        """Get the IndexMaps that describes the parallel distribution of
+           the mesh entities, for each entity type of the dimension.
 
         Args:
             dim: Topological dimension.

From d7047aed3cef098f9b1a795271ab600d60429eb8 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Tue, 17 Jun 2025 14:32:10 +0100
Subject: [PATCH 06/28] clang-format

---
 cpp/dolfinx/refinement/uniform.h | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/cpp/dolfinx/refinement/uniform.h b/cpp/dolfinx/refinement/uniform.h
index 7ddafbbfa79..4b286bfbf46 100644
--- a/cpp/dolfinx/refinement/uniform.h
+++ b/cpp/dolfinx/refinement/uniform.h
@@ -5,7 +5,8 @@
 namespace dolfinx::refinement
 {
 
-/// @brief Uniform refinement of a 3D mesh containing hex, tet, prism or pyramids.
+/// @brief Uniform refinement of a 3D mesh containing hex, tet, prism or
+/// pyramids.
 /// @param mesh Input mesh
 mesh::Mesh<double> uniform_refine(const mesh::Mesh<double>& mesh);
 

From 645d55c2fc605c4ed478261909e4bb5df8ee6462 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Tue, 17 Jun 2025 15:25:16 +0100
Subject: [PATCH 07/28] Fix up span

---
 cpp/dolfinx/refinement/uniform.cpp | 14 +++++---------
 1 file changed, 5 insertions(+), 9 deletions(-)

diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
index 42e9c2f5616..dc41b7a4462 100644
--- a/cpp/dolfinx/refinement/uniform.cpp
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -10,6 +10,8 @@
 
 using namespace dolfinx;
 
+namespace
+{
 std::vector<std::int64_t>
 hex_subdivision(const std::vector<std::int64_t>& entities)
 {
@@ -38,6 +40,7 @@ hex_subdivision(const std::vector<std::int64_t>& entities)
   }
   return topology;
 }
+} // namespace
 
 mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
 {
@@ -183,15 +186,8 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
 
   auto partitioner = mesh::create_cell_partitioner(mesh::GhostMode::none);
 
-  std::vector<std::span<const std::int64_t>> topo_span;
-  for (auto v : mixed_topology)
-  {
-    std::cout << "[";
-    topo_span.push_back(std::span(v));
-    for (auto q : topo_span.back())
-      std::cout << q << " ";
-    std::cout << "]\n";
-  }
+  std::vector<std::span<const std::int64_t>> topo_span(mixed_topology.begin(),
+                                                       mixed_topology.end());
 
   mesh::Mesh new_mesh = mesh::create_mesh(
       mesh.comm(), mesh.comm(), topo_span, mesh.geometry().cmaps(), mesh.comm(),

From 8e69f2e93ad3c169ee3c238dde01d612e27ba6bc Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Tue, 17 Jun 2025 15:35:12 +0100
Subject: [PATCH 08/28] Add refinement for tets

---
 cpp/dolfinx/refinement/uniform.cpp | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
index dc41b7a4462..f30a6c8c087 100644
--- a/cpp/dolfinx/refinement/uniform.cpp
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -12,6 +12,20 @@ using namespace dolfinx;
 
 namespace
 {
+
+std::vector<std::int64_t>
+tet_subdivision(const std::vector<std::int64_t>& entities)
+{
+  std::vector<std::int64_t> topology;
+
+  int cell_list[32] = {0, 7, 8, 9, 1, 5, 6, 9, 2, 4, 6, 8, 3, 4, 5, 7,
+                       9, 4, 6, 8, 9, 4, 8, 7, 9, 4, 7, 5, 9, 4, 5, 6};
+
+  for (int i = 0; i < 32; ++i)
+    topology.push_back(entities[cell_list[i]]);
+  return topology;
+}
+
 std::vector<std::int64_t>
 hex_subdivision(const std::vector<std::int64_t>& entities)
 {
@@ -181,6 +195,12 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
         mixed_topology[k].insert(mixed_topology[k].end(), new_cells.begin(),
                                  new_cells.end());
       }
+      else if (entity_types[3][k] == mesh::CellType::tetrahedron)
+      {
+        auto new_cells = tet_subdivision(entities);
+        mixed_topology[k].insert(mixed_topology[k].end(), new_cells.begin(),
+                                 new_cells.end());
+      }
     }
   }
 

From 60dd5bc88de4a22e4bac4ad59168ee1a999d783f Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Tue, 17 Jun 2025 15:53:07 +0100
Subject: [PATCH 09/28] Add prism cells

---
 cpp/dolfinx/refinement/uniform.cpp | 42 +++++++++++++++++++++---------
 1 file changed, 30 insertions(+), 12 deletions(-)

diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
index f30a6c8c087..d182907d4c9 100644
--- a/cpp/dolfinx/refinement/uniform.cpp
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -26,6 +26,19 @@ tet_subdivision(const std::vector<std::int64_t>& entities)
   return topology;
 }
 
+std::vector<std::int64_t>
+prism_subdivision(const std::vector<std::int64_t>& entities)
+{
+  std::vector<std::int64_t> topology;
+  int cell_list[48]
+      = {0,  6,  7,  8,  15, 16, 6,  1,  9,  15, 10, 17, 7,  9,  2,  16,
+         17, 11, 6,  9,  7,  15, 17, 16, 15, 17, 16, 12, 14, 13, 8,  15,
+         16, 3,  12, 13, 11, 17, 16, 5,  14, 13, 10, 15, 17, 4,  12, 14};
+  for (int i = 0; i < 48; ++i)
+    topology.push_back(entities[cell_list[i]]);
+  return topology;
+}
+
 std::vector<std::int64_t>
 hex_subdivision(const std::vector<std::int64_t>& entities)
 {
@@ -168,6 +181,20 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
   // Create new topology...
 
   std::vector<std::vector<std::int64_t>> mixed_topology(entity_types[3].size());
+
+  std::vector<std::function<std::vector<std::int64_t>(
+      const std::vector<std::int64_t>&)>>
+      subdiv;
+  for (int k = 0; k < static_cast<int>(entity_types[3].size()); ++k)
+  {
+    if (entity_types[3][k] == mesh::CellType::hexahedron)
+      subdiv.push_back(hex_subdivision);
+    if (entity_types[3][k] == mesh::CellType::prism)
+      subdiv.push_back(prism_subdivision);
+    if (entity_types[3][k] == mesh::CellType::tetrahedron)
+      subdiv.push_back(tet_subdivision);
+  }
+
   for (int k = 0; k < static_cast<int>(entity_types[3].size()); ++k)
   {
     auto c_to_v = topology->connectivity({3, k}, {0, 0});
@@ -189,18 +216,9 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
       if (e_index.size() > 3)
         entities.push_back(new_v[3][c]);
 
-      if (entity_types[3][k] == mesh::CellType::hexahedron)
-      {
-        auto new_cells = hex_subdivision(entities);
-        mixed_topology[k].insert(mixed_topology[k].end(), new_cells.begin(),
-                                 new_cells.end());
-      }
-      else if (entity_types[3][k] == mesh::CellType::tetrahedron)
-      {
-        auto new_cells = tet_subdivision(entities);
-        mixed_topology[k].insert(mixed_topology[k].end(), new_cells.begin(),
-                                 new_cells.end());
-      }
+      auto new_cells = subdiv[k](entities);
+      mixed_topology[k].insert(mixed_topology[k].end(), new_cells.begin(),
+                               new_cells.end());
     }
   }
 

From 921aafad8b29baf39858a47d4ffa20618f9d5912 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Tue, 17 Jun 2025 16:52:01 +0100
Subject: [PATCH 10/28] Add pyramid cells

---
 cpp/dolfinx/refinement/uniform.cpp | 31 +++++++++++++++++++++++++++---
 1 file changed, 28 insertions(+), 3 deletions(-)

diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
index d182907d4c9..f3775bfd65e 100644
--- a/cpp/dolfinx/refinement/uniform.cpp
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -13,6 +13,24 @@ using namespace dolfinx;
 namespace
 {
 
+std::vector<std::int64_t>
+pyr_subdivision(const std::vector<std::int64_t>& entities)
+{
+  int cells_pyr[25] = {0,  5,  6, 13, 7,  1,  8,  5, 13, 9,  3,  10, 8,
+                       13, 12, 2, 6,  10, 13, 11, 7, 9,  11, 12, 4};
+
+  //  topo_tet
+  //      = [[edges [0], facets [0], edges [2], edges [4]],
+  //         [edges [1], facets [0], edges [6], edges [2]],
+  //         [edges [5], facets [0], edges [7], edges [6]],
+  //         [edges [3], facets [0], edges [4], edges [7]]]
+
+  std::vector<std::int64_t> topology(25);
+  for (int i = 0; i < 25; ++i)
+    topology[i] = entities[cells_pyr[i]];
+  return topology;
+}
+
 std::vector<std::int64_t>
 tet_subdivision(const std::vector<std::int64_t>& entities)
 {
@@ -193,6 +211,8 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
       subdiv.push_back(prism_subdivision);
     if (entity_types[3][k] == mesh::CellType::tetrahedron)
       subdiv.push_back(tet_subdivision);
+    if (entity_types[3][k] == mesh::CellType::pyramid)
+      subdiv.push_back(pyr_subdivision);
   }
 
   for (int k = 0; k < static_cast<int>(entity_types[3].size()); ++k)
@@ -209,9 +229,14 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
         entities.push_back(new_v[1][i]);
       if (e_index.size() > 2)
       {
-        for (std::int32_t i :
-             topology->connectivity({3, k}, {2, e_index[2]})->links(c))
-          entities.push_back(new_v[2][i]);
+        auto conn = topology->connectivity({3, k}, {2, e_index[2]});
+        if (conn)
+        {
+          spdlog::debug("Get connectivity from (3,{})->(2,{})", k, e_index[2]);
+          for (std::int32_t i :
+               topology->connectivity({3, k}, {2, e_index[2]})->links(c))
+            entities.push_back(new_v[2][i]);
+        }
       }
       if (e_index.size() > 3)
         entities.push_back(new_v[3][c]);

From 7cd7e361c72eeff6ae6f8a28fb83c678cbb976c0 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Thu, 19 Jun 2025 11:30:26 +0200
Subject: [PATCH 11/28] Simplify to lists of int for each cell type

---
 cpp/dolfinx/refinement/uniform.cpp | 143 ++++++++++-------------------
 1 file changed, 50 insertions(+), 93 deletions(-)

diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
index f3775bfd65e..6b9b0248eed 100644
--- a/cpp/dolfinx/refinement/uniform.cpp
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -10,83 +10,6 @@
 
 using namespace dolfinx;
 
-namespace
-{
-
-std::vector<std::int64_t>
-pyr_subdivision(const std::vector<std::int64_t>& entities)
-{
-  int cells_pyr[25] = {0,  5,  6, 13, 7,  1,  8,  5, 13, 9,  3,  10, 8,
-                       13, 12, 2, 6,  10, 13, 11, 7, 9,  11, 12, 4};
-
-  //  topo_tet
-  //      = [[edges [0], facets [0], edges [2], edges [4]],
-  //         [edges [1], facets [0], edges [6], edges [2]],
-  //         [edges [5], facets [0], edges [7], edges [6]],
-  //         [edges [3], facets [0], edges [4], edges [7]]]
-
-  std::vector<std::int64_t> topology(25);
-  for (int i = 0; i < 25; ++i)
-    topology[i] = entities[cells_pyr[i]];
-  return topology;
-}
-
-std::vector<std::int64_t>
-tet_subdivision(const std::vector<std::int64_t>& entities)
-{
-  std::vector<std::int64_t> topology;
-
-  int cell_list[32] = {0, 7, 8, 9, 1, 5, 6, 9, 2, 4, 6, 8, 3, 4, 5, 7,
-                       9, 4, 6, 8, 9, 4, 8, 7, 9, 4, 7, 5, 9, 4, 5, 6};
-
-  for (int i = 0; i < 32; ++i)
-    topology.push_back(entities[cell_list[i]]);
-  return topology;
-}
-
-std::vector<std::int64_t>
-prism_subdivision(const std::vector<std::int64_t>& entities)
-{
-  std::vector<std::int64_t> topology;
-  int cell_list[48]
-      = {0,  6,  7,  8,  15, 16, 6,  1,  9,  15, 10, 17, 7,  9,  2,  16,
-         17, 11, 6,  9,  7,  15, 17, 16, 15, 17, 16, 12, 14, 13, 8,  15,
-         16, 3,  12, 13, 11, 17, 16, 5,  14, 13, 10, 15, 17, 4,  12, 14};
-  for (int i = 0; i < 48; ++i)
-    topology.push_back(entities[cell_list[i]]);
-  return topology;
-}
-
-std::vector<std::int64_t>
-hex_subdivision(const std::vector<std::int64_t>& entities)
-{
-  int facet_list[8][3] = {{0, 1, 2}, {0, 1, 3}, {0, 2, 4}, {0, 3, 4},
-                          {1, 2, 5}, {1, 3, 5}, {2, 4, 5}, {3, 4, 5}};
-
-  int edge_list[8][3] = {{0, 1, 2}, {0, 3, 4},  {1, 5, 6},  {3, 5, 7},
-                         {2, 8, 9}, {4, 8, 10}, {6, 9, 11}, {7, 10, 11}};
-
-  std::vector<std::int64_t> topology;
-  for (int vi = 0; vi < 8; ++vi)
-  {
-    int edge_offset = 8;
-    int facet_offset = edge_offset + 12;
-    auto ee = edge_list[vi];
-    auto ff = facet_list[vi];
-    std::array<std::int64_t, 8> new_cell = {entities[vi],
-                                            entities[ee[1] + edge_offset],
-                                            entities[ee[0] + edge_offset],
-                                            entities[ff[0] + facet_offset],
-                                            entities[ee[2] + edge_offset],
-                                            entities[ff[2] + facet_offset],
-                                            entities[ff[1] + facet_offset],
-                                            entities.back()};
-    topology.insert(topology.end(), new_cell.begin(), new_cell.end());
-  }
-  return topology;
-}
-} // namespace
-
 mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
 {
   // Requires edges and facets to be built already
@@ -200,23 +123,52 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
 
   std::vector<std::vector<std::int64_t>> mixed_topology(entity_types[3].size());
 
-  std::vector<std::function<std::vector<std::int64_t>(
-      const std::vector<std::int64_t>&)>>
-      subdiv;
+  // Find index of tets in topology list, if any
+  int ktet = -1;
+  auto it = std::find(entity_types[3].begin(), entity_types[3].end(),
+                      mesh::CellType::tetrahedron);
+  if (it != entity_types[3].end())
+    ktet = std::distance(entity_types[3].begin(), it);
+  // Topology for tetrahedra which arise from pyramid subdivision
+  std::vector<int> pyr_to_tet_list
+      = {5, 13, 7, 9, 6, 13, 11, 7, 10, 13, 12, 11, 8, 13, 9, 12};
+
+  std::vector<int> refined_cell_list;
   for (int k = 0; k < static_cast<int>(entity_types[3].size()); ++k)
   {
+    // Reserve an estimate of space for the topology of each type
+    mixed_topology[k].reserve(mesh.topology()->index_maps(3)[k]->size_local()
+                              * 8 * 6);
+
+    // Select correct subdivision for celltype
+    // Hex -> 8 hex, Prism -> 8 prism, Tet -> 8 tet, Pyr -> 5 pyr + 4 tet
     if (entity_types[3][k] == mesh::CellType::hexahedron)
-      subdiv.push_back(hex_subdivision);
-    if (entity_types[3][k] == mesh::CellType::prism)
-      subdiv.push_back(prism_subdivision);
-    if (entity_types[3][k] == mesh::CellType::tetrahedron)
-      subdiv.push_back(tet_subdivision);
-    if (entity_types[3][k] == mesh::CellType::pyramid)
-      subdiv.push_back(pyr_subdivision);
-  }
+      refined_cell_list
+          = {0, 9,  8,  20, 10, 22, 21, 26, 1, 11, 8,  20, 12, 23, 21, 26,
+             2, 13, 9,  20, 14, 24, 22, 26, 3, 13, 11, 20, 15, 24, 23, 26,
+             4, 16, 10, 21, 17, 25, 22, 26, 5, 16, 12, 21, 18, 25, 23, 26,
+             6, 17, 14, 22, 19, 25, 24, 26, 7, 18, 15, 23, 19, 25, 24, 26};
+    else if (entity_types[3][k] == mesh::CellType::tetrahedron)
+    {
+      refined_cell_list = {0, 7, 8, 9, 1, 5, 6, 9, 2, 4, 6, 8, 3, 4, 5, 7,
+                           9, 4, 6, 8, 9, 4, 8, 7, 9, 4, 7, 5, 9, 4, 5, 6};
+    }
+    else if (entity_types[3][k] == mesh::CellType::prism)
+    {
+      refined_cell_list
+          = {0,  6,  7,  8,  15, 16, 6,  1,  9,  15, 10, 17, 7,  9,  2,  16,
+             17, 11, 6,  9,  7,  15, 17, 16, 15, 17, 16, 12, 14, 13, 8,  15,
+             16, 3,  12, 13, 11, 17, 16, 5,  14, 13, 10, 15, 17, 4,  12, 14};
+    }
+    else if (entity_types[3][k] == mesh::CellType::pyramid)
+    {
+      refined_cell_list = {0,  5,  6, 13, 7,  1,  8,  5, 13, 9,  3,  10, 8,
+                           13, 12, 2, 6,  10, 13, 11, 7, 9,  11, 12, 4};
+      if (ktet == -1)
+        throw std::runtime_error("Cannot refine mesh with pyramids and no "
+                                 "tetrahedra.");
+    }
 
-  for (int k = 0; k < static_cast<int>(entity_types[3].size()); ++k)
-  {
     auto c_to_v = topology->connectivity({3, k}, {0, 0});
     auto c_to_e = topology->connectivity({3, k}, {1, 0});
 
@@ -241,9 +193,14 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
       if (e_index.size() > 3)
         entities.push_back(new_v[3][c]);
 
-      auto new_cells = subdiv[k](entities);
-      mixed_topology[k].insert(mixed_topology[k].end(), new_cells.begin(),
-                               new_cells.end());
+      for (int i : refined_cell_list)
+        mixed_topology[k].push_back(entities[i]);
+
+      if (entity_types[3][k] == mesh::CellType::pyramid)
+      {
+        for (int i : pyr_to_tet_list)
+          mixed_topology[ktet].push_back(entities[i]);
+      }
     }
   }
 

From 2e1ece01e621bffed7925d432a42568ac62655d7 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Thu, 19 Jun 2025 11:48:30 +0200
Subject: [PATCH 12/28] Throw error for missing topology

---
 cpp/dolfinx/mesh/Topology.cpp      | 4 ++--
 cpp/dolfinx/refinement/uniform.cpp | 4 ++++
 2 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/cpp/dolfinx/mesh/Topology.cpp b/cpp/dolfinx/mesh/Topology.cpp
index c4249e7aa16..88c4c1c2566 100644
--- a/cpp/dolfinx/mesh/Topology.cpp
+++ b/cpp/dolfinx/mesh/Topology.cpp
@@ -814,8 +814,8 @@ Topology::index_maps(int dim) const
   for (std::size_t i = 0; i < _entity_types[dim].size(); ++i)
   {
     auto it = _index_maps.find({dim, int(i)});
-    assert(it != _index_maps.end());
-    maps.push_back(it->second);
+    if (it != _index_maps.end())
+      maps.push_back(it->second);
   }
   return maps;
 }
diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
index 6b9b0248eed..856075ec36b 100644
--- a/cpp/dolfinx/refinement/uniform.cpp
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -42,6 +42,10 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
   std::int64_t nlocal = 0;
   for (std::size_t dim = 0; dim < e_index.size(); ++dim)
   {
+    if (topology->index_maps(dim).empty())
+      throw std::runtime_error(
+          "Missing topology for dimension " + std::to_string(dim)
+          + ", need to call create_entities(" + std::to_string(dim) + ")");
     index_maps.push_back(topology->index_maps(dim)[e_index[dim]]);
     new_v.push_back(std::vector<std::int64_t>(
         index_maps.back()->size_local() + index_maps.back()->num_ghosts()));

From dbfb94c12cd41e3cc00471529bc99c822a4ee13d Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Thu, 19 Jun 2025 12:47:49 +0200
Subject: [PATCH 13/28] Add basic test

---
 python/test/unit/refinement/test_uniform.py | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)
 create mode 100644 python/test/unit/refinement/test_uniform.py

diff --git a/python/test/unit/refinement/test_uniform.py b/python/test/unit/refinement/test_uniform.py
new file mode 100644
index 00000000000..0397fdda273
--- /dev/null
+++ b/python/test/unit/refinement/test_uniform.py
@@ -0,0 +1,20 @@
+from mpi4py import MPI
+
+import pytest
+
+import dolfinx
+from dolfinx.mesh import CellType, Mesh, create_unit_cube
+
+
+@pytest.mark.parametrize("ctype", [CellType.hexahedron, CellType.tetrahedron, CellType.prism])
+def test_uniform_refinement(ctype):
+    mesh = create_unit_cube(MPI.COMM_WORLD, 3, 2, 1, cell_type=ctype)
+
+    ncells0 = mesh.topology.index_map(3).size_local
+    mesh.topology.create_entities(1)
+    mesh.topology.create_entities(2)
+
+    m2 = Mesh(dolfinx.cpp.refinement.uniform_refine(mesh._cpp_object), None)
+    ncells1 = m2.topology.index_map(3).size_local
+
+    assert ncells0 * 8 == ncells1

From 8eb109cebf2c1f01d66337cec683a2228d8c6649 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Thu, 19 Jun 2025 14:45:11 +0200
Subject: [PATCH 14/28] working for 2D too

---
 cpp/dolfinx/refinement/uniform.cpp          | 96 ++++++++++++++-------
 python/test/unit/refinement/test_uniform.py | 19 +++-
 2 files changed, 81 insertions(+), 34 deletions(-)

diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
index 856075ec36b..4765d2effec 100644
--- a/cpp/dolfinx/refinement/uniform.cpp
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -12,12 +12,17 @@ using namespace dolfinx;
 
 mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
 {
-  // Requires edges and facets to be built already
+  // Requires edges (and facets for some 3D meshes) to be built already
   auto topology = mesh.topology();
+  int tdim = topology->dim();
+  if (tdim < 2)
+    throw std::runtime_error("Uniform refinement only for 2D and 3D meshes");
+
+  spdlog::info("Topology dim = {}", tdim);
 
   // Collect up entity types on each dimension
   std::vector<std::vector<mesh::CellType>> entity_types;
-  for (int i = 0; i < 4; ++i)
+  for (int i = 0; i < tdim + 1; ++i)
     entity_types.push_back(topology->entity_types(i));
 
   // Get index maps for vertices and edges and create maps for new vertices
@@ -26,17 +31,21 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
 
   // Indices for vertices and edges on dim 0 and dim 1.
   std::vector<int> e_index = {0, 0};
+
   // Check for quadrilateral faces and get index, if any.
   if (auto it = std::find(entity_types[2].begin(), entity_types[2].end(),
                           mesh::CellType::quadrilateral);
       it != entity_types[2].end())
     e_index.push_back(std::distance(entity_types[2].begin(), it));
 
-  // Check for hexahedral cells, and get index, if any.
-  if (auto it = std::find(entity_types[3].begin(), entity_types[3].end(),
-                          mesh::CellType::hexahedron);
-      it != entity_types[3].end())
-    e_index.push_back(std::distance(entity_types[3].begin(), it));
+  if (tdim == 3)
+  {
+    // In 3D, check for hexahedral cells, and get index, if any.
+    if (auto it = std::find(entity_types[3].begin(), entity_types[3].end(),
+                            mesh::CellType::hexahedron);
+        it != entity_types[3].end())
+      e_index.push_back(std::distance(entity_types[3].begin(), it));
+  }
 
   // Add up all local vertices, edges, quad facets and hex cells.
   std::int64_t nlocal = 0;
@@ -59,19 +68,20 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
                                         + index_maps[0]->num_ghosts());
 
   // Iterate over cells
-  for (int j = 0; j < static_cast<int>(entity_types[3].size()); ++j)
+  for (int j = 0; j < static_cast<int>(entity_types.back().size()); ++j)
   {
     // Get geometry for each cell type
     auto x_dofmap = mesh.geometry().dofmap(j);
-    auto c_to_v = topology->connectivity({3, j}, {0, 0});
+    auto c_to_v = topology->connectivity({tdim, j}, {0, 0});
     auto dof_layout = mesh.geometry().cmaps().at(j).create_dof_layout();
     std::vector<int> entity_dofs(dof_layout.num_dofs());
     for (int k = 0; k < dof_layout.num_dofs(); ++k)
       entity_dofs[k] = dof_layout.entity_dofs(0, k)[0];
 
     // Iterate over cells of this type
-    for (std::int32_t c = 0; c < topology->index_maps(3)[j]->size_local()
-                                     + topology->index_maps(3)[j]->num_ghosts();
+    for (std::int32_t c = 0;
+         c < topology->index_maps(tdim)[j]->size_local()
+                 + topology->index_maps(tdim)[j]->num_ghosts();
          ++c)
     {
       auto vertices = c_to_v->links(c);
@@ -125,58 +135,75 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
 
   // Create new topology...
 
-  std::vector<std::vector<std::int64_t>> mixed_topology(entity_types[3].size());
+  std::vector<std::vector<std::int64_t>> mixed_topology(
+      entity_types.back().size());
 
   // Find index of tets in topology list, if any
   int ktet = -1;
-  auto it = std::find(entity_types[3].begin(), entity_types[3].end(),
+  auto it = std::find(entity_types.back().begin(), entity_types.back().end(),
                       mesh::CellType::tetrahedron);
-  if (it != entity_types[3].end())
-    ktet = std::distance(entity_types[3].begin(), it);
+  if (it != entity_types.back().end())
+    ktet = std::distance(entity_types.back().begin(), it);
   // Topology for tetrahedra which arise from pyramid subdivision
   std::vector<int> pyr_to_tet_list
       = {5, 13, 7, 9, 6, 13, 11, 7, 10, 13, 12, 11, 8, 13, 9, 12};
 
   std::vector<int> refined_cell_list;
-  for (int k = 0; k < static_cast<int>(entity_types[3].size()); ++k)
+  for (int k = 0; k < static_cast<int>(entity_types.back().size()); ++k)
   {
     // Reserve an estimate of space for the topology of each type
-    mixed_topology[k].reserve(mesh.topology()->index_maps(3)[k]->size_local()
+    mixed_topology[k].reserve(mesh.topology()->index_maps(tdim)[k]->size_local()
                               * 8 * 6);
 
     // Select correct subdivision for celltype
     // Hex -> 8 hex, Prism -> 8 prism, Tet -> 8 tet, Pyr -> 5 pyr + 4 tet
-    if (entity_types[3][k] == mesh::CellType::hexahedron)
+    if (entity_types.back()[k] == mesh::CellType::hexahedron)
+    {
+      spdlog::info("Hex subdivision");
       refined_cell_list
           = {0, 9,  8,  20, 10, 22, 21, 26, 1, 11, 8,  20, 12, 23, 21, 26,
              2, 13, 9,  20, 14, 24, 22, 26, 3, 13, 11, 20, 15, 24, 23, 26,
              4, 16, 10, 21, 17, 25, 22, 26, 5, 16, 12, 21, 18, 25, 23, 26,
              6, 17, 14, 22, 19, 25, 24, 26, 7, 18, 15, 23, 19, 25, 24, 26};
-    else if (entity_types[3][k] == mesh::CellType::tetrahedron)
+    }
+    else if (entity_types.back()[k] == mesh::CellType::tetrahedron)
     {
+      spdlog::info("Tet subdivision");
       refined_cell_list = {0, 7, 8, 9, 1, 5, 6, 9, 2, 4, 6, 8, 3, 4, 5, 7,
                            9, 4, 6, 8, 9, 4, 8, 7, 9, 4, 7, 5, 9, 4, 5, 6};
     }
-    else if (entity_types[3][k] == mesh::CellType::prism)
+    else if (entity_types.back()[k] == mesh::CellType::prism)
     {
+      spdlog::info("Prism subdivision");
       refined_cell_list
           = {0,  6,  7,  8,  15, 16, 6,  1,  9,  15, 10, 17, 7,  9,  2,  16,
              17, 11, 6,  9,  7,  15, 17, 16, 15, 17, 16, 12, 14, 13, 8,  15,
              16, 3,  12, 13, 11, 17, 16, 5,  14, 13, 10, 15, 17, 4,  12, 14};
     }
-    else if (entity_types[3][k] == mesh::CellType::pyramid)
+    else if (entity_types.back()[k] == mesh::CellType::pyramid)
     {
+      spdlog::info("Pyramid subdivision");
       refined_cell_list = {0,  5,  6, 13, 7,  1,  8,  5, 13, 9,  3,  10, 8,
                            13, 12, 2, 6,  10, 13, 11, 7, 9,  11, 12, 4};
       if (ktet == -1)
         throw std::runtime_error("Cannot refine mesh with pyramids and no "
                                  "tetrahedra.");
     }
+    else if (entity_types.back()[k] == mesh::CellType::triangle)
+    {
+      spdlog::info("Triangle subdivision");
+      refined_cell_list = {0, 4, 5, 1, 5, 3, 2, 3, 4, 3, 4, 5};
+    }
+    else if (entity_types.back()[k] == mesh::CellType::quadrilateral)
+    {
+      spdlog::info("Quad subdivision");
+      refined_cell_list = {0, 4, 5, 8, 1, 6, 4, 8, 2, 7, 5, 8, 3, 7, 6, 8};
+    }
 
-    auto c_to_v = topology->connectivity({3, k}, {0, 0});
-    auto c_to_e = topology->connectivity({3, k}, {1, 0});
+    auto c_to_v = topology->connectivity({tdim, k}, {0, 0});
+    auto c_to_e = topology->connectivity({tdim, k}, {1, 0});
 
-    for (int c = 0; c < topology->index_maps(3)[k]->size_local(); ++c)
+    for (int c = 0; c < topology->index_maps(tdim)[k]->size_local(); ++c)
     {
       std::vector<std::int64_t> entities;
       for (std::int32_t i : c_to_v->links(c))
@@ -185,22 +212,29 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
         entities.push_back(new_v[1][i]);
       if (e_index.size() > 2)
       {
-        auto conn = topology->connectivity({3, k}, {2, e_index[2]});
-        if (conn)
+        if (tdim == 3)
         {
-          spdlog::debug("Get connectivity from (3,{})->(2,{})", k, e_index[2]);
-          for (std::int32_t i :
-               topology->connectivity({3, k}, {2, e_index[2]})->links(c))
-            entities.push_back(new_v[2][i]);
+          auto conn = topology->connectivity({3, k}, {2, e_index[2]});
+          if (conn)
+          {
+            spdlog::debug("Get connectivity from (3,{})->(2,{})", k,
+                          e_index[2]);
+            for (std::int32_t i :
+                 topology->connectivity({3, k}, {2, e_index[2]})->links(c))
+              entities.push_back(new_v[2][i]);
+          }
         }
+        else
+          entities.push_back(new_v[2][c]);
       }
+
       if (e_index.size() > 3)
         entities.push_back(new_v[3][c]);
 
       for (int i : refined_cell_list)
         mixed_topology[k].push_back(entities[i]);
 
-      if (entity_types[3][k] == mesh::CellType::pyramid)
+      if (entity_types.back()[k] == mesh::CellType::pyramid)
       {
         for (int i : pyr_to_tet_list)
           mixed_topology[ktet].push_back(entities[i]);
diff --git a/python/test/unit/refinement/test_uniform.py b/python/test/unit/refinement/test_uniform.py
index 0397fdda273..1ce55269e0a 100644
--- a/python/test/unit/refinement/test_uniform.py
+++ b/python/test/unit/refinement/test_uniform.py
@@ -3,11 +3,11 @@
 import pytest
 
 import dolfinx
-from dolfinx.mesh import CellType, Mesh, create_unit_cube
+from dolfinx.mesh import CellType, Mesh, create_unit_cube, create_unit_square
 
 
 @pytest.mark.parametrize("ctype", [CellType.hexahedron, CellType.tetrahedron, CellType.prism])
-def test_uniform_refinement(ctype):
+def test_uniform_refinement_3d(ctype):
     mesh = create_unit_cube(MPI.COMM_WORLD, 3, 2, 1, cell_type=ctype)
 
     ncells0 = mesh.topology.index_map(3).size_local
@@ -17,4 +17,17 @@ def test_uniform_refinement(ctype):
     m2 = Mesh(dolfinx.cpp.refinement.uniform_refine(mesh._cpp_object), None)
     ncells1 = m2.topology.index_map(3).size_local
 
-    assert ncells0 * 8 == ncells1
+    assert mesh.comm.allreduce(ncells0) * 8 == mesh.comm.allreduce(ncells1)
+
+
+@pytest.mark.parametrize("ctype", [CellType.triangle, CellType.quadrilateral])
+def test_uniform_refinement_2d(ctype):
+    mesh = create_unit_square(MPI.COMM_WORLD, 12, 11, cell_type=ctype)
+
+    ncells0 = mesh.topology.index_map(2).size_local
+    mesh.topology.create_entities(1)
+
+    m2 = Mesh(dolfinx.cpp.refinement.uniform_refine(mesh._cpp_object), None)
+    ncells1 = m2.topology.index_map(2).size_local
+
+    assert mesh.comm.allreduce(ncells0) * 4 == mesh.comm.allreduce(ncells1)

From b985e1718111252a837585c495ca9c01c6d49425 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Thu, 19 Jun 2025 16:19:41 +0200
Subject: [PATCH 15/28] Fix bug

---
 cpp/dolfinx/mesh/utils.cpp                  |  1 +
 cpp/dolfinx/refinement/uniform.cpp          | 49 +++++++++++----------
 python/test/unit/refinement/test_uniform.py | 25 +++++++++++
 3 files changed, 52 insertions(+), 23 deletions(-)

diff --git a/cpp/dolfinx/mesh/utils.cpp b/cpp/dolfinx/mesh/utils.cpp
index 3491e367120..3054aca08a9 100644
--- a/cpp/dolfinx/mesh/utils.cpp
+++ b/cpp/dolfinx/mesh/utils.cpp
@@ -107,6 +107,7 @@ mesh::create_cell_partitioner(mesh::GhostMode ghost_mode,
   {
     spdlog::info("Compute partition of cells across ranks");
 
+
     // Compute distributed dual graph (for the cells on this process)
     const graph::AdjacencyList dual_graph
         = build_dual_graph(comm, cell_types, cells);
diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
index 4765d2effec..781b6cec45e 100644
--- a/cpp/dolfinx/refinement/uniform.cpp
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -134,22 +134,22 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
   }
 
   // Create new topology...
-
+  const std::vector<mesh::CellType>& cell_entity_types = entity_types.back();
   std::vector<std::vector<std::int64_t>> mixed_topology(
-      entity_types.back().size());
+      cell_entity_types.size());
 
   // Find index of tets in topology list, if any
   int ktet = -1;
-  auto it = std::find(entity_types.back().begin(), entity_types.back().end(),
+  auto it = std::find(cell_entity_types.begin(), cell_entity_types.end(),
                       mesh::CellType::tetrahedron);
-  if (it != entity_types.back().end())
-    ktet = std::distance(entity_types.back().begin(), it);
+  if (it != cell_entity_types.end())
+    ktet = std::distance(cell_entity_types.begin(), it);
   // Topology for tetrahedra which arise from pyramid subdivision
   std::vector<int> pyr_to_tet_list
       = {5, 13, 7, 9, 6, 13, 11, 7, 10, 13, 12, 11, 8, 13, 9, 12};
 
   std::vector<int> refined_cell_list;
-  for (int k = 0; k < static_cast<int>(entity_types.back().size()); ++k)
+  for (int k = 0; k < static_cast<int>(cell_entity_types.size()); ++k)
   {
     // Reserve an estimate of space for the topology of each type
     mixed_topology[k].reserve(mesh.topology()->index_maps(tdim)[k]->size_local()
@@ -157,52 +157,54 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
 
     // Select correct subdivision for celltype
     // Hex -> 8 hex, Prism -> 8 prism, Tet -> 8 tet, Pyr -> 5 pyr + 4 tet
-    if (entity_types.back()[k] == mesh::CellType::hexahedron)
+    if (cell_entity_types[k] == mesh::CellType::hexahedron)
     {
-      spdlog::info("Hex subdivision");
+      spdlog::info("Hex subdivision [{}]", k);
       refined_cell_list
           = {0, 9,  8,  20, 10, 22, 21, 26, 1, 11, 8,  20, 12, 23, 21, 26,
              2, 13, 9,  20, 14, 24, 22, 26, 3, 13, 11, 20, 15, 24, 23, 26,
              4, 16, 10, 21, 17, 25, 22, 26, 5, 16, 12, 21, 18, 25, 23, 26,
              6, 17, 14, 22, 19, 25, 24, 26, 7, 18, 15, 23, 19, 25, 24, 26};
     }
-    else if (entity_types.back()[k] == mesh::CellType::tetrahedron)
+    else if (cell_entity_types[k] == mesh::CellType::tetrahedron)
     {
-      spdlog::info("Tet subdivision");
+      spdlog::info("Tet subdivision [{}]", k);
       refined_cell_list = {0, 7, 8, 9, 1, 5, 6, 9, 2, 4, 6, 8, 3, 4, 5, 7,
                            9, 4, 6, 8, 9, 4, 8, 7, 9, 4, 7, 5, 9, 4, 5, 6};
     }
-    else if (entity_types.back()[k] == mesh::CellType::prism)
+    else if (cell_entity_types[k] == mesh::CellType::prism)
     {
-      spdlog::info("Prism subdivision");
+      spdlog::info("Prism subdivision [{}]", k);
       refined_cell_list
           = {0,  6,  7,  8,  15, 16, 6,  1,  9,  15, 10, 17, 7,  9,  2,  16,
              17, 11, 6,  9,  7,  15, 17, 16, 15, 17, 16, 12, 14, 13, 8,  15,
              16, 3,  12, 13, 11, 17, 16, 5,  14, 13, 10, 15, 17, 4,  12, 14};
     }
-    else if (entity_types.back()[k] == mesh::CellType::pyramid)
+    else if (cell_entity_types[k] == mesh::CellType::pyramid)
     {
-      spdlog::info("Pyramid subdivision");
+      spdlog::info("Pyramid subdivision [{}]", k);
       refined_cell_list = {0,  5,  6, 13, 7,  1,  8,  5, 13, 9,  3,  10, 8,
                            13, 12, 2, 6,  10, 13, 11, 7, 9,  11, 12, 4};
       if (ktet == -1)
         throw std::runtime_error("Cannot refine mesh with pyramids and no "
                                  "tetrahedra.");
     }
-    else if (entity_types.back()[k] == mesh::CellType::triangle)
+    else if (cell_entity_types[k] == mesh::CellType::triangle)
     {
-      spdlog::info("Triangle subdivision");
+      spdlog::info("Triangle subdivision [{}]", k);
       refined_cell_list = {0, 4, 5, 1, 5, 3, 2, 3, 4, 3, 4, 5};
     }
-    else if (entity_types.back()[k] == mesh::CellType::quadrilateral)
+    else if (cell_entity_types[k] == mesh::CellType::quadrilateral)
     {
-      spdlog::info("Quad subdivision");
+      spdlog::info("Quad subdivision [{}]", k);
       refined_cell_list = {0, 4, 5, 8, 1, 6, 4, 8, 2, 7, 5, 8, 3, 7, 6, 8};
     }
 
     auto c_to_v = topology->connectivity({tdim, k}, {0, 0});
     auto c_to_e = topology->connectivity({tdim, k}, {1, 0});
 
+    spdlog::debug("On {}, over {} cells", k,
+                  topology->index_maps(tdim)[k]->size_local());
     for (int c = 0; c < topology->index_maps(tdim)[k]->size_local(); ++c)
     {
       std::vector<std::int64_t> entities;
@@ -217,24 +219,23 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
           auto conn = topology->connectivity({3, k}, {2, e_index[2]});
           if (conn)
           {
-            spdlog::debug("Get connectivity from (3,{})->(2,{})", k,
-                          e_index[2]);
             for (std::int32_t i :
                  topology->connectivity({3, k}, {2, e_index[2]})->links(c))
               entities.push_back(new_v[2][i]);
           }
         }
-        else
+        else if (cell_entity_types[k] == mesh::CellType::quadrilateral)
           entities.push_back(new_v[2][c]);
       }
 
-      if (e_index.size() > 3)
+      if (e_index.size() > 3
+          and cell_entity_types[k] == mesh::CellType::hexahedron)
         entities.push_back(new_v[3][c]);
 
       for (int i : refined_cell_list)
         mixed_topology[k].push_back(entities[i]);
 
-      if (entity_types.back()[k] == mesh::CellType::pyramid)
+      if (cell_entity_types[k] == mesh::CellType::pyramid)
       {
         for (int i : pyr_to_tet_list)
           mixed_topology[ktet].push_back(entities[i]);
@@ -242,11 +243,13 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
     }
   }
 
+  spdlog::debug("Create partitioner");
   auto partitioner = mesh::create_cell_partitioner(mesh::GhostMode::none);
 
   std::vector<std::span<const std::int64_t>> topo_span(mixed_topology.begin(),
                                                        mixed_topology.end());
 
+  spdlog::debug("Create new mesh");
   mesh::Mesh new_mesh = mesh::create_mesh(
       mesh.comm(), mesh.comm(), topo_span, mesh.geometry().cmaps(), mesh.comm(),
       new_x, {new_x.size() / 3, 3}, partitioner);
diff --git a/python/test/unit/refinement/test_uniform.py b/python/test/unit/refinement/test_uniform.py
index 1ce55269e0a..a9d49adf5f4 100644
--- a/python/test/unit/refinement/test_uniform.py
+++ b/python/test/unit/refinement/test_uniform.py
@@ -31,3 +31,28 @@ def test_uniform_refinement_2d(ctype):
     ncells1 = m2.topology.index_map(2).size_local
 
     assert mesh.comm.allreduce(ncells0) * 4 == mesh.comm.allreduce(ncells1)
+
+
+def test_uniform_refine_mixed_mesh(mixed_topology_mesh):
+    mesh = Mesh(mixed_topology_mesh, None)
+
+    ct = mesh.topology.entity_types[3]
+    ncells0 = {ct[j]: mesh.topology.index_maps(3)[j].size_local for j in range(4)}
+    print(ncells0)
+    mesh.topology.create_entities(1)
+    mesh.topology.create_entities(2)
+
+    m2 = Mesh(dolfinx.cpp.refinement.uniform_refine(mesh._cpp_object), None)
+    ncells1 = {ct[j]: m2.topology.index_maps(3)[j].size_local for j in range(4)}
+
+    comm = mesh.comm
+    assert comm.allreduce(ncells0[CellType.hexahedron]) * 8 == comm.allreduce(
+        ncells1[CellType.hexahedron]
+    )
+    assert comm.allreduce(ncells0[CellType.prism]) * 8 == comm.allreduce(ncells1[CellType.prism])
+    assert comm.allreduce(ncells0[CellType.pyramid]) * 5 == comm.allreduce(
+        ncells1[CellType.pyramid]
+    )
+    assert comm.allreduce(ncells0[CellType.tetrahedron]) * 8 + comm.allreduce(
+        ncells0[CellType.pyramid]
+    ) * 4 == comm.allreduce(ncells1[CellType.tetrahedron])

From 398c1bbea727e2e093261f8eeab7c58a61a4cc4a Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Fri, 20 Jun 2025 08:44:22 +0200
Subject: [PATCH 16/28] Update docstring

---
 cpp/dolfinx/refinement/uniform.h |  5 +++--
 python/CMakeLists.txt            |  2 +-
 python/dolfinx/fem/function.py   | 11 +----------
 3 files changed, 5 insertions(+), 13 deletions(-)

diff --git a/cpp/dolfinx/refinement/uniform.h b/cpp/dolfinx/refinement/uniform.h
index 4b286bfbf46..c5b9b869235 100644
--- a/cpp/dolfinx/refinement/uniform.h
+++ b/cpp/dolfinx/refinement/uniform.h
@@ -5,9 +5,10 @@
 namespace dolfinx::refinement
 {
 
-/// @brief Uniform refinement of a 3D mesh containing hex, tet, prism or
-/// pyramids.
+/// @brief Uniform refinement of a 2D or 3D mesh, containing any supported cell
+/// types.
 /// @param mesh Input mesh
+/// @returns Uniformly refined mesh
 mesh::Mesh<double> uniform_refine(const mesh::Mesh<double>& mesh);
 
 } // namespace dolfinx::refinement
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index 68e45482cab..42a5b6ef045 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.21)
 
-project(dolfinx_nanobind LANGUAGES CXX C)
+project(dolfinx_nanobind LANGUAGES CXX)
 
 if(WIN32)
   # Windows requires all symbols to be manually exported. This flag exports all
diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index 4ff1dd9f225..48785e6ed0d 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -37,10 +37,7 @@ class Constant(ufl.Constant):
     ]
 
     def __init__(
-        self,
-        domain,
-        c: typing.Union[np.ndarray, typing.Sequence, np.floating, np.complexfloating],
-        name: typing.Optional[str] = None,
+        self, domain, c: typing.Union[np.ndarray, typing.Sequence, np.floating, np.complexfloating]
     ):
         """A constant with respect to a domain.
 
@@ -64,12 +61,6 @@ def __init__(
         except AttributeError:
             raise AttributeError("Constant value must have a dtype attribute.")
 
-        # Set name
-        if name is None:
-            self.name = "c"
-        else:
-            self.name = name
-
     @property
     def value(self):
         """The value of the constant"""

From 4e8c7f464daf1e671acc9ee5d4092e58425377ce Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Fri, 20 Jun 2025 09:02:05 +0200
Subject: [PATCH 17/28] Templating

---
 cpp/dolfinx/refinement/uniform.cpp     | 16 ++++++++++++----
 cpp/dolfinx/refinement/uniform.h       |  3 ++-
 python/CMakeLists.txt                  |  2 +-
 python/dolfinx/wrappers/refinement.cpp |  2 +-
 4 files changed, 16 insertions(+), 7 deletions(-)

diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
index 781b6cec45e..82303809651 100644
--- a/cpp/dolfinx/refinement/uniform.cpp
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -10,7 +10,8 @@
 
 using namespace dolfinx;
 
-mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
+template <typename T>
+mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
 {
   // Requires edges (and facets for some 3D meshes) to be built already
   auto topology = mesh.topology();
@@ -91,7 +92,7 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
   }
 
   // Copy existing vertices
-  std::vector<double> new_x(nlocal * 3);
+  std::vector<T> new_x(nlocal * 3);
   auto x_g = mesh.geometry().x();
   for (int i = 0; i < index_maps[0]->size_local(); ++i)
     for (int j = 0; j < 3; ++j)
@@ -116,14 +117,14 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
       {
         auto vt = e_to_v->links(w);
         std::size_t nv_ent = vt.size();
-        std::array<double, 3> v = {0, 0, 0};
+        std::array<T, 3> v = {0, 0, 0};
         for (std::size_t i = 0; i < nv_ent; ++i)
         {
           for (int k = 0; k < 3; ++k)
             v[k] += x_g[3 * vertex_to_x[vt[i]] + k];
         }
         for (int k = 0; k < 3; ++k)
-          new_x[(w + w_off) * 3 + k] = v[k] / static_cast<double>(nv_ent);
+          new_x[(w + w_off) * 3 + k] = v[k] / static_cast<T>(nv_ent);
       }
     }
     common::Scatterer sc(*index_maps[j], 1);
@@ -256,3 +257,10 @@ mesh::Mesh<double> refinement::uniform_refine(const mesh::Mesh<double>& mesh)
 
   return new_mesh;
 }
+
+/// @cond Explicit instatiation for float and double
+template mesh::Mesh<double>
+refinement::uniform_refine(const mesh::Mesh<double>& mesh);
+template mesh::Mesh<float>
+refinement::uniform_refine(const mesh::Mesh<float>& mesh);
+/// @endcond
diff --git a/cpp/dolfinx/refinement/uniform.h b/cpp/dolfinx/refinement/uniform.h
index c5b9b869235..9acd93d4f8a 100644
--- a/cpp/dolfinx/refinement/uniform.h
+++ b/cpp/dolfinx/refinement/uniform.h
@@ -9,6 +9,7 @@ namespace dolfinx::refinement
 /// types.
 /// @param mesh Input mesh
 /// @returns Uniformly refined mesh
-mesh::Mesh<double> uniform_refine(const mesh::Mesh<double>& mesh);
+template <typename T>
+mesh::Mesh<T> uniform_refine(const mesh::Mesh<T>& mesh);
 
 } // namespace dolfinx::refinement
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index 42a5b6ef045..68e45482cab 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.21)
 
-project(dolfinx_nanobind LANGUAGES CXX)
+project(dolfinx_nanobind LANGUAGES CXX C)
 
 if(WIN32)
   # Windows requires all symbols to be manually exported. This flag exports all
diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp
index 10b8860eb83..eb6f99e724a 100644
--- a/python/dolfinx/wrappers/refinement.cpp
+++ b/python/dolfinx/wrappers/refinement.cpp
@@ -34,7 +34,7 @@ namespace
 template <std::floating_point T>
 void export_refinement(nb::module_& m)
 {
-  m.def("uniform_refine", dolfinx::refinement::uniform_refine);
+  m.def("uniform_refine", &dolfinx::refinement::uniform_refine<T>);
 
   m.def(
       "refine",

From 92f6155a01c928777af723f0f2d827e1f1a9322d Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Fri, 20 Jun 2025 09:06:31 +0200
Subject: [PATCH 18/28] tparam

---
 cpp/dolfinx/refinement/uniform.h | 1 +
 1 file changed, 1 insertion(+)

diff --git a/cpp/dolfinx/refinement/uniform.h b/cpp/dolfinx/refinement/uniform.h
index 9acd93d4f8a..01c4f268bce 100644
--- a/cpp/dolfinx/refinement/uniform.h
+++ b/cpp/dolfinx/refinement/uniform.h
@@ -7,6 +7,7 @@ namespace dolfinx::refinement
 
 /// @brief Uniform refinement of a 2D or 3D mesh, containing any supported cell
 /// types.
+/// @tparam T Scalar type of the mesh geometry
 /// @param mesh Input mesh
 /// @returns Uniformly refined mesh
 template <typename T>

From da8b150ab09ad645440390f85e0970252b67e070 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Fri, 20 Jun 2025 11:42:35 +0200
Subject: [PATCH 19/28] More description

---
 cpp/dolfinx/refinement/uniform.h | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/cpp/dolfinx/refinement/uniform.h b/cpp/dolfinx/refinement/uniform.h
index 01c4f268bce..835bb83ae15 100644
--- a/cpp/dolfinx/refinement/uniform.h
+++ b/cpp/dolfinx/refinement/uniform.h
@@ -7,6 +7,10 @@ namespace dolfinx::refinement
 
 /// @brief Uniform refinement of a 2D or 3D mesh, containing any supported cell
 /// types.
+/// Hexahedral, tetrahedral and prism cells are subdivided into 8, each being
+/// similar to the original cell. Pyramid cells are subdivided into 5 similar
+/// pyramids, plus 4 tetrahedra. Triangle and quadrilateral cells are subdivided
+/// into 4 similar subcells.
 /// @tparam T Scalar type of the mesh geometry
 /// @param mesh Input mesh
 /// @returns Uniformly refined mesh

From f200ee48978ef0797a0bd899fc5ecf63b3f83d35 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Fri, 20 Jun 2025 12:08:41 +0200
Subject: [PATCH 20/28] Remove blank line

---
 cpp/dolfinx/mesh/utils.cpp | 1 -
 1 file changed, 1 deletion(-)

diff --git a/cpp/dolfinx/mesh/utils.cpp b/cpp/dolfinx/mesh/utils.cpp
index 3054aca08a9..3491e367120 100644
--- a/cpp/dolfinx/mesh/utils.cpp
+++ b/cpp/dolfinx/mesh/utils.cpp
@@ -107,7 +107,6 @@ mesh::create_cell_partitioner(mesh::GhostMode ghost_mode,
   {
     spdlog::info("Compute partition of cells across ranks");
 
-
     // Compute distributed dual graph (for the cells on this process)
     const graph::AdjacencyList dual_graph
         = build_dual_graph(comm, cell_types, cells);

From 4dbce363081e35da77c1a232eb72abc53319e57e Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Sun, 22 Jun 2025 15:18:15 +0100
Subject: [PATCH 21/28] Remove commented line

---
 cpp/dolfinx/refinement/uniform.cpp | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
index 82303809651..5c729a3ba57 100644
--- a/cpp/dolfinx/refinement/uniform.cpp
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -54,7 +54,7 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
   {
     if (topology->index_maps(dim).empty())
       throw std::runtime_error(
-          "Missing topology for dimension " + std::to_string(dim)
+          "Missing entities of dimension " + std::to_string(dim)
           + ", need to call create_entities(" + std::to_string(dim) + ")");
     index_maps.push_back(topology->index_maps(dim)[e_index[dim]]);
     new_v.push_back(std::vector<std::int64_t>(
@@ -62,8 +62,6 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
     nlocal += index_maps.back()->size_local();
   }
 
-  // std::int64_t nv_global = index_maps[0]->size_global();
-
   // Get current geometry and put into new array for vertices
   std::vector<std::int32_t> vertex_to_x(index_maps[0]->size_local()
                                         + index_maps[0]->num_ghosts());

From 4fa4a9c7ce2a2e0bbe8fa1827b739684ff53c2f5 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Mon, 23 Jun 2025 10:57:40 +0100
Subject: [PATCH 22/28] Suggestions from review

---
 cpp/dolfinx/refinement/uniform.cpp | 103 +++++++++++++++++------------
 1 file changed, 62 insertions(+), 41 deletions(-)

diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
index 5c729a3ba57..8585388d781 100644
--- a/cpp/dolfinx/refinement/uniform.cpp
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -25,8 +25,11 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
   std::vector<std::vector<mesh::CellType>> entity_types;
   for (int i = 0; i < tdim + 1; ++i)
     entity_types.push_back(topology->entity_types(i));
+  const std::vector<mesh::CellType>& cell_entity_types = entity_types.back();
 
   // Get index maps for vertices and edges and create maps for new vertices
+  // New vertices are created for each old vertex, and also at the centre of
+  // each edge, each quadrilateral facet, and each hexahedral cell.
   std::vector<std::vector<std::int64_t>> new_v;
   std::vector<std::shared_ptr<const common::IndexMap>> index_maps;
 
@@ -42,10 +45,10 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
   if (tdim == 3)
   {
     // In 3D, check for hexahedral cells, and get index, if any.
-    if (auto it = std::find(entity_types[3].begin(), entity_types[3].end(),
+    if (auto it = std::find(cell_entity_types.begin(), cell_entity_types.end(),
                             mesh::CellType::hexahedron);
         it != entity_types[3].end())
-      e_index.push_back(std::distance(entity_types[3].begin(), it));
+      e_index.push_back(std::distance(cell_entity_types.begin(), it));
   }
 
   // Add up all local vertices, edges, quad facets and hex cells.
@@ -66,8 +69,8 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
   std::vector<std::int32_t> vertex_to_x(index_maps[0]->size_local()
                                         + index_maps[0]->num_ghosts());
 
-  // Iterate over cells
-  for (int j = 0; j < static_cast<int>(entity_types.back().size()); ++j)
+  // Iterate over cells for each cell type
+  for (int j = 0; j < static_cast<int>(cell_entity_types.size()); ++j)
   {
     // Get geometry for each cell type
     auto x_dofmap = mesh.geometry().dofmap(j);
@@ -75,13 +78,11 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
     auto dof_layout = mesh.geometry().cmaps().at(j).create_dof_layout();
     std::vector<int> entity_dofs(dof_layout.num_dofs());
     for (int k = 0; k < dof_layout.num_dofs(); ++k)
-      entity_dofs[k] = dof_layout.entity_dofs(0, k)[0];
+      entity_dofs[k] = dof_layout.entity_dofs(0, k).front();
 
     // Iterate over cells of this type
-    for (std::int32_t c = 0;
-         c < topology->index_maps(tdim)[j]->size_local()
-                 + topology->index_maps(tdim)[j]->num_ghosts();
-         ++c)
+    const auto im = topology->index_maps(tdim)[j];
+    for (std::int32_t c = 0; c < im->size_local() + im->num_ghosts(); ++c)
     {
       auto vertices = c_to_v->links(c);
       for (std::size_t i = 0; i < vertices.size(); ++i)
@@ -89,6 +90,11 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
     }
   }
 
+  // Compute offset for vertices related to each entity type
+  std::vector<std::int32_t> entity_offsets(index_maps.size() + 1, 0);
+  for (std::size_t i = 0; i < index_maps.size(); ++i)
+    entity_offsets[i + 1] = entity_offsets[i] + index_maps[i]->size_local();
+
   // Copy existing vertices
   std::vector<T> new_x(nlocal * 3);
   auto x_g = mesh.geometry().x();
@@ -96,44 +102,51 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
     for (int j = 0; j < 3; ++j)
       new_x[i * 3 + j] = x_g[vertex_to_x[i] * 3 + j];
 
+  // Create vertices on edges, quad facet and hex cell centres
+  for (int j = 1; j < static_cast<int>(index_maps.size()); ++j)
+  {
+    auto e_to_v = topology->connectivity({j, e_index[j]}, {0, 0});
+    std::int32_t w_off = entity_offsets[j];
+    std::int32_t num_entities = index_maps[j]->size_local();
+    assert(num_entities == entity_offsets[j + 1] - entity_offsets[j]);
+
+    for (std::int32_t w = 0; w < num_entities; ++w)
+    {
+      auto vt = e_to_v->links(w);
+      std::size_t nv_ent = vt.size();
+      std::array<T, 3> v = {0, 0, 0};
+      for (std::size_t i = 0; i < nv_ent; ++i)
+      {
+        for (int k = 0; k < 3; ++k)
+          v[k] += x_g[3 * vertex_to_x[vt[i]] + k];
+      }
+      // Get the (linear) center of the entity.
+      // TODO: change for higher-order geometry to use map.
+      for (int k = 0; k < 3; ++k)
+        new_x[(w + w_off) * 3 + k] = v[k] / static_cast<T>(nv_ent);
+    }
+  }
+
   // Get this process global offset and range
   std::int64_t nscan;
   MPI_Scan(&nlocal, &nscan, 1, MPI_INT64_T, MPI_SUM, mesh.comm());
   std::array<std::int64_t, 2> local_range = {nscan - nlocal, nscan};
 
-  std::int32_t w_off = 0;
-  for (int j = 0; j < static_cast<int>(index_maps.size()); ++j)
+  // Compute indices for new vertices and share across processes
+  for (std::size_t j = 0; j < index_maps.size(); ++j)
   {
     std::int32_t num_entities = index_maps[j]->size_local();
+    assert(num_entities == entity_offsets[j + 1] - entity_offsets[j]);
     std::iota(new_v[j].begin(), std::next(new_v[j].begin(), num_entities),
-              local_range[0] + w_off);
+              local_range[0] + entity_offsets[j]);
 
-    if (j > 0)
-    {
-      auto e_to_v = topology->connectivity({j, e_index[j]}, {0, 0});
-      for (std::int32_t w = 0; w < num_entities; ++w)
-      {
-        auto vt = e_to_v->links(w);
-        std::size_t nv_ent = vt.size();
-        std::array<T, 3> v = {0, 0, 0};
-        for (std::size_t i = 0; i < nv_ent; ++i)
-        {
-          for (int k = 0; k < 3; ++k)
-            v[k] += x_g[3 * vertex_to_x[vt[i]] + k];
-        }
-        for (int k = 0; k < 3; ++k)
-          new_x[(w + w_off) * 3 + k] = v[k] / static_cast<T>(nv_ent);
-      }
-    }
     common::Scatterer sc(*index_maps[j], 1);
     sc.scatter_fwd(std::span<const std::int64_t>(new_v[j]),
                    std::span(std::next(new_v[j].begin(), num_entities),
                              index_maps[j]->num_ghosts()));
-    w_off += num_entities;
   }
 
-  // Create new topology...
-  const std::vector<mesh::CellType>& cell_entity_types = entity_types.back();
+  // Create new topology
   std::vector<std::vector<std::int64_t>> mixed_topology(
       cell_entity_types.size());
 
@@ -144,7 +157,7 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
   if (it != cell_entity_types.end())
     ktet = std::distance(cell_entity_types.begin(), it);
   // Topology for tetrahedra which arise from pyramid subdivision
-  std::vector<int> pyr_to_tet_list
+  std::array<int, 16> pyr_to_tet_list
       = {5, 13, 7, 9, 6, 13, 11, 7, 10, 13, 12, 11, 8, 13, 9, 12};
 
   std::vector<int> refined_cell_list;
@@ -156,9 +169,13 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
 
     // Select correct subdivision for celltype
     // Hex -> 8 hex, Prism -> 8 prism, Tet -> 8 tet, Pyr -> 5 pyr + 4 tet
+    // Each `refined_cell_list` is the topology for cells of the same type,
+    // flattened, hence 64 entries for hex, 32 for tet, 48 for prism, 25 for
+    // pyramid. Additionally for pyramid there are 16 tetrahedron entries (see
+    // above).
     if (cell_entity_types[k] == mesh::CellType::hexahedron)
     {
-      spdlog::info("Hex subdivision [{}]", k);
+      spdlog::debug("Hex subdivision [{}]", k);
       refined_cell_list
           = {0, 9,  8,  20, 10, 22, 21, 26, 1, 11, 8,  20, 12, 23, 21, 26,
              2, 13, 9,  20, 14, 24, 22, 26, 3, 13, 11, 20, 15, 24, 23, 26,
@@ -167,13 +184,13 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
     }
     else if (cell_entity_types[k] == mesh::CellType::tetrahedron)
     {
-      spdlog::info("Tet subdivision [{}]", k);
+      spdlog::debug("Tet subdivision [{}]", k);
       refined_cell_list = {0, 7, 8, 9, 1, 5, 6, 9, 2, 4, 6, 8, 3, 4, 5, 7,
                            9, 4, 6, 8, 9, 4, 8, 7, 9, 4, 7, 5, 9, 4, 5, 6};
     }
     else if (cell_entity_types[k] == mesh::CellType::prism)
     {
-      spdlog::info("Prism subdivision [{}]", k);
+      spdlog::debug("Prism subdivision [{}]", k);
       refined_cell_list
           = {0,  6,  7,  8,  15, 16, 6,  1,  9,  15, 10, 17, 7,  9,  2,  16,
              17, 11, 6,  9,  7,  15, 17, 16, 15, 17, 16, 12, 14, 13, 8,  15,
@@ -181,7 +198,7 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
     }
     else if (cell_entity_types[k] == mesh::CellType::pyramid)
     {
-      spdlog::info("Pyramid subdivision [{}]", k);
+      spdlog::debug("Pyramid subdivision [{}]", k);
       refined_cell_list = {0,  5,  6, 13, 7,  1,  8,  5, 13, 9,  3,  10, 8,
                            13, 12, 2, 6,  10, 13, 11, 7, 9,  11, 12, 4};
       if (ktet == -1)
@@ -190,31 +207,33 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
     }
     else if (cell_entity_types[k] == mesh::CellType::triangle)
     {
-      spdlog::info("Triangle subdivision [{}]", k);
+      spdlog::debug("Triangle subdivision [{}]", k);
       refined_cell_list = {0, 4, 5, 1, 5, 3, 2, 3, 4, 3, 4, 5};
     }
     else if (cell_entity_types[k] == mesh::CellType::quadrilateral)
     {
-      spdlog::info("Quad subdivision [{}]", k);
+      spdlog::debug("Quad subdivision [{}]", k);
       refined_cell_list = {0, 4, 5, 8, 1, 6, 4, 8, 2, 7, 5, 8, 3, 7, 6, 8};
     }
 
     auto c_to_v = topology->connectivity({tdim, k}, {0, 0});
     auto c_to_e = topology->connectivity({tdim, k}, {1, 0});
 
-    spdlog::debug("On {}, over {} cells", k,
-                  topology->index_maps(tdim)[k]->size_local());
     for (int c = 0; c < topology->index_maps(tdim)[k]->size_local(); ++c)
     {
+      // Cell topology defined through its globally numbered vertices
       std::vector<std::int64_t> entities;
+      // Extract new global vertex number for existing vertices
       for (std::int32_t i : c_to_v->links(c))
         entities.push_back(new_v[0][i]);
+      // Indices for vertices inserted on edges
       for (std::int32_t i : c_to_e->links(c))
         entities.push_back(new_v[1][i]);
       if (e_index.size() > 2)
       {
         if (tdim == 3)
         {
+          // Add any vertices on quadrilateral facets (3D mesh)
           auto conn = topology->connectivity({3, k}, {2, e_index[2]});
           if (conn)
           {
@@ -223,10 +242,12 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
               entities.push_back(new_v[2][i]);
           }
         }
+        // Add vertices for quadrilateral cells (2D mesh)
         else if (cell_entity_types[k] == mesh::CellType::quadrilateral)
           entities.push_back(new_v[2][c]);
       }
 
+      // Add vertices for hex cell centres
       if (e_index.size() > 3
           and cell_entity_types[k] == mesh::CellType::hexahedron)
         entities.push_back(new_v[3][c]);

From 4a67cc6b83abffd15e7444d94d29380caf86cbb0 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Mon, 23 Jun 2025 11:31:18 +0100
Subject: [PATCH 23/28] Try fixing pyamg

---
 .github/workflows/ccpp.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.github/workflows/ccpp.yml b/.github/workflows/ccpp.yml
index 82657df5479..89d2fd755a3 100644
--- a/.github/workflows/ccpp.yml
+++ b/.github/workflows/ccpp.yml
@@ -147,6 +147,7 @@ jobs:
       - name: Run demos (Python, serial)
         run: |
           pip install pytest-xdist
+          pip install --upgrade pyamg
           python -m pytest -n auto -m serial --durations=10 python/demo/test.py
       - name: Run demos (Python, MPI (np=3))
         run: python -m pytest -m mpi --num-proc=3 python/demo/test.py

From bb0323da7bb4d3f31ef09ef8e5f601f68e9c9884 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Mon, 23 Jun 2025 13:41:50 +0100
Subject: [PATCH 24/28] Pass in partitioner

---
 cpp/dolfinx/refinement/uniform.cpp | 16 ++++++++--------
 cpp/dolfinx/refinement/uniform.h   |  5 ++++-
 2 files changed, 12 insertions(+), 9 deletions(-)

diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
index 8585388d781..30c4223a13c 100644
--- a/cpp/dolfinx/refinement/uniform.cpp
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -11,7 +11,9 @@
 using namespace dolfinx;
 
 template <typename T>
-mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
+mesh::Mesh<T>
+refinement::uniform_refine(const mesh::Mesh<T>& mesh,
+                           const mesh::CellPartitionFunction& partitioner)
 {
   // Requires edges (and facets for some 3D meshes) to be built already
   auto topology = mesh.topology();
@@ -263,13 +265,9 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
     }
   }
 
-  spdlog::debug("Create partitioner");
-  auto partitioner = mesh::create_cell_partitioner(mesh::GhostMode::none);
-
+  spdlog::debug("Create new mesh");
   std::vector<std::span<const std::int64_t>> topo_span(mixed_topology.begin(),
                                                        mixed_topology.end());
-
-  spdlog::debug("Create new mesh");
   mesh::Mesh new_mesh = mesh::create_mesh(
       mesh.comm(), mesh.comm(), topo_span, mesh.geometry().cmaps(), mesh.comm(),
       new_x, {new_x.size() / 3, 3}, partitioner);
@@ -279,7 +277,9 @@ mesh::Mesh<T> refinement::uniform_refine(const mesh::Mesh<T>& mesh)
 
 /// @cond Explicit instatiation for float and double
 template mesh::Mesh<double>
-refinement::uniform_refine(const mesh::Mesh<double>& mesh);
+refinement::uniform_refine(const mesh::Mesh<double>& mesh,
+                           const mesh::CellPartitionFunction&);
 template mesh::Mesh<float>
-refinement::uniform_refine(const mesh::Mesh<float>& mesh);
+refinement::uniform_refine(const mesh::Mesh<float>& mesh,
+                           const mesh::CellPartitionFunction&);
 /// @endcond
diff --git a/cpp/dolfinx/refinement/uniform.h b/cpp/dolfinx/refinement/uniform.h
index 835bb83ae15..a103f9e688d 100644
--- a/cpp/dolfinx/refinement/uniform.h
+++ b/cpp/dolfinx/refinement/uniform.h
@@ -15,6 +15,9 @@ namespace dolfinx::refinement
 /// @param mesh Input mesh
 /// @returns Uniformly refined mesh
 template <typename T>
-mesh::Mesh<T> uniform_refine(const mesh::Mesh<T>& mesh);
+mesh::Mesh<T>
+uniform_refine(const mesh::Mesh<T>& mesh,
+               const mesh::CellPartitionFunction& partitioner
+               = mesh::create_cell_partitioner(mesh::GhostMode::none));
 
 } // namespace dolfinx::refinement

From bbb4186072e30971cba96b62232b34b0290621ac Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Mon, 23 Jun 2025 22:49:58 +0100
Subject: [PATCH 25/28] Update

---
 .github/workflows/ccpp.yml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/.github/workflows/ccpp.yml b/.github/workflows/ccpp.yml
index 89d2fd755a3..82657df5479 100644
--- a/.github/workflows/ccpp.yml
+++ b/.github/workflows/ccpp.yml
@@ -147,7 +147,6 @@ jobs:
       - name: Run demos (Python, serial)
         run: |
           pip install pytest-xdist
-          pip install --upgrade pyamg
           python -m pytest -n auto -m serial --durations=10 python/demo/test.py
       - name: Run demos (Python, MPI (np=3))
         run: python -m pytest -m mpi --num-proc=3 python/demo/test.py

From 2d29f30aed4da9da72394358b875a9aedb5b9c2c Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Mon, 23 Jun 2025 22:54:42 +0100
Subject: [PATCH 26/28] Use switch

---
 cpp/dolfinx/refinement/uniform.cpp | 38 +++++++++++++++++-------------
 1 file changed, 22 insertions(+), 16 deletions(-)

diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp
index 30c4223a13c..ae9cc499ed5 100644
--- a/cpp/dolfinx/refinement/uniform.cpp
+++ b/cpp/dolfinx/refinement/uniform.cpp
@@ -175,47 +175,53 @@ refinement::uniform_refine(const mesh::Mesh<T>& mesh,
     // flattened, hence 64 entries for hex, 32 for tet, 48 for prism, 25 for
     // pyramid. Additionally for pyramid there are 16 tetrahedron entries (see
     // above).
-    if (cell_entity_types[k] == mesh::CellType::hexahedron)
+
+    switch (cell_entity_types[k])
     {
+    case mesh::CellType::hexahedron:
       spdlog::debug("Hex subdivision [{}]", k);
       refined_cell_list
           = {0, 9,  8,  20, 10, 22, 21, 26, 1, 11, 8,  20, 12, 23, 21, 26,
              2, 13, 9,  20, 14, 24, 22, 26, 3, 13, 11, 20, 15, 24, 23, 26,
              4, 16, 10, 21, 17, 25, 22, 26, 5, 16, 12, 21, 18, 25, 23, 26,
              6, 17, 14, 22, 19, 25, 24, 26, 7, 18, 15, 23, 19, 25, 24, 26};
-    }
-    else if (cell_entity_types[k] == mesh::CellType::tetrahedron)
-    {
+      break;
+
+    case mesh::CellType::tetrahedron:
       spdlog::debug("Tet subdivision [{}]", k);
       refined_cell_list = {0, 7, 8, 9, 1, 5, 6, 9, 2, 4, 6, 8, 3, 4, 5, 7,
                            9, 4, 6, 8, 9, 4, 8, 7, 9, 4, 7, 5, 9, 4, 5, 6};
-    }
-    else if (cell_entity_types[k] == mesh::CellType::prism)
-    {
+      break;
+
+    case mesh::CellType::prism:
       spdlog::debug("Prism subdivision [{}]", k);
       refined_cell_list
           = {0,  6,  7,  8,  15, 16, 6,  1,  9,  15, 10, 17, 7,  9,  2,  16,
              17, 11, 6,  9,  7,  15, 17, 16, 15, 17, 16, 12, 14, 13, 8,  15,
              16, 3,  12, 13, 11, 17, 16, 5,  14, 13, 10, 15, 17, 4,  12, 14};
-    }
-    else if (cell_entity_types[k] == mesh::CellType::pyramid)
-    {
+      break;
+
+    case mesh::CellType::pyramid:
       spdlog::debug("Pyramid subdivision [{}]", k);
       refined_cell_list = {0,  5,  6, 13, 7,  1,  8,  5, 13, 9,  3,  10, 8,
                            13, 12, 2, 6,  10, 13, 11, 7, 9,  11, 12, 4};
       if (ktet == -1)
         throw std::runtime_error("Cannot refine mesh with pyramids and no "
                                  "tetrahedra.");
-    }
-    else if (cell_entity_types[k] == mesh::CellType::triangle)
-    {
+      break;
+
+    case mesh::CellType::triangle:
       spdlog::debug("Triangle subdivision [{}]", k);
       refined_cell_list = {0, 4, 5, 1, 5, 3, 2, 3, 4, 3, 4, 5};
-    }
-    else if (cell_entity_types[k] == mesh::CellType::quadrilateral)
-    {
+      break;
+
+    case mesh::CellType::quadrilateral:
       spdlog::debug("Quad subdivision [{}]", k);
       refined_cell_list = {0, 4, 5, 8, 1, 6, 4, 8, 2, 7, 5, 8, 3, 7, 6, 8};
+      break;
+
+    default:
+      throw std::runtime_error("Unhandled cell type");
     }
 
     auto c_to_v = topology->connectivity({tdim, k}, {0, 0});

From 8da67b0ed3401a95e410da3dd7383fc233612500 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Tue, 24 Jun 2025 11:46:25 +0100
Subject: [PATCH 27/28] Try to fix binding

---
 python/dolfinx/wrappers/refinement.cpp | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp
index eb6f99e724a..f79add65eda 100644
--- a/python/dolfinx/wrappers/refinement.cpp
+++ b/python/dolfinx/wrappers/refinement.cpp
@@ -34,7 +34,10 @@ namespace
 template <std::floating_point T>
 void export_refinement(nb::module_& m)
 {
-  m.def("uniform_refine", &dolfinx::refinement::uniform_refine<T>);
+  m.def(
+      "uniform_refine", [](const dolfinx::mesh::Mesh<T>& mesh)
+      { return dolfinx::refinement::uniform_refine<T>(mesh); },
+      nb::arg("mesh"));
 
   m.def(
       "refine",

From 9d70bef55fb7062826d69986675fdedd8c4ed2ba Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Tue, 24 Jun 2025 11:56:02 +0100
Subject: [PATCH 28/28] Doc fix

---
 cpp/dolfinx/refinement/uniform.h | 1 +
 1 file changed, 1 insertion(+)

diff --git a/cpp/dolfinx/refinement/uniform.h b/cpp/dolfinx/refinement/uniform.h
index a103f9e688d..5f70c150e4f 100644
--- a/cpp/dolfinx/refinement/uniform.h
+++ b/cpp/dolfinx/refinement/uniform.h
@@ -13,6 +13,7 @@ namespace dolfinx::refinement
 /// into 4 similar subcells.
 /// @tparam T Scalar type of the mesh geometry
 /// @param mesh Input mesh
+/// @param partitioner Function to partition new mesh across processes.
 /// @returns Uniformly refined mesh
 template <typename T>
 mesh::Mesh<T>