diff --git a/docs/source/expression_tree/unary_operator.rst b/docs/source/expression_tree/unary_operator.rst
index 4848ed6f5..c0f9e04a3 100644
--- a/docs/source/expression_tree/unary_operator.rst
+++ b/docs/source/expression_tree/unary_operator.rst
@@ -33,3 +33,7 @@ Unary Operators
 .. autofunction:: pybamm.div
 
 .. autofunction:: pybamm.surf
+
+.. autofunction:: pybamm.boundary_value
+
+.. autofunction:: pybamm.average
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 4ee569b8e..ec9eabadb 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -80,11 +80,12 @@ from .expression_tree.unary_operators import (
     BoundaryValue,
     Integral,
     IndefiniteIntegral,
+    Diagonal,
     grad,
     div,
     surf,
-    Diagonal,
-    integrate,
+    average,
+    boundary_value,
 )
 from .expression_tree.parameter import Parameter, FunctionParameter
 from .expression_tree.broadcasts import Broadcast, NumpyBroadcast
diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index 42368a58c..b9f4a5cc2 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -399,8 +399,9 @@ def simplify_multiplication_division(myclass, left, right):
             if nonconst_denominator_expr is None:
                 result = constant_numerator_expr * new_numerator
             else:
-                result = constant_numerator_expr * new_numerator \
-                    / nonconst_denominator_expr
+                result = (
+                    constant_numerator_expr * new_numerator / nonconst_denominator_expr
+                )
 
     elif not numerator_has_mat_mul and denominator_has_mat_mul:
         new_denominator = simplify_with_mat_mul(denominator, denominator_types)
@@ -422,8 +423,9 @@ def simplify_multiplication_division(myclass, left, right):
             if nonconst_numerator_expr is None:
                 result = constant_numerator_expr / new_denominator
             else:
-                result = constant_numerator_expr * nonconst_numerator_expr \
-                    / new_denominator
+                result = (
+                    constant_numerator_expr * nonconst_numerator_expr / new_denominator
+                )
 
     else:
         # can reorder the numerator since no matrix multiplies
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 24fac2f89..3d99accbc 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -6,6 +6,7 @@ from __future__ import print_function, unicode_literals
 import pybamm
 
 import autograd
+import copy
 import numpy as np
 import numbers
 from scipy.sparse import csr_matrix, diags
@@ -488,7 +489,8 @@ def div(expression):
 
 
 def surf(variable):
-    """convenience function for creating a :class:`SurfaceValue`
+    """convenience function for creating a right :class:`BoundaryValue`, usually in the
+    spherical geometry
 
     Parameters
     ----------
@@ -503,25 +505,66 @@ def surf(variable):
         the surface value of ``variable``
     """
 
-    return BoundaryValue(variable, "right")
+    return boundary_value(variable, "right")
 
 
-def integrate(expression, variable):
-    """convenience function for creating a :class:`Integral`
+def average(symbol):
+    """convenience function for creating an average
 
     Parameters
     ----------
-
-    expression : :class:`pybamm.Symbol`
-        The function to be integrated
-    integration_variable : :class:`pybamm.IndependentVariable`
-        The variable over which to integrate
+    symbol : :class:`pybamm.Symbol`
+        The function to be averaged
 
     Returns
     -------
+    :class:`Symbol`
+        the new averaged symbol
+    """
+    # If symbol doesn't have a domain, its average value is itself
+    if symbol.domain == []:
+        new_symbol = copy.deepcopy(symbol)
+        new_symbol.parent = None
+        return new_symbol
+    # If symbol is a Broadcast, its average value is its child
+    elif isinstance(symbol, pybamm.Broadcast):
+        return symbol.orphans[0]
+    # Otherwise, use Integral to calculate average value
+    else:
+        if symbol.domain == ["negative electrode"]:
+            x = pybamm.standard_spatial_vars.x_n
+            l = pybamm.geometric_parameters.l_n
+        elif symbol.domain == ["separator"]:
+            x = pybamm.standard_spatial_vars.x_s
+            l = pybamm.geometric_parameters.l_s
+        elif symbol.domain == ["positive electrode"]:
+            x = pybamm.standard_spatial_vars.x_p
+            l = pybamm.geometric_parameters.l_p
+
+        return Integral(symbol, x) / l
+
+
+def boundary_value(symbol, side):
+    """convenience function for creating a :class:`Integral`
+
+    Parameters
+    ----------
+    symbol : `pybamm.Symbol`
+        The symbol whose boundary value to take
+    side : str
+        Which side to take the boundary value on ("left" or "right")
 
-    :class:`Integral`
+    Returns
+    -------
+    :class:`BoundaryValue`
         the new integrated expression tree
     """
