Repository: pybamm-team/PyBaMM
Pull request #5112 (2025-07-12 08:46:21)

3D Thermal Model
# Description

Added a `Basic3DThermalSPM` model with two way coupling.

Fixes # (issue)

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CHANGELOG.md) to document the change (include PR #)

# Important checks:

Please confirm the following before marking the PR as ready for review:
- [x] No style issues: `nox -s pre-commit`
- [x] All tests pass: `nox -s tests`
- [x] The documentation builds: `nox -s doctests`
- [x] Code is commented for hard-to-understand areas
- [x] Tests added that prove fix is effective or that feature works


Rishab87: @rtimms can you review this?

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6f6ef3199..6fc3206f0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,7 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
 ## Features
+- Added a `Basic3DThermalSPM` with two way coupling. ([#5112](https://github.com/pybamm-team/PyBaMM/pull/5112))
 - Enables the passing of `inputs` throughout `set_initial_soc`. ([#5122](https://github.com/pybamm-team/PyBaMM/pull/5122))
 - Adds `on_failure` option to `BaseSolver` with options for `"warn"`, `"ignore"`, and `"raise"` to change behaviour on solver failure. Defaults to "raise" to retain historic functionality. ([#5105](https://github.com/pybamm-team/PyBaMM/pull/5105))
 - Creates a boundary mesh size object that returns the distance from the center of the leftmost/rightmost control volume to the boundary of the domain ([#5108](https://github.com/pybamm-team/PyBaMM/pull/5108))
diff --git a/docs/source/api/models/lithium_ion/spm.rst b/docs/source/api/models/lithium_ion/spm.rst
index fac565455..35df566ee 100644
--- a/docs/source/api/models/lithium_ion/spm.rst
+++ b/docs/source/api/models/lithium_ion/spm.rst
@@ -7,4 +7,6 @@ Single Particle Model (SPM)
 .. autoclass:: pybamm.lithium_ion.BasicSPM
     :members:
 
+.. autoclass:: pybamm.lithium_ion.Basic3DThermalSPM
+
 .. footbibliography::
diff --git a/examples/scripts/3d_examples/3d_basic_spm_vs_lumped.py b/examples/scripts/3d_examples/3d_basic_spm_vs_lumped.py
new file mode 100644
index 000000000..7c1d54414
--- /dev/null
+++ b/examples/scripts/3d_examples/3d_basic_spm_vs_lumped.py
@@ -0,0 +1,109 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+import pybamm
+
+# Load models
+models = {
+    "Lumped": pybamm.lithium_ion.SPM(options={"thermal": "lumped"}),
+    "3D": pybamm.lithium_ion.Basic3DThermalSPM(
+        options={"cell geometry": "pouch", "dimensionality": 3}
+    ),
+}
+
+# Load base parameters
+parameter_values = pybamm.ParameterValues("Marquis2019")
+
+# Define heat transfer coefficients to test
+h_values = [0.1, 1, 10, 100]  # W/m2.K
+
+# Set up experiment
+experiment = pybamm.Experiment(
+    [
+        ("Discharge at 3C until 2.8V", "Rest for 10 minutes"),
+    ]
+)
+
+# Set spatial points
+var_pts = {
+    "x_n": 20,
+    "x_s": 20,
+    "x_p": 20,
+    "r_n": 30,
+    "r_p": 30,
+    "x": None,
+    "y": None,
+    "z": None,
+}
+
+# Dictionary to store solutions for each h value
+all_solutions = {}
+
+# Loop over different h values
+for h in h_values:
+    print(f"\nSolving models for h = {h} W/m2.K...")
+
+    # Update parameter values for this h
+    h_params = parameter_values.copy()
+    h_params.update(
+        {
+            # Used in lumped model
+            "Total heat transfer coefficient [W.m-2.K-1]": h,
+            # Used in 3D model
+            "Left face heat transfer coefficient [W.m-2.K-1]": h,
+            "Right face heat transfer coefficient [W.m-2.K-1]": h,
+            "Front face heat transfer coefficient [W.m-2.K-1]": h,
+            "Back face heat transfer coefficient [W.m-2.K-1]": h,
+            "Bottom face heat transfer coefficient [W.m-2.K-1]": h,
+            "Top face heat transfer coefficient [W.m-2.K-1]": h,
+        },
+        check_already_exists=False,
+    )
+
+    solutions = {}
+    for model_name, model in models.items():
+        print(f"Solving {model_name} model...")
+        sim = pybamm.Simulation(
+            model, parameter_values=h_params, var_pts=var_pts, experiment=experiment
+        )
+        solutions[model_name] = sim.solve()
+        print(f"{model_name} model solved.")
+
+    all_solutions[h] = solutions
+
+# Plot results for all h values on same axes
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))
+# Define color map for h values
+colors = plt.cm.viridis(np.linspace(0, 1, len(h_values)))
+styles = {
+    "3D": {"linestyle": "-", "label": "3D Thermal Model"},
+    "Lumped": {"linestyle": "--", "label": "Lumped Thermal Model"},
+}
+
+for i, h in enumerate(h_values):
+    solutions = all_solutions[h]
+    for name, sol in solutions.items():
+        style = styles[name].copy()
+        style["color"] = colors[i]
+        style["label"] = f"{style['label']} (h={h})"
+
+        ax1.plot(sol["Time [s]"].data, sol["Voltage [V]"].data, **style)
+        ax2.plot(
+            sol["Time [s]"].data,
+            sol["Volume-averaged cell temperature [K]"].data,
+            **style,
+        )
+
+ax1.set_xlabel("Time [s]")
+ax1.set_ylabel("Voltage [V]")
+ax1.set_title("Voltage Comparison")
+ax1.grid(True)
+
+ax2.set_xlabel("Time [s]")
+ax2.set_ylabel("Temperature [K]")
+ax2.set_title("Temperature Comparison")
+ax2.grid(True)
+ax2.legend(bbox_to_anchor=(1.05, 1), loc="upper left")
+fig.suptitle("Model Comparison: 3D Thermal vs. Lumped Thermal", fontsize=16)
+plt.tight_layout()
+plt.show()
diff --git a/examples/scripts/3d_examples/3d_equation_prismatic.py b/examples/scripts/3d_examples/3d_equation_prismatic.py
index ade089f13..ae60db1af 100644
--- a/examples/scripts/3d_examples/3d_equation_prismatic.py
+++ b/examples/scripts/3d_examples/3d_equation_prismatic.py
@@ -48,7 +48,7 @@ geometry = {
 }
 
 submesh_types = {
-    "current collector": pybamm.ScikitFemGenerator3D(geom_type="box", h=0.15)
+    "current collector": pybamm.ScikitFemGenerator3D(geom_type="pouch", h=0.15)
 }
 var_pts = {x: None, y: None, z: None}
 mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
diff --git a/src/pybamm/__init__.py b/src/pybamm/__init__.py
index b33712a96..89b2d370d 100644
--- a/src/pybamm/__init__.py
+++ b/src/pybamm/__init__.py
@@ -164,6 +164,7 @@ from .solvers.solution import Solution, EmptySolution, make_cycle_solution
 from .solvers.processed_variable_time_integral import ProcessedVariableTimeIntegral
 from .solvers.processed_variable import ProcessedVariable, process_variable
 from .solvers.processed_variable_computed import ProcessedVariableComputed
+from .solvers.processed_variable import ProcessedVariableUnstructured
 from .solvers.summary_variable import SummaryVariables
 from .solvers.base_solver import BaseSolver
 from .solvers.dummy_solver import DummySolver
diff --git a/src/pybamm/expression_tree/binary_operators.py b/src/pybamm/expression_tree/binary_operators.py
index e9c20b103..5b271ce4f 100644
--- a/src/pybamm/expression_tree/binary_operators.py
+++ b/src/pybamm/expression_tree/binary_operators.py
@@ -1542,16 +1542,18 @@ def source(
         corresponding to a source term in the bulk.
     """
     # Broadcast if left is number
+    right_domain = right.domain
     if isinstance(left, numbers.Number):
-        left = pybamm.PrimaryBroadcast(left, "current collector")
+        left = pybamm.PrimaryBroadcast(left, right_domain)
 
     # force type cast for mypy
     left = cast(pybamm.Symbol, left)
 
-    if left.domain != ["current collector"] or right.domain != ["current collector"]:
+    allowed_domains = [["cell"], ["current collector"]]
+    if left.domain not in allowed_domains or right_domain not in allowed_domains:
         raise pybamm.DomainError(
-            "'source' only implemented in the 'current collector' domain, "
-            f"but symbols have domains {left.domain} and {right.domain}"
+            "'source' only implemented in the 'cell' or 'current collector' domains, "
+            f"but symbols have domains {left.domain} and {right_domain}"
         )
     if boundary:
         return pybamm.BoundaryMass(right) @ left
diff --git a/src/pybamm/geometry/battery_geometry.py b/src/pybamm/geometry/battery_geometry.py
index 3e1986507..01de12edb 100644
--- a/src/pybamm/geometry/battery_geometry.py
+++ b/src/pybamm/geometry/battery_geometry.py
@@ -150,6 +150,14 @@ def battery_geometry(
                     },
                 },
             }
+        elif current_collector_dimension == 3:
+            geometry["current collector"] = {"z": {"position": 1}}
+            geometry["cell"] = {
+                "x": {"min": 0, "max": geo.L_x},
+                "y": {"min": 0, "max": geo.L_y},
+                "z": {"min": 0, "max": geo.L_z},
+            }
+
     elif form_factor == "cylindrical":
         if current_collector_dimension == 0:
             geometry["current collector"] = {"r_macro": {"position": 1}}
@@ -157,6 +165,12 @@ def battery_geometry(
             geometry["current collector"] = {
                 "r_macro": {"min": geo.r_inner, "max": 1},
             }
+        elif current_collector_dimension == 3:
+            geometry["current collector"] = {"z": {"position": 1}}
+            geometry["cell"] = {
+                "r_macro": {"min": geo.r_inner, "max": geo.r_outer},
+                "z": {"min": 0, "max": geo.L_z},
+            }
         else:
             raise pybamm.GeometryError(
                 f"Invalid current collector dimension '{current_collector_dimension}' (should be 0 or 1 for "
@@ -164,7 +178,7 @@ def battery_geometry(
             )
     else:
         raise pybamm.GeometryError(
-            f"Invalid form factor '{form_factor}' (should be 'pouch' or 'cylindrical'"
+            f"Invalid form factor '{form_factor}' (should be 'pouch' or 'cylindrical')"
         )
 
     return pybamm.Geometry(geometry)
diff --git a/src/pybamm/meshes/meshes.py b/src/pybamm/meshes/meshes.py
index e90afd56c..15152b29a 100644
--- a/src/pybamm/meshes/meshes.py
+++ b/src/pybamm/meshes/meshes.py
@@ -44,7 +44,11 @@ class Mesh(dict):
         for domain, generator in submesh_types.items():
             if isinstance(generator, pybamm.ScikitFemGenerator3D):
                 for var in geometry[domain]:
-                    provided_pts = var_name_pts.get(var.name)
+                    if isinstance(var, str):
+                        var_name = var
+                    else:
+                        var_name = var.name
+                    provided_pts = var_name_pts.get(var_name)
                     if provided_pts is not None:
                         pybamm.logger.warning(
                             f"For the 3D FEM submesh on domain '{domain}', the value "
diff --git a/src/pybamm/meshes/scikit_fem_submeshes_3d.py b/src/pybamm/meshes/scikit_fem_submeshes_3d.py
index 22f865dae..5e41c9efe 100644
--- a/src/pybamm/meshes/scikit_fem_submeshes_3d.py
+++ b/src/pybamm/meshes/scikit_fem_submeshes_3d.py
@@ -11,7 +11,7 @@ class ScikitFemGenerator3D(pybamm.MeshGenerator):
     Parameters
     ----------
     geom_type : str
-        The type of geometry to generate. Must be one of "box" for a rectangular
+        The type of geometry to generate. Must be one of "pouch" for a rectangular
         prism, or "cylinder" for a cylindrical annulus.
     gen_params : dict
         A dictionary of geometry-specific parameters. for eg:
@@ -23,7 +23,7 @@ class ScikitFemGenerator3D(pybamm.MeshGenerator):
     """
 
     def __init__(self, geom_type, **gen_params):
-        supported_geometries = ["box", "cylinder"]
+        supported_geometries = ["pouch", "cylinder"]
         if geom_type not in supported_geometries:
             raise pybamm.GeometryError(
                 f"geom_type must be one of {supported_geometries}, not '{geom_type}'"
@@ -32,9 +32,9 @@ class ScikitFemGenerator3D(pybamm.MeshGenerator):
         self.geom_type = geom_type
         self.gen_params = gen_params
 
-    def _make_box_mesh(self, x_lim, y_lim, z_lim, h):
+    def _make_pouch_mesh(self, x_lim, y_lim, z_lim, h):
         """
-        Create a structured box mesh.
+        Create a structured pouch mesh.
 
         Parameters
         ----------
@@ -46,7 +46,7 @@ class ScikitFemGenerator3D(pybamm.MeshGenerator):
         Returns
         -------
         skfem.MeshTet
-            Box mesh with proper boundary tags
+            Pouch mesh with proper boundary tags
         """
         skfem = import_optional_dependency("skfem")
         x_min, x_max = x_lim
@@ -195,23 +195,29 @@ class ScikitFemGenerator3D(pybamm.MeshGenerator):
         ScikitFemSubMesh3D
             Generated 3D submesh
         """
-        h = self.gen_params.get("h", 0.3)
-
-        if self.geom_type == "box":
-            x_key = next(k for k in lims if k.name == "x")
-            y_key = next(k for k in lims if k.name == "y")
-            z_key = next(k for k in lims if k.name == "z")
-            x_lim = tuple(lims[x_key].values())
-            y_lim = tuple(lims[y_key].values())
-            z_lim = tuple(lims[z_key].values())
+        h = float(self.gen_params.get("h", 0.3))
+
+        lims_dict = {}
+        for var, lim in lims.items():
+            if isinstance(var, pybamm.SpatialVariable):
+                lims_dict[var.name] = (lim["min"], lim["max"])
+            else:
+                lims_dict[var] = (float(lim["min"]), float(lim["max"]))
+
+        if self.geom_type == "pouch":
+            x_lim = lims_dict["x"]
+            y_lim = lims_dict["y"]
+            z_lim = lims_dict["z"]
             coord_sys = "cartesian"
-            mesh = self._make_box_mesh(x_lim, y_lim, z_lim, h)
+            mesh = self._make_pouch_mesh(x_lim, y_lim, z_lim, h)
 
         elif self.geom_type == "cylinder":
-            r_key = next(k for k in lims if k.name == "r")
-            z_key = next(k for k in lims if k.name == "z")
-            r_lim = tuple(lims[r_key].values())
-            z_lim = tuple(lims[z_key].values())
+            if "r" in lims_dict:
+                r_lim = lims_dict["r"]
+            else:
+                r_lim = lims_dict["r_macro"]
+
+            z_lim = lims_dict["z"]
             coord_sys = "cylindrical polar"
             mesh = self._make_cylindrical_mesh(r_lim, z_lim, h)
         else:
diff --git a/src/pybamm/models/full_battery_models/base_battery_model.py b/src/pybamm/models/full_battery_models/base_battery_model.py
index b5c3f48b7..c8283db63 100644
--- a/src/pybamm/models/full_battery_models/base_battery_model.py
+++ b/src/pybamm/models/full_battery_models/base_battery_model.py
@@ -294,7 +294,7 @@ class BatteryModelOptions(pybamm.FuzzyDict):
         self.possible_options = {
             "calculate discharge energy": ["false", "true"],
             "calculate heat source for isothermal models": ["false", "true"],
-            "cell geometry": ["arbitrary", "pouch"],
+            "cell geometry": ["arbitrary", "pouch", "cylindrical"],
             "contact resistance": ["false", "true"],
             "convection": ["none", "uniform transverse", "full transverse"],
             "current collector": [
@@ -303,7 +303,7 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 "potential pair quite conductive",
             ],
             "diffusivity": ["single", "current sigmoid"],
-            "dimensionality": [0, 1, 2],
+            "dimensionality": [0, 1, 2, 3],
             "electrolyte conductivity": [
                 "default",
                 "full",
@@ -662,6 +662,17 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 raise pybamm.OptionError(
                     "cannot have transverse convection in 0D model"
                 )
+        if options["dimensionality"] == 3:
+            if options["cell geometry"] not in ["pouch", "cylindrical"]:
+                raise pybamm.OptionError(
+                    "'cell geometry' must be 'pouch' or 'cylindrical' if 'dimensionality' is '3'"
+                )
+
+        if options["cell geometry"] == "cylindrical":
+            if options["dimensionality"] != 3:
+                raise pybamm.OptionError(
+                    "'dimensionality' must be '3' if 'cell geometry' is 'cylindrical'"
+                )
 
         if (
             options["thermal"] in ["x-lumped", "x-full"]
@@ -936,7 +947,12 @@ class BaseBatteryModel(pybamm.BaseModel):
 
     @property
     def default_geometry(self):
-        return pybamm.battery_geometry(options=self.options)
+        if self.options["cell geometry"] == "cylindrical":
+            return pybamm.battery_geometry(
+                options=self.options, form_factor="cylindrical"
+            )
+        else:
+            return pybamm.battery_geometry(options=self.options)
 
     @property
     def default_var_pts(self):
@@ -990,6 +1006,17 @@ class BaseBatteryModel(pybamm.BaseModel):
 
         elif self.options["dimensionality"] == 2:
             base_submeshes["current collector"] = pybamm.ScikitUniform2DSubMesh
+        elif self.options["dimensionality"] == 3:
+            base_submeshes["current collector"] = pybamm.SubMesh0D
+            geom_type = self.options.get("cell geometry", "pouch")
+            if geom_type == "pouch":
+                base_submeshes["cell"] = pybamm.ScikitFemGenerator3D(
+                    geom_type="pouch", h="0.1"
+                )
+            elif geom_type == "cylindrical":
+                base_submeshes["cell"] = pybamm.ScikitFemGenerator3D(
+                    geom_type="cylinder", h="0.1"
+                )
         return base_submeshes
 
     @property
@@ -1018,6 +1045,11 @@ class BaseBatteryModel(pybamm.BaseModel):
             base_spatial_methods["current collector"] = pybamm.FiniteVolume()
         elif self.options["dimensionality"] == 2:
             base_spatial_methods["current collector"] = pybamm.ScikitFiniteElement()
+        elif self.options["dimensionality"] == 3:
+            base_spatial_methods["current collector"] = (
+                pybamm.ZeroDimensionalSpatialMethod()
+            )
+            base_spatial_methods["cell"] = pybamm.ScikitFiniteElement3D()
         return base_spatial_methods
 
     @property
diff --git a/src/pybamm/models/full_battery_models/lithium_ion/__init__.py b/src/pybamm/models/full_battery_models/lithium_ion/__init__.py
index 556e8de31..7b2f6d9be 100644
--- a/src/pybamm/models/full_battery_models/lithium_ion/__init__.py
+++ b/src/pybamm/models/full_battery_models/lithium_ion/__init__.py
@@ -19,6 +19,7 @@ from .dfn import DFN
 from .newman_tobias import NewmanTobias
 from .basic_dfn import BasicDFN
 from .basic_spm import BasicSPM
+from .basic_spm_with_3d_thermal import Basic3DThermalSPM
 from .basic_dfn_half_cell import BasicDFNHalfCell
 from .basic_dfn_composite import BasicDFNComposite
 from .Yang2017 import Yang2017
@@ -29,4 +30,4 @@ from .basic_splitOCVR import SplitOCVR
 __all__ = ['Yang2017', 'base_lithium_ion_model', 'basic_dfn',
            'basic_dfn_composite', 'basic_dfn_half_cell', 'basic_spm', 'dfn',
            'electrode_soh', 'electrode_soh_half_cell', 'mpm', 'msmr',
-           'newman_tobias', 'spm', 'spme', 'basic_splitOCVR']
+           'newman_tobias', 'spm', 'spme', 'basic_splitOCVR', 'basic_spm_with_3d_thermal']
diff --git a/src/pybamm/models/full_battery_models/lithium_ion/basic_spm_with_3d_thermal.py b/src/pybamm/models/full_battery_models/lithium_ion/basic_spm_with_3d_thermal.py
new file mode 100644
index 000000000..e39b01839
--- /dev/null
+++ b/src/pybamm/models/full_battery_models/lithium_ion/basic_spm_with_3d_thermal.py
@@ -0,0 +1,296 @@
+import pybamm
+
+from .base_lithium_ion_model import BaseModel
+
+
+class Basic3DThermalSPM(BaseModel):
+    """Single Particle Model (SPM) model of a lithium-ion battery, from
+    :footcite:t:`Marquis2019`.
+
+    This class differs from the :class:`pybamm.lithium_ion.SPM` model class in that it
+    shows the whole model in a single class. This comes at the cost of flexibility in
+    combining different physical effects, and in general the main SPM class should be
+    used instead.
+
+    The model is based on the SPM, but with a separate cell domain for the temperature
+    variable. This allows for a more detailed treatment of the thermal effects in the
+    cell, as the temperature can vary independently of the other variables in the
+    model.
+
+    Parameters
+    ----------
+    name : str, optional
+        The name of the model.
+    """
+
+    def __init__(self, options=None, name="SPM with Separate Cell Domain"):
+        super().__init__(options, name)
+        pybamm.citations.register("Marquis2019")
+        # `param` is a class containing all the relevant parameters and functions for
+        # this model. These are purely symbolic at this stage, and will be set by the
+        # `ParameterValues` class when the model is processed.
+
+        ######################
+        # Variables
+        ######################
+        # Variables that depend on time only are created without a domain
+        Q = pybamm.Variable("Discharge capacity [A.h]")
+        # Variables that vary spatially are created with a domain
+        c_s_n = pybamm.Variable(
+            "X-averaged negative particle concentration [mol.m-3]",
+            domain="negative particle",
+        )
+        c_s_p = pybamm.Variable(
+            "X-averaged positive particle concentration [mol.m-3]",
+            domain="positive particle",
+        )
+
+        T = pybamm.Variable("Cell temperature [K]", domain="cell")
+
+        if self.options.get("cell geometry") == "pouch":
+            x = pybamm.SpatialVariable("x", domain="cell")
+            y = pybamm.SpatialVariable("y", domain="cell")
+            z = pybamm.SpatialVariable("z", domain="cell")
+            integration_vars = [x, y, z]
+        elif self.options.get("cell geometry") == "cylindrical":
+            r = pybamm.SpatialVariable("r_macro", domain="cell")
+            z = pybamm.SpatialVariable("z", domain="cell")
+            integration_vars = [r, z]
+        else:
+            raise ValueError(
+                f"Geometry type '{self.options.get('cell geometry')}' is not supported. "
+                "Supported types are 'pouch' and 'cylindrical'."
+            )  # pragma: no cover
+
+        volume = pybamm.Integral(pybamm.PrimaryBroadcast(1.0, "cell"), integration_vars)
+        T_av = pybamm.Integral(T, integration_vars) / volume
+
+        ######################
+        # Other set-up
+        ######################
+
+        # Current density
+        i_cell = self.param.current_density_with_time
+        a_n = 3 * self.param.n.prim.epsilon_s_av / self.param.n.prim.R_typ
+        a_p = 3 * self.param.p.prim.epsilon_s_av / self.param.p.prim.R_typ
+        j_n = i_cell / (self.param.n.L * a_n)
+        j_p = -i_cell / (self.param.p.L * a_p)
+
+        ######################
+        # State of Charge
+        ######################
+        I = self.param.current_with_time
+        # The `rhs` dictionary contains differential equations, with the key being the
+        # variable in the d/dt
+        self.rhs[Q] = I / 3600
+        # Initial conditions must be provided for the ODEs
+        self.initial_conditions[Q] = pybamm.Scalar(0)
+
+        ######################
+        # Particles
+        ######################
+
+        N_s_n = -self.param.n.prim.D(c_s_n, T_av) * pybamm.grad(c_s_n)
+        N_s_p = -self.param.p.prim.D(c_s_p, T_av) * pybamm.grad(c_s_p)
+        self.rhs[c_s_n] = -pybamm.div(N_s_n)
+        self.rhs[c_s_p] = -pybamm.div(N_s_p)
+
+        # Surf takes the surface value of a variable, i.e. its boundary value on the
+        # right side. This is also accessible via `boundary_value(x, "right")`, with
+        # "left" providing the boundary value of the left side
+        c_s_surf_n = pybamm.surf(c_s_n)
+        c_s_surf_p = pybamm.surf(c_s_p)
+
+        # Boundary conditions must be provided for equations with spatial derivatives
+        self.boundary_conditions[c_s_n] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (
+                -j_n / (self.param.F * pybamm.surf(self.param.n.prim.D(c_s_n, T_av))),
+                "Neumann",
+            ),
+        }
+        self.boundary_conditions[c_s_p] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (
+                -j_p / (self.param.F * pybamm.surf(self.param.p.prim.D(c_s_p, T_av))),
+                "Neumann",
+            ),
+        }
+
+        # c_n_init and c_p_init are functions of r and x, but for the SPM we
+        # take the x-averaged value since there is no x-dependence in the particles
+        self.initial_conditions[c_s_n] = pybamm.x_average(self.param.n.prim.c_init)
+        self.initial_conditions[c_s_p] = pybamm.x_average(self.param.p.prim.c_init)
+        # Events specify points at which a solution should terminate
+        sto_surf_n = c_s_surf_n / self.param.n.prim.c_max
+        sto_surf_p = c_s_surf_p / self.param.p.prim.c_max
+
+        self.events += [
+            pybamm.Event(
+                "Minimum negative particle surface stoichiometry",
+                pybamm.min(sto_surf_n) - 0.01,
+            ),
+            pybamm.Event(
+                "Maximum negative particle surface stoichiometry",
+                (1 - 0.01) - pybamm.max(sto_surf_n),
+            ),
+            pybamm.Event(
+                "Minimum positive particle surface stoichiometry",
+                pybamm.min(sto_surf_p) - 0.01,
+            ),
+            pybamm.Event(
+                "Maximum positive particle surface stoichiometry",
+                (1 - 0.01) - pybamm.max(sto_surf_p),
+            ),
+        ]
+
+        # Note that the SPM does not have any algebraic equations, so the `algebraic`
+        # dictionary remains empty
+
+        ######################
+        # (Some) variables
+        ######################
+        # Interfacial reactions
+        RT_F = self.param.R * T_av / self.param.F
+        j0_n = self.param.n.prim.j0(self.param.c_e_init_av, c_s_surf_n, T_av)
+        j0_p = self.param.p.prim.j0(self.param.c_e_init_av, c_s_surf_p, T_av)
+        eta_n = (2 / self.param.n.prim.ne) * RT_F * pybamm.arcsinh(j_n / (2 * j0_n))
+        eta_p = (2 / self.param.p.prim.ne) * RT_F * pybamm.arcsinh(j_p / (2 * j0_p))
+        phi_s_n = 0
+        phi_e = -eta_n - self.param.n.prim.U(sto_surf_n, T_av)
+        phi_s_p = eta_p + phi_e + self.param.p.prim.U(sto_surf_p, T_av)
+        V = phi_s_p
+        num_cells = pybamm.Parameter(
+            "Number of cells connected in series to make a battery"
+        )
+
+        # Compute heat source
+
+        # Reversible (entropic) heating
+        dUdT_n = self.param.n.prim.dUdT(sto_surf_n)
+        dUdT_p = self.param.p.prim.dUdT(sto_surf_p)
+        Q_rev_n = a_n * j_n * T_av * dUdT_n
+        Q_rev_p = a_p * j_p * T_av * dUdT_p
+
+        # Irreversible (Joule) heating
+        Q_irr_n = a_n * j_n * eta_n
+        Q_irr_p = a_p * j_p * eta_p
+
+        # Total heating in each electrode domain
+        Q_total_n = Q_rev_n + Q_irr_n
+        Q_total_p = Q_rev_p + Q_irr_p
+
+        # Average the 1D heat source and broadcast to the 3D 'cell' domain
+        L_n = self.param.n.L
+        L_p = self.param.p.L
+        L_x = self.param.L_x  # Total cell thickness
+
+        # Calculate the true volume-weighted average heat source in W/m^3
+        Q_vol_avg = (Q_total_n * L_n + Q_total_p * L_p) / L_x
+
+        # Broadcast this uniform volumetric heat source to the entire 3D cell domain
+        Q_source = pybamm.PrimaryBroadcast(Q_vol_avg, "cell")
+
+        # Wrap in source term to get the correct mass matrix
+        Q_source = pybamm.source(Q_source, T)
+
+        # Define the 3D heat equation
+        # Effective parameters are functions of temperature
+        rho_c_p_eff = self.param.rho_c_p_eff(T)
+        lambda_eff = self.param.lambda_eff(T)
+
+        # The heat equation is d(rho*cp*T)/dt = div(lambda*grad(T)) + Q
+        term1 = lambda_eff * pybamm.laplacian(T)
+        term2 = pybamm.inner(pybamm.grad(lambda_eff), pybamm.grad(T))
+
+        self.rhs[T] = (term1 + term2 + Q_source) / rho_c_p_eff
+
+        # Cooling taken care of by boundary conditions
+
+        self.initial_conditions[T] = pybamm.PrimaryBroadcast(self.param.T_init, "cell")
+        self.set_thermal_bcs(T)
+
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        # The `variables` dictionary contains all variables that might be useful for
+        # visualising the solution of the model
+        # Primary broadcasts are used to broadcast scalar quantities across a domain
+        # into a vector of the right shape, for multiplying with other vectors
+        self.variables = {
+            "Time [s]": pybamm.t,
+            "Discharge capacity [A.h]": Q,
+            "X-averaged negative particle concentration [mol.m-3]": c_s_n,
+            "Negative particle surface "
+            "concentration [mol.m-3]": pybamm.PrimaryBroadcast(
+                c_s_surf_n, "negative electrode"
+            ),
+            "Electrolyte concentration [mol.m-3]": pybamm.PrimaryBroadcast(
+                self.param.c_e_init_av, whole_cell
+            ),
+            "X-averaged positive particle concentration [mol.m-3]": c_s_p,
+            "Positive particle surface "
+            "concentration [mol.m-3]": pybamm.PrimaryBroadcast(
+                c_s_surf_p, "positive electrode"
+            ),
+            "Current [A]": I,
+            "Current variable [A]": I,  # for compatibility with pybamm.Experiment
+            "Negative electrode potential [V]": pybamm.PrimaryBroadcast(
+                phi_s_n, "negative electrode"
+            ),
+            "Electrolyte potential [V]": pybamm.PrimaryBroadcast(phi_e, whole_cell),
+            "Positive electrode potential [V]": pybamm.PrimaryBroadcast(
+                phi_s_p, "positive electrode"
+            ),
+            "Voltage [V]": V,
+            "Battery voltage [V]": V * num_cells,
+        }
+
+        # Add new thermal variables
+        self.variables.update(
+            {
+                "Cell temperature [K]": T,
+                "Volume-averaged cell temperature [K]": T_av,
+            }
+        )
+
+        # Events specify points at which a solution should terminate
+        self.events += [
+            pybamm.Event("Minimum voltage [V]", V - self.param.voltage_low_cut),
+            pybamm.Event("Maximum voltage [V]", self.param.voltage_high_cut - V),
+        ]
+
+    def set_thermal_bcs(self, T):
+        geometry_type = self.options.get("cell geometry", "pouch")
+        if geometry_type == "pouch":
+            y = pybamm.SpatialVariable("y", "cell")
+            z = pybamm.SpatialVariable("z", "cell")
+            T_amb = self.param.T_amb(y, z, pybamm.t)
+            face_params = {
+                "x_min": self.param.h_edge_x_min,
+                "x_max": self.param.h_edge_x_max,
+                "y_min": self.param.h_edge_y_min,
+                "y_max": self.param.h_edge_y_max,
+                "z_min": self.param.h_edge_z_min,
+                "z_max": self.param.h_edge_z_max,
+            }
+        elif geometry_type == "cylindrical":
+            r = pybamm.SpatialVariable("r_macro", "cell")
+            z = pybamm.SpatialVariable("z", "cell")
+            T_amb = self.param.T_amb(r, z, pybamm.t)
+            face_params = {
+                "r_min": self.param.h_edge_radial_min,
+                "r_max": self.param.h_edge_radial_max,
+                "z_min": self.param.h_edge_z_min,
+                "z_max": self.param.h_edge_z_max,
+            }
+        else:
+            raise ValueError(
+                f"Geometry type '{geometry_type}' is not supported. "
+                "Supported types are 'pouch' and 'cylindrical'."
+            )  # pragma: no cover
+
+        self.boundary_conditions[T] = {}
+        lambda_eff = self.param.lambda_eff(T)
+
+        for face, h_coeff in face_params.items():
+            q_face = -h_coeff * (T - T_amb)
+            self.boundary_conditions[T][face] = (q_face / lambda_eff, "Neumann")
diff --git a/src/pybamm/parameters/lithium_ion_parameters.py b/src/pybamm/parameters/lithium_ion_parameters.py
index edabd4c0c..17e31f6d5 100644
--- a/src/pybamm/parameters/lithium_ion_parameters.py
+++ b/src/pybamm/parameters/lithium_ion_parameters.py
@@ -58,6 +58,18 @@ class LithiumIonParameters(BaseParameters):
         self.lambda_eff = self.therm.lambda_eff
         self.cell_heat_capacity = self.therm.cell_heat_capacity
 
+        # pouch bcs
+        self.h_edge_x_min = self.therm.h_edge_x_min
+        self.h_edge_x_max = self.therm.h_edge_x_max
+        self.h_edge_y_max = self.therm.h_edge_y_max
+        self.h_edge_y_min = self.therm.h_edge_y_min
+        self.h_edge_z_min = self.therm.h_edge_z_min
+        self.h_edge_z_max = self.therm.h_edge_z_max
+
+        # cylinder specific bcs
+        self.h_edge_radial_min = self.therm.h_edge_radial_min
+        self.h_edge_radial_max = self.therm.h_edge_radial_max
+
         # Macroscale geometry
         self.L_x = self.geo.L_x
         self.L = self.geo.L
diff --git a/src/pybamm/parameters/thermal_parameters.py b/src/pybamm/parameters/thermal_parameters.py
index 10db1ba26..0ddd0a3b4 100644
--- a/src/pybamm/parameters/thermal_parameters.py
+++ b/src/pybamm/parameters/thermal_parameters.py
@@ -44,6 +44,34 @@ class ThermalParameters(BaseParameters):
         # Lumped model heat capacity (volumetric)
         self.cell_heat_capacity = pybamm.Parameter("Cell heat capacity [J.K-1.m-3]")
 
+        # Box geometry faces
+        self.h_edge_x_min = pybamm.Parameter(
+            "Left face heat transfer coefficient [W.m-2.K-1]"
+        )
+        self.h_edge_x_max = pybamm.Parameter(
+            "Right face heat transfer coefficient [W.m-2.K-1]"
+        )
+        self.h_edge_y_max = pybamm.Parameter(
+            "Front face heat transfer coefficient [W.m-2.K-1]"
+        )
+        self.h_edge_y_min = pybamm.Parameter(
+            "Back face heat transfer coefficient [W.m-2.K-1]"
+        )
+        self.h_edge_z_min = pybamm.Parameter(
+            "Bottom face heat transfer coefficient [W.m-2.K-1]"
+        )
+        self.h_edge_z_max = pybamm.Parameter(
+            "Top face heat transfer coefficient [W.m-2.K-1]"
+        )
+
+        # Cylinder specific geometry faces
+        self.h_edge_radial_min = pybamm.Parameter(
+            "Inner radius heat transfer coefficient [W.m-2.K-1]"
+        )
+        self.h_edge_radial_max = pybamm.Parameter(
+            "Outer radius heat transfer coefficient [W.m-2.K-1]"
+        )
+
     def T_amb(self, y, z, t):
         """Ambient temperature [K]"""
         return pybamm.FunctionParameter(
diff --git a/src/pybamm/solvers/processed_variable.py b/src/pybamm/solvers/processed_variable.py
index 349f4c6ac..99cc200fe 100644
--- a/src/pybamm/solvers/processed_variable.py
+++ b/src/pybamm/solvers/processed_variable.py
@@ -1162,6 +1162,110 @@ class ProcessedVariable3DSciKitFEM(ProcessedVariable3D):
         return entries
 
 
+class ProcessedVariableUnstructured(ProcessedVariable):
+    """
+    A processed variable for data on an unstructured mesh (e.g., from a FEM solution).
+    This class correctly uses scipy's LinearNDInterpolator for spatial interpolation,
+    which is the required method for scattered data points from a FEM mesh.
+
+    Parameters
+    ----------
+    name : str
+        The name of the variable
+    base_variables : list of :class:`pybamm.Symbol`
+        A list of base variables with a method `evaluate(t,y)`, each entry of which
+        returns the value of that variable for that particular sub-solution.
+        A Solution can be comprised of sub-solutions which are the solutions of
+        different models.
+        base_variables : list of :class:`pybamm.Symbol`
+        A list of base variables with a method `evaluate(t,y)`, each entry of which
+        returns the value of that variable for that particular sub-solution.
+        A Solution can be comprised of sub-solutions which are the solutions of
+        different models.
+        Note that this can be any kind of node in the expression tree, not
+        just a :class:`pybamm.Variable`.
+        When evaluated, returns an array of size (m,n)
+    base_variables_casadi : list of :class:`casadi.Function`
+        A list of casadi functions. When evaluated, returns the same thing as
+        `base_Variable.evaluate` (but more efficiently).
+    solution : :class:`pybamm.Solution`
+        The solution object to be used to create the processed variables
+    time_integral : pybamm.ProcessedVariableTimeIntegral, optional
+        An optional time integral object to handle time integration of the variable.
+        If provided, the processed variable will handle time integration using this object.
+    """
+
+    def __init__(
+        self,
+        name: str,
+        base_variables,
+        base_variables_casadi,
+        solution,
+        time_integral: pybamm.ProcessedVariableTimeIntegral | None = None,
+    ):
+        self.dimensions = base_variables[0].mesh.dimension
+        super().__init__(
+            name,
+            base_variables,
+            base_variables_casadi,
+            solution,
+            time_integral=time_integral,
+        )
+        self._time_interpolator = None
+
+    def initialise(self):
+        if self.entries_raw_initialized:
+            return
+        self._entries_raw = self.observe_raw()
+
+        from scipy.interpolate import interp1d
+
+        self._time_interpolator = interp1d(
+            self.t_pts,
+            self._entries_raw,
+            kind="linear",
+            axis=1,  # Interpolate along the time axis
+            bounds_error=False,
+            fill_value="extrapolate",
+        )
+
+    def _shape(self, t):
+        return [self.mesh.npts, len(t)]
+
+    def __call__(self, t, x=None, y=None, z=None, r=None, R=None, fill_value=np.nan):
+        from scipy.interpolate import LinearNDInterpolator
+
+        self.initialise()
+
+        data_at_t = self._time_interpolator(t)
+
+        spatial_coords_provided = any(c is not None for c in [x, y, z])
+        if not spatial_coords_provided:
+            return data_at_t  # Return all node data if no spatial coords are given
+
+        node_coords = self.mesh.nodes
+
+        eval_points = np.column_stack([x.ravel(), y.ravel(), z.ravel()])
+        output_shape = x.shape
+
+        # If t was a single value (scalar), data_at_t is 1D array of (n_nodes)
+        if isinstance(t, int | float):
+            spatial_interpolator = LinearNDInterpolator(
+                points=node_coords, values=data_at_t, fill_value=fill_value
+            )
+            result = spatial_interpolator(eval_points)
+        else:  # If t was a vector, we must create an interpolator for each time step
+            result = np.empty((len(eval_points), len(t)))
+            for i in range(len(t)):
+                spatial_interpolator = LinearNDInterpolator(
+                    points=node_coords, values=data_at_t[:, i], fill_value=fill_value
+                )
+                result[:, i] = spatial_interpolator(eval_points)
+
+        final_shape = (*output_shape, len(t) if not isinstance(t, int | float) else 1)
+        return result.reshape(final_shape).squeeze()
+
+
 def process_variable(name: str, base_variables, *args, **kwargs):
     mesh = base_variables[0].mesh
     domain = base_variables[0].domain
@@ -1170,6 +1274,9 @@ def process_variable(name: str, base_variables, *args, **kwargs):
     base_eval_shape = base_variables[0].shape
     base_eval_size = base_variables[0].size
 
+    if isinstance(mesh, pybamm.ScikitFemSubMesh3D):
+        return ProcessedVariableUnstructured(name, base_variables, *args, **kwargs)
+
     # handle 2D (in space) finite element variables differently
     if (
         mesh
diff --git a/src/pybamm/spatial_methods/scikit_finite_element_3d.py b/src/pybamm/spatial_methods/scikit_finite_element_3d.py
index 0e0ad4558..d6e972fc5 100644
--- a/src/pybamm/spatial_methods/scikit_finite_element_3d.py
+++ b/src/pybamm/spatial_methods/scikit_finite_element_3d.py
@@ -49,7 +49,7 @@ class ScikitFiniteElement3D(pybamm.SpatialMethod):
             entries = mesh.nodes[:, 1][:, np.newaxis]
         elif symbol.name == "z":
             entries = mesh.nodes[:, 2][:, np.newaxis]
-        elif symbol.name == "r":
+        elif symbol.name == "r" or symbol.name == "r_macro":
             entries = np.sqrt(mesh.nodes[:, 0] ** 2 + mesh.nodes[:, 1] ** 2)[
                 :, np.newaxis
             ]


Test patch:
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_basic_spm_with_3d.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_basic_spm_with_3d.py
new file mode 100644
index 000000000..e7c3e39e5
--- /dev/null
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_basic_spm_with_3d.py
@@ -0,0 +1,158 @@
+import numpy as np
+
+import pybamm
+
+
+class TestBasicSPMWith3DThermal:
+    def test_spm_3d_vs_lumped_pouch(self):
+        models = {
+            "Lumped": pybamm.lithium_ion.SPM(options={"thermal": "lumped"}),
+            "3D": pybamm.lithium_ion.Basic3DThermalSPM(
+                options={"cell geometry": "pouch", "dimensionality": 3}
+            ),
+        }
+
+        parameter_values = pybamm.ParameterValues("Marquis2019")
+        h_values = [0.1, 1, 10, 100]
+
+        experiment = pybamm.Experiment(
+            [
+                ("Discharge at 3C until 2.8V", "Rest for 10 minutes"),
+            ]
+        )
+
+        var_pts = {
+            "x_n": 20,
+            "x_s": 20,
+            "x_p": 20,
+            "r_n": 30,
+            "r_p": 30,
+            "x": None,
+            "y": None,
+            "z": None,
+        }
+
+        all_solutions = {}
+
+        for h in h_values:
+            h_params = parameter_values.copy()
+            h_params.update(
+                {
+                    "Total heat transfer coefficient [W.m-2.K-1]": h,
+                    "Left face heat transfer coefficient [W.m-2.K-1]": h,
+                    "Right face heat transfer coefficient [W.m-2.K-1]": h,
+                    "Front face heat transfer coefficient [W.m-2.K-1]": h,
+                    "Back face heat transfer coefficient [W.m-2.K-1]": h,
+                    "Bottom face heat transfer coefficient [W.m-2.K-1]": h,
+                    "Top face heat transfer coefficient [W.m-2.K-1]": h,
+                },
+                check_already_exists=False,
+            )
+
+            solutions = {}
+            for model_name, model in models.items():
+                sim = pybamm.Simulation(
+                    model,
+                    parameter_values=h_params,
+                    var_pts=var_pts,
+                    experiment=experiment,
+                )
+                solutions[model_name] = sim.solve()
+
+            all_solutions[h] = solutions
+
+        for _h, solutions in all_solutions.items():
+            lumped_sol = solutions["Lumped"]
+            three_d_sol = solutions["3D"]
+
+            np.testing.assert_allclose(lumped_sol.t[-1], three_d_sol.t[-1], rtol=0.01)
+
+            lumped_temp_final = lumped_sol[
+                "Volume-averaged cell temperature [K]"
+            ].entries[-1]
+            three_d_temp_final = three_d_sol[
+                "Volume-averaged cell temperature [K]"
+            ].entries[-1]
+            np.testing.assert_allclose(lumped_temp_final, three_d_temp_final, rtol=0.02)
+
+            lumped_final_voltage = lumped_sol["Voltage [V]"].entries[-1]
+            three_d_final_voltage = three_d_sol["Voltage [V]"].entries[-1]
+            np.testing.assert_allclose(
+                lumped_final_voltage, three_d_final_voltage, rtol=0.01
+            )
+
+    def test_spm_3d_vs_lumped_cylinder(self):
+        models = {
+            "Lumped": pybamm.lithium_ion.SPM(options={"thermal": "lumped"}),
+            "3D": pybamm.lithium_ion.Basic3DThermalSPM(
+                options={"cell geometry": "cylindrical", "dimensionality": 3}
+            ),
+        }
+
+        parameter_values = pybamm.ParameterValues("NCA_Kim2011")
+        h_values = [0.1, 1, 10, 100]
+
+        experiment = pybamm.Experiment(
+            [
+                ("Discharge at 3C until 2.8V", "Rest for 10 minutes"),
+            ]
+        )
+
+        var_pts = {
+            "x_n": 20,
+            "x_s": 20,
+            "x_p": 20,
+            "r_n": 30,
+            "r_p": 30,
+            "r_macro": None,
+            "y": None,
+            "z": None,
+        }
+
+        all_solutions = {}
+
+        for h in h_values:
+            h_params = parameter_values.copy()
+            h_params.update(
+                {
+                    "Inner cell radius [m]": 0.005,
+                    "Outer cell radius [m]": 0.018,
+                    "Total heat transfer coefficient [W.m-2.K-1]": h,
+                    "Outer radius heat transfer coefficient [W.m-2.K-1]": h,
+                    "Inner radius heat transfer coefficient [W.m-2.K-1]": h,
+                    "Bottom face heat transfer coefficient [W.m-2.K-1]": h,
+                    "Top face heat transfer coefficient [W.m-2.K-1]": h,
+                },
+                check_already_exists=False,
+            )
+
+            solutions = {}
+            for model_name, model in models.items():
+                sim = pybamm.Simulation(
+                    model,
+                    parameter_values=h_params,
+                    var_pts=var_pts,
+                    experiment=experiment,
+                )
+                solutions[model_name] = sim.solve()
+
+            all_solutions[h] = solutions
+
+        for _h, solutions in all_solutions.items():
+            lumped_sol = solutions["Lumped"]
+            three_d_sol = solutions["3D"]
+
+            np.testing.assert_allclose(lumped_sol.t[-1], three_d_sol.t[-1], rtol=0.01)
+            lumped_temp_final = lumped_sol[
+                "Volume-averaged cell temperature [K]"
+            ].entries[-1]
+            three_d_temp_final = three_d_sol[
+                "Volume-averaged cell temperature [K]"
+            ].entries[-1]
+            np.testing.assert_allclose(lumped_temp_final, three_d_temp_final, rtol=0.02)
+
+            lumped_final_voltage = lumped_sol["Voltage [V]"].entries[-1]
+            three_d_final_voltage = three_d_sol["Voltage [V]"].entries[-1]
+            np.testing.assert_allclose(
+                lumped_final_voltage, three_d_final_voltage, rtol=0.01
+            )
diff --git a/tests/shared.py b/tests/shared.py
index 317561e22..680c629ea 100644
--- a/tests/shared.py
+++ b/tests/shared.py
@@ -116,8 +116,8 @@ def get_mesh_for_testing(
     return pybamm.Mesh(geometry, submesh_types, var_pts)
 
 
-def get_unit_3d_mesh_for_testing(geom_type="box", **geom_params):
-    if geom_type == "box":
+def get_unit_3d_mesh_for_testing(geom_type="pouch", **geom_params):
+    if geom_type == "pouch":
         x = pybamm.SpatialVariable("x", ["current collector"])
         y = pybamm.SpatialVariable("y", ["current collector"])
         z = pybamm.SpatialVariable("z", ["current collector"])
@@ -160,7 +160,7 @@ def get_unit_3d_mesh_for_testing(geom_type="box", **geom_params):
 
 
 def get_3d_mesh_for_testing(
-    xpts=5, ypts=5, zpts=5, geom_type="box", include_particles=False, **geom_params
+    xpts=5, ypts=5, zpts=5, geom_type="pouch", include_particles=False, **geom_params
 ):
     param = pybamm.ParameterValues(
         {
diff --git a/tests/unit/test_experiments/test_experiment.py b/tests/unit/test_experiments/test_experiment.py
index 181632185..af7bae075 100644
--- a/tests/unit/test_experiments/test_experiment.py
+++ b/tests/unit/test_experiments/test_experiment.py
@@ -360,3 +360,156 @@ class TestExperiment:
         result = np.array(f(test_points)).flatten()
 
         np.testing.assert_allclose(result, expected, rtol=1e-7, atol=1e-6)
+
+    def test_spm_3d_vs_lumped_pouch(self):
+        models = {
+            "Lumped": pybamm.lithium_ion.SPM(options={"thermal": "lumped"}),
+            "3D": pybamm.lithium_ion.Basic3DThermalSPM(
+                options={"cell geometry": "pouch", "dimensionality": 3}
+            ),
+        }
+
+        parameter_values = pybamm.ParameterValues("Marquis2019")
+        h_values = [0.1, 1, 10, 100]
+
+        experiment = pybamm.Experiment(
+            [
+                ("Discharge at 3C until 2.8V", "Rest for 10 minutes"),
+            ]
+        )
+
+        var_pts = {
+            "x_n": 20,
+            "x_s": 20,
+            "x_p": 20,
+            "r_n": 30,
+            "r_p": 30,
+            "x": None,
+            "y": None,
+            "z": None,
+        }
+
+        all_solutions = {}
+
+        for h in h_values:
+            h_params = parameter_values.copy()
+            h_params.update(
+                {
+                    "Total heat transfer coefficient [W.m-2.K-1]": h,
+                    "Left face heat transfer coefficient [W.m-2.K-1]": h,
+                    "Right face heat transfer coefficient [W.m-2.K-1]": h,
+                    "Front face heat transfer coefficient [W.m-2.K-1]": h,
+                    "Back face heat transfer coefficient [W.m-2.K-1]": h,
+                    "Bottom face heat transfer coefficient [W.m-2.K-1]": h,
+                    "Top face heat transfer coefficient [W.m-2.K-1]": h,
+                },
+                check_already_exists=False,
+            )
+
+            solutions = {}
+            for model_name, model in models.items():
+                sim = pybamm.Simulation(
+                    model,
+                    parameter_values=h_params,
+                    var_pts=var_pts,
+                    experiment=experiment,
+                )
+                solutions[model_name] = sim.solve()
+
+            all_solutions[h] = solutions
+
+        for _h, solutions in all_solutions.items():
+            lumped_sol = solutions["Lumped"]
+            three_d_sol = solutions["3D"]
+
+            np.testing.assert_allclose(lumped_sol.t[-1], three_d_sol.t[-1], rtol=0.01)
+
+            lumped_temp_final = lumped_sol[
+                "Volume-averaged cell temperature [K]"
+            ].entries[-1]
+            three_d_temp_final = three_d_sol[
+                "Volume-averaged cell temperature [K]"
+            ].entries[-1]
+            np.testing.assert_allclose(lumped_temp_final, three_d_temp_final, rtol=0.02)
+
+            lumped_final_voltage = lumped_sol["Voltage [V]"].entries[-1]
+            three_d_final_voltage = three_d_sol["Voltage [V]"].entries[-1]
+            np.testing.assert_allclose(
+                lumped_final_voltage, three_d_final_voltage, rtol=0.01
+            )
+
+    def test_spm_3d_vs_lumped_cylinder(self):
+        models = {
+            "Lumped": pybamm.lithium_ion.SPM(options={"thermal": "lumped"}),
+            "3D": pybamm.lithium_ion.Basic3DThermalSPM(
+                options={"cell geometry": "cylindrical", "dimensionality": 3}
+            ),
+        }
+
+        parameter_values = pybamm.ParameterValues("NCA_Kim2011")
+        h_values = [0.1, 1, 10, 100]
+
+        experiment = pybamm.Experiment(
+            [
+                ("Discharge at 3C until 2.8V", "Rest for 10 minutes"),
+            ]
+        )
+
+        var_pts = {
+            "x_n": 20,
+            "x_s": 20,
+            "x_p": 20,
+            "r_n": 30,
+            "r_p": 30,
+            "r_macro": None,
+            "y": None,
+            "z": None,
+        }
+
+        all_solutions = {}
+
+        for h in h_values:
+            h_params = parameter_values.copy()
+            h_params.update(
+                {
+                    "Inner cell radius [m]": 0.005,
+                    "Outer cell radius [m]": 0.018,
+                    "Total heat transfer coefficient [W.m-2.K-1]": h,
+                    "Outer radius heat transfer coefficient [W.m-2.K-1]": h,
+                    "Inner radius heat transfer coefficient [W.m-2.K-1]": h,
+                    "Bottom face heat transfer coefficient [W.m-2.K-1]": h,
+                    "Top face heat transfer coefficient [W.m-2.K-1]": h,
+                },
+                check_already_exists=False,
+            )
+
+            solutions = {}
+            for model_name, model in models.items():
+                sim = pybamm.Simulation(
+                    model,
+                    parameter_values=h_params,
+                    var_pts=var_pts,
+                    experiment=experiment,
+                )
+                solutions[model_name] = sim.solve()
+
+            all_solutions[h] = solutions
+
+        for _h, solutions in all_solutions.items():
+            lumped_sol = solutions["Lumped"]
+            three_d_sol = solutions["3D"]
+
+            np.testing.assert_allclose(lumped_sol.t[-1], three_d_sol.t[-1], rtol=0.01)
+            lumped_temp_final = lumped_sol[
+                "Volume-averaged cell temperature [K]"
+            ].entries[-1]
+            three_d_temp_final = three_d_sol[
+                "Volume-averaged cell temperature [K]"
+            ].entries[-1]
+            np.testing.assert_allclose(lumped_temp_final, three_d_temp_final, rtol=0.02)
+
+            lumped_final_voltage = lumped_sol["Voltage [V]"].entries[-1]
+            three_d_final_voltage = three_d_sol["Voltage [V]"].entries[-1]
+            np.testing.assert_allclose(
+                lumped_final_voltage, three_d_final_voltage, rtol=0.01
+            )
diff --git a/tests/unit/test_meshes/test_scikit_fem_submesh_3d.py b/tests/unit/test_meshes/test_scikit_fem_submesh_3d.py
index b731be0be..ee36b5a5e 100644
--- a/tests/unit/test_meshes/test_scikit_fem_submesh_3d.py
+++ b/tests/unit/test_meshes/test_scikit_fem_submesh_3d.py
@@ -5,7 +5,7 @@ import pybamm
 
 
 class TestScikitFemGenerator3D:
-    def test_box_mesh_creation(self):
+    def test_pouch_mesh_creation(self):
         try:
             from pybamm.meshes.scikit_fem_submeshes_3d import ScikitFemGenerator3D
         except ImportError:
@@ -19,7 +19,7 @@ class TestScikitFemGenerator3D:
             }
         }
 
-        mesh_gen = ScikitFemGenerator3D("box", h=0.3)
+        mesh_gen = ScikitFemGenerator3D("pouch", h=0.3)
 
         mesh = pybamm.Mesh(
             geometry,
@@ -38,7 +38,7 @@ class TestScikitFemGenerator3D:
         assert submesh.dimension == 3, "Mesh should be 3D"
         assert submesh.npts > 0, "Mesh should have nodes"
         assert submesh.nodes.shape[1] == 3, "Nodes should be 3D coordinates"
-        assert submesh.coord_sys == "cartesian", "Box mesh should be cartesian"
+        assert submesh.coord_sys == "cartesian", "Pouch mesh should be cartesian"
 
     def test_cylinder_mesh_creation(self):
         try:
@@ -104,7 +104,7 @@ class TestScikitFemSubMesh3D:
             }
         }
 
-        mesh_gen = ScikitFemGenerator3D("box", h=0.4)
+        mesh_gen = ScikitFemGenerator3D("pouch", h=0.4)
         mesh = pybamm.Mesh(
             geometry,
             {"domain": mesh_gen},
@@ -150,7 +150,7 @@ class TestScikitFemSubMesh3D:
             }
         }
 
-        mesh_gen = ScikitFemGenerator3D("box", h=0.5)
+        mesh_gen = ScikitFemGenerator3D("pouch", h=0.5)
         mesh = pybamm.Mesh(
             geometry,
             {"domain": mesh_gen},
diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
index a73a8aed4..0b2ed8221 100644
--- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
@@ -20,12 +20,12 @@ OPTIONS_DICT = {
 PRINT_OPTIONS_OUTPUT = """\
 'calculate discharge energy': 'false' (possible: ['false', 'true'])
 'calculate heat source for isothermal models': 'false' (possible: ['false', 'true'])
-'cell geometry': 'pouch' (possible: ['arbitrary', 'pouch'])
+'cell geometry': 'pouch' (possible: ['arbitrary', 'pouch', 'cylindrical'])
 'contact resistance': 'false' (possible: ['false', 'true'])
 'convection': 'none' (possible: ['none', 'uniform transverse', 'full transverse'])
 'current collector': 'uniform' (possible: ['uniform', 'potential pair', 'potential pair quite conductive'])
 'diffusivity': 'single' (possible: ['single', 'current sigmoid'])
-'dimensionality': 0 (possible: [0, 1, 2])
+'dimensionality': 0 (possible: [0, 1, 2, 3])
 'electrolyte conductivity': 'default' (possible: ['default', 'full', 'leading order', 'composite', 'integrated'])
 'exchange-current density': 'single' (possible: ['single', 'current sigmoid'])
 'heat of mixing': 'false' (possible: ['false', 'true'])
@@ -127,6 +127,15 @@ class TestBaseBatteryModel:
             model.default_submesh_types["current collector"],
             pybamm.ScikitUniform2DSubMesh,
         )
+        model = pybamm.BaseBatteryModel({"dimensionality": 3, "cell geometry": "pouch"})
+        assert issubclass(
+            model.default_submesh_types["current collector"],
+            pybamm.SubMesh0D,
+        )
+        assert isinstance(
+            model.default_submesh_types["cell"],
+            pybamm.ScikitFemGenerator3D,
+        )
 
     def test_default_var_pts(self):
         var_pts = {
@@ -170,6 +179,15 @@ class TestBaseBatteryModel:
             model.default_spatial_methods["current collector"],
             pybamm.ScikitFiniteElement,
         )
+        model = pybamm.BaseBatteryModel({"dimensionality": 3, "cell geometry": "pouch"})
+        assert isinstance(
+            model.default_spatial_methods["current collector"],
+            pybamm.ZeroDimensionalSpatialMethod,
+        )
+        assert isinstance(
+            model.default_spatial_methods["cell"],
+            pybamm.ScikitFiniteElement3D,
+        )
 
     def test_options(self):
         with pytest.raises(pybamm.OptionError, match="Option"):
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
index 0bf00040a..5ec42fc69 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
@@ -73,3 +73,29 @@ class TestSPM(BaseUnitTestLithiumIon):
         new_model_cs_eqn = list(new_model.rhs.values())[1]
         model_cs_eqn = list(model.rhs.values())[1]
         assert new_model_cs_eqn == model_cs_eqn
+
+    def test_basic_spm_with_3d_thermal_pouch(self):
+        options = {"cell geometry": "pouch", "dimensionality": 3}
+        self.model = pybamm.lithium_ion.Basic3DThermalSPM
+        self.check_well_posedness(options)
+
+    def test_basic_spm_with_3d_thermal_cylinder(self):
+        options = {"cell geometry": "cylindrical", "dimensionality": 3}
+        self.model = pybamm.lithium_ion.Basic3DThermalSPM
+        self.check_well_posedness(options)
+
+    def test_basic_spm_with_3d_thermal_incompatible_options(self):
+        options = {"cell geometry": "cylindrical", "dimensionality": 2}
+        self.model = pybamm.lithium_ion.Basic3DThermalSPM
+        with pytest.raises(
+            pybamm.OptionError,
+            match="'dimensionality' must be '3' if 'cell geometry' is 'cylindrical'",
+        ):
+            self.check_well_posedness(options)
+
+        options = {"cell geometry": "arbitrary", "dimensionality": 3}
+        with pytest.raises(
+            pybamm.OptionError,
+            match="'cell geometry' must be 'pouch' or 'cylindrical' if 'dimensionality' is '3'",
+        ):
+            self.check_well_posedness(options)
diff --git a/tests/unit/test_parameters/test_lithium_ion_parameters.py b/tests/unit/test_parameters/test_lithium_ion_parameters.py
index e9741efcc..f7223f0e2 100644
--- a/tests/unit/test_parameters/test_lithium_ion_parameters.py
+++ b/tests/unit/test_parameters/test_lithium_ion_parameters.py
@@ -109,6 +109,14 @@ class TestLithiumIonParameterValues:
         values.update(
             {
                 "Cell heat capacity [J.K-1.m-3]": 2.5e6,
+                "Left face heat transfer coefficient [W.m-2.K-1]": 5.0,
+                "Right face heat transfer coefficient [W.m-2.K-1]": 5.0,
+                "Front face heat transfer coefficient [W.m-2.K-1]": 5.0,
+                "Back face heat transfer coefficient [W.m-2.K-1]": 5.0,
+                "Bottom face heat transfer coefficient [W.m-2.K-1]": 5.0,
+                "Top face heat transfer coefficient [W.m-2.K-1]": 5.0,
+                "Inner radius heat transfer coefficient [W.m-2.K-1]": 5.0,
+                "Outer radius heat transfer coefficient [W.m-2.K-1]": 5.0,
             },
             check_already_exists=False,
         )
@@ -132,6 +140,14 @@ class TestLithiumIonParameterValues:
         # other thermal parameters
         np.testing.assert_equal(values.evaluate(param.T_init), 298.15)
         np.testing.assert_equal(values.evaluate(param.cell_heat_capacity), 2.5e6)
+        np.testing.assert_equal(values.evaluate(param.h_edge_x_min), 5.0)
+        np.testing.assert_equal(values.evaluate(param.h_edge_x_max), 5.0)
+        np.testing.assert_equal(values.evaluate(param.h_edge_y_min), 5.0)
+        np.testing.assert_equal(values.evaluate(param.h_edge_y_max), 5.0)
+        np.testing.assert_equal(values.evaluate(param.h_edge_z_min), 5.0)
+        np.testing.assert_equal(values.evaluate(param.h_edge_z_max), 5.0)
+        np.testing.assert_equal(values.evaluate(param.h_edge_radial_min), 5.0)
+        np.testing.assert_equal(values.evaluate(param.h_edge_radial_max), 5.0)
 
     def test_parameter_functions(self):
         values = pybamm.lithium_ion.BaseModel().default_parameter_values
diff --git a/tests/unit/test_solvers/test_processed_variable.py b/tests/unit/test_solvers/test_processed_variable.py
index 934c0c63a..3e0f05d13 100644
--- a/tests/unit/test_solvers/test_processed_variable.py
+++ b/tests/unit/test_solvers/test_processed_variable.py
@@ -1638,3 +1638,377 @@ class TestProcessedVariable:
         np.testing.assert_array_equal(
             computed_var(t=t_sol, x=x_sol, R=R_sol, r=r_sol).shape, (6, 7, Nx, 50)
         )
+
+    def test_processed_variable_unstructured_3d_pouch(self):
+        from pybamm.meshes.scikit_fem_submeshes_3d import ScikitFemGenerator3D
+
+        x = pybamm.SpatialVariable(
+            "x", domain=["current collector"], coord_sys="cartesian"
+        )
+        y = pybamm.SpatialVariable(
+            "y", domain=["current collector"], coord_sys="cartesian"
+        )
+        z = pybamm.SpatialVariable(
+            "z", domain=["current collector"], coord_sys="cartesian"
+        )
+
+        geometry = {
+            "current collector": {
+                x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+                y: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+                z: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+            }
+        }
+
+        submesh_types = {
+            "current collector": ScikitFemGenerator3D(geom_type="pouch", h=0.3)
+        }
+        var_pts = {x: None, y: None, z: None}
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+
+        spatial_methods = {"current collector": pybamm.ScikitFiniteElement3D()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        var = pybamm.Variable("var", domain=["current collector"])
+        disc.set_variable_slices([var])
+        var_disc = disc.process_symbol(var)
+
+        t_sol = np.linspace(0, 1, 10)
+        y_sol = np.ones((var_disc.shape[0], len(t_sol))) * np.linspace(0, 5, len(t_sol))
+
+        var_casadi = to_casadi(var_disc, y_sol)
+        model = pybamm.BaseModel()
+        model._geometry = geometry
+
+        solution = pybamm.Solution(t_sol, y_sol, model, {})
+
+        processed_var = pybamm.ProcessedVariableUnstructured(
+            "test", [var_disc], [var_casadi], solution
+        )
+
+        assert processed_var.dimensions == 3
+        assert processed_var._time_interpolator is None
+
+        processed_var.initialise()
+        assert processed_var.entries_raw_initialized
+        assert processed_var._time_interpolator is not None
+
+        nodes = var_disc.mesh.nodes
+        x_test = nodes[:5, 0]
+        y_test = nodes[:5, 1]
+        z_test = nodes[:5, 2]
+
+        result_scalar_t = processed_var(0.5, x=x_test, y=y_test, z=z_test)
+        assert result_scalar_t.shape == (5,)
+
+        result_vector_t = processed_var(t_sol[:3], x=x_test, y=y_test, z=z_test)
+        assert result_vector_t.shape == (5, 3)
+
+        result_no_coords = processed_var(0.5)
+        assert result_no_coords.shape == (var_disc.mesh.npts,)
+
+    def test_processed_variable_unstructured_scalar_vs_vector_time(self):
+        from pybamm.meshes.scikit_fem_submeshes_3d import ScikitFemGenerator3D
+
+        x = pybamm.SpatialVariable(
+            "x", domain=["current collector"], coord_sys="cartesian"
+        )
+        y = pybamm.SpatialVariable(
+            "y", domain=["current collector"], coord_sys="cartesian"
+        )
+        z = pybamm.SpatialVariable(
+            "z", domain=["current collector"], coord_sys="cartesian"
+        )
+
+        geometry = {
+            "current collector": {
+                x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+                y: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+                z: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+            }
+        }
+
+        submesh_types = {
+            "current collector": ScikitFemGenerator3D(geom_type="pouch", h=0.3)
+        }
+        var_pts = {x: None, y: None, z: None}
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+
+        spatial_methods = {"current collector": pybamm.ScikitFiniteElement3D()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        var = pybamm.Variable("var", domain=["current collector"])
+        disc.set_variable_slices([var])
+        var_disc = disc.process_symbol(var)
+
+        t_sol = np.array([0.0, 0.5, 1.0])
+        y_sol = np.ones((var_disc.shape[0], len(t_sol))) * np.array([1, 2, 3])
+
+        var_casadi = to_casadi(var_disc, y_sol)
+        model = pybamm.BaseModel()
+        model._geometry = geometry
+
+        solution = pybamm.Solution(t_sol, y_sol, model, {})
+
+        processed_var = pybamm.ProcessedVariableUnstructured(
+            "test", [var_disc], [var_casadi], solution
+        )
+
+        nodes = var_disc.mesh.nodes
+        x_test = nodes[:3, 0]
+        y_test = nodes[:3, 1]
+        z_test = nodes[:3, 2]
+
+        result_scalar_int = processed_var(0, x=x_test, y=y_test, z=z_test)
+        assert isinstance(0, int)
+        assert result_scalar_int.shape == (3,)
+
+        result_scalar_float = processed_var(0.0, x=x_test, y=y_test, z=z_test)
+        assert isinstance(0.0, float)
+        assert result_scalar_float.shape == (3,)
+
+        result_vector = processed_var(
+            np.array([0.0, 0.5]), x=x_test, y=y_test, z=z_test
+        )
+        assert result_vector.shape == (3, 2)
+
+    def test_processed_variable_unstructured_fill_value(self):
+        from pybamm.meshes.scikit_fem_submeshes_3d import ScikitFemGenerator3D
+
+        x = pybamm.SpatialVariable(
+            "x", domain=["current collector"], coord_sys="cartesian"
+        )
+        y = pybamm.SpatialVariable(
+            "y", domain=["current collector"], coord_sys="cartesian"
+        )
+        z = pybamm.SpatialVariable(
+            "z", domain=["current collector"], coord_sys="cartesian"
+        )
+
+        geometry = {
+            "current collector": {
+                x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+                y: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+                z: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+            }
+        }
+
+        submesh_types = {
+            "current collector": ScikitFemGenerator3D(geom_type="pouch", h=0.3)
+        }
+        var_pts = {x: None, y: None, z: None}
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+
+        spatial_methods = {"current collector": pybamm.ScikitFiniteElement3D()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        var = pybamm.Variable("var", domain=["current collector"])
+        disc.set_variable_slices([var])
+        var_disc = disc.process_symbol(var)
+
+        t_sol = np.array([0])
+        y_sol = np.ones((var_disc.shape[0], 1))
+
+        var_casadi = to_casadi(var_disc, y_sol)
+        model = pybamm.BaseModel()
+        model._geometry = geometry
+
+        solution = pybamm.Solution(t_sol, y_sol, model, {})
+
+        processed_var = pybamm.ProcessedVariableUnstructured(
+            "test", [var_disc], [var_casadi], solution
+        )
+
+        x_outside = np.array([2.0])
+        y_outside = np.array([2.0])
+        z_outside = np.array([2.0])
+
+        result_default = processed_var(0, x=x_outside, y=y_outside, z=z_outside)
+        assert np.isnan(result_default)
+
+        result_custom = processed_var(
+            0, x=x_outside, y=y_outside, z=z_outside, fill_value=999
+        )
+        assert result_custom == 999
+
+    def test_processed_variable_unstructured_shape_method(self):
+        from pybamm.meshes.scikit_fem_submeshes_3d import ScikitFemGenerator3D
+
+        x = pybamm.SpatialVariable(
+            "x", domain=["current collector"], coord_sys="cartesian"
+        )
+        y = pybamm.SpatialVariable(
+            "y", domain=["current collector"], coord_sys="cartesian"
+        )
+        z = pybamm.SpatialVariable(
+            "z", domain=["current collector"], coord_sys="cartesian"
+        )
+
+        geometry = {
+            "current collector": {
+                x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+                y: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+                z: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+            }
+        }
+
+        submesh_types = {
+            "current collector": ScikitFemGenerator3D(geom_type="pouch", h=0.3)
+        }
+        var_pts = {x: None, y: None, z: None}
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+
+        spatial_methods = {"current collector": pybamm.ScikitFiniteElement3D()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        var = pybamm.Variable("var", domain=["current collector"])
+        disc.set_variable_slices([var])
+        var_disc = disc.process_symbol(var)
+
+        t_sol = np.linspace(0, 1, 5)
+        y_sol = np.ones((var_disc.shape[0], len(t_sol)))
+
+        var_casadi = to_casadi(var_disc, y_sol)
+        model = pybamm.BaseModel()
+        model._geometry = geometry
+
+        solution = pybamm.Solution(t_sol, y_sol, model, {})
+
+        processed_var = pybamm.ProcessedVariableUnstructured(
+            "test", [var_disc], [var_casadi], solution
+        )
+
+        shape = processed_var._shape(t_sol)
+        assert shape == [var_disc.mesh.npts, len(t_sol)]
+
+    def test_processed_variable_unstructured_time_integral(self):
+        from pybamm.meshes.scikit_fem_submeshes_3d import ScikitFemGenerator3D
+
+        x = pybamm.SpatialVariable(
+            "x", domain=["current collector"], coord_sys="cartesian"
+        )
+        y = pybamm.SpatialVariable(
+            "y", domain=["current collector"], coord_sys="cartesian"
+        )
+        z = pybamm.SpatialVariable(
+            "z", domain=["current collector"], coord_sys="cartesian"
+        )
+
+        geometry = {
+            "current collector": {
+                x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+                y: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+                z: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+            }
+        }
+
+        submesh_types = {
+            "current collector": ScikitFemGenerator3D(geom_type="pouch", h=0.3)
+        }
+        var_pts = {x: None, y: None, z: None}
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+
+        spatial_methods = {"current collector": pybamm.ScikitFiniteElement3D()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        var = pybamm.Variable("var", domain=["current collector"])
+        disc.set_variable_slices([var])
+        var_disc = disc.process_symbol(var)
+
+        t_sol = np.linspace(0, 1, 5)
+        y_sol = np.ones((var_disc.shape[0], len(t_sol)))
+
+        var_casadi = to_casadi(var_disc, y_sol)
+        model = pybamm.BaseModel()
+        model._geometry = geometry
+
+        solution = pybamm.Solution(t_sol, y_sol, model, {})
+
+        processed_var = pybamm.ProcessedVariableUnstructured(
+            "test", [var_disc], [var_casadi], solution
+        )
+
+        assert processed_var.time_integral is None
+
+        time_integral = pybamm.ProcessedVariableTimeIntegral(
+            method="simpson",
+            sum_node=True,
+            initial_condition=np.zeros(var_disc.shape[0]),
+            discrete_times=t_sol,
+        )
+
+        processed_var_with_integral = pybamm.ProcessedVariableUnstructured(
+            "test", [var_disc], [var_casadi], solution, time_integral=time_integral
+        )
+
+        assert processed_var_with_integral.time_integral is time_integral
+
+    def test_process_variable_regular_mesh_fallback(self):
+        var = pybamm.Variable("var", domain=["negative electrode"])
+
+        disc = tests.get_discretisation_for_testing()
+        disc.set_variable_slices([var])
+        var_sol = disc.process_symbol(var)
+
+        t_sol = np.array([0, 1])
+        y_sol = np.ones((var_sol.shape[0], 2))
+
+        var_casadi = to_casadi(var_sol, y_sol)
+
+        model = tests.get_base_model_with_battery_geometry()
+        solution = pybamm.Solution(t_sol, y_sol, model, {})
+
+        processed_var = pybamm.process_variable(
+            "test", [var_sol], [var_casadi], solution
+        )
+
+        assert not isinstance(processed_var, pybamm.ProcessedVariableUnstructured)
+
+    def test_process_variable_unstructured_detection(self):
+        from pybamm.meshes.scikit_fem_submeshes_3d import ScikitFemGenerator3D
+
+        x = pybamm.SpatialVariable(
+            "x", domain=["current collector"], coord_sys="cartesian"
+        )
+        y = pybamm.SpatialVariable(
+            "y", domain=["current collector"], coord_sys="cartesian"
+        )
+        z = pybamm.SpatialVariable(
+            "z", domain=["current collector"], coord_sys="cartesian"
+        )
+
+        geometry = {
+            "current collector": {
+                x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+                y: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+                z: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)},
+            }
+        }
+
+        submesh_types = {
+            "current collector": ScikitFemGenerator3D(geom_type="pouch", h=0.3)
+        }
+        var_pts = {x: None, y: None, z: None}
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+
+        spatial_methods = {"current collector": pybamm.ScikitFiniteElement3D()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        var = pybamm.Variable("var", domain=["current collector"])
+        disc.set_variable_slices([var])
+        var_disc = disc.process_symbol(var)
+
+        t_sol = np.array([0, 1])
+        y_sol = np.ones((var_disc.shape[0], 2))
+
+        var_casadi = to_casadi(var_disc, y_sol)
+        model = pybamm.BaseModel()
+        model._geometry = geometry
+
+        solution = pybamm.Solution(t_sol, y_sol, model, {})
+
+        processed_var = pybamm.process_variable(
+            "test", [var_disc], [var_casadi], solution
+        )
+
+        assert isinstance(processed_var, pybamm.ProcessedVariableUnstructured)
diff --git a/tests/unit/test_spatial_methods/test_scikit_finite_element_3d.py b/tests/unit/test_spatial_methods/test_scikit_finite_element_3d.py
index 60ab6f672..04e150c79 100644
--- a/tests/unit/test_spatial_methods/test_scikit_finite_element_3d.py
+++ b/tests/unit/test_spatial_methods/test_scikit_finite_element_3d.py
@@ -17,8 +17,8 @@ class TestScikitFiniteElement3D:
         with pytest.raises(NotImplementedError):
             spatial_method.indefinite_integral(None, None, None)
 
-    def test_discretise_equations_box(self):
-        mesh = get_3d_mesh_for_testing(geom_type="box", include_particles=False)
+    def test_discretise_equations_pouch(self):
+        mesh = get_3d_mesh_for_testing(geom_type="pouch", include_particles=False)
         spatial_methods = {
             "negative electrode": pybamm.ScikitFiniteElement3D(),
             "separator": pybamm.ScikitFiniteElement3D(),
@@ -86,7 +86,7 @@ class TestScikitFiniteElement3D:
         np.testing.assert_allclose(np.mean(grad_eval[:, 1]), 3, rtol=1e-7, atol=5e-6)
         np.testing.assert_allclose(np.mean(grad_eval[:, 2]), 4, rtol=1e-7, atol=5e-6)
 
-    def test_laplacian_3d_manufactured_box(self):
+    def test_laplacian_3d_manufactured_pouch(self):
         mesh = get_unit_3d_mesh_for_testing(h=0.4)
         disc = pybamm.Discretisation(
             mesh, {"current collector": pybamm.ScikitFiniteElement3D()}
@@ -185,7 +185,7 @@ class TestScikitFiniteElement3D:
     def test_3d_thermal_equation_analytical(self):
         L_x = 2.0
         mesh = get_unit_3d_mesh_for_testing(
-            h=0.4, geom_type="box", x_max=L_x, y_max=1.0, z_max=1.0
+            h=0.4, geom_type="pouch", x_max=L_x, y_max=1.0, z_max=1.0
         )
         disc = pybamm.Discretisation(
             mesh, {"current collector": pybamm.ScikitFiniteElement3D()}
@@ -222,7 +222,7 @@ class TestScikitFiniteElement3D:
 
     def test_gradient_squared_3d(self):
         mesh = get_unit_3d_mesh_for_testing(
-            xpts=6, ypts=6, zpts=6, geom_type="box", include_particles=False
+            xpts=6, ypts=6, zpts=6, geom_type="pouch", include_particles=False
         )
         spatial_methods = {
             "current collector": pybamm.ScikitFiniteElement3D(),
@@ -352,7 +352,7 @@ class TestScikitFiniteElement3D:
 
     def test_definite_integral_3d(self):
         mesh = get_3d_mesh_for_testing(
-            xpts=6, ypts=6, zpts=6, geom_type="box", include_particles=False
+            xpts=6, ypts=6, zpts=6, geom_type="pouch", include_particles=False
         )
         spatial_methods = {
             "negative electrode": pybamm.ScikitFiniteElement3D(),
@@ -380,7 +380,7 @@ class TestScikitFiniteElement3D:
 
     def test_boundary_value_3d(self):
         mesh = get_3d_mesh_for_testing(
-            xpts=6, ypts=6, zpts=6, geom_type="box", include_particles=False
+            xpts=6, ypts=6, zpts=6, geom_type="pouch", include_particles=False
         )
         spatial_methods = {
             "negative electrode": pybamm.ScikitFiniteElement3D(),
@@ -399,7 +399,7 @@ class TestScikitFiniteElement3D:
             np.testing.assert_allclose(result, 1, rtol=1e-6, atol=1e-6)
 
     def test_boundary_integral_3d(self):
-        mesh = get_3d_mesh_for_testing(geom_type="box")
+        mesh = get_3d_mesh_for_testing(geom_type="pouch")
         disc = pybamm.Discretisation(
             mesh, {"negative electrode": pybamm.ScikitFiniteElement3D()}
         )
@@ -428,7 +428,7 @@ class TestScikitFiniteElement3D:
 
     def test_disc_spatial_var_3d(self):
         mesh = get_unit_3d_mesh_for_testing(
-            xpts=4, ypts=4, zpts=4, geom_type="box", include_particles=False
+            xpts=4, ypts=4, zpts=4, geom_type="pouch", include_particles=False
         )
         spatial_methods = {
             "current collector": pybamm.ScikitFiniteElement3D(),
@@ -704,10 +704,10 @@ class TestScikitFiniteElement3D:
         assert nodes[:, 2].min() > 0.0 - 1e-7
         assert nodes[:, 2].max() < 0.8 + 1e-7
 
-    def test_boundary_integral_3d_box(self):
+    def test_boundary_integral_3d_pouch(self):
         x_max, y_max, z_max = 1.0, 2.0, 3.0
         mesh = get_unit_3d_mesh_for_testing(
-            geom_type="box", h=0.3, x_max=x_max, y_max=y_max, z_max=z_max
+            geom_type="pouch", h=0.3, x_max=x_max, y_max=y_max, z_max=z_max
         )
         domain = "current collector"
         disc = pybamm.Discretisation(mesh, {domain: pybamm.ScikitFiniteElement3D()})