From f82f6677d7352644d198e256b636bf621a1a9c39 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Fri, 20 Mar 2026 15:39:01 +0100
Subject: [PATCH 01/22] Add reproducible mypy make recipe

---
 .github/workflows/ccpp.yml | 10 ++--------
 python/Makefile            | 14 ++++++++++++++
 2 files changed, 16 insertions(+), 8 deletions(-)
 create mode 100644 python/Makefile

diff --git a/.github/workflows/ccpp.yml b/.github/workflows/ccpp.yml
index c009dfd4f0..21fe82f6f9 100644
--- a/.github/workflows/ccpp.yml
+++ b/.github/workflows/ccpp.yml
@@ -132,15 +132,9 @@ jobs:
           pip install --check-build-dependencies --no-build-isolation --config-settings=cmake.build-type="Developer" 'python/[test]'
           python -c "from mpi4py import MPI; import dolfinx; assert not dolfinx.has_petsc; assert not dolfinx.has_petsc4py; assert dolfinx.has_superlu_dist"
 
-      - name: Run mypy  # Use a venv to avoid NumPy upgrades that are incompatible with numba
+      - name: Run mypy
         working-directory: python
-        run: |
-          python -m venv mypy-env
-          . ./mypy-env/bin/activate
-          pip install mypy types-cffi scipy-stubs
-          mypy -p dolfinx
-        # mypy test
-        # mypy demo
+        run: make lint-mypy
 
       - name: Install gmsh and pyvista (and dependencies)
         run: |
diff --git a/python/Makefile b/python/Makefile
new file mode 100644
index 0000000000..4f59ecff41
--- /dev/null
+++ b/python/Makefile
@@ -0,0 +1,14 @@
+# Note: mypy is currenlty on CI tested without petsc4py!
+lint-mypy:
+	python -m venv dolfinx-venv-mypy
+	. dolfinx-venv-mypy/bin/activate && pip install mypy types-cffi scipy-stubs
+	. dolfinx-venv-mypy/bin/activate && mypy -p dolfinx
+# 	. dolfinx-venv-mypy/bin/activate && mypy test
+# 	. dolfinx-venv-mypy/bin/activate && mypy demo
+
+lint-mypy-clean:
+	rm -rf dolfinx-venv-mypy
+
+lint: lint-mypy
+
+clean: lint-mypy-clean

From 452f4ae05d43a74122da83870d06341c8f5c891f Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Fri, 20 Mar 2026 15:46:31 +0100
Subject: [PATCH 02/22] Mypy check test/

---
 python/Makefile | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/Makefile b/python/Makefile
index 4f59ecff41..a8a0c8213f 100644
--- a/python/Makefile
+++ b/python/Makefile
@@ -3,7 +3,7 @@ lint-mypy:
 	python -m venv dolfinx-venv-mypy
 	. dolfinx-venv-mypy/bin/activate && pip install mypy types-cffi scipy-stubs
 	. dolfinx-venv-mypy/bin/activate && mypy -p dolfinx
-# 	. dolfinx-venv-mypy/bin/activate && mypy test
+	. dolfinx-venv-mypy/bin/activate && mypy test
 # 	. dolfinx-venv-mypy/bin/activate && mypy demo
 
 lint-mypy-clean:

From 4cb783a790f19a2aa17729d7b531a900ac585564 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Fri, 20 Mar 2026 16:10:27 +0100
Subject: [PATCH 03/22] Fix: dbc value type

---
 python/dolfinx/fem/bcs.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/python/dolfinx/fem/bcs.py b/python/dolfinx/fem/bcs.py
index 10dce42ba3..d4abcaa8a7 100644
--- a/python/dolfinx/fem/bcs.py
+++ b/python/dolfinx/fem/bcs.py
@@ -12,7 +12,6 @@
 
 from __future__ import annotations
 
-import numbers
 from collections.abc import Callable, Iterable
 
 import numpy as np
@@ -172,7 +171,7 @@ def dof_indices(self) -> tuple[npt.NDArray[np.int32], int]:
 
 
 def dirichletbc(
-    value: Function | Constant | np.ndarray,
+    value: Function | Constant | np.ndarray | float | complex,
     dofs: npt.NDArray[np.int32],
     V: dolfinx.fem.FunctionSpace | None = None,
 ) -> DirichletBC:
