diff --git a/DG_Approximation.py b/DG_Approximation.py
index 1727f6e6d3a6fad4bdffa31682e55e188483cfdf..61d6741d8b1a757c0361ced33fe26f47ab3a35a2 100644
--- a/DG_Approximation.py
+++ b/DG_Approximation.py
@@ -21,10 +21,8 @@ TODO: Contain number of grid cells in mesh -> Done
 TODO: Create data dict for mesh separately -> Done
 TODO: Create data dict for basis separately -> Done
 TODO: Remove unnecessary instance variables from TCD -> Done
-TODO: Refactor eval_points in do_initial_projection()
+TODO: Replace eval_points with mesh in do_initial_projection() -> Done
 TODO: Replace getter with property attributes for quadrature
-TODO: Check whether ghost cells are handled/set correctly
-TODO: Ensure uniform use of mesh and grid
 TODO: Remove use of DGScheme from ANN_Data_Generator
 TODO: Find error in centering for ANN training
 TODO: Adapt TCD from Soraya
@@ -34,6 +32,8 @@ TODO: Add TC condition to only flag cell if left-adjacent one is flagged as
 TODO: Move plot_approximation_results() into plotting script
 TODO: Add verbose output
 TODO: Improve file naming (e.g. use '.' instead of '__')
+TODO: Check whether ghost cells are handled/set correctly
+TODO: Ensure uniform use of mesh and grid
 
 Critical, but not urgent:
 TODO: Combine ANN workflows
@@ -51,6 +51,7 @@ TODO: Use cfl_number for updating, not just time
 Currently not critical:
 TODO: Unify use of 'length' and 'len' in naming
 TODO: Replace loops with list comprehension if feasible
+TODO: Replace loops/list comprehension with vectorization if feasible
 TODO: Check whether 'projection' is always a ndarray
 TODO: Check whether all instance variables are sensible
 TODO: Rename files according to standard
@@ -357,17 +358,15 @@ def do_initial_projection(initial_condition, mesh, basis, quadrature,
     """
     # Initialize matrix and set first entry to accommodate for ghost cell
     output_matrix = [0]
-    basis_vector = basis.basis
 
-    for cell in range(mesh.num_grid_cells):
+    for eval_point in mesh.non_ghost_cells:
         new_row = []
-        eval_point = mesh.bounds[0] + (cell+0.5)*mesh.cell_len
 
         for degree in range(basis.polynomial_degree + 1):
             new_row.append(np.float64(sum(initial_condition.calculate(
                 mesh, eval_point + mesh.cell_len/2
                 * quadrature.get_eval_points()[point] - adjustment)
-                * basis_vector[degree].subs(
+                * basis.basis[degree].subs(
                     x, quadrature.get_eval_points()[point])
                 * quadrature.get_weights()[point]
                 for point in range(quadrature.get_num_points()))))