diff --git a/python/demo/demo_biharmonic.py b/python/demo/demo_biharmonic.py
index bf8f0d9c9b..4ed8595132 100644
--- a/python/demo/demo_biharmonic.py
+++ b/python/demo/demo_biharmonic.py
@@ -131,7 +131,6 @@
 from dolfinx import fem, io, mesh, plot
 from dolfinx.fem.petsc import LinearProblem
 from dolfinx.mesh import CellType, GhostMode
-from ufl import CellDiameter, FacetNormal, avg, div, dS, dx, grad, inner, jump, pi, sin
 
 # -
 
@@ -192,8 +191,8 @@
 # sides of a facet.
 
 alpha = ScalarType(8.0)
-h = CellDiameter(msh)
-n = FacetNormal(msh)
+h = ufl.CellDiameter(msh)
+n = ufl.FacetNormal(msh)
 h_avg = (h("+") + h("-")) / 2.0
 
 # After that, we can define the variational problem consisting of the bilinear
@@ -210,15 +209,15 @@
 u = ufl.TrialFunction(V)
 v = ufl.TestFunction(V)
 x = ufl.SpatialCoordinate(msh)
-f = 4.0 * pi**4 * sin(pi * x[0]) * sin(pi * x[1])
+f = 4.0 * ufl.pi**4 * ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1])
 
 a = (
-    inner(div(grad(u)), div(grad(v))) * dx
-    - inner(avg(div(grad(u))), jump(grad(v), n)) * dS
-    - inner(jump(grad(u), n), avg(div(grad(v)))) * dS
-    + alpha / h_avg * inner(jump(grad(u), n), jump(grad(v), n)) * dS
+    ufl.inner(ufl.div(ufl.grad(u)), ufl.div(ufl.grad(v))) * ufl.dx
+    - ufl.inner(ufl.avg(ufl.div(ufl.grad(u))), ufl.jump(ufl.grad(v), n)) * ufl.dS
+    - ufl.inner(ufl.jump(ufl.grad(u), n), ufl.avg(ufl.div(ufl.grad(v)))) * ufl.dS
+    + alpha / h_avg * ufl.inner(ufl.jump(ufl.grad(u), n), ufl.jump(ufl.grad(v), n)) * ufl.dS
 )
-L = inner(f, v) * dx
+L = ufl.inner(f, v) * ufl.dx
 # -
 
 # We create a {py:class}`LinearProblem <dolfinx.fem.petsc.LinearProblem>`
diff --git a/python/demo/demo_cahn-hilliard.py b/python/demo/demo_cahn-hilliard.py
index 3cfed7d48e..85e7aa45e9 100644
--- a/python/demo/demo_cahn-hilliard.py
+++ b/python/demo/demo_cahn-hilliard.py
@@ -145,7 +145,6 @@
 from dolfinx.io import XDMFFile
 from dolfinx.mesh import CellType, create_unit_square
 from dolfinx.nls.petsc import NewtonSolver
-from ufl import dx, grad, inner
 
 try:
     import pyvista as pv
@@ -200,7 +199,7 @@
 c0, mu0 = ufl.split(u0)
 # -
 
-# The line `c, mu = ufl.split(u)` permits direct access to the
+# The line `c, mu = split(u)` permits direct access to the
 # components of a mixed function. Note that `c` and `mu` are references
 # for components of `u`, and not copies.
 #
@@ -249,8 +248,16 @@
 # which is then used in the definition of the variational forms:
 
 # Weak statement of the equations
-F0 = inner(c, q) * dx - inner(c0, q) * dx + dt * inner(grad(mu_mid), grad(q)) * dx
-F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
+F0 = (
+    ufl.inner(c, q) * ufl.dx
+    - ufl.inner(c0, q) * ufl.dx
+    + dt * ufl.inner(ufl.grad(mu_mid), ufl.grad(q)) * ufl.dx
+)
+F1 = (
+    ufl.inner(mu, v) * ufl.dx
+    - ufl.inner(dfdc, v) * ufl.dx
+    - lmbda * ufl.inner(ufl.grad(c), ufl.grad(v)) * ufl.dx
+)
 F = F0 + F1
 
 # This is a statement of the time-discrete equations presented as part
