Repository: pybamm-team/PyBaMM

Issue:
Remove Outer and Kron
These are now only used for broadcasting and can be replaced by matrix multiplication

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9a15d3259..a3ab4302a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -53,6 +53,7 @@
 
 ## Breaking changes
 
+-   Removed `Outer` and `Kron` nodes as no longer used ([#777](https://github.com/pybamm-team/PyBaMM/pull/777))
 -   Moved `results` to separate repositories ([#761](https://github.com/pybamm-team/PyBaMM/pull/761))
 -   The parameters "Bruggeman coefficient" must now be specified separately as "Bruggeman coefficient (electrolyte)" and "Bruggeman coefficient (electrode)"
 -   The current classes (`GetConstantCurrent`, `GetUserCurrent` and `GetUserData`) have now been removed. Please refer to the [`change-input-current` notebook](https://github.com/pybamm-team/PyBaMM/blob/master/examples/notebooks/change-input-current.ipynb) for information on how to specify an input current
diff --git a/docs/source/expression_tree/binary_operator.rst b/docs/source/expression_tree/binary_operator.rst
index 822864bbc..d8082d8c4 100644
--- a/docs/source/expression_tree/binary_operator.rst
+++ b/docs/source/expression_tree/binary_operator.rst
@@ -25,15 +25,7 @@ Binary Operators
 .. autoclass:: pybamm.Inner
   :members:
 
-.. autoclass:: pybamm.Outer
-  :members:
-
-.. autoclass:: pybamm.Kron
-  :members:
-
 .. autoclass:: pybamm.Heaviside
   :members:
 
-.. autofunction:: pybamm.outer
-
 .. autofunction:: pybamm.source
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index b47efcfa2..663540e3c 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -81,10 +81,7 @@ from .expression_tree.binary_operators import (
     Division,
     Inner,
     inner,
-    Outer,
-    Kron,
     Heaviside,
-    outer,
     source,
 )
 from .expression_tree.concatenations import (
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 1503e3c94..f99234378 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -989,7 +989,7 @@ class Discretisation(object):
         """
         Check variables in variable list against rhs
         Be lenient with size check if the variable in model.variables is broadcasted, or
-        a concatenation, or an outer product
+        a concatenation
         (if broadcasted, variable is a multiplication with a vector of ones)
         """
         for rhs_var in model.rhs.keys():
@@ -1001,7 +1001,6 @@ class Discretisation(object):
                 )
 
                 not_concatenation = not isinstance(var, pybamm.Concatenation)
-                not_outer = not isinstance(var, pybamm.Outer)
 
                 not_mult_by_one_vec = not (
                     isinstance(var, pybamm.Multiplication)
@@ -1012,7 +1011,6 @@ class Discretisation(object):
                 if (
                     different_shapes
                     and not_concatenation
-                    and not_outer
                     and not_mult_by_one_vec
                 ):
                     raise pybamm.ModelError(
diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index f5f2cb024..7765b3a95 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -5,7 +5,7 @@ import pybamm
 
 import numpy as np
 import numbers
-from scipy.sparse import issparse, csr_matrix, kron
+from scipy.sparse import issparse, csr_matrix
 
 
 def is_scalar_zero(expr):
@@ -79,17 +79,8 @@ class BinaryOperator(pybamm.Symbol):
     def __init__(self, name, left, right):
         left, right = self.format(left, right)
 
-        # Check and process domains, except for Outer symbol which takes the outer
-        # product of two smbols in different domains, and gives it the domain of the
-        # right child.
-        if isinstance(self, (pybamm.Outer, pybamm.Kron)):
-            domain = right.domain
-            auxiliary_domains = {}
-            if domain != []:
-                auxiliary_domains["secondary"] = left.domain
-        else:
-            domain = self.get_children_domains(left.domain, right.domain)
-            auxiliary_domains = self.get_children_auxiliary_domains([left, right])
+        domain = self.get_children_domains(left.domain, right.domain)
+        auxiliary_domains = self.get_children_auxiliary_domains([left, right])
         super().__init__(
             name,
             children=[left, right],
@@ -116,11 +107,7 @@ class BinaryOperator(pybamm.Symbol):
             )
 
         # Do some broadcasting in special cases, to avoid having to do this manually
-        if (
-            not isinstance(self, (Outer, Kron))
-            and left.domain != []
-            and right.domain != []
-        ):
+        if left.domain != [] and right.domain != []:
             if (
                 left.domain != right.domain
                 and "secondary" in right.auxiliary_domains
@@ -654,86 +641,6 @@ def inner(left, right):
     return pybamm.Inner(left, right)
 
 
-class Outer(BinaryOperator):
-    """A node in the expression tree representing an outer product.
-    This takes a 1D vector in the current collector domain of size (n,1) and a 1D
-    variable of size (m,1), takes their outer product, and reshapes this into a vector
-    of size (nm,1). It can also take in a vector in a single particle and a vector
-    of the electrolyte domain to repeat that particle.
-    Note: this class might be a bit dangerous, so at the moment it is very restrictive
-    in what symbols can be passed to it
-
-    **Extends:** :class:`BinaryOperator`
-    """
-
-    def __init__(self, left, right):
-        """ See :meth:`pybamm.BinaryOperator.__init__()`. """
-        # cannot have certain types of objects in the right symbol, as these
-        # can already be 2D objects (so we can't take an outer product with them)
-        if right.has_symbol_of_classes(
-            (pybamm.Variable, pybamm.StateVector, pybamm.Matrix, pybamm.SpatialVariable)
-        ):
-            raise TypeError("right child must only contain Vectors and Scalars" "")
-
-        super().__init__("outer product", left, right)
-
-    def __str__(self):
-        """ See :meth:`pybamm.Symbol.__str__()`. """
-        return "outer({!s}, {!s})".format(self.left, self.right)
-
-    def diff(self, variable):
-        """ See :meth:`pybamm.Symbol.diff()`. """
-        raise NotImplementedError("diff not implemented for symbol of type 'Outer'")
-
-    def _outer_jac(self, left_jac, right_jac, variable):
-        """
-        Calculate jacobian of outer product.
-        See :meth:`pybamm.Jacobian._jac()`.
-        """
-        # right cannot be a StateVector, so no need for product rule
-        left, right = self.orphans
-        if left.evaluates_to_number():
-            # Return zeros of correct size
-            return pybamm.Matrix(
-                csr_matrix((self.size, variable.evaluation_array.count(True)))
-            )
-        else:
-            return pybamm.Kron(left_jac, right)
-
-    def _binary_evaluate(self, left, right):
-        """ See :meth:`pybamm.BinaryOperator._binary_evaluate()`. """
-
-        return np.outer(left, right).reshape(-1, 1)
-
-
-class Kron(BinaryOperator):
-    """A node in the expression tree representing a (sparse) kronecker product operator
-
-    **Extends:** :class:`BinaryOperator`
-    """
-
-    def __init__(self, left, right):
-        """ See :meth:`pybamm.BinaryOperator.__init__()`. """
-
-        super().__init__("kronecker product", left, right)
-
-    def __str__(self):
-        """ See :meth:`pybamm.Symbol.__str__()`. """
-        return "kron({!s}, {!s})".format(self.left, self.right)
-
-    def diff(self, variable):
-        """ See :meth:`pybamm.Symbol.diff()`. """
-        raise NotImplementedError("diff not implemented for symbol of type 'Kron'")
-
-    def _binary_jac(self, left_jac, right_jac):
-        """ See :meth:`pybamm.BinaryOperator._binary_jac()`. """
-        raise NotImplementedError("jac not implemented for symbol of type 'Kron'")
-
-    def _binary_evaluate(self, left, right):
-        """ See :meth:`pybamm.BinaryOperator._binary_evaluate()`. """
-        return csr_matrix(kron(left, right))
-
-
 class Heaviside(BinaryOperator):
     """A node in the expression tree representing a heaviside step function.
 
@@ -783,18 +690,6 @@ class Heaviside(BinaryOperator):
         return Heaviside(left, right, self.equal)
 
 
-def outer(left, right):
-    """
-    Return outer product of two symbols. If the symbols have the same domain, the outer
-    product is just a multiplication. If they have different domains, make a copy of the
-    left child with same domain as right child, and then take outer product.
-    """
-    try:
-        return left * right
-    except pybamm.DomainError:
-        return pybamm.Outer(left, right)
-
-
 def source(left, right, boundary=False):
     """A convinience function for creating (part of) an expression tree representing
     a source term. This is necessary for spatial methods where the mass matrix
diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py
index 983008d70..706960d3c 100644
--- a/pybamm/expression_tree/operations/convert_to_casadi.py
+++ b/pybamm/expression_tree/operations/convert_to_casadi.py
@@ -59,11 +59,8 @@ class CasadiConverter(object):
             # process children
             converted_left = self.convert(left, t, y, u)
             converted_right = self.convert(right, t, y, u)
-            if isinstance(symbol, pybamm.Outer):
-                return casadi.kron(converted_left, converted_right)
-            else:
-                # _binary_evaluate defined in derived classes for specific rules
-                return symbol._binary_evaluate(converted_left, converted_right)
+            # _binary_evaluate defined in derived classes for specific rules
+            return symbol._binary_evaluate(converted_left, converted_right)
 
         elif isinstance(symbol, pybamm.UnaryOperator):
             converted_child = self.convert(symbol.child, t, y, u)
diff --git a/pybamm/expression_tree/operations/evaluate.py b/pybamm/expression_tree/operations/evaluate.py
index 6558357c3..a11769b06 100644
--- a/pybamm/expression_tree/operations/evaluate.py
+++ b/pybamm/expression_tree/operations/evaluate.py
@@ -92,14 +92,6 @@ def find_symbols(symbol, constant_symbols, variable_symbols):
                 "if scipy.sparse.issparse({1}) else "
                 "{0} * {1}".format(children_vars[0], children_vars[1])
             )
-        elif isinstance(symbol, pybamm.Outer):
-            symbol_str = "np.outer({}, {}).reshape(-1, 1)".format(
-                children_vars[0], children_vars[1]
-            )
-        elif isinstance(symbol, pybamm.Kron):
-            symbol_str = "scipy.sparse.csr_matrix(scipy.sparse.kron({}, {}))".format(
-                children_vars[0], children_vars[1]
-            )
         else:
             symbol_str = children_vars[0] + " " + symbol.name + " " + children_vars[1]
 
diff --git a/pybamm/expression_tree/operations/jacobian.py b/pybamm/expression_tree/operations/jacobian.py
index 629d498b9..42bd3683d 100644
--- a/pybamm/expression_tree/operations/jacobian.py
+++ b/pybamm/expression_tree/operations/jacobian.py
@@ -46,15 +46,8 @@ class Jacobian(object):
             # process children
             left_jac = self.jac(left, variable)
             right_jac = self.jac(right, variable)
-            # Need to treat outer differently. If the left child of an Outer
-            # evaluates to number then we need to return a matrix of zeros
-            # of the correct size, which requires variable.evaluation_array
-            if isinstance(symbol, pybamm.Outer):
-                # _outer_jac defined in pybamm.Outer
-                jac = symbol._outer_jac(left_jac, right_jac, variable)
-            else:
-                # _binary_jac defined in derived classes for specific rules
-                jac = symbol._binary_jac(left_jac, right_jac)
+            # _binary_jac defined in derived classes for specific rules
+            jac = symbol._binary_jac(left_jac, right_jac)
 
         elif isinstance(symbol, pybamm.UnaryOperator):
             child_jac = self.jac(symbol.child, variable)
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index ec9ee9433..1976b806c 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -98,9 +98,15 @@ class SpatialMethod:
         )
 
         if broadcast_type == "primary":
-            out = pybamm.Outer(
-                symbol, pybamm.Vector(np.ones(primary_domain_size), domain=domain)
-            )
+            # Make copies of the child stacked on top of each other
+            sub_vector = np.ones((primary_domain_size, 1))
+            if symbol.shape_for_testing == ():
+                out = symbol * pybamm.Vector(sub_vector)
+            else:
+                # Repeat for secondary points
+                matrix = csr_matrix(kron(eye(symbol.shape_for_testing[0]), sub_vector))
+                out = pybamm.Matrix(matrix) @ symbol
+            out.domain = domain
         elif broadcast_type == "secondary":
             secondary_domain_size = sum(
                 self.mesh[dom][0].npts_for_broadcast


Test patch:
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 8ff690c17..5e823c2a0 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -781,13 +781,18 @@ class TestDiscretise(unittest.TestCase):
 
         disc.set_variable_slices([var])
         broad_disc = disc.process_symbol(broad)
-        self.assertIsInstance(broad_disc, pybamm.Outer)
-        self.assertIsInstance(broad_disc.children[0], pybamm.StateVector)
-        self.assertIsInstance(broad_disc.children[1], pybamm.Vector)
+        self.assertIsInstance(broad_disc, pybamm.MatrixMultiplication)
+        self.assertIsInstance(broad_disc.children[0], pybamm.Matrix)
+        self.assertIsInstance(broad_disc.children[1], pybamm.StateVector)
         self.assertEqual(
             broad_disc.shape,
             (mesh["separator"][0].npts * mesh["current collector"][0].npts, 1),
         )
+        y_test = np.linspace(0, 1, mesh["current collector"][0].npts)
+        np.testing.assert_array_equal(
+            broad_disc.evaluate(y=y_test),
+            np.outer(y_test, np.ones(mesh["separator"][0].npts)).reshape(-1, 1),
+        )
 
     def test_secondary_broadcast_2D(self):
         # secondary broadcast in 2D --> Matrix multiplication
@@ -808,29 +813,6 @@ class TestDiscretise(unittest.TestCase):
             (mesh["negative particle"][0].npts * mesh["negative electrode"][0].npts, 1),
         )
 
-    def test_outer(self):
-
-        # create discretisation
-        disc = get_1p1d_discretisation_for_testing()
-        mesh = disc.mesh
-
-        var_z = pybamm.Variable("var_z", ["current collector"])
-        var_x = pybamm.Vector(
-            np.linspace(0, 1, mesh["separator"][0].npts), domain="separator"
-        )
-
-        # process Outer variable
-        disc.set_variable_slices([var_z, var_x])
-        outer = pybamm.outer(var_z, var_x)
-        outer_disc = disc.process_symbol(outer)
-        self.assertIsInstance(outer_disc, pybamm.Outer)
-        self.assertIsInstance(outer_disc.children[0], pybamm.StateVector)
-        self.assertIsInstance(outer_disc.children[1], pybamm.Vector)
-        self.assertEqual(
-            outer_disc.shape,
-            (mesh["separator"][0].npts * mesh["current collector"][0].npts, 1),
-        )
-
     def test_concatenation(self):
         a = pybamm.Symbol("a")
         b = pybamm.Symbol("b")
diff --git a/tests/unit/test_expression_tree/test_binary_operators.py b/tests/unit/test_expression_tree/test_binary_operators.py
index 2381bcab6..8127b41b2 100644
--- a/tests/unit/test_expression_tree/test_binary_operators.py
+++ b/tests/unit/test_expression_tree/test_binary_operators.py
@@ -58,54 +58,6 @@ class TestBinaryOperators(unittest.TestCase):
         pow2 = pybamm.Power(a, b)
         self.assertEqual(pow2.evaluate(), 16)
 
-    def test_outer(self):
-        # Outer class
-        v = pybamm.Vector(np.ones(5), domain="current collector")
-        w = pybamm.Vector(2 * np.ones(3), domain="test")
-        outer = pybamm.Outer(v, w)
-        np.testing.assert_array_equal(outer.evaluate(), 2 * np.ones((15, 1)))
-        self.assertEqual(outer.domain, w.domain)
-        self.assertEqual(
-            str(outer), "outer(Column vector of length 5, Column vector of length 3)"
-        )
-
-        # outer function
-        # if there is no domain clash, normal multiplication is retured
-        u = pybamm.Vector(np.linspace(0, 1, 5))
-        outer = pybamm.outer(u, v)
-        self.assertIsInstance(outer, pybamm.Multiplication)
-        np.testing.assert_array_equal(outer.evaluate(), u.evaluate())
-        # otherwise, Outer class is returned
-        outer_fun = pybamm.outer(v, w)
-        outer_class = pybamm.Outer(v, w)
-        self.assertEqual(outer_fun.id, outer_class.id)
-
-        # failures
-        y = pybamm.StateVector(slice(10))
-        with self.assertRaisesRegex(
-            TypeError, "right child must only contain Vectors and Scalars"
-        ):
-            pybamm.Outer(v, y)
-        with self.assertRaises(NotImplementedError):
-            outer_fun.diff(None)
-
-    def test_kron(self):
-        # Kron class
-        A = pybamm.Matrix(np.eye(2))
-        b = pybamm.Vector(np.array([[4], [5]]))
-        kron = pybamm.Kron(A, b)
-        np.testing.assert_array_equal(
-            kron.evaluate().toarray(), np.kron(A.entries, b.entries)
-        )
-
-        # failures
-        with self.assertRaises(NotImplementedError):
-            kron.diff(None)
-
-        y = pybamm.StateVector(slice(0, 2))
-        with self.assertRaises(NotImplementedError):
-            kron.jac(y)
-
     def test_known_eval(self):
         # Scalars
         a = pybamm.Scalar(4)
diff --git a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
index d0ee4fcf1..a5fc04748 100644
--- a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
+++ b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
@@ -71,12 +71,6 @@ class TestCasadiConverter(unittest.TestCase):
         # State Vector
         self.assert_casadi_equal(pybamm_y.to_casadi(casadi_t, casadi_y), casadi_y)
 
-        # outer product
-        outer = pybamm.Outer(pybamm_a, pybamm_a)
-        self.assert_casadi_equal(
-            outer.to_casadi(), casadi.MX(outer.evaluate()), evalf=True
-        )
-
     def test_special_functions(self):
         a = pybamm.Array(np.array([1, 2, 3, 4, 5]))
         self.assert_casadi_equal(pybamm.max(a).to_casadi(), casadi.MX(5), evalf=True)
diff --git a/tests/unit/test_expression_tree/test_operations/test_evaluate.py b/tests/unit/test_expression_tree/test_operations/test_evaluate.py
index 52addcb63..cd6aae7dd 100644
--- a/tests/unit/test_expression_tree/test_operations/test_evaluate.py
+++ b/tests/unit/test_expression_tree/test_operations/test_evaluate.py
@@ -391,14 +391,6 @@ class TestEvaluate(unittest.TestCase):
             result = evaluator.evaluate(t=t, y=y).toarray()
             np.testing.assert_allclose(result, expr.evaluate(t=t, y=y).toarray())
 
-        # test Outer
-        v = pybamm.Vector(np.ones(5), domain="current collector")
-        w = pybamm.Vector(2 * np.ones(3), domain="test")
-        expr = pybamm.Outer(v, w)
-        evaluator = pybamm.EvaluatorPython(expr)
-        result = evaluator.evaluate()
-        np.testing.assert_allclose(result, expr.evaluate())
-
         # test Inner
         v = pybamm.Vector(np.ones(5), domain="test")
         w = pybamm.Vector(2 * np.ones(5), domain="test")
diff --git a/tests/unit/test_expression_tree/test_operations/test_jac.py b/tests/unit/test_expression_tree/test_operations/test_jac.py
index 526f31390..270b030c4 100644
--- a/tests/unit/test_expression_tree/test_operations/test_jac.py
+++ b/tests/unit/test_expression_tree/test_operations/test_jac.py
@@ -56,32 +56,6 @@ class TestJacobian(unittest.TestCase):
         du_dv = u.jac(v).evaluate().toarray()
         np.testing.assert_array_equal(du_dv, jacobian)
 
-        # test Jacobian of Outer (must set domain to be 'current collector')
-        u.domain = ["current collector"]
-        func = pybamm.Outer(u, pybamm.Scalar(4))
-        jacobian = np.array([[4, 0, 0, 0], [0, 4, 0, 0]])
-        dfunc_dy = func.jac(y).evaluate(y=y0)
-        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
-
-        func = pybamm.Outer(u, pybamm.Vector(np.array([1, 2, 3])))
-        jacobian = np.array(
-            [
-                [1, 0, 0, 0],
-                [2, 0, 0, 0],
-                [3, 0, 0, 0],
-                [0, 1, 0, 0],
-                [0, 2, 0, 0],
-                [0, 3, 0, 0],
-            ]
-        )
-        dfunc_dy = func.jac(y).evaluate(y=y0)
-        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
-
-        # test jac of outer if left evaluates to number
-        func = pybamm.Outer(pybamm.Scalar(1), pybamm.Scalar(4))
-        dfunc_dy = func.jac(y).evaluate(y=y0)
-        np.testing.assert_array_equal(0, dfunc_dy.toarray())
-
     def test_nonlinear(self):
         y = pybamm.StateVector(slice(0, 4))
         u = pybamm.StateVector(slice(0, 2))
diff --git a/tests/unit/test_expression_tree/test_operations/test_jac_2D.py b/tests/unit/test_expression_tree/test_operations/test_jac_2D.py
index 0d3679715..d08e2b313 100644
--- a/tests/unit/test_expression_tree/test_operations/test_jac_2D.py
+++ b/tests/unit/test_expression_tree/test_operations/test_jac_2D.py
@@ -86,40 +86,6 @@ class TestJacobian(unittest.TestCase):
         with self.assertRaises(NotImplementedError):
             u.jac(v)
 
-        # test Jacobian of Outer (must set domain to be 'current collector')
-        u.domain = ["current collector"]
-        func = pybamm.Outer(u, pybamm.Scalar(4))
-        jacobian = np.array(
-            [
-                [4, 0, 0, 0, 0, 0, 0, 0],
-                [0, 4, 0, 0, 0, 0, 0, 0],
-                [0, 0, 0, 0, 4, 0, 0, 0],
-                [0, 0, 0, 0, 0, 4, 0, 0],
-            ]
-        )
-        dfunc_dy = func.jac(y).evaluate(y=y0)
-        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
-
-        func = pybamm.Outer(u, pybamm.Vector(np.array([1, 2, 3])))
-        jacobian = np.array(
-            [
-                [1, 0, 0, 0, 0, 0, 0, 0],
-                [2, 0, 0, 0, 0, 0, 0, 0],
-                [3, 0, 0, 0, 0, 0, 0, 0],
-                [0, 1, 0, 0, 0, 0, 0, 0],
-                [0, 2, 0, 0, 0, 0, 0, 0],
-                [0, 3, 0, 0, 0, 0, 0, 0],
-                [0, 0, 0, 0, 1, 0, 0, 0],
-                [0, 0, 0, 0, 2, 0, 0, 0],
-                [0, 0, 0, 0, 3, 0, 0, 0],
-                [0, 0, 0, 0, 0, 1, 0, 0],
-                [0, 0, 0, 0, 0, 2, 0, 0],
-                [0, 0, 0, 0, 0, 3, 0, 0],
-            ]
-        )
-        dfunc_dy = func.jac(y).evaluate(y=y0)
-        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
-
     def test_nonlinear(self):
         y = pybamm.StateVector(slice(0, 8))
         u = pybamm.StateVector(slice(0, 2), slice(4, 6))
@@ -251,21 +217,16 @@ class TestJacobian(unittest.TestCase):
         a_dom = ["negative electrode"]
         b_dom = ["separator"]
         c_dom = ["positive electrode"]
-        a_npts = mesh[a_dom[0]][0].npts
-        b_npts = mesh[b_dom[0]][0].npts
-        c_npts = mesh[c_dom[0]][0].npts
         cc_npts = mesh["current collector"][0].npts
         curr_coll_vector = pybamm.Vector(np.ones(cc_npts), domain="current collector")
-        a = 2 * pybamm.Outer(
-            curr_coll_vector, pybamm.Vector(np.ones(a_npts), domain=a_dom)
-        )
-        b = pybamm.Outer(curr_coll_vector, pybamm.Vector(np.ones(b_npts), domain=b_dom))
-        c = 3 * pybamm.Outer(
-            curr_coll_vector, pybamm.Vector(np.ones(c_npts), domain=c_dom)
-        )
+        a = 2 * pybamm.PrimaryBroadcast(curr_coll_vector, a_dom)
+        b = pybamm.PrimaryBroadcast(curr_coll_vector, b_dom)
+        c = 3 * pybamm.PrimaryBroadcast(curr_coll_vector, c_dom)
 
-        conc = pybamm.DomainConcatenation([a, b, c], mesh)
-        jac = conc.jac(y).evaluate().toarray()
+        conc = pybamm.Concatenation(a, b, c)
+        disc.set_variable_slices([conc])
+        conc_disc = disc.process_symbol(conc)
+        jac = conc_disc.jac(y).evaluate().toarray()
         np.testing.assert_array_equal(jac, np.zeros((1500, 1500)))
 
         # Jacobian of a DomainConcatenation of StateVectors
diff --git a/tests/unit/test_expression_tree/test_operations/test_simplify.py b/tests/unit/test_expression_tree/test_operations/test_simplify.py
index 9625a090d..9c47da440 100644
--- a/tests/unit/test_expression_tree/test_operations/test_simplify.py
+++ b/tests/unit/test_expression_tree/test_operations/test_simplify.py
@@ -570,40 +570,6 @@ class TestSimplify(unittest.TestCase):
         broad_simp = broad.simplify()
         self.assertEqual(broad_simp.id, broad.id)
 
-    def test_simplify_outer(self):
-        v = pybamm.Vector(np.ones(5), domain="current collector")
-        w = pybamm.Vector(2 * np.ones(3), domain="test")
-        outer_simp = pybamm.Outer(v, w).simplify()
-        self.assertIsInstance(outer_simp, pybamm.Vector)
-        np.testing.assert_array_equal(outer_simp.evaluate(), 2 * np.ones((15, 1)))
-
-    def test_simplify_divide_outer(self):
-        u = pybamm.Scalar(1)
-        v = pybamm.StateVector(slice(0, 5), domain="current collector")
-        outer = pybamm.Outer(v, u)
-
-        exp1 = pybamm.Division(pybamm.Division(outer, u), u)
-        self.assertIsInstance(exp1.simplify(), pybamm.Outer)
-
-        exp2 = pybamm.Division(pybamm.Division(outer, 2 * u), u)
-        self.assertIsInstance(exp2.simplify(), pybamm.Multiplication)
-
-        exp3 = pybamm.Division(pybamm.Division(outer, u), 2 * u)
-        self.assertIsInstance(exp3.simplify(), pybamm.Multiplication)
-
-        exp4 = pybamm.Division(pybamm.Division(outer, 2 * u), 2 * u)
-        self.assertIsInstance(exp4.simplify(), pybamm.Multiplication)
-
-    def test_simplify_kron(self):
-        A = pybamm.Matrix(np.eye(2))
-        b = pybamm.Vector(np.array([[4], [5]]))
-        kron = pybamm.Kron(A, b)
-        kron_simp = kron.simplify()
-        self.assertIsInstance(kron_simp, pybamm.Matrix)
-        np.testing.assert_array_equal(
-            kron_simp.evaluate().toarray(), np.kron(A.entries, b.entries)
-        )
-
     def test_simplify_heaviside(self):
         a = pybamm.Scalar(1)
         b = pybamm.Scalar(2)
diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py
index 75bc3a4a7..a41b3e5ea 100644
--- a/tests/unit/test_parameters/test_parameter_values.py
+++ b/tests/unit/test_parameters/test_parameter_values.py
@@ -235,12 +235,6 @@ class TestParameterValues(unittest.TestCase):
             processed_g.evaluate(y=np.ones(10)), np.ones((10, 1))
         )
 
-        # process outer
-        c = pybamm.Parameter("c", domain="current collector")
-        outer = pybamm.Outer(c, b)
-        processed_outer = parameter_values.process_symbol(outer)
-        self.assertIsInstance(processed_outer, pybamm.Outer)
-
         # not implemented
         sym = pybamm.Symbol("sym")
         with self.assertRaises(NotImplementedError):