From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Mon, 24 Aug 2015 00:30:23 +0000 (-0500)
Subject: Make the Mapping and FE output objects members of FEValues.
X-Git-Tag: v8.4.0-rc2~555^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=639fcbb1b167799c9b0595a6f7a3f742e0d54634;p=dealii.git

Make the Mapping and FE output objects members of FEValues.

Historically, the fields in the internal::FEValues::MappingRelatedData and
internal::FEValues::FiniteElementRelatedData classes were part of the (now
removed) base class FEValuesData. These two classes neatly split the fields of
this base class into separate categories, but they continued to enter the
FEValuesBase class via protected inheritance. This is not only slightly awkward,
but also not the easiest approach to understand if you start looking at stuff.

This patch is in essence incredibly boring: instead of having two protected
base classes, it introduces two protected member variables. The remainder
of the patch is then simply an exercise in making sure every use of the
henceforth member variables now accessese the members of the two new
class-type member variables. The only interesting aspect (and that's where
everything becomes much clearer with this design) is that when we call
Mapping::fill_fe_values() or the same function in FiniteElement, we no
longer need to implicitly cast down '*this' to the base class reference,
but can instead just use the new member variable.

Other than that the only noteworthy part of the patch is the introduction
of memory_consumption() functions for the two classes previously split out
from FEValuesData.
---

diff --git a/include/deal.II/fe/fe_update_flags.h b/include/deal.II/fe/fe_update_flags.h
index 9cf01649f3..5a93350dec 100644
--- a/include/deal.II/fe/fe_update_flags.h
+++ b/include/deal.II/fe/fe_update_flags.h
@@ -372,6 +372,12 @@ namespace internal
       void initialize (const unsigned int n_quadrature_points,
                        const UpdateFlags  flags);
 
+      /**
+       * Compute and return an estimate for the memory consumption (in
+       * bytes) of this object.
+       */
+      std::size_t memory_consumption () const;
+
       /**
        * Store an array of weights times the Jacobi determinant at the quadrature
        * points. This function is reset each time reinit() is called. The Jacobi
@@ -442,6 +448,12 @@ namespace internal
                        const FiniteElement<dim,spacedim> &fe,
                        const UpdateFlags         flags);
 
+      /**
+       * Compute and return an estimate for the memory consumption (in
+       * bytes) of this object.
+       */
+      std::size_t memory_consumption () const;
+
       /**
        * Storage type for shape values. Each row in the matrix denotes the values
        * of a single shape function at the different points, columns are for a
diff --git a/include/deal.II/fe/fe_values.h b/include/deal.II/fe/fe_values.h
index c73190e914..16aff9940b 100644
--- a/include/deal.II/fe/fe_values.h
+++ b/include/deal.II/fe/fe_values.h
@@ -1345,15 +1345,19 @@ namespace internal
  * FiniteElement to schedule auxiliary data fields for updating. Still, it is
  * recommended to give <b>all</b> needed update flags to FEValues.
  *
- * The mechanisms by which this class works is also discussed on the page on
- * @ref UpdateFlags.
+ *
+ * <h3>Internals about the implementation</h3>
+ *
+ * The mechanisms by which this class work are discussed on the page on
+ * @ref UpdateFlags "Update flags" and about the
+ * @ref FE_vs_Mapping_vs_FEValues "How Mapping, FiniteElement, and FEValues work together".
+ *
  *
  * @ingroup feaccess
  * @author Wolfgang Bangerth, 1998, 2003, Guido Kanschat, 2001
  */
 template <int dim, int spacedim>
-class FEValuesBase : protected dealii::internal::FEValues::MappingRelatedData<dim, spacedim>,
-  protected dealii::internal::FEValues::FiniteElementRelatedData<dim, spacedim>,
+class FEValuesBase :
   public Subscriptor
 {
 public:
@@ -1398,7 +1402,9 @@ public:
    * Destructor.
    */
   ~FEValuesBase ();
-  /// @name ShapeAccess Access to shape function values. These fields are filled by the finite element
+
+
+  /// @name ShapeAccess Access to shape function values. These fields are filled by the finite element.
   //@{
 
   /**
@@ -2344,25 +2350,43 @@ protected:
   maybe_invalidate_previous_present_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell);
 
   /**
-   * Storage for the mapping object.
+   * A pointer to the mapping object associated with this FEValues object.
    */
   const SmartPointer<const Mapping<dim,spacedim>,FEValuesBase<dim,spacedim> > mapping;
 
   /**
-   * Store the finite element for later use.
+   * A pointer to the internal data object of mapping, obtained from
+   * Mapping::get_data(), Mapping::get_face_data(), or
+   * Mapping::get_subface_data().
    */
-  const SmartPointer<const FiniteElement<dim,spacedim>,FEValuesBase<dim,spacedim> > fe;
+  std_cxx11::unique_ptr<typename Mapping<dim,spacedim>::InternalDataBase> mapping_data;
 
   /**
-   * Internal data of mapping.
+   * An object into which the Mapping::fill_fe_values() and similar
+   * functions place their output.
    */
-  std_cxx11::unique_ptr<typename Mapping<dim,spacedim>::InternalDataBase> mapping_data;
+  dealii::internal::FEValues::MappingRelatedData<dim, spacedim> mapping_output;
+
 
   /**
-   * Internal data of finite element.
+   * A pointer to the finite element object associated with this FEValues object.
+   */
+  const SmartPointer<const FiniteElement<dim,spacedim>,FEValuesBase<dim,spacedim> > fe;
+
+  /**
+   * A pointer to the internal data object of finite element, obtained from
+   * FiniteElement::get_data(), Mapping::get_face_data(), or
+   * FiniteElement::get_subface_data().
    */
   std_cxx11::unique_ptr<typename FiniteElement<dim,spacedim>::InternalDataBase> fe_data;
 
+  /**
+   * An object into which the FiniteElement::fill_fe_values() and similar
+   * functions place their output.
+   */
+  dealii::internal::FEValues::FiniteElementRelatedData<dim, spacedim> finite_element_output;
+
+
   /**
    * Original update flags handed to the constructor of FEValues.
    */
@@ -2871,9 +2895,9 @@ namespace FEValuesViews
     // except that here we know the component as fixed and we have
     // pre-computed and cached a bunch of information. See the comments there.
     if (shape_function_data[shape_function].is_nonzero_shape_function_component)
-      return fe_values.shape_values(shape_function_data[shape_function]
-                                    .row_index,
-                                    q_point);
+      return fe_values.finite_element_output.shape_values(shape_function_data[shape_function]
+                                                          .row_index,
+                                                          q_point);
     else
       return 0;
   }
@@ -2900,8 +2924,8 @@ namespace FEValuesViews
     // pre-computed and cached a bunch of
     // information. See the comments there.
     if (shape_function_data[shape_function].is_nonzero_shape_function_component)