diff --git a/python/demo/demo_elasticity.py b/python/demo/demo_elasticity.py
index 8301667b09..fbca3dbdad 100644
--- a/python/demo/demo_elasticity.py
+++ b/python/demo/demo_elasticity.py
@@ -53,7 +53,6 @@
 from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector
 from dolfinx.io import XDMFFile
 from dolfinx.mesh import CellType, GhostMode, create_box, locate_entities_boundary
-from ufl import dx, grad, inner
 
 dtype = PETSc.ScalarType  # type: ignore
 # -
@@ -135,7 +134,7 @@ def build_nullspace(V: FunctionSpace):
 
 def σ(v):
     """Return an expression for the stress σ given a displacement field"""
-    return 2.0 * μ * ufl.sym(grad(v)) + λ * ufl.tr(ufl.sym(grad(v))) * ufl.Identity(len(v))
+    return 2.0 * μ * ufl.sym(ufl.grad(v)) + λ * ufl.tr(ufl.sym(ufl.grad(v))) * ufl.Identity(len(v))
 
 
 # -
@@ -146,8 +145,8 @@ def σ(v):
 
 V = functionspace(msh, ("Lagrange", 1, (msh.geometry.dim,)))
 u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
-a = form(inner(σ(u), grad(v)) * dx)
-L = form(inner(f, v) * dx)
+a = form(ufl.inner(σ(u), ufl.grad(v)) * ufl.dx)
+L = form(ufl.inner(f, v) * ufl.dx)
 
 # A homogeneous (zero) boundary condition is created on $x_0 = 0$ and
 # $x_1 = 1$ by finding all facets on these boundaries, and then creating
@@ -240,7 +239,7 @@ def σ(v):
 
 # +
 sigma_dev = σ(uh) - (1 / 3) * ufl.tr(σ(uh)) * ufl.Identity(len(uh))
-sigma_vm = ufl.sqrt((3 / 2) * inner(sigma_dev, sigma_dev))
+sigma_vm = ufl.sqrt((3 / 2) * ufl.inner(sigma_dev, sigma_dev))
 # -
 
 # Next, the Von Mises stress is interpolated in a piecewise-constant
diff --git a/python/demo/demo_hdg.py b/python/demo/demo_hdg.py
index 8e0e54a686..c40ef304d3 100644
--- a/python/demo/demo_hdg.py
+++ b/python/demo/demo_hdg.py
@@ -41,7 +41,6 @@
 import ufl
 from dolfinx import fem, mesh
 from dolfinx.cpp.mesh import cell_num_entities
-from ufl import div, dot, grad, inner
 
 
 def par_print(comm, string):
@@ -147,17 +146,17 @@ def u_e(x):
 x = ufl.SpatialCoordinate(msh)
 c = 1.0 + 0.1 * ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1])
 a = (
-    inner(c * grad(u), grad(v)) * dx_c
-    - inner(c * (u - ubar), dot(grad(v), n)) * ds_c(cell_boundaries)
-    - inner(dot(grad(u), n), c * (v - vbar)) * ds_c(cell_boundaries)
-    + gamma * inner(c * (u - ubar), v - vbar) * ds_c(cell_boundaries)
+    ufl.inner(c * ufl.grad(u), ufl.grad(v)) * dx_c
+    - ufl.inner(c * (u - ubar), ufl.dot(ufl.grad(v), n)) * ds_c(cell_boundaries)
+    - ufl.inner(ufl.dot(ufl.grad(u), n), c * (v - vbar)) * ds_c(cell_boundaries)
+    + gamma * ufl.inner(c * (u - ubar), v - vbar) * ds_c(cell_boundaries)
 )
 
 # Manufacture a source term
-f = -div(c * grad(u_e(x)))
+f = -ufl.div(c * ufl.grad(u_e(x)))
 
-L = inner(f, v) * dx_c
-L += inner(fem.Constant(facet_mesh, dtype(0.0)), vbar) * dx_f
+L = ufl.inner(f, v) * dx_c
+L += ufl.inner(fem.Constant(facet_mesh, dtype(0.0)), vbar) * dx_f
 
 # Define block structure
 a_blocked = dolfinx.fem.form(ufl.extract_blocks(a), entity_maps=entity_maps)
