From bba2c791a07ccfee4e0867d196de9cdaf9ce47e4 Mon Sep 17 00:00:00 2001
From: kronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 19 Feb 2014 09:56:31 +0000
Subject: [PATCH] Make it possible to allocate FEEvaluation on the heap using
 an AlignedVector type. It turns out that the revision 32498 must be reverted
 again: a shared ptr is needed to keep track of the allocated shape info in
 case it is copied around

git-svn-id: https://svn.dealii.org/trunk@32511 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/include/deal.II/base/aligned_vector.h |  62 ++--
 .../deal.II/matrix_free/fe_evaluation.h       | 239 ++++++++++--
 tests/matrix_free/copy_feevaluation.cc        | 349 ++++++++++++++++++
 tests/matrix_free/copy_feevaluation.output    |  12 +
 4 files changed, 592 insertions(+), 70 deletions(-)
 create mode 100644 tests/matrix_free/copy_feevaluation.cc
 create mode 100644 tests/matrix_free/copy_feevaluation.output

diff --git a/deal.II/include/deal.II/base/aligned_vector.h b/deal.II/include/deal.II/base/aligned_vector.h
index c013786d69..2916bcd0ba 100644
--- a/deal.II/include/deal.II/base/aligned_vector.h
+++ b/deal.II/include/deal.II/base/aligned_vector.h
@@ -23,7 +23,6 @@
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/memory_consumption.h>
 #include <deal.II/base/parallel.h>
-#include <deal.II/base/vectorization.h>
 
 
 #if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0
@@ -76,7 +75,8 @@ public:
    * Sets the vector size to the given size and initializes all elements with
    * T().
    */
-  AlignedVector (const size_type size);
+  AlignedVector (const size_type size,
+                 const T        &init = T());
 
   /**
    * Destructor.
@@ -259,13 +259,10 @@ namespace internal
     static const std::size_t minimum_parallel_grain_size = 160000/sizeof(T)+1;
   public:
     /**
-     * Constructor. Issues a parallel call if
-     * there are sufficiently many elements,
-     * otherwise work in serial. Copies the data
-     * from source to destination and then calls
-     * destructor on the source. If the optional
-     * argument is set to true, the source is left
-     * untouched instead.
+     * Constructor. Issues a parallel call if there are sufficiently many
+     * elements, otherwise work in serial. Copies the data from source to
+     * destination and then calls destructor on the source. If the optional
+     * argument is set to true, the source is left untouched instead.
      */
     AlignedVectorMove (T *source_begin,
                        T *source_end,
@@ -285,30 +282,26 @@ namespace internal
     }
 
     /**
-     * This method moves elements from the source
-     * to the destination given in the constructor
-     * on a subrange given by two integers.
+     * This method moves elements from the source to the destination given in
+     * the constructor on a subrange given by two integers.
      */
     virtual void apply_to_subrange (const std::size_t begin,
                                     const std::size_t end) const
     {
-      // for classes trivial assignment can use
-      // memcpy
+      // for classes trivial assignment can use memcpy
       if (std_cxx1x::is_trivial<T>::value == true)
         std::memcpy (destination_+begin, source_+begin, (end-begin)*sizeof(T));
       else if (copy_only_ == false)
         for (std::size_t i=begin; i<end; ++i)
           {
-            // initialize memory, copy, and destruct
-            new (&destination_[i]) T;
-            destination_[i] = source_[i];
+            // initialize memory (copy construct), and destruct
+            new (&destination_[i]) T(source_[i]);
             source_[i].~T();
           }
       else
         for (std::size_t i=begin; i<end; ++i)
           {
-            new (&destination_[i]) T;
-            destination_[i] = source_[i];
+            new (&destination_[i]) T(source_[i]);
           }
     }
 
