Repository: pybamm-team/PyBaMM (python)

Issue:
Allow for different tab locations in 1+1D
At present the 1plus1D current collector model assumes both tabs are at the top (regardless of the tab locations set in the parameters). Also, in 2plus1D we use "left" and "right" to correspond to the negative and positive tab, respectively. It would be cleaner to be able to set BCs on "negative tab" and "positive tab". It would also be useful to then have something like `pybamm.BoundaryValue(var, "negative tab")` which can be used to get e.g. terminal voltage.


Pull request description:
604 allow BCs set on neg and pos tabs
# Description
Allows setting boundary conditions and taking boundary value on tabs. Currently this is implemented for 1D meshes by having `discretisation` change "negative tab" and "positive tab" to "left" or "right" depending on the tab location, which is now stored in the mesh. The benefit of doing it this way is that when implementing any new 1D spatial methods you just deal with left and right, and don't worry about catching tabs too. There may be benefits to doing this in a different way though? In 2D the spatial method deals directly with the strings  "negative tab" and "positive tab" since "left" and "right" don't really make sense in this context, and were always taken to mean the tabs. 

What are your thoughts on this @Scottmar93 @tinosulzer @martinjrobins ?

Fixes #604 

## Type of change

Please delete options that are not relevant.

- [ ] Bug fix (non-breaking change which fixes an issue)
- [ ] New feature (non-breaking change which adds functionality)
- [ ] Breaking change (fix or feature that would cause existing functionality to not work as expected)


# Key checklist:

- [ ] No style issues: `$ flake8`
- [ ] All tests pass: `$ python run-tests.py --unit`
- [ ] The documentation builds: `$ cd docs` and then `$ make clean; make html`

You can run all three at once, using `$ python run-tests.py --quick`.

## Further checks: 

- [ ] Code is commented, particularly in hard-to-understand areas
- [ ] Tests added that prove fix is effective or that feature works
- [ ] Any dependent changes have been merged and published in downstream modules


Patch:
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 1e63ed9fe..65c8f71d8 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -324,6 +324,12 @@ class Discretisation(object):
         # in discrisation of other boundary conditions
         for key, bcs in model.boundary_conditions.items():
             processed_bcs[key.id] = {}
+
+            # Handle any boundary conditions applied on the tabs
+            if any("tab" in side for side in list(bcs.keys())):
+                bcs = self.check_tab_conditions(key, bcs)
+
+            # Process boundary conditions
             for side, bc in bcs.items():
                 eqn, typ = bc
                 pybamm.logger.debug("Discretise {} ({} bc)".format(key, side))
@@ -332,15 +338,55 @@ class Discretisation(object):
 
         return processed_bcs
 
-    def _process_bc_entry(self, key, bcs):
-        processed_entry = {key.id: {}}
-        for side, bc in bcs.items():
-            eqn, typ = bc
-            pybamm.logger.debug("Discretise {} ({} bc)".format(key, side))
-            processed_eqn = self.process_symbol(eqn)
-            processed_entry[key.id][side] = (processed_eqn, typ)
+    def check_tab_conditions(self, symbol, bcs):
+        """
+        Check any boundary conditions applied on "negative tab", "positive tab"
+        and "no tab". For 1D current collector meshes, these conditions are
+        converted into boundary conditions on "left" (tab at z=0) or "right"
+        (tab at z=l_z) depending on the tab location stored in the mesh. For 2D
+        current collector meshes, the boundary conditions can be applied on the
+        tabs directly.
+
+        Parameters
+        ----------
+        symbol : :class:`pybamm.expression_tree.symbol.Symbol`
+            The symbol on which the boundary conditions are applied.
+        bcs : dict
+            The dictionary of boundary conditions (a dict of {side: equation}).
 
-        return processed_entry
+        Returns
+        -------
+        dict
+            The dictionary of boundary conditions, with the keys changed to
+            "left" and "right" where necessary.
+
+        """
+        # Check symbol domain
+        domain = symbol.domain[0]
+        mesh = self.mesh[domain][0]
+
+        if domain != "current collector":
+            raise pybamm.ModelError(
+                """Boundary conditions can only be applied on the tabs in the domain
+            'current collector', but {} has domain {}""".format(
+                    symbol, domain
+                )
+            )
+
+        # Replace keys with "left" and "right" as appropriate for 1D meshes
+        if isinstance(mesh, pybamm.SubMesh1D):
+            # replace negative and/or positive tab
+            for tab in ["negative tab", "positive tab"]:
+                if any(tab in side for side in list(bcs.keys())):
+                    bcs[mesh.tabs[tab]] = bcs.pop(tab)
+            # replace no tab
+            if any("no tab" in side for side in list(bcs.keys())):
+                if "left" in list(bcs.keys()):
+                    bcs["right"] = bcs.pop("no tab")  # tab at bottom
+                else:
+                    bcs["left"] = bcs.pop("no tab")  # tab at top
+
+        return bcs
 
     def process_rhs_and_algebraic(self, model):
         """Discretise model equations - differential ('rhs') and algebraic.
@@ -502,6 +548,13 @@ class Discretisation(object):
 
         if symbol.domain != []:
             spatial_method = self.spatial_methods[symbol.domain[0]]
+            # If boundary conditions are provided, need to check for BCs on tabs
+            if self.bcs:
+                key_id = list(self.bcs.keys())[0]
+                if any("tab" in side for side in list(self.bcs[key_id].keys())):
+                    self.bcs[key_id] = self.check_tab_conditions(
+                        symbol, self.bcs[key_id]
+                    )
 
         if isinstance(symbol, pybamm.BinaryOperator):
             # Pre-process children
@@ -520,6 +573,7 @@ class Discretisation(object):
             disc_child = self.process_symbol(child)
             if child.domain != []:
                 child_spatial_method = self.spatial_methods[child.domain[0]]
+
             if isinstance(symbol, pybamm.Gradient):
                 return child_spatial_method.gradient(child, disc_child, self.bcs)
 
@@ -574,6 +628,13 @@ class Discretisation(object):
                 return symbol
 
             elif isinstance(symbol, pybamm.BoundaryOperator):
+                # if boundary operator applied on "negative tab" or
+                # "positive tab" *and* the mesh is 1D then change side to
+                # "left" or "right" as appropriate
+                if symbol.side in ["negative tab", "positive tab"]:
+                    mesh = self.mesh[symbol.children[0].domain[0]][0]
+                    if isinstance(mesh, pybamm.SubMesh1D):
+                        symbol.side = mesh.tabs[symbol.side]
                 return child_spatial_method.boundary_value_or_flux(symbol, disc_child)
 
             else:
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index fd3be4c75..e8b24dc8f 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -635,8 +635,18 @@ class BoundaryOperator(SpatialOperator):
     """
 
     def __init__(self, name, child, side):