diff --git a/python/demo/demo_helmholtz.py b/python/demo/demo_helmholtz.py
index 0c9ae774bd..051ba18b2b 100644
--- a/python/demo/demo_helmholtz.py
+++ b/python/demo/demo_helmholtz.py
@@ -51,7 +51,6 @@
 from dolfinx.fem.petsc import LinearProblem
 from dolfinx.io import XDMFFile
 from dolfinx.mesh import create_unit_square
-from ufl import FacetNormal, dot, ds, dx, grad, inner
 
 # Wavenumber
 k0 = 4 * np.pi
@@ -79,14 +78,14 @@
 
 # Define variational problem:
 u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
-a = inner(grad(u), grad(v)) * dx - k0**2 * inner(u, v) * dx
+a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx - k0**2 * ufl.inner(u, v) * ufl.dx
 
 # solve for plane wave with mixed Dirichlet and Neumann BCs
 theta = np.pi / 4
 u_exact.interpolate(lambda x: A * np.exp(1j * k0 * (np.cos(theta) * x[0] + np.sin(theta) * x[1])))
-n = FacetNormal(msh)
-g = -dot(n, grad(u_exact))
-L = -inner(g, v) * ds
+n = ufl.FacetNormal(msh)
+g = -ufl.dot(n, ufl.grad(u_exact))
+L = -ufl.inner(g, v) * ufl.ds
 
 
 dofs_D = locate_dofs_geometrical(
@@ -122,13 +121,17 @@
 
 # H1 errors
 diff = uh - u_exact
-H1_diff = msh.comm.allreduce(assemble_scalar(form(inner(grad(diff), grad(diff)) * dx)), op=MPI.SUM)
+H1_diff = msh.comm.allreduce(
+    assemble_scalar(form(ufl.inner(ufl.grad(diff), ufl.grad(diff)) * ufl.dx)), op=MPI.SUM
+)
 H1_exact = msh.comm.allreduce(
-    assemble_scalar(form(inner(grad(u_exact), grad(u_exact)) * dx)), op=MPI.SUM
+    assemble_scalar(form(ufl.inner(ufl.grad(u_exact), ufl.grad(u_exact)) * ufl.dx)), op=MPI.SUM
 )
 print("Relative H1 error of FEM solution:", abs(np.sqrt(H1_diff) / np.sqrt(H1_exact)))
 
 # L2 errors
-L2_diff = msh.comm.allreduce(assemble_scalar(form(inner(diff, diff) * dx)), op=MPI.SUM)
-L2_exact = msh.comm.allreduce(assemble_scalar(form(inner(u_exact, u_exact) * dx)), op=MPI.SUM)
+L2_diff = msh.comm.allreduce(assemble_scalar(form(ufl.inner(diff, diff) * ufl.dx)), op=MPI.SUM)
+L2_exact = msh.comm.allreduce(
+    assemble_scalar(form(ufl.inner(u_exact, u_exact) * ufl.dx)), op=MPI.SUM
+)
 print("Relative L2 error of FEM solution:", abs(np.sqrt(L2_diff) / np.sqrt(L2_exact)))
diff --git a/python/demo/demo_lagrange_variants.py b/python/demo/demo_lagrange_variants.py
index 9348aa86d2..ef7153be18 100644
--- a/python/demo/demo_lagrange_variants.py
+++ b/python/demo/demo_lagrange_variants.py
@@ -24,9 +24,8 @@
 
 import basix
 import basix.ufl
-import ufl  # type: ignore
+import ufl
 from dolfinx import default_real_type, fem, mesh
-from ufl import dx
 
 # -
 
@@ -175,7 +174,7 @@ def saw_tooth(x):
     V = fem.functionspace(msh, ufl_element)
     uh = fem.Function(V)
     uh.interpolate(lambda x: saw_tooth(x[0]))
-    M = fem.form((u_exact - uh) ** 2 * dx)
+    M = fem.form((u_exact - uh) ** 2 * ufl.dx)
     error = msh.comm.allreduce(fem.assemble_scalar(M), op=MPI.SUM)
     print(f"Computed L2 interpolation error ({variant.name}):", error**0.5)
 # -
diff --git a/python/demo/demo_mixed-poisson.py b/python/demo/demo_mixed-poisson.py
index 51086346b1..44eb2a9103 100644
--- a/python/demo/demo_mixed-poisson.py
+++ b/python/demo/demo_mixed-poisson.py
@@ -110,20 +110,11 @@
 import numpy as np
 
 import dolfinx.fem.petsc
+import ufl
 from basix.ufl import element
 from dolfinx import fem, la, mesh
 from dolfinx.fem.petsc import discrete_gradient, interpolation_matrix
 from dolfinx.mesh import CellType, create_unit_square
-from ufl import (
-    Measure,
-    SpatialCoordinate,
-    TestFunction,
-    TrialFunction,
-    ZeroBaseForm,
-    div,
-    exp,
-    inner,
-)
 
 # Solution scalar (e.g., float32, complex128) and geometry (float32/64)
 # types
@@ -161,16 +152,16 @@
 # $W$, with corresponding test functions $\tau$ and $v$:
 
 # +
-(sigma, u) = TrialFunction(V), TrialFunction(W)
-(tau, v) = TestFunction(V), TestFunction(W)
+(sigma, u) = ufl.TrialFunction(V), ufl.TrialFunction(W)
+(tau, v) = ufl.TestFunction(V), ufl.TestFunction(W)
 # -
 
 # The source function is set to be $f = 10\exp(-((x_{0} - 0.5)^2 +
 # (x_{1} - 0.5)^2) / 0.02)$:
 
 # +
-x = SpatialCoordinate(msh)
-f = 10 * exp(-((x[0] - 0.5) * (x[0] - 0.5) + (x[1] - 0.5) * (x[1] - 0.5)) / 0.02)
+x = ufl.SpatialCoordinate(msh)
+f = 10 * ufl.exp(-((x[0] - 0.5) * (x[0] - 0.5) + (x[1] - 0.5) * (x[1] - 0.5)) / 0.02)
 # -
 
 # We now declare the blocked bilinear and linear forms. The rows of `a`
@@ -184,9 +175,12 @@
 # corresponds to $u_{0} = 0$ on $\Gamma_{D}$.*
 
 # +
-dx = Measure("dx", msh)
-a = [[inner(sigma, tau) * dx, inner(u, div(tau)) * dx], [inner(div(sigma), v) * dx, None]]
-L = [ZeroBaseForm((tau,)), -inner(f, v) * dx]
+dx = ufl.Measure("dx", msh)
+a = [
+    [ufl.inner(sigma, tau) * dx, ufl.inner(u, ufl.div(tau)) * dx],
+    [ufl.inner(ufl.div(sigma), v) * dx, None],
+]
+L = [ufl.ZeroBaseForm((tau,)), -ufl.inner(f, v) * dx]
 # -
 
 # We now compile the abstract/symbolic forms in `a` and `L` into
@@ -263,7 +257,10 @@
 
 # +
 a_p = fem.form(
-    [[inner(sigma, tau) * dx + inner(div(sigma), div(tau)) * dx, None], [None, inner(u, v) * dx]],
+    [
+        [ufl.inner(sigma, tau) * dx + ufl.inner(ufl.div(sigma), ufl.div(tau)) * dx, None],
+        [None, ufl.inner(u, v) * dx],
+    ],
     dtype=dtype,
 )
 P = fem.petsc.assemble_matrix_nest(a_p, bcs=bcs)
diff --git a/python/demo/demo_navier-stokes.py b/python/demo/demo_navier-stokes.py
index a48b9aa3a7..ec5bd5d3ce 100644
--- a/python/demo/demo_navier-stokes.py
+++ b/python/demo/demo_navier-stokes.py
@@ -179,27 +179,9 @@
 # +
 import numpy as np
 
+import ufl
 from dolfinx import default_real_type, fem, io, mesh
 from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
-from ufl import (
-    CellDiameter,
-    FacetNormal,
-    MixedFunctionSpace,
-    TestFunctions,
-    TrialFunctions,
-    avg,
-    conditional,
-    div,
-    dot,
-    dS,
-    ds,
-    dx,
-    extract_blocks,
-    grad,
-    gt,
-    inner,
-    outer,
-)
 
 try:
     from petsc4py import PETSc
@@ -225,15 +207,18 @@
 # +
 def norm_L2(comm, v):
     """Compute the L2(Ω)-norm of v"""
-    return np.sqrt(comm.allreduce(fem.assemble_scalar(fem.form(inner(v, v) * dx)), op=MPI.SUM))
+    return np.sqrt(
+        comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(v, v) * ufl.dx)), op=MPI.SUM)
+    )
 
 
 def domain_average(msh, v):
     """Compute the average of a function over the domain"""
     vol = msh.comm.allreduce(
-        fem.assemble_scalar(fem.form(fem.Constant(msh, default_real_type(1.0)) * dx)), op=MPI.SUM
+        fem.assemble_scalar(fem.form(fem.Constant(msh, default_real_type(1.0)) * ufl.dx)),
+        op=MPI.SUM,
     )
