diff --git a/source/mesh_deformation/diffusion.cc b/source/mesh_deformation/diffusion.cc
index a39cfb075f9..0993292bb56 100644
--- a/source/mesh_deformation/diffusion.cc
+++ b/source/mesh_deformation/diffusion.cc
@@ -281,10 +281,10 @@ namespace aspect
                 cell_matrix = 0;
 
                 // Loop over the quadrature points of the current face
-                for (unsigned int point=0; point<n_fs_face_q_points; ++point)
+                for (unsigned int q=0; q<n_fs_face_q_points; ++q)
                   {
                     // Get the gravity vector to compute the outward direction of displacement
-                    Tensor<1,dim> direction = -(this->get_gravity_model().gravity_vector(fs_fe_face_values.quadrature_point(point)));
+                    Tensor<1,dim> direction = -(this->get_gravity_model().gravity_vector(fs_fe_face_values.quadrature_point(q)));
                     // Normalize direction vector
                     if (direction.norm() > 0.0)
                       direction *= 1./direction.norm();
@@ -299,15 +299,16 @@ namespace aspect
 
                     // Compute the total displacement in the gravity direction,
                     // i.e. the initial topography + any additional mesh displacement.
-                    const double displacement = direction * (displacement_values[point] + initial_topography_values[point]);
+                    const double displacement = direction * (displacement_values[q] + initial_topography_values[q]);
 
                     // To project onto the tangent space of the surface,
                     // we define the projection P:= I- n x n,
                     // with I the unit tensor and n the unit normal to the surface.
                     // The surface gradient then is P times the usual gradient of the shape functions.
                     const Tensor<2, dim, double> projection = unit_symmetric_tensor<dim>() -
-                                                              outer_product(fs_fe_face_values.normal_vector(point), fs_fe_face_values.normal_vector(point));
+                                                              outer_product(fs_fe_face_values.normal_vector(q), fs_fe_face_values.normal_vector(q));
 
+                    const double JxW = fs_fe_face_values.JxW(q);
 
                     // The shape values for the i-loop
                     std::vector<double> phi(dofs_per_cell);
@@ -327,12 +328,12 @@ namespace aspect
                         if (mesh_deformation_dof_handler.get_fe().system_to_component_index(i).first == dim-1)
                           {
                             // Precompute shape values and projected shape value gradients
-                            phi[i] = fs_fe_face_values.shape_value (i, point);
-                            projected_grad_phi[i] = projection * fs_fe_face_values.shape_grad(i, point);
+                            phi[i] = fs_fe_face_values.shape_value (i, q);
+                            projected_grad_phi[i] = projection * fs_fe_face_values.shape_grad(i, q);
 
                             // Assemble the RHS
                             // RHS = M*H_old
-                            cell_vector(i) += phi[i] * displacement * fs_fe_face_values.JxW(point);
+                            cell_vector(i) += phi[i] * displacement * JxW;
 
                             for (unsigned int j=0; j<dofs_per_cell; ++j)
                               {
@@ -344,12 +345,12 @@ namespace aspect
                                     // Matrix := (M+dt*K) = (M+dt*B^T*kappa*B)
                                     cell_matrix(i,j) +=
                                       (
-                                        phi[i] * fs_fe_face_values.shape_value (j, point) +
+                                        phi[i] * fs_fe_face_values.shape_value (j, q) +
                                         this->get_timestep() * diffusivity *
                                         projected_grad_phi[i] *
-                                        (projection * fs_fe_face_values.shape_grad(j, point))
+                                        (projection * fs_fe_face_values.shape_grad(j, q))
                                       )
-                                      * fs_fe_face_values.JxW(point);
+                                      * JxW;
                                   }
                               }
                           }
diff --git a/source/mesh_deformation/free_surface.cc b/source/mesh_deformation/free_surface.cc
index 06e5edfc17c..43b04db1f91 100644
--- a/source/mesh_deformation/free_surface.cc
+++ b/source/mesh_deformation/free_surface.cc
@@ -156,31 +156,37 @@ namespace aspect
 
                 cell_vector = 0;
                 cell_matrix = 0;