@@ -329,9 +322,8 @@ namespace internal
     static const std::size_t minimum_parallel_grain_size = 160000/sizeof(T)+1;
   public:
     /**
-     * Constructor. Issues a parallel call if
-     * there are sufficiently many elements,
-     * otherwise work in serial.
+     * Constructor. Issues a parallel call if there are sufficiently many
+     * elements, otherwise work in serial.
      */
     AlignedVectorSet (const std::size_t size,
                       const T &element,
@@ -359,23 +351,18 @@ namespace internal
   private:
 
     /**
-     * This sets elements on a subrange given by
-     * two integers.
+     * This sets elements on a subrange given by two integers.
      */
     virtual void apply_to_subrange (const std::size_t begin,
                                     const std::size_t end) const
     {
-      // for classes with trivial assignment of zero
-      // can use memset
+      // for classes with trivial assignment of zero can use memset
       if (std_cxx1x::is_trivial<T>::value == true && trivial_element)
         std::memset (destination_+begin, 0, (end-begin)*sizeof(T));
       else
+        // initialize memory and set
         for (std::size_t i=begin; i<end; ++i)
-          {
-            // initialize memory and set
-            new (&destination_[i]) T;
-            destination_[i] = element_;
-          }
+          new (&destination_[i]) T(element_);
     }
 
     const T &element_;
@@ -401,14 +388,15 @@ AlignedVector<T>::AlignedVector ()
 
 template < class T >
 inline
-AlignedVector<T>::AlignedVector (const size_type size)
+AlignedVector<T>::AlignedVector (const size_type size,
+                                 const T &init)
   :
   _data (0),
   _end_data (0),
   _end_allocated (0)
 {
   if (size > 0)
-    resize (size);
+    resize (size, init);
 }
 
 
@@ -504,11 +492,9 @@ AlignedVector<T>::reserve (const size_type size_alloc)
 
 #if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0
 
-      // allocate and align along boundaries of the size of
-      // VectorizedArray<double>, which is 16 bytes for SSE and 32 bytes for
-      // AVX
-      T *new_data = static_cast<T *>(_mm_malloc (size_actual_allocate,
-                                                 sizeof(VectorizedArray<double>)));
+      // allocate and align along 64-byte boundaries (this is enough for all
+      // levels of vectorization currently supported by deal.II)
+      T *new_data = static_cast<T *>(_mm_malloc (size_actual_allocate, 64));
 #else
       T *new_data = static_cast<T *>(malloc (size_actual_allocate));
 #endif
diff --git a/deal.II/include/deal.II/matrix_free/fe_evaluation.h b/deal.II/include/deal.II/matrix_free/fe_evaluation.h
index 92af89a2f7..141ebe579e 100644
--- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h
+++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h
@@ -88,11 +88,6 @@ public:
    * @name 1: General operations
    */
   //@{
-  /**
-   * Destructor.
-   */
-  ~FEEvaluationBase();
-
   /**
    * Initializes the operation pointer to the current cell. Unlike the
    * FEValues::reinit function, where the information related to a particular
@@ -535,6 +530,11 @@ protected:
                     const DoFHandler<dim>               &dof_handler,
                     const unsigned int                   base_element = 0);
 
+  /**
+   * Copy constructor
+   */
+  FEEvaluationBase (const FEEvaluationBase &other);
+
   /**
    * A unified function to read from and write into vectors based on the given
    * template operation. It can perform the operation for @p read_dof_values,
@@ -641,10 +641,10 @@ protected:
    * Stores a reference to the unit cell data, i.e., values, gradients and
    * Hessians in 1D at the quadrature points that constitute the tensor
    * product. Also contained in matrix_info, but it simplifies code if we
-   * store a reference to it. If initialized without MatrixFree object, this
-   * call actually initializes the evaluation.
+   * store a reference to it. If the object is initialized without MatrixFree
+   * object, the constructor creates this data structure.
    */
-  const internal::MatrixFreeFunctions::ShapeInfo<Number> *data;
+  std_cxx1x::shared_ptr<const internal::MatrixFreeFunctions::ShapeInfo<Number> > data;
 
   /**
    * A pointer to the Cartesian Jacobian information of the present cell. Only
@@ -797,7 +797,6 @@ protected:
                       const unsigned int            dofs_per_cell,
                       const unsigned int            n_q_points);
 
-
   /**
    * Constructor that comes with reduced functionality and works similar as
    * FEValues. The user has to provide a structure of type EvaluatedGeometry
@@ -817,6 +816,11 @@ protected:
   FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
                       const DoFHandler<dim>               &dof_handler,
                       const unsigned int                   base_element = 0);
+
+  /**
+   * Copy constructor
+   */
+  FEEvaluationAccess (const FEEvaluationAccess &other);
 };
 
 
