diff --git a/DG_Approximation.py b/DG_Approximation.py
index 63a9cbd6c693c80462584f1bed491ba4ea6cf67b..7aa8b94bf6dd17564113255ee9d9622cf0256cf7 100644
--- a/DG_Approximation.py
+++ b/DG_Approximation.py
@@ -16,7 +16,7 @@ TODO: Put basis initialization for plots in function -> Done
 TODO: Contain cell length in mesh -> Done
 TODO: Contain bounds in mesh -> Done
 TODO: Contain interval length in mesh -> Done
-TODO: Contain number of grid cells in mesh
+TODO: Contain number of grid cells in mesh -> Done
 TODO: Create data dict for mesh separately
 TODO: Create data dict for basis separately
 TODO: Refactor eval_points in do_initial_projection()
@@ -161,7 +161,6 @@ class DGScheme:
         # Unpack keyword arguments
         self._wave_speed = kwargs.pop('wave_speed', 1)
         self._cfl_number = kwargs.pop('cfl_number', 0.2)
-        self._num_grid_cells = kwargs.pop('num_grid_cells', 64)
         self._final_time = kwargs.pop('final_time', 1)
         self._verbose = kwargs.pop('verbose', False)
         self._history_threshold = kwargs.pop('history_threshold',
@@ -175,16 +174,13 @@ class DGScheme:
         self._quadrature = kwargs.pop('quadrature', 'Gauss')
         self._quadrature_config = kwargs.pop('quadrature_config', {})
         self._update_scheme = kwargs.pop('update_scheme', 'SSPRK3')
-
-        polynomial_degree = kwargs.pop('polynomial_degree', 2)
-        self._basis = OrthonormalLegendre(polynomial_degree)
+        self._basis = OrthonormalLegendre(kwargs.pop('polynomial_degree', 2))
 
         # Initialize mesh with two ghost cells on each side
-        left_bound = kwargs.pop('left_bound', -1)
-        right_bound = kwargs.pop('right_bound', 1)
-        self._mesh = Mesh(num_grid_cells=self._num_grid_cells,
-                          num_ghost_cells=2, left_bound=left_bound,
-                          right_bound=right_bound)
+        self._mesh = Mesh(num_grid_cells=kwargs.pop('num_grid_cells', 64),
+                          left_bound=kwargs.pop('left_bound', -1),
+                          right_bound=kwargs.pop('right_bound', 1),
+                          num_ghost_cells=2)
         print(len(self._mesh.cells))
         print(type(self._mesh.cells))
 
@@ -219,12 +215,12 @@ class DGScheme:
             config=self._quadrature_config)
         self._detector = getattr(Troubled_Cell_Detector, self._detector)(
             config=self._detector_config, mesh=self._mesh,
-            wave_speed=self._wave_speed, num_grid_cells=self._num_grid_cells,
-            final_time=self._final_time, basis=self._basis,
-            init_cond=self._init_cond, quadrature=self._quadrature)
+            wave_speed=self._wave_speed, final_time=self._final_time,
+            basis=self._basis, init_cond=self._init_cond,
+            quadrature=self._quadrature)
         self._update_scheme = getattr(Update_Scheme, self._update_scheme)(
-            polynomial_degree=polynomial_degree,
-            num_grid_cells=self._num_grid_cells, detector=self._detector,
+            polynomial_degree=self._basis.polynomial_degree,
+            num_grid_cells=self._mesh.num_grid_cells, detector=self._detector,
             limiter=self._limiter)
 
     def approximate(self, data_file):
