Repository: pybamm-team/PyBaMM

Issue:
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)
```

Discussion:

Patch:
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)


Test patch:
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")