@@ -193,7 +192,7 @@ def dirichletbc(
         A representation of the boundary condition for modifying linear
         systems.
     """
-    if isinstance(value, numbers.Number):
+    if isinstance(value, float | complex):
         value = np.asarray(value)
 
     try:

From a75f18d34bee6f37da15f1c10e4ee6bd44d7d15b Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Fri, 20 Mar 2026 15:51:55 +0100
Subject: [PATCH 04/22] Mypy check (and fix) demo/

---
 python/Makefile                                    |  2 +-
 python/demo/demo_pml.py                            |  4 ++--
 python/demo/demo_pyamg.py                          | 10 +++++-----
 python/demo/demo_scattering_boundary_conditions.py |  4 ++--
 4 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/python/Makefile b/python/Makefile
index a8a0c8213f..88332551f5 100644
--- a/python/Makefile
+++ b/python/Makefile
@@ -4,7 +4,7 @@ lint-mypy:
 	. dolfinx-venv-mypy/bin/activate && pip install mypy types-cffi scipy-stubs
 	. dolfinx-venv-mypy/bin/activate && mypy -p dolfinx
 	. dolfinx-venv-mypy/bin/activate && mypy test
-# 	. dolfinx-venv-mypy/bin/activate && mypy demo
+	. dolfinx-venv-mypy/bin/activate && mypy demo
 
 lint-mypy-clean:
 	rm -rf dolfinx-venv-mypy
diff --git a/python/demo/demo_pml.py b/python/demo/demo_pml.py
index 27d4fb9464..c5634bf72d 100644
--- a/python/demo/demo_pml.py
+++ b/python/demo/demo_pml.py
@@ -214,8 +214,8 @@ def generate_mesh_wire(
 
 def compute_a(nu: int, m: complex, alpha: float) -> float:
     """Compute the Mie coefficient a_nu for a cylinder."""
-    J_nu_alpha = jv(nu, alpha)
-    J_nu_malpha = jv(nu, m * alpha)
+    J_nu_alpha = jv(nu, alpha)  # type: ignore
+    J_nu_malpha = jv(nu, m * alpha)  # type: ignore
     J_nu_alpha_p = jvp(nu, alpha, 1)
     J_nu_malpha_p = jvp(nu, m * alpha, 1)
 
diff --git a/python/demo/demo_pyamg.py b/python/demo/demo_pyamg.py
index bf6e4e68c8..6027e028ce 100644
--- a/python/demo/demo_pyamg.py
+++ b/python/demo/demo_pyamg.py
@@ -63,7 +63,7 @@ def poisson_problem(dtype: npt.DTypeLike, solver_type: str) -> None:
         solver_type: pyamg solver type, either "ruge_stuben" or
             "smoothed_aggregation"
     """
-    real_type = np.real(dtype(0)).dtype
+    real_type = np.real(dtype(0)).dtype  # type: ignore
     mesh = create_box(
         comm=MPI.COMM_WORLD,
         points=[(0.0, 0.0, 0.0), (3.0, 2.0, 1.0)],
@@ -84,7 +84,7 @@ def poisson_problem(dtype: npt.DTypeLike, solver_type: str) -> None:
 
     dofs = locate_dofs_topological(V=V, entity_dim=fdim, entities=facets)
 
-    bc = dirichletbc(value=dtype(0), dofs=dofs, V=V)
+    bc = dirichletbc(value=0.0, dofs=dofs, V=V)
 
     u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
     x = ufl.SpatialCoordinate(mesh)
@@ -110,14 +110,14 @@ def poisson_problem(dtype: npt.DTypeLike, solver_type: str) -> None:
     print(ml)
 
     # Solve linear systems
-    print(f"\nSolve Poisson equation: {dtype.__name__}")
+    print(f"\nSolve Poisson equation: {dtype!s}")
     res: list[float] = []
     tol = 1e-10 if real_type == np.float64 else 1e-6
     uh.x.array[:] = ml.solve(b.array, tol=tol, residuals=res, accel="cg")
     for i, q in enumerate(res):
         print(f"Convergence history: iteration {i}, residual= {q}")
 
-    with io.XDMFFile(mesh.comm, f"out_pyamg/poisson_{dtype.__name__}.xdmf", "w") as file:
+    with io.XDMFFile(mesh.comm, f"out_pyamg/poisson_{dtype!s}.xdmf", "w") as file:
         file.write_mesh(mesh)
         file.write_function(uh)
 
@@ -126,7 +126,7 @@ def poisson_problem(dtype: npt.DTypeLike, solver_type: str) -> None:
 
 
 # +
-def nullspace_elasticty(Q: fem.FunctionSpace) -> list[np.ndarray]:
+def nullspace_elasticty(Q: fem.FunctionSpace) -> npt.NDArray:
     """Create the elasticity (near)nulspace.
 
     Args:
diff --git a/python/demo/demo_scattering_boundary_conditions.py b/python/demo/demo_scattering_boundary_conditions.py
index f3bd735f6d..3c40b83574 100644
--- a/python/demo/demo_scattering_boundary_conditions.py
+++ b/python/demo/demo_scattering_boundary_conditions.py
@@ -220,8 +220,8 @@ def generate_mesh_wire(
 # +
 def compute_a(nu: int, m: complex, alpha: float) -> float:
     """Compute the a_nu coefficient."""
-    J_nu_alpha = jv(nu, alpha)
-    J_nu_malpha = jv(nu, m * alpha)
+    J_nu_alpha = jv(nu, alpha)  # type: ignore
+    J_nu_malpha = jv(nu, m * alpha)  # type: ignore
     J_nu_alpha_p = jvp(nu, alpha, 1)
     J_nu_malpha_p = jvp(nu, m * alpha, 1)
 

From ce1b91135aa93b5115dc653cfbd65250b072adbb Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Mon, 23 Mar 2026 13:26:28 +0100
Subject: [PATCH 05/22] Fix

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

diff --git a/python/demo/demo_pyamg.py b/python/demo/demo_pyamg.py
index 6027e028ce..0da47e25ed 100644
--- a/python/demo/demo_pyamg.py
+++ b/python/demo/demo_pyamg.py
@@ -63,7 +63,7 @@ def poisson_problem(dtype: npt.DTypeLike, solver_type: str) -> None:
         solver_type: pyamg solver type, either "ruge_stuben" or
             "smoothed_aggregation"
     """
-    real_type = np.real(dtype(0)).dtype  # type: ignore
+    real_type = np.real(np.zeros(0, dtype=dtype)).dtype
     mesh = create_box(
         comm=MPI.COMM_WORLD,
         points=[(0.0, 0.0, 0.0), (3.0, 2.0, 1.0)],
@@ -84,7 +84,7 @@ def poisson_problem(dtype: npt.DTypeLike, solver_type: str) -> None:
 
     dofs = locate_dofs_topological(V=V, entity_dim=fdim, entities=facets)
 
-    bc = dirichletbc(value=0.0, dofs=dofs, V=V)
+    bc = dirichletbc(value=dtype(0.0), dofs=dofs, V=V)  # type: ignore
 
     u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
     x = ufl.SpatialCoordinate(mesh)

From 0270c2ca23ceeeabc38f3fc7bcc7df76da7d41af Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 19 Mar 2026 10:50:48 +0100
Subject: [PATCH 06/22] Add typing generics to vector

---
 python/dolfinx/la/__init__.py | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/python/dolfinx/la/__init__.py b/python/dolfinx/la/__init__.py
index 8c3f40aaf4..283d2e3dbc 100644
--- a/python/dolfinx/la/__init__.py
+++ b/python/dolfinx/la/__init__.py
@@ -5,6 +5,8 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 """Linear algebra functionality."""
 
+from typing import Generic, TypeVar
+
 import numpy as np
 import numpy.typing as npt
 
@@ -26,7 +28,10 @@
 ]
 
 
-class Vector:
+_T = TypeVar("_T", np.float32, np.float64, np.complex64, np.complex128, np.int8, np.int32, np.int64)
+
+
+class Vector(Generic[_T]):
     """Distributed vector object."""
 
     _cpp_object: (
@@ -79,7 +84,7 @@ def block_size(self) -> int:
         return self._cpp_object.bs
 
     @property
-    def array(self) -> np.ndarray:
+    def array(self) -> npt.NDArray[_T]:
         """Local representation of the vector."""
         return self._cpp_object.array
 
@@ -335,17 +340,17 @@ def vector(map, bs=1, dtype: npt.DTypeLike = np.float64) -> Vector:
     return Vector(vtype(map, bs))
 
 
-def orthonormalize(basis: list[Vector]):
+def orthonormalize(basis: list[Vector[_T]]) -> None:
     """Orthogonalise set of vectors in-place."""
     _cpp.la.orthonormalize([x._cpp_object for x in basis])
 
 
-def is_orthonormal(basis: list[Vector], eps: float = 1.0e-12) -> bool:
+def is_orthonormal(basis: list[Vector[_T]], eps: float = 1.0e-12) -> bool:
     """Check that list of vectors are orthonormal."""
     return _cpp.la.is_orthonormal([x._cpp_object for x in basis], eps)
 
 
-def norm(x: Vector, type: _cpp.la.Norm = _cpp.la.Norm.l2) -> np.floating:
+def norm(x: Vector[_T], type: _cpp.la.Norm = _cpp.la.Norm.l2) -> np.floating:
     """Compute a norm of the vector.
 
     Args:

From 2f731b2d5e3d4cbd543656a8a80b132ac82de7ec Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 19 Mar 2026 10:58:20 +0100
Subject: [PATCH 07/22] Add typing generics to matrix

---
 python/dolfinx/la/__init__.py | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/python/dolfinx/la/__init__.py b/python/dolfinx/la/__init__.py
index 283d2e3dbc..ab008d6077 100644
--- a/python/dolfinx/la/__init__.py
+++ b/python/dolfinx/la/__init__.py
@@ -122,7 +122,10 @@ def scatter_reverse(self, mode: InsertMode) -> None:
         self._cpp_object.scatter_reverse(mode)
 
 
-class MatrixCSR:
+_MT = TypeVar("_MT", np.float32, np.float64, np.complex64, np.complex128)
+
+
+class MatrixCSR(Generic[_MT]):
     """Distributed compressed sparse row matrix."""
 
     _cpp_object: (
@@ -160,7 +163,7 @@ def index_map(self, i: int) -> IndexMap:
         """
         return self._cpp_object.index_map(i)
 
-    def mult(self, x: Vector, y: Vector) -> None:
+    def mult(self, x: Vector[_MT], y: Vector[_MT]) -> None:
         """Compute ``y += Ax``.
 
         Args:
@@ -176,7 +179,7 @@ def block_size(self) -> list:
 
     def add(
         self,
-        x: npt.NDArray[np.floating],
+        x: npt.NDArray[_MT],
         rows: npt.NDArray[np.int32],
         cols: npt.NDArray[np.int32],
         bs: int = 1,
@@ -186,7 +189,7 @@ def add(
 
     def set(
         self,
-        x: npt.NDArray[np.floating],
+        x: npt.NDArray[_MT],
         rows: npt.NDArray[np.int32],
         cols: npt.NDArray[np.int32],
         bs: int = 1,
@@ -194,7 +197,7 @@ def set(
         """Set a block of values in the matrix."""
         self._cpp_object.set(x, rows, cols, bs)
 
-    def set_value(self, x: np.floating) -> None:
+    def set_value(self, x: _MT) -> None:
         """Set all non-zero entries to a value.
 
         Args:
@@ -215,7 +218,7 @@ def squared_norm(self) -> np.floating:
         return self._cpp_object.squared_norm()
 
     @property
-    def data(self) -> npt.NDArray[np.floating]:
+    def data(self) -> npt.NDArray[_MT]:
         """Underlying matrix entry data."""
         return self._cpp_object.data
 
@@ -229,7 +232,7 @@ def indptr(self) -> npt.NDArray[np.int64]:
         """Local row pointers."""
         return self._cpp_object.indptr
 
-    def to_dense(self) -> npt.NDArray[np.floating]:
+    def to_dense(self) -> npt.NDArray[_MT]:
         """Copy to a dense 2D array.
 
         Note:

From 545902073ab29c9d91b551ec572589ff92c49046 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 19 Mar 2026 11:40:06 +0100
Subject: [PATCH 08/22] Add typing generics to superludist{matrix, solver}

---
 python/dolfinx/la/superlu_dist.py | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

diff --git a/python/dolfinx/la/superlu_dist.py b/python/dolfinx/la/superlu_dist.py
index 72ea934c43..4a5d8e9c8b 100644
--- a/python/dolfinx/la/superlu_dist.py
+++ b/python/dolfinx/la/superlu_dist.py
@@ -12,6 +12,8 @@
   Users with advanced linear solver requirements should use PETSc/petsc4py.
 """
 
+from typing import Generic, TypeVar
+
 import numpy as np
 import numpy.typing as npt
 
@@ -22,8 +24,10 @@
 
 __all__ = ["SuperLUDistMatrix", "SuperLUDistSolver", "superlu_dist_matrix", "superlu_dist_solver"]
 
+_T = TypeVar("_T", np.float32, np.float64, np.complex128)
+
 
-class SuperLUDistMatrix:
+class SuperLUDistMatrix(Generic[_T]):
     """SuperLU_DIST matrix."""
 
     _cpp_object: (
@@ -52,7 +56,7 @@ def dtype(self) -> npt.DTypeLike:
         return self._cpp_object.dtype
 
 
-def superlu_dist_matrix(A: dolfinx.la.MatrixCSR) -> SuperLUDistMatrix:
+def superlu_dist_matrix(A: dolfinx.la.MatrixCSR[_T]) -> SuperLUDistMatrix[_T]:
     """Create a SuperLU_DIST matrix.
 
     Deep copies all required data from ``A``.
@@ -75,7 +79,7 @@ def superlu_dist_matrix(A: dolfinx.la.MatrixCSR) -> SuperLUDistMatrix:
     return SuperLUDistMatrix(stype(A._cpp_object))
 
 
-class SuperLUDistSolver:
+class SuperLUDistSolver(Generic[_T]):
     """SuperLU_DIST solver."""
 
     _cpp_object: (
@@ -112,7 +116,7 @@ def set_option(self, name: str, value: str):
         """
         self._cpp_object.set_option(name, value)
 
-    def set_A(self, A: SuperLUDistMatrix):
+    def set_A(self, A: SuperLUDistMatrix[_T]):
         """Set assembled left-hand side matrix.
 
         For advanced use with SuperLU_DIST option `Factor` allowing use of
@@ -123,7 +127,7 @@ def set_A(self, A: SuperLUDistMatrix):
         """
         self._cpp_object.set_A(A._cpp_object)
 
-    def solve(self, b: dolfinx.la.Vector, u: dolfinx.la.Vector) -> int:
+    def solve(self, b: dolfinx.la.Vector[_T], u: dolfinx.la.Vector[_T]) -> int:
         """Solve linear system :math:`Au = b`.
 
         Note:
@@ -155,7 +159,7 @@ def solve(self, b: dolfinx.la.Vector, u: dolfinx.la.Vector) -> int:
         return self._cpp_object.solve(b._cpp_object, u._cpp_object)
 
 
-def superlu_dist_solver(A: SuperLUDistMatrix) -> SuperLUDistSolver:
+def superlu_dist_solver(A: SuperLUDistMatrix[_T]) -> SuperLUDistSolver[_T]:
     """Create a SuperLU_DIST linear solver.
 
     Solve linear system :math:`Au = b` via LU decomposition.

From f1333e9f2d94e0e8282fe879e5ffa2c92aaffca5 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 19 Mar 2026 14:49:52 +0100
Subject: [PATCH 09/22] Add typing generics to CoordinateElement

---
 python/dolfinx/fem/element.py | 21 ++++++++++-----------
 1 file changed, 10 insertions(+), 11 deletions(-)

diff --git a/python/dolfinx/fem/element.py b/python/dolfinx/fem/element.py
index 55ad6626f1..d0d1479f81 100644
--- a/python/dolfinx/fem/element.py
+++ b/python/dolfinx/fem/element.py
@@ -6,6 +6,7 @@
 """Finite elements."""
 
 from functools import singledispatch
+from typing import Generic, TypeVar
 
 import numpy as np
 import numpy.typing as npt
@@ -14,8 +15,10 @@
 import basix.ufl
 from dolfinx import cpp as _cpp
 
+_T = TypeVar("_T", np.float32, np.float64)
 
-class CoordinateElement:
+
+class CoordinateElement(Generic[_T]):
     """Coordinate element describing the geometry map for mesh cells."""
 
     _cpp_object: _cpp.fem.CoordinateElement_float32 | _cpp.fem.CoordinateElement_float64
@@ -60,11 +63,7 @@ def create_dof_layout(self) -> _cpp.fem.ElementDofLayout:
         """Compute and return the dof layout."""
         return self._cpp_object.create_dof_layout()
 
-    def push_forward(
-        self,
-        X: npt.NDArray[np.float32] | npt.NDArray[np.float64],
-        cell_geometry: npt.NDArray[np.float32] | npt.NDArray[np.float64],
-    ) -> npt.NDArray[np.float32] | npt.NDArray[np.float64]:
+    def push_forward(self, X: npt.NDArray[_T], cell_geometry: npt.NDArray[_T]) -> npt.NDArray[_T]:
         """Push points on the reference cell forward to the physical cell.
 
         Args:
@@ -82,11 +81,11 @@ def push_forward(
 
     def pull_back(
         self,
-        x: npt.NDArray[np.float32] | npt.NDArray[np.float64],
-        cell_geometry: npt.NDArray[np.float32] | npt.NDArray[np.float64],
+        x: npt.NDArray[_T],
+        cell_geometry: npt.NDArray[_T],
         tol: float = 1.0e-6,
         maxit: int = 15,
-    ) -> npt.NDArray[np.float32] | npt.NDArray[np.float64]:
+    ) -> npt.NDArray[_T]:
         """Pull points on the physical cell back to the reference cell.
 
         For non-affine cells, the pull-back is a nonlinear operation.
@@ -130,7 +129,7 @@ def coordinate_element(
     degree: int,
     variant=int(basix.LagrangeVariant.unset),
     dtype: npt.DTypeLike = np.float64,
-):
+) -> CoordinateElement:
     """Create a Lagrange CoordinateElement from element metadata.
 
     Coordinate elements are typically used to create meshes.
@@ -153,7 +152,7 @@ def coordinate_element(
 
 
 @coordinate_element.register(basix.finite_element.FiniteElement)
-def _(e: basix.finite_element.FiniteElement):
+def _(e: basix.finite_element.FiniteElement) -> CoordinateElement:
     """Create a Lagrange CoordinateElement from a Basix finite element.
 
     Coordinate elements are typically used when creating meshes.

From f9e3e9b3545f459bc2348ae5e523a58a8a2fe0e2 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 19 Mar 2026 14:52:32 +0100
Subject: [PATCH 10/22] Add typing generics to finiteelement

---
 python/dolfinx/fem/element.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/python/dolfinx/fem/element.py b/python/dolfinx/fem/element.py
index d0d1479f81..9e13b11f8f 100644
--- a/python/dolfinx/fem/element.py
+++ b/python/dolfinx/fem/element.py
@@ -169,7 +169,7 @@ def _(e: basix.finite_element.FiniteElement) -> CoordinateElement:
         return CoordinateElement(_cpp.fem.CoordinateElement_float64(e._e))
 
 
-class FiniteElement:
+class FiniteElement(Generic[_T]):
     """A finite element."""
 
     _cpp_object: _cpp.fem.FiniteElement_float32 | _cpp.fem.FiniteElement_float64
@@ -223,7 +223,7 @@ def value_shape(self) -> npt.NDArray[np.integer]:
         return self._cpp_object.value_shape
 
     @property
-    def interpolation_points(self) -> npt.NDArray[np.floating]:
+    def interpolation_points(self) -> npt.NDArray[_T]:
         """Points at which to evaluate the function to be interpolated.
 
         Interpolation point coordinates on the reference cell, returning
@@ -281,7 +281,7 @@ def signature(self) -> str:
         return self._cpp_object.signature
 
     def T_apply(
-        self, x: npt.NDArray[np.floating], cell_permutations: npt.NDArray[np.uint32], dim: int
+        self, x: npt.NDArray[_T], cell_permutations: npt.NDArray[np.uint32], dim: int
     ) -> None:
         """Transform basis from reference to physical ordering/orientation.
 
@@ -304,7 +304,7 @@ def T_apply(
         self._cpp_object.T_apply(x, cell_permutations, dim)
 
     def Tt_apply(
-        self, x: npt.NDArray[np.floating], cell_permutations: npt.NDArray[np.uint32], dim: int
+        self, x: npt.NDArray[_T], cell_permutations: npt.NDArray[np.uint32], dim: int
     ) -> None:
         """Apply the transpose of the operator applied by T_apply().
 
@@ -318,7 +318,7 @@ def Tt_apply(
         self._cpp_object.Tt_apply(x, cell_permutations, dim)
 
     def Tt_inv_apply(
-        self, x: npt.NDArray[np.floating], cell_permutations: npt.NDArray[np.uint32], dim: int
+        self, x: npt.NDArray[_T], cell_permutations: npt.NDArray[np.uint32], dim: int
     ) -> None:
         """Apply the inverse transpose of T_apply().
 
@@ -335,7 +335,7 @@ def Tt_inv_apply(
 def finiteelement(
     cell_type: _cpp.mesh.CellType,
     ufl_e: basix.ufl._ElementBase,
-    FiniteElement_dtype: np.dtype,
+    FiniteElement_dtype: npt.DTypeLike,
 ) -> FiniteElement:
     """Create a DOLFINx element from a basix.ufl element.
 

From 8d1c31e06fb024d33c87eed739e0a37e3828a928 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 19 Mar 2026 14:59:05 +0100
Subject: [PATCH 11/22] Add typing generics to constant

---
 python/dolfinx/fem/function.py | 21 ++++++++++++---------
 1 file changed, 12 insertions(+), 9 deletions(-)

diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index 69c9cce7f2..fc15227f50 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -11,6 +11,7 @@
 import typing
 from collections.abc import Callable, Sequence
 from functools import cached_property, singledispatch
+from typing import Generic, TypeVar
 
 import numpy as np
 import numpy.typing as npt
@@ -28,8 +29,10 @@
 
     from dolfinx.mesh import Mesh
 
+_S = TypeVar("_S", np.float32, np.float64, np.complex64, np.complex128)  # scalar
 
-class Constant(ufl.Constant):
+
+class Constant(ufl.Constant, Generic[_S]):
     """A constant with respect to a domain."""
 
     _cpp_object: (
@@ -72,27 +75,27 @@ def value(self):
         return self._cpp_object.value
 
     @value.setter
-    def value(self, v):
+    def value(self, v: npt.NDArray[_S]) -> None:
         np.copyto(self._cpp_object.value, np.asarray(v))
 
     @property
-    def dtype(self) -> np.dtype:
+    def dtype(self) -> npt.DTypeLike:
         """Value dtype of the constant."""
         return np.dtype(self._cpp_object.dtype)
 
-    def __float__(self):
+    def __float__(self) -> float:
         """Real representation of the constant."""
         if self.ufl_shape or self.ufl_free_indices:
             raise TypeError("Cannot evaluate a nonscalar expression to a scalar value.")
-        else:
-            return float(self.value)
 
-    def __complex__(self):
+        return float(self.value)
+
+    def __complex__(self) -> complex:
         """Complex representation of the constant."""
         if self.ufl_shape or self.ufl_free_indices:
             raise TypeError("Cannot evaluate a nonscalar expression to a scalar value.")
-        else:
-            return complex(self.value)
+
+        return complex(self.value)
 
 
 class Expression:

From 222800ea5b55cfd1ae7ff62c9fae30f142941d91 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 19 Mar 2026 16:25:08 +0100
Subject: [PATCH 12/22] Add typing generics dirichletbc

---
 python/dolfinx/fem/bcs.py | 21 ++++++++++++---------
 1 file changed, 12 insertions(+), 9 deletions(-)

diff --git a/python/dolfinx/fem/bcs.py b/python/dolfinx/fem/bcs.py
index d4abcaa8a7..5d32f7e015 100644
--- a/python/dolfinx/fem/bcs.py
+++ b/python/dolfinx/fem/bcs.py
@@ -13,6 +13,7 @@
 from __future__ import annotations
 
 from collections.abc import Callable, Iterable
+from typing import Generic, TypeVar
 
 import numpy as np
 import numpy.typing as npt
@@ -90,7 +91,10 @@ def locate_dofs_topological(
     return _cpp.fem.locate_dofs_topological(_V, entity_dim, _entities, remote)
 
 
-class DirichletBC:
+_T = TypeVar("_T", np.float32, np.float64, np.complex64, np.complex128)
+
+
+class DirichletBC(Generic[_T]):
     """Representation of Dirichlet boundary conditions.
 
     The conditions are imposed on a linear system.
@@ -119,8 +123,9 @@ class initialiser. This class is combined with different
         self._cpp_object = bc
 
     @property
-    def g(self) -> Function | Constant | np.ndarray:
+    def g(self) -> Function | Constant:
         """The boundary condition value(s)."""
+        # TODO: needs to be wrapped into Function or Constant
         return self._cpp_object.value
 
     @property
@@ -128,9 +133,7 @@ def function_space(self) -> dolfinx.fem.FunctionSpace:
         """Function space on which the boundary condition is defined."""
         return self._cpp_object.function_space
 
-    def set(
-        self, x: npt.NDArray, x0: npt.NDArray[np.int32] | None = None, alpha: float = 1
-    ) -> None:
+    def set(self, x: npt.NDArray[_T], x0: npt.NDArray[_T] | None = None, alpha: float = 1) -> None:
         """Set array entries that are constrained by a Dirichlet condition.
 
         Entries in ``x`` that are constrained by a Dirichlet boundary
@@ -171,10 +174,10 @@ def dof_indices(self) -> tuple[npt.NDArray[np.int32], int]:
 
 
 def dirichletbc(
-    value: Function | Constant | np.ndarray | float | complex,
+    value: Function | Constant | npt.NDArray[_T] | float | complex,
     dofs: npt.NDArray[np.int32],
     V: dolfinx.fem.FunctionSpace | None = None,
-) -> DirichletBC:
+) -> DirichletBC[_T]:
     """Representation of Dirichlet boundary condition.
 
     Args:
@@ -231,8 +234,8 @@ def dirichletbc(
 
 
 def bcs_by_block(
-    spaces: Iterable[FunctionSpace | None], bcs: Iterable[DirichletBC]
-) -> list[list[DirichletBC]]:
+    spaces: Iterable[FunctionSpace | None], bcs: Iterable[DirichletBC[_T]]
+) -> list[list[DirichletBC[_T]]]:
     """Arrange boundary conditions by the space that they constrain.
 
     Given a sequence of function spaces ``spaces`` and a sequence of

From 8fdf0ec77cd0ff51df47da94f9e810c028f0023b Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 19 Mar 2026 16:35:32 +0100
Subject: [PATCH 13/22] Add typing generics Expression

---
 python/dolfinx/fem/function.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index fc15227f50..5b654b7044 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -98,7 +98,7 @@ def __complex__(self) -> complex:
         return complex(self.value)
 
 
-class Expression:
+class Expression(Generic[_S]):
     """An object for evaluating functions of finite element functions.
 
     Represents a mathematical expression evaluated at a pre-defined set
@@ -213,9 +213,9 @@ def _create_expression(dtype):
     def eval(
         self,
         mesh: Mesh,
-        entities: np.ndarray,
-        values: np.ndarray | None = None,
-    ) -> np.ndarray:
+        entities: npt.NDArray[np.int32],
+        values: npt.NDArray[_S] | None = None,
+    ) -> npt.NDArray[_S]:
         """Evaluate Expression on entities.
 
         Args:
@@ -271,12 +271,12 @@ def eval(
         )
         return values
 
-    def X(self) -> np.ndarray:
+    def X(self) -> npt.NDArray:
         """Evaluation points on the reference cell."""
         return self._cpp_object.X()
 
     @property
-    def ufl_expression(self):
+    def ufl_expression(self) -> ufl.core.expr.Expr:
         """Original UFL Expression."""
         return self._ufl_expression
 
@@ -306,7 +306,7 @@ def code(self) -> str:
         return self._code
 
     @property
-    def dtype(self) -> np.dtype:
+    def dtype(self) -> npt.DTypeLike:
         """Expression value dtype."""
         return np.dtype(self._cpp_object.dtype)
 

From 275ba55cb81b099bb5d78eed65999eac5eb9b6b4 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 19 Mar 2026 16:40:13 +0100
Subject: [PATCH 14/22] Add typing generics Function

---
 python/dolfinx/fem/function.py | 34 ++++++++++++++++++----------------
 1 file changed, 18 insertions(+), 16 deletions(-)

diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index 5b654b7044..50ee0d9583 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -311,7 +311,7 @@ def dtype(self) -> npt.DTypeLike:
         return np.dtype(self._cpp_object.dtype)
 
 
-class Function(ufl.Coefficient):
+class Function(ufl.Coefficient, Generic[_S]):
     """A finite element function.
 
     A finite element function is represented by a function space
@@ -327,10 +327,12 @@ class Function(ufl.Coefficient):
         | _cpp.fem.Function_float64
     )
 
+    _x: la.Vector[_S]
+
     def __init__(
         self,
         V: FunctionSpace,
-        x: la.Vector | None = None,
+        x: la.Vector[_S] | None = None,
         name: str | None = None,
         dtype: npt.DTypeLike | None = None,
     ):
@@ -398,11 +400,11 @@ def function_space(self) -> FunctionSpace:
     def eval(
         self,
         x: npt.ArrayLike,
-        cells: npt.ArrayLike,
-        u: None | npt.NDArray[np.float32 | np.float64 | np.complex128 | np.complex64] = None,
+        cells: npt.NDArray[np.int32],
+        u: None | npt.NDArray[_S] = None,
         tol: float = 1.0e-6,
         maxit: int = 15,
-    ) -> np.ndarray:
+    ) -> npt.NDArray[_S]:
         """Evaluate Function at points x.
 
         Args:
@@ -446,7 +448,7 @@ def eval(
 
     def interpolate_nonmatching(
         self,
-        u0: Function,
+        u0: Function[_S],
         cells: npt.NDArray[np.int32],
         interpolation_data: PointOwnershipData,
         tol: float = 1e-6,
@@ -472,9 +474,9 @@ def interpolate_nonmatching(
 
     def interpolate(
         self,
-        u0: Callable | Expression | Function,
-        cells0: np.ndarray | None = None,
-        cells1: np.ndarray | None = None,
+        u0: Callable | Expression[_S] | Function[_S],
+        cells0: npt.NDArray[np.int32] | None = None,
+        cells1: npt.NDArray[np.int32] | None = None,
     ) -> None:
         """Interpolate an expression.
 
@@ -520,7 +522,7 @@ def _(e0: Expression):
             )
             self._cpp_object.interpolate_f(np.asarray(u0(x), dtype=self.dtype), cells0)
 
-    def copy(self) -> Function:
+    def copy(self) -> Function[_S]:
         """Create a copy of the Function.
 
         The function space is shared and the degree-of-freedom vector is
@@ -536,12 +538,12 @@ def copy(self) -> Function:
         )
 
     @property
-    def x(self) -> la.Vector:
+    def x(self) -> la.Vector[_S]:
         """Vector holding the degrees-of-freedom."""
         return self._x
 
     @property
-    def dtype(self) -> np.dtype:
+    def dtype(self) -> npt.DTypeLike:
         """Function value dtype."""
         return np.dtype(self._cpp_object.x.array.dtype)
 
@@ -554,11 +556,11 @@ def name(self) -> str:
     def name(self, name):
         self._cpp_object.name = name
 
-    def __str__(self):
+    def __str__(self) -> str:
         """Pretty print representation."""
         return self.name
 
-    def sub(self, i: int) -> Function:
+    def sub(self, i: int) -> Function[_S]:
         """Return a sub-function (a view into the ``Function``).
 
         Sub-functions are indexed ``i = 0, ..., N-1``, where ``N`` is
@@ -577,7 +579,7 @@ def sub(self, i: int) -> Function:
         """
         return Function(self._V.sub(i), self.x, name=f"{self!s}_{i}")
 
-    def split(self) -> tuple[Function, ...]:
+    def split(self) -> tuple[Function[_S], ...]:
         """Extract (any) sub-functions.
 
         A sub-function can be extracted from a discrete function that is
@@ -592,7 +594,7 @@ def split(self) -> tuple[Function, ...]:
             raise RuntimeError("No subfunctions to extract")
         return tuple(self.sub(i) for i in range(num_sub_spaces))
 
-    def collapse(self) -> Function:
+    def collapse(self) -> Function[_S]:
         """Create a collapsed version of this Function."""
         u_collapsed = self._cpp_object.collapse()  # type: ignore
         V_collapsed = FunctionSpace(

From 791870c61f61bcf5ceda84b6ba74bd33e6715249 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Mon, 23 Mar 2026 16:05:04 +0100
Subject: [PATCH 15/22] [tmp] ci setup

---
 .github/workflows/ccpp.yml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/ccpp.yml b/.github/workflows/ccpp.yml
index 21fe82f6f9..c2d56ae74d 100644
--- a/.github/workflows/ccpp.yml
+++ b/.github/workflows/ccpp.yml
@@ -7,8 +7,8 @@ on:
     tags:
       - "v*"
   pull_request:
-    branches:
-      - main
+    # branches:
+    #   - main
   merge_group:
     branches:
       - main

From 800f79694b9aca0075f8e1ebab6a77fe2f3f776c Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Tue, 24 Mar 2026 18:30:13 +0100
Subject: [PATCH 16/22] Fix: demo_lagrange

---
 python/demo/demo_lagrange_variants.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/python/demo/demo_lagrange_variants.py b/python/demo/demo_lagrange_variants.py
index 209e7312c6..135f12a863 100644
--- a/python/demo/demo_lagrange_variants.py
+++ b/python/demo/demo_lagrange_variants.py
@@ -26,6 +26,7 @@
 from mpi4py import MPI
 
 import matplotlib.pylab as plt
+import numpy as np
 
 import basix
 import basix.ufl
@@ -142,11 +143,11 @@ def saw_tooth(x):
     uh.interpolate(lambda x: saw_tooth(x[0]))
     if MPI.COMM_WORLD.size == 1:  # Skip this plotting in parallel
         pts: list[list[float]] = []
-        cells: list[int] = []
+        cells = np.empty((0,), dtype=np.int32)
         for cell in range(N):
             for i in range(51):
                 pts.append([cell / N + i / 50 / N, 0, 0])
-                cells.append(cell)
+                cells = np.append(cells, [cell])
         values = uh.eval(pts, cells)
         plt.plot(pts, [saw_tooth(i[0]) for i in pts], "k--")
         plt.plot(pts, values, "r-")

From 69420739fe78ecb36af5f46973cb31acfda1f512 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Tue, 24 Mar 2026 18:37:20 +0100
Subject: [PATCH 17/22] Add typing generics Geometry

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

diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index c326acc04b..1ae21404cf 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -82,6 +82,9 @@
 ]
 
 
+_T = typing.TypeVar("_T", np.float32, np.float64)
+
+
 @singledispatch
 def create_cell_partitioner(
     part: Callable, mode: GhostMode, max_facet_to_cell_links: int
@@ -267,7 +270,7 @@ def cell_type(self) -> CellType:
         return self._cpp_object.cell_type
 
 
-class Geometry:
+class Geometry(typing.Generic[_T]):
     """The geometry of a :class:`dolfinx.mesh.Mesh`."""
 
     _cpp_object: _cpp.mesh.Geometry_float32 | _cpp.mesh.Geometry_float64
@@ -314,7 +317,7 @@ def input_global_indices(self) -> npt.NDArray[np.int64]:
         return self._cpp_object.input_global_indices
 
     @property
-    def x(self) -> npt.NDArray[np.float32] | npt.NDArray[np.float64]:
+    def x(self) -> npt.NDArray[_T]:
         """Geometry coordinate points.
 
         Shape is ``shape=(num_points, 3)``.

From b279c81d6c44c65045cc121c41a08a6e6cab365b Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Tue, 24 Mar 2026 18:38:54 +0100
Subject: [PATCH 18/22] Add typing generics Mesh

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

diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index 1ae21404cf..1b5fbf7fcb 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -325,12 +325,12 @@ def x(self) -> npt.NDArray[_T]:
         return self._cpp_object.x
 
 
-class Mesh:
+class Mesh(typing.Generic[_T]):
     """A mesh."""
 
     _mesh: _cpp.mesh.Mesh_float32 | _cpp.mesh.Mesh_float64
     _topology: Topology
-    _geometry: Geometry
+    _geometry: Geometry[_T]
     _ufl_domain: ufl.Mesh | None
 
     def __init__(
@@ -414,7 +414,7 @@ def topology(self) -> Topology:
         return self._topology
 
     @property
-    def geometry(self) -> Geometry:
+    def geometry(self) -> Geometry[_T]:
         """Mesh geometry."""
         return self._geometry
 

From 1a13cc89f73818c90a7407a8650e01d111c85b24 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Tue, 24 Mar 2026 18:41:21 +0100
Subject: [PATCH 19/22] Add typing generics PointOwnershipData

---
 python/dolfinx/geometry.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/python/dolfinx/geometry.py b/python/dolfinx/geometry.py
index 9c9d734797..9600551042 100644
--- a/python/dolfinx/geometry.py
+++ b/python/dolfinx/geometry.py
@@ -34,8 +34,10 @@
     "squared_distance",
 ]
 
+_T = typing.TypeVar("_T", np.float32, np.float64)
 
-class PointOwnershipData:
+
+class PointOwnershipData(typing.Generic[_T]):
     """Class for storing data related to the ownership of points."""
 
     _cpp_object: _cpp.geometry.PointOwnershipData_float32 | _cpp.geometry.PointOwnershipData_float64
@@ -55,7 +57,7 @@ def dest_owner(self) -> npt.NDArray[np.int32]:
         return self._cpp_object.dest_owners
 
     @property
-    def dest_points(self) -> npt.NDArray[np.floating]:
+    def dest_points(self) -> npt.NDArray[_T]:
         """Points owned by current rank."""
         return self._cpp_object.dest_points
 
@@ -321,10 +323,10 @@ def compute_distances_gjk(
 
 def determine_point_ownership(
     mesh: Mesh,
-    points: npt.NDArray[np.floating],
+    points: npt.NDArray[_T],
     padding: float,
     cells: npt.NDArray[np.int32] | None = None,
-) -> PointOwnershipData:
+) -> PointOwnershipData[_T]:
     """Build point ownership data for a mesh-points pair.
 
     First, potential collisions are found by computing intersections

From 8e7da84a528c9a109bd1538d1b4e9d60308a8632 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Tue, 24 Mar 2026 18:45:39 +0100
Subject: [PATCH 20/22] Add typing generics BoundingBoxTree

---
 python/dolfinx/geometry.py | 42 +++++++++++++++++++-------------------
 1 file changed, 21 insertions(+), 21 deletions(-)

diff --git a/python/dolfinx/geometry.py b/python/dolfinx/geometry.py
index 9600551042..900bb5b1f9 100644
--- a/python/dolfinx/geometry.py
+++ b/python/dolfinx/geometry.py
@@ -67,7 +67,7 @@ def dest_cells(self) -> npt.NDArray[np.int32]:
         return self._cpp_object.dest_cells
 
 
-class BoundingBoxTree:
+class BoundingBoxTree(typing.Generic[_T]):
     """Bounding box trees used in collision detection."""
 
     _cpp_object: _cpp.geometry.BoundingBoxTree_float32 | _cpp.geometry.BoundingBoxTree_float64
@@ -88,7 +88,7 @@ def num_bboxes(self) -> int:
         return self._cpp_object.num_bboxes
 
     @property
-    def bbox_coordinates(self) -> npt.NDArray[np.float32] | npt.NDArray[np.float64]:
+    def bbox_coordinates(self) -> npt.NDArray[_T]:
         """Coordinates of lower and upper corners of bounding boxes.
 
         Note:
@@ -97,7 +97,7 @@ def bbox_coordinates(self) -> npt.NDArray[np.float32] | npt.NDArray[np.float64]:
         """
         return self._cpp_object.bbox_coordinates
 
-    def get_bbox(self, i) -> npt.NDArray[np.floating]:
+    def get_bbox(self, i) -> npt.NDArray[_T]:
         """Get lower and upper corners of the ith bounding box.
 
         Args:
@@ -110,18 +110,18 @@ def get_bbox(self, i) -> npt.NDArray[np.floating]:
         """
         return self._cpp_object.get_bbox(i)
 
-    def create_global_tree(self, comm) -> BoundingBoxTree:
+    def create_global_tree(self, comm) -> BoundingBoxTree[_T]:
         """Create a global bounding box tree."""
         return BoundingBoxTree(self._cpp_object.create_global_tree(comm))
 
 
 def bb_tree(
-    mesh: Mesh,
+    mesh: Mesh[_T],
     dim: int,
     *,
     padding: float = 0.0,
     entities: npt.NDArray[np.int32] | None = None,
-) -> BoundingBoxTree:
+) -> BoundingBoxTree[_T]:
     """Create a bounding box tree for use in collision detection.
 
     Args:
@@ -153,7 +153,7 @@ def bb_tree(
 
 
 def compute_collisions_trees(
-    tree0: BoundingBoxTree, tree1: BoundingBoxTree
+    tree0: BoundingBoxTree[_T], tree1: BoundingBoxTree[_T]
 ) -> npt.NDArray[np.int32]:
     """Compute all collisions between two bounding box trees.
 
@@ -169,7 +169,7 @@ def compute_collisions_trees(
     return _cpp.geometry.compute_collisions_trees(tree0._cpp_object, tree1._cpp_object)
 
 
-def compute_collisions_points(tree: BoundingBoxTree, x: npt.NDArray[np.floating]) -> AdjacencyList:
+def compute_collisions_points(tree: BoundingBoxTree[_T], x: npt.NDArray[_T]) -> AdjacencyList:
     """Compute collisions between points and leaf bounding boxes.
 
     Bounding boxes can overlap, therefore points can collide with more
@@ -188,10 +188,10 @@ def compute_collisions_points(tree: BoundingBoxTree, x: npt.NDArray[np.floating]
 
 
 def compute_closest_entity(
-    tree: BoundingBoxTree,
-    midpoint_tree: BoundingBoxTree,
-    mesh: Mesh,
-    points: npt.NDArray[np.floating],
+    tree: BoundingBoxTree[_T],
+    midpoint_tree: BoundingBoxTree[_T],
+    mesh: Mesh[_T],
+    points: npt.NDArray[_T],
 ) -> npt.NDArray[np.int32]:
     """Compute closest mesh entity to a point.
 
@@ -213,7 +213,9 @@ def compute_closest_entity(
     )
 
 
-def create_midpoint_tree(mesh: Mesh, dim: int, entities: npt.NDArray[np.int32]) -> BoundingBoxTree:
+def create_midpoint_tree(
+    mesh: Mesh[_T], dim: int, entities: npt.NDArray[np.int32]
+) -> BoundingBoxTree[_T]:
     """Create bounding box tree for the midpoints of a subset of entities.
 
     Args:
@@ -228,7 +230,7 @@ def create_midpoint_tree(mesh: Mesh, dim: int, entities: npt.NDArray[np.int32])
 
 
 def compute_colliding_cells(
-    msh: Mesh, candidates: AdjacencyList, x: npt.NDArray[np.floating]
+    msh: Mesh[_T], candidates: AdjacencyList, x: npt.NDArray[_T]
 ) -> AdjacencyList:
     """From a mesh, find which cells collide with a set of points.
 
@@ -249,8 +251,8 @@ def compute_colliding_cells(
 
 
 def squared_distance(
-    mesh: Mesh, dim: int, entities: npt.NDArray[np.int32], points: npt.NDArray[np.floating]
-) -> npt.NDArray[np.floating]:
+    mesh: Mesh[_T], dim: int, entities: npt.NDArray[np.int32], points: npt.NDArray[_T]
+) -> npt.NDArray[_T]:
     """Compute the squared distance between a point and a mesh entity.
 
     The distance is computed between the ith input points and the ith
@@ -270,9 +272,7 @@ def squared_distance(
     return _cpp.geometry.squared_distance(mesh._cpp_object, dim, entities, points)
 
 
-def compute_distance_gjk(
-    p: npt.NDArray[np.floating], q: npt.NDArray[np.floating]
-) -> npt.NDArray[np.floating]:
+def compute_distance_gjk(p: npt.NDArray[_T], q: npt.NDArray[_T]) -> npt.NDArray[_T]:
     """Compute the distance between two convex bodies.
 
     Each body is defined by a set of points. Uses the
@@ -294,8 +294,8 @@ def compute_distance_gjk(
 
 
 def compute_distances_gjk(
-    bodies: list[npt.NDArray[np.floating]], q: npt.NDArray[np.floating], num_threads: int
-) -> npt.NDArray[np.floating]:
+    bodies: list[npt.NDArray[_T]], q: npt.NDArray[_T], num_threads: int
+) -> npt.NDArray[_T]:
     """Compute the distance between a set of convex bodies.
 
     For each convex body defined in `bodies`;

From de6166cc81559454c7860600c9f551ecd0a504f6 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Tue, 24 Mar 2026 18:50:11 +0100
Subject: [PATCH 21/22] Add typing generics FunctionSpace

---
 python/dolfinx/fem/function.py | 21 +++++++++++----------
 1 file changed, 11 insertions(+), 10 deletions(-)

diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index 50ee0d9583..d20cd6ea9b 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -29,6 +29,7 @@
 
     from dolfinx.mesh import Mesh
 
+_T = TypeVar("_T", np.float32, np.float64)
 _S = TypeVar("_S", np.float32, np.float64, np.complex64, np.complex128)  # scalar
 
 
@@ -674,11 +675,11 @@ def functionspace(
     return FunctionSpace(mesh, ufl_e, cppV)
 
 
-class FunctionSpace(ufl.FunctionSpace):
+class FunctionSpace(ufl.FunctionSpace, Generic[_T]):
     """A space on which Functions (fields) can be defined."""
 
     _cpp_object: _cpp.fem.FunctionSpace_float32 | _cpp.fem.FunctionSpace_float64
-    _mesh: Mesh
+    _mesh: Mesh[_T]
 
     def __init__(
         self,
@@ -705,7 +706,7 @@ def __init__(
         self._mesh = mesh
         super().__init__(ufl_domain, element)
 
-    def clone(self) -> FunctionSpace:
+    def clone(self) -> FunctionSpace[_T]:
         """Create a FunctionSpace which shares data with this space.
 
         The new space has a different unique integer ID.
@@ -740,7 +741,7 @@ def num_sub_spaces(self) -> int:
         """Number of sub spaces."""
         return self.element.num_sub_elements
 
-    def sub(self, i: int) -> FunctionSpace:
+    def sub(self, i: int) -> FunctionSpace[_T]:
         """Return the i-th sub space.
 
         Args:
@@ -788,7 +789,7 @@ def ufl_function_space(self) -> ufl.FunctionSpace:
         return self
 
     @cached_property
-    def element(self) -> FiniteElement:
+    def element(self) -> FiniteElement[_T]:
         """Function space finite element."""
         return FiniteElement(self._cpp_object.element)
 
@@ -802,22 +803,22 @@ def dofmaps(self, idx: int) -> DofMap:
         return DofMap(self._cpp_object.dofmaps(idx))
 
     @property
-    def mesh(self) -> Mesh:
+    def mesh(self) -> Mesh[_T]:
         """Mesh on which the function space is defined."""
         return self._mesh
 
-    def collapse(self) -> tuple[FunctionSpace, np.ndarray]:
+    def collapse(self) -> tuple[FunctionSpace[_T], list[npt.NDArray[np.int32]]]:
         """Create a new function space by collapsing a subspace.
 
         Returns:
             A new function space and the map from new to old
             degrees-of-freedom.
         """
-        cpp_space, dofs = self._cpp_object.collapse()  # type: ignore
+        cpp_space, dofs = self._cpp_object.collapse()
         V = FunctionSpace(self._mesh, self.ufl_element(), cpp_space)
         return V, dofs
 
-    def tabulate_dof_coordinates(self) -> npt.NDArray[np.float64]:
+    def tabulate_dof_coordinates(self) -> npt.NDArray[_T]:
         """Tabulate coordinates of function space degrees-of-freedom.
 
         Returns:
@@ -827,4 +828,4 @@ def tabulate_dof_coordinates(self) -> npt.NDArray[np.float64]:
             This method is only for elements with point evaluation
             degrees-of-freedom.
         """
-        return self._cpp_object.tabulate_dof_coordinates()  # type: ignore
+        return self._cpp_object.tabulate_dof_coordinates()

From 59448366fcb563b7f684157f933511c252b3fa82 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Tue, 24 Mar 2026 18:51:45 +0100
Subject: [PATCH 22/22] Add typing generics Form

---
 python/dolfinx/fem/forms.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py
index 242c88ded1..fbcc6377aa 100644
--- a/python/dolfinx/fem/forms.py
+++ b/python/dolfinx/fem/forms.py
@@ -32,8 +32,10 @@
     from dolfinx.mesh import EntityMap as _EntityMap
     from dolfinx.mesh import Mesh, MeshTags
 
+_S = typing.TypeVar("_S", np.float32, np.float64, np.complex64, np.complex128)  # scalar
 
-class Form:
+
+class Form(typing.Generic[_S]):
     """A finite element form."""
 
     _cpp_object: (