-
-    return Integral(expression, variable)
+    # If symbol doesn't have a domain, its boundary value is itself
+    if symbol.domain == []:
+        return symbol
+    # If symbol is a Broadcast, its boundary value is its child
+    elif isinstance(symbol, pybamm.Broadcast):
+        return symbol.children[0]
+    # Otherwise, calculate boundary value
+    else:
+        return BoundaryValue(symbol, side)

Boundary Value and Average Value of a constant / broadcasted variable
Each of these should just return the variable itself, since it is constant. Possibly should be left to after simplifying, but I think it is a trivial enough simplification to do straight away via helper functions. e.g.
```python3
def boundary(symbol, side):
   if symbol.domain == []:
       return symbol 
   elif isinstance(symbol, pybamm.Broadcast):
       return symbol.orphans[0] 
   else: 
       return BoundaryValue(symbol, side)
```
diff --git a/tests/unit/test_expression_tree/test_matrix.py b/tests/unit/test_expression_tree/test_matrix.py
index 5d8128c3d..dbaf50f8e 100644
--- a/tests/unit/test_expression_tree/test_matrix.py
+++ b/tests/unit/test_expression_tree/test_matrix.py
@@ -157,11 +157,11 @@ class TestMatrix(unittest.TestCase):
         # we expect simplified solution to be much faster
         timer = pybamm.Timer()
         start = timer.time()
-        for _ in range(20):
+        for _ in range(100):
             expr3.evaluate(y=np.ones(300))
         end = timer.time()
         start_simp = timer.time()
-        for _ in range(20):
+        for _ in range(100):
             expr3simp.evaluate(y=np.ones(300))
         end_simp = timer.time()
         self.assertLess(end_simp - start_simp, 1.5 * (end - start))
diff --git a/tests/unit/test_expression_tree/test_symbol.py b/tests/unit/test_expression_tree/test_symbol.py
index 8a843415e..3b2845722 100644
--- a/tests/unit/test_expression_tree/test_symbol.py
+++ b/tests/unit/test_expression_tree/test_symbol.py
@@ -56,11 +56,12 @@ class TestSymbol(unittest.TestCase):
 
         # Integral
         self.assertIsInstance(
-            (pybamm.integrate(a, pybamm.t)).simplify(), pybamm.Integral
+            (pybamm.Integral(a, pybamm.t)).simplify(), pybamm.Integral
         )
 
         # BoundaryValue
-        self.assertIsInstance((pybamm.surf(v)).simplify(), pybamm.BoundaryValue)
+        v_neg = pybamm.Variable("v", domain=["negative electrode"])
+        self.assertIsInstance((pybamm.surf(v_neg)).simplify(), pybamm.BoundaryValue)
 
         # addition
         self.assertIsInstance((a + b).simplify(), pybamm.Scalar)
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index aa2f2a78c..8efb80df7 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -100,18 +100,12 @@ class TestUnaryOperators(unittest.TestCase):
         self.assertEqual(inta.integration_variable, x)
         self.assertEqual(inta.domain, [])
 
-        # # Indefinite
-        # for inta in [
-        #     pybamm.Integral(a, x, definite=False),
-        #     pybamm.IndefiniteIntegral(a, x),
-        # ]:
-        #     self.assertEqual(
-        #         inta.name, "indefinite integral dspace (['negative electrode'])"
-        #     )
-        #     self.assertFalse(inta.definite)
-        #     self.assertEqual(inta.children[0].name, a.name)
-        #     self.assertEqual(inta.integration_variable, x)
-        #     self.assertEqual(inta.domain, [])
+        # Indefinite
+        inta = pybamm.IndefiniteIntegral(a, x)
+        self.assertEqual(inta.name, "a integrated w.r.t x on ['negative electrode']")
+        self.assertEqual(inta.children[0].name, a.name)
+        self.assertEqual(inta.integration_variable, x)
+        self.assertEqual(inta.domain, ["negative electrode"])
 
         # expected errors
         a = pybamm.Symbol("a", domain=["negative electrode"])
@@ -172,6 +166,39 @@ class TestUnaryOperators(unittest.TestCase):
         un5 = pybamm.UnaryOperator("test", d)
         self.assertNotEqual(un1.id, un5.id)
 
