Repository: pybamm-team/PyBaMM
Pull request #5009 (2025-05-09 17:38:14)

[GSOC 2025] 3D FEM
# Description

This PR will add a support for 3D FEM, 3D meshing along with a mesh generator with support for prismatic and cylindrical geometry

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


ig-raiderler: The documentation in this repository is excellent! Very clear and comprehensive.
aabills: Hey Rishab, yeah Ghost nodes are a real pain. Let me wrap mine up this week and I'll get back to you ASAP
Rishab87: thanks @aabills, in meantime I'll try to fix these errors myself.
Rishab87: I've removed spiral added heat diffusion examples everything seems to be working now, would love to get some feedback on this
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/5009?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) Report
Attention: Patch coverage is `90.00000%` with `45 lines` in your changes missing coverage. Please review.
> Project coverage is 98.42%. Comparing base [(`c11ee32`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/c11ee32daaf7d9a6bfdeb810a6e77e758c15a91c?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) to head [(`01dbe2a`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/01dbe2ad38063c41cbbd3c238b4f1b8804d26513?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).
> Report is 1 commits behind head on develop.

| [Files with missing lines](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/5009?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) | Patch % | Lines |
|---|---|---|
| [...pybamm/spatial\_methods/scikit\_finite\_element\_3d.py](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/5009?src=pr&el=tree&filepath=src%2Fpybamm%2Fspatial_methods%2Fscikit_finite_element_3d.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team#diff-c3JjL3B5YmFtbS9zcGF0aWFsX21ldGhvZHMvc2Npa2l0X2Zpbml0ZV9lbGVtZW50XzNkLnB5) | 89.08% | [25 Missing :warning: ](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/5009?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) |
| [src/pybamm/expression\_tree/unary\_operators.py](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/5009?src=pr&el=tree&filepath=src%2Fpybamm%2Fexpression_tree%2Funary_operators.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team#diff-c3JjL3B5YmFtbS9leHByZXNzaW9uX3RyZWUvdW5hcnlfb3BlcmF0b3JzLnB5) | 14.28% | [18 Missing :warning: ](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/5009?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) |
| [src/pybamm/meshes/scikit\_fem\_submeshes\_3d.py](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/5009?src=pr&el=tree&filepath=src%2Fpybamm%2Fmeshes%2Fscikit_fem_submeshes_3d.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team#diff-c3JjL3B5YmFtbS9tZXNoZXMvc2Npa2l0X2ZlbV9zdWJtZXNoZXNfM2QucHk=) | 98.92% | [2 Missing :warning: ](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/5009?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) |

<details><summary>Additional details and impacted files</summary>


```diff
@@             Coverage Diff             @@
##           develop    #5009      +/-   ##
===========================================
- Coverage    98.58%   98.42%   -0.16%     
===========================================
  Files          304      306       +2     
  Lines        23725    24172     +447     
===========================================
+ Hits         23389    23792     +403     
- Misses         336      380      +44     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/5009?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).

<details><summary> :rocket: New features to boost your workflow: </summary>

- :snowflake: [Test Analytics](https://docs.codecov.com/docs/test-analytics): Detect flaky tests, report on failures, and find test suite problems.
</details>

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 98e65c7f6..f85ad6c80 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,7 @@
 - Adds `all` key-value pair to `output_variables` sensitivity dictionaries, accessible through `solution[var].sensitivities['all']`. Aligns shape with conventional solution sensitivities object. ([#5067](https://github.com/pybamm-team/PyBaMM/pull/5067))
 - Added a new `BaseHysteresisOpenCircuitPotential` class that sets variables for the lithiation and delithiation OCP and the hysteresis voltage (`H = U_lith - U_delith`). Allow the initial hysteresis state to be a function of position through the electrode. Allow the hysteresis decay rates of the Axen and Wycisk models to be a function of stoichiometry and temperature. Added a heat source term in each active material phase `Q_hys = i_vol * (U - U_eq)` where `i_vol` is the volumetric interfacial current density, `U` is the OCP (i.e. includes hysteresis), and `U_eq` is the "equilibrium OCP". Renamed the open-circuit potential models to be more descriptive. The options "Axen" and "Wycisk" are now "one-state hysteresis" and "one-state differential capacity hysteresis". The old option names still work but will raise a warning. ([#4893](https://github.com/pybamm-team/PyBaMM/pull/4893))
 - Add support for `output_variables` to `pybamm.DiscreteTimeSum` and `pybamm.ExplicitTimeIntegral` expressions. ([#5071](https://github.com/pybamm-team/PyBaMM/pull/5071))
+- Added 3D FEM and meshes supporting rectangular and cylindrical geometries ([#5009](https://github.com/pybamm-team/PyBaMM/pull/5009))
 
 ## Bug fixes
 
diff --git a/docs/source/api/meshes/index.rst b/docs/source/api/meshes/index.rst
index 968b78de0..143e93638 100644
--- a/docs/source/api/meshes/index.rst
+++ b/docs/source/api/meshes/index.rst
@@ -7,3 +7,4 @@ Meshes
   zero_dimensional_submeshes
   one_dimensional_submeshes
   two_dimensional_submeshes
+  three_dimensional_submeshes
diff --git a/docs/source/api/meshes/three_dimensional_submeshes.rst b/docs/source/api/meshes/three_dimensional_submeshes.rst
new file mode 100644
index 000000000..575cf33d9
--- /dev/null
+++ b/docs/source/api/meshes/three_dimensional_submeshes.rst
@@ -0,0 +1,8 @@
+3D Sub Meshes
+==================
+
+.. autoclass:: pybamm.ScikitFemSubMesh3D
+    :members:
+
+.. autoclass:: pybamm.ScikitFemGenerator3D
+    :members:
diff --git a/docs/source/api/spatial_methods/index.rst b/docs/source/api/spatial_methods/index.rst
index 060eb91c0..f9ccacd0d 100644
--- a/docs/source/api/spatial_methods/index.rst
+++ b/docs/source/api/spatial_methods/index.rst
@@ -9,3 +9,4 @@ Discretisation and spatial methods
   spectral_volume
   scikit_finite_element
   zero_dimensional_method
+  scikit_finite_element_3d
diff --git a/docs/source/api/spatial_methods/scikit_finite_element_3d.rst b/docs/source/api/spatial_methods/scikit_finite_element_3d.rst
new file mode 100644
index 000000000..d6e4ff3db
--- /dev/null
+++ b/docs/source/api/spatial_methods/scikit_finite_element_3d.rst
@@ -0,0 +1,5 @@
+Scikit Finite Elements 3D
+=========================
+
+.. autoclass:: pybamm.ScikitFiniteElement3D
+    :members:
diff --git a/examples/scripts/3d_examples/3d_eqation_cylindric_non_zero_flux.py b/examples/scripts/3d_examples/3d_eqation_cylindric_non_zero_flux.py
new file mode 100644
index 000000000..4d5f1f406
--- /dev/null
+++ b/examples/scripts/3d_examples/3d_eqation_cylindric_non_zero_flux.py
@@ -0,0 +1,210 @@
+import matplotlib.pyplot as plt
+import matplotlib.tri as tri
+import numpy as np
+from scipy.interpolate import griddata
+
+import pybamm
+from pybamm import (
+    BaseModel,
+    Discretisation,
+    Scalar,
+    ScikitFiniteElement3D,
+)
+from pybamm.meshes.scikit_fem_submeshes_3d import ScikitFemGenerator3D
+
+R_outer = 0.6
+R_inner = 0.1
+H = 1.0
+alpha = 0.1
+
+heat_generation = 1000.0
+heat_flux_top = 500.0
+heat_flux_bottom = -300.0
+
+model = BaseModel()
+
+r = pybamm.SpatialVariable("r", ["current collector"], coord_sys="cylindrical polar")
+theta = pybamm.SpatialVariable(
+    "theta", ["current collector"], coord_sys="cylindrical polar"
+)
+z = pybamm.SpatialVariable("z", ["current collector"], coord_sys="cylindrical polar")
+T = pybamm.Variable("T", domain="current collector")
+
+source_term = pybamm.source(heat_generation, T)
+
+model.algebraic = {T: alpha * pybamm.laplacian(T) + source_term}
+model.initial_conditions = {T: Scalar(25)}
+
+model.boundary_conditions = {
+    T: {
+        "r_min": (Scalar(100), "Dirichlet"),
+        "r_max": (Scalar(20), "Dirichlet"),
+        "z_min": (Scalar(heat_flux_bottom), "Neumann"),
+        "z_max": (Scalar(heat_flux_top), "Neumann"),
+    }
+}
+
+geometry = {
+    "current collector": {
+        r: {"min": pybamm.Scalar(R_inner), "max": pybamm.Scalar(R_outer)},
+        z: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(H)},
+    }
+}
+
+submesh_types = {
+    "current collector": ScikitFemGenerator3D(geom_type="cylinder", h=0.08)
+}
+
+var_pts = {r: None, theta: None, z: None}
+mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+spatial_methods = {"current collector": ScikitFiniteElement3D()}
+disc = Discretisation(mesh, spatial_methods)
+disc.process_model(model)
+
+solver = pybamm.AlgebraicSolver()
+solution = solver.solve(model)
+
+nodes = mesh["current collector"].nodes
+T_solution = solution.y.flatten()
+
+print(f"Number of FEM nodes: {nodes.shape[0]}")
+print(f"Temperature range: [{T_solution.min():.2f}, {T_solution.max():.2f}]°C")
+
+Nr_vis, Ntheta_vis, Nz_vis = 30, 40, 25
+r_grid = np.linspace(R_inner, R_outer, Nr_vis)
+theta_grid = np.linspace(0, 2 * np.pi, Ntheta_vis)
+z_grid = np.linspace(0, H, Nz_vis)
+R_grid, Theta_grid, Z_grid = np.meshgrid(r_grid, theta_grid, z_grid, indexing="ij")
+
+X_grid = R_grid * np.cos(Theta_grid)
+Y_grid = R_grid * np.sin(Theta_grid)
+grid_points = np.column_stack([X_grid.ravel(), Y_grid.ravel(), Z_grid.ravel()])
+
+T_interpolated = griddata(
+    nodes, T_solution, grid_points, method="linear", fill_value=np.nan
+).reshape((Nr_vis, Ntheta_vis, Nz_vis))
+
+fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
+fig.suptitle("Heat Distribution with Source Term and Neumann BCs", fontsize=16)
+
+mid_k = Nz_vis // 2
+R_plane = R_grid[:, :, mid_k]
+Theta_plane = Theta_grid[:, :, mid_k]
+T_plane = T_interpolated[:, :, mid_k]
+ax1.remove()
+ax1 = fig.add_subplot(2, 2, 1, projection="polar")
+pcm1 = ax1.contourf(Theta_plane, R_plane, T_plane, levels=20, cmap="plasma")
+ax1.set_title(f"T(r,θ) at z={H / 2:.1f}m")
+ax1.set_ylim(R_inner, R_outer)
+fig.colorbar(pcm1, ax=ax1)
+
+mid_j = 0
+R_axial = R_grid[:, mid_j, :]
+Z_axial = Z_grid[:, mid_j, :]
+T_axial = T_interpolated[:, mid_j, :]
+pcm2 = ax2.contourf(R_axial, Z_axial, T_axial, levels=20, cmap="plasma")
+ax2.set_xlabel("r [m]")
+ax2.set_ylabel("z [m]")
+ax2.set_title("T(r,z) - Shows effect of Neumann BCs")
+ax2.set_aspect("equal")
+fig.colorbar(pcm2, ax=ax2)
+
+unique_z_coords = np.unique(nodes[:, 2])
+mid_plane_z = unique_z_coords[np.argmin(np.abs(unique_z_coords - H / 2))]
+
+z_mid_indices = nodes[:, 2] == mid_plane_z
+z_mid_nodes = nodes[z_mid_indices]
+T_mid = T_solution[z_mid_indices]
+
+triang = tri.Triangulation(z_mid_nodes[:, 0], z_mid_nodes[:, 1])
+pcm3 = ax3.tricontourf(triang, T_mid, levels=20, cmap="plasma")
+ax3.plot(
+    R_outer * np.cos(np.linspace(0, 2 * np.pi, 100)),
+    R_outer * np.sin(np.linspace(0, 2 * np.pi, 100)),
+    "k-",
+    linewidth=2,
+)
+ax3.plot(
+    R_inner * np.cos(np.linspace(0, 2 * np.pi, 100)),
+    R_inner * np.sin(np.linspace(0, 2 * np.pi, 100)),
+    "w-",
+    linewidth=2,
+)
+ax3.set_xlabel("x [m]")
+ax3.set_ylabel("y [m]")
+ax3.set_title("FEM Solution at Mid-Height")
+ax3.set_aspect("equal")
+fig.colorbar(pcm3, ax=ax3, label="Temperature [°C]")
+
+z_line = np.linspace(0, H, 100)
+mid_radius = (R_inner + R_outer) / 2
+axial_points = np.column_stack(
+    [np.full_like(z_line, mid_radius), np.zeros_like(z_line), z_line]
+)
+T_axial_profile = griddata(nodes, T_solution, axial_points, method="linear")
+ax4.plot(z_line, T_axial_profile, "b-", linewidth=3, label="FEM solution")
+ax4.axhline(100, color="r", linestyle="--", alpha=0.7, label="Inner BC (100°C)")
+ax4.axhline(20, color="g", linestyle="--", alpha=0.7, label="Outer BC (20°C)")
+ax4.set_xlabel("z [m]")
+ax4.set_ylabel("Temperature [°C]")
+ax4.set_title("Axial Temperature Profile\n(shows Neumann BC effects)")
+ax4.legend()
+ax4.grid(True, alpha=0.3)
+
+ax4.text(
+    0.05,
+    T_axial_profile[5],
+    f"q = {heat_flux_bottom} W/m²\n(cooling)",
+    bbox=dict(boxstyle="round,pad=0.3", facecolor="lightblue"),
+)
+ax4.text(
+    0.85,
+    T_axial_profile[-5],
+    f"q = {heat_flux_top} W/m²\n(heating)",
+    bbox=dict(boxstyle="round,pad=0.3", facecolor="lightcoral"),
+)
+
+plt.tight_layout(rect=[0, 0.03, 1, 0.95])
+plt.show()
+
+
+def find_closest_node(target_point, nodes):
+    distances = np.linalg.norm(nodes - target_point, axis=1)
+    return np.argmin(distances)
+
+
+inner_point = np.array([R_inner, 0, H / 2])
+outer_point = np.array([R_outer, 0, H / 2])
+mid_point = np.array([(R_inner + R_outer) / 2, 0, H / 2])
+bottom_point = np.array([(R_inner + R_outer) / 2, 0, 0])
+top_point = np.array([(R_inner + R_outer) / 2, 0, H])
+
+inner_idx = find_closest_node(inner_point, nodes)
+outer_idx = find_closest_node(outer_point, nodes)
+mid_idx = find_closest_node(mid_point, nodes)
+bottom_idx = find_closest_node(bottom_point, nodes)
+top_idx = find_closest_node(top_point, nodes)
+
+print("\n" + "=" * 60)
+print("HEAT DISTRIBUTION WITH SOURCE TERM AND NEUMANN BCs")
+print("=" * 60)
+print("Boundary Conditions:")
+print(f"  Inner radius (r={R_inner}m): T = 100°C (Dirichlet)")
+print(f"  Outer radius (r={R_outer}m): T = 20°C (Dirichlet)")
+print(f"  Bottom surface (z=0): q = {heat_flux_bottom} W/m² (Neumann - cooling)")
+print(f"  Top surface (z={H}m): q = {heat_flux_top} W/m² (Neumann - heating)")
+print(f"  Volumetric heat source: Q = {heat_generation} W/m³")
+print("\nTemperature Results:")
+print(f"  Inner surface: {T_solution[inner_idx]:.1f}°C")
+print(f"  Outer surface: {T_solution[outer_idx]:.1f}°C")
+print(f"  Mid-radius center: {T_solution[mid_idx]:.1f}°C")
+print(f"  Bottom center: {T_solution[bottom_idx]:.1f}°C")
+print(f"  Top center: {T_solution[top_idx]:.1f}°C")
+print(f"  Overall range: [{T_solution.min():.1f}, {T_solution.max():.1f}]°C")
+
+temp_rise_axial = T_solution[top_idx] - T_solution[bottom_idx]
+print(f"\nAxial temperature rise due to Neumann BCs: {temp_rise_axial:.1f}°C")
+print("This demonstrates how:")
+print("  1. Positive Neumann flux (top) heats the surface")
+print("  2. Negative Neumann flux (bottom) cools the surface")
+print("  3. Source term raises overall temperature level")
diff --git a/examples/scripts/3d_examples/3d_equation_cylindric.py b/examples/scripts/3d_examples/3d_equation_cylindric.py
new file mode 100644
index 000000000..486a31c17
--- /dev/null
+++ b/examples/scripts/3d_examples/3d_equation_cylindric.py
@@ -0,0 +1,201 @@
+import matplotlib.pyplot as plt
+import matplotlib.tri as tri
+import numpy as np
+from scipy.interpolate import griddata
+
+import pybamm
+from pybamm import (
+    BaseModel,
+    Discretisation,
+    Scalar,
+    ScikitFiniteElement3D,
+)
+from pybamm.meshes.scikit_fem_submeshes_3d import ScikitFemGenerator3D
+
+R_outer = 0.6
+R_inner = 0.1
+H = 1.0
+alpha = 0.1
+
+model = BaseModel()
+
+r = pybamm.SpatialVariable("r", ["current collector"], coord_sys="cylindrical polar")
+theta = pybamm.SpatialVariable(
+    "theta", ["current collector"], coord_sys="cylindrical polar"
+)
+z = pybamm.SpatialVariable("z", ["current collector"], coord_sys="cylindrical polar")
+T = pybamm.Variable("T", domain="current collector")
+
+model.algebraic = {T: alpha * pybamm.laplacian(T)}
+model.initial_conditions = {T: Scalar(25)}
+model.boundary_conditions = {
+    T: {
+        "r_min": (Scalar(100), "Dirichlet"),
+        "r_max": (Scalar(0), "Dirichlet"),
+        "z_min": (Scalar(0), "Neumann"),
+        "z_max": (Scalar(0), "Neumann"),
+    }
+}
+
+geometry = {
+    "current collector": {
+        r: {"min": pybamm.Scalar(R_inner), "max": pybamm.Scalar(R_outer)},
+        z: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(H)},
+    }
+}
+
+# For unstructured 3D FEM meshes, we control the density and quality using a
+# characteristic length parameter 'h'. The generator will try to create
+# tetrahedral elements with edge lengths around this value.
+submesh_types = {
+    "current collector": ScikitFemGenerator3D(geom_type="cylinder", h=0.08)
+}
+
+# The 'var_pts' dictionary is required by the pybamm.Mesh class for validation,
+# but our ScikitFemGenerator3D completely IGNORES it. The number of points
+# is determined by the meshing algorithm based on 'h', not by 'var_pts'.
+# We pass placeholder values to satisfy the class constructor.
+# We may want to pass var_pts anyway since other domains may not be in 3D in a simulation
+var_pts = {r: None, theta: None, z: None}
+mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+spatial_methods = {"current collector": ScikitFiniteElement3D()}
+disc = Discretisation(mesh, spatial_methods)
+disc.process_model(model)
+
+solver = pybamm.AlgebraicSolver()
+solution = solver.solve(model)
+
+nodes = mesh["current collector"].nodes
+T_solution = solution.y.flatten()
+
+print(f"Number of FEM nodes: {nodes.shape[0]}")
+print(f"Temperature range: [{T_solution.min():.2f}, {T_solution.max():.2f}]°C")
+
+# Create a regular grid for visualization purposes.
+Nr_vis, Ntheta_vis, Nz_vis = 30, 40, 25
+r_grid = np.linspace(R_inner, R_outer, Nr_vis)
+theta_grid = np.linspace(0, 2 * np.pi, Ntheta_vis)
+z_grid = np.linspace(0, H, Nz_vis)
+R_grid, Theta_grid, Z_grid = np.meshgrid(r_grid, theta_grid, z_grid, indexing="ij")
+
+X_grid = R_grid * np.cos(Theta_grid)
+Y_grid = R_grid * np.sin(Theta_grid)
+grid_points = np.column_stack([X_grid.ravel(), Y_grid.ravel(), Z_grid.ravel()])
+
+# Interpolate the unstructured FEM solution onto the regular visualization grid.
+# This can sometimes make the plot look "lumpy" if the mesh is coarse.
+T_interpolated = griddata(
+    nodes, T_solution, grid_points, method="linear", fill_value=np.nan
+).reshape((Nr_vis, Ntheta_vis, Nz_vis))
+
+# Create the 4-panel plot
+fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
+fig.suptitle("Steady-State Heat Distribution Analysis", fontsize=16)
+
+# Plot 1: Polar plot (r-theta plane) at mid-height
+mid_k = Nz_vis // 2
+R_plane = R_grid[:, :, mid_k]
+Theta_plane = Theta_grid[:, :, mid_k]
+T_plane = T_interpolated[:, :, mid_k]
+ax1.remove()
+ax1 = fig.add_subplot(2, 2, 1, projection="polar")
+pcm1 = ax1.contourf(Theta_plane, R_plane, T_plane, levels=20, cmap="viridis")
+ax1.set_title(f"T(r,θ) at z={H / 2:.1f}m (Interpolated)")
+ax1.set_ylim(R_inner, R_outer)
+fig.colorbar(pcm1, ax=ax1)
+
+# Plot 2: Cross-section (r-z plane) at theta=0
+mid_j = 0
+R_axial = R_grid[:, mid_j, :]
+Z_axial = Z_grid[:, mid_j, :]
+T_axial = T_interpolated[:, mid_j, :]
+pcm2 = ax2.contourf(R_axial, Z_axial, T_axial, levels=20, cmap="viridis")
+ax2.set_xlabel("r [m]")
+ax2.set_ylabel("z [m]")
+ax2.set_title("T(r,z) at θ=0 (Interpolated)")
+ax2.set_aspect("equal")
+fig.colorbar(pcm2, ax=ax2)
+
+# Plot 3: Direct FEM solution plot at mid-height (no interpolation)
+unique_z_coords = np.unique(nodes[:, 2])
+mid_plane_z = unique_z_coords[np.argmin(np.abs(unique_z_coords - H / 2))]
+print(
+    f"\nPlotting direct FEM solution on the mesh plane closest to mid-height (z={mid_plane_z:.3f}m)."
+)
+
+z_mid_indices = nodes[:, 2] == mid_plane_z
+z_mid_nodes = nodes[z_mid_indices]
+T_mid = T_solution[z_mid_indices]
+
+triang = tri.Triangulation(z_mid_nodes[:, 0], z_mid_nodes[:, 1])
+pcm3 = ax3.tricontourf(triang, T_mid, levels=20, cmap="viridis")
+ax3.plot(
+    R_outer * np.cos(np.linspace(0, 2 * np.pi, 100)),
+    R_outer * np.sin(np.linspace(0, 2 * np.pi, 100)),
+    "k-",
+)
+ax3.plot(
+    R_inner * np.cos(np.linspace(0, 2 * np.pi, 100)),
+    R_inner * np.sin(np.linspace(0, 2 * np.pi, 100)),
+    "r-",
+)
+ax3.set_xlabel("x [m]")
+ax3.set_ylabel("y [m]")
+ax3.set_title("FEM Solution at Mid-Height (Direct)")
+ax3.set_aspect("equal")
+fig.colorbar(pcm3, ax=ax3, label="Temperature [°C]")
+
+# Plot 4: Radial temperature profile vs. analytical solution
+r_line = np.linspace(R_inner, R_outer, 100)
+radial_points = np.column_stack(
+    [r_line, np.zeros_like(r_line), np.full_like(r_line, H / 2)]
+)
+T_radial = griddata(nodes, T_solution, radial_points, method="linear")
+ax4.plot(r_line, T_radial, "b-", linewidth=3, label="FEM solution")
+A = -100 / np.log(R_outer / R_inner)
+B = 100 - A * np.log(R_inner)
+T_analytical = A * np.log(r_line) + B
+ax4.plot(r_line, T_analytical, "r--", linewidth=2, label="Analytical solution")
+ax4.set_xlabel("r [m]")
+ax4.set_ylabel("Temperature [°C]")
+ax4.set_title("Radial Temperature Profile")
+ax4.legend()
+ax4.grid(True, alpha=0.5)
+
+plt.tight_layout(rect=[0, 0.03, 1, 0.95])
+plt.show()
+
+# Mesh Visualisation
+print("\nGenerating mesh visualization...")
+fig_mesh, ax_mesh = plt.subplots(figsize=(8, 8))
+ax_mesh.triplot(triang, "k-", lw=0.5, alpha=0.7)
+ax_mesh.plot(z_mid_nodes[:, 0], z_mid_nodes[:, 1], "o", markersize=2, color="blue")
+ax_mesh.set_title(f"FEM Mesh Visualization (Slice at z={mid_plane_z:.3f}m)")
+ax_mesh.set_xlabel("x [m]")
+ax_mesh.set_ylabel("y [m]")
+ax_mesh.set_aspect("equal")
+plt.show()
+
+
+def find_closest_node(target_point, nodes):
+    distances = np.linalg.norm(nodes - target_point, axis=1)
+    return np.argmin(distances)
+
+
+inner_point = np.array([R_inner, 0, H / 2])
+outer_point = np.array([R_outer, 0, H / 2])
+mid_point = np.array([(R_inner + R_outer) / 2, 0, H / 2])
+inner_idx = find_closest_node(inner_point, nodes)
+outer_idx = find_closest_node(outer_point, nodes)
+mid_idx = find_closest_node(mid_point, nodes)
+
+print("\n" + "=" * 50)
+print("STEADY-STATE HEAT DISTRIBUTION ANALYSIS")
+print("=" * 50)
+print(f"FEM Inner surface temperature: {T_solution[inner_idx]:.2f}°C")
+print(f"FEM Outer surface temperature: {T_solution[outer_idx]:.2f}°C")
+print(f"FEM Mid-radius temperature: {T_solution[mid_idx]:.2f}°C")
+T_mid_analytical = A * np.log((R_inner + R_outer) / 2) + B
+print(f"Analytical mid-radius temperature: {T_mid_analytical:.2f}°C")
+error = abs(T_solution[mid_idx] - T_mid_analytical)
+print(f"Error at mid-radius: {error:.3f}°C")
diff --git a/examples/scripts/3d_examples/3d_equation_prismatic.py b/examples/scripts/3d_examples/3d_equation_prismatic.py
new file mode 100644
index 000000000..ade089f13
--- /dev/null
+++ b/examples/scripts/3d_examples/3d_equation_prismatic.py
@@ -0,0 +1,176 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.interpolate import griddata
+
+import pybamm
+from pybamm import (
+    BaseModel,
+    Discretisation,
+    Scalar,
+    ScikitFiniteElement3D,
+    SpatialVariable,
+)
+
+Lx, Ly, Lz = 1.0, 0.8, 0.6
+alpha = 0.1
+t_max = 0.2
+
+model = BaseModel()
+
+x = SpatialVariable("x", ["current collector"], coord_sys="cartesian")
+y = SpatialVariable("y", ["current collector"], coord_sys="cartesian")
+z = SpatialVariable("z", ["current collector"], coord_sys="cartesian")
+
+T = pybamm.Variable("T", domain="current collector")
+model.variables = {"T": T}
+
+model.algebraic = {T: alpha * pybamm.laplacian(T) - pybamm.source(0, T)}
+
+model.initial_conditions = {T: Scalar(0)}
+
+model.boundary_conditions = {
+    T: {
+        "x_min": (Scalar(100), "Dirichlet"),
+        "x_max": (Scalar(0), "Dirichlet"),
+        "y_min": (Scalar(0), "Neumann"),
+        "y_max": (Scalar(0), "Neumann"),
+        "z_min": (Scalar(0), "Neumann"),
+        "z_max": (Scalar(0), "Neumann"),
+    }
+}
+
+geometry = {
+    "current collector": {
+        x: {"min": Scalar(0), "max": Scalar(Lx)},
+        y: {"min": Scalar(0), "max": Scalar(Ly)},
+        z: {"min": Scalar(0), "max": Scalar(Lz)},
+    }
+}
+
+submesh_types = {
+    "current collector": pybamm.ScikitFemGenerator3D(geom_type="box", h=0.15)
+}
+var_pts = {x: None, y: None, z: None}
+mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+
+spatial_methods = {"current collector": ScikitFiniteElement3D()}
+disc = Discretisation(mesh, spatial_methods)
+disc.process_model(model)
+
+solver = pybamm.AlgebraicSolver()
+solution = solver.solve(model)
+
+nodes = mesh["current collector"].nodes
+T_solution = solution.y.flatten()
+
+print(f"Number of FEM nodes: {nodes.shape[0]}")
+print(f"Temperature range: [{T_solution.min():.2f}, {T_solution.max():.2f}]°C")
+
+Nx_vis, Ny_vis, Nz_vis = 25, 20, 15
+x_grid = np.linspace(0, Lx, Nx_vis)
+y_grid = np.linspace(0, Ly, Ny_vis)
+z_grid = np.linspace(0, Lz, Nz_vis)
+X, Y, Z = np.meshgrid(x_grid, y_grid, z_grid, indexing="ij")
+grid_points = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()])
+T_interp_3d = griddata(
+    nodes, T_solution, grid_points, method="linear", fill_value=0
+).reshape((Nx_vis, Ny_vis, Nz_vis))
+
+fig, axes = plt.subplots(1, 3, figsize=(18, 5))
+mid_y = Ny_vis // 2
+mid_z = Nz_vis // 2
+mid_x = Nx_vis // 2
+
+im0 = axes[0].contourf(
+    X[:, mid_y, :], Z[:, mid_y, :], T_interp_3d[:, mid_y, :], levels=20, cmap="inferno"
+)
+axes[0].set_xlabel("x [m]")
+axes[0].set_ylabel("z [m]")
+axes[0].set_title(f"T(x,z) at y={y_grid[mid_y]:.2f}")
+plt.colorbar(im0, ax=axes[0])
+
+im1 = axes[1].contourf(
+    X[:, :, mid_z], Y[:, :, mid_z], T_interp_3d[:, :, mid_z], levels=20, cmap="inferno"
+)
+axes[1].set_xlabel("x [m]")
+axes[1].set_ylabel("y [m]")
+axes[1].set_title(f"T(x,y) at z={z_grid[mid_z]:.2f}")
+plt.colorbar(im1, ax=axes[1])
+
+im2 = axes[2].contourf(
+    Y[mid_x, :, :], Z[mid_x, :, :], T_interp_3d[mid_x, :, :], levels=20, cmap="inferno"
+)
+axes[2].set_xlabel("y [m]")
+axes[2].set_ylabel("z [m]")
+axes[2].set_title(f"T(y,z) at x={x_grid[mid_x]:.2f}")
+plt.colorbar(im2, ax=axes[2])
+
+plt.tight_layout()
+plt.show()
+
+fig = plt.figure(figsize=(8, 6))
+ax = fig.add_subplot(111, projection="3d")
+surf = ax.plot_surface(
+    X[:, mid_y, :],
+    Z[:, mid_y, :],
+    T_interp_3d[:, mid_y, :],
+    cmap="inferno",
+    edgecolor="none",
+    alpha=0.9,
+)
+ax.set_xlabel("x [m]")
+ax.set_ylabel("z [m]")
+ax.set_zlabel("T [°C]")
+ax.set_title(f"T(x,z) at y={y_grid[mid_y]:.2f} (3D Surface)")
+fig.colorbar(surf, ax=ax, shrink=0.5, label="Temperature [°C]")
+plt.tight_layout()
+plt.show()
+
+fig = plt.figure(figsize=(10, 8))
+ax = fig.add_subplot(111, projection="3d")
+scatter = ax.scatter(
+    nodes[:, 0], nodes[:, 1], nodes[:, 2], c=T_solution, cmap="inferno", s=15, alpha=0.6
+)
+ax.set_xlabel("x [m]")
+ax.set_ylabel("y [m]")
+ax.set_zlabel("z [m]")
+ax.set_title("3D Temperature Distribution (FEM Nodes)")
+fig.colorbar(scatter, ax=ax, shrink=0.5, label="Temperature [°C]")
+plt.tight_layout()
+plt.show()
+
+
+def find_closest_node(target_point, nodes):
+    distances = np.linalg.norm(nodes - target_point, axis=1)
+    return np.argmin(distances)
+
+
+center_point = np.array([Lx / 2, Ly / 2, Lz / 2])
+left_edge = np.array([0.1, Ly / 2, Lz / 2])
+right_edge = np.array([Lx - 0.1, Ly / 2, Lz / 2])
+
+center_idx = find_closest_node(center_point, nodes)
+left_idx = find_closest_node(left_edge, nodes)
+right_idx = find_closest_node(right_edge, nodes)
+
+x_line_points = np.column_stack(
+    [np.linspace(0, Lx, 50), np.full(50, Ly / 2), np.full(50, Lz / 2)]
+)
+T_line = griddata(nodes, T_solution, x_line_points, method="linear")
+
+plt.figure(figsize=(8, 5))
+plt.plot(x_line_points[:, 0], T_line, "b-", linewidth=2, label="FEM solution")
+plt.axhline(100, color="r", linestyle="--", alpha=0.7, label="Left BC (100°C)")
+plt.axhline(0, color="g", linestyle="--", alpha=0.7, label="Right BC (0°C)")
+plt.xlabel("x [m]")
+plt.ylabel("Temperature [°C]")
+plt.title("Temperature Profile Along x-axis (centerline)")
+plt.legend()
+plt.grid(True, alpha=0.3)
+plt.show()
+
+print("\nSTEADY-STATE HEAT DISTRIBUTION")
+print(f"Center temperature: {T_solution[center_idx]:.1f}°C")
+print(f"Near left (hot): {T_solution[left_idx]:.1f}°C")
+print(f"Near right (cold): {T_solution[right_idx]:.1f}°C")
+print("Heat flows from left (100°C) to right (0°C)")
diff --git a/src/pybamm/__init__.py b/src/pybamm/__init__.py
index cccfe1538..b33712a96 100644
--- a/src/pybamm/__init__.py
+++ b/src/pybamm/__init__.py
@@ -134,6 +134,7 @@ from .meshes.one_dimensional_submeshes import (
     SpectralVolume1DSubMesh,
     SymbolicUniform1DSubMesh,
 )
+
 from .meshes.scikit_fem_submeshes import (
     ScikitSubMesh2D,
     ScikitUniform2DSubMesh,
@@ -142,6 +143,11 @@ from .meshes.scikit_fem_submeshes import (
     UserSupplied2DSubMesh,
 )
 
+from .meshes.scikit_fem_submeshes_3d import (
+    ScikitFemSubMesh3D,
+    ScikitFemGenerator3D,
+)
+
 # Serialisation
 from .models.base_model import load_model
 
@@ -151,6 +157,7 @@ from .spatial_methods.zero_dimensional_method import ZeroDimensionalSpatialMetho
 from .spatial_methods.finite_volume import FiniteVolume
 from .spatial_methods.spectral_volume import SpectralVolume
 from .spatial_methods.scikit_finite_element import ScikitFiniteElement
+from .spatial_methods.scikit_finite_element_3d import ScikitFiniteElement3D
 
 # Solver classes
 from .solvers.solution import Solution, EmptySolution, make_cycle_solution
diff --git a/src/pybamm/expression_tree/__init__.py b/src/pybamm/expression_tree/__init__.py
index 7ac80e535..84f74626b 100644
--- a/src/pybamm/expression_tree/__init__.py
+++ b/src/pybamm/expression_tree/__init__.py
@@ -2,4 +2,4 @@ __all__ = ['array', 'averages', 'binary_operators', 'broadcasts',
            'concatenations', 'exceptions', 'functions', 'independent_variable',
            'input_parameter', 'interpolant', 'matrix', 'operations',
            'parameter', 'printing', 'scalar', 'state_vector', 'symbol',
-           'unary_operators', 'variable', 'vector', 'discrete_time_sum' ]
+           'unary_operators', 'variable', 'vector', 'discrete_time_sum']
diff --git a/src/pybamm/expression_tree/unary_operators.py b/src/pybamm/expression_tree/unary_operators.py
index 8a9f21359..828bc4bdf 100644
--- a/src/pybamm/expression_tree/unary_operators.py
+++ b/src/pybamm/expression_tree/unary_operators.py
@@ -905,6 +905,22 @@ class BoundaryIntegral(SpatialOperator):
             name += "negative tab"
         elif region == "positive tab":
             name += "positive tab"
+        elif region == "x_min":
+            name += "x_min"
+        elif region == "x_max":
+            name += "x_max"
+        elif region == "y_min":
+            name += "y_min"
+        elif region == "y_max":
+            name += "y_max"
+        elif region == "z_min":
+            name += "z_min"
+        elif region == "z_max":
+            name += "z_max"
+        elif region == "r_min":  # pragma: no cover
+            name += "r_min"
+        elif region == "r_max":  # pragma: no cover
+            name += "r_max"
         self.region = region
         super().__init__(name, child, domains)
 
diff --git a/src/pybamm/meshes/__init__.py b/src/pybamm/meshes/__init__.py
index 0bd3f8d78..a539e40e4 100644
--- a/src/pybamm/meshes/__init__.py
+++ b/src/pybamm/meshes/__init__.py
@@ -1,2 +1,2 @@
 __all__ = ['meshes', 'one_dimensional_submeshes', 'scikit_fem_submeshes',
-           'zero_dimensional_submesh']
+           'zero_dimensional_submesh', 'scikit_fem_submeshes_3d']
diff --git a/src/pybamm/meshes/meshes.py b/src/pybamm/meshes/meshes.py
index 8a296b7f3..e90afd56c 100644
--- a/src/pybamm/meshes/meshes.py
+++ b/src/pybamm/meshes/meshes.py
@@ -41,6 +41,21 @@ class Mesh(dict):
         # convert var_pts to an id dict
         var_name_pts = {var.name: pts for var, pts in var_pts.items()}
 
+        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 provided_pts is not None:
+                        pybamm.logger.warning(
+                            f"For the 3D FEM submesh on domain '{domain}', the value "
+                            f"of all 'var_pts' are ignored. "
+                            "The mesh is generated to satisfy the 'h' parameter. "
+                            "It is recommended to pass 'None' for spatial variables "
+                            "on a 3D FEM domain to avoid this warning."
+                        )
+                        # We only need to warn once per domain.
+                        break
+
         # create submesh_pts from var_pts
         submesh_pts = {}
         for domain in geometry:
@@ -231,7 +246,12 @@ class Mesh(dict):
             for domain, submesh in self.items()
             if (
                 len(domain) == 1
-                and not isinstance(submesh, pybamm.SubMesh0D | pybamm.ScikitSubMesh2D)
+                and not isinstance(
+                    submesh,
+                    pybamm.SubMesh0D
+                    | pybamm.ScikitSubMesh2D
+                    | pybamm.ScikitFemSubMesh3D,
+                )
             )
         ]
         for domain, submesh in submeshes:
diff --git a/src/pybamm/meshes/scikit_fem_submeshes_3d.py b/src/pybamm/meshes/scikit_fem_submeshes_3d.py
new file mode 100644
index 000000000..22f865dae
--- /dev/null
+++ b/src/pybamm/meshes/scikit_fem_submeshes_3d.py
@@ -0,0 +1,319 @@
+import numpy as np
+
+import pybamm
+from pybamm.util import import_optional_dependency
+
+
+class ScikitFemGenerator3D(pybamm.MeshGenerator):
+    """
+    A mesh generator that creates 3D tetrahedral meshes using scikit-fem.
+
+    Parameters
+    ----------
+    geom_type : str
+        The type of geometry to generate. Must be one of "box" for a rectangular
+        prism, or "cylinder" for a cylindrical annulus.
+    gen_params : dict
+        A dictionary of geometry-specific parameters. for eg:
+
+        - "h" : float, optional
+            The target characteristic length of the mesh elements. A smaller 'h'
+            results in a finer mesh and a more accurate solution, but
+            increases computation time. Default is 0.3.
+    """
+
+    def __init__(self, geom_type, **gen_params):
+        supported_geometries = ["box", "cylinder"]
+        if geom_type not in supported_geometries:
+            raise pybamm.GeometryError(
+                f"geom_type must be one of {supported_geometries}, not '{geom_type}'"
+            )
+        super().__init__(ScikitFemSubMesh3D)
+        self.geom_type = geom_type
+        self.gen_params = gen_params
+
+    def _make_box_mesh(self, x_lim, y_lim, z_lim, h):
+        """
+        Create a structured box mesh.
+
+        Parameters
+        ----------
+        x_lim, y_lim, z_lim : tuple
+            Domain limits for each dimension
+        h : float
+            Target mesh size
+
+        Returns
+        -------
+        skfem.MeshTet
+            Box mesh with proper boundary tags
+        """
+        skfem = import_optional_dependency("skfem")
+        x_min, x_max = x_lim
+        y_min, y_max = y_lim
+        z_min, z_max = z_lim
+
+        nx = max(2, int(np.ceil((x_max - x_min) / h)) + 1)
+        ny = max(2, int(np.ceil((y_max - y_min) / h)) + 1)
+        nz = max(2, int(np.ceil((z_max - z_min) / h)) + 1)
+
+        mesh = skfem.MeshTet.init_tensor(
+            np.linspace(x_min, x_max, nx),
+            np.linspace(y_min, y_max, ny),
+            np.linspace(z_min, z_max, nz),
+        )
+
+        bnd_facets = mesh.boundary_facets()
+        midpoints = mesh.p[:, mesh.facets[:, bnd_facets]].mean(axis=1)
+
+        boundaries = {
+            "x_min": bnd_facets[np.isclose(midpoints[0], x_min)],
+            "x_max": bnd_facets[np.isclose(midpoints[0], x_max)],
+            "y_min": bnd_facets[np.isclose(midpoints[1], y_min)],
+            "y_max": bnd_facets[np.isclose(midpoints[1], y_max)],
+            "z_min": bnd_facets[np.isclose(midpoints[2], z_min)],
+            "z_max": bnd_facets[np.isclose(midpoints[2], z_max)],
+        }
+
+        subdomains = {"default": np.arange(mesh.nelements)}
+        return mesh.with_boundaries(boundaries).with_subdomains(subdomains)
+
+    def _make_cylindrical_mesh(self, r_lim, z_lim, h):
+        """
+        Create a cylindrical annulus mesh.
+
+        Parameters
+        ----------
+        r_lim : tuple
+            Radial limits (inner radius, outer radius)
+        z_lim : tuple
+            Vertical limits (z_min, z_max)
+        h : float
+            Target mesh size
+
+        Returns
+        -------
+        skfem.MeshTet
+            Cylindrical mesh with proper boundary tags
+        """
+        skfem = import_optional_dependency("skfem")
+        from scipy.spatial import Delaunay
+
+        r_inner, r_outer = r_lim
+        z_min, z_max = z_lim
+        height = z_max - z_min
+
+        n_radial = max(4, int((r_outer - r_inner) / h) + 1)
+        n_theta_base = max(16, int(2 * np.pi * r_outer / h))
+        r_coords = np.linspace(r_inner, r_outer, n_radial)
+        points_list = []
+        for r in r_coords:
+            if r > 1e-9:
+                n_theta = max(8, int(n_theta_base * (r / r_outer)))
+                thetas = np.linspace(0, 2 * np.pi, n_theta, endpoint=False)
+                for theta in thetas:
+                    points_list.append([r * np.cos(theta), r * np.sin(theta)])
+        points_2d = np.array(points_list)
+        tri_2d = Delaunay(points_2d)
+        triangles_base = tri_2d.simplices
+
+        n_z = max(5, int(height / h))
+        z_coords = np.linspace(z_min, z_max, n_z)
+        n_nodes_per_layer = points_2d.shape[0]
+        nodes_3d = np.zeros((n_nodes_per_layer * n_z, 3))
+        for i, z in enumerate(z_coords):
+            start, end = i * n_nodes_per_layer, (i + 1) * n_nodes_per_layer
+            nodes_3d[start:end, :2] = points_2d
+            nodes_3d[start:end, 2] = z
+        tetrahedra = []
+        for i in range(n_z - 1):
+            for tri in triangles_base:
+                v = [node_idx + (i * n_nodes_per_layer) for node_idx in tri]
+                v_top = [node_idx + ((i + 1) * n_nodes_per_layer) for node_idx in tri]
+                tetrahedra.extend(
+                    [
+                        [v[0], v[1], v[2], v_top[2]],
+                        [v[0], v[1], v_top[1], v_top[2]],
+                        [v[0], v_top[0], v_top[1], v_top[2]],
+                    ]
+                )
+        mesh = skfem.MeshTet(nodes_3d.T, np.array(tetrahedra).T)
+
+        bottom_nodes = set(range(n_nodes_per_layer))
+        top_nodes = set(range(n_nodes_per_layer * (n_z - 1), n_nodes_per_layer * n_z))
+
+        points_2d_r = np.linalg.norm(points_2d, axis=1)
+        tol = h * 0.1
+        inner_nodes_2d = np.where(np.isclose(points_2d_r, r_inner, atol=tol))[0]
+        outer_nodes_2d = np.where(np.isclose(points_2d_r, r_outer, atol=tol))[0]
+
+        inner_wall_nodes = {
+            n + (i * n_nodes_per_layer) for i in range(n_z) for n in inner_nodes_2d
+        }
+        outer_wall_nodes = {
+            n + (i * n_nodes_per_layer) for i in range(n_z) for n in outer_nodes_2d
+        }
+
+        bottom_facets, top_facets, inner_facets, outer_facets = [], [], [], []
+
+        for facet_idx in mesh.boundary_facets():
+            facet_nodes = set(mesh.facets[:, facet_idx])
+            if facet_nodes.issubset(bottom_nodes):
+                bottom_facets.append(facet_idx)
+            elif facet_nodes.issubset(top_nodes):
+                top_facets.append(facet_idx)
+            elif facet_nodes.issubset(inner_wall_nodes):
+                inner_facets.append(facet_idx)
+            elif facet_nodes.issubset(outer_wall_nodes):
+                outer_facets.append(facet_idx)
+
+        boundaries = {}
+        if len(bottom_facets) > 0:
+            boundaries["z_min"] = np.array(bottom_facets)
+        if len(top_facets) > 0:
+            boundaries["z_max"] = np.array(top_facets)
+        if len(inner_facets) > 0:
+            boundaries["r_min"] = np.array(inner_facets)
+        if len(outer_facets) > 0:
+            boundaries["r_max"] = np.array(outer_facets)
+
+        return mesh.with_boundaries(boundaries)
+
+    def __call__(self, lims, npts):
+        """
+        Main entry point called by PyBaMM's Discretisation.
+
+        Parameters
+        ----------
+        lims : dict
+            Dictionary of domain limits
+        npts : dict
+            Dictionary of target number of points
+
+        Returns
+        -------
+        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())
+            coord_sys = "cartesian"
+            mesh = self._make_box_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())
+            coord_sys = "cylindrical polar"
+            mesh = self._make_cylindrical_mesh(r_lim, z_lim, h)
+        else:
+            raise ValueError(f"Unknown geom_type: {self.geom_type}")  # pragma: no cover
+
+        if mesh is None:
+            raise pybamm.DiscretisationError(
+                f"Mesh generation failed for {self.geom_type}"
+            )  # pragma: no cover
+
+        nodes = mesh.p.T
+        elements = mesh.t.T
+        submesh = self.submesh_type(mesh, nodes, elements, coord_sys)
+        return submesh
+
+
+class ScikitFemSubMesh3D(pybamm.SubMesh):
+    """
+    A 3D submesh class for unstructured tetrahedral meshes generated by scikit-fem.
+
+    Parameters
+    ----------
+    skfem_mesh : skfem.MeshTet
+        The scikit-fem mesh object
+    nodes : array_like
+        Array of node coordinates (npts, 3)
+    elements : array_like
+        Array of element connectivity (nelements, 4)
+    coord_sys : str
+        The coordinate system of the submesh
+    tabs : dict, optional
+        Information about tabs (unused in 3D)
+    """
+
+    def __init__(self, skfem_mesh, nodes, elements, coord_sys, tabs=None):
+        super().__init__()
+        skfem = import_optional_dependency("skfem")
+
+        self._skfem_mesh = skfem_mesh
+        self.nodes = nodes
+        self.elements = elements
+        self.coord_sys = coord_sys
+        self.tabs = tabs
+        self.dimension = 3
+        self.npts = nodes.shape[0]
+        self.nelements = elements.shape[0]
+        self.edges = self._skfem_mesh.edges.T
+
+        self.basis = skfem.InteriorBasis(self._skfem_mesh, skfem.ElementTetP1())
+        self.facet_basis = skfem.FacetBasis(self._skfem_mesh, self.basis.elem)
+        self.element = skfem.ElementTetP1()
+
+        if (
+            hasattr(self._skfem_mesh, "boundaries")
+            and self._skfem_mesh.boundaries is not None
+        ):
+            for name in self._skfem_mesh.boundaries:
+                facet_basis = skfem.FacetBasis(
+                    self._skfem_mesh,
+                    self.basis.elem,
+                    facets=self._skfem_mesh.boundaries[name],
+                )
+                setattr(self, f"{name}_basis", facet_basis)
+
+                dofs = self.basis.get_dofs(name).all()
+                setattr(self, f"{name}_dofs", dofs)
+
+                normals = facet_basis.normals
+                setattr(self, f"{name}_normals", normals)
+
+        self._compute_element_volumes()
+        self._compute_element_centroids()
+        self.internal_boundaries = []
+
+    def _compute_element_volumes(self):
+        """Compute volumes of all tetrahedral elements."""
+        vertices = self.nodes[self.elements]
+        v1 = vertices[:, 1] - vertices[:, 0]
+        v2 = vertices[:, 2] - vertices[:, 0]
+        v3 = vertices[:, 3] - vertices[:, 0]
+        self.element_volumes = np.abs(np.einsum("ij,ij->i", np.cross(v1, v2), v3)) / 6
+
+    def _compute_element_centroids(self):
+        """Compute centroids of all tetrahedral elements."""
+        self.element_centroids = np.mean(self.nodes[self.elements], axis=1)
+
+    def to_json(self):
+        """Convert mesh to JSON format."""
+        return {
+            "mesh_type": self.__class__.__name__,
+            "nodes": self.nodes.tolist(),
+            "elements": self.elements.tolist(),
+            "coord_sys": self.coord_sys,
+            "element_volumes": self.element_volumes.tolist(),
+            "element_centroids": self.element_centroids.tolist(),
+        }
+
+    @classmethod
+    def _from_json(cls, json_dict):
+        """Create mesh from JSON format."""
+        skfem = import_optional_dependency("skfem")
+        nodes = np.array(json_dict["nodes"])
+        elements = np.array(json_dict["elements"])
+        mesh = skfem.MeshTet(nodes.T, elements.T)
+        return cls(mesh, nodes, elements, json_dict["coord_sys"])
diff --git a/src/pybamm/parameters/parameter_values.py b/src/pybamm/parameters/parameter_values.py
index 4d6f2970a..19b41b530 100644
--- a/src/pybamm/parameters/parameter_values.py
+++ b/src/pybamm/parameters/parameter_values.py
@@ -584,7 +584,21 @@ class ParameterValues:
         "positive tab": pos. tab bc "no tab": no tab bc}.
         """
         new_boundary_conditions = {}
-        sides = ["left", "right", "negative tab", "positive tab", "no tab"]
+        sides = [
+            "left",
+            "right",
+            "negative tab",
+            "positive tab",
+            "no tab",
+            "x_min",
+            "x_max",
+            "y_min",
+            "y_max",
+            "z_min",
+            "z_max",
+            "r_min",
+            "r_max",
+        ]
         for variable, bcs in model.boundary_conditions.items():
             processed_variable = self.process_symbol(variable)
             new_boundary_conditions[processed_variable] = {}
diff --git a/src/pybamm/spatial_methods/__init__.py b/src/pybamm/spatial_methods/__init__.py
index 56b998853..e141c7eac 100644
--- a/src/pybamm/spatial_methods/__init__.py
+++ b/src/pybamm/spatial_methods/__init__.py
@@ -1,2 +1,2 @@
 __all__ = ['finite_volume', 'scikit_finite_element', 'spatial_method',
-           'spectral_volume', 'zero_dimensional_method']
+           'spectral_volume', 'zero_dimensional_method', 'scikit_finite_element_3d']
diff --git a/src/pybamm/spatial_methods/scikit_finite_element_3d.py b/src/pybamm/spatial_methods/scikit_finite_element_3d.py
new file mode 100644
index 000000000..0e0ad4558
--- /dev/null
+++ b/src/pybamm/spatial_methods/scikit_finite_element_3d.py
@@ -0,0 +1,624 @@
+import numpy as np
+from scipy.sparse import csc_matrix
+from scipy.sparse.linalg import inv
+
+import pybamm
+from pybamm.util import import_optional_dependency
+
+
+class ScikitFiniteElement3D(pybamm.SpatialMethod):
+    """
+    A class which implements the steps specific to the 3D finite element method
+    during discretisation using scikit-fem.
+
+    Parameters
+    ----------
+    options : dict-like, optional
+        A dictionary of options to be passed to the spatial method.
+    """
+
+    def __init__(self, options=None):
+        super().__init__(options)
+        pybamm.citations.register("Gustafsson2020")
+
+    def build(self, mesh):
+        super().build(mesh)
+        for dom in mesh.keys():
+            mesh[dom].npts_for_broadcast_to_nodes = mesh[dom].npts
+
+    def spatial_variable(self, symbol):
+        """
+        Creates a discretised spatial variable for x, y, z, r or theta.
+
+        Parameters
+        ----------
+        symbol : :class:`pybamm.SpatialVariable`
+            The spatial variable to discretise (must be 'x', 'y', 'z', 'r' or 'theta')
+
+        Returns
+        -------
+        :class:`pybamm.Vector`
+            The discretised spatial variable as a vector of nodal values
+        """
+        domain = symbol.domain[0]
+        mesh = self.mesh[domain]
+
+        if symbol.name == "x":
+            entries = mesh.nodes[:, 0][:, np.newaxis]
+        elif symbol.name == "y":
+            entries = mesh.nodes[:, 1][:, np.newaxis]
+        elif symbol.name == "z":
+            entries = mesh.nodes[:, 2][:, np.newaxis]
+        elif symbol.name == "r":
+            entries = np.sqrt(mesh.nodes[:, 0] ** 2 + mesh.nodes[:, 1] ** 2)[
+                :, np.newaxis
+            ]
+        elif symbol.name == "theta":
+            entries = np.arctan2(mesh.nodes[:, 1], mesh.nodes[:, 0])[
+                :, np.newaxis
+            ]  # pragma: no cover
+        else:
+            raise pybamm.GeometryError(
+                f"Spatial variable must be 'x', 'y', 'z', 'r' or 'theta', not {symbol.name}"
+            )  # pragma: no cover
+        return pybamm.Vector(entries, domains=symbol.domains)
+
+    def gradient(self, symbol, discretised_symbol, boundary_conditions):
+        """
+        Matrix-vector multiplication to implement the 3D gradient operator.
+        Returns a Concatenation of [grad_x, grad_y, grad_z] or [grad_r, grad_theta, grad_z].
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Symbol`
+            The symbol that we will take the gradient of.
+        discretised_symbol: :class:`pybamm.Symbol`
+            The discretised symbol of the correct size
+        boundary_conditions : dict
+            The boundary conditions of the model
+
+        Returns
+        -------
+        :class:`pybamm.Concatenation`
+            The 3D gradient as concatenation of x, y, z or z, r, theta components
+        """
+        skfem = import_optional_dependency("skfem")
+        domain = symbol.domain[0]
+        mesh = self.mesh[domain]
+        coord_sys = mesh.coord_sys
+
+        grad_x_matrix, grad_y_matrix, grad_z_matrix = self.gradient_matrix(
+            symbol, boundary_conditions
+        )
+
+        @skfem.BilinearForm
+        def mass_form(u, v, w):
+            return u * v
+
+        mass = skfem.asm(mass_form, mesh.basis)
+        mass_inv = pybamm.Matrix(inv(csc_matrix(mass)))
+
+        grad_x = mass_inv @ (grad_x_matrix @ discretised_symbol)
+        grad_y = mass_inv @ (grad_y_matrix @ discretised_symbol)
+        grad_z = mass_inv @ (grad_z_matrix @ discretised_symbol)
+
+        if coord_sys == "cartesian":
+            grad = pybamm.Concatenation(
+                grad_x, grad_y, grad_z, check_domain=False, concat_fun=np.hstack
+            )
+        elif coord_sys == "cylindrical polar":
+            x_nodes = pybamm.Vector(mesh.nodes[:, 0][:, np.newaxis])
+            y_nodes = pybamm.Vector(mesh.nodes[:, 1][:, np.newaxis])
+            r_nodes = pybamm.sqrt(x_nodes**2 + y_nodes**2)
+            r_nodes.entries[r_nodes.entries < 1e-12] = 1e-12
+
+            cos_theta = x_nodes / r_nodes
+            sin_theta = y_nodes / r_nodes
+
+            grad_r = grad_x * cos_theta + grad_y * sin_theta
+            grad_theta = -grad_x * sin_theta + grad_y * cos_theta
+
+            grad = pybamm.Concatenation(
+                grad_r, grad_theta, grad_z, check_domain=False, concat_fun=np.hstack
+            )
+        else:
+            raise pybamm.DiscretisationError(
+                f"Unknown coordinate system '{coord_sys}'"
+            )  # pragma: no cover
+
+        grad.copy_domains(symbol)
+
+        return grad
+
+    def gradient_squared(self, symbol, discretised_symbol, boundary_conditions):
+        """
+        Multiplication to implement the inner product of the gradient operator
+        with itself in 3D.
+        """
+        grad = self.gradient(symbol, discretised_symbol, boundary_conditions)
+        grad_x, grad_y, grad_z = grad.orphans
+        return grad_x**2 + grad_y**2 + grad_z**2
+
+    def gradient_matrix(self, symbol, boundary_conditions):
+        """
+        Gradient matrices for finite elements in 3D.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Symbol`
+            The symbol for which we want to calculate the gradient matrix
+        boundary_conditions : dict
+            The boundary conditions of the model
+
+        Returns
+        -------
+        tuple of :class:`pybamm.Matrix`
+            The (sparse) finite element gradient matrices.
+        """
+        skfem = import_optional_dependency("skfem")
+        domain = symbol.domain[0]
+        mesh = self.mesh[domain]
+
+        @skfem.BilinearForm
+        def gradient_dx(u, v, w):
+            return u.grad[0] * v
+
+        @skfem.BilinearForm
+        def gradient_dy(u, v, w):
+            return u.grad[1] * v
+
+        @skfem.BilinearForm
+        def gradient_dz(u, v, w):
+            return u.grad[2] * v
+
+        grad_x = skfem.asm(gradient_dx, mesh.basis)
+        grad_y = skfem.asm(gradient_dy, mesh.basis)
+        grad_z = skfem.asm(gradient_dz, mesh.basis)
+
+        return pybamm.Matrix(grad_x), pybamm.Matrix(grad_y), pybamm.Matrix(grad_z)
+
+    def laplacian(self, symbol, discretised_symbol, boundary_conditions):
+        """
+        Matrix-vector multiplication to implement the 3D Laplacian operator.
+        This should be called only after boundary conditions have been properly
+        discretised by the PyBaMM discretisation process.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Symbol`
+            The symbol for which we want to calculate the Laplacian
+        discretised_symbol: :class:`pybamm.Symbol`
+            The discretised symbol of the correct size
+        boundary_conditions : dict
+            The boundary conditions of the model
+
+        Returns
+        -------
+        :class:`pybamm.Symbol`
+            The Laplacian of the symbol
+        """
+        stiffness_matrix = self.stiffness_matrix(symbol, boundary_conditions)
+        boundary_load = self.laplacian_boundary_load(symbol, boundary_conditions)
+
+        return -stiffness_matrix @ discretised_symbol + boundary_load
+
+    def divergence(self, symbol, discretised_symbol, boundary_conditions):
+        """
+        Matrix-vector multiplication to implement the 3D divergence operator.
+        Expects discretised_symbol to be a Concatenation of [Fx, Fy, Fz].
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Symbol`
+            The symbol representing the divergence operation
+        discretised_symbol: :class:`pybamm.Concatenation`
+            The discretised vector field as concatenation [Fx, Fy, Fz]
+        boundary_conditions : dict
+            The boundary conditions of the model
+
+        Returns
+        -------
+        :class:`pybamm.Symbol`
+            The divergence of the vector field
+        """
+        domain_key = symbol.domain[0]
+
+        if not (
+            hasattr(discretised_symbol, "orphans")
+            and len(discretised_symbol.orphans) == 3
+        ):
+            raise ValueError(
+                "divergence expects a concatenation of 3 vector components"
+            )  # pragma: no cover
+
+        Fx_var = pybamm.Variable("Fx", domain=[domain_key])
+        grad_M = self.gradient_matrix(Fx_var, boundary_conditions)
+
+        skfem = import_optional_dependency("skfem")
+        mesh = self.mesh[domain_key]
+        coord_sys = mesh.coord_sys
+
+        if coord_sys == "cartesian":
+            Fx, Fy, Fz = discretised_symbol.orphans
+        elif coord_sys == "cylindrical polar":
+            Fr, Ftheta, Fz_cyl = discretised_symbol.orphans
+            x_nodes = pybamm.Vector(mesh.nodes[:, 0][:, np.newaxis])
+            y_nodes = pybamm.Vector(mesh.nodes[:, 1][:, np.newaxis])
+            r_nodes = pybamm.sqrt(x_nodes**2 + y_nodes**2)
+            r_nodes.entries[r_nodes.entries < 1e-12] = 1e-12
+
+            cos_theta = x_nodes / r_nodes
+            sin_theta = y_nodes / r_nodes
+
+            Fx = Fr * cos_theta - Ftheta * sin_theta
+            Fy = Fr * sin_theta + Ftheta * cos_theta
+            Fz = Fz_cyl
+        else:
+            raise pybamm.DiscretisationError(
+                f"Unknown coordinate system '{coord_sys}'"
+            )  # pragma: no cover
+
+        @skfem.BilinearForm
+        def mass_form(u, v, w):
+            return u * v
+
+        mass = skfem.asm(mass_form, mesh.basis)
+        mass_inv = pybamm.Matrix(inv(csc_matrix(mass)))
+
+        grad_x_T = pybamm.Matrix(grad_M[0].entries)
+        grad_y_T = pybamm.Matrix(grad_M[1].entries)
+        grad_z_T = pybamm.Matrix(grad_M[2].entries)
+
+        div_x = mass_inv @ (grad_x_T @ Fx)
+        div_y = mass_inv @ (grad_y_T @ Fy)
+        div_z = mass_inv @ (grad_z_T @ Fz)
+
+        div_result = div_x + div_y + div_z
+        div_result.clear_domains()
+        return div_result
+
+    def stiffness_matrix(self, symbol, boundary_conditions):
+        """
+        Laplacian (stiffness) matrix for finite elements in 3D.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Symbol`
+            The symbol for which we want to calculate the stiffness matrix
+        boundary_conditions : dict
+            The boundary conditions of the model
+
+        Returns
+        -------
+        :class:`pybamm.Matrix`
+            The (sparse) finite element stiffness matrix
+        """
+        skfem = import_optional_dependency("skfem")
+        domain = symbol.domain[0]
+        mesh = self.mesh[domain]
+
+        @skfem.BilinearForm
+        def stiffness_form(u, v, w):
+            return u.grad[0] * v.grad[0] + u.grad[1] * v.grad[1] + u.grad[2] * v.grad[2]
+
+        stiffness = skfem.asm(stiffness_form, mesh.basis)
+
+        if symbol in boundary_conditions:
+            bcs = boundary_conditions[symbol]
+            for name, (_, bc_type) in bcs.items():
+                if bc_type == "Dirichlet":
+                    boundary_dofs = getattr(mesh, f"{name}_dofs", None)
+                    if boundary_dofs is not None and len(boundary_dofs) > 0:
+                        self.bc_apply(stiffness, boundary_dofs)
+
+        return pybamm.Matrix(stiffness)
+
+    def laplacian_boundary_load(self, symbol_for_laplacian, boundary_conditions_dict):
+        """
+        Assembles the boundary load vector for the 3D Laplacian.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Symbol`
+            The symbol for which we want to calculate the boundary load
+        boundary_conditions : dict
+            The boundary conditions of the model
+
+        Returns
+        -------
+        :class:`pybamm.Vector`
+            The boundary load vector
+        """
+        skfem = import_optional_dependency("skfem")
+        domain = symbol_for_laplacian.domain[0]
+        fem_mesh = self.mesh[domain]
+
+        current_boundary_load_symbol = pybamm.Vector(
+            np.zeros((fem_mesh.npts, 1)), domains=symbol_for_laplacian.domains
+        )
+
+        bcs_for_symbol = boundary_conditions_dict.get(symbol_for_laplacian, {})
+
+        unit_bc_load_form = None
+        if any(bc_type == "Neumann" for _, (_, bc_type) in bcs_for_symbol.items()):
+
+            @skfem.LinearForm
+            def _unit_bc_load_form(v, w):
+                return v
+
+            unit_bc_load_form = _unit_bc_load_form
+        if symbol_for_laplacian in boundary_conditions_dict:
+            for name, (bc_value_symbol, bc_type) in bcs_for_symbol.items():
+                term_contribution = None
+                if bc_type == "Neumann":
+                    numeric_coeffs = skfem.asm(
+                        unit_bc_load_form, getattr(fem_mesh, f"{name}_basis")
+                    )
+                    if numeric_coeffs.ndim == 1:
+                        numeric_coeffs = numeric_coeffs[:, np.newaxis]
+                    term_contribution = bc_value_symbol * pybamm.Vector(numeric_coeffs)
+                    current_boundary_load_symbol += term_contribution
+                elif bc_type == "Dirichlet":
+                    boundary_dofs = getattr(fem_mesh, f"{name}_dofs")
+                    numeric_mask = np.zeros((fem_mesh.npts, 1))
+                    numeric_mask[boundary_dofs, 0] = 1.0
+                    term_contribution = bc_value_symbol * pybamm.Vector(numeric_mask)
+                    current_boundary_load_symbol += term_contribution
+
+        return current_boundary_load_symbol
+
+    def integral(self, child, discretised_child, integration_dimension):
+        """
+        Vector-vector dot product to implement the 3D integral operator.
+
+        Parameters
+        ----------
+        child : :class:`pybamm.Symbol`
+            The symbol being integrated
+        discretised_child : :class:`pybamm.Symbol`
+            The discretised symbol being integrated (vector of nodal values)
+        integration_dimension : str
+            The dimension over which to integrate (e.g. "primary" for volume)
+
+        Returns
+        -------
+        :class:`pybamm.Symbol`
+            The result of the integration (a scalar value).
+        """
+        integration_vector = self.definite_integral_matrix(child)
+        out = integration_vector @ discretised_child
+        return out
+
+    def definite_integral_matrix(self, child, vector_type="row"):
+        """
+        Matrix for definite integral *vector* (vector of ones).
+        The child is used to determine the domain and size.
+
+        Parameters
+        ----------
+        child : :class:`pybamm.Symbol`
+            The symbol whose domain/mesh determines the size.
+        vector_type : str, optional
+            "row" or "column" for the shape of the resulting vector of ones.
+
+        Returns
+        -------
+        :class:`pybamm.Matrix`
+            A matrix representing a vector of ones (either row or column).
+        """
+        skfem = import_optional_dependency("skfem")
+        # get primary domain mesh
+        domain = child.domain[0]
+        mesh = self.mesh[domain]
+
+        # make form for the integral
+        @skfem.LinearForm
+        def integral_form(v, w):
+            return v
+
+        vector = skfem.asm(integral_form, mesh.basis)
+
+        if vector_type == "row":
+            return pybamm.Matrix(vector[np.newaxis, :])
+        elif vector_type == "column":  # pragma: no cover
+            return pybamm.Matrix(vector[:, np.newaxis])
+
+    def indefinite_integral(self, child, discretised_child, direction):
+        """
+        Implementation of the indefinite integral operator in 3D.
+        """
+        raise NotImplementedError(
+            "Indefinite integral not implemented for 3D finite elements"
+        )
+
+    def boundary_integral(self, child, discretised_child, region):
+        """
+        Implementation of the 3D boundary integral operator.
+
+        Parameters
+        ----------
+        child : :class:`pybamm.Symbol`
+            The symbol being integrated
+        discretised_child : :class:`pybamm.Symbol`
+            The discretised symbol being integrated
+        region : str
+            The boundary region over which to integrate
+
+        Returns
+        -------
+        :class:`pybamm.Symbol`
+            The result of the boundary integration
+        """
+        integration_vector = self.boundary_integral_vector(child.domain, region=region)
+        out = integration_vector @ discretised_child
+        out.clear_domains()
+        return out
+
+    def boundary_integral_vector(self, domain, region):
+        """
+        A vector representing an integral operator over the boundary.
+
+        Parameters
+        ----------
+        domain : list
+            The domain(s) of the variable in the integrand
+        region : str
+            The region of the boundary over which to integrate
+
+        Returns
+        -------
+        :class:`pybamm.Matrix`
+            The finite element boundary integral vector
+        """
+        skfem = import_optional_dependency("skfem")
+        mesh = self.mesh[domain[0]]
+
+        @skfem.LinearForm
+        def integral_form(v, w):
+            return v
+
+        if region == "entire":
+            integration_vector = skfem.asm(integral_form, mesh.facet_basis)
+        elif hasattr(mesh, f"{region}_basis"):
+            integration_vector = skfem.asm(
+                integral_form, getattr(mesh, f"{region}_basis")
+            )
+
+        return pybamm.Matrix(integration_vector[np.newaxis, :])
+
+    def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
+        """
+        Returns the boundary value or flux of a variable in 3D.
+
+        Parameters
+        ----------
+        symbol : :class:`pybamm.Symbol`
+            The boundary symbol
+        discretised_child : :class:`pybamm.Symbol`
+            The discretised variable
+        bcs : dict, optional
+            Boundary conditions
+
+        Returns
+        -------
+        :class:`pybamm.Symbol`
+            The boundary value or flux
+        """
+        domain = symbol.children[0].domain
+        region = symbol.side
+
+        if isinstance(symbol, pybamm.BoundaryValue):
+            integration_vector = self.boundary_integral_vector(domain, region=region)
+            boundary_val_vector = integration_vector / (
+                integration_vector @ pybamm.Vector(np.ones(integration_vector.shape[1]))
+            )
+            boundary_value = boundary_val_vector @ discretised_child
+            boundary_value.copy_domains(symbol)
+            return boundary_value
+
+        elif isinstance(symbol, pybamm.BoundaryGradient):
+            raise NotImplementedError(
+                "Boundary gradient not implemented for 3D finite elements"
+            )  # pragma: no cover
+
+    def mass_matrix(self, symbol, boundary_conditions):
+        """
+        Calculates the mass matrix for the 3D finite element method.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Variable`
+            The variable corresponding to the equation
+        boundary_conditions : dict
+            The boundary conditions of the model
+
+        Returns
+        -------
+        :class:`pybamm.Matrix`
+            The (sparse) mass matrix for the spatial method.
+        """
+        return self.assemble_mass_form(symbol, boundary_conditions)
+
+    def boundary_mass_matrix(self, symbol, boundary_conditions):
+        """
+        Calculates the mass matrix assembled over the 3D boundary.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Variable`
+            The variable corresponding to the equation
+        boundary_conditions : dict
+            The boundary conditions of the model
+
+        Returns
+        -------
+        :class:`pybamm.Matrix`
+            The (sparse) boundary mass matrix
+        """
+        return self.assemble_mass_form(symbol, boundary_conditions, region="boundary")
+
+    def assemble_mass_form(self, symbol, boundary_conditions, region="interior"):
+        """
+        Assembles the form of the 3D finite element mass matrix.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Variable`
+            The variable corresponding to the equation
+        boundary_conditions : dict
+            The boundary conditions of the model
+        region: str, optional
+            The domain over which to assemble
+
+        Returns
+        -------
+        :class:`pybamm.Matrix`
+            The (sparse) mass matrix
+        """
+        skfem = import_optional_dependency("skfem")
+        domain = symbol.domain[0]
+        mesh = self.mesh[domain]
+
+        @skfem.BilinearForm
+        def mass_form(u, v, w):
+            return u * v
+
+        if region == "interior":
+            mass = skfem.asm(mass_form, mesh.basis)
+        elif region == "boundary":  # pragma: no cover
+            mass = skfem.asm(mass_form, mesh.facet_basis)
+
+        if symbol in boundary_conditions:
+            bcs = boundary_conditions[symbol]
+            for name, (_, bc_type) in bcs.items():
+                if bc_type == "Dirichlet":
+                    boundary_dofs = getattr(mesh, f"{name}_dofs", None)
+                    if boundary_dofs is not None and len(boundary_dofs) > 0:
+                        self.bc_apply(mass, boundary_dofs, zero=True)
+        return pybamm.Matrix(mass)
+
+    def bc_apply(self, M, boundary_dofs, zero=False):
+        """
+        Adjusts matrices for 3D boundary conditions.
+
+        Parameters
+        ----------
+        M : scipy.sparse matrix
+            Matrix to modify
+        boundary_dofs : array_like
+            Boundary degrees of freedom
+        zero : bool, optional
+            Whether to zero the diagonal entries
+        """
+        M_lil = M.tolil()
+
+        for dof in boundary_dofs:
+            M_lil[dof, :] = 0
+            if not zero:
+                M_lil[dof, dof] = 1.0
+
+        M_csr = M_lil.tocsr()
+        M.data = M_csr.data
+        M.indices = M_csr.indices
+        M.indptr = M_csr.indptr
+        M._shape = M_csr._shape


Test patch:
diff --git a/tests/__init__.py b/tests/__init__.py
index 8227cd065..571db5728 100644
--- a/tests/__init__.py
+++ b/tests/__init__.py
@@ -24,6 +24,8 @@ from .integration.test_models.test_full_battery_models.test_lithium_ion.base_lit
 
 from .shared import (
     get_mesh_for_testing,
+    get_3d_mesh_for_testing,
+    get_unit_3d_mesh_for_testing,
     get_mesh_for_testing_symbolic,
     get_p2d_mesh_for_testing,
     get_size_distribution_mesh_for_testing,
diff --git a/tests/shared.py b/tests/shared.py
index be15d0bf8..317561e22 100644
--- a/tests/shared.py
+++ b/tests/shared.py
@@ -116,6 +116,140 @@ 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":
+        x = pybamm.SpatialVariable("x", ["current collector"])
+        y = pybamm.SpatialVariable("y", ["current collector"])
+        z = pybamm.SpatialVariable("z", ["current collector"])
+        x_max = geom_params.get("x_max", 1.0)
+        y_max = geom_params.get("y_max", 2.0)
+        z_max = geom_params.get("z_max", 3.0)
+        geometry = {
+            "current collector": {
+                x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(x_max)},
+                y: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(y_max)},
+                z: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(z_max)},
+            }
+        }
+        var_pts = {x: 5, y: 5, z: 5}
+    elif geom_type == "cylinder":
+        r = pybamm.SpatialVariable(
+            "r", ["current collector"], coord_sys="cylindrical polar"
+        )
+        z = pybamm.SpatialVariable(
+            "z", ["current collector"], coord_sys="cylindrical polar"
+        )
+        radius = geom_params.get("radius", 1.0)
+        height = geom_params.get("height", 1.0)
+        r_inner = geom_params.get("r_inner", 0.0)
+        geometry = {
+            "current collector": {
+                r: {"min": pybamm.Scalar(r_inner), "max": pybamm.Scalar(radius)},
+                z: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(height)},
+            }
+        }
+        var_pts = {r: 5, z: 5}
+    else:
+        raise ValueError(f"geom_type '{geom_type}' not recognised")
+
+    generator_params = {"h": 0.2}
+    generator_params.update(geom_params)
+    generator = pybamm.ScikitFemGenerator3D(geom_type, **generator_params)
+    submesh_types = {"current collector": generator}
+    return pybamm.Mesh(geometry, submesh_types, var_pts)
+
+
+def get_3d_mesh_for_testing(
+    xpts=5, ypts=5, zpts=5, geom_type="box", include_particles=False, **geom_params
+):
+    param = pybamm.ParameterValues(
+        {
+            "Electrode width [m]": 1.0,
+            "Electrode height [m]": 1.0,
+            "Electrode depth [m]": 1.0,
+            "Negative electrode thickness [m]": 1 / 3,
+            "Separator thickness [m]": 1 / 3,
+            "Positive electrode thickness [m]": 1 / 3,
+        }
+    )
+
+    x = pybamm.SpatialVariable(
+        "x", ["negative electrode", "separator", "positive electrode"]
+    )
+    y = pybamm.SpatialVariable(
+        "y", ["negative electrode", "separator", "positive electrode"]
+    )
+    z = pybamm.SpatialVariable(
+        "z", ["negative electrode", "separator", "positive electrode"]
+    )
+
+    geometry = {
+        "negative electrode": {
+            x: {
+                "min": pybamm.Scalar(0),
+                "max": pybamm.Parameter("Negative electrode thickness [m]"),
+            },
+            y: {
+                "min": pybamm.Scalar(0),
+                "max": pybamm.Parameter("Electrode width [m]"),
+            },
+            z: {
+                "min": pybamm.Scalar(0),
+                "max": pybamm.Parameter("Electrode height [m]"),
+            },
+        },
+        "separator": {
+            x: {
+                "min": pybamm.Parameter("Negative electrode thickness [m]"),
+                "max": pybamm.Parameter("Separator thickness [m]")
+                + pybamm.Parameter("Negative electrode thickness [m]"),
+            },
+            y: {
+                "min": pybamm.Scalar(0),
+                "max": pybamm.Parameter("Electrode width [m]"),
+            },
+            z: {
+                "min": pybamm.Scalar(0),
+                "max": pybamm.Parameter("Electrode height [m]"),
+            },
+        },
+        "positive electrode": {
+            x: {
+                "min": pybamm.Parameter("Separator thickness [m]")
+                + pybamm.Parameter("Negative electrode thickness [m]"),
+                "max": pybamm.Parameter("Positive electrode thickness [m]")
+                + pybamm.Parameter("Separator thickness [m]")
+                + pybamm.Parameter("Negative electrode thickness [m]"),
+            },
+            y: {
+                "min": pybamm.Scalar(0),
+                "max": pybamm.Parameter("Electrode width [m]"),
+            },
+            z: {
+                "min": pybamm.Scalar(0),
+                "max": pybamm.Parameter("Electrode height [m]"),
+            },
+        },
+    }
+    param.process_geometry(geometry)
+
+    # Create generator with parameters
+    generator_params = {"h": 0.2}
+    generator_params.update(geom_params)
+
+    generator = pybamm.ScikitFemGenerator3D(geom_type, **generator_params)
+
+    submesh_types = {
+        "negative electrode": generator,
+        "separator": generator,
+        "positive electrode": generator,
+    }
+
+    var_pts = {x: xpts, y: ypts, z: zpts}
+
+    return pybamm.Mesh(geometry, submesh_types, var_pts)
+
+
 def get_p2d_mesh_for_testing(xpts=None, rpts=10):
     geometry = pybamm.battery_geometry()
     return get_mesh_for_testing(xpts=xpts, rpts=rpts, geometry=geometry)
diff --git a/tests/unit/test_meshes/test_scikit_fem_submesh_3d.py b/tests/unit/test_meshes/test_scikit_fem_submesh_3d.py
new file mode 100644
index 000000000..b731be0be
--- /dev/null
+++ b/tests/unit/test_meshes/test_scikit_fem_submesh_3d.py
@@ -0,0 +1,187 @@
+import numpy as np
+import pytest
+
+import pybamm
+
+
+class TestScikitFemGenerator3D:
+    def test_box_mesh_creation(self):
+        try:
+            from pybamm.meshes.scikit_fem_submeshes_3d import ScikitFemGenerator3D
+        except ImportError:
+            pytest.skip("scikit-fem not available")
+
+        geometry = {
+            "domain": {
+                pybamm.standard_spatial_vars.x: {"min": 0, "max": 1},
+                pybamm.standard_spatial_vars.y: {"min": 0, "max": 2},
+                pybamm.standard_spatial_vars.z: {"min": 0, "max": 3},
+            }
+        }
+
+        mesh_gen = ScikitFemGenerator3D("box", h=0.3)
+
+        mesh = pybamm.Mesh(
+            geometry,
+            {"domain": mesh_gen},
+            {
+                pybamm.standard_spatial_vars.x: None,
+                pybamm.standard_spatial_vars.y: None,
+                pybamm.standard_spatial_vars.z: None,
+            },
+        )
+
+        submesh = mesh["domain"]
+        assert hasattr(submesh, "_skfem_mesh"), "Mesh should have _skfem_mesh attribute"
+        assert hasattr(submesh, "nodes"), "Mesh should have nodes attribute"
+        assert hasattr(submesh, "elements"), "Mesh should have elements attribute"
+        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"
+
+    def test_cylinder_mesh_creation(self):
+        try:
+            from pybamm.meshes.scikit_fem_submeshes_3d import ScikitFemGenerator3D
+        except ImportError:
+            pytest.skip("scikit-fem not available")
+
+        r = pybamm.SpatialVariable("r", ["domain"], coord_sys="cylindrical polar")
+        z = pybamm.SpatialVariable("z", ["domain"], coord_sys="cylindrical polar")
+
+        geometry = {
+            "domain": {
+                r: {"min": 0.1, "max": 1.0},
+                z: {"min": 0, "max": 2.0},
+            }
+        }
+
+        mesh_gen = ScikitFemGenerator3D("cylinder", h=0.5)
+
+        var_pts = {r: None, z: None}
+
+        mesh = pybamm.Mesh(geometry, {"domain": mesh_gen}, var_pts)
+
+        submesh = mesh["domain"]
+        assert hasattr(submesh, "_skfem_mesh")
+        assert submesh.dimension == 3
+        assert submesh.coord_sys == "cylindrical polar"
+
+        nodes = submesh.nodes
+        x_coords, y_coords, z_coords = nodes.T
+        radii = np.sqrt(x_coords**2 + y_coords**2)
+
+        assert np.max(radii) < 1.0 + 1e-7
+        assert np.min(radii) > 0.1 - 1e-7
+        assert np.min(z_coords) >= 0.0 - 1e-7
+        assert np.max(z_coords) <= 2.0 + 1e-7
+
+    def test_invalid_geometry_type(self):
+        try:
+            from pybamm.meshes.scikit_fem_submeshes_3d import ScikitFemGenerator3D
+        except ImportError:
+            pytest.skip("scikit-fem not available")
+
+        with pytest.raises(pybamm.GeometryError, match="geom_type must be one of"):
+            ScikitFemGenerator3D("invalid_type", h=0.3)
+
+
+class TestScikitFemSubMesh3D:
+    def test_submesh_properties(self):
+        try:
+            from pybamm.meshes.scikit_fem_submeshes_3d import (
+                ScikitFemGenerator3D,
+                ScikitFemSubMesh3D,
+            )
+        except ImportError:
+            pytest.skip("scikit-fem not available")
+
+        geometry = {
+            "domain": {
+                pybamm.standard_spatial_vars.x: {"min": 0, "max": 1},
+                pybamm.standard_spatial_vars.y: {"min": 0, "max": 2},
+                pybamm.standard_spatial_vars.z: {"min": 0, "max": 3},
+            }
+        }
+
+        mesh_gen = ScikitFemGenerator3D("box", h=0.4)
+        mesh = pybamm.Mesh(
+            geometry,
+            {"domain": mesh_gen},
+            {
+                pybamm.standard_spatial_vars.x: None,
+                pybamm.standard_spatial_vars.y: None,
+                pybamm.standard_spatial_vars.z: None,
+            },
+        )
+
+        submesh = mesh["domain"]
+
+        assert isinstance(submesh, ScikitFemSubMesh3D), (
+            "Should be ScikitFemSubMesh3D instance"
+        )
+
+        assert hasattr(submesh, "_skfem_mesh"), "Should have _skfem_mesh"
+        assert hasattr(submesh, "nodes"), "Should have nodes"
+        assert hasattr(submesh, "elements"), "Should have elements"
+        assert hasattr(submesh, "npts"), "Should have npts"
+        assert hasattr(submesh, "dimension"), "Should have dimension"
+        assert hasattr(submesh, "coord_sys"), "Should have coord_sys"
+
+        assert submesh.dimension == 3, "Should be 3D"
+        assert submesh.npts == len(submesh.nodes), "npts should equal number of nodes"
+        assert submesh.nodes.shape[1] == 3, "Nodes should be 3D"
+        assert len(submesh.elements) > 0, "Should have elements"
+
+    def test_serialization(self):
+        try:
+            from pybamm.meshes.scikit_fem_submeshes_3d import (
+                ScikitFemGenerator3D,
+                ScikitFemSubMesh3D,
+            )
+        except ImportError:
+            pytest.skip("scikit-fem not available")
+
+        geometry = {
+            "domain": {
+                pybamm.standard_spatial_vars.x: {"min": 0, "max": 1},
+                pybamm.standard_spatial_vars.y: {"min": 0, "max": 2},
+                pybamm.standard_spatial_vars.z: {"min": 0, "max": 3},
+            }
+        }
+
+        mesh_gen = ScikitFemGenerator3D("box", h=0.5)
+        mesh = pybamm.Mesh(
+            geometry,
+            {"domain": mesh_gen},
+            {
+                pybamm.standard_spatial_vars.x: None,
+                pybamm.standard_spatial_vars.y: None,
+                pybamm.standard_spatial_vars.z: None,
+            },
+        )
+
+        submesh = mesh["domain"]
+
+        json_data = submesh.to_json()
+        assert isinstance(json_data, dict), "to_json should return dict"
+        assert "mesh_type" in json_data, "Should have mesh_type"
+        assert json_data["mesh_type"] == "ScikitFemSubMesh3D", (
+            "Should identify as ScikitFemSubMesh3D"
+        )
+        assert "coord_sys" in json_data, "Should have coord_sys"
+        assert "nodes" in json_data, "Should have nodes"
+        assert "elements" in json_data, "Should have elements"
+
+        restored_submesh = ScikitFemSubMesh3D._from_json(json_data)
+        assert isinstance(restored_submesh, ScikitFemSubMesh3D), (
+            "Should restore to ScikitFemSubMesh3D"
+        )
+        assert restored_submesh.coord_sys == submesh.coord_sys, (
+            "Should preserve coord_sys"
+        )
+        assert restored_submesh.npts == submesh.npts, "Should preserve npts"
+
+    def test_invalid_geometry_type(self):
+        with pytest.raises(pybamm.GeometryError, match="geom_type must be one of"):
+            pybamm.ScikitFemGenerator3D("invalid_type", h=0.3)
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
new file mode 100644
index 000000000..60ab6f672
--- /dev/null
+++ b/tests/unit/test_spatial_methods/test_scikit_finite_element_3d.py
@@ -0,0 +1,784 @@
+import numpy as np
+import pytest
+
+import pybamm
+from tests import (
+    get_3d_mesh_for_testing,
+    get_unit_3d_mesh_for_testing,
+)
+
+
+class TestScikitFiniteElement3D:
+    def test_not_implemented(self):
+        mesh = get_3d_mesh_for_testing(include_particles=False)
+        spatial_method = pybamm.ScikitFiniteElement3D()
+        spatial_method.build(mesh)
+        assert spatial_method.mesh == mesh
+        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)
+        spatial_methods = {
+            "negative electrode": pybamm.ScikitFiniteElement3D(),
+            "separator": pybamm.ScikitFiniteElement3D(),
+            "positive electrode": pybamm.ScikitFiniteElement3D(),
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        var = pybamm.Variable("var", domain="negative electrode")
+        x = pybamm.SpatialVariable("x", ["negative electrode"])
+        y = pybamm.SpatialVariable("y", ["negative electrode"])
+        z = pybamm.SpatialVariable("z", ["negative electrode"])
+        disc.set_variable_slices([var])
+        y_test = np.ones(mesh["negative electrode"].npts)
+        boundary_pairs = [("x_min", "x_max"), ("z_min", "z_max"), ("y_min", "y_max")]
+        for side1, side2 in boundary_pairs:
+            for eqn in [
+                pybamm.laplacian(var),
+                pybamm.grad_squared(var),
+                pybamm.laplacian(var) + 2 * var,
+                var * var,
+                pybamm.Integral(var, [x, y, z]) - 1,
+            ]:
+                disc.bcs = {
+                    var: {
+                        side1: (pybamm.Scalar(0), "Dirichlet"),
+                        side2: (pybamm.Scalar(1), "Dirichlet"),
+                    }
+                }
+                eqn_disc = disc.process_symbol(eqn)
+                result = eqn_disc.evaluate(None, y_test)
+                assert result is not None
+                disc.bcs = {
+                    var: {
+                        side1: (pybamm.Scalar(0), "Neumann"),
+                        side2: (pybamm.Scalar(1), "Neumann"),
+                    }
+                }
+                eqn_disc = disc.process_symbol(eqn)
+                result = eqn_disc.evaluate(None, y_test)
+                assert result is not None
+                disc.bcs = {
+                    var: {
+                        side1: (pybamm.Scalar(0), "Neumann"),
+                        side2: (pybamm.Scalar(1), "Dirichlet"),
+                    }
+                }
+                eqn_disc = disc.process_symbol(eqn)
+                result = eqn_disc.evaluate(None, y_test)
+                assert result is not None
+
+    def test_gradient_3d_manufactured(self):
+        mesh = get_unit_3d_mesh_for_testing(h=0.3)
+        disc = pybamm.Discretisation(
+            mesh, {"current collector": pybamm.ScikitFiniteElement3D()}
+        )
+        var = pybamm.Variable("var", domain="current collector")
+        disc.set_variable_slices([var])
+
+        grad_disc = disc.process_symbol(pybamm.grad(var))
+        x, y, z = mesh["current collector"].nodes.T
+        test_func = 2 * x + 3 * y + 4 * z
+
+        grad_eval = grad_disc.evaluate(None, test_func)
+
+        np.testing.assert_allclose(np.mean(grad_eval[:, 0]), 2, rtol=1e-7, atol=1e-6)
+        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):
+        mesh = get_unit_3d_mesh_for_testing(h=0.4)
+        disc = pybamm.Discretisation(
+            mesh, {"current collector": pybamm.ScikitFiniteElement3D()}
+        )
+        var = pybamm.Variable("var", domain="current collector")
+
+        eqn = pybamm.laplacian(var)
+
+        processed_x = disc.process_symbol(
+            pybamm.SpatialVariable("x", ["current collector"])
+        )
+        processed_y = disc.process_symbol(
+            pybamm.SpatialVariable("y", ["current collector"])
+        )
+        processed_z = disc.process_symbol(
+            pybamm.SpatialVariable("z", ["current collector"])
+        )
+
+        u_analytical_sym = processed_x * processed_y * processed_z
+
+        all_boundaries = ["x_min", "x_max", "y_min", "y_max", "z_min", "z_max"]
+        disc.bcs = {
+            var: {name: (u_analytical_sym, "Dirichlet") for name in all_boundaries}
+        }
+        disc.set_variable_slices([var])
+        eqn_disc = disc.process_symbol(eqn)
+
+        x_num, y_num, z_num = mesh["current collector"].nodes.T
+        u_analytical_num = x_num * y_num * z_num
+
+        result = eqn_disc.evaluate(None, u_analytical_num)
+
+        submesh = mesh["current collector"]
+        boundary_dofs = np.unique(
+            np.concatenate(
+                [getattr(submesh, f"{name}_dofs") for name in all_boundaries]
+            )
+        )
+        interior_mask = np.ones(submesh.npts, dtype=bool)
+        interior_mask[boundary_dofs] = False
+        l2_error = np.sqrt(np.mean(result[interior_mask] ** 2))
+        assert l2_error < 1e-10
+
+    def test_divergence_3d_manufactured(self):
+        mesh = get_unit_3d_mesh_for_testing(h=0.3)
+        disc = pybamm.Discretisation(
+            mesh, {"current collector": pybamm.ScikitFiniteElement3D()}
+        )
+
+        u = pybamm.Variable("u", domain="current collector")
+        disc.set_variable_slices([u])
+
+        x, y, z = mesh["current collector"].nodes.T
+        u_exact = 4.5 * x**2
+
+        grad_u = pybamm.grad(u)
+        div_grad_u = pybamm.div(grad_u)
+
+        div_disc = disc.process_symbol(div_grad_u)
+        div_result = div_disc.evaluate(None, u_exact)
+
+        submesh = mesh["current collector"]
+
+        all_boundary_names = ["x_min", "x_max", "y_min", "y_max", "z_min", "z_max"]
+        boundary_dofs = np.unique(
+            np.concatenate(
+                [getattr(submesh, f"{name}_dofs") for name in all_boundary_names]
+            )
+        )
+
+        interior_mask = np.ones(submesh.npts, dtype=bool)
+        interior_mask[boundary_dofs] = False
+
+        mean_on_interior = np.mean(div_result[interior_mask])
+
+        np.testing.assert_allclose(mean_on_interior, 9, rtol=3e-2, atol=1e-2)
+
+    def test_neumann_boundary_conditions(self):
+        mesh = get_unit_3d_mesh_for_testing(h=0.8)
+        spatial_method = pybamm.ScikitFiniteElement3D()
+        spatial_method.build(mesh)
+
+        var = pybamm.Variable("var", domain="current collector")
+        bcs = {var: {"x_max": (pybamm.Scalar(5.0), "Neumann")}}
+
+        boundary_load_disc = spatial_method.laplacian_boundary_load(var, bcs)
+        boundary_load_eval = boundary_load_disc.evaluate()
+
+        x_max_dofs = mesh["current collector"].x_max_dofs
+        assert np.sum(boundary_load_eval) > 0
+        assert np.count_nonzero(boundary_load_eval[x_max_dofs]) > 0
+
+        non_boundary_load = np.delete(boundary_load_eval, x_max_dofs)
+        assert np.count_nonzero(np.abs(non_boundary_load) > 1e-12) == 0
+
+    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
+        )
+        disc = pybamm.Discretisation(
+            mesh, {"current collector": pybamm.ScikitFiniteElement3D()}
+        )
+        T = pybamm.Variable("T", domain="current collector")
+
+        model = pybamm.BaseModel()
+        # Equation: ∇²T + 1 = 0
+        model.algebraic = {T: pybamm.laplacian(T) + pybamm.source(1, T)}
+
+        model.boundary_conditions = {
+            T: {
+                "x_min": (pybamm.Scalar(0), "Dirichlet"),
+                "x_max": (pybamm.Scalar(0), "Dirichlet"),
+                "y_min": (pybamm.Scalar(0), "Neumann"),
+                "y_max": (pybamm.Scalar(0), "Neumann"),
+                "z_min": (pybamm.Scalar(0), "Neumann"),
+                "z_max": (pybamm.Scalar(0), "Neumann"),
+            }
+        }
+        model.initial_conditions = {T: pybamm.Scalar(0)}
+        model.variables = {"T": T}
+
+        disc.process_model(model)
+        solver = pybamm.AlgebraicSolver()
+        solution = solver.solve(model)
+        T_numerical = solution.y.flatten()
+
+        # Define the analytical solution: T(x) = 0.5 * x * (L - x)
+        x_coords = mesh["current collector"].nodes[:, 0]
+        T_analytical = 0.5 * x_coords * (L_x - x_coords)
+
+        np.testing.assert_allclose(T_numerical, T_analytical, atol=1e-2)
+
+    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
+        )
+        spatial_methods = {
+            "current collector": pybamm.ScikitFiniteElement3D(),
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        var = pybamm.Variable("var", domain="current collector")
+        disc.set_variable_slices([var])
+        x = mesh["current collector"].nodes[:, 0]
+        y = mesh["current collector"].nodes[:, 1]
+        z = mesh["current collector"].nodes[:, 2]
+        eqn = pybamm.grad_squared(var)
+        eqn_disc = disc.process_symbol(eqn)
+        test_function = x**2 + y**2 + z**2
+        ans = eqn_disc.evaluate(None, test_function)
+        np.testing.assert_array_less(0, ans)
+
+    def test_manufactured_solution_3d(self):
+        mesh = get_unit_3d_mesh_for_testing(h=0.4)
+        disc = pybamm.Discretisation(
+            mesh, {"current collector": pybamm.ScikitFiniteElement3D()}
+        )
+        submesh = mesh["current collector"]
+        x_vertices, y_vertices, z_vertices = submesh.nodes.T
+
+        Lx = np.max(x_vertices) - np.min(x_vertices)
+        Ly = np.max(y_vertices) - np.min(y_vertices)
+        Lz = np.max(z_vertices) - np.min(z_vertices)
+
+        var = pybamm.Variable("var", domain="current collector")
+        disc.set_variable_slices([var])
+        var_disc = disc.process_symbol(var)
+        np.testing.assert_allclose(
+            var_disc.evaluate(None, z_vertices),
+            z_vertices[:, np.newaxis],
+            rtol=1e-7,
+            atol=1e-6,
+        )
+        np.testing.assert_allclose(
+            var_disc.evaluate(None, 6 * y_vertices),
+            6 * y_vertices[:, np.newaxis],
+            rtol=1e-7,
+            atol=1e-6,
+        )
+        np.testing.assert_allclose(
+            var_disc.evaluate(None, 3 * x_vertices),
+            3 * x_vertices[:, np.newaxis],
+            rtol=1e-7,
+            atol=1e-6,
+        )
+        np.testing.assert_allclose(
+            var_disc.evaluate(None, x_vertices * y_vertices * z_vertices),
+            (x_vertices * y_vertices * z_vertices)[:, np.newaxis],
+            rtol=1e-7,
+            atol=1e-6,
+        )
+
+        var1 = pybamm.Variable("var1", domain="current collector")
+        eqn1 = pybamm.laplacian(var1)
+        disc.bcs = {
+            var1: {
+                "z_min": (pybamm.Scalar(0), "Dirichlet"),
+                "z_max": (pybamm.Scalar(0), "Dirichlet"),
+            }
+        }
+        disc.set_variable_slices([var1])
+        eqn1_disc = disc.process_symbol(eqn1)
+        u1 = np.sin(np.pi * z_vertices / Lz)  # Use Lz
+        mass1 = pybamm.Mass(var1)
+        mass1_disc = disc.process_symbol(mass1)
+        soln1 = -((np.pi / Lz) ** 2) * u1  # Use Lz
+
+        np.testing.assert_allclose(
+            eqn1_disc.evaluate(None, u1).flatten(),
+            mass1_disc.entries @ soln1,
+            rtol=1e-2,  # Looser tolerance for coarse mesh
+            atol=1e-2,
+        )
+
+        var2 = pybamm.Variable("var2", domain="current collector")
+        eqn2 = pybamm.laplacian(var2)
+        disc.bcs = {
+            var2: {
+                "y_min": (pybamm.Scalar(0), "Dirichlet"),
+                "y_max": (pybamm.Scalar(0), "Dirichlet"),
+            }
+        }
+        disc.set_variable_slices([var2])
+        eqn2_disc = disc.process_symbol(eqn2)
+        u2 = (
+            np.cos(np.pi * x_vertices / Lx)
+            * np.sin(np.pi * y_vertices / Ly)
+            * np.cos(np.pi * z_vertices / Lz)
+        )
+        mass2 = pybamm.Mass(var2)
+        mass2_disc = disc.process_symbol(mass2)
+        soln2 = -((np.pi / Lx) ** 2 + (np.pi / Ly) ** 2 + (np.pi / Lz) ** 2) * u2
+
+        np.testing.assert_allclose(
+            eqn2_disc.evaluate(None, u2).flatten(),
+            mass2_disc.entries @ soln2,
+            rtol=1e-2,
+            atol=1e-1,
+        )
+
+        var3 = pybamm.Variable("var3", domain="current collector")
+        eqn3 = pybamm.laplacian(var3)
+        all_boundaries = ["x_min", "x_max", "y_min", "y_max", "z_min", "z_max"]
+        disc.bcs = {
+            var3: {name: (pybamm.Scalar(0), "Dirichlet") for name in all_boundaries}
+        }
+        disc.set_variable_slices([var3])
+        eqn3_disc = disc.process_symbol(eqn3)
+        u3 = (
+            np.sin(np.pi * x_vertices / Lx)
+            * np.sin(np.pi * y_vertices / Ly)
+            * np.sin(np.pi * z_vertices / Lz)
+        )
+        mass3 = pybamm.Mass(var3)
+        mass3_disc = disc.process_symbol(mass3)
+        soln3 = -((np.pi / Lx) ** 2 + (np.pi / Ly) ** 2 + (np.pi / Lz) ** 2) * u3
+        np.testing.assert_allclose(
+            eqn3_disc.evaluate(None, u3).flatten(),
+            mass3_disc.entries @ soln3,
+            rtol=1e-2,
+            atol=1e-1,
+        )
+
+    def test_definite_integral_3d(self):
+        mesh = get_3d_mesh_for_testing(
+            xpts=6, ypts=6, zpts=6, geom_type="box", include_particles=False
+        )
+        spatial_methods = {
+            "negative electrode": pybamm.ScikitFiniteElement3D(),
+            "separator": pybamm.ScikitFiniteElement3D(),
+            "positive electrode": pybamm.ScikitFiniteElement3D(),
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        var = pybamm.Variable("var", domain="negative electrode")
+        x = pybamm.SpatialVariable("x", ["negative electrode"])
+        y = pybamm.SpatialVariable("y", ["negative electrode"])
+        z = pybamm.SpatialVariable("z", ["negative electrode"])
+        integral_eqn = pybamm.Integral(var, [x, y, z])
+        disc.set_variable_slices([var])
+        integral_eqn_disc = disc.process_symbol(integral_eqn)
+        constant_val = 5
+        y_test = constant_val * np.ones(mesh["negative electrode"].npts)
+        fem_mesh = mesh["negative electrode"]
+        lx = fem_mesh.nodes[:, 0].max() - fem_mesh.nodes[:, 0].min()
+        ly = fem_mesh.nodes[:, 1].max() - fem_mesh.nodes[:, 1].min()
+        lz = fem_mesh.nodes[:, 2].max() - fem_mesh.nodes[:, 2].min()
+        expected_volume = lx * ly * lz
+        result = integral_eqn_disc.evaluate(None, y_test)
+        expected = constant_val * expected_volume
+        np.testing.assert_allclose(result, expected, rtol=1e-1, atol=1e-1)
+
+    def test_boundary_value_3d(self):
+        mesh = get_3d_mesh_for_testing(
+            xpts=6, ypts=6, zpts=6, geom_type="box", include_particles=False
+        )
+        spatial_methods = {
+            "negative electrode": pybamm.ScikitFiniteElement3D(),
+            "separator": pybamm.ScikitFiniteElement3D(),
+            "positive electrode": pybamm.ScikitFiniteElement3D(),
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        var = pybamm.Variable("var", domain="negative electrode")
+        disc.set_variable_slices([var])
+        boundary_faces = ["x_min", "x_max", "z_min", "z_max", "y_min", "y_max"]
+        for face in boundary_faces:
+            boundary_val = pybamm.BoundaryValue(var, face)
+            boundary_val_disc = disc.process_symbol(boundary_val)
+            constant_y = np.ones(mesh["negative electrode"].npts)[:, np.newaxis]
+            result = boundary_val_disc.evaluate(None, constant_y)
+            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")
+        disc = pybamm.Discretisation(
+            mesh, {"negative electrode": pybamm.ScikitFiniteElement3D()}
+        )
+
+        var = pybamm.Variable("var", domain="negative electrode")
+        disc.set_variable_slices([var])
+
+        const_value = 5
+        var_with_const_value = pybamm.Vector(
+            np.full(mesh["negative electrode"].npts, const_value)
+        )
+
+        full_boundary_integral = pybamm.BoundaryIntegral(var, region="entire")
+        integral_disc = disc.process_symbol(full_boundary_integral)
+
+        numerical_result = integral_disc.evaluate(y=var_with_const_value.entries)
+
+        submesh = mesh["negative electrode"]
+        lx = submesh.nodes[:, 0].max() - submesh.nodes[:, 0].min()
+        ly = submesh.nodes[:, 1].max() - submesh.nodes[:, 1].min()
+        lz = submesh.nodes[:, 2].max() - submesh.nodes[:, 2].min()
+        analytical_area = 2 * (lx * ly + lx * lz + ly * lz)
+
+        expected_result = const_value * analytical_area
+        np.testing.assert_allclose(numerical_result, expected_result, rtol=1e-2)
+
+    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
+        )
+        spatial_methods = {
+            "current collector": pybamm.ScikitFiniteElement3D(),
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        x = pybamm.SpatialVariable("x", ["current collector"])
+        y = pybamm.SpatialVariable("y", ["current collector"])
+        z = pybamm.SpatialVariable("z", ["current collector"])
+        x_disc = disc.process_symbol(x)
+        y_disc = disc.process_symbol(y)
+        z_disc = disc.process_symbol(z)
+        x_actual = mesh["current collector"].nodes[:, 0][:, np.newaxis]
+        y_actual = mesh["current collector"].nodes[:, 1][:, np.newaxis]
+        z_actual = mesh["current collector"].nodes[:, 2][:, np.newaxis]
+        np.testing.assert_array_equal(x_disc.evaluate(), x_actual)
+        np.testing.assert_array_equal(y_disc.evaluate(), y_actual)
+        np.testing.assert_array_equal(z_disc.evaluate(), z_actual)
+
+    def test_pure_neumann_poisson_3d(self):
+        # Solves grad^2 u = 1, with du/dz = Lz at z=Lz, and du/dn=0 elsewhere.
+        # The solution is constrained to have zero average value.
+        # Exact solution for this case: u(z) = z^2 / 2 - Lz^2 / 6
+
+        u = pybamm.Variable("u", domain="current collector")
+        c = pybamm.Variable("c")
+
+        x = pybamm.SpatialVariable("x", ["current collector"])
+        y = pybamm.SpatialVariable("y", ["current collector"])
+        z = pybamm.SpatialVariable("z", ["current collector"])
+
+        model = pybamm.BaseModel()
+        model.algebraic = {
+            u: pybamm.laplacian(u)
+            - pybamm.source(1, u)
+            + c * pybamm.DefiniteIntegralVector(u, vector_type="column"),
+            c: pybamm.Integral(u, [x, y, z]) + pybamm.Multiplication(0, c),
+        }
+        model.initial_conditions = {u: pybamm.Scalar(0), c: pybamm.Scalar(0)}
+        Lz = 3.0
+        model.boundary_conditions = {
+            u: {
+                "z_min": (pybamm.Scalar(0), "Neumann"),
+                "z_max": (pybamm.Scalar(Lz), "Neumann"),
+                "x_min": (pybamm.Scalar(0), "Neumann"),
+                "x_max": (pybamm.Scalar(0), "Neumann"),
+                "y_min": (pybamm.Scalar(0), "Neumann"),
+                "y_max": (pybamm.Scalar(0), "Neumann"),
+            }
+        }
+        model.variables = {"c": c, "u": u}
+
+        mesh = get_unit_3d_mesh_for_testing(h=0.4, z_max=Lz, x_max=1.0, y_max=2.0)
+
+        spatial_methods = {"current collector": pybamm.ScikitFiniteElement3D()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+
+        solver = pybamm.AlgebraicSolver()
+        solution = solver.solve(model)
+        u_numerical = solution.y[:-1]
+
+        z_coords = mesh["current collector"].nodes[:, 2]
+
+        u_exact_at_nodes = (z_coords**2 / 2) - (Lz**2 / 6)
+
+        np.testing.assert_allclose(
+            u_numerical.flatten(), u_exact_at_nodes.flatten(), rtol=1e-2, atol=5e-2
+        )
+
+    def test_dirichlet_bcs_3d(self):
+        # Manufactured solution u(x) = a*x^2 + b*x + k
+        # Equation: -laplacian(u) + source(2*a, u) = 0  => laplacian(u) = 2*a
+        a_val = 3
+        b_val = 4
+        k_val = 5
+
+        u = pybamm.Variable("u", domain="current collector")
+        model = pybamm.BaseModel()
+
+        model.algebraic = {u: -pybamm.laplacian(u) + pybamm.source(2 * a_val, u)}
+
+        model.boundary_conditions = {
+            u: {
+                "x_min": (pybamm.Scalar(k_val), "Dirichlet"),
+                "x_max": (pybamm.Scalar(a_val + b_val + k_val), "Dirichlet"),
+                "z_min": (pybamm.Scalar(0), "Neumann"),
+                "z_max": (pybamm.Scalar(0), "Neumann"),
+                "y_max": (pybamm.Scalar(0), "Neumann"),
+                "y_min": (pybamm.Scalar(0), "Neumann"),
+            }
+        }
+        model.initial_conditions = {u: pybamm.Scalar(1)}
+        model.variables = {"u": u}
+
+        mesh = get_unit_3d_mesh_for_testing(h=0.4)
+        spatial_methods = {"current collector": pybamm.ScikitFiniteElement3D()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+
+        solver = pybamm.AlgebraicSolver()
+        solution = solver.solve(model)
+
+        x_coords = mesh["current collector"].nodes[:, 0]
+        u_exact_at_nodes = a_val * x_coords**2 + b_val * x_coords + k_val
+
+        np.testing.assert_allclose(
+            solution.y.flatten(), u_exact_at_nodes.flatten(), rtol=1e-2, atol=1e-2
+        )
+
+    def test_scalar_field_discretization(self):
+        mesh = get_unit_3d_mesh_for_testing(h=0.2)
+        disc = pybamm.Discretisation(
+            mesh, {"current collector": pybamm.ScikitFiniteElement3D()}
+        )
+
+        phi = pybamm.Variable("phi", domain="current collector")
+
+        disc.set_variable_slices([phi])
+        phi_disc = disc.process_symbol(phi)
+
+        nodes = mesh["current collector"].nodes
+        x_vals, y_vals, z_vals = nodes.T
+
+        constant_field = np.full_like(x_vals, 5.0)
+        result_1 = phi_disc.evaluate(None, constant_field)
+        np.testing.assert_allclose(
+            result_1.flatten(),
+            constant_field,
+            rtol=1e-12,
+            atol=1e-12,
+            err_msg="Constant scalar field should be preserved exactly",
+        )
+
+        linear_field = 2 * x_vals + 3 * y_vals + 4 * z_vals
+        result_2 = phi_disc.evaluate(None, linear_field)
+        np.testing.assert_allclose(
+            result_2.flatten(),
+            linear_field,
+            rtol=1e-10,
+            atol=1e-10,
+            err_msg="Linear scalar field should be preserved exactly",
+        )
+
+        quadratic_field = x_vals**2 + y_vals**2 + z_vals**2
+        result_3 = phi_disc.evaluate(None, quadratic_field)
+        np.testing.assert_allclose(
+            result_3.flatten(),
+            quadratic_field,
+            rtol=1e-8,
+            atol=1e-8,
+            err_msg="Quadratic scalar field should be well-preserved",
+        )
+
+        mixed_field = x_vals * y_vals + y_vals * z_vals + x_vals * z_vals
+        result_4 = phi_disc.evaluate(None, mixed_field)
+        np.testing.assert_allclose(
+            result_4.flatten(),
+            mixed_field,
+            rtol=1e-8,
+            atol=1e-8,
+            err_msg="Mixed polynomial scalar field should be well-preserved",
+        )
+
+        assert np.abs(np.max(result_3) - np.max(quadratic_field)) < 1e-6
+        assert np.abs(np.min(result_3) - np.min(quadratic_field)) < 1e-6
+
+        mid_idx = len(x_vals) // 2
+        x_mid, y_mid, z_mid = x_vals[mid_idx], y_vals[mid_idx], z_vals[mid_idx]
+        expected_at_mid = x_mid**2 + y_mid**2 + z_mid**2
+        actual_at_mid = result_3[mid_idx, 0]
+
+        np.testing.assert_allclose(
+            actual_at_mid,
+            expected_at_mid,
+            rtol=1e-8,
+            atol=1e-8,
+            err_msg="Scalar field should be accurate at specific points",
+        )
+
+    def test_gradient_cylinder_manufactured(self):
+        mesh = get_unit_3d_mesh_for_testing(
+            geom_type="cylinder", r_inner=0.1, radius=0.5, height=1.0, h=0.15
+        )
+        disc = pybamm.Discretisation(
+            mesh, {"current collector": pybamm.ScikitFiniteElement3D()}
+        )
+        var = pybamm.Variable("var", domain="current collector")
+        disc.set_variable_slices([var])
+        grad_disc = disc.process_symbol(pybamm.grad(var))
+        nodes = mesh["current collector"].nodes
+        x, y, z = nodes.T
+        r = np.sqrt(x**2 + y**2)
+        u_exact_num = r**2 * z
+        grad_u_exact_r = 2 * r * z
+        grad_u_exact_theta = np.zeros_like(r)
+        grad_u_exact_z = r**2
+        grad_eval = grad_disc.evaluate(None, u_exact_num)
+        l2_error_r = np.sqrt(np.mean((grad_eval[:, 0] - grad_u_exact_r) ** 2))
+        l2_error_theta = np.sqrt(np.mean((grad_eval[:, 1] - grad_u_exact_theta) ** 2))
+        l2_error_z = np.sqrt(np.mean((grad_eval[:, 2] - grad_u_exact_z) ** 2))
+        assert l2_error_r < 0.1
+        assert l2_error_theta < 5e-2
+        assert l2_error_z < 0.1
+
+    def test_laplacian_cylinder_manufactured(self):
+        mesh = get_unit_3d_mesh_for_testing(
+            geom_type="cylinder", r_inner=0.2, radius=0.5, height=1.0, h=0.4
+        )
+        model = pybamm.BaseModel()
+        u = pybamm.Variable("u", "current collector")
+        r = pybamm.SpatialVariable("r", "current collector")
+        z = pybamm.SpatialVariable("z", "current collector")
+        u_exact = r**2 * z
+        f = -4 * z
+        model.algebraic = {u: pybamm.laplacian(u) - pybamm.source(f, u)}
+        model.boundary_conditions = {
+            u: {
+                "z_min": (u_exact, "Dirichlet"),
+                "z_max": (u_exact, "Dirichlet"),
+                "r_min": (u_exact, "Dirichlet"),
+                "r_max": (u_exact, "Dirichlet"),
+            }
+        }
+        model.initial_conditions = {u: pybamm.Scalar(0)}
+        model.variables = {"u": u}
+        spatial_methods = {"current collector": pybamm.ScikitFiniteElement3D()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+        solver = pybamm.AlgebraicSolver()
+        solution = solver.solve(model)
+        u_numerical = solution.y.flatten()
+        nodes = mesh["current collector"].nodes
+        x_nodes, y_nodes, z_nodes = nodes.T
+        r_nodes = np.sqrt(x_nodes**2 + y_nodes**2)
+        u_analytical = r_nodes**2 * z_nodes
+        l2_error = np.sqrt(np.mean((u_numerical - u_analytical) ** 2))
+        assert l2_error < 0.08
+
+    def test_divergence_cylinder_manufactured(self):
+        mesh = get_unit_3d_mesh_for_testing(
+            geom_type="cylinder", r_inner=0.1, radius=0.5, height=1.0, h=0.2
+        )
+        disc = pybamm.Discretisation(
+            mesh, {"current collector": pybamm.ScikitFiniteElement3D()}
+        )
+        u = pybamm.Variable("u", domain="current collector")
+        disc.set_variable_slices([u])
+
+        nodes = mesh["current collector"].nodes
+        x, y, z = nodes.T
+        r = np.sqrt(x**2 + y**2)
+
+        # Manufactured solution u = r^3 * z
+        # div(grad(u)) = lap(u) = 9*r*z
+        u_exact_num = r**3 * z
+        div_grad_u_exact_num = 9 * r * z
+
+        div_grad_u_sym = pybamm.div(pybamm.grad(u))
+        div_grad_u_disc = disc.process_symbol(div_grad_u_sym)
+        div_grad_u_eval = div_grad_u_disc.evaluate(None, u_exact_num).flatten()
+
+        l2_error = np.sqrt(np.mean((div_grad_u_eval - div_grad_u_exact_num) ** 2))
+        assert l2_error < 1.5
+
+    def test_cylinder_mesh_properties(self):
+        mesh = get_unit_3d_mesh_for_testing(
+            geom_type="cylinder", radius=0.4, height=0.8
+        )
+        nodes = mesh["current collector"].nodes
+        radii = np.sqrt(nodes[:, 0] ** 2 + nodes[:, 1] ** 2)
+        assert radii.max() < 0.4 + 1e-7
+        assert radii.min() > 0.0 - 1e-7
+        assert nodes[:, 2].min() > 0.0 - 1e-7
+        assert nodes[:, 2].max() < 0.8 + 1e-7
+
+    def test_boundary_integral_3d_box(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
+        )
+        domain = "current collector"
+        disc = pybamm.Discretisation(mesh, {domain: pybamm.ScikitFiniteElement3D()})
+
+        var = pybamm.Variable("var", domain=domain)
+        disc.set_variable_slices([var])
+
+        y_test = np.ones(mesh[domain].npts)
+
+        regions_and_areas = {
+            "x_min": y_max * z_max,
+            "x_max": y_max * z_max,
+            "y_min": x_max * z_max,
+            "y_max": x_max * z_max,
+            "z_min": x_max * y_max,
+            "z_max": x_max * y_max,
+        }
+
+        for region, analytical_area in regions_and_areas.items():
+            integral_op = pybamm.BoundaryIntegral(var, region=region)
+            integral_disc = disc.process_symbol(integral_op)
+            numerical_area = integral_disc.evaluate(y=y_test)
+            np.testing.assert_allclose(
+                numerical_area,
+                analytical_area,
+                rtol=5e-2,
+                err_msg=f"Failed for region: '{region}'",
+            )
+
+    def test_discretise_expression_cylinder(self):
+        mesh = get_unit_3d_mesh_for_testing(
+            geom_type="cylinder", r_inner=0.2, radius=0.5, height=1.0, h=0.1
+        )
+        disc = pybamm.Discretisation(
+            mesh, {"current collector": pybamm.ScikitFiniteElement3D()}
+        )
+
+        r_sym = pybamm.SpatialVariable("r", ["current collector"])
+        z_sym = pybamm.SpatialVariable("z", ["current collector"])
+        expression_sym = r_sym**2 + z_sym
+
+        expression_disc = disc.process_symbol(expression_sym)
+        result_numerical = expression_disc.evaluate().flatten()
+
+        nodes = mesh["current collector"].nodes
+        x, y, z = nodes.T
+        r_analytical = np.sqrt(x**2 + y**2)
+        result_analytical = r_analytical**2 + z
+
+        np.testing.assert_allclose(result_numerical, result_analytical, atol=1e-12)
+
+    def test_scalar_times_gradient_3d(self):
+        mesh = get_unit_3d_mesh_for_testing(h=0.3, x_max=1.0, y_max=1.0, z_max=1.0)
+        disc = pybamm.Discretisation(
+            mesh, {"current collector": pybamm.ScikitFiniteElement3D()}
+        )
+        T = pybamm.Variable("T", domain="current collector")
+        disc.set_variable_slices([T])
+
+        expression_sym = T * pybamm.grad(T)
+        expression_disc = disc.process_symbol(expression_sym)
+
+        nodes = mesh["current collector"].nodes
+        x, y, z = nodes.T
+        T_analytical_func = 2 * x + 3 * y + 4 * z
+
+        numerical_result = expression_disc.evaluate(None, T_analytical_func)
+
+        grad_T_analytical = np.array([2.0, 3.0, 4.0])
+        analytical_result = (
+            T_analytical_func[:, np.newaxis] * grad_T_analytical[np.newaxis, :]
+        )
+
+        np.testing.assert_allclose(numerical_result, analytical_result)