-      return fe_values.shape_gradients[shape_function_data[shape_function]
-                                       .row_index][q_point];
+      return fe_values.finite_element_output.shape_gradients[shape_function_data[shape_function]
+                                                             .row_index][q_point];
     else
       return gradient_type();
   }
@@ -2927,7 +2951,7 @@ namespace FEValuesViews
     // pre-computed and cached a bunch of
     // information. See the comments there.
     if (shape_function_data[shape_function].is_nonzero_shape_function_component)
-      return fe_values.shape_hessians[shape_function_data[shape_function].row_index][q_point];
+      return fe_values.finite_element_output.shape_hessians[shape_function_data[shape_function].row_index][q_point];
     else
       return hessian_type();
   }
@@ -2955,7 +2979,7 @@ namespace FEValuesViews
       {
         value_type return_value;
         return_value[shape_function_data[shape_function].single_nonzero_component_index]
-          = fe_values.shape_values(snc,q_point);
+          = fe_values.finite_element_output.shape_values(snc,q_point);
         return return_value;
       }
     else
@@ -2964,7 +2988,7 @@ namespace FEValuesViews
         for (unsigned int d=0; d<dim; ++d)
           if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
             return_value[d]
-              = fe_values.shape_values(shape_function_data[shape_function].row_index[d],q_point);
+              = fe_values.finite_element_output.shape_values(shape_function_data[shape_function].row_index[d],q_point);
 
         return return_value;
       }
@@ -2993,7 +3017,7 @@ namespace FEValuesViews
       {
         gradient_type return_value;
         return_value[shape_function_data[shape_function].single_nonzero_component_index]
-          = fe_values.shape_gradients[snc][q_point];
+          = fe_values.finite_element_output.shape_gradients[snc][q_point];
         return return_value;
       }
     else
@@ -3002,7 +3026,7 @@ namespace FEValuesViews
         for (unsigned int d=0; d<dim; ++d)
           if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
             return_value[d]
-              = fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point];
+              = fe_values.finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point];
 
         return return_value;
       }
@@ -3031,14 +3055,14 @@ namespace FEValuesViews
       return divergence_type();
     else if (snc != -1)
       return
-        fe_values.shape_gradients[snc][q_point][shape_function_data[shape_function].single_nonzero_component_index];
+        fe_values.finite_element_output.shape_gradients[snc][q_point][shape_function_data[shape_function].single_nonzero_component_index];
     else
       {
         divergence_type return_value = 0;
         for (unsigned int d=0; d<dim; ++d)
           if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
             return_value
-            += fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point][d];
+            += fe_values.finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point][d];
 
         return return_value;
       }
@@ -3084,9 +3108,9 @@ namespace FEValuesViews
               // can only be zero
               // or one in 2d
               if (shape_function_data[shape_function].single_nonzero_component_index == 0)
-                return_value[0] = -1.0 * fe_values.shape_gradients[snc][q_point][1];
+                return_value[0] = -1.0 * fe_values.finite_element_output.shape_gradients[snc][q_point][1];
               else
-                return_value[0] = fe_values.shape_gradients[snc][q_point][0];
+                return_value[0] = fe_values.finite_element_output.shape_gradients[snc][q_point][0];
 
               return return_value;
             }
@@ -3099,11 +3123,11 @@ namespace FEValuesViews
 
               if (shape_function_data[shape_function].is_nonzero_shape_function_component[0])
                 return_value[0]
-                -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][1];
+                -= fe_values.finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][1];
 
               if (shape_function_data[shape_function].is_nonzero_shape_function_component[1])
                 return_value[0]
-                += fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][0];
+                += fe_values.finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][0];
 
               return return_value;
             }
@@ -3120,23 +3144,23 @@ namespace FEValuesViews
                 case 0:
                 {
                   return_value[0] = 0;
-                  return_value[1] = fe_values.shape_gradients[snc][q_point][2];
-                  return_value[2] = -1.0 * fe_values.shape_gradients[snc][q_point][1];
+                  return_value[1] = fe_values.finite_element_output.shape_gradients[snc][q_point][2];
+                  return_value[2] = -1.0 * fe_values.finite_element_output.shape_gradients[snc][q_point][1];
                   return return_value;
                 }
 
                 case 1:
                 {
-                  return_value[0] = -1.0 * fe_values.shape_gradients[snc][q_point][2];
+                  return_value[0] = -1.0 * fe_values.finite_element_output.shape_gradients[snc][q_point][2];
                   return_value[1] = 0;
-                  return_value[2] = fe_values.shape_gradients[snc][q_point][0];
+                  return_value[2] = fe_values.finite_element_output.shape_gradients[snc][q_point][0];
                   return return_value;
                 }
 
                 default:
                 {
-                  return_value[0] = fe_values.shape_gradients[snc][q_point][1];
-                  return_value[1] = -1.0 * fe_values.shape_gradients[snc][q_point][0];
+                  return_value[0] = fe_values.finite_element_output.shape_gradients[snc][q_point][1];
+                  return_value[1] = -1.0 * fe_values.finite_element_output.shape_gradients[snc][q_point][0];
                   return_value[2] = 0;
                   return return_value;
                 }
@@ -3153,25 +3177,25 @@ namespace FEValuesViews
               if (shape_function_data[shape_function].is_nonzero_shape_function_component[0])
                 {
                   return_value[1]
-                  += fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][2];
+                  += fe_values.finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][2];
                   return_value[2]
-                  -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][1];
+                  -= fe_values.finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][1];
                 }
 
               if (shape_function_data[shape_function].is_nonzero_shape_function_component[1])
                 {
                   return_value[0]
-                  -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][2];
+                  -= fe_values.finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][2];
                   return_value[2]
-                  += fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][0];
+                  += fe_values.finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][0];
                 }
 
               if (shape_function_data[shape_function].is_nonzero_shape_function_component[2])
                 {
                   return_value[0]
-                  += fe_values.shape_gradients[shape_function_data[shape_function].row_index[2]][q_point][1];
+                  += fe_values.finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[2]][q_point][1];
                   return_value[1]
-                  -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[2]][q_point][0];
+                  -= fe_values.finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[2]][q_point][0];
                 }
 
               return return_value;
@@ -3206,7 +3230,7 @@ namespace FEValuesViews
       {
         hessian_type return_value;
         return_value[shape_function_data[shape_function].single_nonzero_component_index]
-          = fe_values.shape_hessians[snc][q_point];
+          = fe_values.finite_element_output.shape_hessians[snc][q_point];
         return return_value;
       }
     else
@@ -3215,7 +3239,7 @@ namespace FEValuesViews
         for (unsigned int d=0; d<dim; ++d)
           if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
             return_value[d]
-              = fe_values.shape_hessians[shape_function_data[shape_function].row_index[d]][q_point];
+              = fe_values.finite_element_output.shape_hessians[shape_function_data[shape_function].row_index[d]][q_point];
 
         return return_value;
       }
