diff --git a/ANN_Data_Generator.py b/ANN_Data_Generator.py
index ab96ba03515a78b9f4a0e8eb63b35db8298977c5..bcadffa7ed04015742ec98dd8408d4427734af5a 100644
--- a/ANN_Data_Generator.py
+++ b/ANN_Data_Generator.py
@@ -49,7 +49,7 @@ class TrainingDataGenerator:
         self._quadrature_list = [Gauss({'num_nodes': pol_deg+1})
                                  for pol_deg in range(5)]
         self._mesh_list = [Mesh(left_bound=left_bound, right_bound=right_bound,
-                                num_ghost_cells=0, num_grid_cells=2**exp)
+                                num_ghost_cells=0, num_cells=2**exp)
                            for exp in range(3, 12)]
 
     def build_training_data(self, init_cond_list, num_samples, balance=0.5,
diff --git a/DG_Approximation.py b/DG_Approximation.py
index fdd06fa343a089879f28a68f82ea3c976d940099..b6b06bf977b847d7f6f1453e52075df830bfd57a 100644
--- a/DG_Approximation.py
+++ b/DG_Approximation.py
@@ -19,11 +19,12 @@ TODO: Discuss adding kwargs to attributes in documentation
 TODO: Discuss descriptions (matrices, cfl number, right-hand side,
     limiting slope, basis, wavelet, etc.)
 TODO: Discuss referencing info on SSPRK3
+TODO: Discuss name for quadrature mesh (now: grid)
 
 Urgent:
 TODO: Unify use of 'length' and 'len' in naming -> Done
 TODO: Unify use of 'initial_condition' and 'init_cond' in naming -> Done
-TODO: Unify use of 'mesh' and 'grid' in naming
+TODO: Unify use of 'mesh' and 'grid' in naming -> Done
 TODO: Check sign change in stiffness matrix to accommodate negative wave speed
 TODO: Force input_size for each ANN model to be stencil length
 TODO: Induce shift in IC class
@@ -45,6 +46,7 @@ TODO: Add images to report
 TODO: Add verbose output
 
 Critical, but not urgent:
+TODO: Restructure 'calculate_approximate_solution()'
 TODO: Rework Theoretical TC for efficiency
 TODO: Extract object initialization from DGScheme
 TODO: Replace loops with list comprehension if feasible
@@ -141,7 +143,7 @@ class DGScheme:
             Polynomial degree. Default: 2.
         cfl_number : float, optional
             CFL number to ensure stability. Default: 0.2.
-        num_grid_cells : int, optional
+        num_mesh_cells : int, optional
             Number of cells in the mesh. Usually exponential of 2. Default: 64.
         final_time : float, optional
             Final time for which approximation is calculated. Default: 1.
@@ -191,7 +193,7 @@ class DGScheme:
         self._basis = OrthonormalLegendre(kwargs.pop('polynomial_degree', 2))
 
         # Initialize mesh with two ghost cells on each side
-        self._mesh = Mesh(num_grid_cells=kwargs.pop('num_grid_cells', 64),
+        self._mesh = Mesh(num_cells=kwargs.pop('num_mesh_cells', 64),
                           left_bound=kwargs.pop('left_bound', -1),
                           right_bound=kwargs.pop('right_bound', 1),
                           num_ghost_cells=2)
@@ -231,7 +233,7 @@ class DGScheme:
             config=self._detector_config, mesh=self._mesh, basis=self._basis)
         self._update_scheme = getattr(Update_Scheme, self._update_scheme)(
             polynomial_degree=self._basis.polynomial_degree,
-            num_grid_cells=self._mesh.num_grid_cells, detector=self._detector,
+            num_cells=self._mesh.num_cells, detector=self._detector,
             limiter=self._limiter)
 
     def approximate(self, data_file):
@@ -326,7 +328,7 @@ def do_initial_projection(init_cond, mesh, basis, quadrature,
     -------
     ndarray
         Matrix containing projection of size (N+2, p+1) with N being the
-        number of grid cells and p being the polynomial degree of the basis.
+        number of mesh cells and p being the polynomial degree of the basis.
 
     """
     # Initialize matrix and set first entry to accommodate for ghost cell
@@ -347,7 +349,7 @@ def do_initial_projection(init_cond, mesh, basis, quadrature,
         output_matrix.append(basis.inverse_mass_matrix @ new_row)
 
     # Set ghost cells to respective value
-    output_matrix[0] = output_matrix[mesh.num_grid_cells]
+    output_matrix[0] = output_matrix[mesh.num_cells]
     output_matrix.append(output_matrix[1])
 
     # print(np.array(output_matrix).shape)
diff --git a/Plotting.py b/Plotting.py
index 687a37a899251e2b4bd6bf4bed666e1897fef9c5..9ec5354c97fa0f041218bb6d274caf094ff3b62c 100644
--- a/Plotting.py
+++ b/Plotting.py
@@ -94,7 +94,7 @@ def plot_error(grid: ndarray, exact: ndarray, approx: ndarray) -> None:
     plt.title('Errors')
 
 
-def plot_shock_tube(num_grid_cells: int, troubled_cell_history: list,
+def plot_shock_tube(num_cells: int, troubled_cell_history: list,
                     time_history: list) -> None:
     """Plots shock tube.
 
@@ -103,7 +103,7 @@ def plot_shock_tube(num_grid_cells: int, troubled_cell_history: list,
 
     Parameters
     ----------
-    num_grid_cells : int
+    num_cells : int
         Number of cells in the mesh. Usually exponential of 2.
     troubled_cell_history : list
         List of detected troubled cells for each time step.
@@ -116,7 +116,7 @@ def plot_shock_tube(num_grid_cells: int, troubled_cell_history: list,
         current_cells = troubled_cell_history[pos]
         for cell in current_cells:
             plt.plot(cell, time_history[pos], 'k.')
-    plt.xlim((0, num_grid_cells))
+    plt.xlim((0, num_cells))
     plt.xlabel('Cell')
     plt.ylabel('Time')
     plt.title('Shock Tubes')
@@ -139,16 +139,16 @@ def plot_details(fine_projection: ndarray, fine_mesh: Mesh, basis: Basis,
         Matrix of multiwavelet coefficients.
 
     """
-    num_coarse_grid_cells = len(coarse_projection[0])
+    num_coarse_mesh_cells = len(coarse_projection[0])
     averaged_projection = [[coarse_projection[degree][cell]
                             * basis.basis[degree].subs(x, value)
-                            for cell in range(num_coarse_grid_cells)
+                            for cell in range(num_coarse_mesh_cells)
                             for value in [-0.5, 0.5]]
                            for degree in range(basis.polynomial_degree + 1)]
 
     wavelet_projection = [[multiwavelet_coeffs[degree][cell]
                            * basis.wavelet[degree].subs(z, 0.5) * value
-                           for cell in range(num_coarse_grid_cells)
+                           for cell in range(num_coarse_mesh_cells)
                            for value in [(-1) ** (basis.polynomial_degree
                                                   + degree + 1), 1]]
                           for degree in range(basis.polynomial_degree + 1)]
@@ -375,8 +375,8 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
     Plots exact and approximate solution, errors, and troubled cells of a
     projection given its evaluation history.
 
-    If coarse grid and projection are given, solutions are displayed for
-    both coarse and fine grid. Additionally, coefficient details are plotted.
+    If coarse mesh and projection are given, solutions are displayed for
+    both coarse and fine mesh. Additionally, coefficient details are plotted.
 
     Parameters
     ----------
@@ -401,7 +401,7 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
     colors: dict
         Dictionary of colors used for plots.
     coarse_projection: ndarray, optional
-        Matrix of projection on coarse grid for each polynomial degree.
+        Matrix of projection on coarse mesh for each polynomial degree.
         Default: None.
     multiwavelet_coeffs: ndarray, optional
         Matrix of wavelet coefficients. Default: None.
@@ -413,7 +413,7 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
     colors = _check_colors(colors)
 
     # Plot troubled cells
-    plot_shock_tube(mesh.num_grid_cells, troubled_cell_history, time_history)
+    plot_shock_tube(mesh.num_cells, troubled_cell_history, time_history)
 
     # Determine exact and approximate solution
     grid, exact = calculate_exact_solution(
@@ -422,9 +422,9 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
         projection[:, 1:-1], quadrature.nodes,
         basis.polynomial_degree, basis.basis)
 
-    # Plot multiwavelet solution (fine and coarse grid)
+    # Plot multiwavelet solution (fine and coarse mesh)
     if coarse_projection is not None:
-        coarse_mesh = Mesh(num_grid_cells=mesh.num_grid_cells//2,
+        coarse_mesh = Mesh(num_cells=mesh.num_cells//2,
                            num_ghost_cells=1, left_bound=mesh.bounds[0],
                            right_bound=mesh.bounds[1])
 
@@ -462,7 +462,7 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
     plot_error(grid, exact, approx)
 
     print('p =', basis.polynomial_degree)
-    print('N =', mesh.num_grid_cells)
+    print('N =', mesh.num_cells)
     print('maximum error =', max_error)
 
 
diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py
index 347c06c8ea9ccb7487bcb988f8ad393b53d06e6f..cd7f98acdda727b7f956966cf8f7877bf4569b99 100644
--- a/Troubled_Cell_Detector.py
+++ b/Troubled_Cell_Detector.py
@@ -224,7 +224,7 @@ class WaveletDetector(TroubledCellDetector):
             = self._basis.multiwavelet_projection
 
     def _calculate_wavelet_coeffs(self, projection):
-        """Calculates wavelet coefficients used for projection to coarser grid.
+        """Calculates wavelet coefficients used for projection to coarser mesh.
 
         Parameters
         ----------
@@ -238,7 +238,7 @@ class WaveletDetector(TroubledCellDetector):
 
         """
         output_matrix = []
-        for i in range(self._mesh.num_grid_cells):
+        for i in range(self._mesh.num_cells):
             new_entry = 0.5*(
                     projection[:, i+1] @ self._wavelet_projection_left
                     + projection[:, i+2] @ self._wavelet_projection_right)
@@ -282,7 +282,7 @@ class WaveletDetector(TroubledCellDetector):
         Returns
         -------
         ndarray
-            Matrix of projection on coarse grid for each polynomial degree.
+            Matrix of projection on coarse mesh for each polynomial degree.
 
         """
         basis_projection_left, basis_projection_right\
@@ -293,7 +293,7 @@ class WaveletDetector(TroubledCellDetector):
 
         # Calculate projection on coarse mesh
         output_matrix = []
-        for i in range(self._mesh.num_grid_cells//2):
+        for i in range(self._mesh.num_cells//2):
             new_entry = 0.5 * (
                     projection[:, 2 * i] @ basis_projection_left
                     + projection[:, 2 * i + 1] @ basis_projection_right)
@@ -353,18 +353,18 @@ class Boxplot(WaveletDetector):
         self._adjust_outer_fences = config.pop('adjust_outer_fences', True)
         self._extreme_outlier_only = config.pop('extreme_outlier_only', True)
 
-        if self._mesh.num_grid_cells < self._fold_len:
-            self._fold_len = self._mesh.num_grid_cells
+        if self._mesh.num_cells < self._fold_len:
+            self._fold_len = self._mesh.num_cells
 
         self._quantile_method = config.pop('quantile_method', 'weibull')
         num_overlapping_cells = config.pop('num_overlapping_cells', 1)
-        num_folds = self._mesh.num_grid_cells//self._fold_len
+        num_folds = self._mesh.num_cells//self._fold_len
         self._fold_indices = np.zeros([num_folds,
                                        self._fold_len + 2 *
                                        num_overlapping_cells]).astype(np.int32)
         for fold in range(num_folds):
             self._fold_indices[fold] = np.array(
-                [i % self._mesh.num_grid_cells for i in range(
+                [i % self._mesh.num_cells for i in range(
                     fold * self._fold_len - num_overlapping_cells,
                     (fold+1) * self._fold_len + num_overlapping_cells)])
 
@@ -476,9 +476,9 @@ class Theoretical(WaveletDetector):
         troubled_cells = []
         max_avg = np.sqrt(0.5) \
             * max(1, max(abs(projection[0][cell+1])
-                         for cell in range(self._mesh.num_grid_cells)))
+                         for cell in range(self._mesh.num_cells)))
 
-        for cell in range(self._mesh.num_grid_cells):
+        for cell in range(self._mesh.num_cells):
             if self._is_troubled_cell(coeffs, cell, max_avg):
                 troubled_cells.append(cell)
 
@@ -504,6 +504,6 @@ class Theoretical(WaveletDetector):
         """
         max_value = abs(coeffs[cell])/max_avg
         eps = self._cutoff_factor\
-            / (self._mesh.cell_len*self._mesh.num_grid_cells)
+            / (self._mesh.cell_len*self._mesh.num_cells)
 
         return max_value > eps
diff --git a/Update_Scheme.py b/Update_Scheme.py
index 2bc72674c89172271a408600ec66340fe7b37152..d7adb2aa1382537a379f4dc66a9277208171511b 100644
--- a/Update_Scheme.py
+++ b/Update_Scheme.py
@@ -26,14 +26,14 @@ class UpdateScheme(ABC):
         Performs time step.
 
     """
-    def __init__(self, polynomial_degree, num_grid_cells, detector, limiter):
+    def __init__(self, polynomial_degree, num_cells, detector, limiter):
         """Initializes UpdateScheme.
 
         Parameters
         ----------
         polynomial_degree : int
             Polynomial degree.
-        num_grid_cells : int
+        num_cells : int
             Number of cells in the mesh. Usually exponential of 2.
         detector : TroubledCellDetector object
             Troubled cell detector for evaluation.
@@ -43,7 +43,7 @@ class UpdateScheme(ABC):
         """
         # Unpack positional arguments
         self._polynomial_degree = polynomial_degree
-        self._num_grid_cells = num_grid_cells
+        self._num_cells = num_cells
         self._detector = detector
         self._limiter = limiter
 
@@ -169,8 +169,8 @@ class UpdateScheme(ABC):
             degree.
 
         """
-        current_projection[:, 0] = current_projection[:, self._num_grid_cells]
-        current_projection[:, self._num_grid_cells+1] \
+        current_projection[:, 0] = current_projection[:, self._num_cells]
+        current_projection[:, self._num_cells+1] \
             = current_projection[:, 1]
         return current_projection
 
@@ -313,14 +313,14 @@ class SSPRK3(UpdateScheme):
         # Initialize vector and set first entry to accommodate for ghost cell
         right_hand_side = [0]
 
-        for j in range(self._num_grid_cells):
+        for j in range(self._num_cells):
             right_hand_side.append(2*(self._stiffness_matrix
                                       @ current_projection[:, j+1]
                                       + self._boundary_matrix
                                       @ current_projection[:, j]))
 
         # Set ghost cells to respective value
-        right_hand_side[0] = right_hand_side[self._num_grid_cells]
+        right_hand_side[0] = right_hand_side[self._num_cells]
         right_hand_side.append(right_hand_side[1])
 
         return np.transpose(right_hand_side)
diff --git a/config.yaml b/config.yaml
index 2792e37f5c175a77bf2e3b3fde2876c6ead985ae..8cc9d56dad83a9314ccfb0eb9e76270725cbdf88 100644
--- a/config.yaml
+++ b/config.yaml
@@ -1,4 +1,4 @@
-data_dir: 'Sep29'
+data_dir: 'Oct03'
 random_seed: 1234
 
 # Parameter for Approximation with Troubled Cell Detection
@@ -8,7 +8,7 @@ Approximation:
       wave_speed: 1
       polynomial_degree: 2
       cfl_number: 0.2
-      num_grid_cells: 32    # 40 elements work well for Condition 3
+      num_mesh_cells: 32    # 40 elements work well for Condition 3
       final_time: 1
       left_bound: -1
       right_bound: 1
diff --git a/projection_utils.py b/projection_utils.py
index 540a197fe541143e3f8cad9b19821dadc8ce971e..9df58307473561a3bb0e899bc2f877b34409b6d2 100644
--- a/projection_utils.py
+++ b/projection_utils.py
@@ -20,7 +20,7 @@ x = Symbol('x')
 
 
 class Mesh:
-    """Class for mesh/grid.
+    """Class for mesh.
 
     Each cell is characterized by its center.
 
@@ -28,7 +28,7 @@ class Mesh:
     ----------
     mode : str
         Mode for mesh use. Either 'training' or 'evaluation'.
-    num_grid_cells : int
+    num_cells : int
         Number of cells in the mesh (ghost cells notwithstanding). Usually
         exponential of 2.
     bounds : Tuple[float, float]
@@ -42,14 +42,14 @@ class Mesh:
 
     """
 
-    def __init__(self, num_grid_cells: int, num_ghost_cells: int,
+    def __init__(self, num_cells: int, num_ghost_cells: int,
                  left_bound: float, right_bound: float,
                  training_data_mode: bool = False) -> None:
         """Initialize Mesh.
 
         Parameters
         ----------
-        num_grid_cells : int
+        num_cells : int
             Number of cells in the mesh (ghost cells notwithstanding). Has
             to be an exponential of 2.
         num_ghost_cells : int
@@ -63,10 +63,10 @@ class Mesh:
             generation. Default: False.
 
         """
-        self._num_grid_cells = num_grid_cells
+        self._num_cells = num_cells
         self._mode = 'training' if training_data_mode else 'evaluation'
         if not training_data_mode:
-            if not math.log(self._num_grid_cells, 2).is_integer():
+            if not math.log(self._num_cells, 2).is_integer():
                 raise ValueError('The number of cells in the mesh has to be '
                                  'an exponential of 2')
         self._num_ghost_cells = num_ghost_cells
@@ -79,9 +79,9 @@ class Mesh:
         return self._mode
 
     @property
-    def num_grid_cells(self) -> int:
-        """Return number of grid cells."""
-        return self._num_grid_cells
+    def num_cells(self) -> int:
+        """Return number of mesh cells."""
+        return self._num_cells
 
     @property
     def bounds(self) -> Tuple[float, float]:
@@ -96,7 +96,7 @@ class Mesh:
     @property
     def cell_len(self) -> float:
         """Return the length of a mesh cell."""
-        return self.interval_len/self.num_grid_cells
+        return self.interval_len/self.num_cells
 
     @property
     @cache
@@ -114,7 +114,7 @@ class Mesh:
 
     def create_data_dict(self) -> dict:
         """Return dictionary with data necessary to construct mesh."""
-        return {'num_grid_cells': self._num_grid_cells,
+        return {'num_cells': self._num_cells,
                 'left_bound': self._left_bound,
                 'right_bound': self._right_bound,
                 'num_ghost_cells': self._num_ghost_cells}
@@ -134,17 +134,17 @@ class Mesh:
         # Pick random point between left and right bound
         point = np.random.uniform(self._left_bound, self._right_bound)
 
-        # Adjust grid spacing to be within interval if necessary
-        grid_spacing = self.cell_len
+        # Adjust mesh spacing to be within interval if necessary
+        mesh_spacing = self.cell_len
         max_spacing = 2/stencil_len*min(point-self._left_bound,
-                                           self._right_bound-point)
-        while grid_spacing > max_spacing:
-            grid_spacing /= 2
+                                        self._right_bound-point)
+        while mesh_spacing > max_spacing:
+            mesh_spacing /= 2
 
         # Return new mesh instance
-        return Mesh(left_bound=point - stencil_len/2 * grid_spacing,
-                    right_bound=point + stencil_len/2 * grid_spacing,
-                    num_grid_cells=stencil_len, num_ghost_cells=2,
+        return Mesh(left_bound=point - stencil_len/2 * mesh_spacing,
+                    right_bound=point + stencil_len/2 * mesh_spacing,
+                    num_cells=stencil_len, num_ghost_cells=2,
                     training_data_mode=True)