-    return (1 / vol) * msh.comm.allreduce(fem.assemble_scalar(fem.form(v * dx)), op=MPI.SUM)
+    return (1 / vol) * msh.comm.allreduce(fem.assemble_scalar(fem.form(v * ufl.dx)), op=MPI.SUM)
 
 
 def u_e_expr(x):
@@ -282,7 +267,7 @@ def f_expr(x):
 # Function spaces for the velocity and for the pressure
 V = fem.functionspace(msh, ("Raviart-Thomas", k + 1))
 Q = fem.functionspace(msh, ("Discontinuous Lagrange", k))
-VQ = MixedFunctionSpace(V, Q)
+VQ = ufl.MixedFunctionSpace(V, Q)
 
 # Function space for visualising the velocity field
 gdim = msh.geometry.dim
@@ -290,18 +275,18 @@ def f_expr(x):
 
 # Define trial and test functions
 
-u, p = TrialFunctions(VQ)
-v, q = TestFunctions(VQ)
+u, p = ufl.TrialFunctions(VQ)
+v, q = ufl.TestFunctions(VQ)
 
 delta_t = fem.Constant(msh, default_real_type(t_end / num_time_steps))
 alpha = fem.Constant(msh, default_real_type(6.0 * k**2))
 
-h = CellDiameter(msh)
-n = FacetNormal(msh)
+h = ufl.CellDiameter(msh)
+n = ufl.FacetNormal(msh)
 
 
 def jump(phi, n):