@@ -3318,14 +3342,14 @@ namespace FEValuesViews
       return symmetric_gradient_type();
     else if (snc != -1)
       return symmetrize_single_row (shape_function_data[shape_function].single_nonzero_component_index,
-                                    fe_values.shape_gradients[snc][q_point]);
+                                    fe_values.finite_element_output.shape_gradients[snc][q_point]);
     else
       {
         gradient_type return_value;
         for (unsigned int d=0; d<dim; ++d)
           if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
             return_value[d]
-              = fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point];
+              = fe_values.finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point];
 
         return symmetrize(return_value);
       }
@@ -3365,7 +3389,7 @@ namespace FEValuesViews
         const unsigned int comp =
           shape_function_data[shape_function].single_nonzero_component_index;
         return_value[value_type::unrolled_to_component_indices(comp)]
-          = fe_values.shape_values(snc,q_point);
+          = fe_values.finite_element_output.shape_values(snc,q_point);
         return return_value;
       }
     else
@@ -3374,7 +3398,7 @@ namespace FEValuesViews
         for (unsigned int d = 0; d < value_type::n_independent_components; ++d)
           if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
             return_value[value_type::unrolled_to_component_indices(d)]
-              = fe_values.shape_values(shape_function_data[shape_function].row_index[d],q_point);
+              = fe_values.finite_element_output.shape_values(shape_function_data[shape_function].row_index[d],q_point);
         return return_value;
       }
   }
@@ -3451,7 +3475,7 @@ namespace FEValuesViews
         // b_jj := \dfrac{\partial phi_{ii,jj}}{\partial x_jj}.
         // again, all other entries of 'b' are
         // zero
-        const dealii::Tensor<1, spacedim> phi_grad = fe_values.shape_gradients[snc][q_point];
+        const dealii::Tensor<1, spacedim> phi_grad = fe_values.finite_element_output.shape_gradients[snc][q_point];
 
         divergence_type return_value;
         return_value[ii] = phi_grad[jj];
@@ -3502,7 +3526,7 @@ namespace FEValuesViews
         const unsigned int comp =
           shape_function_data[shape_function].single_nonzero_component_index;
         const TableIndices<2> indices = dealii::Tensor<2,spacedim>::unrolled_to_component_indices(comp);
-        return_value[indices] = fe_values.shape_values(snc,q_point);
+        return_value[indices] = fe_values.finite_element_output.shape_values(snc,q_point);
         return return_value;
       }
     else
@@ -3513,7 +3537,7 @@ namespace FEValuesViews
             {
               const TableIndices<2> indices = dealii::Tensor<2,spacedim>::unrolled_to_component_indices(d);
               return_value[indices]
-                = fe_values.shape_values(shape_function_data[shape_function].row_index[d],q_point);
+                = fe_values.finite_element_output.shape_values(shape_function_data[shape_function].row_index[d],q_point);
             }
         return return_value;
       }
@@ -3564,7 +3588,7 @@ namespace FEValuesViews
         const unsigned int ii = indices[0];
         const unsigned int jj = indices[1];
 
-        const dealii::Tensor<1, spacedim> phi_grad = fe_values.shape_gradients[snc][q_point];
+        const dealii::Tensor<1, spacedim> phi_grad = fe_values.finite_element_output.shape_gradients[snc][q_point];
 
         divergence_type return_value;
         return_value[jj] = phi_grad[ii];
@@ -3663,7 +3687,7 @@ FEValuesBase<dim,spacedim>::shape_value (const unsigned int i,
   // if the entire FE is primitive,
   // then we can take a short-cut:
   if (fe->is_primitive())
-    return this->shape_values(i,j);
+    return this->finite_element_output.shape_values(i,j);
   else
     {
       // otherwise, use the mapping
@@ -3675,8 +3699,8 @@ FEValuesBase<dim,spacedim>::shape_value (const unsigned int i,
       // so we can call
       // system_to_component_index
       const unsigned int
-      row = this->shape_function_to_row_table[i * fe->n_components() + fe->system_to_component_index(i).first];
-      return this->shape_values(row, j);
+      row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + fe->system_to_component_index(i).first];
+      return this->finite_element_output.shape_values(row, j);
     }
 }
 
@@ -3706,8 +3730,8 @@ FEValuesBase<dim,spacedim>::shape_value_component (const unsigned int i,
   // table and take the data from
   // there
   const unsigned int
-  row = this->shape_function_to_row_table[i * fe->n_components() + component];
-  return this->shape_values(row, j);
+  row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + component];
+  return this->finite_element_output.shape_values(row, j);
 }
 
 
@@ -3724,15 +3748,15 @@ FEValuesBase<dim,spacedim>::shape_grad (const unsigned int i,
           ExcAccessToUninitializedField("update_gradients"));
   Assert (fe->is_primitive (i),
           ExcShapeFunctionNotPrimitive(i));
-  Assert (i<this->shape_gradients.size(),
-          ExcIndexRange (i, 0, this->shape_gradients.size()));
-  Assert (j<this->shape_gradients[0].size(),
-          ExcIndexRange (j, 0, this->shape_gradients[0].size()));
+  Assert (i<this->finite_element_output.shape_gradients.size(),
+          ExcIndexRange (i, 0, this->finite_element_output.shape_gradients.size()));
+  Assert (j<this->finite_element_output.shape_gradients[0].size(),
+          ExcIndexRange (j, 0, this->finite_element_output.shape_gradients[0].size()));
 
   // if the entire FE is primitive,
   // then we can take a short-cut:
   if (fe->is_primitive())
-    return this->shape_gradients[i][j];
+    return this->finite_element_output.shape_gradients[i][j];
   else
     {
       // otherwise, use the mapping
@@ -3744,8 +3768,8 @@ FEValuesBase<dim,spacedim>::shape_grad (const unsigned int i,
       // so we can call
       // system_to_component_index
       const unsigned int
-      row = this->shape_function_to_row_table[i * fe->n_components() + fe->system_to_component_index(i).first];
-      return this->shape_gradients[row][j];
+      row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + fe->system_to_component_index(i).first];
+      return this->finite_element_output.shape_gradients[row][j];
     }
 }
 
@@ -3775,8 +3799,8 @@ FEValuesBase<dim,spacedim>::shape_grad_component (const unsigned int i,
   // table and take the data from
   // there
   const unsigned int
-  row = this->shape_function_to_row_table[i * fe->n_components() + component];
-  return this->shape_gradients[row][j];
+  row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + component];
+  return this->finite_element_output.shape_gradients[row][j];
 }
 
 
@@ -3793,15 +3817,15 @@ FEValuesBase<dim,spacedim>::shape_hessian (const unsigned int i,
           ExcAccessToUninitializedField("update_hessians"));
   Assert (fe->is_primitive (i),
           ExcShapeFunctionNotPrimitive(i));