-                for (unsigned int point=0; point<n_face_q_points; ++point)
+                for (unsigned int q=0; q<n_face_q_points; ++q)
                   {
                     // Select the direction onto which to project the velocity solution
                     Tensor<1,dim> direction;
                     if ( advection_direction == SurfaceAdvection::normal ) // project onto normal vector
-                      direction = fs_fe_face_values.normal_vector(point);
+                      direction = fs_fe_face_values.normal_vector(q);
                     else if ( advection_direction == SurfaceAdvection::vertical ) // project onto local gravity
-                      direction = this->get_gravity_model().gravity_vector(fs_fe_face_values.quadrature_point(point));
+                      direction = this->get_gravity_model().gravity_vector(fs_fe_face_values.quadrature_point(q));
                     else
                       AssertThrow(false, ExcInternalError());
 
                     direction *= ( direction.norm() > 0.0 ? 1./direction.norm() : 0.0 );
 
+                    const double JxW = fs_fe_face_values.JxW(q);
+
+                    small_vector<Tensor<1,dim>> phi_u(dofs_per_cell);
+                    for (unsigned int i=0; i<dofs_per_cell; ++i)
+                      phi_u[i] = fs_fe_face_values[extract_vel].value(i,q);
+
                     for (unsigned int i=0; i<dofs_per_cell; ++i)
                       {
                         for (unsigned int j=0; j<dofs_per_cell; ++j)
                           {
-                            cell_matrix(i,j) += (fs_fe_face_values[extract_vel].value(j,point) *
-                                                 fs_fe_face_values[extract_vel].value(i,point) ) *
-                                                fs_fe_face_values.JxW(point);
+                            cell_matrix(i,j) += (phi_u[j] *
+                                                 phi_u[i]) *
+                                                JxW;
                           }
 
-                        cell_vector(i) += (fs_fe_face_values[extract_vel].value(i,point) * direction)
-                                          * (velocity_values[point] * direction)
-                                          * fs_fe_face_values.JxW(point);
+                        cell_vector(i) += (phi_u[i] * direction)
+                                          * (velocity_values[q] * direction)
+                                          * JxW;
                       }
                   }
 
diff --git a/source/mesh_deformation/interface.cc b/source/mesh_deformation/interface.cc
index cd15360c4e0..2fb7974f0fc 100644
--- a/source/mesh_deformation/interface.cc
+++ b/source/mesh_deformation/interface.cc
@@ -109,39 +109,49 @@ namespace aspect
 
               this->get_material_model().evaluate(scratch.face_material_model_inputs, scratch.face_material_model_outputs);
 
-              for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point)
+              for (unsigned int q = 0; q < n_face_q_points; ++q)
                 {
                   for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/)
                     {
                       if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
                         {
-                          scratch.phi_u[i_stokes] = scratch.face_finite_element_values[introspection.extractors.velocities].value(i, q_point);
+                          scratch.phi_u[i_stokes] = scratch.face_finite_element_values[introspection.extractors.velocities].value(i, q);
                           ++i_stokes;
                         }
                       ++i;
                     }
 
                   const Tensor<1,dim>
-                  gravity = this->get_gravity_model().gravity_vector(scratch.face_finite_element_values.quadrature_point(q_point));
+                  gravity = this->get_gravity_model().gravity_vector(scratch.face_finite_element_values.quadrature_point(q));
                   const double g_norm = gravity.norm();
 
                   // construct the relevant vectors
-                  const Tensor<1,dim> n_hat = scratch.face_finite_element_values.normal_vector(q_point);
+                  const Tensor<1,dim> n_hat = scratch.face_finite_element_values.normal_vector(q);
                   const Tensor<1,dim> g_hat = (g_norm == 0.0 ? Tensor<1,dim>() : gravity/g_norm);
 
-                  const double pressure_perturbation = scratch.face_material_model_outputs.densities[q_point] *
+                  small_vector<double> phi_u_times_g_hat(stokes_dofs_per_cell);
+                  small_vector<double> phi_u_times_n_hat(stokes_dofs_per_cell);
+                  for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
+                    {
+                      phi_u_times_g_hat[i] = scratch.phi_u[i] * g_hat;
+                      phi_u_times_n_hat[i] = scratch.phi_u[i] * n_hat;
+                    }
+
+                  const double pressure_perturbation = scratch.face_material_model_outputs.densities[q] *
                                                        this->get_timestep() *
                                                        free_surface_theta *
                                                        g_norm;
 