+    def test_boundary_value(self):
+        a = pybamm.Symbol("a")
+        boundary_a = pybamm.boundary_value(a, "right")
+        self.assertEqual(boundary_a.id, a.id)
+
+        boundary_broad_a = pybamm.boundary_value(
+            pybamm.Broadcast(a, ["negative electrode"]), "left"
+        )
+        self.assertEqual(boundary_broad_a.id, a.id)
+
+        a = pybamm.Symbol("a", domain=["separator"])
+        boundary_a = pybamm.boundary_value(a, "right")
+        self.assertIsInstance(boundary_a, pybamm.BoundaryValue)
+        self.assertEqual(boundary_a.side, "right")
+        self.assertEqual(boundary_a.domain, [])
+
+    def test_average(self):
+        a = pybamm.Symbol("a")
+        average_a = pybamm.average(a)
+        self.assertEqual(average_a.id, a.id)
+
+        average_broad_a = pybamm.average(pybamm.Broadcast(a, ["negative electrode"]))
+        self.assertEqual(average_broad_a.id, a.id)
+
+        for domain in [["negative electrode"], ["separator"], ["positive electrode"]]:
+            a = pybamm.Symbol("a", domain=domain)
+            x = pybamm.SpatialVariable("x", domain)
+            av_a = pybamm.average(a)
+            self.assertIsInstance(av_a, pybamm.Division)
+            self.assertIsInstance(av_a.children[0], pybamm.Integral)
+            self.assertEqual(av_a.children[0].integration_variable.domain, x.domain)
+            self.assertEqual(av_a.domain, [])
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

pip install -e .[dev,docs]
parse_log_pytest
requirements.txt
pytest
apt-get update apt-get install -y gfortran gcc libopenblas-dev liblapack-dev
3.6
requirements-docs.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
alabaster==0.7.13
anyio==3.6.2
anytree==2.8.0
argon2-cffi==21.3.0
argon2-cffi-bindings==21.2.0
async-generator==1.10
attrs==22.2.0
autograd==1.6.2
Babel==2.11.0
backcall==0.2.0
bleach==4.1.0
certifi==2021.5.30
cffi==1.15.1
charset-normalizer==2.0.12
comm==0.1.4
contextvars==2.4
cycler==0.11.0
dataclasses==0.8
decorator==5.1.1
defusedxml==0.7.1
docutils==0.17.1
entrypoints==0.4
flake8==5.0.4
future==1.0.0
guzzle-sphinx-theme==0.7.11
idna==3.10
imagesize==1.4.1
immutables==0.19
importlib-metadata==4.2.0
iniconfig==1.1.1
ipykernel==5.5.6
ipython==7.16.3
ipython-genutils==0.2.0
ipywidgets==7.8.5
jedi==0.17.2
Jinja2==3.0.3
json5==0.9.16
jsonschema==3.2.0
jupyter==1.1.1
jupyter-client==7.1.2
jupyter-console==6.4.3
jupyter-core==4.9.2
jupyter-server==1.13.1
jupyterlab==3.2.9
jupyterlab-pygments==0.1.2
jupyterlab-server==2.10.3
jupyterlab_widgets==1.1.11
kiwisolver==1.3.1
MarkupSafe==2.0.1
matplotlib==3.3.4
mccabe==0.7.0
mistune==0.8.4
nbclassic==0.3.5
nbclient==0.5.9
nbconvert==6.0.7
nbformat==5.1.3
nest-asyncio==1.6.0
notebook==6.4.10
numpy==1.19.5
packaging==21.3
pandas==1.1.5
pandocfilters==1.5.1
parso==0.7.1
pexpect==4.9.0
pickleshare==0.7.5
Pillow==8.4.0
pluggy==1.0.0
prometheus-client==0.17.1
prompt-toolkit==3.0.36
ptyprocess==0.7.0
py==1.11.0
-e git+https://github.com/pybamm-team/PyBaMM.git@afebdce35b6c8ef9c0f21e7009e36cae18aa85cb#egg=pybamm
pycodestyle==2.9.1
pycparser==2.21
pyflakes==2.5.0
Pygments==2.14.0
pyparsing==3.1.4
pyrsistent==0.18.0
pytest==7.0.1
python-dateutil==2.9.0.post0
pytz==2025.2
pyzmq==25.1.2
requests==2.27.1
scipy==1.5.4
Send2Trash==1.8.3
six==1.17.0
sniffio==1.2.0
snowballstemmer==2.2.0
Sphinx==4.3.2
sphinxcontrib-applehelp==1.0.2
sphinxcontrib-devhelp==1.0.2
sphinxcontrib-htmlhelp==2.0.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==1.0.3
sphinxcontrib-serializinghtml==1.1.5
terminado==0.12.1
testpath==0.6.0
tomli==1.2.3
tornado==6.1
traitlets==4.3.3
typing_extensions==4.1.1
urllib3==1.26.20
wcwidth==0.2.13
webencodings==0.5.1
websocket-client==1.3.1
widgetsnbextension==3.6.10
zipp==3.6.0