-  Assert (i<this->shape_hessians.size(),
-          ExcIndexRange (i, 0, this->shape_hessians.size()));
-  Assert (j<this->shape_hessians[0].size(),
-          ExcIndexRange (j, 0, this->shape_hessians[0].size()));
+  Assert (i<this->finite_element_output.shape_hessians.size(),
+          ExcIndexRange (i, 0, this->finite_element_output.shape_hessians.size()));
+  Assert (j<this->finite_element_output.shape_hessians[0].size(),
+          ExcIndexRange (j, 0, this->finite_element_output.shape_hessians[0].size()));
 
   // if the entire FE is primitive,
   // then we can take a short-cut:
   if (fe->is_primitive())
-    return this->shape_hessians[i][j];
+    return this->finite_element_output.shape_hessians[i][j];
   else
     {
       // otherwise, use the mapping
@@ -3813,8 +3837,8 @@ FEValuesBase<dim,spacedim>::shape_hessian (const unsigned int i,
       // so we can call
       // system_to_component_index
       const unsigned int
-      row = this->shape_function_to_row_table[i * fe->n_components() + fe->system_to_component_index(i).first];
-      return this->shape_hessians[row][j];
+      row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + fe->system_to_component_index(i).first];
+      return this->finite_element_output.shape_hessians[row][j];
     }
 }
 
@@ -3844,8 +3868,8 @@ FEValuesBase<dim,spacedim>::shape_hessian_component (const unsigned int i,
   // table and take the data from
   // there
   const unsigned int
-  row = this->shape_function_to_row_table[i * fe->n_components() + component];
-  return this->shape_hessians[row][j];
+  row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + component];
+  return this->finite_element_output.shape_hessians[row][j];
 }
 
 
@@ -3886,7 +3910,7 @@ FEValuesBase<dim,spacedim>::get_quadrature_points () const
 {
   Assert (this->update_flags & update_quadrature_points,
           ExcAccessToUninitializedField("update_quadrature_points"));
-  return this->quadrature_points;
+  return this->mapping_output.quadrature_points;
 }
 
 
@@ -3898,7 +3922,7 @@ FEValuesBase<dim,spacedim>::get_JxW_values () const
 {
   Assert (this->update_flags & update_JxW_values,
           ExcAccessToUninitializedField("update_JxW_values"));
-  return this->JxW_values;
+  return this->mapping_output.JxW_values;
 }
 
 
@@ -3910,7 +3934,7 @@ FEValuesBase<dim,spacedim>::get_jacobians () const
 {
   Assert (this->update_flags & update_jacobians,
           ExcAccessToUninitializedField("update_jacobians"));
-  return this->jacobians;
+  return this->mapping_output.jacobians;
 }
 
 
@@ -3922,7 +3946,7 @@ FEValuesBase<dim,spacedim>::get_jacobian_grads () const
 {
   Assert (this->update_flags & update_jacobian_grads,
           ExcAccessToUninitializedField("update_jacobians_grads"));
-  return this->jacobian_grads;
+  return this->mapping_output.jacobian_grads;
 }
 
 
@@ -3934,7 +3958,7 @@ FEValuesBase<dim,spacedim>::get_inverse_jacobians () const
 {
   Assert (this->update_flags & update_inverse_jacobians,
           ExcAccessToUninitializedField("update_inverse_jacobians"));
-  return this->inverse_jacobians;
+  return this->mapping_output.inverse_jacobians;
 }
 
 
@@ -3946,9 +3970,10 @@ FEValuesBase<dim,spacedim>::quadrature_point (const unsigned int i) const
 {
   Assert (this->update_flags & update_quadrature_points,
           ExcAccessToUninitializedField("update_quadrature_points"));
-  Assert (i<this->quadrature_points.size(), ExcIndexRange(i, 0, this->quadrature_points.size()));
+  Assert (i<this->mapping_output.quadrature_points.size(),
+          ExcIndexRange(i, 0, this->mapping_output.quadrature_points.size()));
 
-  return this->quadrature_points[i];
+  return this->mapping_output.quadrature_points[i];
 }
 
 
@@ -3961,9 +3986,10 @@ FEValuesBase<dim,spacedim>::JxW (const unsigned int i) const
 {
   Assert (this->update_flags & update_JxW_values,
           ExcAccessToUninitializedField("update_JxW_values"));
-  Assert (i<this->JxW_values.size(), ExcIndexRange(i, 0, this->JxW_values.size()));
+  Assert (i<this->mapping_output.JxW_values.size(),
+          ExcIndexRange(i, 0, this->mapping_output.JxW_values.size()));
 
-  return this->JxW_values[i];
+  return this->mapping_output.JxW_values[i];
 }
 
 
@@ -3975,9 +4001,10 @@ FEValuesBase<dim,spacedim>::jacobian (const unsigned int i) const
 {
   Assert (this->update_flags & update_jacobians,
           ExcAccessToUninitializedField("update_jacobians"));
-  Assert (i<this->jacobians.size(), ExcIndexRange(i, 0, this->jacobians.size()));
+  Assert (i<this->mapping_output.jacobians.size(),
+          ExcIndexRange(i, 0, this->mapping_output.jacobians.size()));
 
-  return this->jacobians[i];
+  return this->mapping_output.jacobians[i];
 }
 
 
@@ -3989,9 +4016,10 @@ FEValuesBase<dim,spacedim>::jacobian_grad (const unsigned int i) const
 {
   Assert (this->update_flags & update_jacobian_grads,
           ExcAccessToUninitializedField("update_jacobians_grads"));
-  Assert (i<this->jacobian_grads.size(), ExcIndexRange(i, 0, this->jacobian_grads.size()));
+  Assert (i<this->mapping_output.jacobian_grads.size(),
+          ExcIndexRange(i, 0, this->mapping_output.jacobian_grads.size()));
 
-  return this->jacobian_grads[i];
+  return this->mapping_output.jacobian_grads[i];
 }
 
 
@@ -4003,9 +4031,10 @@ FEValuesBase<dim,spacedim>::inverse_jacobian (const unsigned int i) const
 {
   Assert (this->update_flags & update_inverse_jacobians,
           ExcAccessToUninitializedField("update_inverse_jacobians"));
-  Assert (i<this->inverse_jacobians.size(), ExcIndexRange(i, 0, this->inverse_jacobians.size()));
+  Assert (i<this->mapping_output.inverse_jacobians.size(),
+          ExcIndexRange(i, 0, this->mapping_output.inverse_jacobians.size()));
 
-  return this->inverse_jacobians[i];
+  return this->mapping_output.inverse_jacobians[i];
 }
 
 
