diff --git a/source/mesh_refinement/composition_gradient.cc b/source/mesh_refinement/composition_gradient.cc
index af6105b1dc1..36cc2a0717d 100644
--- a/source/mesh_refinement/composition_gradient.cc
+++ b/source/mesh_refinement/composition_gradient.cc
@@ -43,7 +43,7 @@ namespace aspect
 
       for (const unsigned int k : this->introspection().get_compositional_field_base_element_indices())
         {
-          const Quadrature<dim> quadrature (composition.base_element(k).get_unit_support_points());
+          const Quadrature<dim> quadrature (this->get_fe().base_element(k).get_unit_support_points());
           const unsigned int dofs_per_cell = quadrature.size();
           FEValues<dim> fe_values (this->get_mapping(),
                                    this->get_fe(),
@@ -54,16 +54,18 @@ namespace aspect
           // we have to extract them in this structure
           std::vector<Tensor<1,dim>> composition_gradients (quadrature.size());
 
-          for (const auto &cell : this->get_dof_handler().active_cell_iterators())
-            if (cell->is_locally_owned())
-              {
-                const unsigned int idx = cell->active_cell_index();
-                fe_values.reinit(cell);
-
-                for (const unsigned int c : this->introspection().get_compositional_field_indices_with_base_element(k))
+          for (const unsigned int c : this->introspection().get_compositional_field_indices_with_base_element(k))
+            {
+              const typename Simulator<dim>::AdvectionField composition = Simulator<dim>::AdvectionField::composition(c);
+              for (const auto &cell : this->get_dof_handler().active_cell_iterators())
+                if (cell->is_locally_owned())
                   {
-                    fe_values[composition.scalar_extractor(c)].get_function_gradients (this->get_solution(),
-                                                                                       composition_gradients);
+                    const unsigned int idx = cell->active_cell_index();
+                    fe_values.reinit(cell);
+
+
+                    fe_values[composition.scalar_extractor(this->introspection())].get_function_gradients (this->get_solution(),
+                        composition_gradients);
 
                     // For each composition dof, write into the output vector the
                     // composition gradient. Note that quadrature points and dofs
@@ -79,8 +81,9 @@ namespace aspect
                     // will yield convergence of the error indicators to zero as h->0)
                     this_indicator[idx] *= std::pow(cell->diameter(), power);
                   }
-                indicators.add(composition_scaling_factors[c], this_indicator);
-              }
+
+              indicators.add(composition_scaling_factors[c], this_indicator);
+            }
         }
     }