@@ -959,6 +963,11 @@ protected:
   FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
                       const DoFHandler<dim>               &dof_handler,
                       const unsigned int                   base_element = 0);
+
+  /**
+   * Copy constructor
+   */
+  FEEvaluationAccess (const FEEvaluationAccess &other);
 };
 
 
@@ -1110,6 +1119,11 @@ protected:
   FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
                       const DoFHandler<dim>               &dof_handler,
                       const unsigned int                   base_element = 0);
+
+  /**
+   * Copy constructor
+   */
+  FEEvaluationAccess (const FEEvaluationAccess &other);
 };
 
 
@@ -1197,6 +1211,11 @@ public:
                        const DoFHandler<dim>               &dof_handler,
                        const unsigned int                   base_element = 0);
 
+  /**
+   * Copy constructor
+   */
+  FEEvaluationGeneral (const FEEvaluationGeneral &other);
+
   /**
    * Evaluates the function values, the gradients, and the Laplacians of the
    * FE function given at the DoF values in the input vector at the quadrature
@@ -1267,6 +1286,12 @@ protected:
    * Internally stored variables for the different data fields.
    */
   VectorizedArray<Number> my_data_array[n_components*(dofs_per_cell+(dim*dim+2*dim+1)*n_q_points)];
+
+private:
+  /**
+   * Sets the pointers from the data array to values_dof, etc.
+   */
+  void set_data_pointers();
 };
 
 
@@ -1347,6 +1372,11 @@ public:
                 const unsigned int            fe_no   = 0,
                 const unsigned int            quad_no = 0);
 
+  /**
+   * Copy constructor
+   */
+  FEEvaluation (const FEEvaluation &other);
+
   /**
    * Constructor that comes with reduced functionality and works similar as
    * FEValues. The user has to provide a structure of type EvaluatedGeometry
@@ -1521,6 +1551,11 @@ public:
                   const DoFHandler<dim>               &dof_handler,
                   const unsigned int                   base_element = 0);
 
+  /**
+   * Copy constructor
+   */
+  FEEvaluationGL (const FEEvaluationGL &other);
+
   /**
    * Evaluates the function values, the gradients, and the Hessians of the FE
    * function given at the DoF values in the input vector at the quadrature
@@ -1643,6 +1678,11 @@ public:
                    const DoFHandler<dim>               &dof_handler,
                    const unsigned int                   base_element = 0);
 
+  /**
+   * Copy constructor
+   */
+  FEEvaluationDGP (const FEEvaluationDGP &other);
+
   /**
    * Evaluates the function values, the gradients, and the Hessians of the FE
    * function given at the DoF values in the input vector at the quadrature
@@ -1672,6 +1712,34 @@ public:
 #ifndef DOXYGEN
 
 
+namespace internal
+{
+  namespace MatrixFreeFunctions
+  {
+    // a small class that gives control over the delete behavior of
+    // std::shared_ptr: we need to disable it when we initialize a pointer
+    // from another structure.
+    template <typename CLASS>
+    struct DummyDeleter
+    {
+      DummyDeleter (const bool do_delete = false)
+        :
+        do_delete(do_delete)
+      {}
+
+      void operator () (CLASS *pointer)
+      {
+        if (do_delete)
+          delete pointer;
+      }
+
+      const bool do_delete;
+    };
+  }
+}
+
+
+
 /*----------------------- FEEvaluationBase ----------------------------------*/
 
 template <int dim, int n_components_, typename Number>