+        # side can only be "negative tab" or "positive tab" if domain is
+        # "current collector"
+        if side in ["negative tab", "positive tab"]:
+            if child.domain[0] != "current collector":
+                raise pybamm.ModelError(
+                    """Can only take boundary value on the tabs in the domain
+                'current collector', but {} has domain {}""".format(
+                        child, child.domain[0]
+                    )
+                )
         self.side = side
-        # integral of a child takes the domain from auxiliary domain of the child
+        # boundary value of a child takes the domain from auxiliary domain of the child
         if child.auxiliary_domains != {}:
             domain = child.auxiliary_domains["secondary"]
         # if child has no auxiliary domain, integral removes domain
diff --git a/pybamm/meshes/meshes.py b/pybamm/meshes/meshes.py
index 4b144f33a..e4f79ff33 100644
--- a/pybamm/meshes/meshes.py
+++ b/pybamm/meshes/meshes.py
@@ -75,7 +75,7 @@ class Mesh(dict):
         # evaluate any expressions in geometry
         for domain in geometry:
             for prim_sec_tabs, variables in geometry[domain].items():
-                # process tab information if using 2D current collectors
+                # process tab information if using 1 or 2D current collectors
                 if prim_sec_tabs == "tabs":
                     for tab, position_size in variables.items():
                         for position_size, sym in position_size.items():
@@ -95,10 +95,11 @@ class Mesh(dict):
 
         # Create submeshes
         for domain in geometry:
-            # need to pass tab information if primary domian is 2D current collector
+            # need to pass tab information if primary domain is 1 or 2D
+            # current collector
             if (
                 domain == "current collector"
-                and submesh_types[domain] == pybamm.Scikit2DSubMesh
+                and submesh_types[domain] != pybamm.SubMesh0D
             ):
                 self[domain] = [
                     submesh_types[domain](
diff --git a/pybamm/meshes/one_dimensional_submeshes.py b/pybamm/meshes/one_dimensional_submeshes.py
index 7f1906b4d..e26906457 100644
--- a/pybamm/meshes/one_dimensional_submeshes.py
+++ b/pybamm/meshes/one_dimensional_submeshes.py
@@ -9,10 +9,11 @@ import numpy as np
 class SubMesh1D:
     """
     1D submesh class.
-    Contains the position of the nodes and the number of mesh points.
+    Contains the position of the nodes, the number of mesh points, and
+    (optionally) information about the tab locations.
     """
 
-    def __init__(self, edges, coord_sys):
+    def __init__(self, edges, coord_sys, tabs=None):
         self.edges = edges
         self.nodes = (self.edges[1:] + self.edges[:-1]) / 2
         self.d_edges = np.diff(self.edges)
@@ -20,6 +21,20 @@ class SubMesh1D:
         self.npts = self.nodes.size
         self.coord_sys = coord_sys
 
+        # Add tab locations in terms of "left" and "right"
+        if tabs:
+            self.tabs = {}
+            l_z = self.edges[-1]
+
+            def near(x, point, tol=3e-16):
+                return abs(x - point) < tol
+
+            for tab in ["negative", "positive"]:
+                if near(tabs[tab]["z_centre"], 0):
+                    self.tabs[tab + " tab"] = "left"
+                elif near(tabs[tab]["z_centre"], l_z):
+                    self.tabs[tab + " tab"] = "right"
+
 
 class Uniform1DSubMesh(SubMesh1D):
     """
@@ -32,9 +47,12 @@ class Uniform1DSubMesh(SubMesh1D):
     npts : dict
         A dictionary that contains the number of points to be used on each
         spatial variable
+    tabs : dict, optional
+        A dictionary that contains information about the size and location of
+        the tabs
     """
 
-    def __init__(self, lims, npts):
+    def __init__(self, lims, npts, tabs=None):
 
         # currently accept lims and npts as dicts. This may get changed at a future
         # date depending on the form of mesh we desire for 2D/3D
@@ -51,4 +69,4 @@ class Uniform1DSubMesh(SubMesh1D):
 
         coord_sys = spatial_var.coord_sys
 
-        super().__init__(edges, coord_sys=coord_sys)
+        super().__init__(edges, coord_sys=coord_sys, tabs=tabs)
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 1cb138ff7..bec0eb93a 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -639,26 +639,15 @@ class BaseBatteryModel(pybamm.BaseModel):
         eta_r_av = eta_r_p_av - eta_r_n_av
         eta_r_av_dim = eta_r_p_av_dim - eta_r_n_av_dim
 
-        # terminal voltage
-        phi_s_cn = self.variables["Negative current collector potential"]
+        # terminal voltage (Note: phi_s_cn is zero at the negative tab)
         phi_s_cp = self.variables["Positive current collector potential"]
-        phi_s_cn_dim = self.variables["Negative current collector potential [V]"]
         phi_s_cp_dim = self.variables["Positive current collector potential [V]"]
         if self.options["dimensionality"] == 0:
             V = phi_s_cp
             V_dim = phi_s_cp_dim
-        elif self.options["dimensionality"] == 1:
-            # In 1D both tabs are at "right"
-            V = pybamm.BoundaryValue(phi_s_cp, "right")
-            V_dim = pybamm.BoundaryValue(phi_s_cp_dim, "right")
-        elif self.options["dimensionality"] == 2:
-            # In 2D left corresponds to the negative tab and right the positive tab
-            V = pybamm.BoundaryValue(phi_s_cp, "right") - pybamm.BoundaryValue(
-                phi_s_cn, "left"
-            )
-            V_dim = pybamm.BoundaryValue(phi_s_cp_dim, "right") - pybamm.BoundaryValue(
-                phi_s_cn_dim, "left"
-            )
+        elif self.options["dimensionality"] in [1, 2]:
+            V = pybamm.BoundaryValue(phi_s_cp, "positive tab")
+            V_dim = pybamm.BoundaryValue(phi_s_cp_dim, "positive tab")
 
         # TODO: add current collector losses to the voltage in 3D
 
diff --git a/pybamm/models/submodels/current_collector/effective_resistance_current_collector.py b/pybamm/models/submodels/current_collector/effective_resistance_current_collector.py
index a49872a3f..87f234939 100644
--- a/pybamm/models/submodels/current_collector/effective_resistance_current_collector.py
+++ b/pybamm/models/submodels/current_collector/effective_resistance_current_collector.py
@@ -68,10 +68,13 @@ class EffectiveResistance2D(pybamm.BaseModel):
 
         self.boundary_conditions = {
             psi: {
-                "left": (psi_neg_tab_bc, "Neumann"),
-                "right": (psi_pos_tab_bc, "Neumann"),
+                "negative tab": (psi_neg_tab_bc, "Neumann"),
+                "positive tab": (psi_pos_tab_bc, "Neumann"),
+            },
+            W: {
+                "negative tab": (W_neg_tab_bc, "Neumann"),
+                "positive tab": (W_pos_tab_bc, "Neumann")
             },
-            W: {"left": (W_neg_tab_bc, "Neumann"), "right": (W_pos_tab_bc, "Neumann")},
         }
 
         # "Initial conditions" provides initial guess for solver
@@ -84,10 +87,10 @@ class EffectiveResistance2D(pybamm.BaseModel):
         }
 
         # Define effective current collector resistance
-        psi_neg_tab = pybamm.BoundaryValue(psi, "left")
-        psi_pos_tab = pybamm.BoundaryValue(psi, "right")
-        W_neg_tab = pybamm.BoundaryValue(W, "left")
-        W_pos_tab = pybamm.BoundaryValue(W, "right")
+        psi_neg_tab = pybamm.BoundaryValue(psi, "negative tab")
+        psi_pos_tab = pybamm.BoundaryValue(psi, "positive tab")
+        W_neg_tab = pybamm.BoundaryValue(W, "negative tab")
+        W_pos_tab = pybamm.BoundaryValue(W, "positive tab")
 
         R_cc = (
             (alpha_prime / l_y)
diff --git a/pybamm/models/submodels/current_collector/potential_pair.py b/pybamm/models/submodels/current_collector/potential_pair.py
index efcd7b8b9..08e9304e1 100644
--- a/pybamm/models/submodels/current_collector/potential_pair.py
+++ b/pybamm/models/submodels/current_collector/potential_pair.py
@@ -98,12 +98,12 @@ class PotentialPair1plus1D(BasePotentialPair):
         # even though phi_s_cp isn't a pybamm.Variable object
         self.boundary_conditions = {
             phi_s_cn: {
-                "left": (pybamm.Scalar(0), "Neumann"),
-                "right": (pybamm.Scalar(0), "Dirichlet"),
+                "negative tab": (pybamm.Scalar(0), "Dirichlet"),
+                "no tab": (pybamm.Scalar(0), "Neumann"),
             },
             phi_s_cp: {
-                "left": (pybamm.Scalar(0), "Neumann"),
-                "right": (pos_tab_bc, "Neumann"),
+                "no tab": (pybamm.Scalar(0), "Neumann"),
+                "positive tab": (pos_tab_bc, "Neumann"),
             },
         }
 
@@ -145,14 +145,29 @@ class PotentialPair2plus1D(BasePotentialPair):
 
         # Boundary condition needs to be on the variables that go into the Laplacian,
         # even though phi_s_cp isn't a pybamm.Variable object
+        # In the 2+1D model, the equations for the current collector potentials
+        # are solved on a 2D domain and the regions "negative tab" and "positive tab"
+        # are the projections of the tabs onto this 2D domain.
+        # In the 2D formulation it is assumed that no flux boundary conditions
+        # are applied everywhere apart from the tabs.
+        # The reference potenital is taken to be zero on the negative tab,
+        # giving the zero Dirichlet condition on phi_s_cn. Elsewhere, the boundary
+        # is insulated, giving no flux conditions on phi_s_cn. This is automatically
+        # applied everywhere, apart from the region corresponding to the projection
+        # of the positive tab, so we need to explititly apply a zero-flux boundary
+        # condition on the region "positive tab" for phi_s_cn.
+        # A current is drawn from the positive tab, giving the non-zero Neumann
+        # boundary condition on phi_s_cp at "positive tab". Elsewhere, the boundary is
+        # insulated, so, as with phi_s_cn, we need to explicitly give the zero-flux
+        # condition on the region "negative tab" for phi_s_cp.
         self.boundary_conditions = {
             phi_s_cn: {
-                "left": (pybamm.Scalar(0), "Dirichlet"),
-                "right": (pybamm.Scalar(0), "Neumann"),
+                "negative tab": (pybamm.Scalar(0), "Dirichlet"),
+                "positive tab": (pybamm.Scalar(0), "Neumann"),
             },
             phi_s_cp: {
-                "left": (pybamm.Scalar(0), "Neumann"),
-                "right": (pos_tab_bc, "Neumann"),
+                "negative tab": (pybamm.Scalar(0), "Neumann"),
+                "positive tab": (pos_tab_bc, "Neumann"),
             },
         }
 
diff --git a/pybamm/models/submodels/thermal/x_lumped/x_lumped_1D_current_collectors.py b/pybamm/models/submodels/thermal/x_lumped/x_lumped_1D_current_collectors.py
index 6fc9a3d18..3c0fe26b6 100644
--- a/pybamm/models/submodels/thermal/x_lumped/x_lumped_1D_current_collectors.py
+++ b/pybamm/models/submodels/thermal/x_lumped/x_lumped_1D_current_collectors.py
@@ -27,19 +27,27 @@ class CurrentCollector1D(BaseModel):
 
     def set_boundary_conditions(self, variables):
         T_av = variables["X-averaged cell temperature"]
-        T_av_left = pybamm.boundary_value(T_av, "left")
-        T_av_right = pybamm.boundary_value(T_av, "right")
+        T_av_left = pybamm.boundary_value(T_av, "negative tab")
+        T_av_right = pybamm.boundary_value(T_av, "positive tab")
 
+        # Three boundary conditions here to handle the cases of both tabs at
+        # the same side (top or bottom), or one either side. For both tabs on the
+        # same side, T_av_left and T_av_right are equal, and the boundary condition
+        # "no tab" is used on the other side.
         self.boundary_conditions = {
             T_av: {
-                "left": (
+                "negative tab": (
                     self.param.h * T_av_left / self.param.lambda_k / self.param.delta,
                     "Neumann",
                 ),
-                "right": (
+                "positive tab": (
                     -self.param.h * T_av_right / self.param.lambda_k / self.param.delta,
                     "Neumann",
                 ),
+                "no tab": (
+                    pybamm.Scalar(0),
+                    "Neumann",
+                ),
             }
         }
 
diff --git a/pybamm/models/submodels/thermal/x_lumped/x_lumped_2D_current_collectors.py b/pybamm/models/submodels/thermal/x_lumped/x_lumped_2D_current_collectors.py
index 85ff431fe..1fb0dd3be 100644
--- a/pybamm/models/submodels/thermal/x_lumped/x_lumped_2D_current_collectors.py
+++ b/pybamm/models/submodels/thermal/x_lumped/x_lumped_2D_current_collectors.py
@@ -36,8 +36,8 @@ class CurrentCollector2D(BaseModel):
         # TODO: update to allow different cooling conditions at the tabs
         self.boundary_conditions = {
             T_av: {
-                "left": (pybamm.Scalar(0), "Neumann"),
-                "right": (pybamm.Scalar(0), "Neumann"),
+                "negative tab": (pybamm.Scalar(0), "Neumann"),
+                "positive tab": (pybamm.Scalar(0), "Neumann"),
             }
         }
 
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index 115dae34c..a6fb433aa 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -124,19 +124,26 @@ class ParameterValues(dict):
             model.initial_conditions[variable] = processing_function(equation)
 
         # Boundary conditions are dictionaries {"left": left bc, "right": right bc}
+        # in general, but may be imposed on the tabs (or *not* on the tab) for a
+        # small number of variables, e.g. {"negative tab": neg. tab bc,
+        # "positive tab": pos. tab bc "no tab": no tab bc}.
         new_boundary_conditions = {}
         for variable, bcs in model.boundary_conditions.items():
             processed_variable = processing_function(variable)
             new_boundary_conditions[processed_variable] = {}
-            for side in ["left", "right"]:
-                bc, typ = bcs[side]
-                pybamm.logger.debug(
-                    "{} parameters for {!r} ({} bc)".format(
-                        processing.capitalize(), variable, side
+            for side in ["left", "right", "negative tab", "positive tab", "no tab"]:
+                try:
+                    bc, typ = bcs[side]
+                    pybamm.logger.debug(
+                        "{} parameters for {!r} ({} bc)".format(
+                            processing.capitalize(), variable, side
+                        )
                     )
-                )
-                processed_bc = (processing_function(bc), typ)
-                new_boundary_conditions[processed_variable][side] = processed_bc
+                    processed_bc = (processing_function(bc), typ)
+                    new_boundary_conditions[processed_variable][side] = processed_bc
+                except KeyError:
+                    pass
+
         model.boundary_conditions = new_boundary_conditions
 
         for variable, equation in model.variables.items():
@@ -187,7 +194,7 @@ class ParameterValues(dict):
         """
         for domain in geometry:
             for prim_sec_tabs, variables in geometry[domain].items():
-                # process tab information if using 2D current collectors
+                # process tab information if using 1 or 2D current collectors
                 if prim_sec_tabs == "tabs":
                     for tab, position_size in variables.items():
                         for position_size, sym in position_size.items():
diff --git a/pybamm/spatial_methods/scikit_finite_element.py b/pybamm/spatial_methods/scikit_finite_element.py
index 4ab75dcaa..410c07c88 100644
--- a/pybamm/spatial_methods/scikit_finite_element.py
+++ b/pybamm/spatial_methods/scikit_finite_element.py
@@ -91,7 +91,7 @@ class ScikitFiniteElement(pybamm.SpatialMethod):
             The discretised symbol of the correct size
         boundary_conditions : dict
             The boundary conditions of the model
-            ({symbol.id: {"left": left bc, "right": right bc}})
+            ({symbol.id: {"negative tab": neg. tab bc, "positive tab": pos. tab bc}})
 
         Returns
         -------
@@ -104,50 +104,50 @@ class ScikitFiniteElement(pybamm.SpatialMethod):
 
         stiffness_matrix = self.stiffness_matrix(symbol, boundary_conditions)
 
-        # get boundary conditions and type, here lbc: negative tab, rbc: positive tab
-        lbc_value, lbc_type = boundary_conditions[symbol.id]["left"]
-        rbc_value, rbc_type = boundary_conditions[symbol.id]["right"]
+        # get boundary conditions and type
+        neg_bc_value, neg_bc_type = boundary_conditions[symbol.id]["negative tab"]
+        pos_bc_value, pos_bc_type = boundary_conditions[symbol.id]["positive tab"]
         # boundary load vector is adjusted to account for boundary conditions below
         boundary_load = pybamm.Vector(np.zeros(mesh.npts))
 
         # assemble boundary load if Neumann boundary conditions
-        if "Neumann" in [lbc_type, rbc_type]:
+        if "Neumann" in [neg_bc_type, pos_bc_type]:
             # make form for unit load over the boundary
             @skfem.linear_form
             def unit_bc_load_form(v, dv, w):
                 return v
 
-        if lbc_type == "Neumann":
+        if neg_bc_type == "Neumann":
             # assemble unit load over tab
-            lbc_load = skfem.asm(unit_bc_load_form, mesh.negative_tab_basis)
+            neg_bc_load = skfem.asm(unit_bc_load_form, mesh.negative_tab_basis)
             # value multiplied by weights
-            boundary_load = boundary_load + lbc_value * pybamm.Vector(lbc_load)
-        elif lbc_type == "Dirichlet":
+            boundary_load = boundary_load + neg_bc_value * pybamm.Vector(neg_bc_load)
+        elif neg_bc_type == "Dirichlet":
             # set Dirichlet value at facets corresponding to tab
-            lbc_load = np.zeros(mesh.npts)
-            lbc_load[mesh.negative_tab_dofs] = 1
-            boundary_load = boundary_load - lbc_value * pybamm.Vector(lbc_load)
+            neg_bc_load = np.zeros(mesh.npts)
+            neg_bc_load[mesh.negative_tab_dofs] = 1
+            boundary_load = boundary_load - neg_bc_value * pybamm.Vector(neg_bc_load)
         else:
             raise ValueError(
                 "boundary condition must be Dirichlet or Neumann, not '{}'".format(
-                    lbc_type
+                    neg_bc_type
                 )
             )
 
-        if rbc_type == "Neumann":
+        if pos_bc_type == "Neumann":
             # assemble unit load over tab
-            rbc_load = skfem.asm(unit_bc_load_form, mesh.positive_tab_basis)
+            pos_bc_load = skfem.asm(unit_bc_load_form, mesh.positive_tab_basis)
             # value multiplied by weights
-            boundary_load = boundary_load + rbc_value * pybamm.Vector(rbc_load)
-        elif rbc_type == "Dirichlet":
+            boundary_load = boundary_load + pos_bc_value * pybamm.Vector(pos_bc_load)
+        elif pos_bc_type == "Dirichlet":
             # set Dirichlet value at facets corresponding to tab
-            rbc_load = np.zeros(mesh.npts)
-            rbc_load[mesh.positive_tab_dofs] = 1
-            boundary_load = boundary_load - rbc_value * pybamm.Vector(rbc_load)
+            pos_bc_load = np.zeros(mesh.npts)
+            pos_bc_load[mesh.positive_tab_dofs] = 1
+            boundary_load = boundary_load - pos_bc_value * pybamm.Vector(pos_bc_load)
         else:
             raise ValueError(
                 "boundary condition must be Dirichlet or Neumann, not '{}'".format(
-                    rbc_type
+                    pos_bc_type
                 )
             )
 
@@ -172,7 +172,7 @@ class ScikitFiniteElement(pybamm.SpatialMethod):
             The symbol for which we want to calculate the laplacian matrix
         boundary_conditions : dict
             The boundary conditions of the model
-            ({symbol.id: {"left": left bc, "right": right bc}})
+            ({symbol.id: {"negative tab": neg. tab bc, "positive tab": pos. tab bc}})
 
         Returns
         -------
@@ -191,19 +191,19 @@ class ScikitFiniteElement(pybamm.SpatialMethod):
         # assemble the stifnness matrix
         stiffness = skfem.asm(stiffness_form, mesh.basis)
 
-        # get boundary conditions and type, here lbc: negative tab, rbc: positive tab
+        # get boundary conditions and type
         try:
-            _, lbc_type = boundary_conditions[symbol.id]["left"]
-            _, rbc_type = boundary_conditions[symbol.id]["right"]
+            _, neg_bc_type = boundary_conditions[symbol.id]["negative tab"]
+            _, pos_bc_type = boundary_conditions[symbol.id]["positive tab"]
         except KeyError:
             raise pybamm.ModelError(
                 "No boundary conditions provided for symbol `{}``".format(symbol)
             )
 
         # adjust matrix for Dirichlet boundary conditions
-        if lbc_type == "Dirichlet":
+        if neg_bc_type == "Dirichlet":
             self.bc_apply(stiffness, mesh.negative_tab_dofs)
-        if rbc_type == "Dirichlet":
+        if pos_bc_type == "Dirichlet":
             self.bc_apply(stiffness, mesh.positive_tab_dofs)
 
         return pybamm.Matrix(stiffness)
@@ -329,19 +329,19 @@ class ScikitFiniteElement(pybamm.SpatialMethod):
 
     def boundary_value_or_flux(self, symbol, discretised_child):
         """
-        Returns the average value of the symbol over the negative tab ("left")
-        or the positive tab ("right") in the Finite Element Method.
+        Returns the average value of the symbol over the negative tab ("negative tab")
+        or the positive tab ("positive tab") in the Finite Element Method.
 
         Overwrites the default :meth:`pybamm.SpatialMethod.boundary_value`
         """
 
-        # Return average value on the negative tab for "left" and positive tab
-        # for "right"
+        # Return average value on the negative tab for "negative tab" and positive tab
+        # for "positive tab"
         if isinstance(symbol, pybamm.BoundaryValue):
             # get integration_vector
-            if symbol.side == "left":
+            if symbol.side == "negative tab":
                 region = "negative tab"
-            elif symbol.side == "right":
+            elif symbol.side == "positive tab":
                 region = "positive tab"
             domain = symbol.children[0].domain[0]
             integration_vector = self.boundary_integral_vector(domain, region=region)
@@ -372,7 +372,7 @@ class ScikitFiniteElement(pybamm.SpatialMethod):
             calculating the mass matrix.
         boundary_conditions : dict
             The boundary conditions of the model
-            ({symbol.id: {"left": left bc, "right": right bc}})
+            ({symbol.id: {"negative tab": neg. tab bc, "positive tab": pos. tab bc}})
 
         Returns
         -------
@@ -393,7 +393,7 @@ class ScikitFiniteElement(pybamm.SpatialMethod):
             calculating the mass matrix.
         boundary_conditions : dict
             The boundary conditions of the model
-            ({symbol.id: {"left": left bc, "right": right bc}})
+            ({symbol.id: {"negative tab": neg. tab bc, "positive tab": pos. tab bc}})
 
         Returns
         -------
@@ -414,7 +414,7 @@ class ScikitFiniteElement(pybamm.SpatialMethod):
             calculating the mass matrix.
         boundary_conditions : dict
             The boundary conditions of the model
-            ({symbol.id: {"left": left bc, "right": right bc}})
+            ({symbol.id: {"negative tab": neg. tab bc, "positive tab": pos. tab bc}})
         region: str, optional
             The domain over which to assemble the mass matrix form. Can be "interior"
             (default) or "boundary".
@@ -439,15 +439,15 @@ class ScikitFiniteElement(pybamm.SpatialMethod):
         if region == "boundary":
             mass = skfem.asm(mass_form, mesh.facet_basis)
 
-        # get boundary conditions and type, here lbc: negative tab, rbc: positive tab
+        # get boundary conditions and type
         if symbol.id in boundary_conditions:
-            _, lbc_type = boundary_conditions[symbol.id]["left"]
-            _, rbc_type = boundary_conditions[symbol.id]["right"]
+            _, neg_bc_type = boundary_conditions[symbol.id]["negative tab"]
+            _, pos_bc_type = boundary_conditions[symbol.id]["positive tab"]
 
-            if lbc_type == "Dirichlet":
+            if neg_bc_type == "Dirichlet":
                 # set source terms to zero on boundary by zeroing out mass matrix
                 self.bc_apply(mass, mesh.negative_tab_dofs, zero=True)
-            if rbc_type == "Dirichlet":
+            if pos_bc_type == "Dirichlet":
                 # set source terms to zero on boundary by zeroing out mass matrix
                 self.bc_apply(mass, mesh.positive_tab_dofs, zero=True)
 
diff --git a/results/2plus1D/dfn_2plus1D.py b/results/2plus1D/dfn_2plus1D.py
index 1f8efc831..af05b8770 100644
--- a/results/2plus1D/dfn_2plus1D.py
+++ b/results/2plus1D/dfn_2plus1D.py
@@ -10,7 +10,7 @@ pybamm.set_logging_level("INFO")
 options = {
     "current collector": "potential pair",
     "dimensionality": 2,
-    "thermal": "full",
+    "thermal": "x-lumped",
 }
 model = pybamm.lithium_ion.DFN(options)
 


Test patch:
diff --git a/tests/shared.py b/tests/shared.py
index d2d4a94b1..1e55ce4d4 100644
--- a/tests/shared.py
+++ b/tests/shared.py
@@ -47,7 +47,7 @@ def get_mesh_for_testing(
             "Electrode height [m]": 0.5,
             "Negative tab width [m]": 0.1,
             "Negative tab centre y-coordinate [m]": 0.1,
-            "Negative tab centre z-coordinate [m]": 0.5,
+            "Negative tab centre z-coordinate [m]": 0.0,
             "Positive tab width [m]": 0.1,
             "Positive tab centre y-coordinate [m]": 0.3,
             "Positive tab centre z-coordinate [m]": 0.5,
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 5f3bb7de0..820cde672 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -821,6 +821,21 @@ class TestDiscretise(unittest.TestCase):
         }
         disc.process_model(model)
 
+    def test_check_tab_bcs_error(self):
+        a = pybamm.Variable("a", domain=["current collector"])
+        b = pybamm.Variable("b", domain=["negative electrode"])
+        bcs = {"negative tab": (0, "Dirichlet"), "positive tab": (0, "Neumann")}
+
+        disc = get_discretisation_for_testing()
+
+        # for 0D bcs keys should be unchanged
+        new_bcs = disc.check_tab_conditions(a, bcs)
+        self.assertListEqual(list(bcs.keys()), list(new_bcs.keys()))
+
+        # error if domain not "current collector"
+        with self.assertRaisesRegex(pybamm.ModelError, "Boundary conditions"):
+            disc.check_tab_conditions(b, bcs)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index 5308eebe3..156812525 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -227,6 +227,12 @@ class TestUnaryOperators(unittest.TestCase):
         self.assertEqual(boundary_a_tert.domain, ["current collector"])
         self.assertEqual(boundary_a_tert.auxiliary_domains, {"secondary": ["bla"]})
 
+        # error if boundary value on tabs and domain is not "current collector"
+        var = pybamm.Variable("var", domain=["negative electrode"])
+        with self.assertRaisesRegex(pybamm.ModelError, "Can only take boundary"):
+            pybamm.boundary_value(var, "negative tab")
+            pybamm.boundary_value(var, "positive tab")
+
     def test_average(self):
         a = pybamm.Scalar(1)
         average_a = pybamm.x_average(a)
diff --git a/tests/unit/test_meshes/test_meshes.py b/tests/unit/test_meshes/test_meshes.py
index c0d76faae..26366ef4b 100644
--- a/tests/unit/test_meshes/test_meshes.py
+++ b/tests/unit/test_meshes/test_meshes.py
@@ -378,6 +378,80 @@ class TestMesh(unittest.TestCase):
         with self.assertRaises(pybamm.GeometryError):
             pybamm.Mesh(geometry, submesh_types, var_pts)
 
+    def test_1plus1D_tabs_left_right(self):
+        param = pybamm.ParameterValues(
+            base_parameters={
+                "Electrode width [m]": 0.4,
+                "Electrode height [m]": 0.5,
+                "Negative tab centre z-coordinate [m]": 0.0,
+                "Positive tab centre z-coordinate [m]": 0.5,
+                "Negative electrode thickness [m]": 0.3,
+                "Separator thickness [m]": 0.3,
+                "Positive electrode thickness [m]": 0.3,
+            }
+        )
+
+        geometry = pybamm.Geometryxp1DMacro(cc_dimension=1)
+        param.process_geometry(geometry)
+
+        var = pybamm.standard_spatial_vars
+        var_pts = {var.x_n: 10, var.x_s: 7, var.x_p: 12, var.z: 24}
+
+        submesh_types = {
+            "negative electrode": pybamm.Uniform1DSubMesh,
+            "separator": pybamm.Uniform1DSubMesh,
+            "positive electrode": pybamm.Uniform1DSubMesh,
+            "current collector": pybamm.Uniform1DSubMesh,
+        }
+
+        mesh_type = pybamm.Mesh
+
+        # create mesh
+        mesh = mesh_type(geometry, submesh_types, var_pts)
+
+        # negative tab should be "left"
+        self.assertEqual(mesh["current collector"][0].tabs["negative tab"], "left")
+
+        # positive tab should be "right"
+        self.assertEqual(mesh["current collector"][0].tabs["positive tab"], "right")
+
+    def test_1plus1D_tabs_right_left(self):
+        param = pybamm.ParameterValues(
+            base_parameters={
+                "Electrode width [m]": 0.4,
+                "Electrode height [m]": 0.5,
+                "Negative tab centre z-coordinate [m]": 0.5,
+                "Positive tab centre z-coordinate [m]": 0.0,
+                "Negative electrode thickness [m]": 0.3,
+                "Separator thickness [m]": 0.3,
+                "Positive electrode thickness [m]": 0.3,
+            }
+        )
+
+        geometry = pybamm.Geometryxp1DMacro(cc_dimension=1)
+        param.process_geometry(geometry)
+
+        var = pybamm.standard_spatial_vars
+        var_pts = {var.x_n: 10, var.x_s: 7, var.x_p: 12, var.z: 24}
+
+        submesh_types = {
+            "negative electrode": pybamm.Uniform1DSubMesh,
+            "separator": pybamm.Uniform1DSubMesh,
+            "positive electrode": pybamm.Uniform1DSubMesh,
+            "current collector": pybamm.Uniform1DSubMesh,
+        }
+
+        mesh_type = pybamm.Mesh
+
+        # create mesh
+        mesh = mesh_type(geometry, submesh_types, var_pts)
+
+        # negative tab should be "right"
+        self.assertEqual(mesh["current collector"][0].tabs["negative tab"], "right")
+
+        # positive tab should be "left"
+        self.assertEqual(mesh["current collector"][0].tabs["positive tab"], "left")
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_spatial_methods/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume.py
index 441617acb..fbe5ba1a7 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume.py
@@ -173,6 +173,31 @@ class TestFiniteVolume(unittest.TestCase):
         extrap_right_disc = disc.process_symbol(extrap_right)
         self.assertEqual(extrap_right_disc.domain, [])
 
+        # test extrapolate to "negative tab" gives same as "left" and
+        # "positive tab" gives same "right" (see get_mesh_for_testing)
+        var = pybamm.Variable("var", domain="current collector")
+        disc.set_variable_slices([var])
+        submesh = mesh["current collector"]
+        constant_y = np.ones_like(submesh[0].nodes[:, np.newaxis])
+
+        extrap_neg = pybamm.BoundaryValue(var, "negative tab")
+        extrap_neg_disc = disc.process_symbol(extrap_neg)
+        extrap_left = pybamm.BoundaryValue(var, "left")
+        extrap_left_disc = disc.process_symbol(extrap_left)
+        np.testing.assert_array_equal(
+            extrap_neg_disc.evaluate(None, constant_y),
+            extrap_left_disc.evaluate(None, constant_y),
+        )
+
+        extrap_pos = pybamm.BoundaryValue(var, "positive tab")
+        extrap_pos_disc = disc.process_symbol(extrap_pos)
+        extrap_right = pybamm.BoundaryValue(var, "right")
+        extrap_right_disc = disc.process_symbol(extrap_right)
+        np.testing.assert_array_equal(
+            extrap_pos_disc.evaluate(None, constant_y),
+            extrap_right_disc.evaluate(None, constant_y),
+        )
+
     def test_discretise_diffusivity_times_spatial_operator(self):
         # Set up
         whole_cell = ["negative electrode", "separator", "positive electrode"]
@@ -1332,6 +1357,66 @@ class TestFiniteVolume(unittest.TestCase):
         self.assertEqual(c_s_n_surf_disc.domain, ["negative electrode"])
         self.assertEqual(c_s_p_surf_disc.domain, ["positive electrode"])
 
+    def test_grad_div_with_bcs_on_tab(self):
+        # 2d macroscale
+        mesh = get_1p1d_mesh_for_testing()
+        spatial_methods = {
+            "macroscale": pybamm.FiniteVolume,
+            "negative particle": pybamm.FiniteVolume,
+            "positive particle": pybamm.FiniteVolume,
+            "current collector": pybamm.FiniteVolume,
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        y_test = np.ones(mesh["current collector"][0].npts)
+
+        # var
+        var = pybamm.Variable("var", domain="current collector")
+        disc.set_variable_slices([var])
+        # grad
+        grad_eqn = pybamm.grad(var)
+        # div
+        N = pybamm.grad(var)
+        div_eqn = pybamm.div(N)
+
+        # bcs (on each tab)
+        boundary_conditions = {
+            var.id: {
+                "negative tab": (pybamm.Scalar(1), "Dirichlet"),
+                "positive tab": (pybamm.Scalar(0), "Neumann"),
+            }
+        }
+        disc.bcs = boundary_conditions
+        grad_eqn_disc = disc.process_symbol(grad_eqn)
+        grad_eqn_disc.evaluate(None, y_test)
+        div_eqn_disc = disc.process_symbol(div_eqn)
+        div_eqn_disc.evaluate(None, y_test)
+
+        # bcs (one pos, one not tab)
+        boundary_conditions = {
+            var.id: {
+                "no tab": (pybamm.Scalar(1), "Dirichlet"),
+                "positive tab": (pybamm.Scalar(0), "Dirichlet"),
+            }
+        }
+        disc.bcs = boundary_conditions
+        grad_eqn_disc = disc.process_symbol(grad_eqn)
+        grad_eqn_disc.evaluate(None, y_test)
+        div_eqn_disc = disc.process_symbol(div_eqn)
+        div_eqn_disc.evaluate(None, y_test)
+
+        # bcs (one neg, one not tab)
+        boundary_conditions = {
+            var.id: {
+                "negative tab": (pybamm.Scalar(1), "Neumann"),
+                "no tab": (pybamm.Scalar(0), "Neumann"),
+            }
+        }
+        disc.bcs = boundary_conditions
+        grad_eqn_disc = disc.process_symbol(grad_eqn)
+        grad_eqn_disc.evaluate(None, y_test)
+        div_eqn_disc = disc.process_symbol(div_eqn)
+        div_eqn_disc.evaluate(None, y_test)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_spatial_methods/test_scikit_finite_element.py b/tests/unit/test_spatial_methods/test_scikit_finite_element.py
index 69af77c83..f2765a88e 100644
--- a/tests/unit/test_spatial_methods/test_scikit_finite_element.py
+++ b/tests/unit/test_spatial_methods/test_scikit_finite_element.py
@@ -37,8 +37,8 @@ class TestScikitFiniteElement(unittest.TestCase):
         unit_source = pybamm.Broadcast(1, "current collector")
         disc.bcs = {
             var.id: {
-                "left": (pybamm.Scalar(0), "Neumann"),
-                "right": (pybamm.Scalar(0), "Neumann"),
+                "negative tab": (pybamm.Scalar(0), "Neumann"),
+                "positive tab": (pybamm.Scalar(0), "Neumann"),
             }
         }
 
@@ -60,8 +60,8 @@ class TestScikitFiniteElement(unittest.TestCase):
             # Dirichlet
             disc.bcs = {
                 var.id: {
-                    "left": (pybamm.Scalar(0), "Dirichlet"),
-                    "right": (pybamm.Scalar(1), "Dirichlet"),
+                    "negative tab": (pybamm.Scalar(0), "Dirichlet"),
+                    "positive tab": (pybamm.Scalar(1), "Dirichlet"),
                 }
             }
             eqn_disc = disc.process_symbol(eqn)
@@ -69,8 +69,8 @@ class TestScikitFiniteElement(unittest.TestCase):
             # Neumann
             disc.bcs = {
                 var.id: {
-                    "left": (pybamm.Scalar(0), "Neumann"),
-                    "right": (pybamm.Scalar(1), "Neumann"),
+                    "negative tab": (pybamm.Scalar(0), "Neumann"),
+                    "positive tab": (pybamm.Scalar(1), "Neumann"),
                 }
             }
             eqn_disc = disc.process_symbol(eqn)
@@ -78,8 +78,8 @@ class TestScikitFiniteElement(unittest.TestCase):
             # One of each
             disc.bcs = {
                 var.id: {
-                    "left": (pybamm.Scalar(0), "Neumann"),
-                    "right": (pybamm.Scalar(1), "Dirichlet"),
+                    "negative tab": (pybamm.Scalar(0), "Neumann"),
+                    "positive tab": (pybamm.Scalar(1), "Dirichlet"),
                 }
             }
             eqn_disc = disc.process_symbol(eqn)
@@ -87,8 +87,8 @@ class TestScikitFiniteElement(unittest.TestCase):
             # One of each
             disc.bcs = {
                 var.id: {
-                    "left": (pybamm.Scalar(0), "Dirichlet"),
-                    "right": (pybamm.Scalar(1), "Neumann"),
+                    "negative tab": (pybamm.Scalar(0), "Dirichlet"),
+                    "positive tab": (pybamm.Scalar(1), "Neumann"),
                 }
             }
             eqn_disc = disc.process_symbol(eqn)
@@ -98,16 +98,16 @@ class TestScikitFiniteElement(unittest.TestCase):
         eqn = pybamm.laplacian(var) - pybamm.source(unit_source, var)
         disc.bcs = {
             var.id: {
-                "left": (pybamm.Scalar(0), "Dirichlet"),
-                "right": (pybamm.Scalar(1), "Other BC"),
+                "negative tab": (pybamm.Scalar(0), "Dirichlet"),
+                "positive tab": (pybamm.Scalar(1), "Other BC"),
             }
         }
         with self.assertRaises(ValueError):
             eqn_disc = disc.process_symbol(eqn)
         disc.bcs = {
             var.id: {
-                "left": (pybamm.Scalar(0), "Other BC"),
-                "right": (pybamm.Scalar(1), "Neumann"),
+                "negative tab": (pybamm.Scalar(0), "Other BC"),
+                "positive tab": (pybamm.Scalar(1), "Neumann"),
             }
         }
         with self.assertRaises(ValueError):
@@ -157,12 +157,12 @@ class TestScikitFiniteElement(unittest.TestCase):
         # laplace of u = sin(pi*z)
         var = pybamm.Variable("var", domain="current collector")
         eqn_zz = pybamm.laplacian(var)
-        # set boundary conditions ("left" = bottom of unit square, "right" = top
-        # of unit square, elsewhere normal derivative is zero)
+        # set boundary conditions ("negative tab" = bottom of unit square,
+        # "positive tab" = top of unit square, elsewhere normal derivative is zero)
         disc.bcs = {
             var.id: {
-                "left": (pybamm.Scalar(0), "Dirichlet"),
-                "right": (pybamm.Scalar(0), "Dirichlet"),
+                "negative tab": (pybamm.Scalar(0), "Dirichlet"),
+                "positive tab": (pybamm.Scalar(0), "Dirichlet"),
             }
         }
         disc.set_variable_slices([var])
@@ -179,12 +179,12 @@ class TestScikitFiniteElement(unittest.TestCase):
         # laplace of u = cos(pi*y)*sin(pi*z)
         var = pybamm.Variable("var", domain="current collector")
         laplace_eqn = pybamm.laplacian(var)
-        # set boundary conditions ("left" = bottom of unit square, "right" = top
-        # of unit square, elsewhere normal derivative is zero)
+        # set boundary conditions ("negative tab" = bottom of unit square,
+        # "positive tab" = top of unit square, elsewhere normal derivative is zero)
         disc.bcs = {
             var.id: {
-                "left": (pybamm.Scalar(0), "Dirichlet"),
-                "right": (pybamm.Scalar(0), "Dirichlet"),
+                "negative tab": (pybamm.Scalar(0), "Dirichlet"),
+                "positive tab": (pybamm.Scalar(0), "Dirichlet"),
             }
         }
         disc.set_variable_slices([var])
@@ -242,7 +242,7 @@ class TestScikitFiniteElement(unittest.TestCase):
         self.assertEqual(vec_disc.shape[0], mesh["current collector"][0].npts)
         self.assertEqual(vec_disc.shape[1], 1)
 
-    def test_left_right(self):
+    def test_neg_pos(self):
         mesh = get_2p1d_mesh_for_testing()
         spatial_methods = {
             "macroscale": pybamm.FiniteVolume,
@@ -252,15 +252,18 @@ class TestScikitFiniteElement(unittest.TestCase):
         var = pybamm.Variable("var", domain="current collector")
         disc.set_variable_slices([var])
 
-        extrap_left = pybamm.BoundaryValue(var, "left")
-        extrap_right = pybamm.BoundaryValue(var, "right")
-        extrap_left_disc = disc.process_symbol(extrap_left)
-        extrap_right_disc = disc.process_symbol(extrap_right)
-
+        extrap_neg = pybamm.BoundaryValue(var, "negative tab")
+        extrap_pos = pybamm.BoundaryValue(var, "positive tab")
+        extrap_neg_disc = disc.process_symbol(extrap_neg)
+        extrap_pos_disc = disc.process_symbol(extrap_pos)
         # check constant returns constant at tab
-        constant_y = np.ones(mesh["current collector"][0].npts)
-        self.assertEqual(extrap_left_disc.evaluate(None, constant_y), 1)
-        self.assertEqual(extrap_right_disc.evaluate(None, constant_y), 1)
+        constant_y = np.ones(mesh["current collector"][0].npts)[:, np.newaxis]
+        np.testing.assert_array_almost_equal(
+            extrap_neg_disc.evaluate(None, constant_y), 1
+        )
+        np.testing.assert_array_almost_equal(
+            extrap_pos_disc.evaluate(None, constant_y), 1
+        )
 
     def test_boundary_integral(self):
         mesh = get_2p1d_mesh_for_testing()
@@ -314,10 +317,10 @@ class TestScikitFiniteElement(unittest.TestCase):
             c: pybamm.Integral(u, [y, z]) + 0 * c,
         }
         model.initial_conditions = {u: pybamm.Scalar(0), c: pybamm.Scalar(0)}
-        # set boundary conditions ("left" = bottom of unit square, "right" = top
-        # of unit square, elsewhere normal derivative is zero)
+        # set boundary conditions ("negative tab" = bottom of unit square,
+        # "positive tab" = top of unit square, elsewhere normal derivative is zero)
         model.boundary_conditions = {
-            u: {"left": (0, "Neumann"), "right": (1, "Neumann")}
+            u: {"negative tab": (0, "Neumann"), "positive tab": (1, "Neumann")}
         }
         model.variables = {"c": c, "u": u}
         # create discretisation