-    return outer(phi("+"), n("+")) + outer(phi("-"), n("-"))
+    return ufl.outer(phi("+"), n("+")) + ufl.outer(phi("-"), n("-"))
 
 
 # -
@@ -311,27 +296,28 @@ def jump(phi, n):
 
 # +
 a = (1.0 / Re) * (
-    inner(grad(u), grad(v)) * dx
-    - inner(avg(grad(u)), jump(v, n)) * dS
-    - inner(jump(u, n), avg(grad(v))) * dS
-    + (alpha / avg(h)) * inner(jump(u, n), jump(v, n)) * dS
-    - inner(grad(u), outer(v, n)) * ds
-    - inner(outer(u, n), grad(v)) * ds
-    + (alpha / h) * inner(outer(u, n), outer(v, n)) * ds
+    ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
+    - ufl.inner(ufl.avg(ufl.grad(u)), jump(v, n)) * ufl.dS
+    - ufl.inner(jump(u, n), ufl.avg(ufl.grad(v))) * ufl.dS
+    + (alpha / ufl.avg(h)) * ufl.inner(jump(u, n), jump(v, n)) * ufl.dS
+    - ufl.inner(ufl.grad(u), ufl.outer(v, n)) * ufl.ds
+    - ufl.inner(ufl.outer(u, n), ufl.grad(v)) * ufl.ds
+    + (alpha / h) * ufl.inner(ufl.outer(u, n), ufl.outer(v, n)) * ufl.ds
 )
-a -= inner(p, div(v)) * dx
-a -= inner(div(u), q) * dx
+a -= ufl.inner(p, ufl.div(v)) * ufl.dx
+a -= ufl.inner(ufl.div(u), q) * ufl.dx
 
-a_blocked = fem.form(extract_blocks(a))
+a_blocked = fem.form(ufl.extract_blocks(a))
 
 f = fem.Function(W)
 u_D = fem.Function(V)
 u_D.interpolate(u_e_expr)
