diff --git a/scripts/tcd/projection_utils.py b/scripts/tcd/projection_utils.py
index 6d3c0bf00e021344d36a7096863703e36274f37f..a7c061ea21d4e16efc2874b35d0e42f6aa6e3a08 100644
--- a/scripts/tcd/projection_utils.py
+++ b/scripts/tcd/projection_utils.py
@@ -79,25 +79,15 @@ def calculate_exact_solution(
         Array containing exact evaluation of a function.
 
     """
-    grid = []
-    exact = []
     num_periods = np.floor(wave_speed * final_time / mesh.interval_len)
 
-    for cell_center in mesh.non_ghost_cells:
-        eval_points = cell_center+mesh.cell_len / 2 * \
-                      quadrature.nodes
+    grid = np.repeat(mesh.non_ghost_cells, quadrature.num_nodes) + \
+        mesh.cell_len/2 * np.tile(quadrature.nodes, mesh.num_cells)
+    exact = np.array([init_cond.calculate(
+        mesh=mesh, x=point-wave_speed*final_time+num_periods*mesh.interval_len)
+        for point in grid])
 
-        eval_values = []
-        for eval_point in eval_points:
-            new_entry = init_cond.calculate(mesh=mesh, x=eval_point
-                                            - wave_speed * final_time
-                                            + num_periods * mesh.interval_len)
-            eval_values.append(new_entry)
-
-        grid.append(eval_points)
-        exact.append(eval_values)
-
-    exact = np.reshape(np.array(exact), (1, len(exact) * len(exact[0])))
-    grid = np.reshape(np.array(grid), (1, len(grid) * len(grid[0])))
+    grid = np.reshape(grid, (1, grid.size))
+    exact = np.reshape(exact, (1, exact.size))
 
     return grid, exact