@@ -244,8 +240,7 @@ class DGScheme:
         """
         projection = do_initial_projection(
             initial_condition=self._init_cond, mesh=self._mesh,
-            basis=self._basis, quadrature=self._quadrature,
-            num_grid_cells=self._num_grid_cells)
+            basis=self._basis, quadrature=self._quadrature)
 
         time_step = abs(self._cfl_number * self._mesh.cell_len /
                         self._wave_speed)
@@ -324,7 +319,7 @@ class DGScheme:
         projection = do_initial_projection(
             initial_condition=initial_condition, mesh=self._mesh,
             basis=self._basis, quadrature=self._quadrature,
-            num_grid_cells=self._num_grid_cells, adjustment=adjustment)
+            adjustment=adjustment)
 
         return self._basis.calculate_cell_average(
             projection=projection[:, 1:-1], stencil_length=stencil_length,
@@ -332,7 +327,7 @@ class DGScheme:
 
 
 def do_initial_projection(initial_condition, mesh, basis, quadrature,
-                          num_grid_cells, adjustment=0):
+                          adjustment=0):
     """Calculates initial projection.
 
     Calculates a projection at time step 0 and adds ghost cells on both
@@ -348,8 +343,6 @@ def do_initial_projection(initial_condition, mesh, basis, quadrature,
         Basis used for calculation.
     quadrature: Quadrature object
         Quadrature used for evaluation.
-    num_grid_cells : int
-        Number of cells in the mesh. Usually exponential of 2.
     adjustment: float, optional
         Extent of adjustment of each evaluation point in x-direction.
         Default: 0.
@@ -365,7 +358,7 @@ def do_initial_projection(initial_condition, mesh, basis, quadrature,
     output_matrix = [0]
     basis_vector = basis.basis
 
-    for cell in range(num_grid_cells):
+    for cell in range(mesh.num_grid_cells):
         new_row = []
         eval_point = mesh.bounds[0] + (cell+0.5)*mesh.cell_len
 
@@ -382,7 +375,7 @@ def do_initial_projection(initial_condition, mesh, basis, quadrature,
         output_matrix.append(basis.inverse_mass_matrix @ new_row)
 
     # Set ghost cells to respective value
-    output_matrix[0] = output_matrix[num_grid_cells]
+    output_matrix[0] = output_matrix[mesh.num_grid_cells]
     output_matrix.append(output_matrix[1])
 
     # print(np.array(output_matrix).shape)
diff --git a/Plotting.py b/Plotting.py
index d7c5119a1b87e2f9c4bbbfb8392b85d23733b0b8..02d65555e871d29610f92e9f4e0992e1fe40e1a5 100644
--- a/Plotting.py
+++ b/Plotting.py
@@ -125,8 +125,8 @@ def plot_shock_tube(num_grid_cells: int, troubled_cell_history: list,
 
 
 def plot_details(fine_projection: ndarray, fine_mesh: Mesh, basis: Basis,
-                 coarse_projection: ndarray, multiwavelet_coeffs: ndarray,
-                 num_coarse_grid_cells: int) -> None:
+                 coarse_projection: ndarray,
+                 multiwavelet_coeffs: ndarray) -> None:
     """Plots details of projection to coarser mesh.
 
     Parameters
@@ -139,11 +139,9 @@ def plot_details(fine_projection: ndarray, fine_mesh: Mesh, basis: Basis,
         Basis used for calculation.
     multiwavelet_coeffs : ndarray
         Matrix of multiwavelet coefficients.
-    num_coarse_grid_cells : int
-        Number of cells in the coarse mesh (half the cells of the fine mesh).
-        Usually exponential of 2.
 
     """
+    num_coarse_grid_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)
@@ -369,8 +367,8 @@ def plot_approximation_results(data_file: str, directory: str, plot_name: str,
 
 
 def plot_results(projection: ndarray, troubled_cell_history: list,
-                 time_history: list, mesh: Mesh, num_grid_cells: int,
-                 wave_speed: float, final_time: float, basis: Basis,
+                 time_history: list, mesh: Mesh, wave_speed: float,
+                 final_time: float, basis: Basis,
                  quadrature: Quadrature, init_cond: InitialCondition,
                  colors: dict = None, coarse_projection: ndarray = None,
                  multiwavelet_coeffs: ndarray = None) -> None:
@@ -392,8 +390,6 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
         List of value of each time step.
     mesh : Mesh
         Mesh for calculation.
-    num_grid_cells : int
-        Number of cells in the mesh. Usually exponential of 2.
     wave_speed : float
         Speed of wave in rightward direction.
     final_time : float
@@ -419,7 +415,7 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
     colors = _check_colors(colors)
 
     # Plot troubled cells
-    plot_shock_tube(num_grid_cells, troubled_cell_history, time_history)
+    plot_shock_tube(mesh.num_grid_cells, troubled_cell_history, time_history)
 
     # Determine exact and approximate solution
     grid, exact = calculate_exact_solution(
@@ -430,7 +426,7 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
 
     # Plot multiwavelet solution (fine and coarse grid)
     if coarse_projection is not None:
-        coarse_mesh = Mesh(num_grid_cells=num_grid_cells//2,
+        coarse_mesh = Mesh(num_grid_cells=mesh.num_grid_cells//2,
                            num_ghost_cells=1, left_bound=mesh.bounds[0],
                            right_bound=mesh.bounds[1])
 
@@ -445,9 +441,8 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
             colors['coarse_approx'])
 
         # Plot multiwavelet details
-        num_coarse_grid_cells = num_grid_cells//2
         plot_details(projection[:, 1:-1], mesh, basis, coarse_projection,
-                     multiwavelet_coeffs, num_coarse_grid_cells)
+                     multiwavelet_coeffs)
 
         plot_solution_and_approx(grid, exact, approx,
                                  colors['fine_exact'],
@@ -469,7 +464,7 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
     plot_error(grid, exact, approx)
 
     print('p =', basis.polynomial_degree)
-    print('N =', num_grid_cells)
+    print('N =', mesh.num_grid_cells)
     print('maximum error =', max_error)
 
 
diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py
index 24d9a8b4213c9b459ae32d4c9430397d94d0da72..32038c46fc2f7c349f7292d81c850200c8746ec0 100644
--- a/Troubled_Cell_Detector.py
+++ b/Troubled_Cell_Detector.py
@@ -30,7 +30,7 @@ class TroubledCellDetector(ABC):
 
     """
     def __init__(self, config, init_cond, quadrature, basis, mesh,
-                 wave_speed=1, num_grid_cells=64, final_time=1):
+                 wave_speed=1, final_time=1):
         """Initializes TroubledCellDetector.
 
         Parameters
@@ -47,15 +47,12 @@ class TroubledCellDetector(ABC):
             Mesh for calculation.
         wave_speed : float, optional
             Speed of wave in rightward direction. Default: 1.
-        num_grid_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.
 
         """
         self._mesh = mesh
         self._wave_speed = wave_speed
-        self._num_grid_cells = num_grid_cells
         self._final_time = final_time
         self._basis = basis
         self._init_cond = init_cond
@@ -91,14 +88,12 @@ class TroubledCellDetector(ABC):
         pass
 
     def create_data_dict(self, projection):
-        left_bound, right_bound = self._mesh.bounds
         return {'projection': projection, 'wave_speed': self._wave_speed,
-                'num_grid_cells': self._num_grid_cells,
                 'final_time': self._final_time,
                 'basis': {'polynomial_degree': self._basis.polynomial_degree},
-                'mesh': {'num_grid_cells': self._num_grid_cells,
-                         'left_bound': left_bound,
-                         'right_bound': right_bound,
+                'mesh': {'num_grid_cells': self._mesh.num_grid_cells,
+                         'left_bound': self._mesh.bounds[0],
+                         'right_bound': self._mesh.bounds[1],
                          'num_ghost_cells': 2}
                 }
 
@@ -234,7 +229,7 @@ class WaveletDetector(TroubledCellDetector):
         super()._reset(config)
 
         # Set additional necessary parameter
-        self._num_coarse_grid_cells = self._num_grid_cells//2
+        self._num_coarse_grid_cells = self._mesh.num_grid_cells//2
         self._wavelet_projection_left, self._wavelet_projection_right \
             = self._basis.multiwavelet_projection