diff --git a/Basis_Function.py b/Basis_Function.py
index af43bdfbbc32f2b01fccbb35067bc20ae51e9c7d..832355fbba93d6ccdd9d5e57886ecf6af4af310a 100644
--- a/Basis_Function.py
+++ b/Basis_Function.py
@@ -22,6 +22,8 @@ class Vector:
         Array of basis.
     wavelet : ndarray
         Array of wavelet.
+    inv_mass: ndarray
+        Inverse mass matrix.
 
     Methods
     -------
@@ -29,6 +31,8 @@ class Vector:
         Returns basis vector.
     get_wavelet_vector
         Returns wavelet vector.
+    get_inverse_mass_matrix()
+        Returns inverse mass matrix.
     get_basis_projections()
         Returns basis projections.
     get_wavelet_projections()
@@ -47,6 +51,7 @@ class Vector:
         self._polynomial_degree = polynomial_degree
         self._basis = self._build_basis_vector(x)
         self._wavelet = self._build_wavelet_vector(z)
+        self._inv_mass = self._build_inverse_mass_matrix()
 
     def get_basis_vector(self):
         """Returns basis vector."""
@@ -88,6 +93,21 @@ class Vector:
         """
         return []
 
+    def get_inverse_mass_matrix(self):
+        """Returns inverse mass matrix."""
+        return self._inv_mass
+
+    def _build_inverse_mass_matrix(self):
+        """Constructs inverse mass matrix.
+
+        Returns
+        -------
+        ndarray
+            Inverse mass matrix.
+
+        """
+        pass
+
     def get_basis_projections(self):
         """Returns basis projection."""
         pass
@@ -238,6 +258,18 @@ class OrthonormalLegendre(Legendre):
         raise ValueError('Invalid value: Alpert\'s wavelet is only available \
                          up to degree 4 for this application')
 
+    def _build_inverse_mass_matrix(self):
+        mass_matrix = []
+        for i in range(self._polynomial_degree+1):
+            new_row = []
+            for j in range(self._polynomial_degree+1):
+                new_entry = 0.0
+                if i == j:
+                    new_entry = 1.0
+                new_row.append(new_entry)
+            mass_matrix.append(new_row)
+        return np.array(mass_matrix)
+
     def get_basis_projections(self):
         """Returns basis projection.
 
diff --git a/DG_Approximation.py b/DG_Approximation.py
index e8487ee58702b264c596d1156dfc848cd7878f9f..661a6776cd76d0c725a403a1e93af247c5cbcaab 100644
--- a/DG_Approximation.py
+++ b/DG_Approximation.py
@@ -4,7 +4,7 @@
 
 Urgent:
 TODO: Extract do_initial_projection() from DGScheme -> Done
-TODO: Move inverse mass matrix to basis
+TODO: Move inverse mass matrix to basis -> Done
 TODO: Extract calculate_cell_average() from TCD
 TODO: Extract calculate_[...]_solution() from Plotting
 TODO: Extract plotting from TCD completely
@@ -364,18 +364,6 @@ def do_initial_projection(initial_condition, basis, quadrature,
         number of grid cells and p being the polynomial degree.
 
     """
-    # Set inverse mass matrix
-    mass_matrix = []
-    for i in range(polynomial_degree+1):
-        new_row = []
-        for j in range(polynomial_degree+1):
-            new_entry = 0.0
-            if i == j:
-                new_entry = 1.0
-            new_row.append(new_entry)
-        mass_matrix.append(new_row)
-    inv_mass = np.array(mass_matrix)
-
     # Initialize matrix and set first entry to accommodate for ghost cell
     output_matrix = [0]
     basis_vector = basis.get_basis_vector()
@@ -386,19 +374,16 @@ def do_initial_projection(initial_condition, basis, quadrature,
         eval_point = left_bound + (cell+0.5)*cell_len
 
         for degree in range(polynomial_degree + 1):
-            new_entry = sum(
-                initial_condition.calculate(
-                    eval_point + cell_len/2
-                    * quadrature.get_eval_points()[point]
-                    - adjustment)
+            new_row.append(np.float64(sum(initial_condition.calculate(
+                eval_point + cell_len/2
+                * quadrature.get_eval_points()[point] - adjustment)
                 * basis_vector[degree].subs(
                     x, quadrature.get_eval_points()[point])
                 * quadrature.get_weights()[point]
-                for point in range(quadrature.get_num_points()))
-            new_row.append(np.float64(new_entry))
+                for point in range(quadrature.get_num_points()))))
 
         new_row = np.array(new_row)
-        output_matrix.append(inv_mass @ new_row)
+        output_matrix.append(basis.get_inverse_mass_matrix() @ new_row)
 
     # Set ghost cells to respective value
     output_matrix[0] = output_matrix[num_grid_cells]