@@ -1695,7 +1763,9 @@ FEEvaluationBase<dim,n_components_,Number>
   mapping_info       (&data_in.get_mapping_info()),
   data               (&data_in.get_shape_info
                       (fe_no_in, quad_no_in, active_fe_index,
-                       active_quad_index)),
+                       active_quad_index),
+                      internal::MatrixFreeFunctions::DummyDeleter
+                      <const internal::MatrixFreeFunctions::ShapeInfo<Number> >(false)),
   cartesian_data     (0),
   jacobian           (0),
   J_value            (0),
@@ -1786,16 +1856,44 @@ FEEvaluationBase<dim,n_components_,Number>
 
 template <int dim, int n_components_, typename Number>
 inline
-FEEvaluationBase<dim,n_components_,Number>::~FEEvaluationBase()
+FEEvaluationBase<dim,n_components_,Number>
+::FEEvaluationBase (const FEEvaluationBase<dim,n_components_,Number> &other)
+  :
+  quad_no            (other.quad_no),
+  n_fe_components    (other.n_fe_components),
+  active_fe_index    (other.active_fe_index),
+  active_quad_index  (other.active_quad_index),
+  matrix_info        (other.matrix_info),
+  dof_info           (other.dof_info),
+  mapping_info       (other.mapping_info),
+  data               (other.data),
+  cartesian_data     (other.cartesian_data),
+  jacobian           (other.jacobian),
+  J_value            (other.J_value),
+  quadrature_weights (other.quadrature_weights),
+  quadrature_points  (other.quadrature_points),
+  jacobian_grad      (other.jacobian_grad),
+  jacobian_grad_upper(other.jacobian_grad_upper),
+  cell               (other.cell),
+  cell_type          (other.cell_type),
+  cell_data_number   (other.cell_data_number),
+  evaluated_geometry (other.evaluated_geometry),
+  dof_handler        (other.dof_handler)
 {
-  // delete memory held by data in case this structure was initialized without
-  // a MatrixFree object which could hold the data structure.
-  if (matrix_info == 0)
-    delete data;
+  for (unsigned int c=0; c<n_components_; ++c)
+    {
+      values_dofs[c] = 0;
+      values_quad[c] = 0;
+      for (unsigned int d=0; d<dim; ++d)
+        gradients_quad[c][d] = 0;
+      for (unsigned int d=0; d<(dim*dim+dim)/2; ++d)
+        hessians_quad[c][d] = 0;
+    }
 }
 
 
 
+
 template <int dim, int n_components_, typename Number>
 inline
 void
@@ -3460,6 +3558,16 @@ FEEvaluationAccess<dim,n_components_,Number>
 
 
 
+template <int dim, int n_components_, typename Number>
+inline
+FEEvaluationAccess<dim,n_components_,Number>
+::FEEvaluationAccess (const FEEvaluationAccess<dim,n_components_,Number> &other)
+  :
+  FEEvaluationBase <dim,n_components_,Number>(other)
+{}
+
+
+
 
 /*-------------------- FEEvaluationAccess scalar ----------------------------*/
 
@@ -3491,6 +3599,16 @@ FEEvaluationAccess<dim,1,Number>
 
 
 
+template <int dim, typename Number>
+inline
+FEEvaluationAccess<dim,1,Number>
+::FEEvaluationAccess (const FEEvaluationAccess<dim,1,Number>&other)
+  :
+  FEEvaluationBase <dim,1,Number>(other)
+{}
+
+
+
 template <int dim, typename Number>
 inline
 VectorizedArray<Number>
@@ -3715,6 +3833,16 @@ FEEvaluationAccess<dim,dim,Number>
 
 
 
+template <int dim, typename Number>
+inline
+FEEvaluationAccess<dim,dim,Number>
+::FEEvaluationAccess (const FEEvaluationAccess<dim,dim,Number>&other)
+  :
+  FEEvaluationBase <dim,dim,Number>(other)
+{}
+
+
+
 template <int dim, typename Number>
 inline
 Tensor<2,dim,VectorizedArray<Number> >