-L = inner(f, v) * dx + (1 / Re) * (
-    -inner(outer(u_D, n), grad(v)) * ds + (alpha / h) * inner(outer(u_D, n), outer(v, n)) * ds
+L = ufl.inner(f, v) * ufl.dx + (1 / Re) * (
+    -ufl.inner(ufl.outer(u_D, n), ufl.grad(v)) * ufl.ds
+    + (alpha / h) * ufl.inner(ufl.outer(u_D, n), ufl.outer(v, n)) * ufl.ds
 )
-L += inner(fem.Constant(msh, default_real_type(0.0)), q) * dx
-L_blocked = fem.form(extract_blocks(L))
+L += ufl.inner(fem.Constant(msh, default_real_type(0.0)), q) * ufl.dx
+L_blocked = fem.form(ufl.extract_blocks(L))
 
 # Boundary conditions
 boundary_facets = mesh.exterior_facet_indices(msh.topology)
@@ -404,19 +390,22 @@ def jump(phi, n):
 # Now we add the time stepping and convective terms
 
 # +
-lmbda = conditional(gt(dot(u_n, n), 0), 1, 0)
+lmbda = ufl.conditional(ufl.gt(ufl.dot(u_n, n), 0), 1, 0)
 u_uw = lmbda("+") * u("+") + lmbda("-") * u("-")
 a += (
-    inner(u / delta_t, v) * dx
-    - inner(u, div(outer(v, u_n))) * dx
-    + inner((dot(u_n, n))("+") * u_uw, v("+")) * dS
-    + inner((dot(u_n, n))("-") * u_uw, v("-")) * dS
-    + inner(dot(u_n, n) * lmbda * u, v) * ds
+    ufl.inner(u / delta_t, v) * ufl.dx
+    - ufl.inner(u, ufl.div(ufl.outer(v, u_n))) * ufl.dx
+    + ufl.inner((ufl.dot(u_n, n))("+") * u_uw, v("+")) * ufl.dS
+    + ufl.inner((ufl.dot(u_n, n))("-") * u_uw, v("-")) * ufl.dS
+    + ufl.inner(ufl.dot(u_n, n) * lmbda * u, v) * ufl.ds
 )
-a_blocked = fem.form(extract_blocks(a))
+a_blocked = fem.form(ufl.extract_blocks(a))
 
-L += inner(u_n / delta_t, v) * dx - inner(dot(u_n, n) * (1 - lmbda) * u_D, v) * ds
-L_blocked = fem.form(extract_blocks(L))
+L += (
+    ufl.inner(u_n / delta_t, v) * ufl.dx
+    - ufl.inner(ufl.dot(u_n, n) * (1 - lmbda) * u_D, v) * ufl.ds
+)
+L_blocked = fem.form(ufl.extract_blocks(L))
 
 # Time stepping loop
 for n in range(num_time_steps):
@@ -473,7 +462,7 @@ def jump(phi, n):
 
 # Compute errors
 e_u = norm_L2(msh.comm, u_h - u_e)
-e_div_u = norm_L2(msh.comm, div(u_h))
+e_div_u = norm_L2(msh.comm, ufl.div(u_h))
 
 # This scheme conserves mass exactly, so check this
 assert np.isclose(e_div_u, 0.0, atol=float(1.0e5 * np.finfo(default_real_type).eps))
diff --git a/python/demo/demo_poisson.py b/python/demo/demo_poisson.py
index 38c072211e..0ee822c719 100644
--- a/python/demo/demo_poisson.py
+++ b/python/demo/demo_poisson.py
@@ -85,7 +85,6 @@
 import ufl
 from dolfinx import fem, io, mesh, plot
 from dolfinx.fem.petsc import LinearProblem
-from ufl import ds, dx, grad, inner
 
 # -
 
@@ -146,8 +145,8 @@
 x = ufl.SpatialCoordinate(msh)
 f = 10 * ufl.exp(-((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) / 0.02)
 g = ufl.sin(5 * x[0])
-a = inner(grad(u), grad(v)) * dx
-L = inner(f, v) * dx + inner(g, v) * ds
+a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
+L = ufl.inner(f, v) * ufl.dx + ufl.inner(g, v) * ufl.ds
 # -
 
 # A {py:class}`LinearProblem <dolfinx.fem.petsc.LinearProblem>` object is
diff --git a/python/demo/demo_poisson_matrix_free.py b/python/demo/demo_poisson_matrix_free.py
index b29d998040..fa5fe103bd 100644
--- a/python/demo/demo_poisson_matrix_free.py
+++ b/python/demo/demo_poisson_matrix_free.py
@@ -83,7 +83,6 @@
 import dolfinx
 import ufl
 from dolfinx import fem, la
-from ufl import action, dx, grad, inner
 
 # We begin by using {py:func}`create_rectangle
 # <dolfinx.mesh.create_rectangle>` to create a rectangular
@@ -138,8 +137,8 @@
 u = ufl.TrialFunction(V)
 v = ufl.TestFunction(V)
 f = fem.Constant(mesh, dtype(-6.0))
-a = inner(grad(u), grad(v)) * dx
-L = inner(f, v) * dx
+a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
+L = ufl.inner(f, v) * ufl.dx
 L_fem = fem.form(L, dtype=dtype)
 
 # For the matrix-free solvers we also define a second linear form `M` as
@@ -152,7 +151,7 @@
 # $$
 
 ui = fem.Function(V, dtype=dtype)
-M = action(a, ui)
+M = ufl.action(a, ui)
 M_fem = fem.form(M, dtype=dtype)
 
 # ### Matrix-free conjugate gradient solver
@@ -254,7 +253,7 @@ def _global_dot(comm, v0, v1):
 
 
 def L2Norm(u):
-    val = fem.assemble_scalar(fem.form(inner(u, u) * dx, dtype=dtype))
+    val = fem.assemble_scalar(fem.form(ufl.inner(u, u) * ufl.dx, dtype=dtype))
     return np.sqrt(comm.allreduce(val, op=MPI.SUM))
 
 
diff --git a/python/demo/demo_pyamg.py b/python/demo/demo_pyamg.py
index 69de02486f..dab612f9a1 100644
--- a/python/demo/demo_pyamg.py
+++ b/python/demo/demo_pyamg.py
@@ -48,7 +48,6 @@
     locate_dofs_topological,
 )
 from dolfinx.mesh import CellType, create_box, locate_entities_boundary
-from ufl import ds, dx, grad, inner
 
 # -
 
@@ -88,8 +87,8 @@ def poisson_problem(dtype: npt.DTypeLike, solver_type: str):
     x = ufl.SpatialCoordinate(mesh)
     f = 10 * ufl.exp(-((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) / 0.02)
     g = ufl.sin(5 * x[0])
-    a = form(inner(grad(u), grad(v)) * dx, dtype=dtype)
-    L = form(inner(f, v) * dx + inner(g, v) * ds, dtype=dtype)
+    a = form(ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx, dtype=dtype)
+    L = form(ufl.inner(f, v) * ufl.dx + ufl.inner(g, v) * ufl.ds, dtype=dtype)
 
     A = assemble_matrix(a, [bc]).to_scipy()
     b = assemble_vector(L)
@@ -193,12 +192,14 @@ def elasticity_problem(dtype):
 
     def σ(v):
         """Return an expression for the stress σ given a displacement field"""
-        return 2.0 * μ * ufl.sym(grad(v)) + λ * ufl.tr(ufl.sym(grad(v))) * ufl.Identity(len(v))
+        return 2.0 * μ * ufl.sym(ufl.grad(v)) + λ * ufl.tr(ufl.sym(ufl.grad(v))) * ufl.Identity(
+            len(v)
+        )
 
     V = functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim,)))
     u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