@@ -4017,10 +4046,10 @@ FEValuesBase<dim,spacedim>::normal_vector (const unsigned int i) const
   typedef FEValuesBase<dim,spacedim> FVB;
   Assert (this->update_flags & update_normal_vectors,
           typename FVB::ExcAccessToUninitializedField("update_normal_vectors"));
-  Assert (i<this->normal_vectors.size(),
-          ExcIndexRange(i, 0, this->normal_vectors.size()));
+  Assert (i<this->mapping_output.normal_vectors.size(),
+          ExcIndexRange(i, 0, this->mapping_output.normal_vectors.size()));
 
-  return this->normal_vectors[i];
+  return this->mapping_output.normal_vectors[i];
 }
 
 
@@ -4097,12 +4126,12 @@ const Tensor<1,spacedim> &
 FEFaceValuesBase<dim,spacedim>::boundary_form (const unsigned int i) const
 {
   typedef FEValuesBase<dim,spacedim> FVB;
-  Assert (i<this->boundary_forms.size(),
-          ExcIndexRange(i, 0, this->boundary_forms.size()));
+  Assert (i<this->mapping_output.boundary_forms.size(),
+          ExcIndexRange(i, 0, this->mapping_output.boundary_forms.size()));
   Assert (this->update_flags & update_boundary_forms,
           typename FVB::ExcAccessToUninitializedField("update_boundary_forms"));
 
-  return this->boundary_forms[i];
+  return this->mapping_output.boundary_forms[i];
 }
 
 #endif // DOXYGEN
diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc
index 7003914a6d..e9916001f2 100644
--- a/source/fe/fe_values.cc
+++ b/source/fe/fe_values.cc
@@ -1279,7 +1279,7 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values(fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
     internal::do_function_values<dim,spacedim>
-    (dof_values, fe_values.shape_values, shape_function_data, values);
+    (dof_values, fe_values.finite_element_output.shape_values, shape_function_data, values);
   }
 
 
@@ -1303,7 +1303,7 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
     internal::do_function_derivatives<1,dim,spacedim>
-    (dof_values, fe_values.shape_gradients, shape_function_data, gradients);
+    (dof_values, fe_values.finite_element_output.shape_gradients, shape_function_data, gradients);
   }
 
 
@@ -1327,7 +1327,7 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
     internal::do_function_derivatives<2,dim,spacedim>
-    (dof_values, fe_values.shape_hessians, shape_function_data, hessians);
+    (dof_values, fe_values.finite_element_output.shape_hessians, shape_function_data, hessians);
   }
 
 
@@ -1351,7 +1351,7 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
     internal::do_function_laplacians<dim,spacedim>
-    (dof_values, fe_values.shape_hessians, shape_function_data, laplacians);
+    (dof_values, fe_values.finite_element_output.shape_hessians, shape_function_data, laplacians);
   }
 
 
@@ -1375,7 +1375,7 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
     internal::do_function_values<dim,spacedim>
-    (dof_values, fe_values.shape_values, shape_function_data, values);
+    (dof_values, fe_values.finite_element_output.shape_values, shape_function_data, values);
   }
 
 
@@ -1400,7 +1400,7 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
     internal::do_function_derivatives<1,dim,spacedim>
-    (dof_values, fe_values.shape_gradients, shape_function_data, gradients);
+    (dof_values, fe_values.finite_element_output.shape_gradients, shape_function_data, gradients);
   }
 
 
@@ -1424,7 +1424,7 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
     internal::do_function_symmetric_gradients<dim,spacedim>
-    (dof_values, fe_values.shape_gradients, shape_function_data,
+    (dof_values, fe_values.finite_element_output.shape_gradients, shape_function_data,
      symmetric_gradients);
   }
 
@@ -1450,7 +1450,7 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
     internal::do_function_divergences<dim,spacedim>
-    (dof_values, fe_values.shape_gradients, shape_function_data, divergences);
+    (dof_values, fe_values.finite_element_output.shape_gradients, shape_function_data, divergences);
   }
 
   template <int dim, int spacedim>
@@ -1473,7 +1473,7 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values (fe_function, dof_values);
     internal::do_function_curls<dim,spacedim>
-    (dof_values, fe_values.shape_gradients, shape_function_data, curls);
+    (dof_values, fe_values.finite_element_output.shape_gradients, shape_function_data, curls);
   }
 
 
@@ -1496,7 +1496,7 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
     internal::do_function_derivatives<2,dim,spacedim>
-    (dof_values, fe_values.shape_hessians, shape_function_data, hessians);
+    (dof_values, fe_values.finite_element_output.shape_hessians, shape_function_data, hessians);
   }
 
 
@@ -1523,7 +1523,7 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
     internal::do_function_laplacians<dim,spacedim>
-    (dof_values, fe_values.shape_hessians, shape_function_data, laplacians);
+    (dof_values, fe_values.finite_element_output.shape_hessians, shape_function_data, laplacians);
   }
 
 
@@ -1547,7 +1547,7 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values(fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
     internal::do_function_values<dim,spacedim>
-    (dof_values, fe_values.shape_values, shape_function_data, values);
+    (dof_values, fe_values.finite_element_output.shape_values, shape_function_data, values);
   }
 
 
@@ -1572,7 +1572,7 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values(fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
     internal::do_function_divergences<dim,spacedim>
-    (dof_values, fe_values.shape_gradients, shape_function_data, divergences);
+    (dof_values, fe_values.finite_element_output.shape_gradients, shape_function_data, divergences);
   }
 
   template <int dim, int spacedim>
@@ -1594,7 +1594,7 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values(fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
     internal::do_function_values<dim,spacedim>
-    (dof_values, fe_values.shape_values, shape_function_data, values);
+    (dof_values, fe_values.finite_element_output.shape_values, shape_function_data, values);
   }
 
 
@@ -1619,7 +1619,7 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values(fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
     internal::do_function_divergences<dim,spacedim>
-    (dof_values, fe_values.shape_gradients, shape_function_data, divergences);
+    (dof_values, fe_values.finite_element_output.shape_gradients, shape_function_data, divergences);
   }
 }
 
@@ -2102,6 +2102,23 @@ namespace internal
     }
 
 
+
+    template <int dim, int spacedim>
+    std::size_t
+    MappingRelatedData<dim,spacedim>::memory_consumption () const
+    {
+      return (MemoryConsumption::memory_consumption (JxW_values) +
+              MemoryConsumption::memory_consumption (jacobians) +
+              MemoryConsumption::memory_consumption (jacobian_grads) +
+              MemoryConsumption::memory_consumption (inverse_jacobians) +
+              MemoryConsumption::memory_consumption (quadrature_points) +
+              MemoryConsumption::memory_consumption (normal_vectors) +
+              MemoryConsumption::memory_consumption (boundary_forms));
+    }
+
+
+
+
     template <int dim, int spacedim>
     void
     FiniteElementRelatedData<dim,spacedim>::initialize (const unsigned int        n_quadrature_points,
@@ -2142,6 +2159,19 @@ namespace internal
         this->shape_hessians.resize (n_nonzero_shape_components,
                                      std::vector<Tensor<2,spacedim> > (n_quadrature_points));
     }