@@ -4043,20 +4171,8 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
   :
   BaseClass (data_in, fe_no, quad_no, dofs_per_cell, n_q_points)
 {
-  // set the pointers to the correct position in the data array
-  for (unsigned int c=0; c<n_components_; ++c)
-    {
-      this->values_dofs[c] = &my_data_array[c*dofs_per_cell];
-      this->values_quad[c] = &my_data_array[n_components*dofs_per_cell+c*n_q_points];
-      for (unsigned int d=0; d<dim; ++d)
-        this->gradients_quad[c][d] = &my_data_array[n_components*(dofs_per_cell+n_q_points)
-                                                    +
-                                                    (c*dim+d)*n_q_points];
-      for (unsigned int d=0; d<(dim*dim+dim)/2; ++d)
-        this->hessians_quad[c][d] = &my_data_array[n_components*((dim+1)*n_q_points+dofs_per_cell)
-                                                   +
-                                                   (c*(dim*dim+dim)+d)*n_q_points];
-    }
+  set_data_pointers();
+
 #ifdef DEBUG
   // print error message when the dimensions do not match. Propose a possible
   // fix
@@ -4173,6 +4289,30 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
                        const unsigned int                   base_element)
   :
   BaseClass (geometry, dof_handler, base_element)
+{
+  set_data_pointers();
+}
+
+
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
+          typename Number>
+inline
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluationGeneral (const FEEvaluationGeneral &other)
+  :
+  BaseClass (other)
+{
+  set_data_pointers();
+}
+
+
+
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
+          typename Number>
+inline
+void
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::set_data_pointers()
 {
   // set the pointers to the correct position in the data array
   for (unsigned int c=0; c<n_components_; ++c)
@@ -5788,6 +5928,19 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 
 
 
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
+          typename Number>
+inline
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluation (const FEEvaluation &other)
+  :
+  BaseClass (other)
+{
+  compute_even_odd_factors();
+}
+
+
+
 template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
           typename Number>
 inline
@@ -5957,7 +6110,7 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
   if (fe_degree > 1 || n_q_points_1d > 3)
     internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
              VectorizedArray<Number>, direction, dof_to_quad, add, 0>
-             (shape_val_evenodd, in, out);
+             (&shape_val_evenodd[0], in, out);
   else
     internal::apply_tensor_product_values<dim,fe_degree,n_q_points_1d,
              VectorizedArray<Number>, direction, dof_to_quad, add>
@@ -5978,7 +6131,7 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
   if (fe_degree > 1 || n_q_points_1d > 3)
     internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
              VectorizedArray<Number>, direction, dof_to_quad, add, 1>
-             (shape_gra_evenodd, in, out);
+             (&shape_gra_evenodd[0], in, out);
   else
     internal::apply_tensor_product_gradients<dim,fe_degree,n_q_points_1d,
              VectorizedArray<Number>, direction, dof_to_quad, add>
@@ -6002,7 +6155,7 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
   if (fe_degree > 1 || n_q_points_1d > 3)
     internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
              VectorizedArray<Number>, direction, dof_to_quad, add, 2>
-             (shape_hes_evenodd, in, out);
+             (&shape_hes_evenodd[0], in, out);
   else
     internal::apply_tensor_product_hessians<dim,fe_degree,n_q_points_1d,
              VectorizedArray<Number>, direction, dof_to_quad, add>
@@ -6069,6 +6222,16 @@ FEEvaluationGL<dim,fe_degree,n_components_,Number>
 
 
 
+template <int dim, int fe_degree, int n_components_, typename Number>
+inline
+FEEvaluationGL<dim,fe_degree,n_components_,Number>
+::FEEvaluationGL (const FEEvaluationGL &other)
+  :
+  BaseClass (other)
+{}
+
+
+
 template <int dim, int fe_degree, int n_components_, typename Number>
 inline
 void
@@ -6293,6 +6456,18 @@ FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
 {}
 
 
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+          typename Number>
+inline
+FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluationDGP (const FEEvaluationDGP &other)
+  :
+  BaseClass (other)
+{}
+
+
+
 template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
           typename Number>
 inline
diff --git a/tests/matrix_free/copy_feevaluation.cc b/tests/matrix_free/copy_feevaluation.cc
new file mode 100644
index 0000000000..cb0d126b22
--- /dev/null
+++ b/tests/matrix_free/copy_feevaluation.cc
@@ -0,0 +1,349 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// same test as matrix_vector_stokes_noflux, but allocating FEEvaluation on
+// the heap (using AlignedVector) instead of allocating it on the stack. Tests
+// also copy constructors of FEEvaluation.
+
+#include "../tests.h"
+
+std::ofstream logfile("output");
+
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/aligned_vector.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_sparsity_pattern.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+
+
+template <int dim, int degree_p, typename VectorType>
+class MatrixFreeTest
+{
+public:
+  typedef typename DoFHandler<dim>::active_cell_iterator CellIterator;
+  typedef double Number;
+
+  MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+    data (data_in)
+  {};
+
+  void
+  local_apply (const MatrixFree<dim,Number> &data,
+               VectorType          &dst,
+               const VectorType    &src,
+               const std::pair<unsigned int,unsigned int> &cell_range) const
+  {
+    typedef VectorizedArray<Number> vector_t;
+    // allocate FEEvaluation. This test will test proper alignment
+    AlignedVector<FEEvaluation<dim,degree_p+1,degree_p+2,dim,Number> > velocity
+      (1, FEEvaluation<dim,degree_p+1,degree_p+2,dim,Number>(data, 0));
+    AlignedVector<FEEvaluation<dim,degree_p,  degree_p+2,1,  Number> > pressure
+      (1, FEEvaluation<dim,degree_p,  degree_p+2,1,  Number>(data, 1));
+
+    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+      {
+        velocity[0].reinit (cell);
+        velocity[0].read_dof_values (src.block(0));
+        velocity[0].evaluate (false,true,false);
+        pressure[0].reinit (cell);
+        pressure[0].read_dof_values (src.block(1));
+        pressure[0].evaluate (true,false,false);
+
+        for (unsigned int q=0; q<velocity[0].n_q_points; ++q)
+          {
+            SymmetricTensor<2,dim,vector_t> sym_grad_u =
+              velocity[0].get_symmetric_gradient (q);
+            vector_t pres = pressure[0].get_value(q);
+            vector_t div = -velocity[0].get_divergence(q);
+            pressure[0].submit_value   (div, q);
+
+            // subtract p * I
+            for (unsigned int d=0; d<dim; ++d)
+              sym_grad_u[d][d] -= pres;
+
+            velocity[0].submit_symmetric_gradient(sym_grad_u, q);
+          }
+
+        velocity[0].integrate (false,true);
+        velocity[0].distribute_local_to_global (dst.block(0));
+        pressure[0].integrate (true,false);
+        pressure[0].distribute_local_to_global (dst.block(1));
+      }
+  }
+
+
+  void vmult (VectorType &dst,
+              const VectorType &src) const
+  {
+    dst = 0;
+    data.cell_loop (&MatrixFreeTest<dim,degree_p,VectorType>::local_apply,
+                    this, dst, src);
+  };
+
+private:
+  const MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+  Triangulation<dim>   triangulation;
+  GridGenerator::hyper_shell (triangulation, Point<dim>(),
+                              0.5, 1., 96, true);
+  static HyperShellBoundary<dim> boundary;
+  triangulation.set_boundary (0, boundary);
+  triangulation.set_boundary (1, boundary);
+  triangulation.begin_active()->set_refine_flag();
+  triangulation.last()->set_refine_flag();
+  triangulation.execute_coarsening_and_refinement();
+  triangulation.refine_global (3-dim);
+  triangulation.last()->set_refine_flag();
+  triangulation.execute_coarsening_and_refinement();
+
+  MappingQ<dim>        mapping (3);
+  FE_Q<dim>            fe_u_scal (fe_degree+1);
+  FESystem<dim>        fe_u (fe_u_scal,dim);
+  FE_Q<dim>            fe_p (fe_degree);
+  FESystem<dim>        fe (fe_u_scal, dim, fe_p, 1);
+  DoFHandler<dim>      dof_handler_u (triangulation);
+  DoFHandler<dim>      dof_handler_p (triangulation);
+  DoFHandler<dim>      dof_handler (triangulation);
+
+  MatrixFree<dim,double> mf_data;
+
+  ConstraintMatrix     constraints, constraints_u, constraints_p;
+
+  BlockSparsityPattern      sparsity_pattern;
+  BlockSparseMatrix<double> system_matrix;
+
+  BlockVector<double> solution;
+  BlockVector<double> system_rhs;
+  BlockVector<double> mf_solution;
+
+  dof_handler.distribute_dofs (fe);
+  dof_handler_u.distribute_dofs (fe_u);
+  dof_handler_p.distribute_dofs (fe_p);
+  std::vector<unsigned int> stokes_sub_blocks (dim+1,0);
+  stokes_sub_blocks[dim] = 1;
+  DoFRenumbering::component_wise (dof_handler, stokes_sub_blocks);
+
+  std::set<unsigned char> no_normal_flux_boundaries;
+  no_normal_flux_boundaries.insert (0);
+  no_normal_flux_boundaries.insert (1);
+  DoFTools::make_hanging_node_constraints (dof_handler,
+                                           constraints);
+  VectorTools::compute_no_normal_flux_constraints (dof_handler, 0,
+                                                   no_normal_flux_boundaries,
+                                                   constraints, mapping);
+  constraints.close ();
+  DoFTools::make_hanging_node_constraints (dof_handler_u,
+                                           constraints_u);
+  VectorTools::compute_no_normal_flux_constraints (dof_handler_u, 0,
+                                                   no_normal_flux_boundaries,
+                                                   constraints_u, mapping);
+  constraints_u.close ();
+  DoFTools::make_hanging_node_constraints (dof_handler_p,
+                                           constraints_p);
+  constraints_p.close ();
+
+  std::vector<types::global_dof_index> dofs_per_block (2);
+  DoFTools::count_dofs_per_block (dof_handler, dofs_per_block,
+                                  stokes_sub_blocks);
+
+  //std::cout << "Number of active cells: "
+  //          << triangulation.n_active_cells()
+  //          << std::endl
+  //          << "Number of degrees of freedom: "
+  //          << dof_handler.n_dofs()
+  //          << " (" << n_u << '+' << n_p << ')'
+  //          << std::endl;
+
+  {
+    BlockCompressedSimpleSparsityPattern csp (2,2);
+
+    for (unsigned int d=0; d<2; ++d)
+      for (unsigned int e=0; e<2; ++e)
+        csp.block(d,e).reinit (dofs_per_block[d], dofs_per_block[e]);
+
+    csp.collect_sizes();
+
+    DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false);
+    sparsity_pattern.copy_from (csp);
+  }
+
+  system_matrix.reinit (sparsity_pattern);
+
+  // this is from step-22
+  {
+    QGauss<dim>   quadrature_formula(fe_degree+2);
+
+    FEValues<dim> fe_values (mapping, fe, quadrature_formula,
+                             update_values    |
+                             update_JxW_values |
+                             update_gradients);
+
+    const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+    const unsigned int   n_q_points      = quadrature_formula.size();
+
+    FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+    const FEValuesExtractors::Vector velocities (0);
+    const FEValuesExtractors::Scalar pressure (dim);
+
+    std::vector<SymmetricTensor<2,dim> > phi_grads_u (dofs_per_cell);
+    std::vector<double>                  div_phi_u   (dofs_per_cell);
+    std::vector<double>                  phi_p       (dofs_per_cell);
+
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+        fe_values.reinit (cell);
+        local_matrix = 0;
+
+        for (unsigned int q=0; q<n_q_points; ++q)
+          {
+            for (unsigned int k=0; k<dofs_per_cell; ++k)
+              {
+                phi_grads_u[k] = fe_values[velocities].symmetric_gradient (k, q);
+                div_phi_u[k]   = fe_values[velocities].divergence (k, q);
+                phi_p[k]       = fe_values[pressure].value (k, q);
+              }
+
+            for (unsigned int i=0; i<dofs_per_cell; ++i)
+              {
+                for (unsigned int j=0; j<=i; ++j)
+                  {
+                    local_matrix(i,j) += (phi_grads_u[i] * phi_grads_u[j]
+                                          - div_phi_u[i] * phi_p[j]
+                                          - phi_p[i] * div_phi_u[j])
+                                         * fe_values.JxW(q);
+                  }
+              }
+          }
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          for (unsigned int j=i+1; j<dofs_per_cell; ++j)
+            local_matrix(i,j) = local_matrix(j,i);
+
+        cell->get_dof_indices (local_dof_indices);
+        constraints.distribute_local_to_global (local_matrix,
+                                                local_dof_indices,
+                                                system_matrix);
+      }
+  }
+
+
+  solution.reinit (2);
+  for (unsigned int d=0; d<2; ++d)
+    solution.block(d).reinit (dofs_per_block[d]);
+  solution.collect_sizes ();
+
+  system_rhs.reinit (solution);
+  mf_solution.reinit (solution);
+
+  // fill system_rhs with random numbers
+  for (unsigned int j=0; j<system_rhs.block(0).size(); ++j)
+    if (constraints_u.is_constrained(j) == false)
+      {
+        const double val = -1 + 2.*(double)Testing::rand()/double(RAND_MAX);
+        system_rhs.block(0)(j) = val;
+      }
+  for (unsigned int j=0; j<system_rhs.block(1).size(); ++j)
+    if (constraints_p.is_constrained(j) == false)
+      {
+        const double val = -1 + 2.*(double)Testing::rand()/double(RAND_MAX);
+        system_rhs.block(1)(j) = val;
+      }
+
+  // setup matrix-free structure
+  {
+    std::vector<const DoFHandler<dim>*> dofs;
+    dofs.push_back(&dof_handler_u);
+    dofs.push_back(&dof_handler_p);
+    std::vector<const ConstraintMatrix *> constraints;
+    constraints.push_back (&constraints_u);
+    constraints.push_back (&constraints_p);
+    QGauss<1> quad(fe_degree+2);
+    // no parallelism
+    mf_data.reinit (mapping, dofs, constraints, quad,
+                    typename MatrixFree<dim>::AdditionalData
+                    (MPI_COMM_WORLD,
+                     MatrixFree<dim>::AdditionalData::none));
+  }
+
+  system_matrix.vmult (solution, system_rhs);
+
+  MatrixFreeTest<dim,fe_degree,BlockVector<double> > mf (mf_data);
+  mf.vmult (mf_solution, system_rhs);
+
+  // Verification
+  mf_solution -= solution;
+  const double error = mf_solution.linfty_norm();
+  const double relative = solution.linfty_norm();
+  deallog << "Verification fe degree " << fe_degree  <<  ": "
+          << error/relative << std::endl << std::endl;
+}
+
+
+
+int main ()
+{
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  deallog << std::setprecision (3);
+
+  {
+    deallog << std::endl << "Test with doubles" << std::endl << std::endl;
+    deallog.threshold_double(1.e-12);
+    deallog.push("2d");
+    test<2,1>();
+    test<2,2>();
+    test<2,3>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3,1>();
+    deallog.pop();
+  }
+}
+
diff --git a/tests/matrix_free/copy_feevaluation.output b/tests/matrix_free/copy_feevaluation.output
new file mode 100644
index 0000000000..0434059632
--- /dev/null
+++ b/tests/matrix_free/copy_feevaluation.output
@@ -0,0 +1,12 @@
+
+DEAL::
+DEAL::Test with doubles
+DEAL::
+DEAL:2d::Verification fe degree 1: 0
+DEAL:2d::
+DEAL:2d::Verification fe degree 2: 0
+DEAL:2d::
+DEAL:2d::Verification fe degree 3: 0
+DEAL:2d::
+DEAL:3d::Verification fe degree 1: 0
+DEAL:3d::
-- 
2.39.5