-    a = form(inner(σ(u), grad(v)) * dx, dtype=dtype)
-    L = form(inner(f, v) * dx, dtype=dtype)
+    a = form(ufl.inner(σ(u), ufl.grad(v)) * ufl.dx, dtype=dtype)
+    L = form(ufl.inner(f, v) * ufl.dx, dtype=dtype)
 
     tdim = mesh.topology.dim
     dofs = locate_dofs_topological(V=V, entity_dim=tdim - 1, entities=facets)
diff --git a/python/demo/demo_stokes.py b/python/demo/demo_stokes.py
index 599063f4fc..5b5df67b4c 100644
--- a/python/demo/demo_stokes.py
+++ b/python/demo/demo_stokes.py
@@ -116,7 +116,6 @@
 from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
 from dolfinx.io import XDMFFile
 from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
-from ufl import ZeroBaseForm, div, dx, grad, inner
 
 # We create a {py:class}`Mesh <dolfinx.mesh.Mesh>`, define functions for
 # locating geometrically subsets of the boundary, and define a function
@@ -185,14 +184,19 @@ def lid_velocity_expression(x):
 (v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
 f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))  # type: ignore
 
-a = form([[inner(grad(u), grad(v)) * dx, inner(p, div(v)) * dx], [inner(div(u), q) * dx, None]])
-L = form([inner(f, v) * dx, ZeroBaseForm((q,))])
+a = form(
+    [
+        [ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx, ufl.inner(p, ufl.div(v)) * ufl.dx],
+        [ufl.inner(ufl.div(u), q) * ufl.dx, None],
+    ]
+)
+L = form([ufl.inner(f, v) * ufl.dx, ufl.ZeroBaseForm((q,))])
 # -
 
 # A block-diagonal preconditioner will be used with the iterative
 # solvers for this problem:
 
