diff --git a/python/demo/demo_axis.py b/python/demo/demo_axis.py
index f1a8385fb8f..072ee977879 100644
--- a/python/demo/demo_axis.py
+++ b/python/demo/demo_axis.py
@@ -154,8 +154,8 @@ def generate_mesh_sphere_axis(
 # \cdot (\nabla \times \bar{\mathbf{v}})+\varepsilon_{r} k_{0}^{2}
 # \mathbf{E}_s \cdot \bar{\mathbf{v}}+k_{0}^{2}\left(\varepsilon_{r}
 # -\varepsilon_b\right)\mathbf{E}_b \cdot \bar{\mathbf{v}}~\mathrm{d} x\\
-# +\int_{\Omega_{pml}}\left[\boldsymbol{\mu}^{-1}_{pml} \nabla \times \mathbf{E}_s
-# \right]\cdot \nabla \times \bar{\mathbf{v}}-k_{0}^{2}
+# +\int_{\Omega_{pml}}\left[\boldsymbol{\mu}^{-1}_{pml} \nabla \times
+# \mathbf{E}_s \right]\cdot \nabla \times \bar{\mathbf{v}}-k_{0}^{2}
 # \left[\boldsymbol{\varepsilon}_{pml} \mathbf{E}_s \right]\cdot
 # \bar{\mathbf{v}}~ d x=0
 # \end{split}
@@ -172,8 +172,8 @@ def generate_mesh_sphere_axis(
 # -\varepsilon_b\right)\mathbf{E}_b \cdot
 # \bar{\mathbf{v}}~ \rho d\rho dz d \phi\\
 # +\int_{\Omega_{cs}}
-# \int_{0}^{2\pi}\left[\boldsymbol{\mu}^{-1}_{pml} \nabla \times \mathbf{E}_s
-# \right]\cdot \nabla \times \bar{\mathbf{v}}-k_{0}^{2}
+# \int_{0}^{2\pi}\left[\boldsymbol{\mu}^{-1}_{pml} \nabla \times
+# \mathbf{E}_s \right]\cdot \nabla \times \bar{\mathbf{v}}-k_{0}^{2}
 # \left[\boldsymbol{\varepsilon}_{pml} \mathbf{E}_s \right]\cdot
 # \bar{\mathbf{v}}~ \rho d\rho dz d \phi=0
 # \end{split}
@@ -184,8 +184,10 @@ def generate_mesh_sphere_axis(
 #
 # $$
 # \begin{align}
-# \mathbf{E}_s(\rho, z, \phi) &= \sum_m\mathbf{E}^{(m)}_s(\rho, z)e^{-jm\phi} \\
-# \mathbf{E}_b(\rho, z, \phi) &= \sum_m\mathbf{E}^{(m)}_b(\rho, z)e^{-jm\phi} \\
+# \mathbf{E}_s(\rho, z, \phi) &= \sum_m\mathbf{E}^{(m)}_s(\rho, z)
+#   e^{-jm\phi} \\
+# \mathbf{E}_b(\rho, z, \phi) &= \sum_m\mathbf{E}^{(m)}_b(\rho, z)
+#   e^{-jm\phi} \\
 # \bar{\mathbf{v}}(\rho, z, \phi) &=
 # \sum_m\bar{\mathbf{v}}^{(m)}(\rho, z)e^{+jm\phi}
 # \end{align}
@@ -748,8 +750,10 @@ def create_eps_mu(pml, rho, eps_bkg, mu_bkg):
 # m = np.sqrt(eps_au)/n_bkg
 # x = 2*np.pi*radius_sph/wl0*n_bkg
 #
-# q_sca_analyt, q_abs_analyt = scattnlay(np.array([x], dtype=np.complex128),
-#                                        np.array([m], dtype=np.complex128))[2:4]
+# q_sca_analyt, q_abs_analyt = scattnlay(
+#   np.array([x], dtype=np.complex128),
+#   np.array([m], dtype=np.complex128)
+# )[2:4]
 # ```
 #
 # The numerical values are reported here below:
diff --git a/python/demo/demo_biharmonic.py b/python/demo/demo_biharmonic.py
index 4ed85951321..cce2f6cdb83 100644
--- a/python/demo/demo_biharmonic.py
+++ b/python/demo/demo_biharmonic.py
@@ -33,23 +33,23 @@
 # \nabla^{4} u = f \quad {\rm in} \ \Omega,
 # $$
 #
-# where $\nabla^{4} \equiv \nabla^{2} \nabla^{2}$ is the biharmonic operator
-# and $f$ is a prescribed source term.
+# where $\nabla^{4} \equiv \nabla^{2} \nabla^{2}$ is the biharmonic
+# operator and $f$ is a prescribed source term.
 # To formulate a complete boundary value problem, the biharmonic equation
 # must be complemented by suitable boundary conditions.
 #
 # ### Weak formulation
 #
 # Multiplying the biharmonic equation by a test function and integrating
-# by parts twice leads to a problem of second-order derivatives, which would
-# require $H^{2}$ conforming (roughly $C^{1}$ continuous) basis functions.
-# To solve the biharmonic equation using Lagrange finite element basis
-# functions, the biharmonic equation can be split into two second-order
-# equations (see the Mixed Poisson demo for a mixed method for the Poisson
-# equation), or a variational formulation can be constructed that imposes
-# weak continuity of normal derivatives between finite element cells.
-# This demo uses a discontinuous Galerkin approach to impose continuity
-# of the normal derivative weakly.
+# by parts twice leads to a problem of second-order derivatives, which
+# would require $H^{2}$ conforming (roughly $C^{1}$ continuous) basis
+# functions. To solve the biharmonic equation using Lagrange finite element
+# basis functions, the biharmonic equation can be split into two second-
+# order equations (see the Mixed Poisson demo for a mixed method for the
+# Poisson equation), or a variational formulation can be constructed that
+# imposes weak continuity of normal derivatives between finite element
+# cells. This demo uses a discontinuous Galerkin approach to impose
+# continuity of the normal derivative weakly.
 #
 # Consider a triangulation $\mathcal{T}$ of the domain $\Omega$, where
 # the set of interior facets is denoted by $\mathcal{E}_h^{\rm int}$.
@@ -71,7 +71,8 @@
 # \end{align}
 # $$
 #
-# a weak formulation of the biharmonic problem reads: find $u \in V$ such that
+# a weak formulation of the biharmonic problem reads: find $u \in V$ such
+# that
 #
 # $$
 # a(u,v)=L(v) \quad \forall \ v \in V,
@@ -158,8 +159,9 @@
 # Next, we locate the mesh facets that lie on the boundary
 # $\Gamma_D = \partial\Omega$.
 # We do this using using {py:func}`exterior_facet_indices
-# <dolfinx.mesh.exterior_facet_indices>` which returns all mesh boundary facets
-# (Note: if we are only interested in a subset of those, consider {py:func}`locate_entities_boundary
+# <dolfinx.mesh.exterior_facet_indices>` which returns all mesh boundary
+# facets (Note: if we are only interested in a subset of those, consider
+# {py:func}`locate_entities_boundary
 # <dolfinx.mesh.locate_entities_boundary>`).
 
 tdim = msh.topology.dim
@@ -182,8 +184,8 @@
 
 # Next, we express the variational problem using UFL.
 #
-# First, the penalty parameter $\alpha$ is defined. In addition, we define a
-# variable `h` for the cell diameter $h_E$, a variable `n`for the
+# First, the penalty parameter $\alpha$ is defined. In addition, we define
+# a variable `h` for the cell diameter $h_E$, a variable `n`for the
 # outward-facing normal vector $n$ and a variable `h_avg` for the
 # average size of cells sharing a facet
 # $\left< h \right> = \frac{1}{2} (h_{+} + h_{-})$. Here, the UFL syntax
@@ -195,12 +197,12 @@
 n = ufl.FacetNormal(msh)
 h_avg = (h("+") + h("-")) / 2.0
 
-# After that, we can define the variational problem consisting of the bilinear
-# form $a$ and the linear form $L$. The source term is prescribed as
-# $f = 4.0 \pi^4\sin(\pi x)\sin(\pi y)$. Note that with `dS`, integration is
-# carried out over all the interior facets $\mathcal{E}_h^{\rm int}$, whereas
-# with `ds` it would be only the facets on the boundary of the domain, i.e.
-# $\partial\Omega$. The jump operator
+# After that, we can define the variational problem consisting of the
+# bilinear form $a$ and the linear form $L$. The source term is prescribed
+# as $f = 4.0 \pi^4\sin(\pi x)\sin(\pi y)$. Note that with `dS`,
+# integration is carried out over all the interior facets
+# $\mathcal{E}_h^{\rm int}$, whereas with `ds` it would be only the facets
+# on the boundary of the domain, i.e. $\partial\Omega$. The jump operator
 # $[\!\![ w ]\!\!] = w_{+} \cdot n_{+} + w_{-} \cdot n_{-}$ w.r.t. the
 # outward-facing normal vector $n$ is in UFL available as `jump(w, n)`.
 
diff --git a/python/demo/demo_cahn-hilliard.py b/python/demo/demo_cahn-hilliard.py
index 85e7aa45e91..8f0cde6ad4d 100644
--- a/python/demo/demo_cahn-hilliard.py
+++ b/python/demo/demo_cahn-hilliard.py
@@ -62,7 +62,8 @@
 # \begin{align}
 # \frac{\partial c}{\partial t} - \nabla \cdot M \nabla\mu
 #     &= 0 \quad {\rm in} \ \Omega, \\
-# \mu -  \frac{d f}{d c} + \lambda \nabla^{2}c &= 0 \quad {\rm in} \ \Omega.
+# \mu -  \frac{d f}{d c} + \lambda \nabla^{2}c &= 0 \quad {\rm in}
+#   \ \Omega.
 # \end{align}
 # $$
 #
@@ -74,8 +75,8 @@
 # \int_{\Omega} \frac{\partial c}{\partial t} q \, {\rm d} x +
 #     \int_{\Omega} M \nabla\mu \cdot \nabla q \, {\rm d} x
 #     &= 0 \quad \forall \ q \in V,  \\
-# \int_{\Omega} \mu v \, {\rm d} x - \int_{\Omega} \frac{d f}{d c} v \, {\rm d} x
-#   - \int_{\Omega} \lambda \nabla c \cdot \nabla v \, {\rm d} x
+# \int_{\Omega} \mu v \, {\rm d} x - \int_{\Omega} \frac{d f}{d c} v \,
+#   {\rm d} x - \int_{\Omega} \lambda \nabla c \cdot \nabla v \, {\rm d} x
 #    &= 0 \quad \forall \ v \in V.
 # \end{align}
 # $$
@@ -89,15 +90,16 @@
 # $$
 # \begin{align}
 # \int_{\Omega} \frac{c_{n+1} - c_{n}}{dt} q \, {\rm d} x
-# + \int_{\Omega} M \nabla \mu_{n+\theta} \cdot \nabla q \, {\rm d} x
-#        &= 0 \quad \forall \ q \in V  \\
-# \int_{\Omega} \mu_{n+1} v  \, {\rm d} x - \int_{\Omega} \frac{d f_{n+1}}{d c} v  \, {\rm d} x
-# - \int_{\Omega} \lambda \nabla c_{n+1} \cdot \nabla v \, {\rm d} x
-#        &= 0 \quad \forall \ v \in V
+#   + \int_{\Omega} M \nabla \mu_{n+\theta} \cdot \nabla q \, {\rm d} x
+#   &= 0 \quad \forall \ q \in V  \\
+# \int_{\Omega} \mu_{n+1} v  \, {\rm d} x - \int_{\Omega}
+#   \frac{d f_{n+1}}{d c} v  \, {\rm d} x - \int_{\Omega} \lambda \nabla
+#   c_{n+1} \cdot \nabla v \, {\rm d} x &= 0 \quad \forall \ v \in V
 # \end{align}
 # $$
 #
-# where $dt = t_{n+1} - t_{n}$ and $\mu_{n+\theta} = (1-\theta) \mu_{n} + \theta \mu_{n+1}$.
+# where $dt = t_{n+1} - t_{n}$ and $\mu_{n+\theta} = (1-\theta) \mu_{n} +
+# \theta \mu_{n+1}$.
 # The task is: given $c_{n}$ and $\mu_{n}$, solve the above equation to
 # find $c_{n+1}$ and $\mu_{n+1}$.
 #
diff --git a/python/demo/demo_hdg.py b/python/demo/demo_hdg.py
index 5fec78cc4e3..49fcfa6401c 100644
--- a/python/demo/demo_hdg.py
+++ b/python/demo/demo_hdg.py
@@ -103,8 +103,8 @@ def u_e(x):
 facets = np.arange(num_facets, dtype=np.int32)
 
 # Create the sub-mesh
-# NOTE Despite all facets being present in the submesh, the entity map isn't
-# necessarily the identity in parallel
+# NOTE Despite all facets being present in the submesh, the entity map
+# isn't necessarily the identity in parallel
 facet_mesh, facet_mesh_to_mesh = mesh.create_submesh(msh, fdim, facets)[:2]
 
 # Define function spaces
diff --git a/python/demo/demo_lagrange_variants.py b/python/demo/demo_lagrange_variants.py
index ef7153be181..7fd89096655 100644
--- a/python/demo/demo_lagrange_variants.py
+++ b/python/demo/demo_lagrange_variants.py
@@ -37,8 +37,9 @@
 # influence of interpolation point position, we create a degree 10
 # element on an interval using equally spaced points, and plot the basis
 # functions. We create this element using `basix.ufl`'s
-# `element` function. The function `element.tabulate` returns a 3-dimensional
-# array with shape (derivatives, points, (value size) * (basis functions)).
+# `element` function. The function `element.tabulate` returns a 3-
+# dimensional array with shape (derivatives, points, (value size) *
+# (basis functions)).
 # In this example, we only tabulate the 0th derivative and the value
 # size is 1, so we take the slice `[0, :, :]` to get a 2-dimensional
 # array.
@@ -158,7 +159,8 @@ def saw_tooth(x):
 # Lagrange compared to the GLL variant. To quantify the error, we
 # compute the interpolation error in the $L_2$ norm,
 #
-# $$\left\|u - u_h\right\|_2 = \left(\int_0^1 (u - u_h)^2\right)^{\frac{1}{2}},$$
+# $$\left\|u - u_h\right\|_2 =
+#   \left(\int_0^1 (u - u_h)^2\right)^{\frac{1}{2}},$$
 #
 # where $u$ is the function and $u_h$ is its interpolation in the finite
 # element space. The following code uses UFL to compute the $L_2$ error
diff --git a/python/demo/demo_mixed-poisson.py b/python/demo/demo_mixed-poisson.py
index 6ad1238c68b..c8a99b4df87 100644
--- a/python/demo/demo_mixed-poisson.py
+++ b/python/demo/demo_mixed-poisson.py
@@ -8,7 +8,7 @@
 #       jupytext_version: 1.14.1
 # ---
 
-# # Mixed formulation of the Poisson equation with a block-preconditioner/solver
+# # Mixed formulation of the Poisson equation with a block-preconditioner/solver # noqa
 #
 # This demo illustrates how to solve the Poisson equation using a mixed
 # (two-field) formulation and a block-preconditioned iterative solver.
diff --git a/python/demo/demo_mixed-topology.py b/python/demo/demo_mixed-topology.py
index adec8ba5de5..292d7584060 100644
--- a/python/demo/demo_mixed-topology.py
+++ b/python/demo/demo_mixed-topology.py
@@ -110,8 +110,8 @@
 V_cpp = _cpp.fem.FunctionSpace_float64(mesh, elements_cpp, dofmaps)
 
 # Create forms for each cell type.
-# FIXME This hack is required at the moment because UFL does not yet know about
-# mixed topology meshes.
+# FIXME This hack is required at the moment because UFL does not yet know
+# about mixed topology meshes.
 a = []
 L = []
 for i, cell_name in enumerate(["hexahedron", "prism"]):
diff --git a/python/demo/demo_navier-stokes.py b/python/demo/demo_navier-stokes.py
index 7424528750b..f25f7b268d6 100644
--- a/python/demo/demo_navier-stokes.py
+++ b/python/demo/demo_navier-stokes.py
@@ -8,7 +8,7 @@
 #       jupytext_version: 1.13.6
 # ---
 
-# # Divergence conforming discontinuous Galerkin method for the Navier--Stokes equations
+# # Divergence conforming discontinuous Galerkin method for the Navier--Stokes equations # noqa
 #
 # This demo ({download}`demo_navier-stokes.py`) illustrates how to
 # implement a divergence conforming discontinuous Galerkin method for
@@ -29,7 +29,8 @@
 #
 # $$
 # \begin{align}
-# \partial_t u - \nu \Delta u + (u \cdot \nabla)u + \nabla p &= f \text{ in } \Omega_t, \\
+# \partial_t u - \nu \Delta u + (u \cdot \nabla)u + \nabla p &= f
+# \text{ in } \Omega_t, \\
 # \nabla \cdot u &= 0 \text{ in } \Omega_t,
 # \end{align}
 # $$
@@ -124,7 +125,8 @@
 #     \end{cases}
 # $$
 #
-# where $\Gamma^\psi = \left\{x \in \Gamma; \; \psi(x) \cdot n(x) < 0\right\}$.
+# where $\Gamma^\psi = \left\{x \in \Gamma; \; \psi(x) \cdot n(x) < 0
+# \right\}$.
 #
 # The semi-discrete version problem (in dimensionless form) is: find
 # $(u_h, p_h) \in V_h^{u_D} \times Q_h$ such that
@@ -132,8 +134,8 @@
 # $$
 # \begin{align}
 #   \int_\Omega \partial_t u_h \cdot v + a_h(u_h, v_h) + c_h(u_h; u_h, v_h)
-#   + b_h(v_h, p_h) &= \int_\Omega f \cdot v_h + L_{a_h}(v_h) + L_{c_h}(v_h)
-#   \quad \forall v_h \in V_h^0, \\
+#   + b_h(v_h, p_h) &= \int_\Omega f \cdot v_h + L_{a_h}(v_h) +
+#   L_{c_h}(v_h) \quad \forall v_h \in V_h^0, \\
 #   b_h(u_h, q_h) &= 0 \quad \forall q_h \in Q_h,
 # \end{align}
 # $$
@@ -151,8 +153,10 @@
 #   c_h(w; u, v) &= - \sumK \int_K u \cdot \nabla \cdot (v \otimes w)
 #   + \sumK \int_{\partial_K} w \cdot n \hat{u}^{w} \cdot v, \\
 #   L_{a_h}(v_h) &= Re^{-1} \left(- \int_{\partial \Omega} u_D \otimes n :
-#   \nabla_h v_h + \frac{\alpha}{h} u_D \otimes n : v_h \otimes n \right), \\
-#   L_{c_h}(v_h) &= - \int_{\partial \Omega} u_D \cdot n \hat{u}_D \cdot v_h, \\
+#   \nabla_h v_h + \frac{\alpha}{h} u_D \otimes n : v_h \otimes n \right),
+#   \\
+#   L_{c_h}(v_h) &= - \int_{\partial \Omega} u_D \cdot n \hat{u}_D \cdot
+#   v_h, \\
 #   b_h(v, q) &= - \int_K \nabla \cdot v q.
 # \end{align}
 # $$
diff --git a/python/demo/demo_pml.py b/python/demo/demo_pml.py
index 1a879837250..4703662b814 100644
--- a/python/demo/demo_pml.py
+++ b/python/demo/demo_pml.py
@@ -1,4 +1,4 @@
-# # Electromagnetic scattering from a wire with perfectly matched layer condition
+# # Electromagnetic scattering from a wire with perfectly matched layer condition # noqa
 #
 # Copyright (C) 2022 Michele Castriotta, Igor Baratta, Jørgen S. Dokken
 #
@@ -221,7 +221,8 @@ def generate_mesh_wire(
 # being hit normally by a TM-polarized electromagnetic wave.
 #
 # The formula are taken from:
-# Milton Kerker, "The Scattering of Light and Other Electromagnetic Radiation",
+# Milton Kerker, "The Scattering of Light and Other Electromagnetic
+# Radiation",
 # Chapter 6, Elsevier, 1969.
 #
 # ## Implementation
@@ -489,11 +490,12 @@ def pml_coordinates(x: ufl.indexed.Indexed, alpha: float, k0: complex, l_dom: fl
 #
 # $$
 # \begin{align}
-# \text{PML}_\text{corners} \rightarrow \mathbf{r}^\prime &= (x^\prime, y^\prime) \\
+# \text{PML}_\text{corners} \rightarrow \mathbf{r}^\prime &= (x^\prime,
+#   y^\prime) \\
 # \text{PML}_\text{rectangles along x} \rightarrow
-#                                       \mathbf{r}^\prime &= (x^\prime, y) \\
+#   \mathbf{r}^\prime &= (x^\prime, y) \\
 # \text{PML}_\text{rectangles along y} \rightarrow
-#                                       \mathbf{r}^\prime &= (x, y^\prime).
+#   \mathbf{r}^\prime &= (x, y^\prime).
 # \end{align}
 # $$
 #
@@ -662,8 +664,8 @@ def create_eps_mu(
 # The final weak form in the PML region is:
 #
 # $$
-# \int_{\Omega_{pml}}\left[\boldsymbol{\mu}^{-1}_{pml} \nabla \times \mathbf{E}
-# \right]\cdot \nabla \times \bar{\mathbf{v}}-k_{0}^{2}
+# \int_{\Omega_{pml}}\left[\boldsymbol{\mu}^{-1}_{pml} \nabla \times
+# \mathbf{E} \right]\cdot \nabla \times \bar{\mathbf{v}}-k_{0}^{2}
 # \left[\boldsymbol{\varepsilon}_{pml} \mathbf{E} \right]\cdot
 # \bar{\mathbf{v}}~ d x=0,
 # $$
diff --git a/python/demo/demo_poisson.py b/python/demo/demo_poisson.py
index 0ee822c719d..ae6bd040dc8 100644
--- a/python/demo/demo_poisson.py
+++ b/python/demo/demo_poisson.py
@@ -43,7 +43,8 @@
 # $$
 # \begin{align}
 #   a(u, v) &:= \int_{\Omega} \nabla u \cdot \nabla v \, {\rm d} x, \\
-#   L(v)    &:= \int_{\Omega} f v \, {\rm d} x + \int_{\Gamma_{N}} g v \, {\rm d} s.
+#   L(v)    &:= \int_{\Omega} f v \, {\rm d} x + \int_{\Gamma_{N}} g v \,
+#               {\rm d} s.
 # \end{align}
 # $$
 #
diff --git a/python/demo/demo_poisson_matrix_free.py b/python/demo/demo_poisson_matrix_free.py
index fa5fe103bdd..194a335bc9d 100644
--- a/python/demo/demo_poisson_matrix_free.py
+++ b/python/demo/demo_poisson_matrix_free.py
@@ -107,9 +107,10 @@
 #
 # Next, we locate the mesh facets that lie on the domain boundary
 # $\partial\Omega$. We do this by first calling
-# {py:func}`create_connectivity <dolfinx.mesh.topology.create_connectivity>`
-# and then retrieving all facets on the boundary using
-# {py:func}`exterior_facet_indices <dolfinx.mesh.exterior_facet_indices>`.
+# {py:func}`create_connectivity
+# <dolfinx.mesh.topology.create_connectivity>`  and then retrieving all
+# facets on the boundary using {py:func}`exterior_facet_indices
+# <dolfinx.mesh.exterior_facet_indices>`.
 
 tdim = mesh.topology.dim
 mesh.topology.create_connectivity(tdim - 1, tdim)
@@ -156,10 +157,11 @@
 
 # ### Matrix-free conjugate gradient solver
 #
-# The right hand side vector $b - A x_{\rm bc}$ is the assembly of the linear
-# form $L$ where the essential Dirichlet boundary conditions are implemented
-# using lifting. Since we want to avoid assembling the matrix `A`, we compute
-# the necessary matrix-vector product using the linear form `M` explicitly.
+# The right hand side vector $b - A x_{\rm bc}$ is the assembly of the
+# linear form $L$ where the essential Dirichlet boundary conditions are
+# implemented using lifting. Since we want to avoid assembling the matrix
+# `A`, we compute the necessary matrix-vector product using the linear form
+# `M` explicitly.
 
 # Apply lifting: b <- b - A * x_bc
 b = fem.assemble_vector(L_fem)
diff --git a/python/demo/demo_pyamg.py b/python/demo/demo_pyamg.py
index dab612f9a1f..3a7f057eab7 100644
--- a/python/demo/demo_pyamg.py
+++ b/python/demo/demo_pyamg.py
@@ -58,7 +58,8 @@ def poisson_problem(dtype: npt.DTypeLike, solver_type: str):
 
     Args:
         dtype: Scalar type to use.
-        solver_type: pyamg solver type, either "ruge_stuben" or "smoothed_aggregation"
+        solver_type: pyamg solver type, either "ruge_stuben" or
+            "smoothed_aggregation"
     """
 
     real_type = np.real(dtype(0)).dtype
@@ -191,7 +192,8 @@ def elasticity_problem(dtype):
     λ = E * ν / ((1.0 + ν) * (1.0 - 2.0 * ν))
 
     def σ(v):
-        """Return an expression for the stress σ given a displacement field"""
+        """Return an expression for the stress σ given a displacement
+        field"""
         return 2.0 * μ * ufl.sym(ufl.grad(v)) + λ * ufl.tr(ufl.sym(ufl.grad(v))) * ufl.Identity(
             len(v)
         )
diff --git a/python/demo/demo_pyvista.py b/python/demo/demo_pyvista.py
index c1b0cd81dc8..fdf2b3d9bbf 100644
--- a/python/demo/demo_pyvista.py
+++ b/python/demo/demo_pyvista.py
@@ -256,8 +256,8 @@ def plot_nedelec():
     # Add this grid (as a wireframe) to the plotter
     plotter.add_mesh(grid, style="wireframe", line_width=2, color="black")
 
-    # Create a function space consisting of first order Nédélec (first kind)
-    # elements and interpolate a vector-valued expression
+    # Create a function space consisting of first order Nédélec (first
+    # kind) elements and interpolate a vector-valued expression
     V = functionspace(msh, ("N1curl", 2))
     u = Function(V, dtype=np.float64)
     u.interpolate(lambda x: (x[2] ** 2, np.zeros(x.shape[1]), -x[0] * x[2]))
@@ -276,7 +276,8 @@ def plot_nedelec():
     cells, cell_types, x = plot.vtk_mesh(V0)
     grid = pyvista.UnstructuredGrid(cells, cell_types, x)
 
-    # Create point cloud of vertices, and add the vertex values to the cloud
+    # Create point cloud of vertices, and add the vertex values to the
+    # cloud
     grid.point_data["u"] = u0.x.array.reshape(x.shape[0], V0.dofmap.index_map_bs)
     glyphs = grid.glyph(orient="u", factor=0.1)
 
diff --git a/python/demo/demo_scattering_boundary_conditions.py b/python/demo/demo_scattering_boundary_conditions.py
index 5332e45cb74..f3f5baf5c53 100644
--- a/python/demo/demo_scattering_boundary_conditions.py
+++ b/python/demo/demo_scattering_boundary_conditions.py
@@ -13,7 +13,7 @@
 #     name: python3-complex
 # ---
 
-# # Electromagnetic scattering from a wire with scattering boundary conditions
+# # Electromagnetic scattering from a wire with scattering boundary conditions # noqa
 #
 # Copyright (C) 2022 Michele Castriotta, Igor Baratta, Jørgen S. Dokken
 #
@@ -163,8 +163,8 @@ def generate_mesh_wire(
 # being hit normally by a TM-polarized electromagnetic wave.
 #
 # The formula are taken from:
-# Milton Kerker, "The Scattering of Light and Other Electromagnetic Radiation",
-# Chapter 6, Elsevier, 1969.
+# Milton Kerker, "The Scattering of Light and Other Electromagnetic
+# Radiation", Chapter 6, Elsevier, 1969.
 #
 # ## Implementation
 # First of all, let's define the parameters of the problem:
@@ -562,8 +562,8 @@ def curl_2d(f: fem.Function):
 # & \int_{\Omega}-\nabla \cdot(\nabla\times\mathbf{E}_s \times
 # \bar{\mathbf{v}})-\nabla \times \mathbf{E}_s \cdot \nabla
 # \times\bar{\mathbf{v}}+\varepsilon_{r} k_{0}^{2} \mathbf{E}_s
-# \cdot \bar{\mathbf{v}}+k_{0}^{2}\left(\varepsilon_{r}-\varepsilon_b\right)
-# \mathbf{E}_b \cdot \bar{\mathbf{v}}~\mathrm{dx} \\
+# \cdot \bar{\mathbf{v}}+k_{0}^{2}\left(\varepsilon_{r}-\varepsilon_b
+# \right) \mathbf{E}_b \cdot \bar{\mathbf{v}}~\mathrm{dx} \\
 # +&\int_{\partial \Omega}
 # (\mathbf{n} \times \nabla \times \mathbf{E}_s) \cdot \bar{\mathbf{v}}
 # +\left(j n_bk_{0}+\frac{1}{2r}\right) (\mathbf{n} \times \mathbf{E}_s
@@ -591,10 +591,11 @@ def curl_2d(f: fem.Function):
 #
 # Cancelling $-(\nabla\times\mathbf{E}_s \times \bar{\mathbf{V}})
 # \cdot\mathbf{n}$  and $\mathbf{n} \times \nabla \times \mathbf{E}_s
-# \cdot \bar{\mathbf{V}}$ and rearrange $\left((\mathbf{n} \times \mathbf{E}_s)
-# \times \mathbf{n}\right) \cdot \bar{\mathbf{v}}$ to $ (\mathbf{E}_s \times\mathbf{n})
-# \cdot (\bar{\mathbf{v}} \times \mathbf{n})$ using the triple product rule $\mathbf{A}
-# \cdot(\mathbf{B} \times \mathbf{C})=\mathbf{B} \cdot(\mathbf{C} \times
+# \cdot \bar{\mathbf{V}}$ and rearrange $\left((\mathbf{n} \times
+# \mathbf{E}_s) \times \mathbf{n}\right) \cdot \bar{\mathbf{v}}$ to
+# $(\mathbf{E}_s \times\mathbf{n}) \cdot (\bar{\mathbf{v}} \times
+# \mathbf{n})$ using the triple product rule $\mathbf{A} \cdot(\mathbf{B}
+# \times \mathbf{C})=\mathbf{B} \cdot(\mathbf{C} \times
 # \mathbf{A})=\mathbf{C} \cdot(\mathbf{A} \times \mathbf{B})$, we get:
 #
 # $$
@@ -604,8 +605,8 @@ def curl_2d(f: fem.Function):
 # \bar{\mathbf{v}}+k_{0}^{2}\left(\varepsilon_{r}-\varepsilon_b\right)
 # \mathbf{E}_b \cdot \bar{\mathbf{v}}~\mathrm{d}x \\
 # +&\int_{\partial \Omega}
-# \left(j n_bk_{0}+\frac{1}{2r}\right)( \mathbf{n} \times \mathbf{E}_s \times
-# \mathbf{n}) \cdot \bar{\mathbf{v}} ~\mathrm{d} s = 0.
+# \left(j n_bk_{0}+\frac{1}{2r}\right)( \mathbf{n} \times \mathbf{E}_s
+# \times \mathbf{n}) \cdot \bar{\mathbf{v}} ~\mathrm{d} s = 0.
 # \end{align}
 # $$
 #
diff --git a/python/demo/demo_static-condensation.py b/python/demo/demo_static-condensation.py
index 70dd5d52f36..53114a8f6df 100644
--- a/python/demo/demo_static-condensation.py
+++ b/python/demo/demo_static-condensation.py
@@ -107,7 +107,8 @@
 
 
 def sigma_u(u):
-    """Constitutive relation for stress-strain. Assuming plane-stress in XY"""
+    """Constitutive relation for stress-strain. Assuming plane-stress in
+    XY"""
     eps = 0.5 * (ufl.grad(u) + ufl.grad(u).T)
     sigma = E / (1.0 - nu**2) * ((1.0 - nu) * eps + nu * ufl.Identity(2) * ufl.tr(eps))
     return sigma
diff --git a/python/demo/demo_stokes.py b/python/demo/demo_stokes.py
index 018fd83e486..eec543f63a4 100644
--- a/python/demo/demo_stokes.py
+++ b/python/demo/demo_stokes.py
@@ -75,7 +75,8 @@
 #
 # The Stokes problem using Taylor-Hood elements is solved using:
 # 1. [Block preconditioner using PETSc MatNest and VecNest data
-#    structures. Each 'block' is a standalone object.](#nested-matrix-solver)
+#    structures. Each 'block' is a standalone object.]
+#    (#nested-matrix-solver)
 # 1. [Block preconditioner with the `u` and `p` fields stored block-wise
 #    in a single matrix](#monolithic-block-iterative-solver)
 # 1. [Direct solver with the `u` and `p` fields stored block-wise in a
@@ -210,7 +211,8 @@ def lid_velocity_expression(x):
 
 
 def nested_iterative_solver():
-    """Solve the Stokes problem using nest matrices and an iterative solver."""
+    """Solve the Stokes problem using nest matrices and an iterative
+    solver."""
 
     # Assemble nested matrix operators
     A = fem.petsc.assemble_matrix(a, bcs=bcs, kind="nest")
@@ -331,7 +333,8 @@ def nested_iterative_solver():
 
 
 def block_operators():
-    """Return block operators and block RHS vector for the Stokes problem."""
+    """Return block operators and block RHS vector for the Stokes
+    problem."""
 
     # Assembler matrix operator, preconditioner and RHS vector into
     # single objects but preserving block structure
diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index 0f98298310d..d2b5be272a8 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -74,10 +74,10 @@
 # For elements where the coefficients matrix is not an identity, we can
 # use the properties of orthonormal polynomials to compute `wcoeffs`.
 # Let $\{q_0, q_1,\dots\}$ be the orthonormal polynomials of a given
-# degree for a given cell, and suppose that we're trying to represent a function
-# $f_i\in\operatorname{span}\{q_1, q_2,\dots\}$ (as $\{f_0, f_1,\dots\}$ is a
-# basis of the polynomial space for our element). Using the properties of
-# orthonormal polynomials, we see that
+# degree for a given cell, and suppose that we're trying to represent a
+# function $f_i\in\operatorname{span}\{q_1, q_2,\dots\}$ (as
+# $\{f_0, f_1,\dots\}$ is a basis of the polynomial space for our element).
+# Using the properties of orthonormal polynomials, we see that
 # $f_i = \sum_j\left(\int_R f_iq_j\,\mathrm{d}\mathbf{x}\right)q_j$,
 # and so the coefficients are given by
 # $a_{ij}=\int_R f_iq_j\,\mathrm{d}\mathbf{x}$.
diff --git a/python/demo/demo_types.py b/python/demo/demo_types.py
index 7d647ecfde9..a801028a38d 100644
--- a/python/demo/demo_types.py
+++ b/python/demo/demo_types.py
@@ -182,7 +182,8 @@ def elasticity(dtype) -> fem.Function:
     μ, λ = E / (2.0 * (1.0 + ν)), E * ν / ((1.0 + ν) * (1.0 - 2.0 * ν))
 
     def σ(v):
-        """Return an expression for the stress σ given a displacement field"""
+        """Return an expression for the stress σ given a displacement
+        field"""
         return 2.0 * μ * ufl.sym(ufl.grad(v)) + λ * ufl.tr(ufl.sym(ufl.grad(v))) * ufl.Identity(
             len(v)
         )
diff --git a/python/demo/pyproject.toml b/python/demo/pyproject.toml
index 28e7749be23..55d86c22d7d 100644
--- a/python/demo/pyproject.toml
+++ b/python/demo/pyproject.toml
@@ -3,3 +3,6 @@ extend = "../pyproject.toml"
 
 [tool.ruff.lint.flake8-import-conventions]
 banned-from = ["ufl"]
+
+[tool.ruff.lint.pycodestyle]
+max-doc-length = 75
diff --git a/python/dolfinx/fem/__init__.py b/python/dolfinx/fem/__init__.py
index a28892ec1cd..f08516eb093 100644
--- a/python/dolfinx/fem/__init__.py
+++ b/python/dolfinx/fem/__init__.py
@@ -99,7 +99,8 @@ def create_interpolation_data(
     cells: npt.NDArray[np.int32],
     padding: float = 1e-14,
 ) -> _PointOwnershipData:
-    """Generate data needed to interpolate discrete functions across different meshes.
+    """Generate data needed to interpolate discrete functions across
+    different meshes.
 
     Args:
         V_to: Function space to interpolate into
@@ -159,7 +160,8 @@ def discrete_gradient(space0: FunctionSpace, space1: FunctionSpace) -> _MatrixCS
 
 
 def interpolation_matrix(space0: FunctionSpace, space1: FunctionSpace) -> _MatrixCSR:
-    """Assemble an interpolation matrix for two function spaces on the same mesh.
+    """Assemble an interpolation matrix for two function spaces on the same
+    mesh.
 
     Args:
         space0: space to interpolate from
@@ -174,7 +176,8 @@ def interpolation_matrix(space0: FunctionSpace, space1: FunctionSpace) -> _Matri
 def compute_integration_domains(
     integral_type: IntegralType, topology: Topology, entities: np.ndarray
 ):
-    """Given an integral type and a set of entities compute integration entities.
+    """Given an integral type and a set of entities compute integration
+    entities.
 
     This function returns a list ``[(id, entities)]``. For cell
     integrals ``entities`` are the cell indices. For exterior facet
diff --git a/python/dolfinx/fem/assemble.py b/python/dolfinx/fem/assemble.py
index 4aa9d47b17f..c7662c66ee0 100644
--- a/python/dolfinx/fem/assemble.py
+++ b/python/dolfinx/fem/assemble.py
@@ -92,7 +92,7 @@ def _pack(form):
     return _pack(form)
 
 
-# -- Vector and matrix instantiation -----------------------------------------
+# -- Vector and matrix instantiation --------------------------------------
 
 
 def create_vector(L: Form) -> la.Vector:
@@ -111,7 +111,8 @@ def create_vector(L: Form) -> la.Vector:
 
 
 def create_matrix(a: Form, block_mode: typing.Optional[la.BlockMode] = None) -> la.MatrixCSR:
-    """Create a sparse matrix that is compatible with a given bilinear form.
+    """Create a sparse matrix that is compatible with a given bilinear
+    form.
 
     Args:
         a: Bilinear form.
@@ -129,7 +130,7 @@ def create_matrix(a: Form, block_mode: typing.Optional[la.BlockMode] = None) ->
         return la.matrix_csr(sp, dtype=a.dtype)
 
 
-# -- Scalar assembly ---------------------------------------------------------
+# -- Scalar assembly ------------------------------------------------------
 
 
 def assemble_scalar(M: Form, constants: typing.Optional[np.ndarray] = None, coeffs=None):
@@ -159,7 +160,7 @@ def assemble_scalar(M: Form, constants: typing.Optional[np.ndarray] = None, coef
     return _cpp.fem.assemble_scalar(M._cpp_object, constants, coeffs)
 
 
-# -- Vector assembly ---------------------------------------------------------
+# -- Vector assembly ------------------------------------------------------
 
 
 @functools.singledispatch
@@ -234,7 +235,7 @@ def _assemble_vector_array(
     return b
 
 
-# -- Matrix assembly ---------------------------------------------------------
+# -- Matrix assembly ------------------------------------------------------
 
 
 @functools.singledispatch
@@ -321,7 +322,7 @@ def _assemble_matrix_csr(
     return A
 
 
-# -- Modifiers for Dirichlet conditions ---------------------------------------
+# -- Modifiers for Dirichlet conditions -----------------------------------
 
 
 def apply_lifting(
@@ -333,7 +334,8 @@ def apply_lifting(
     constants: typing.Optional[Iterable[np.ndarray]] = None,
     coeffs=None,
 ) -> None:
-    """Modify right-hand side vector ``b`` for lifting of Dirichlet boundary conditions.
+    """Modify right-hand side vector ``b`` for lifting of Dirichlet
+    boundary conditions.
 
     Consider the discrete algebraic system:
 
@@ -349,9 +351,12 @@ def apply_lifting(
 
     .. math::
 
-        \\begin{bmatrix} A_{0}^{(0)} & A_{0}^{(1)} & A_{1}^{(0)} & A_{1}^{(1)}\\end{bmatrix}
-        \\begin{bmatrix}u_{0}^{(0)} \\\\ u_{0}^{(1)}
-        \\\\ u_{1}^{(0)} \\\\ u_{1}^{(1)}\\end{bmatrix}
+        \\begin{bmatrix}
+            A_{0}^{(0)} & A_{0}^{(1)} & A_{1}^{(0)} & A_{1}^{(1)}
+        \\end{bmatrix}
+        \\begin{bmatrix}
+            u_{0}^{(0)} \\\\ u_{0}^{(1)} \\\\ u_{1}^{(0)} \\\\ u_{1}^{(1)}
+        \\end{bmatrix}
         = b.
 
     If :math:`u_{i}^{(1)} = \\alpha(g_{i} - x_{i})`, where :math:`g_{i}`
@@ -360,7 +365,9 @@ def apply_lifting(
 
     .. math::
 
-        \\begin{bmatrix}A_{0}^{(0)} & A_{0}^{(1)} & A_{1}^{(0)} & A_{1}^{(1)} \\end{bmatrix}
+        \\begin{bmatrix}
+            A_{0}^{(0)} & A_{0}^{(1)} & A_{1}^{(0)} & A_{1}^{(1)}
+        \\end{bmatrix}
         \\begin{bmatrix}u_{0}^{(0)} \\\\ \\alpha(g_{0} - x_{0})
         \\\\ u_{1}^{(0)} \\\\ \\alpha(g_{1} - x_{1})\\end{bmatrix}
         = b.
@@ -371,7 +378,8 @@ def apply_lifting(
 
         \\begin{bmatrix}A_{0}^{(0)} & A_{1}^{(0)}\\end{bmatrix}
         \\begin{bmatrix}u_{0}^{(0)} \\\\ u_{1}^{(0)}\\end{bmatrix}
-        = b - \\alpha A_{0}^{(1)} (g_{0} - x_{0}) - \\alpha A_{1}^{(1)} (g_{1} - x_{1}).
+        = b - \\alpha A_{0}^{(1)} (g_{0} - x_{0})
+        - \\alpha A_{1}^{(1)} (g_{1} - x_{1}).
 
     The modified  :math:`b` vector is
 
@@ -397,7 +405,10 @@ def apply_lifting(
             list-of-lists of `DirichletBC` from a list of forms ``a``
             and a flat list of `DirichletBC` objects ``bcs``::
 
-                bcs1 = fem.bcs_by_block(fem.extract_function_spaces([a], 1), bcs)
+                bcs1 = fem.bcs_by_block(
+                    fem.extract_function_spaces([a], 1),
+                    bcs
+                )
 
         x0: The array :math:`x_{i}` above. If ``None`` it is set to
             zero.
diff --git a/python/dolfinx/fem/bcs.py b/python/dolfinx/fem/bcs.py
index e433bc81d95..d5656204da5 100644
--- a/python/dolfinx/fem/bcs.py
+++ b/python/dolfinx/fem/bcs.py
@@ -1,4 +1,5 @@
-# Copyright (C) 2017-2021 Chris N. Richardson, Garth N. Wells and Jørgen S. Dokken
+# Copyright (C) 2017-2021 Chris N. Richardson, Garth N. Wells and
+# Jørgen S. Dokken
 #
 # This file is part of DOLFINx (https://www.fenicsproject.org)
 #
@@ -99,8 +100,8 @@ class DirichletBC:
     ]
 
     def __init__(self, bc):
-        """Representation of Dirichlet boundary condition which is imposed on
-        a linear system.
+        """Representation of Dirichlet boundary condition which is imposed
+        on a linear system.
 
         Note:
             Dirichlet boundary conditions  should normally be
@@ -135,7 +136,8 @@ def function_space(self) -> dolfinx.fem.FunctionSpace:
     def set(
         self, x: npt.NDArray, x0: typing.Optional[npt.NDArray[np.int32]] = None, alpha: float = 1
     ) -> None:
-        """Set entries in an array that are constrained by Dirichlet boundary conditions.
+        """Set entries in an array that are constrained by Dirichlet
+        boundary conditions.
 
         Entries in ``x`` that are constrained by a Dirichlet boundary
         conditions are set to ``alpha * (x_bc - x0)``, where ``x_bc`` is
@@ -160,9 +162,9 @@ def set(
         self._cpp_object.set(x, x0, alpha)
 
     def dof_indices(self) -> tuple[npt.NDArray[np.int32], int]:
-        """Access dof indices `(local indices, unrolled)`, including ghosts, to
-        which a Dirichlet condition is applied, and the index to the first
-        non-owned (ghost) index. The array of indices is sorted.
+        """Access dof indices `(local indices, unrolled)`, including
+        ghosts, to which a Dirichlet condition is applied, and the index to
+        the first non-owned (ghost) index. The array of indices is sorted.
 
         Note:
             The returned array is read-only.
diff --git a/python/dolfinx/fem/dofmap.py b/python/dolfinx/fem/dofmap.py
index b8a848170dc..d8421032cfc 100644
--- a/python/dolfinx/fem/dofmap.py
+++ b/python/dolfinx/fem/dofmap.py
@@ -43,7 +43,8 @@ def dof_layout(self):
 
     @property
     def index_map(self):
-        """Index map that described the parallel distribution of the dofmap."""
+        """Index map that described the parallel distribution of the
+        dofmap."""
         return self._cpp_object.index_map
 
     @property
diff --git a/python/dolfinx/fem/element.py b/python/dolfinx/fem/element.py
index 178d19f925b..78a36d207f7 100644
--- a/python/dolfinx/fem/element.py
+++ b/python/dolfinx/fem/element.py
@@ -98,8 +98,9 @@ def pull_back(
             x: Physical coordinates to pull back to the reference cells,
                 ``shape=(num_points, geometrical_dimension)``.
             cell_geometry: Physical coordinates describing the cell,
-                shape ``(num_of_geometry_basis_functions, geometrical_dimension)``
-                They can be created by accessing ``geometry.x[geometry.dofmap.cell_dofs(i)]``,
+                shape ``(num_of_geometry_basis_functions,
+                geometrical_dimension)``. They can be created by accessing
+                ``geometry.x[geometry.dofmap.cell_dofs(i)]``,
 
         Returns:
             Reference coordinates of the physical points ``x``.
@@ -111,8 +112,8 @@ def variant(self) -> int:
         """Lagrange variant of the coordinate element.
 
         Note:
-            The return type is an integer. A Basix enum can be created using
-            ``basix.LagrangeVariant(value)``.
+            The return type is an integer. A Basix enum can be
+            created using ``basix.LagrangeVariant(value)``.
         """
         return self._cpp_object.variant
 
@@ -178,7 +179,8 @@ def __init__(
         """Creates a Python wrapper for the exported finite element class.
 
         Note:
-            Do not use this constructor directly. Instead use :func:``finiteelement``.
+            Do not use this constructor directly. Instead use
+            :func:``finiteelement``.
 
         Args:
             The underlying cpp instance that this object will wrap.
@@ -190,7 +192,8 @@ def __eq__(self, other):
 
     @property
     def dtype(self) -> np.dtype:
-        """Geometry type of the Mesh that the FunctionSpace is defined on."""
+        """Geometry type of the Mesh that the FunctionSpace is defined
+        on."""
         return self._cpp_object.dtype
 
     @property
@@ -211,40 +214,45 @@ def num_sub_elements(self) -> int:
     def value_shape(self) -> npt.NDArray[np.integer]:
         """Value shape of the finite element field.
 
-        The value shape describes the shape of the finite element field, e.g. ``{}`` for a scalar,
-        ``{2}`` for a vector in 2D, ``{3, 3}`` for a rank-2 tensor in 3D, etc.
+        The value shape describes the shape of the finite element field,
+        e.g. ``{}`` for a scalar, ``{2}`` for a vector in 2D, ``{3, 3}``
+        for a rank-2 tensor in 3D, etc.
         """
         return self._cpp_object.value_shape
 
     @property
     def interpolation_points(self) -> npt.NDArray[np.floating]:
-        """Points on the reference cell at which an expression needs to be evaluated in order to
-        interpolate the expression in the finite element space.
+        """Points on the reference cell at which an expression needs to be
+        evaluated in order to interpolate the expression in the finite
+        element space.
 
-        Interpolation point coordinates on the reference cell, returning the coordinates data
-        (row-major) storage with shape ``(num_points, tdim)``.
+        Interpolation point coordinates on the reference cell, returning
+        the coordinates data (row-major) storage with shape
+        ``(num_points, tdim)``.
 
         Note:
-            For Lagrange elements the points will just be the nodal positions. For other elements
-            the points will typically be the quadrature points used to evaluate moment degrees of
-            freedom.
+            For Lagrange elements the points will just be the nodal
+            positions. For other elements the points will typically be the
+            quadrature points used to evaluate moment degrees of freedom.
         """
         return self._cpp_object.interpolation_points()
 
     @property
     def interpolation_ident(self) -> bool:
-        """Check if interpolation into the finite element space is an identity operation given the
-        evaluation on an expression at specific points, i.e. the degree-of-freedom are equal to
-        point evaluations. The function will return `true` for Lagrange elements."""
+        """Check if interpolation into the finite element space is an
+        identity operation given the evaluation on an expression at
+        specific points, i.e. the degree-of-freedom are equal to point
+        evaluations. The function will return `true` for Lagrange
+        elements."""
         return self._cpp_object.interpolation_ident
 
     @property
     def space_dimension(self) -> int:
-        """Dimension of the finite element function space (the number of degrees-of-freedom for the
-        element).
+        """Dimension of the finite element function space (the number of
+        degrees-of-freedom for the element).
 
-        For 'blocked' elements, this function returns the dimension of the full element rather than
-        the dimension of the base element.
+        For 'blocked' elements, this function returns the dimension of the
+        full element rather than the dimension of the base element.
         """
         return self._cpp_object.space_dimension
 
@@ -252,13 +260,15 @@ def space_dimension(self) -> int:
     def needs_dof_transformations(self) -> bool:
         """Check if DOF transformations are needed for this element.
 
-        DOF transformations will be needed for elements which might not be continuous when two
-        neighbouring cells disagree on the orientation of a shared sub-entity, and when this cannot
-        be corrected for by permuting the DOF numbering in the dofmap.
+        DOF transformations will be needed for elements which might not be
+        continuous when two neighbouring cells disagree on the orientation
+        of a shared sub-entity, and when this cannot be corrected for by
+        permuting the DOF numbering in the dofmap.
 
-        For example, Raviart-Thomas elements will need DOF transformations, as the neighbouring
-        cells may disagree on the orientation of a basis function, and this orientation cannot be
-        corrected for by permuting the DOF numbers on each cell.
+        For example, Raviart-Thomas elements will need DOF transformations,
+        as the neighbouring cells may disagree on the orientation of a
+        basis function, and this orientation cannot be corrected for by
+        permuting the DOF numbers on each cell.
         """
         return self._cpp_object.needs_dof_transformations
 
@@ -270,18 +280,21 @@ def signature(self) -> str:
     def T_apply(
         self, x: npt.NDArray[np.floating], cell_permutations: npt.NDArray[np.uint32], dim: int
     ) -> None:
-        """Transform basis functions from the reference element ordering and orientation to the
-        globally consistent physical element ordering and orientation.
+        """Transform basis functions from the reference element ordering
+        and orientation to the globally consistent physical element
+        ordering and orientation.
 
         Args:
-            x: Data to transform (in place). The shape is ``(num_cells, n, dim)``, where ``n`` is
-                the number degrees-of-freedom and the data is flattened (row-major).
+            x: Data to transform (in place). The shape is
+                ``(num_cells, n, dim)``, where ``n`` is the number degrees-
+                of-freedom and the data is flattened (row-major).
             cell_permutations: Permutation data for the cell.
             dim: Number of columns in ``data``.
 
         Note:
-            Exposed for testing. Function is not vectorised across multiple cells. Please see
-            `basix.numba_helpers` for performant versions.
+            Exposed for testing. Function is not vectorised across multiple
+            cells. Please see `basix.numba_helpers` for performant
+            versions.
         """
         self._cpp_object.T_apply(x, cell_permutations, dim)
 
@@ -291,8 +304,9 @@ def Tt_apply(
         """Apply the transpose of the operator applied by T_apply().
 
         Args:
-            x: Data to transform (in place). The shape is ``(num_cells, n, dim)``, where ``n`` is
-                the number degrees-of-freedom and the data is flattened (row-major).
+            x: Data to transform (in place). The shape is
+                ``(num_cells, n, dim)``, where ``n`` is the number degrees-
+                of-freedom and the data is flattened (row-major).
             cell_permutations: Permutation data for the cells
             dim: Number of columns in ``data``.
         """
@@ -301,11 +315,13 @@ def Tt_apply(
     def Tt_inv_apply(
         self, x: npt.NDArray[np.floating], cell_permutations: npt.NDArray[np.uint32], dim: int
     ) -> None:
-        """Apply the inverse transpose of the operator applied by T_apply().
+        """Apply the inverse transpose of the operator applied by
+        T_apply().
 
         Args:
-            x: Data to transform (in place). The shape is ``(num_cells, n, dim)``, where ``n`` is
-                the number degrees-of-freedom and the data is flattened (row-major).
+            x: Data to transform (in place). The shape is
+                ``(num_cells, n, dim)``, where ``n`` is the number degrees-
+                of-freedom and the data is flattened (row-major).
             cell_permutations: Permutation data for the cells
             dim: Number of columns in ``data``.
         """
@@ -321,7 +337,8 @@ def finiteelement(
 
     Args:
         cell_type: Element cell type, see ``mesh.CellType``
-        ufl_e: UFL element, holding quadrature rule and other properties of the selected element.
+        ufl_e: UFL element, holding quadrature rule and other properties of
+            the selected element.
         FiniteElement_dtype: Geometry type of the element.
     """
     if np.issubdtype(FiniteElement_dtype, np.float32):
diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py
index 750eaacf245..6db45890cc2 100644
--- a/python/dolfinx/fem/forms.py
+++ b/python/dolfinx/fem/forms.py
@@ -165,7 +165,8 @@ def form_cpp_class(
 ) -> typing.Union[
     _cpp.fem.Form_float32, _cpp.fem.Form_float64, _cpp.fem.Form_complex64, _cpp.fem.Form_complex128
 ]:
-    """Return the wrapped C++ class of a variational form of a specific scalar type.
+    """Return the wrapped C++ class of a variational form of a specific
+    scalar type.
 
     Args:
         dtype: Scalar type of the required form class.
@@ -207,9 +208,9 @@ def mixed_topology_form(
     """
     Create a mixed-topology from from an array of Forms.
 
-    # FIXME: This function is a temporary hack for mixed-topology meshes. It is needed
-    # because UFL does not know about mixed-topology meshes, so we need
-    # to pass a list of forms for each cell type.
+    # FIXME: This function is a temporary hack for mixed-topology meshes.
+    # It is needed because UFL does not know about mixed-topology meshes,
+    # so we need to pass a list of forms for each cell type.
 
     Args:
         form: A list of UFL forms. Each form should be the same, just
@@ -395,7 +396,8 @@ def _form(form):
         return Form(f, ufcx_form, code, module)
 
     def _zero_form(form):
-        """Compile a single 'zero' UFL form, i.e. a form with no integrals."""
+        """Compile a single 'zero' UFL form, i.e. a form with no
+        integrals."""
         V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]
         assert len(V) > 0
         msh = V[0].mesh
@@ -540,7 +542,8 @@ def form_cpp_creator(
     _cpp.fem.Form_complex64,
     _cpp.fem.Form_complex128,
 ]:
-    """Return the wrapped C++ constructor for creating a variational form of a specific scalar type.
+    """Return the wrapped C++ constructor for creating a variational form
+    of a specific scalar type.
 
     Args:
         dtype: Scalar type of the required form class.
diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index 0374a6fa24d..48785e6ed0d 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -1,4 +1,5 @@
-# Copyright (C) 2009-2024 Chris N. Richardson, Garth N. Wells, Michal Habera and Jørgen S. Dokken
+# Copyright (C) 2009-2024 Chris N. Richardson, Garth N. Wells,
+# Michal Habera and Jørgen S. Dokken
 #
 # This file is part of DOLFINx (https://www.fenicsproject.org)
 #
@@ -759,7 +760,8 @@ def collapse(self) -> tuple[FunctionSpace, np.ndarray]:
         return V, dofs
 
     def tabulate_dof_coordinates(self) -> npt.NDArray[np.float64]:
-        """Tabulate the coordinates of the degrees-of-freedom in the function space.
+        """Tabulate the coordinates of the degrees-of-freedom in the
+        function space.
 
         Returns:
             Coordinates of the degrees-of-freedom.
diff --git a/python/dolfinx/fem/petsc.py b/python/dolfinx/fem/petsc.py
index e98e3abb390..5de40699c71 100644
--- a/python/dolfinx/fem/petsc.py
+++ b/python/dolfinx/fem/petsc.py
@@ -102,7 +102,7 @@ def fn(form):
     return rows, cols
 
 
-# -- Vector instantiation ----------------------------------------------------
+# -- Vector instantiation -------------------------------------------------
 
 
 def create_vector(
@@ -186,14 +186,15 @@ def create_vector(
             )
 
 
-# -- Matrix instantiation ----------------------------------------------------
+# -- Matrix instantiation -------------------------------------------------
 
 
 def create_matrix(
     a: typing.Union[Form, Iterable[Iterable[Form]]],
     kind: typing.Optional[typing.Union[str, Iterable[Iterable[str]]]] = None,
 ) -> PETSc.Mat:
-    """Create a PETSc matrix that is compatible with the (sequence) of bilinear form(s).
+    """Create a PETSc matrix that is compatible with the (sequence) of
+    bilinear form(s).
 
     Three cases are supported:
 
@@ -236,7 +237,7 @@ def create_matrix(
                 return _cpp.fem.petsc.create_matrix_nest(_a, kind)  # Array of 'kind' types
 
 
-# -- Vector assembly ---------------------------------------------------------
+# -- Vector assembly ------------------------------------------------------
 
 
 @functools.singledispatch
@@ -376,7 +377,7 @@ def _assemble_vector_vec(
     return b
 
 
-# -- Matrix assembly ---------------------------------------------------------
+# -- Matrix assembly ------------------------------------------------------
 @functools.singledispatch
 def assemble_matrix(
     a: typing.Union[Form, Iterable[Iterable[Form]]],
@@ -560,7 +561,7 @@ def assemble_matrix_mat(
     return A
 
 
-# -- Modifiers for Dirichlet conditions ---------------------------------------
+# -- Modifiers for Dirichlet conditions -----------------------------------
 
 
 def apply_lifting(
@@ -596,7 +597,9 @@ def apply_lifting(
                of `DirichletBC` from the 2D ``a`` and a flat list of
                `DirichletBC` objects ``bcs``::
 
-                   bcs1 = fem.bcs_by_block(fem.extract_function_spaces(a, 1), bcs)
+                   bcs1 = fem.bcs_by_block(
+                    fem.extract_function_spaces(a, 1), bcs
+                   )
 
             2. ``bcs`` is a sequence of :class:`dolfinx.fem.DirichletBC`
                objects. The function deduces which `DiricletBC` objects
@@ -774,10 +777,11 @@ def __init__(
 
         Example::
 
-            problem = LinearProblem(a, L, [bc0, bc1], petsc_options={"ksp_type": "preonly",
-                                                                     "pc_type": "lu",
-                                                                     "pc_factor_mat_solver_type":
-                                                                       "mumps"})
+            problem = LinearProblem(a, L, [bc0, bc1], petsc_options={
+                "ksp_type": "preonly",
+                "pc_type": "lu",
+                "pc_factor_mat_solver_type": "mumps"
+            })
         """
         self._a = _create_form(
             a,
@@ -890,8 +894,8 @@ def solver(self) -> PETSc.KSP:
 class NonlinearProblem:
     """Nonlinear problem class for solving the non-linear problems.
 
-    Solves problems of the form :math:`F(u, v) = 0 \\ \\forall v \\in V` using
-    PETSc as the linear algebra backend.
+    Solves problems of the form :math:`F(u, v) = 0 \\ \\forall v \\in V`
+    using PETSc as the linear algebra backend.
     """
 
     def __init__(
@@ -903,7 +907,8 @@ def __init__(
         form_compiler_options: typing.Optional[dict] = None,
         jit_options: typing.Optional[dict] = None,
     ):
-        """Initialize solver for solving a non-linear problem using Newton's method`.
+        """Initialize solver for solving a non-linear problem using
+        Newton's method`.
 
         Args:
             F: The PDE residual F(u, v).
diff --git a/python/dolfinx/geometry.py b/python/dolfinx/geometry.py
index 0960e45f7fa..b43940a677a 100644
--- a/python/dolfinx/geometry.py
+++ b/python/dolfinx/geometry.py
@@ -1,4 +1,5 @@
-# Copyright (C) 2018-2021 Michal Habera, Garth N. Wells and Jørgen S. Dokken
+# Copyright (C) 2018-2021 Michal Habera, Garth N. Wells and
+# Jørgen S. Dokken
 #
 # This file is part of DOLFINx (https://www.fenicsproject.org)
 #
@@ -33,7 +34,8 @@
 
 
 class PointOwnershipData:
-    """Convenience class for storing data related to the ownership of points."""
+    """Convenience class for storing data related to the ownership of
+    points."""
 
     _cpp_object: typing.Union[
         _cpp.geometry.PointOwnershipData_float32, _cpp.geometry.PointOwnershipData_float64
@@ -44,7 +46,8 @@ def __init__(self, ownership_data):
         self._cpp_object = ownership_data
 
     def src_owner(self) -> npt.NDArray[np.int32]:
-        """Ranks owning each point sent into ownership determination for current process."""
+        """Ranks owning each point sent into ownership determination for
+        current process."""
         return self._cpp_object.src_owner
 
     def dest_owner(self) -> npt.NDArray[np.int32]:
@@ -56,7 +59,8 @@ def dest_points(self) -> npt.NDArray[np.floating]:
         return self._cpp_object.dest_points
 
     def dest_cells(self) -> npt.NDArray[np.int32]:
-        """Cell indices (local to process) where each entry of ``dest_points`` is located."""
+        """Cell indices (local to process) where each entry of
+        ``dest_points`` is located."""
         return self._cpp_object.dest_cells
 
 
@@ -193,7 +197,8 @@ def compute_closest_entity(
         midpoint_tree: A bounding box tree with the midpoints of all
             the mesh entities. This is used to accelerate the search.
         mesh: The mesh.
-        points: The points to check for collision, ``shape=(num_points,3)``.
+        points: The points to check for collision,
+            ``shape=(num_points,3)``.
 
     Returns:
         Mesh entity index for each point in ``points``. Returns -1 for a
@@ -206,7 +211,8 @@ def compute_closest_entity(
 
 
 def create_midpoint_tree(mesh: Mesh, dim: int, entities: npt.NDArray[np.int32]) -> BoundingBoxTree:
-    """Create a bounding box tree for the midpoints of a subset of entities.
+    """Create a bounding box tree for the midpoints of a subset of
+    entities.
 
     Args:
         mesh: The mesh.
@@ -229,7 +235,8 @@ def compute_colliding_cells(
         mesh: The mesh.
         candidate_cells: Adjacency list of candidate colliding cells for
             the ith point in ``x``.
-        points: The points to check for collision ``shape=(num_points, 3)``,
+        points: The points to check for collision
+            ``shape=(num_points, 3)``,
 
     Returns:
         Adjacency list where the ith node is the list of entities that
@@ -266,7 +273,8 @@ def squared_distance(
 def compute_distance_gjk(
     p: npt.NDArray[np.floating], q: npt.NDArray[np.floating]
 ) -> npt.NDArray[np.floating]:
-    """Compute the distance between two convex bodies p and q, each defined by a set of points.
+    """Compute the distance between two convex bodies p and q, each defined
+    by a set of points.
 
     Uses the Gilbert-Johnson-Keerthi (GJK) distance algorithm.
 
diff --git a/python/dolfinx/graph.py b/python/dolfinx/graph.py
index 5ad919a1fa2..a4d7d3469b9 100644
--- a/python/dolfinx/graph.py
+++ b/python/dolfinx/graph.py
@@ -41,7 +41,8 @@ def __init__(self, cpp_object: Union[_cpp.la.AdjacencyList_int32, _cpp.la.Adjace
         """Creates a Python wrapper for the exported adjacency list class.
 
         Note:
-            Do not use this constructor directly. Instead use :func:`adjacencylist`.
+            Do not use this constructor directly. Instead use
+            :func:`adjacencylist`.
 
         Args:
             The underlying cpp instance that this object will wrap.
@@ -104,7 +105,8 @@ def adjacencylist(
     Returns:
         An adjacency list.
     """
-    # TODO: Switch to np.isdtype(data.dtype, np.int32) once numpy >= 2.0 is enforced
+    # TODO: Switch to np.isdtype(data.dtype, np.int32) once numpy >= 2.0 is
+    # enforced
     if data.dtype == np.int32:
         cpp_t = _cpp.graph.AdjacencyList_int32
     elif data.dtype == np.int64:
diff --git a/python/dolfinx/io/gmshio.py b/python/dolfinx/io/gmshio.py
index f17e8002102..ef55f078076 100644
--- a/python/dolfinx/io/gmshio.py
+++ b/python/dolfinx/io/gmshio.py
@@ -34,7 +34,8 @@
 
 
 class TopologyDict(typing.TypedDict):
-    """TopologyDict is a TypedDict for storing the topology of the marked cell.
+    """TopologyDict is a TypedDict for storing the topology of the marked
+    cell.
 
     Args:
         topology: 2D array containing the topology of the marked cell.
@@ -95,7 +96,8 @@ class MeshData(typing.NamedTuple):
 
 
 def ufl_mesh(gmsh_cell: int, gdim: int, dtype: npt.DTypeLike) -> ufl.Mesh:
-    """Create a UFL mesh from a Gmsh cell identifier and geometric dimension.
+    """Create a UFL mesh from a Gmsh cell identifier and geometric
+    dimension.
 
     See https://gmsh.info//doc/texinfo/gmsh.html#MSH-file-format.
 
@@ -126,7 +128,8 @@ def ufl_mesh(gmsh_cell: int, gdim: int, dtype: npt.DTypeLike) -> ufl.Mesh:
 
 
 def cell_perm_array(cell_type: CellType, num_nodes: int) -> list[int]:
-    """The permutation array for permuting Gmsh ordering to DOLFINx ordering.
+    """The permutation array for permuting Gmsh ordering to DOLFINx
+    ordering.
 
     Args:
         cell_type: DOLFINx cell type.
@@ -142,7 +145,8 @@ def cell_perm_array(cell_type: CellType, num_nodes: int) -> list[int]:
 def extract_topology_and_markers(
     model, name: typing.Optional[str] = None
 ) -> tuple[dict[int, TopologyDict], dict[str, tuple[int, int]]]:
-    """Extract all entities tagged with a physical marker in the gmsh model.
+    """Extract all entities tagged with a physical marker in the gmsh
+    model.
 
     Returns a nested dictionary where the first key is the gmsh MSH
     element type integer. Each element type present in the model
@@ -155,13 +159,12 @@ def extract_topology_and_markers(
             model will be used.
 
     Returns:
-        A tuple ``(topologies, physical_groups)``, where ``topologies`` is a
-        nested dictionary where each key corresponds to a gmsh cell
-        type. Each cell type found in the mesh has a 2D array containing
-        the topology of the marked cell and a list with the
-        corresponding markers. ``physical_groups`` is a dictionary where the key
-        is the physical name and the value is a tuple with the dimension
-        and tag.
+        A tuple ``(topologies, physical_groups)``, where ``topologies`` is
+        a nested dictionary where each key corresponds to a gmsh cell type.
+        Each cell type found in the mesh has a 2D array containing the
+        topology of the marked cell and a list with the corresponding
+        markers. ``physical_groups`` is a dictionary where the key is the
+        physical name and the value is a tuple with the dimension and tag.
 
     """
     if name is not None:
@@ -311,7 +314,8 @@ def model_to_mesh(
         num_nodes = cell_information[perm_sort[-1]]["num_nodes"]
         cell_id, num_nodes = comm.bcast([cell_id, num_nodes], root=rank)
 
-        # Check for facet, edge and vertex data and broadcast relevant info if True
+        # Check for facet, edge and vertex data and broadcast relevant info
+        # if True
         has_facet_data = (tdim - 1) in cell_dimensions
         has_edge_data = (tdim - 2) in cell_dimensions
         has_vertex_data = (tdim - 3) in cell_dimensions
@@ -455,7 +459,8 @@ def read_from_msh(
         typing.Callable[[_MPI.Comm, int, int, AdjacencyList], AdjacencyList_int32]
     ] = None,
 ) -> MeshData:
-    """Read a Gmsh .msh file and return a :class:`dolfinx.mesh.Mesh` and cell facet markers.
+    """Read a Gmsh .msh file and return a :class:`dolfinx.mesh.Mesh` and
+    cell facet markers.
 
     Note:
         This function requires the Gmsh Python module.
diff --git a/python/dolfinx/io/utils.py b/python/dolfinx/io/utils.py
index b50eb1acf30..5387391e759 100644
--- a/python/dolfinx/io/utils.py
+++ b/python/dolfinx/io/utils.py
@@ -1,5 +1,5 @@
-# Copyright (C) 2017-2022 Chris N. Richardson, Garth N. Wells, Michal Habera
-# and Jørgen S. Dokken
+# Copyright (C) 2017-2022 Chris N. Richardson, Garth N. Wells,
+# Michal Habera and Jørgen S. Dokken
 #
 # This file is part of DOLFINx (https://www.fenicsproject.org)
 #
@@ -135,7 +135,8 @@ def write_mesh(self, mesh: Mesh, t: float = 0.0) -> None:
         self.write(mesh._cpp_object, t)
 
     def write_function(self, u: typing.Union[list[Function], Function], t: float = 0.0) -> None:
-        """Write a single function or a list of functions to file for a given time (default 0.0)"""
+        """Write a single function or a list of functions to file for a
+        given time (default 0.0)"""
         super().write(_extract_cpp_objects(u), t)
 
 
@@ -187,7 +188,8 @@ def read_mesh(
         cells = super().read_topology_data(name, xpath)
         x = super().read_geometry_data(name, xpath)
 
-        # Get coordinate element, special handling for second order serendipity.
+        # Get coordinate element, special handling for second order
+        # serendipity.
         num_nodes_per_cell = cells.shape[1]
         if (cell_shape == CellType.quadrilateral and num_nodes_per_cell == 8) or (
             cell_shape == CellType.hexahedron and num_nodes_per_cell == 20
@@ -224,7 +226,8 @@ def read_meshtags(
         attribute_name: typing.Optional[str] = None,
         xpath: str = "/Xdmf/Domain",
     ) -> MeshTags:
-        """Read MeshTags with a specific name as specified in the XMDF file.
+        """Read MeshTags with a specific name as specified in the XMDF
+        file.
 
         Args:
             mesh: Mesh that the input data is defined on.
@@ -249,7 +252,8 @@ def read_meshtags(
 def distribute_entity_data(
     mesh: Mesh, entity_dim: int, entities: npt.NDArray[np.int64], values: np.ndarray
 ) -> tuple[npt.NDArray[np.int64], np.ndarray]:
-    """Given a set of mesh entities and values, distribute them to the process that owns the entity.
+    """Given a set of mesh entities and values, distribute them to the
+    process that owns the entity.
 
     The entities are described by the global vertex indices of the mesh.
     These entity indices are using the original input ordering.
diff --git a/python/dolfinx/io/vtkhdf.py b/python/dolfinx/io/vtkhdf.py
index 826f2579afe..c36af6cd358 100644
--- a/python/dolfinx/io/vtkhdf.py
+++ b/python/dolfinx/io/vtkhdf.py
@@ -26,10 +26,10 @@ def read_mesh(
 ):
     """Read a mesh from a VTKHDF format file
     Args:
-           comm: An MPI communicator.
-           filename: File to read from.
-           dtype: Scalar type of mesh geometry (need not match dtype in file)
-           gdim: Geometric dimension of the mesh.
+        comm: An MPI communicator.
+        filename: File to read from.
+        dtype: Scalar type of mesh geometry (need not match dtype in file)
+        gdim: Geometric dimension of the mesh.
     """
     if dtype == np.float64:
         mesh_cpp = read_vtkhdf_mesh_float64(comm, filename, gdim)
@@ -58,7 +58,7 @@ def read_mesh(
 def write_mesh(filename: typing.Union[str, Path], mesh: Mesh):
     """Write a mesh to file in VTKHDF format
     Args:
-           filename: File to write to.
-           mesh: Mesh.
+        filename: File to write to.
+        mesh: Mesh.
     """
     write_vtkhdf_mesh(filename, mesh._cpp_object)
diff --git a/python/dolfinx/jit.py b/python/dolfinx/jit.py
index da4574b0358..26c3cd16e3c 100644
--- a/python/dolfinx/jit.py
+++ b/python/dolfinx/jit.py
@@ -193,7 +193,8 @@ def ffcx_jit(
         -  **$XDG_CONFIG_HOME/ffcx/ffcx_options.json** (user options)
         -  **FFCX_DEFAULT_OPTIONS** in `ffcx.options`
 
-        `$XDG_CONFIG_HOME` is `~/.config/` if the environment variable is not set.
+        `$XDG_CONFIG_HOME` is `~/.config/` if the environment variable is
+        not set.
 
         The contents of the `dolfinx_options.json` files are cached
         on the first call. Subsequent calls to this function use this
@@ -201,7 +202,8 @@ def ffcx_jit(
 
         Example `dolfinx_jit_options.json` file:
 
-            **{ "cffi_extra_compile_args": ["-O2", "-march=native" ],  "cffi_verbose": True }**
+            **{ "cffi_extra_compile_args": ["-O2", "-march=native" ],
+            "cffi_verbose": True }**
 
     """
     p_ffcx = ffcx.get_options(form_compiler_options)
diff --git a/python/dolfinx/la/__init__.py b/python/dolfinx/la/__init__.py
index 8fe2e31a553..1a83610c70b 100644
--- a/python/dolfinx/la/__init__.py
+++ b/python/dolfinx/la/__init__.py
@@ -91,8 +91,8 @@ def petsc_vec(self):
         cached and the cached ``Vec`` is returned upon subsequent calls.
 
         Note:
-          When the object is destroyed it will destroy the underlying petsc4py
-          vector automatically.
+          When the object is destroyed it will destroy the underlying
+          petsc4py vector automatically.
         """
         assert dolfinx.has_petsc4py
 
@@ -133,7 +133,8 @@ def __init__(
             _cpp.la.MatrixCSR_complex128,
         ],
     ):
-        """A distributed sparse matrix that uses compressed sparse row storage.
+        """A distributed sparse matrix that uses compressed sparse row
+        storage.
 
         Note:
             Objects of this type should be created using
diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index 46385f0d345..28c5eacdb55 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -1,4 +1,5 @@
-# Copyright (C) 2017-2024 Chris N. Richardson, Garth N. Wells and Jørgen S. Dokken
+# Copyright (C) 2017-2024 Chris N. Richardson, Garth N. Wells and
+# Jørgen S. Dokken
 #
 # This file is part of DOLFINx (https://www.fenicsproject.org)
 #
@@ -81,7 +82,8 @@ def __init__(self, topology: _cpp.mesh.Topology):
 
         Note:
             Topology objects should usually be constructed with the
-            :func:`dolfinx.cpp.mesh.create_topology` and not this class initializer.
+            :func:`dolfinx.cpp.mesh.create_topology` and not this class
+            initializer.
         """
         self._cpp_object = topology
 
@@ -90,7 +92,8 @@ def cell_name(self) -> str:
         return to_string(self._cpp_object.cell_type)
 
     def connectivity(self, d0: int, d1: int) -> _cpp.graph.AdjacencyList_int32:
-        """Return connectivity from entities of dimension ``d0`` to entities of dimension ``d1``.
+        """Return connectivity from entities of dimension ``d0`` to
+        entities of dimension ``d1``.
 
         Args:
             d0: Dimension of entity one is mapping from
@@ -169,7 +172,8 @@ def get_facet_permutations(self) -> npt.NDArray[np.uint8]:
         return self._cpp_object.get_facet_permutations()
 
     def index_map(self, dim: int) -> _cpp.common.IndexMap:
-        """Get the IndexMap that describes the parallel distribution of the mesh entities.
+        """Get the IndexMap that describes the parallel distribution of the
+        mesh entities.
 
         Args:
             dim: Topological dimension.
@@ -186,7 +190,8 @@ def index_map(self, dim: int) -> _cpp.common.IndexMap:
             )
 
     def interprocess_facets(self) -> npt.NDArray[np.int32]:
-        """List of inter-process facets, if facet topology has been computed."""
+        """List of inter-process facets, if facet topology has been
+        computed."""
         return self._cpp_object.interprocess_facets()
 
     @property
@@ -233,11 +238,13 @@ def dim(self):
 
     @property
     def dofmap(self) -> npt.NDArray[np.int32]:
-        """Dofmap for the geometry, shape ``(num_cells, dofs_per_cell)``."""
+        """Dofmap for the geometry, shape
+        ``(num_cells, dofs_per_cell)``."""
         return self._cpp_object.dofmap
 
     def index_map(self) -> _IndexMap:
-        """Index map describing the layout of the geometry points (nodes)."""
+        """Index map describing the layout of the geometry points
+        (nodes)."""
         return self._cpp_object.index_map()
 
     @property
@@ -340,7 +347,8 @@ def geometry(self) -> Geometry:
 
 
 class MeshTags:
-    """Mesh tags associate data (markers) with a subset of mesh entities of a given dimension."""
+    """Mesh tags associate data (markers) with a subset of mesh entities of
+    a given dimension."""
 
     def __init__(self, meshtags):
         """Initialize tags from a C++ MeshTags object.
@@ -408,7 +416,8 @@ def find(self, value) -> npt.NDArray[np.int32]:
 def compute_incident_entities(
     topology: Topology, entities: npt.NDArray[np.int32], d0: int, d1: int
 ) -> npt.NDArray[np.int32]:
-    """Compute all entities of ``d1`` connected to ``entities`` of dimension ``d0``.
+    """Compute all entities of ``d1`` connected to ``entities`` of
+    dimension ``d0``.
 
     Args:
         topology: The topology.
@@ -486,7 +495,8 @@ def transfer_meshtag(
     parent_cell: npt.NDArray[np.int32],
     parent_facet: typing.Optional[npt.NDArray[np.int8]] = None,
 ) -> MeshTags:
-    """Generate cell mesh tags on a refined mesh from the mesh tags on the coarse parent mesh.
+    """Generate cell mesh tags on a refined mesh from the mesh tags on the
+    coarse parent mesh.
 
     Args:
         meshtag: Mesh tags on the coarse, parent mesh.
@@ -576,12 +586,14 @@ def create_mesh(
 
     Args:
         comm: MPI communicator to define the mesh on.
-        cells: Cells of the mesh. ``cells[i]`` are the 'nodes' of cell ``i``.
-        x: Mesh geometry ('node' coordinates), with shape ``(num_nodes, gdim)``.
+            cells: Cells of the mesh. ``cells[i]`` are the 'nodes' of cell
+            ``i``.
+        x: Mesh geometry ('node' coordinates), with shape
+            ``(num_nodes, gdim)``.
         e: UFL mesh. The mesh scalar type is determined by the scalar
             type of ``e``.
-        partitioner: Function that determines the parallel distribution
-            of cells across MPI ranks.
+        partitioner: Function that determines the parallel distribution of
+            cells across MPI ranks.
 
     Note:
         If required, the coordinates ``x`` will be cast to the same
@@ -1024,7 +1036,8 @@ def entities_to_geometry(
 
 
 def exterior_facet_indices(topology: Topology) -> npt.NDArray[np.int32]:
-    """Compute the indices of all exterior facets that are owned by the caller.
+    """Compute the indices of all exterior facets that are owned by the
+    caller.
 
     An exterior facet (co-dimension 1) is one that is connected globally
     to only one cell of co-dimension 0).
@@ -1059,7 +1072,8 @@ def create_geometry(
             row in the point coordinates ``x`` of each local geometry
             node. ``shape=(num_cells, num_dofs_per_cell)``.
         element: Element that describes the cell geometry map.
-        x: The point coordinates. The shape is ``(num_points, geometric_dimension).``
+        x: The point coordinates. The shape is
+            ``(num_points, geometric_dimension).``
         input_global_indices: The 'global' input index of each point,
             commonly from a mesh input file.
     """
diff --git a/python/dolfinx/pyproject.toml b/python/dolfinx/pyproject.toml
new file mode 100644
index 00000000000..945f316fa4f
--- /dev/null
+++ b/python/dolfinx/pyproject.toml
@@ -0,0 +1,5 @@
+[tool.ruff]
+extend = "../pyproject.toml"
+
+[tool.ruff.lint.pycodestyle]
+max-doc-length = 75
diff --git a/python/dolfinx/utils.py b/python/dolfinx/utils.py
index af1fd2eb056..ef1d43b95a0 100644
--- a/python/dolfinx/utils.py
+++ b/python/dolfinx/utils.py
@@ -69,7 +69,8 @@ def set_vals(A: int,
                          m: int, rows: npt.NDArray[PETSc.IntType],
                          n: int, cols: npt.NDArray[PETSc.IntType],
                          data: npt.NDArray[PETSc.ScalarTYpe], mode: int):
-                MatSetValuesLocal(A, m, rows.ctypes, n, cols.ctypes, data.ctypes, mode)
+                MatSetValuesLocal(A, m, rows.ctypes, n, cols.ctypes,
+                                  data.ctypes, mode)
     """
 
     try:
@@ -127,7 +128,8 @@ def set_vals(A: int,
                          m: int, rows: npt.NDArray[PETSc.IntType],
                          n: int, cols: npt.NDArray[PETSc.IntType],
                          data: npt.NDArray[PETSc.ScalarTYpe], mode: int):
-                MatSetValuesLocal(A, m, rows.ctypes, n, cols.ctypes, data.ctypes, mode)
+                MatSetValuesLocal(A, m, rows.ctypes, n, cols.ctypes,
+                                  data.ctypes, mode)
     """
 
     try:
@@ -176,12 +178,13 @@ class cffi_utils:
     Registers Numba's complex types with CFFI.
 
     If PETSc is available, CFFI convenience functions for calling PETSc C
-    functions are also created. These are typically called from within Numba
-    functions.
+    functions are also created. These are typically called from within
+    Numba functions.
 
     Note:
         `CFFI <https://cffi.readthedocs.io/>`_ and  `Numba
-        <https://numba.pydata.org/>`_ must be available to use these utilities.
+        <https://numba.pydata.org/>`_ must be available to use these
+        utilities.
 
     Examples:
         A typical use of these utility functions is::
@@ -192,8 +195,9 @@ def set_vals(A: int,
                          m: int, rows: npt.NDArray[PETSc.IntType],
                          n: int, cols: npt.NDArray[PETSc.IntType],
                          data: npt.NDArray[PETSc.ScalarType], mode: int):
-                MatSetValuesLocal(A, m, ffi.from_buffer(rows), n, ffi.from_buffer(cols),
-                                ffi.from_buffer(rows(data), mode)
+                MatSetValuesLocal(A, m, ffi.from_buffer(rows), n,
+                                  ffi.from_buffer(cols),
+                                  ffi.from_buffer(rows(data), mode)
     """
 
     import cffi as _cffi