+                  const double JxW = scratch.face_finite_element_values.JxW(q);
+
                   // see Kaus et al 2010 for details of the stabilization term
                   for (unsigned int i=0; i< stokes_dofs_per_cell; ++i)
                     for (unsigned int j=0; j< stokes_dofs_per_cell; ++j)
                       {
                         // The fictive stabilization stress is (phi_u[i].g)*(phi_u[j].n)
                         const double stress_value = -pressure_perturbation*
-                                                    (scratch.phi_u[i]*g_hat) * (scratch.phi_u[j]*n_hat)
-                                                    *scratch.face_finite_element_values.JxW(q_point);
+                                                    phi_u_times_g_hat[i] * phi_u_times_n_hat[j]
+                                                    * JxW;
 
                         data.local_matrix(i,j) += stress_value;
                       }
@@ -867,14 +877,17 @@ namespace aspect
 
             cell_vector = 0;
             cell_matrix = 0;
-            for (unsigned int point=0; point<n_q_points; ++point)
-              for (unsigned int i=0; i<dofs_per_cell; ++i)
-                {
-                  for (unsigned int j=0; j<dofs_per_cell; ++j)
-                    cell_matrix(i,j) += scalar_product( fe_values[extract_vel].gradient(i,point),
-                                                        fe_values[extract_vel].gradient(j,point) ) *
-                                        fe_values.JxW(point);
-                }
+            for (unsigned int q=0; q<n_q_points; ++q)
+              {
+                const double JxW = fe_values.JxW(q);
+                for (unsigned int i=0; i<dofs_per_cell; ++i)
+                  {
+                    for (unsigned int j=0; j<dofs_per_cell; ++j)
+                      cell_matrix(i,j) += scalar_product( fe_values[extract_vel].gradient(i,q),
+                                                          fe_values[extract_vel].gradient(j,q) ) *
+                                          JxW;
+                  }
+              }
 
             mesh_velocity_constraints.distribute_local_to_global (cell_matrix, cell_vector,
                                                                   cell_dof_indices, mesh_matrix, rhs, false);
diff --git a/source/postprocess/dynamic_topography.cc b/source/postprocess/dynamic_topography.cc
index 9e26682f877..aafb2eadb4a 100644
--- a/source/postprocess/dynamic_topography.cc
+++ b/source/postprocess/dynamic_topography.cc
@@ -160,6 +160,7 @@ namespace aspect
                 const double density = out_volume.densities[q];
                 const bool is_compressible = this->get_material_model().is_compressible();
                 const Tensor<1,dim> gravity = this->get_gravity_model().gravity_vector(in_volume.position[q]);
+                const double JxW = fe_volume_values.JxW(q);
 
                 // Set up shape function values
                 for (unsigned int k=0; k<dofs_per_cell; ++k)
@@ -173,20 +174,26 @@ namespace aspect
                   {
                     // Viscous stress part
                     local_vector(i) += 2.0 * eta * ( epsilon_phi_u[i] * in_volume.strain_rate[q]
-                                                     - (is_compressible ? 1./3. * div_phi_u[i] * div_solution[q] : 0.0) ) * fe_volume_values.JxW(q);
+                                                     - (is_compressible ? 1./3. * div_phi_u[i] * div_solution[q] : 0.0) ) * JxW;
                     // Pressure and compressibility parts
-                    local_vector(i) -= div_phi_u[i] * in_volume.pressure[q] * fe_volume_values.JxW(q);
+                    local_vector(i) -= div_phi_u[i] * in_volume.pressure[q] * JxW;
                     // Force part
-                    local_vector(i) -= density * gravity * phi_u[i] * fe_volume_values.JxW(q);
+                    local_vector(i) -= density * gravity * phi_u[i] * JxW;
                   }
               }
             // Assemble the mass matrix for cell face. Since we are using GLL
             // quadrature, the mass matrix will be diagonal, and we can just assemble it into a vector.
             for (unsigned int q=0; q < n_face_q_points; ++q)