-a_p11 = form(inner(p, q) * dx)
+a_p11 = form(ufl.inner(p, q) * ufl.dx)
 a_p = [[a[0][0], None], [None, a_p11]]
 
 # ### Nested matrix solver
@@ -503,8 +507,11 @@ def mixed_direct():
     (u, p) = ufl.TrialFunctions(W)
     (v, q) = ufl.TestFunctions(W)
     f = Function(Q)
-    a = form((inner(grad(u), grad(v)) + inner(p, div(v)) + inner(div(u), q)) * dx)
-    L = form(inner(f, v) * dx)
+    a = form(
+        (ufl.inner(ufl.grad(u), ufl.grad(v)) + ufl.inner(p, ufl.div(v)) + ufl.inner(ufl.div(u), q))
+        * ufl.dx
+    )
+    L = form(ufl.inner(f, v) * ufl.dx)
 
     # Assemble LHS matrix and RHS vector
     A = fem.petsc.assemble_matrix(a, bcs=bcs)
diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index 06eadf560b..810d18533a 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -41,9 +41,9 @@
 
 import basix
 import basix.ufl
+import ufl
 from dolfinx import default_real_type, fem, mesh
 from dolfinx.fem.petsc import LinearProblem
-from ufl import SpatialCoordinate, TestFunction, TrialFunction, cos, div, dx, grad, inner, sin
 
 mpl.use("agg")
 # -
@@ -256,14 +256,14 @@ def create_tnt_quad(degree):
 
 def poisson_error(V: fem.FunctionSpace):
     msh = V.mesh
-    u, v = TrialFunction(V), TestFunction(V)
+    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
 
-    x = SpatialCoordinate(msh)
-    u_exact = sin(10 * x[1]) * cos(15 * x[0])
-    f = -div(grad(u_exact))
+    x = ufl.SpatialCoordinate(msh)
+    u_exact = ufl.sin(10 * x[1]) * ufl.cos(15 * x[0])
+    f = -ufl.div(ufl.grad(u_exact))
 
-    a = inner(grad(u), grad(v)) * dx
-    L = inner(f, v) * dx
+    a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
+    L = ufl.inner(f, v) * ufl.dx
 
     # Create Dirichlet boundary condition
     u_bc = fem.Function(V)
@@ -278,7 +278,7 @@ def poisson_error(V: fem.FunctionSpace):
     problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_rtol": 1e-12})
     uh = problem.solve()
 
-    M = (u_exact - uh) ** 2 * dx
+    M = (u_exact - uh) ** 2 * ufl.dx
     M = fem.form(M)
     error = msh.comm.allreduce(fem.assemble_scalar(M), op=MPI.SUM)
     return error**0.5
diff --git a/python/demo/pyproject.toml b/python/demo/pyproject.toml
new file mode 100644
index 0000000000..28e7749be2
--- /dev/null
+++ b/python/demo/pyproject.toml
@@ -0,0 +1,5 @@
+[tool.ruff]
+extend = "../pyproject.toml"
+
+[tool.ruff.lint.flake8-import-conventions]
+banned-from = ["ufl"]