+
+
+
+
+    template <int dim, int spacedim>
+    std::size_t
+    FiniteElementRelatedData<dim,spacedim>::memory_consumption () const
+    {
+      return (MemoryConsumption::memory_consumption (shape_values) +
+              MemoryConsumption::memory_consumption (shape_gradients) +
+              MemoryConsumption::memory_consumption (shape_hessians) +
+              MemoryConsumption::memory_consumption (shape_function_to_row_table));
+    }
   }
 }
 
@@ -2601,7 +2631,7 @@ void FEValuesBase<dim,spacedim>::get_function_values (
   // get function values of dofs on this cell
   Vector<Number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
-  internal::do_function_values (dof_values.begin(), this->shape_values,
+  internal::do_function_values (dof_values.begin(), this->finite_element_output.shape_values,
                                 values);
 }
 
@@ -2626,14 +2656,14 @@ void FEValuesBase<dim,spacedim>::get_function_values (
       Number dof_values[100];
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_values(&dof_values[0], this->shape_values, values);
+      internal::do_function_values(&dof_values[0], this->finite_element_output.shape_values, values);
     }
   else
     {
       Vector<Number> dof_values(dofs_per_cell);
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_values(dof_values.begin(), this->shape_values,
+      internal::do_function_values(dof_values.begin(), this->finite_element_output.shape_values,
                                    values);
     }
 }
@@ -2658,8 +2688,8 @@ void FEValuesBase<dim,spacedim>::get_function_values (
   Vector<Number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
   VectorSlice<std::vector<Vector<Number> > > val(values);
-  internal::do_function_values(dof_values.begin(), this->shape_values, *fe,
-                               this->shape_function_to_row_table, val);
+  internal::do_function_values(dof_values.begin(), this->finite_element_output.shape_values, *fe,
+                               this->finite_element_output.shape_function_to_row_table, val);
 }
 
 
@@ -2685,8 +2715,8 @@ void FEValuesBase<dim,spacedim>::get_function_values (
       Number dof_values[100];
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_values(&dof_values[0], this->shape_values, *fe,
-                                   this->shape_function_to_row_table, val,
+      internal::do_function_values(&dof_values[0], this->finite_element_output.shape_values, *fe,
+                                   this->finite_element_output.shape_function_to_row_table, val,
                                    false, indices.size()/dofs_per_cell);
     }
   else
@@ -2694,8 +2724,8 @@ void FEValuesBase<dim,spacedim>::get_function_values (
       Vector<Number> dof_values(100);
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_values(dof_values.begin(), this->shape_values, *fe,
-                                   this->shape_function_to_row_table, val,
+      internal::do_function_values(dof_values.begin(), this->finite_element_output.shape_values, *fe,
+                                   this->finite_element_output.shape_function_to_row_table, val,
                                    false, indices.size()/dofs_per_cell);
     }
 }
@@ -2724,8 +2754,8 @@ void FEValuesBase<dim,spacedim>::get_function_values (
       Number dof_values[100];
       for (unsigned int i=0; i<indices.size(); ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_values(&dof_values[0], this->shape_values, *fe,
-                                   this->shape_function_to_row_table, values,
+      internal::do_function_values(&dof_values[0], this->finite_element_output.shape_values, *fe,
+                                   this->finite_element_output.shape_function_to_row_table, values,
                                    quadrature_points_fastest,
                                    indices.size()/dofs_per_cell);
     }
@@ -2734,8 +2764,8 @@ void FEValuesBase<dim,spacedim>::get_function_values (
       Vector<Number> dof_values(indices.size());
       for (unsigned int i=0; i<indices.size(); ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_values(dof_values.begin(), this->shape_values, *fe,
-                                   this->shape_function_to_row_table, values,
+      internal::do_function_values(dof_values.begin(), this->finite_element_output.shape_values, *fe,
+                                   this->finite_element_output.shape_function_to_row_table, values,
                                    quadrature_points_fastest,
                                    indices.size()/dofs_per_cell);
     }
@@ -2761,7 +2791,7 @@ FEValuesBase<dim,spacedim>::get_function_gradients (
   // get function values of dofs on this cell
   Vector<Number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
-  internal::do_function_derivatives(dof_values.begin(), this->shape_gradients,
+  internal::do_function_derivatives(dof_values.begin(), this->finite_element_output.shape_gradients,
                                     gradients);
 }
 
@@ -2784,7 +2814,7 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
       Number dof_values[100];
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_derivatives(&dof_values[0], this->shape_gradients,
+      internal::do_function_derivatives(&dof_values[0], this->finite_element_output.shape_gradients,
                                         gradients);
     }
   else
@@ -2792,7 +2822,7 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
       Vector<Number> dof_values(dofs_per_cell);
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_derivatives(dof_values.begin(), this->shape_gradients,
+      internal::do_function_derivatives(dof_values.begin(), this->finite_element_output.shape_gradients,
                                         gradients);
     }
 }
@@ -2818,8 +2848,8 @@ FEValuesBase<dim,spacedim>::get_function_gradients (
   Vector<Number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
   VectorSlice<std::vector<std::vector<Tensor<1,spacedim,Number> > > > grads(gradients);
-  internal::do_function_derivatives(dof_values.begin(), this->shape_gradients,
-                                    *fe, this->shape_function_to_row_table,
+  internal::do_function_derivatives(dof_values.begin(), this->finite_element_output.shape_gradients,
+                                    *fe, this->finite_element_output.shape_function_to_row_table,
                                     grads);
 }
 
@@ -2846,8 +2876,8 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
       Number dof_values[100];
       for (unsigned int i=0; i<indices.size(); ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_derivatives(&dof_values[0], this->shape_gradients,
-                                        *fe, this->shape_function_to_row_table,
+      internal::do_function_derivatives(&dof_values[0], this->finite_element_output.shape_gradients,
+                                        *fe, this->finite_element_output.shape_function_to_row_table,
                                         gradients, quadrature_points_fastest,
                                         indices.size()/dofs_per_cell);
     }
@@ -2856,8 +2886,8 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
       Vector<Number> dof_values(indices.size());
       for (unsigned int i=0; i<indices.size(); ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_derivatives(dof_values.begin(),this->shape_gradients,
-                                        *fe, this->shape_function_to_row_table,
+      internal::do_function_derivatives(dof_values.begin(),this->finite_element_output.shape_gradients,
+                                        *fe, this->finite_element_output.shape_function_to_row_table,
                                         gradients, quadrature_points_fastest,
                                         indices.size()/dofs_per_cell);
     }
@@ -2883,7 +2913,7 @@ get_function_hessians (const InputVector                &fe_function,
   // get function values of dofs on this cell
   Vector<Number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
-  internal::do_function_derivatives(dof_values.begin(), this->shape_hessians,
+  internal::do_function_derivatives(dof_values.begin(), this->finite_element_output.shape_hessians,
                                     hessians);
 }
 
@@ -2906,7 +2936,7 @@ void FEValuesBase<dim,spacedim>::get_function_hessians (
       Number dof_values[100];
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_derivatives(&dof_values[0], this->shape_hessians,
+      internal::do_function_derivatives(&dof_values[0], this->finite_element_output.shape_hessians,
                                         hessians);
     }
   else
@@ -2914,7 +2944,7 @@ void FEValuesBase<dim,spacedim>::get_function_hessians (
       Vector<Number> dof_values(dofs_per_cell);
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_derivatives(dof_values.begin(), this->shape_hessians,
+      internal::do_function_derivatives(dof_values.begin(), this->finite_element_output.shape_hessians,
                                         hessians);
     }
 }
@@ -2941,8 +2971,8 @@ get_function_hessians (const InputVector                         &fe_function,
   Vector<Number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
   VectorSlice<std::vector<std::vector<Tensor<2,spacedim,Number> > > > hes(hessians);
-  internal::do_function_derivatives(dof_values.begin(), this->shape_hessians,
-                                    *fe, this->shape_function_to_row_table,
+  internal::do_function_derivatives(dof_values.begin(), this->finite_element_output.shape_hessians,
+                                    *fe, this->finite_element_output.shape_function_to_row_table,
                                     hes, quadrature_points_fastest);
 }
 
@@ -2966,8 +2996,8 @@ void FEValuesBase<dim, spacedim>::get_function_hessians (
       Number dof_values[100];
       for (unsigned int i=0; i<indices.size(); ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_derivatives(&dof_values[0], this->shape_hessians,
-                                        *fe, this->shape_function_to_row_table,
+      internal::do_function_derivatives(&dof_values[0], this->finite_element_output.shape_hessians,
+                                        *fe, this->finite_element_output.shape_function_to_row_table,
                                         hessians, quadrature_points_fastest,
                                         indices.size()/dofs_per_cell);
     }
@@ -2976,8 +3006,8 @@ void FEValuesBase<dim, spacedim>::get_function_hessians (
       Vector<Number> dof_values(indices.size());
       for (unsigned int i=0; i<indices.size(); ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_derivatives(dof_values.begin(),this->shape_hessians,
-                                        *fe, this->shape_function_to_row_table,
+      internal::do_function_derivatives(dof_values.begin(),this->finite_element_output.shape_hessians,
+                                        *fe, this->finite_element_output.shape_function_to_row_table,
                                         hessians, quadrature_points_fastest,
                                         indices.size()/dofs_per_cell);
     }
@@ -3002,7 +3032,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
   // get function values of dofs on this cell
   Vector<Number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
-  internal::do_function_laplacians(dof_values.begin(), this->shape_hessians,
+  internal::do_function_laplacians(dof_values.begin(), this->finite_element_output.shape_hessians,
                                    laplacians);
 }
 
@@ -3025,7 +3055,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
       Number dof_values[100];
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_laplacians(&dof_values[0], this->shape_hessians,
+      internal::do_function_laplacians(&dof_values[0], this->finite_element_output.shape_hessians,
                                        laplacians);
     }
   else
@@ -3033,7 +3063,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
       Vector<Number> dof_values(dofs_per_cell);
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_laplacians(dof_values.begin(), this->shape_hessians,
+      internal::do_function_laplacians(dof_values.begin(), this->finite_element_output.shape_hessians,
                                        laplacians);
     }
 }
@@ -3056,8 +3086,8 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
   // get function values of dofs on this cell
   Vector<Number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
-  internal::do_function_laplacians(dof_values.begin(), this->shape_hessians,
-                                   *fe, this->shape_function_to_row_table,
+  internal::do_function_laplacians(dof_values.begin(), this->finite_element_output.shape_hessians,
+                                   *fe, this->finite_element_output.shape_function_to_row_table,
                                    laplacians);
 }
 
@@ -3082,8 +3112,8 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
       Number dof_values[100];
       for (unsigned int i=0; i<indices.size(); ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_laplacians(&dof_values[0], this->shape_hessians,
-                                       *fe, this->shape_function_to_row_table,
+      internal::do_function_laplacians(&dof_values[0], this->finite_element_output.shape_hessians,
+                                       *fe, this->finite_element_output.shape_function_to_row_table,
                                        laplacians, false,
                                        indices.size()/dofs_per_cell);
     }
@@ -3092,8 +3122,8 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
       Vector<Number> dof_values(indices.size());
       for (unsigned int i=0; i<indices.size(); ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_laplacians(dof_values.begin(),this->shape_hessians,
-                                       *fe, this->shape_function_to_row_table,
+      internal::do_function_laplacians(dof_values.begin(),this->finite_element_output.shape_hessians,
+                                       *fe, this->finite_element_output.shape_function_to_row_table,
                                        laplacians, false,
                                        indices.size()/dofs_per_cell);
     }
@@ -3119,8 +3149,8 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
       Number dof_values[100];
       for (unsigned int i=0; i<indices.size(); ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_laplacians(&dof_values[0], this->shape_hessians,
-                                       *fe, this->shape_function_to_row_table,
+      internal::do_function_laplacians(&dof_values[0], this->finite_element_output.shape_hessians,
+                                       *fe, this->finite_element_output.shape_function_to_row_table,
                                        laplacians, quadrature_points_fastest,
                                        indices.size()/dofs_per_cell);
     }
@@ -3129,8 +3159,8 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
       Vector<Number> dof_values(indices.size());
       for (unsigned int i=0; i<indices.size(); ++i)
         dof_values[i] = get_vector_element (fe_function, indices[i]);
-      internal::do_function_laplacians(dof_values.begin(),this->shape_hessians,
-                                       *fe, this->shape_function_to_row_table,
+      internal::do_function_laplacians(dof_values.begin(),this->finite_element_output.shape_hessians,
+                                       *fe, this->finite_element_output.shape_function_to_row_table,
                                        laplacians, quadrature_points_fastest,
                                        indices.size()/dofs_per_cell);
     }
@@ -3154,7 +3184,7 @@ FEValuesBase<dim,spacedim>::get_all_normal_vectors () const
   typedef FEValuesBase<dim,spacedim> FEVB;
   Assert (this->update_flags & update_normal_vectors,
           typename FEVB::ExcAccessToUninitializedField("update_normal_vectors"));
-  return this->normal_vectors;
+  return this->mapping_output.normal_vectors;
 }
 
 
@@ -3168,9 +3198,9 @@ FEValuesBase<dim,spacedim>::get_normal_vectors () const
           typename FEVB::ExcAccessToUninitializedField("update_normal_vectors"));
 
   // copy things into a vector of Points, then return that
-  std::vector<Point<spacedim> > tmp (this->normal_vectors.size());
-  for (unsigned int q=0; q<this->normal_vectors.size(); ++q)
-    tmp[q] = Point<spacedim>(this->normal_vectors[q]);
+  std::vector<Point<spacedim> > tmp (this->mapping_output.normal_vectors.size());
+  for (unsigned int q=0; q<this->mapping_output.normal_vectors.size(); ++q)
+    tmp[q] = Point<spacedim>(this->mapping_output.normal_vectors[q]);
 
   return tmp;
 }
@@ -3194,26 +3224,18 @@ template <int dim, int spacedim>
 std::size_t
 FEValuesBase<dim,spacedim>::memory_consumption () const
 {
-  return (MemoryConsumption::memory_consumption (this->shape_values) +
-          MemoryConsumption::memory_consumption (this->shape_gradients) +
-          MemoryConsumption::memory_consumption (this->shape_hessians) +
-          MemoryConsumption::memory_consumption (this->JxW_values) +
-          MemoryConsumption::memory_consumption (this->jacobians) +
-          MemoryConsumption::memory_consumption (this->jacobian_grads) +
-          MemoryConsumption::memory_consumption (this->inverse_jacobians) +
-          MemoryConsumption::memory_consumption (this->quadrature_points) +
-          MemoryConsumption::memory_consumption (this->normal_vectors) +
-          MemoryConsumption::memory_consumption (this->boundary_forms) +
-          sizeof(this->update_flags) +
+  return (sizeof(this->update_flags) +
           MemoryConsumption::memory_consumption (n_quadrature_points) +
+          sizeof (cell_similarity) +
           MemoryConsumption::memory_consumption (dofs_per_cell) +
           MemoryConsumption::memory_consumption (mapping) +
-          MemoryConsumption::memory_consumption (fe) +
           MemoryConsumption::memory_consumption (mapping_data) +
           MemoryConsumption::memory_consumption (*mapping_data) +
+          MemoryConsumption::memory_consumption (mapping_output) +
+          MemoryConsumption::memory_consumption (fe) +
           MemoryConsumption::memory_consumption (fe_data) +
           MemoryConsumption::memory_consumption (*fe_data) +
-          MemoryConsumption::memory_consumption (this->shape_function_to_row_table));
+          MemoryConsumption::memory_consumption (finite_element_output));
 }
 
 
@@ -3434,8 +3456,8 @@ FEValues<dim,spacedim>::initialize (const UpdateFlags update_flags)
                                         quadrature);
 
   // initialize the base classes
-  internal::FEValues::MappingRelatedData<dim,spacedim>::initialize(this->n_quadrature_points, flags);
-  internal::FEValues::FiniteElementRelatedData<dim,spacedim>::initialize(this->n_quadrature_points, *this->fe, flags);
+  this->mapping_output.initialize(this->n_quadrature_points, flags);
+  this->finite_element_output.initialize(this->n_quadrature_points, *this->fe, flags);
 
   this->update_flags = flags;
 
@@ -3540,7 +3562,7 @@ void FEValues<dim,spacedim>::do_reinit ()
                                          this->cell_similarity,
                                          quadrature,
                                          *this->mapping_data,
-                                         *this);
+                                         this->mapping_output);
 
   // then call the finite element and, with the data
   // already filled by the mapping, let it compute the
@@ -3551,8 +3573,8 @@ void FEValues<dim,spacedim>::do_reinit ()
                                 quadrature,
                                 *this->mapping_data,
                                 *this->fe_data,
-                                *this,
-                                *this,
+                                this->mapping_output,
+                                this->finite_element_output,
                                 this->cell_similarity);
 
   this->fe_data->clear_first_cell ();
@@ -3597,7 +3619,7 @@ FEFaceValuesBase<dim,spacedim>::get_boundary_forms () const
   typedef FEValuesBase<dim,spacedim> FEVB;
   Assert (this->update_flags & update_boundary_forms,
           typename FEVB::ExcAccessToUninitializedField("update_boundary_forms"));
-  return this->boundary_forms;
+  return this->mapping_output.boundary_forms;
 }
 
 
@@ -3674,8 +3696,8 @@ FEFaceValues<dim,spacedim>::initialize (const UpdateFlags update_flags)
                                         this->quadrature);
 
   // initialize the base classes
-  internal::FEValues::MappingRelatedData<dim,spacedim>::initialize(this->n_quadrature_points, flags);
-  internal::FEValues::FiniteElementRelatedData<dim,spacedim>::initialize(this->n_quadrature_points, *this->fe, flags);
+  this->mapping_output.initialize(this->n_quadrature_points, flags);
+  this->finite_element_output.initialize(this->n_quadrature_points, *this->fe, flags);
 
   this->update_flags = flags;
 
@@ -3751,15 +3773,15 @@ void FEFaceValues<dim,spacedim>::do_reinit (const unsigned int face_no)
                                           face_no,
                                           this->quadrature,
                                           *this->mapping_data,
-                                          *this);
+                                          this->mapping_output);
 
   this->get_fe().fill_fe_face_values(this->get_mapping(),
                                      *this->present_cell, face_no,
                                      this->quadrature,
                                      *this->mapping_data,
                                      *this->fe_data,
-                                     *this,
-                                     *this);
+                                     this->mapping_output,
+                                     this->finite_element_output);
 
   this->fe_data->clear_first_cell ();
 }
@@ -3831,8 +3853,8 @@ FESubfaceValues<dim,spacedim>::initialize (const UpdateFlags update_flags)
                                         this->quadrature);
 
   // initialize the base classes
-  internal::FEValues::MappingRelatedData<dim,spacedim>::initialize(this->n_quadrature_points, flags);
-  internal::FEValues::FiniteElementRelatedData<dim,spacedim>::initialize(this->n_quadrature_points, *this->fe, flags);
+  this->mapping_output.initialize(this->n_quadrature_points, flags);
+  this->finite_element_output.initialize(this->n_quadrature_points, *this->fe, flags);
 
   this->update_flags = flags;
 
@@ -3993,7 +4015,7 @@ void FESubfaceValues<dim,spacedim>::do_reinit (const unsigned int face_no,
                                              subface_no,
                                              this->quadrature,
                                              *this->mapping_data,
-                                             *this);
+                                             this->mapping_output);
 
   this->get_fe().fill_fe_subface_values(this->get_mapping(),
                                         *this->present_cell,
@@ -4001,8 +4023,8 @@ void FESubfaceValues<dim,spacedim>::do_reinit (const unsigned int face_no,
                                         this->quadrature,
                                         *this->mapping_data,
                                         *this->fe_data,
-                                        *this,
-                                        *this);
+                                        this->mapping_output,
+                                        this->finite_element_output);
 
   this->fe_data->clear_first_cell ();
 }