-              for (unsigned int i=0; i<dofs_per_cell; ++i)
-                local_mass_matrix(i) += fe_face_values[this->introspection().extractors.velocities].value(i,q) *
-                                        fe_face_values[this->introspection().extractors.velocities].value(i,q) *
-                                        fe_face_values.JxW(q);
+              {
+                const double JxW = fe_face_values.JxW(q);
+                for (unsigned int i=0; i<dofs_per_cell; ++i)
+                  {
+                    const Tensor<1,dim> phi_u = fe_face_values[this->introspection().extractors.velocities].value(i,q);
+                    local_mass_matrix(i) += phi_u *
+                                            phi_u *
+                                            JxW;
+                  }
+              }
 
             cell->distribute_local_to_global(local_vector, rhs_vector);
             cell->distribute_local_to_global(local_mass_matrix, mass_matrix);
@@ -330,21 +337,22 @@ namespace aspect
                 const Point<dim> point = fe_output_values.quadrature_point(q);
                 const Tensor<1,dim> normal = fe_output_values.normal_vector(q);
                 const double gravity_norm = this->get_gravity_model().gravity_vector(point).norm();
+                const double JxW = fe_output_values.JxW(q);
 
                 if (at_upper_surface)
                   {
                     const double delta_rho = out_output.densities[q] - density_above;
                     dynamic_topography += (-stress_output_values[q]*normal - surface_pressure)
-                                          / delta_rho / gravity_norm * fe_output_values.JxW(q);
+                                          / delta_rho / gravity_norm * JxW;
                   }
                 else
                   {
                     const double delta_rho = out_output.densities[q] - density_below;
 
                     dynamic_topography += (-stress_output_values[q]*normal - bottom_pressure)
-                                          / delta_rho / gravity_norm * fe_output_values.JxW(q);
+                                          / delta_rho / gravity_norm * JxW;
                   }
-                face_area += fe_output_values.JxW(q);
+                face_area += JxW;
               }
             // Get the average dynamic topography for the cell
             dynamic_topography = dynamic_topography * (backward_advection ? -1. : 1.) / face_area;
diff --git a/source/simulator/assemblers/stokes.cc b/source/simulator/assemblers/stokes.cc
index 83f5661a68a..de8d3b10ced 100644
--- a/source/simulator/assemblers/stokes.cc
+++ b/source/simulator/assemblers/stokes.cc
@@ -908,16 +908,19 @@ namespace aspect
       const unsigned int n_q_points    = scratch.finite_element_values.n_quadrature_points;
 
       for (unsigned int q=0; q<n_q_points; ++q)
-        for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
-          {
-            if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
-              {
-                scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection.extractors.pressure].value (i, q);
-                data.local_pressure_shape_function_integrals(i_stokes) += scratch.phi_p[i_stokes] * scratch.finite_element_values.JxW(q);
-                ++i_stokes;
-              }
-            ++i;
-          }
+        {
+          const double JxW = scratch.finite_element_values.JxW(q);
+          for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
+            {
+              if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
+                {
+                  scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection.extractors.pressure].value (i, q);
+                  data.local_pressure_shape_function_integrals(i_stokes) += scratch.phi_p[i_stokes] * JxW;
+                  ++i_stokes;
+                }
+              ++i;
+            }
+        }
     }
 
 
@@ -951,13 +954,14 @@ namespace aspect
                                      scratch.face_finite_element_values.quadrature_point(q),
                                      scratch.face_finite_element_values.normal_vector(q));
 