name: PyBaMM
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2021.5.30=py36h06a4308_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.3=he6710b0_2
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - pip=21.2.2=py36h06a4308_0
  - python=3.6.13=h12debd9_1
  - readline=8.2=h5eee18b_0
  - setuptools=58.0.4=py36h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - wheel=0.37.1=pyhd3eb1b0_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - alabaster==0.7.13
      - anyio==3.6.2
      - anytree==2.8.0
      - argon2-cffi==21.3.0
      - argon2-cffi-bindings==21.2.0
      - async-generator==1.10
      - attrs==22.2.0
      - autograd==1.6.2
      - babel==2.11.0
      - backcall==0.2.0
      - bleach==4.1.0
      - cffi==1.15.1
      - charset-normalizer==2.0.12
      - comm==0.1.4
      - contextvars==2.4
      - cycler==0.11.0
      - dataclasses==0.8
      - decorator==5.1.1
      - defusedxml==0.7.1
      - docutils==0.17.1
      - entrypoints==0.4
      - flake8==5.0.4
      - future==1.0.0
      - guzzle-sphinx-theme==0.7.11
      - idna==3.10
      - imagesize==1.4.1
      - immutables==0.19
      - importlib-metadata==4.2.0
      - iniconfig==1.1.1
      - ipykernel==5.5.6
      - ipython==7.16.3
      - ipython-genutils==0.2.0
      - ipywidgets==7.8.5
      - jedi==0.17.2
      - jinja2==3.0.3
      - json5==0.9.16
      - jsonschema==3.2.0
      - jupyter==1.1.1
      - jupyter-client==7.1.2
      - jupyter-console==6.4.3
      - jupyter-core==4.9.2
      - jupyter-server==1.13.1
      - jupyterlab==3.2.9
      - jupyterlab-pygments==0.1.2
      - jupyterlab-server==2.10.3
      - jupyterlab-widgets==1.1.11
      - kiwisolver==1.3.1
      - markupsafe==2.0.1
      - matplotlib==3.3.4
      - mccabe==0.7.0
      - mistune==0.8.4
      - nbclassic==0.3.5
      - nbclient==0.5.9
      - nbconvert==6.0.7
      - nbformat==5.1.3
      - nest-asyncio==1.6.0
      - notebook==6.4.10
      - numpy==1.19.5
      - packaging==21.3
      - pandas==1.1.5
      - pandocfilters==1.5.1
      - parso==0.7.1
      - pexpect==4.9.0
      - pickleshare==0.7.5
      - pillow==8.4.0
      - pluggy==1.0.0
      - prometheus-client==0.17.1
      - prompt-toolkit==3.0.36
      - ptyprocess==0.7.0
      - py==1.11.0
      - pycodestyle==2.9.1
      - pycparser==2.21
      - pyflakes==2.5.0
      - pygments==2.14.0
      - pyparsing==3.1.4
      - pyrsistent==0.18.0
      - pytest==7.0.1
      - python-dateutil==2.9.0.post0
      - pytz==2025.2
      - pyzmq==25.1.2
      - requests==2.27.1
      - scipy==1.5.4
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.2.0
      - snowballstemmer==2.2.0
      - sphinx==4.3.2
      - sphinxcontrib-applehelp==1.0.2
      - sphinxcontrib-devhelp==1.0.2
      - sphinxcontrib-htmlhelp==2.0.0
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==1.0.3
      - sphinxcontrib-serializinghtml==1.1.5
      - terminado==0.12.1
      - testpath==0.6.0
      - tomli==1.2.3
      - tornado==6.1
      - traitlets==4.3.3
      - typing-extensions==4.1.1
      - urllib3==1.26.20
      - wcwidth==0.2.13
      - webencodings==0.5.1
      - websocket-client==1.3.1
      - widgetsnbextension==3.6.10
      - zipp==3.6.0
prefix: /opt/conda/envs/PyBaMM

tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_average tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_boundary_value
tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_visualise
tests/unit/test_expression_tree/test_matrix.py::TestMatrix::test_array_wrapper tests/unit/test_expression_tree/test_matrix.py::TestMatrix::test_matrix_evaluate tests/unit/test_expression_tree/test_matrix.py::TestMatrix::test_matrix_modification tests/unit/test_expression_tree/test_matrix.py::TestMatrix::test_matrix_operations tests/unit/test_expression_tree/test_matrix.py::TestMatrix::test_matrix_simplifications tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_ghost_cell_flaf tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_has_spatial_derivatives tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_multiple_symbols tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_orphans tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_diff tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_domains tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_evaluates_to_number tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_evaluation tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_init tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_is_constant tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_methods tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_repr tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_simplify tests/unit/test_expression_tree/test_unary_operators.py::test_const_function tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_absolute tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_diff tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_function tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_function_simplify tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_gradient tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_id tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_integral tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_negation tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_printing tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_unary_operator
docs/source/expression_tree/unary_operator.rst pybamm/__init__.py pybamm/expression_tree/unary_operators.py pybamm/expression_tree/binary_operators.py
docs/source/expression_tree/unary_operator.rst pybamm/__init__.py pybamm/expression_tree/unary_operators.py pybamm/expression_tree/binary_operators.py