Fix non-deterministic ShapeError in 3D FEM gradient method (#5143)

* adding patch to fix 3D FEM CI errors

* updating changelog

Author: Rishab87

CHANGELOG.md

@@ -17,6 +17,7 @@
 
 ## Bug fixes
 
+- Fix non-deterministic ShapeError in 3D FEM gradient method ([#5143](https://github.com/pybamm-team/PyBaMM/pull/5143))
 - Fixes negative electrode boundary values for half-cell voltage contributions. ([#5139](https://github.com/pybamm-team/PyBaMM/pull/5139))
 - Makes `A_cc` L_z * L_y * number of layers ([#5138](https://github.com/pybamm-team/PyBaMM/pull/5138))
 - Fixes `TimeIntegral` expression node summation when dependent on an input parameter. ([#5119](https://github.com/pybamm-team/PyBaMM/pull/5119))

src/pybamm/spatial_methods/scikit_finite_element_3d.py

@@ -82,6 +82,7 @@ def gradient(self, symbol, discretised_symbol, boundary_conditions):
         :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]
@@ -98,6 +99,12 @@ def mass_form(u, v, w):
         mass = skfem.asm(mass_form, mesh.basis)
         mass_inv = pybamm.Matrix(inv(csc_matrix(mass)))
 
+        if isinstance(discretised_symbol, pybamm.Scalar):  # pragma: no cover
+            zeros = pybamm.Vector(np.zeros((mesh.npts, 1)))
+            grad = pybamm.Concatenation(zeros, zeros, zeros, check_domain=False)
+            grad.copy_domains(symbol)
+            return grad
+
         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)