+              const double JxW = scratch.face_finite_element_values.JxW(q);
+
               for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
                 {
                   if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
                     {
                       data.local_rhs(i_stokes) += scratch.face_finite_element_values[introspection.extractors.velocities].value(i,q) *
-                                                  traction *
-                                                  scratch.face_finite_element_values.JxW(q);
+                                                  traction * JxW;
                       ++i_stokes;
                     }
                   ++i;
diff --git a/source/simulator/melt.cc b/source/simulator/melt.cc
index 78e9224a76b..87222eb6b05 100644
--- a/source/simulator/melt.cc
+++ b/source/simulator/melt.cc
@@ -606,6 +606,9 @@ namespace aspect
           const double K_D = this->get_melt_handler().limited_darcy_coefficient(melt_outputs->permeabilities[q] / melt_outputs->fluid_viscosities[q],
                                                                                 is_melt_cell);
 
+          const double JxW = scratch.face_finite_element_values.JxW(q);
+          const Tensor<1,dim> normal_vector = scratch.face_finite_element_values.normal_vector(q);
+
           for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
             {
               const unsigned int component_index_i = fe.system_to_component_index(i).first;
@@ -616,9 +619,9 @@ namespace aspect
                   data.local_rhs(i_stokes) += (scratch.face_finite_element_values[ex_p_f].value(i, q)
                                                * this->get_pressure_scaling() * K_D *
                                                (density_f
-                                                * (scratch.face_finite_element_values.normal_vector(q) * gravity)
+                                                * (normal_vector * gravity)
                                                 - grad_p_f[q])
-                                               * scratch.face_finite_element_values.JxW(q));
+                                               * JxW);
                   ++i_stokes;
                 }
               ++i;
@@ -842,6 +845,8 @@ namespace aspect
               density_c_P_melt = 1.0;
             }
 
+          const double JxW = scratch.finite_element_values.JxW(q);
+
           // do the actual assembly. note that we only need to loop over the advection
           // shape functions because these are the only contributions we compute here
           for (unsigned int i=0; i<advection_dofs_per_cell; ++i)
@@ -853,8 +858,7 @@ namespace aspect
                   * (gamma + melt_transport_RHS)
                   + scratch.phi_field[i]
                   * reaction_term)
-                 *
-                 scratch.finite_element_values.JxW(q);
+                 * JxW;
 
               for (unsigned int j=0; j<advection_dofs_per_cell; ++j)
                 {
@@ -869,8 +873,7 @@ namespace aspect
                           + (factor * scratch.phi_field[i] * scratch.phi_field[j])) *
                        (density_c_P_melt)
                        + time_step * scratch.phi_field[i] * scratch.phi_field[j] * melt_transport_LHS
-                     )
-                     * scratch.finite_element_values.JxW(q);
+                     ) * JxW;
                 }
             }
         }
@@ -990,18 +993,21 @@ namespace aspect
       const FEValuesExtractors::Scalar ex_p_f = introspection.variable("fluid pressure").extractor_scalar();
 
       for (unsigned int q=0; q<n_q_points; ++q)
-        for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
-          {
-            const unsigned int component_index_i = fe.system_to_component_index(i).first;
+        {
+          const double JxW = scratch.finite_element_values.JxW(q);
+          for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
+            {
+              const unsigned int component_index_i = fe.system_to_component_index(i).first;
 
-            if (is_velocity_or_pressures(introspection,p_c_component_index,p_f_component_index,component_index_i))
-              {
-                scratch.phi_p[i_stokes] = scratch.finite_element_values[ex_p_f].value (i, q);
-                data.local_pressure_shape_function_integrals(i_stokes) += scratch.phi_p[i_stokes] * scratch.finite_element_values.JxW(q);
-                ++i_stokes;
-              }
-            ++i;
-          }
+              if (is_velocity_or_pressures(introspection,p_c_component_index,p_f_component_index,component_index_i))
+                {
+                  scratch.phi_p[i_stokes] = scratch.finite_element_values[ex_p_f].value (i, q);
+                  data.local_pressure_shape_function_integrals(i_stokes) += scratch.phi_p[i_stokes] * JxW;
+                  ++i_stokes;
+                }
+              ++i;
+            }
+        }
     }
 
     template <int dim>
@@ -1050,6 +1056,8 @@ namespace aspect
                                      scratch.face_finite_element_values.quadrature_point(q),
                                      scratch.face_finite_element_values.normal_vector(q));
 
