diff --git a/source/particle/generator/quadrature_points.cc b/source/particle/generator/quadrature_points.cc
index 73fadc9be07..cebf81aa9df 100644
--- a/source/particle/generator/quadrature_points.cc
+++ b/source/particle/generator/quadrature_points.cc
@@ -32,7 +32,8 @@ namespace aspect
       void
       QuadraturePoints<dim>::generate_particles(std::multimap<Particles::internal::LevelInd, Particle<dim>> &particles)
       {
-        const QGauss<dim> quadrature_formula(this->get_parameters().stokes_velocity_degree + 1);
+        const Quadrature<dim> quadrature_formula
+          = this->introspection().quadratures.velocities;
 
         types::particle_index n_particles_to_generate = quadrature_formula.size() * this->get_triangulation().n_locally_owned_active_cells();
         types::particle_index prefix_sum = 0;
diff --git a/source/postprocess/melt_statistics.cc b/source/postprocess/melt_statistics.cc
index 05c4785c562..15b4e2ba107 100644
--- a/source/postprocess/melt_statistics.cc
+++ b/source/postprocess/melt_statistics.cc
@@ -36,7 +36,8 @@ namespace aspect
     MeltStatistics<dim>::execute (TableHandler &statistics)
     {
       // create a quadrature formula based on the temperature element alone.
-      const QGauss<dim> quadrature_formula (this->get_fe().base_element(this->introspection().base_elements.temperature).degree+1);
+      const Quadrature<dim> quadrature_formula
+        = this->introspection().quadratures.temperature;
       const unsigned int n_q_points = quadrature_formula.size();
       std::vector<std::vector<double>> composition_values (this->n_compositional_fields(),
                                                            std::vector<double> (n_q_points));
diff --git a/source/simulator/entropy_viscosity.cc b/source/simulator/entropy_viscosity.cc
index f2ea4875dbd..b8c69f505a8 100644
--- a/source/simulator/entropy_viscosity.cc
+++ b/source/simulator/entropy_viscosity.cc
@@ -44,7 +44,8 @@ namespace aspect
       return numbers::signaling_nan<double>();
 
     // record maximal entropy on Gauss quadrature points
-    const QGauss<dim> quadrature_formula (advection_field.polynomial_degree(introspection)+1);
+    const Quadrature<dim> quadrature_formula
+      = introspection.quadratures.compositional_fields;
     const unsigned int n_q_points = quadrature_formula.size();
 
     const FEValuesExtractors::Scalar field = advection_field.scalar_extractor(introspection);
diff --git a/source/simulator/nullspace.cc b/source/simulator/nullspace.cc
index a97a65cc557..9afe1619a42 100644
--- a/source/simulator/nullspace.cc
+++ b/source/simulator/nullspace.cc
@@ -254,10 +254,10 @@ namespace aspect
 
     // compute and remove net linear momentum from velocity field, by computing
     // \int \rho (v + v_const) = 0
-
-    QGauss<dim> quadrature(parameters.stokes_velocity_degree+1);
+    const Quadrature<dim> quadrature = introspection.quadratures.velocities;
     const unsigned int n_q_points = quadrature.size();
-    FEValues<dim> fe(*mapping, finite_element, quadrature,
+    FEValues<dim> fe(*mapping, finite_element,
+                     quadrature,
                      UpdateFlags(update_quadrature_points | update_JxW_values | update_values | update_gradients));
 
     Tensor<1,dim> local_momentum;
@@ -364,8 +364,9 @@ namespace aspect
     // \int \rho u \cdot r_orth = \omega  * \int \rho x^2    ( 2 dimensions)
     // \int \rho r \times u =  I \cdot \omega  (3 dimensions)
 
-    QGauss<dim> quadrature(parameters.stokes_velocity_degree+1);
-    QGauss<dim-1> surface_quadrature(parameters.stokes_velocity_degree+1);
+    const Quadrature<dim> quadrature
+      = introspection.quadratures.velocities;
+    const QGauss<dim-1> surface_quadrature(parameters.stokes_velocity_degree+1);
 
     const unsigned int n_q_points = (limit_to_top_faces == false) ? quadrature.size() : surface_quadrature.size();
     UpdateFlags flags = update_quadrature_points | update_JxW_values | update_values | update_gradients;