+              const double JxW = scratch.face_finite_element_values.JxW(q);
+
               for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
                 {
                   const unsigned int component_index_i = fe.system_to_component_index(i).first;
@@ -1058,7 +1066,7 @@ namespace aspect
                     {
                       data.local_rhs(i_stokes) += scratch.face_finite_element_values[introspection.extractors.velocities].value(i,q) *
                                                   traction *
-                                                  scratch.face_finite_element_values.JxW(q);
+                                                  JxW;
                       ++i_stokes;
                     }
                   ++i;
@@ -1910,8 +1918,17 @@ namespace aspect
                 const Tensor<1,dim> n_hat = scratch.face_finite_element_values.normal_vector(q_point);
                 const Tensor<1,dim> g_hat = (g_norm == 0.0 ? Tensor<1,dim>() : gravity/g_norm);
 
+                small_vector<double> phi_u_times_g_hat(stokes_dofs_per_cell);
+                small_vector<double> phi_u_times_n_hat(stokes_dofs_per_cell);
+                for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
+                  {
+                    phi_u_times_g_hat[i] = scratch.phi_u[i] * g_hat;
+                    phi_u_times_n_hat[i] = scratch.phi_u[i] * n_hat;
+                  }
+
                 const double pressure_perturbation = scratch.face_material_model_outputs.densities[q_point] *
                                                      this->get_timestep() * free_surface_theta * g_norm;
+                const double JxW = scratch.face_finite_element_values.JxW(q_point);
 
                 // see Kaus et al 2010 for details of the stabilization term
                 for (unsigned int i=0; i< stokes_dofs_per_cell; ++i)
@@ -1919,8 +1936,8 @@ namespace aspect
                     {
                       // The fictive stabilization stress is (phi_u[i].g)*(phi_u[j].n)
                       const double stress_value = - pressure_perturbation
-                                                  * (scratch.phi_u[i] * g_hat) * (scratch.phi_u[j] * n_hat)
-                                                  * scratch.face_finite_element_values.JxW(q_point);
+                                                  * phi_u_times_g_hat[i] * phi_u_times_n_hat[j]
+                                                  * JxW;
 
                       data.local_matrix(i,j) += stress_value;
                     }
diff --git a/source/utilities.cc b/source/utilities.cc
index 059705ee6fc..3a586b63b07 100644
--- a/source/utilities.cc
+++ b/source/utilities.cc
@@ -2891,26 +2891,29 @@ namespace aspect
 
             cell_vector = 0;
             local_mass_matrix = 0;
-            for (unsigned int point=0; point<n_q_points; ++point)
-              for (unsigned int i=0; i<dofs_per_cell; ++i)
-                {
-                  if (dof_handler.get_fe().system_to_component_index(i).first == component_index)
-                    cell_vector(i) +=
-                      rhs_values[point] *
-                      fe_values[extractor].value(i,point) *
-                      fe_values.JxW(point);
-
-                  for (unsigned int j=0; j<dofs_per_cell; ++j)
-                    if ((dof_handler.get_fe().system_to_component_index(i).first ==
-                         component_index)
-                        &&
-                        (dof_handler.get_fe().system_to_component_index(j).first ==
-                         component_index))
-                      local_mass_matrix(j,i) += (fe_values[extractor].value(i,point) * fe_values[extractor].value(j,point) *
-                                                 fe_values.JxW(point));
-                    else if (i == j)
-                      local_mass_matrix(i,j) = 1.;
-                }
+            for (unsigned int q=0; q<n_q_points; ++q)
+              {
+                const double JxW = fe_values.JxW(q);
+                for (unsigned int i=0; i<dofs_per_cell; ++i)
+                  {
+                    if (dof_handler.get_fe().system_to_component_index(i).first == component_index)
+                      cell_vector(i) +=
+                        rhs_values[q] *
+                        fe_values[extractor].value(i,q) *
+                        JxW;
+
+                    for (unsigned int j=0; j<dofs_per_cell; ++j)
+                      if ((dof_handler.get_fe().system_to_component_index(i).first ==
+                           component_index)
+                          &&
+                          (dof_handler.get_fe().system_to_component_index(j).first ==
+                           component_index))
+                        local_mass_matrix(j,i) += (fe_values[extractor].value(i,q) * fe_values[extractor].value(j,q) *
+                                                   JxW);
+                      else if (i == j)
+                        local_mass_matrix(i,j) = 1.;
+                  }
+              }
 
             // now invert the local mass matrix and multiply it with the rhs
             local_mass_matrix.gauss_jordan();
diff --git a/source/volume_of_fluid/assembler.cc b/source/volume_of_fluid/assembler.cc
index 1d92fb49bfc..4e9e8afd043 100644
--- a/source/volume_of_fluid/assembler.cc
+++ b/source/volume_of_fluid/assembler.cc
@@ -97,6 +97,8 @@ namespace aspect
 
       for (unsigned int q = 0; q< n_q_points; ++q)
         {
+          const double JxW = scratch.finite_element_values.JxW(q);
+
           // Init FE field vals
           for (unsigned int k=0; k<volume_of_fluid_dofs_per_cell; ++k)
             scratch.phi_field[k] = scratch.finite_element_values[solution_field].value(main_fe.component_to_system_index(solution_component, k), q);
@@ -104,13 +106,13 @@ namespace aspect
           for (unsigned int i = 0; i<volume_of_fluid_dofs_per_cell; ++i)
             {
               data.local_rhs[i] += scratch.old_field_values[q] *
-                                   scratch.finite_element_values.JxW(q);
+                                   JxW;
               for (unsigned int j=0; j<volume_of_fluid_dofs_per_cell; ++j)
                 data.local_matrix (i, j) += scratch.phi_field[i] *
                                             scratch.phi_field[j] *
-                                            scratch.finite_element_values.JxW(q);
+                                            JxW;
             }
-          scratch.volume += scratch.finite_element_values.JxW(q);
+          scratch.volume += JxW;
         }
 
       for (const unsigned int face_no : cell->face_indices())
diff --git a/source/volume_of_fluid/setup_initial_conditions.cc b/source/volume_of_fluid/setup_initial_conditions.cc
index 41209956832..96a2371a119 100644
--- a/source/volume_of_fluid/setup_initial_conditions.cc
+++ b/source/volume_of_fluid/setup_initial_conditions.cc
@@ -100,12 +100,13 @@ namespace aspect
         double volume_of_fluid_val = 0.0;
         double cell_vol = 0.0;
 
-        for (unsigned int i = 0; i < fe_init.n_quadrature_points; ++i)
+        for (unsigned int q = 0; q < fe_init.n_quadrature_points; ++q)
           {
-            const double fraction_at_point = this->get_initial_composition_manager().initial_composition(fe_init.quadrature_point(i),
+            const double fraction_at_point = this->get_initial_composition_manager().initial_composition(fe_init.quadrature_point(q),
                                              field.composition_index);
-            volume_of_fluid_val += fraction_at_point * fe_init.JxW (i);
-            cell_vol += fe_init.JxW(i);
+            const double JxW = fe_init.JxW(q);
+            volume_of_fluid_val += fraction_at_point * JxW;
+            cell_vol += JxW;
           }
 
         volume_of_fluid_val /= cell_vol;
@@ -184,11 +185,11 @@ namespace aspect
           {
 
             // For each quadrature point compute an approximation to the fluid fraction in the surrounding region
-            for (unsigned int i = 0; i < fe_init.n_quadrature_points; ++i)
+            for (unsigned int q = 0; q < fe_init.n_quadrature_points; ++q)
               {
                 double d = 0.0;
                 Tensor<1, dim, double> grad;
-                Point<dim> xU = quadrature.point (i);
+                Point<dim> xU = quadrature.point (q);
 
                 // Get an approximation to local normal at the closest interface (level set gradient)
                 // and the distance to the closest interface (value of level set function)
@@ -208,8 +209,9 @@ namespace aspect
                   }
                 // Use the basic fluid fraction formula to compute an approximation to the fluid fraction
                 const double fraction_at_point = VolumeOfFluid::Utilities::compute_fluid_fraction (grad, d);
-                volume_of_fluid_val += fraction_at_point * fe_init.JxW (i);
-                cell_vol += fe_init.JxW (i);
+                const double JxW = fe_init.JxW(q);
+                volume_of_fluid_val += fraction_at_point * JxW;
+                cell_vol += JxW;
               }
             volume_of_fluid_val /= cell_vol;
           }