]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fig bug and augment FEEValuation. Rename MatrixFree::get_quadrature_formula into...
authorkormann <kormann@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 13 Jun 2012 19:52:16 +0000 (19:52 +0000)
committerkormann <kormann@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 13 Jun 2012 19:52:16 +0000 (19:52 +0000)
git-svn-id: https://svn.dealii.org/trunk@25626 0785d39b-7218-0410-832d-ea1e28bc413d

37 files changed:
deal.II/examples/step-37/doc/intro.dox
deal.II/examples/step-37/doc/results.dox
deal.II/examples/step-37/step-37.cc
deal.II/examples/step-48/doc/intro.dox
deal.II/examples/step-48/doc/results.dox
deal.II/examples/step-48/step-48.cc
deal.II/include/deal.II/matrix_free/dof_info.templates.h
deal.II/include/deal.II/matrix_free/fe_evaluation.h
deal.II/include/deal.II/matrix_free/mapping_info.h
deal.II/include/deal.II/matrix_free/mapping_info.templates.h
deal.II/include/deal.II/matrix_free/matrix_free.h
deal.II/include/deal.II/matrix_free/matrix_free.templates.h
deal.II/source/numerics/matrix_free.inst.in
tests/matrix_free/get_functions_common.h
tests/matrix_free/get_functions_gl.cc
tests/matrix_free/get_functions_multife.cc
tests/matrix_free/get_functions_multife2.cc
tests/matrix_free/get_functions_q_hierarchical.cc
tests/matrix_free/get_functions_rect.cc
tests/matrix_free/get_functions_variants.cc
tests/matrix_free/get_values_plain.cc
tests/matrix_free/integrate_functions.cc
tests/matrix_free/integrate_functions_multife.cc
tests/matrix_free/integrate_functions_multife2.cc
tests/matrix_free/matrix_vector_06.cc
tests/matrix_free/matrix_vector_06/cmp/generic
tests/matrix_free/matrix_vector_08.cc
tests/matrix_free/matrix_vector_08/cmp/generic
tests/matrix_free/matrix_vector_09.cc
tests/matrix_free/matrix_vector_09/cmp/generic
tests/matrix_free/matrix_vector_hp.cc
tests/matrix_free/matrix_vector_stokes.cc
tests/matrix_free/matrix_vector_stokes_noflux.cc
tests/matrix_free/matrix_vector_stokes_noflux/cmp/generic
tests/matrix_free/quadrature_points.cc
tests/matrix_free/thread_correctness_hp.cc
tests/matrix_free/thread_correctness_hp/cmp/generic

index fd54c93359d3cd5578586e488ded15d0be7ae0cb..19e8a07e9f28a8ac27bc58e8005b9d39631d8c60 100644 (file)
@@ -8,7 +8,7 @@ The algorithm for the matrix-vector product is based on the article <a
 href="http://dx.doi.org/10.1016/j.compfluid.2012.04.012">A generic
 interface for parallel cell-based finite element operator
 application</a> by Martin Kronbichler and Katharina Kormann, Computers
-and Fluids, 2012, and the paper &quot;Parallel finite element operator
+and Fluids 63:135&ndash;147, 2012, and the paper &quot;Parallel finite element operator
 application: Graph partitioning and coloring&quot; by Katharina
 Kormann and Martin Kronbichler in: Proceedings of the 7th IEEE
 International Conference on e-Science, 2011.  </i>
index 8c269d575b46eeb75a14f4f6921fcb7778757fd6..ab9a9b0648f2d923476115474d58fad6fe961c00 100644 (file)
@@ -204,8 +204,9 @@ factor 12. This is the threshold when the cache in the processor can no longer
 hold all data necessary for the matrix-vector products and all matrix elements
 must be fetched from main memory. The second deviation is the times for the
 matrix-free solve which increase by less than a factor 8. This is because of
-more parallelism from more cells, exploited by the (involved) dynamic tasks
-scheduling approach CellFEOperator. Note that about 30% of the time in the
+more parallelism from more cells, exploited by the (involved) dynamic task
+scheduling approach taken in the cell loop of the MatrixFree class. Note
+that about 30% of the time in the
 matrix-free solver is spent on restriction and prolongation, which use sparse
 matrices. So the speedup could be even better if all parts where done
 efficiently.
index bc24a553e65226ce70017c3d46dd0e5ad9ff5829..9f6e4b898771716ea0229cd20323f33effd4b582 100644 (file)
@@ -1,8 +1,9 @@
-/* Author: Katharina Kormann, Martin Kronbichler, Uppsala University, 2009-2011 */
+/* $Id$ */
+/* Author: Katharina Kormann, Martin Kronbichler, Uppsala University, 2009-2012 */
 
-/*    $Id$       */
+/*    $Id$    */
 /*                                                                */
-/*    Copyright (C) 2009, 2010, 2011, 2012 by the deal.II authors       */
+/*    Copyright (C) 2009, 2010, 2011, 2012 by the deal.II authors */
 /*                                                                */
 /*    This file is subject to QPL and may not be  distributed     */
 /*    without copyright and license information. Please refer     */
@@ -189,7 +190,7 @@ namespace Step37
                                    // writing to vectors), and hence SIMD can
                                    // be used to process several cells at
                                    // once. In all what follows, you can think
-                                   // of an AlignedVector to hold data from
+                                   // of a VectorizedArray to hold data from
                                    // several cells. For example, we evaluate
                                    // the coefficient shown here not on a
                                    // simple point as usually done, but we
@@ -206,7 +207,7 @@ namespace Step37
                                    // contains the x-coordinate for the first
                                    // point and the second point. You may
                                    // access the coordinates individually
-                                   // using e.g. <code>p[0][j]</code>, j=1,2,
+                                   // using e.g. <code>p[0][j]</code>, j=0,1,
                                    // but it is recommended to define
                                    // operations on a VectorizedArray as much
                                    // as possible in order to make use of
index b1ceccc5535d557e609cb0d9c96c894adc14fc71..04d9ee7e4ec34479e28550dad712cc5c0db4b19f 100644 (file)
@@ -7,7 +7,7 @@ The algorithm for the matrix-vector product is based on the article <a
 href="http://dx.doi.org/10.1016/j.compfluid.2012.04.012">A generic
 interface for parallel cell-based finite element operator
 application</a> by Martin Kronbichler and Katharina Kormann, Computers
-and Fluids, 2012, and the paper &quot;Parallel finite element operator
+and Fluids 63:135&ndash;147, 2012, and the paper &quot;Parallel finite element operator
 application: Graph partitioning and coloring&quot; by Katharina
 Kormann and Martin Kronbichler in: Proceedings of the 7th IEEE
 International Conference on e-Science, 2011.  </i>
@@ -23,10 +23,12 @@ MatrixFree class handles constraints and how it can be
 parallelized over distributed nodes. Finally, we will use an explicit
 time-stepping method to solve the problem and introduce Gauss-Lobatto
 finite elements that are very convenient in this case since they have
-a diagonally, and thus trivially invertable, mass matrix. Moreover,
+a diagonal, and thus trivially invertible, mass matrix. Moreover,
 this type of elements clusters the nodes towards the element
-boundaries which is why it has good properties for high-order
-discretization methods.
+boundaries which is why they have good properties for high-order
+discretization methods. Indeed, the condition number of standard FE_Q
+elements with equidistant nodes grows exponentially with the degree,
+which destroys any benefit for orders of about five and higher.
 
 <h3> Problem statement and discretization </h3>
 
index 34e7cf6dbdb56581d789a106a2eb8b9c77ed3365..52308a17545e15d662d3a7d2b5dfbab3ee73a25c 100644 (file)
@@ -2,7 +2,7 @@
 
 <h3>Comparison with a sparse matrix</h3>
 
-In order to demonstrate the gain in using the CellFEOperator class instead of
+In order to demonstrate the gain in using the MatrixFree class instead of
 the standard <code>deal.II</code> assembly routines for evaluating the
 information from old time steps, we study a simple serial run of the code on a
 nonadaptive mesh. Since much time is spent on evaluating the sine function, we
@@ -59,7 +59,7 @@ following table.
    </tr>
 </table>
 
-It is apparent that the CellFEOperator outperforms the standard assembly
+It is apparent that the matrix-free code outperforms the standard assembly
 routines in deal.II by far. In 3D and for fourth order elements, one operator
 application is also almost ten times as fast as a sparse matrix-vector
 product.
index ddd3174156a8c2c41707e86932589f864f4edd7f..451281dc28532359610579b6c8d56af061d098e6 100644 (file)
@@ -1,9 +1,9 @@
 /* $Id$ */
-/* Author: Katharina Kormann, Martin Kronbichler, Uppsala University, 2011 */
+/* Author: Katharina Kormann, Martin Kronbichler, Uppsala University, 2011-2012 */
 
 /*    $Id$       */
 /*                                                                */
-/*    Copyright (C) 2011, 2012 by the deal.II authors       */
+/*    Copyright (C) 2011, 2012 by the deal.II authors             */
 /*                                                                */
 /*    This file is subject to QPL and may not be  distributed     */
 /*    without copyright and license information. Please refer     */
@@ -96,8 +96,6 @@ namespace Step48
   class SineGordonOperation
   {
     public:
-      typedef VectorizedArray<double> vector_t;
-
       SineGordonOperation(const MatrixFree<dim,double> &data_in,
                           const double                      time_step);
 
@@ -105,12 +103,12 @@ namespace Step48
                   const std::vector<parallel::distributed::Vector<double>*> &src) const;
 
     private:
-      const MatrixFree<dim,double>     &data;
-      const vector_t                        delta_t_sqr;
+      const MatrixFree<dim,double>         &data;
+      const VectorizedArray<double>         delta_t_sqr;
       parallel::distributed::Vector<double> inv_mass_matrix;
 
-      void local_apply (const MatrixFree<dim,double>                  &data,
-                        parallel::distributed::Vector<double>             &dst,
+      void local_apply (const MatrixFree<dim,double>               &data,
+                        parallel::distributed::Vector<double>      &dst,
                         const std::vector<parallel::distributed::Vector<double>*>&src,
                         const std::pair<unsigned int,unsigned int> &cell_range) const;
   };
@@ -141,12 +139,12 @@ namespace Step48
   template <int dim, int fe_degree>
   SineGordonOperation<dim,fe_degree>::
   SineGordonOperation(const MatrixFree<dim,double> &data_in,
-                      const double                      time_step)
+                      const double                  time_step)
                   :
                   data(data_in),
                   delta_t_sqr(make_vectorized_array(time_step*time_step))
   {
-    vector_t one = make_vectorized_array (1.);
+    VectorizedArray<double> one = make_vectorized_array (1.);
 
     data.initialize_dof_vector (inv_mass_matrix);
 
@@ -242,8 +240,8 @@ namespace Step48
 
         for (unsigned int q=0; q<current.n_q_points; ++q)
           {
-            const vector_t current_value = current.get_value(q);
-            const vector_t old_value     = old.get_value(q);
+            const VectorizedArray<double> current_value = current.get_value(q);
+            const VectorizedArray<double> old_value     = old.get_value(q);
 
             current.submit_value (2.*current_value - old_value -
                                   delta_t_sqr * std::sin(current_value),q);
@@ -269,9 +267,9 @@ namespace Step48
                                    // the cell loop is implemented in the cell
                                    // finite element operator class. On each cell
                                    // it applies the routine defined as the
-                                   // <code>operator ()</code> method of the
+                                   // <code>local_apply()</code> method of the
                                    // class <code>SineGordonOperation</code>,
-                                   // i.e., <code>*this</code>. One could also
+                                   // i.e., <code>this</code>. One could also
                                    // provide a function with the same signature
                                    // that is not part of a class.
   template <int dim, int fe_degree>
@@ -468,7 +466,7 @@ namespace Step48
                                      // to use shared-memory parallelization (hence
                                      // one would use multithreading for intra-node
                                      // parallelism and not MPI; note that we here
-                                     // choose the standard option &mdash if we
+                                     // choose the standard option &mdash; if we
                                      // wanted to disable shared memory
                                      // parallelization, we would choose @p
                                      // none). Finally, three solution vectors are
index bba63780168be2345303bed96ba043d0a9057756..8c680c96464649e855ad7cb311a990f57f745e97 100644 (file)
@@ -566,6 +566,8 @@ namespace MatrixFreeFunctions
                                        std::vector<unsigned int> &renumbering,
                                        std::vector<unsigned int> &irregular_cells)
   {
+    if (max_fe_index < 2)
+      return;
     const unsigned int n_active_cells = size_info.n_active_cells;
     const unsigned int vectorization_length = size_info.vectorization_length;
     irregular_cells.resize (0);
index 7e8a42ed09f84ac5db481b7022b4bc85a598aa67..580042ef367bb1f74635e0442ce9781fb22edcb7 100644 (file)
@@ -38,6 +38,11 @@ namespace parallel
 namespace internal
 {
   DeclException0 (ExcAccessToUninitializedField);
+
+  template <typename FEEval>
+  void do_evaluate (FEEval &, const bool, const bool, const bool);
+  template <typename FEEval>
+  void do_integrate (FEEval &, const bool, const bool);
 }
 
 
@@ -74,13 +79,15 @@ namespace internal
  * @author Katharina Kormann and Martin Kronbichler, 2010, 2011
  */
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 class FEEvaluationBase
 {
 public:
-  typedef Tensor<1,n_components,VectorizedArray<Number> > value_type;
-  typedef Tensor<1,n_components,Tensor<1,dim,VectorizedArray<Number> > > gradient_type;
+  typedef Number                            number_type;
+  typedef Tensor<1,n_components_,VectorizedArray<Number> > value_type;
+  typedef Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > > gradient_type;
   static const unsigned int dimension     = dim;
+  static const unsigned int n_components  = n_components_;
   static const unsigned int dofs_per_cell = dofs_per_cell_;
   static const unsigned int n_q_points    = n_q_points_;
 
@@ -175,7 +182,7 @@ public:
                                  */
   template <typename VectorType>
   void read_dof_values (const std::vector<VectorType> &src,
-                       const unsigned int             first_index=0);
+                        const unsigned int             first_index=0);
 
                                 /**
                                  * Reads data from several vectors. Same as
@@ -184,21 +191,7 @@ public:
                                  */
   template <typename VectorType>
   void read_dof_values (const std::vector<VectorType*> &src,
-                       const unsigned int              first_index=0);
-
-                                /**
-                                 * For a collection of several vector @p src,
-                                 * read out the values on the degrees of
-                                 * freedom of the current cell for @p
-                                 * n_components (template argument), and store
-                                 * them internally. Similar functionality as
-                                 * the function
-                                 * ConstraintMatrix::read_dof_values. Note
-                                 * that if vectorization is enabled, the DoF
-                                 * values for several cells are set.
-                                 */
-  template<typename VectorType>
-  void read_dof_values (const VectorType * src_data[]);
+                        const unsigned int              first_index=0);
 
                                 /**
                                  * For the vector @p src, read out the values
@@ -248,20 +241,6 @@ public:
   void read_dof_values_plain (const std::vector<VectorType*> &src,
                              const unsigned int              first_index=0);
 
-                                /**
-                                 * For a collection of several vector @p src,
-                                 * read out the values on the degrees of
-                                 * freedom of the current cell for @p
-                                 * n_components (template argument), and store
-                                 * them internally. Similar functionality as
-                                 * the function
-                                 * DoFAccessor::read_dof_values. Note
-                                 * that if vectorization is enabled, the DoF
-                                 * values for several cells are set.
-                                 */
-  template<typename VectorType>
-  void read_dof_values_plain (const VectorType * src_data[]);
-
                                 /**
                                  * Takes the values stored internally on dof
                                  * values of the current cell and sums them
@@ -304,23 +283,6 @@ public:
   void distribute_local_to_global (std::vector<VectorType*> &dst,
                                    const unsigned int       first_index=0) const;
 
-                                /**
-                                 * Takes the values stored internally on dof
-                                 * values of the current cell for a
-                                 * vector-valued problem consisting of @p
-                                 * n_components (template argument) and sums
-                                 * them into the collection of vectors vector
-                                 * @p dst. The function also applies
-                                 * constraints during the write operation. The
-                                 * functionality is hence similar to the
-                                 * function
-                                 * ConstraintMatrix::distribute_local_to_global.
-                                 * Note that if vectorization is enabled, the
-                                 * DoF values for several cells are used.
-                                 */
-  template<typename VectorType>
-  void distribute_local_to_global (VectorType * dst_data[]) const;
-
                                 /**
                                  * Takes the values stored internally on dof
                                  * values of the current cell and sums them
@@ -363,23 +325,6 @@ public:
   void set_dof_values (std::vector<VectorType*> &dst,
                        const unsigned int        first_index=0) const;
 
-                                /**
-                                 * Takes the values stored internally on dof
-                                 * values of the current cell for a
-                                 * vector-valued problem consisting of @p
-                                 * n_components (template argument) and sums
-                                 * them into the collection of vectors vector
-                                 * @p dst. The function also applies
-                                 * constraints during the write operation. The
-                                 * functionality is hence similar to the
-                                 * function
-                                 * ConstraintMatrix::distribute_local_to_global.
-                                 * Note that if vectorization is enabled, the
-                                 * DoF values for several cells are used.
-                                 */
-  template<typename VectorType>
-  void set_dof_values (VectorType * dst_data[]) const;
-
                                 //@}
 
                                 /**
@@ -518,7 +463,7 @@ public:
                                  * (n_components == 1) and for the
                                  * vector-valued case (n_components == dim).
                                  */
-  Tensor<1,n_components,Tensor<2,dim,VectorizedArray<Number> > >
+  Tensor<1,n_components_,Tensor<2,dim,VectorizedArray<Number> > >
   get_hessian (const unsigned int q_point) const;
 
                                 /**
@@ -569,6 +514,132 @@ public:
 
                                 //@}
 
+                                /**
+                                 * @name 4: Access to internal data
+                                 */
+                                //@{
+                                /**
+                                 * Returns a read-only pointer to the first
+                                 * field of function values on quadrature
+                                 * points. First come the function values on
+                                 * all quadrature points for the first
+                                 * component, then all values for the second
+                                 * component, and so on. This is related to
+                                 * the internal data structures used in this
+                                 * class. The raw data after a call to @p
+                                 * evaluate only contains unit cell
+                                 * operations, so possible transformations,
+                                 * quadrature weights etc. must be applied
+                                 * manually. In general, it is safer to use
+                                 * the get_value() function instead, which
+                                 * does all the transformation internally.
+                                 */
+  const VectorizedArray<Number> * begin_values () const;
+
+                                /**
+                                 * Returns a read and write pointer to the
+                                 * first field of function values on
+                                 * quadrature points. First come the function
+                                 * values on all quadrature points for the
+                                 * first component, then all values for the
+                                 * second component, and so on. This is
+                                 * related to the internal data structures
+                                 * used in this class. The raw data after a
+                                 * call to @p evaluate only contains unit
+                                 * cell operations, so possible
+                                 * transformations, quadrature weights
+                                 * etc. must be applied manually. In general,
+                                 * it is safer to use the get_value() function
+                                 * instead, which does all the transformation
+                                 * internally.
+                                 */
+  VectorizedArray<Number> * begin_values ();
+
+                                /**
+                                 * Returns a read-only pointer to the first
+                                 * field of function gradients on quadrature
+                                 * points. First comes the x-component of the
+                                 * gradient for the first component on all
+                                 * quadrature points, then the y-component,
+                                 * and so on. Next comes the x-component of
+                                 * the second component, and so on. This is
+                                 * related to the internal data structures
+                                 * used in this class. The raw data after a
+                                 * call to @p evaluate only contains unit
+                                 * cell operations, so possible
+                                 * transformations, quadrature weights
+                                 * etc. must be applied manually. In general,
+                                 * it is safer to use the get_gradient() function
+                                 * instead, which does all the transformation
+                                 * internally.
+                                 */
+  const VectorizedArray<Number> * begin_gradients () const;
+
+                                /**
+                                 * Returns a read and write pointer to the
+                                 * first field of function gradients on
+                                 * quadrature points. First comes the
+                                 * x-component of the gradient for the first
+                                 * component on all quadrature points, then
+                                 * the y-component, and so on. Next comes the
+                                 * x-component of the second component, and so
+                                 * on. This is related to the internal data
+                                 * structures used in this class. The raw data
+                                 * after a call to @p evaluate only
+                                 * contains unit cell operations, so possible
+                                 * transformations, quadrature weights
+                                 * etc. must be applied manually. In general,
+                                 * it is safer to use the get_gradient()
+                                 * function instead, which does all the
+                                 * transformation internally.
+                                 */
+  VectorizedArray<Number> * begin_gradients ();
+
+                                /**
+                                 * Returns a read-only pointer to the first
+                                 * field of function hessians on quadrature
+                                 * points. First comes the xx-component of the
+                                 * hessian for the first component on all
+                                 * quadrature points, then the yy-component,
+                                 * zz-component in (3D), then the
+                                 * xy-component, and so on. Next comes the
+                                 * xx-component of the second component, and
+                                 * so on. This is related to the internal data
+                                 * structures used in this class. The raw data
+                                 * after a call to @p evaluate only
+                                 * contains unit cell operations, so possible
+                                 * transformations, quadrature weights
+                                 * etc. must be applied manually. In general,
+                                 * it is safer to use the get_laplacian() or
+                                 * get_hessian() functions instead, which does
+                                 * all the transformation internally.
+                                 */
+  const VectorizedArray<Number> * begin_hessians () const;
+
+                                /**
+                                 * Returns a read and write pointer to the
+                                 * first field of function hessians on
+                                 * quadrature points. First comes the
+                                 * xx-component of the hessian for the first
+                                 * component on all quadrature points, then
+                                 * the yy-component, zz-component in (3D),
+                                 * then the xy-component, and so on. Next
+                                 * comes the xx-component of the second
+                                 * component, and so on. This is related to
+                                 * the internal data structures used in this
+                                 * class. The raw data after a call to @p
+                                 * evaluate only contains unit cell
+                                 * operations, so possible transformations,
+                                 * quadrature weights etc. must be applied
+                                 * manually. In general, it is safer to use
+                                 * the get_laplacian() or get_hessian()
+                                 * functions instead, which does all the
+                                 * transformation internally.
+                                 */
+  VectorizedArray<Number> * begin_hessians ();
+
+                                //@}
+
 protected:
 
                                 /**
@@ -585,6 +656,33 @@ protected:
                     const unsigned int            fe_no   = 0,
                     const unsigned int            quad_no = 0);
 
+                                /**
+                                 * 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, @p
+                                 * distribute_local_to_global, and @p
+                                 * set_dof_values. It performs the operation
+                                 * for several vectors at a time.
+                                 */
+  template<typename VectorType, typename VectorOperation>
+  void read_write_operation (const VectorOperation &operation,
+                             VectorType            *vectors[]) const;
+
+                                /**
+                                 * For a collection of several vector @p src,
+                                 * read out the values on the degrees of
+                                 * freedom of the current cell for @p
+                                 * n_components (template argument), and store
+                                 * them internally. Similar functionality as
+                                 * the function
+                                 * DoFAccessor::read_dof_values. Note
+                                 * that if vectorization is enabled, the DoF
+                                 * values for several cells are set.
+                                 */
+  template<typename VectorType>
+  void read_dof_values_plain (const VectorType * src_data[]);
+
                                 /**
                                  * Internal data fields that store the
                                  * values. Since all array lengths are known
@@ -698,52 +796,6 @@ protected:
                                  */
   const internal::MatrixFreeFunctions::ShapeInfo<Number> &data;
 
-                                /**
-                                 * After a call to reinit(), stores the number
-                                 * of the cell we are currently working with.
-                                 */
-  unsigned int cell;
-
-                                /**
-                                 * Stores the type of the cell we are
-                                 * currently working with after a call to
-                                 * reinit(). Valid values are @p cartesian, @p
-                                 * affine and @p general, which have different
-                                 * implications on how the Jacobian
-                                 * transformations are stored internally in
-                                 * MappingInfo.
-                                 */
-  internal::MatrixFreeFunctions::CellType cell_type;
-
-                                /**
-                                 * The stride to access the correct data in
-                                 * MappingInfo.
-                                 */
-  unsigned int cell_data_number;
-
-                                /**
-                                 * Stores whether the present cell chunk used
-                                 * in vectorization is not completely filled
-                                 * up with physical cells. E.g. if
-                                 * vectorization dictates that four cells
-                                 * should be worked with but only three
-                                 * physical cells are left, this flag will be
-                                 * set to true, otherwise to false. Mainly
-                                 * used for internal checking when reading
-                                 * from vectors or writing to vectors.
-                                 */
-  bool         at_irregular_cell;
-
-                                /**
-                                 * If the present cell chunk for vectorization
-                                 * is not completely filled up with data, this
-                                 * field stores how many physical cells are
-                                 * underlying. Is between 1 and
-                                 * VectorizedArray<Number>::n_array_elements-1
-                                 * (inclusive).
-                                 */
-  unsigned int n_irreg_components_filled;
-
                                 /**
                                  * A pointer to the Cartesian Jacobian
                                  * information of the present cell. Only set
@@ -797,13 +849,59 @@ protected:
                                  */
   const Tensor<1,(dim>1?dim*(dim-1)/2:1),Tensor<1,dim,VectorizedArray<Number> > > * jacobian_grad_upper;
 
+                                /**
+                                 * After a call to reinit(), stores the number
+                                 * of the cell we are currently working with.
+                                 */
+  unsigned int cell;
+
+                                /**
+                                 * Stores the type of the cell we are
+                                 * currently working with after a call to
+                                 * reinit(). Valid values are @p cartesian, @p
+                                 * affine and @p general, which have different
+                                 * implications on how the Jacobian
+                                 * transformations are stored internally in
+                                 * MappingInfo.
+                                 */
+  internal::MatrixFreeFunctions::CellType cell_type;
+
+                                /**
+                                 * The stride to access the correct data in
+                                 * MappingInfo.
+                                 */
+  unsigned int cell_data_number;
+
+                                /**
+                                 * If the present cell chunk for vectorization
+                                 * is not completely filled up with data, this
+                                 * field stores how many physical cells are
+                                 * underlying. Is between 1 and
+                                 * VectorizedArray<Number>::n_array_elements-1
+                                 * (inclusive).
+                                 */
+  unsigned int n_irreg_components_filled;
+
+                                /**
+                                 * Stores whether the present cell chunk used
+                                 * in vectorization is not completely filled
+                                 * up with physical cells. E.g. if
+                                 * vectorization dictates that four cells
+                                 * should be worked with but only three
+                                 * physical cells are left, this flag will be
+                                 * set to true, otherwise to false. Mainly
+                                 * used for internal checking when reading
+                                 * from vectors or writing to vectors.
+                                 */
+  bool at_irregular_cell;
+
                                 /**
                                  * Debug information to track whether dof
                                  * values have been initialized before
                                  * accessed. Used to control exceptions when
                                  * uninitialized data is used.
                                  */
-  bool     dof_values_initialized;
+  bool dof_values_initialized;
 
                                 /**
                                  * Debug information to track whether values
@@ -811,7 +909,7 @@ protected:
                                  * before accessed. Used to control exceptions
                                  * when uninitialized data is used.
                                  */
-  bool     values_quad_initialized;
+  bool values_quad_initialized;
 
                                 /**
                                  * Debug information to track whether
@@ -820,7 +918,7 @@ protected:
                                  * control exceptions when uninitialized data
                                  * is used.
                                  */
-  bool     gradients_quad_initialized;
+  bool gradients_quad_initialized;
 
                                 /**
                                  * Debug information to track whether
@@ -829,7 +927,7 @@ protected:
                                  * control exceptions when uninitialized data
                                  * is used.
                                  */
-  bool     hessians_quad_initialized;
+  bool hessians_quad_initialized;
 
                                 /**
                                  * Debug information to track whether values
@@ -838,7 +936,7 @@ protected:
                                  * actually stared. Used to control exceptions
                                  * when uninitialized data is used.
                                  */
-  bool     values_quad_submitted;
+  bool values_quad_submitted;
 
                                 /**
                                  * Debug information to track whether
@@ -848,7 +946,7 @@ protected:
                                  * control exceptions when uninitialized data
                                  * is used.
                                  */
-  bool     gradients_quad_submitted;
+  bool gradients_quad_submitted;
 };
 
 
@@ -862,17 +960,19 @@ protected:
  * @author Katharina Kormann and Martin Kronbichler, 2010, 2011
  */
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 class FEEvaluationAccess :
-  public FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>
+  public FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
 {
 public:
-  typedef Tensor<1,n_components,VectorizedArray<Number> > value_type;
-  typedef Tensor<1,n_components,Tensor<1,dim,VectorizedArray<Number> > > gradient_type;
+  typedef Number                            number_type;
+  typedef Tensor<1,n_components_,VectorizedArray<Number> > value_type;
+  typedef Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > > gradient_type;
   static const unsigned int dimension     = dim;
+  static const unsigned int n_components  = n_components_;
   static const unsigned int dofs_per_cell = dofs_per_cell_;
   static const unsigned int n_q_points    = n_q_points_;
-  typedef FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,
+  typedef FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,
                            Number> BaseClass;
 
 protected:
@@ -907,6 +1007,7 @@ class FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number> :
   public FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,1,Number>
 {
  public:
+  typedef Number                                 number_type;
   typedef VectorizedArray<Number>                value_type;
   typedef Tensor<1,dim,VectorizedArray<Number> > gradient_type;
   static const unsigned int dimension          = dim;
@@ -1072,9 +1173,11 @@ class FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number> :
   public FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,dim,Number>
 {
  public:
-  typedef Tensor<1,dim,VectorizedArray<Number> >           value_type;
-  typedef Tensor<2,dim,VectorizedArray<Number> >           gradient_type;
+  typedef Number                            number_type;
+  typedef Tensor<1,dim,VectorizedArray<Number> > value_type;
+  typedef Tensor<2,dim,VectorizedArray<Number> > gradient_type;
   static const unsigned int dimension     = dim;
+  static const unsigned int n_components  = dim;
   static const unsigned int dofs_per_cell = dofs_per_cell_;
   static const unsigned int n_q_points    = n_q_points_;
   typedef FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,dim,Number> BaseClass;
@@ -1246,21 +1349,24 @@ protected:
  *
  * @author Katharina Kormann and Martin Kronbichler, 2010, 2011
  */
-template <int dim, int fe_degree, int n_q_points_1d=fe_degree+1,
-          int n_components=1, typename Number=double >
+template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
+          int n_components_ = 1, typename Number = double >
 class FEEvaluationGeneral :
   public FEEvaluationAccess<dim,
                             Utilities::fixed_int_power<fe_degree+1,dim>::value,
                             Utilities::fixed_int_power<n_q_points_1d,dim>::value,
-                            n_components,Number>
+                            n_components_,Number>
 {
  public:
   typedef FEEvaluationAccess<dim,
                              Utilities::fixed_int_power<fe_degree+1,dim>::value,
                              Utilities::fixed_int_power<n_q_points_1d,dim>::value,
-                             n_components, Number> BaseClass;
+                             n_components_, Number> BaseClass;
+  typedef Number                            number_type;
   typedef typename BaseClass::value_type    value_type;
   typedef typename BaseClass::gradient_type gradient_type;
+  static const unsigned int dimension     = dim;
+  static const unsigned int n_components  = n_components_;
   static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell;
   static const unsigned int n_q_points    = BaseClass::n_q_points;
 
@@ -1288,8 +1394,9 @@ class FEEvaluationGeneral :
                                  * @p get_gradient() or @p get_laplacian
                                  * return useful information.
                                  */
-  void evaluate (bool evaluate_val, bool evaluate_grad,
-                 bool evaluate_hess=false);
+  void evaluate (const bool evaluate_val,
+                 const bool evaluate_grad,
+                 const bool evaluate_hess = false);
 
                                 /**
                                  * This function takes the values and/or
@@ -1301,7 +1408,8 @@ class FEEvaluationGeneral :
                                  * integrate_grad are used to enable/disable
                                  * some of values or gradients.
                                  */
-  void integrate (bool integrate_val, bool integrate_grad);
+  void integrate (const bool integrate_val,
+                  const bool integrate_grad);
 
                                 /**
                                  * Returns the q-th quadrature point stored in
@@ -1313,21 +1421,63 @@ class FEEvaluationGeneral :
 protected:
 
                                 /**
-                                 * Internal function that applies the shape
-                                 * function data of the tensor product in a
-                                 * given coordinate direction (first template
+                                 * Internal function that applies the function
+                                 * values of the tensor product in a given
+                                 * coordinate direction (first template
+                                 * argument), from polynomials to values on
+                                 * quadrature points (second flag set to true)
+                                 * or in an integration loop from values on
+                                 * quadrature points to values tested by
+                                 * different test function (second flag set to
+                                 * false), and if the result is to be added to
+                                 * previous content in the data fields or
+                                 * not.
+                                 */
+  template <int direction, bool dof_to_quad, bool add>
+  void apply_values (const VectorizedArray<Number> in [],
+                     VectorizedArray<Number> out []);
+
+                                /**
+                                 * Internal function that applies the gradient
+                                 * operation of the tensor product in a given
+                                 * coordinate direction (first template
                                  * argument), from polynomials to values on
                                  * quadrature points (second flag set to true)
                                  * or in an integration loop from values on
                                  * quadrature points to values tested by
                                  * different test function (second flag set to
                                  * false), and if the result is to be added to
-                                 * some previous results or not.
+                                 * previous content in the data fields or
+                                 * not.
+                                 */
+  template <int direction, bool dof_to_quad, bool add>
+  void apply_gradients (const VectorizedArray<Number> in [],
+                        VectorizedArray<Number> out []);
+
+                                /**
+                                 * Internal function that applies the second
+                                 * derivative operation (Hessian) of the
+                                 * tensor product in a given coordinate
+                                 * direction (first template argument), from
+                                 * polynomials to values on quadrature points
+                                 * (second flag set to true) or in an
+                                 * integration loop from values on quadrature
+                                 * points to values tested by different test
+                                 * function (second flag set to false), and if
+                                 * the result is to be added to previous
+                                 * content in the data fields or not.
                                  */
   template <int direction, bool dof_to_quad, bool add>
-  void apply_tensor_prod (const VectorizedArray<Number> * shape_data,
-                          const VectorizedArray<Number> in [],
-                          VectorizedArray<Number> out []);
+  void apply_hessians (const VectorizedArray<Number> in [],
+                       VectorizedArray<Number> out []);
+
+                                /**
+                                 * Friend declaration.
+                                 */
+  template <typename FEEval> friend void
+  internal::do_evaluate (FEEval &, const bool, const bool, const bool);
+  template <typename FEEval> friend void
+  internal::do_integrate (FEEval &, const bool, const bool);
 };
 
 
@@ -1373,15 +1523,18 @@ protected:
  *
  * @author Katharina Kormann and Martin Kronbichler, 2010, 2011
  */
-template <int dim, int fe_degree, int n_q_points_1d=fe_degree+1,
-          int n_components=1, typename Number=double >
+template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
+          int n_components_ = 1, typename Number = double >
 class FEEvaluation :
-  public FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components,Number>
+  public FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
 {
  public:
-  typedef FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components,Number> BaseClass;
+  typedef FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number> BaseClass;
+  typedef Number                            number_type;
   typedef typename BaseClass::value_type    value_type;
   typedef typename BaseClass::gradient_type gradient_type;
+  static const unsigned int dimension     = dim;
+  static const unsigned int n_components  = n_components_;
   static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell;
   static const unsigned int n_q_points    = BaseClass::n_q_points;
 
@@ -1411,8 +1564,9 @@ class FEEvaluation :
                                  * (unless these values have been set
                                  * manually).
                                  */
-  void evaluate (bool evaluate_val, bool evaluate_grad,
-                 bool evaluate_hess=false);
+  void evaluate (const bool evaluate_val,
+                 const bool evaluate_grad,
+                 const bool evaluate_hess = false);
 
                                 /**
                                  * This function takes the values and/or
@@ -1424,9 +1578,11 @@ class FEEvaluation :
                                  * integrate_grad are used to enable/disable
                                  * some of values or gradients.
                                  */
-  void integrate (bool integrate_val, bool integrate_grad);
+  void integrate (const bool integrate_val,
+                  const bool integrate_grad);
 
 protected:
+
                                 /**
                                  * Internal function that applies the function
                                  * values of the tensor product in a given
@@ -1477,6 +1633,14 @@ protected:
   template <int direction, bool dof_to_quad, bool add>
   void apply_hessians (const VectorizedArray<Number> in [],
                        VectorizedArray<Number> out []);
+
+                                /**
+                                 * Friend declarations.
+                                 */
+  template <typename FEEval> friend void
+  internal::do_evaluate (FEEval &, const bool, const bool, const bool);
+  template <typename FEEval> friend void
+  internal::do_integrate (FEEval &, const bool, const bool);
 };
 
 
@@ -1518,14 +1682,17 @@ protected:
  *
  * @author Katharina Kormann and Martin Kronbichler, 2010, 2011
  */
-template <int dim, int fe_degree, int n_components=1, typename Number=double >
+template <int dim, int fe_degree, int n_components_ = 1, typename Number = double >
 class FEEvaluationGL :
-  public FEEvaluation<dim,fe_degree,fe_degree+1,n_components,Number>
+  public FEEvaluation<dim,fe_degree,fe_degree+1,n_components_,Number>
 {
  public:
-  typedef FEEvaluation<dim,fe_degree,fe_degree+1,n_components,Number> BaseClass;
+  typedef FEEvaluation<dim,fe_degree,fe_degree+1,n_components_,Number> BaseClass;
+  typedef Number                            number_type;
   typedef typename BaseClass::value_type    value_type;
   typedef typename BaseClass::gradient_type gradient_type;
+  static const unsigned int dimension     = dim;
+  static const unsigned int n_components  = n_components_;
   static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell;
   static const unsigned int n_q_points    = BaseClass::n_q_points;
 
@@ -1555,8 +1722,9 @@ class FEEvaluationGL :
                                  * (unless these values have been set
                                  * manually).
                                  */
-  void evaluate (bool evaluate_val, bool evaluate_grad,
-                 bool evaluate_lapl=false);
+  void evaluate (const bool evaluate_val,
+                 const bool evaluate_grad,
+                 const bool evaluate_lapl = false);
 
                                 /**
                                  * This function takes the values and/or
@@ -1568,7 +1736,8 @@ class FEEvaluationGL :
                                  * integrate_grad are used to enable/disable
                                  * some of values or gradients.
                                  */
-  void integrate (bool integrate_val, bool integrate_grad);
+  void integrate (const bool integrate_val,
+                  const bool integrate_grad);
 
 protected:
                                 /**
@@ -1596,15 +1765,15 @@ protected:
 #ifndef DOXYGEN
 
 
-/*----------------------- FEEvaluationBase -------------------------------*/
+/*----------------------- FEEvaluationBase ----------------------------------*/
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 inline
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-FEEvaluationBase (const MatrixFree<dim,Number> &data_in,
-                  const unsigned int fe_no_in,
-                  const unsigned int quad_no_in)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::FEEvaluationBase (const MatrixFree<dim,Number> &data_in,
+                    const unsigned int fe_no_in,
+                    const unsigned int quad_no_in)
   :
   quad_no            (quad_no_in),
   n_fe_components    (data_in.get_dof_info(fe_no_in).n_components),
@@ -1619,11 +1788,6 @@ FEEvaluationBase (const MatrixFree<dim,Number> &data_in,
   data               (data_in.get_shape_info
                       (fe_no_in, quad_no_in, active_fe_index,
                        active_quad_index)),
-  cell               (numbers::invalid_unsigned_int),
-  cell_type          (internal::MatrixFreeFunctions::undefined),
-  cell_data_number   (0),
-  at_irregular_cell  (false),
-  n_irreg_components_filled (0),
   cartesian_data     (0),
   jacobian           (0),
   J_value            (0),
@@ -1631,7 +1795,12 @@ FEEvaluationBase (const MatrixFree<dim,Number> &data_in,
                       quadrature_weights[active_quad_index].begin()),
   quadrature_points  (0),
   jacobian_grad      (0),
-  jacobian_grad_upper(0)
+  jacobian_grad_upper(0),
+  cell               (numbers::invalid_unsigned_int),
+  cell_type          (internal::MatrixFreeFunctions::undefined),
+  cell_data_number   (0),
+  n_irreg_components_filled (0),
+  at_irregular_cell  (false)
 {
   Assert (matrix_info.indices_initialized() == true,
           ExcNotInitialized());
@@ -1654,11 +1823,11 @@ FEEvaluationBase (const MatrixFree<dim,Number> &data_in,
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 inline
 void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-reinit (const unsigned int cell_in)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::reinit (const unsigned int cell_in)
 {
   AssertIndexRange (cell_in, dof_info.row_starts.size()-1);
   AssertDimension (((dof_info.cell_active_fe_index.size() > 0) ?
@@ -1732,11 +1901,11 @@ reinit (const unsigned int cell_in)
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 inline
 unsigned int
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_cell_data_number () const
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_cell_data_number () const
 {
   Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
   return cell_data_number;
@@ -1745,11 +1914,11 @@ get_cell_data_number () const
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 inline
 internal::MatrixFreeFunctions::CellType
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_cell_type () const
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_cell_type () const
 {
   Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
   return cell_type;
@@ -1839,82 +2008,148 @@ namespace internal
                        "Use MatrixFree::initialize_dof_vector to get a "
                        "compatible vector."));
   }
-}
 
+                                // A class to use the same code to read from
+                                // and write to vector
+  template <typename Number>
+  struct VectorReader
+  {
+    template <typename VectorType>
+    void process_dof (const unsigned int  index,
+                      VectorType         &vec,
+                      Number             &res) const
+    {
+      res = vector_access (vec, index);
+    }
 
+    void pre_constraints (const Number &,
+                          Number       &res) const
+    {
+      res = Number();
+    }
 
-template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
-template<typename VectorType>
-inline
-void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values (const VectorType &src)
-{
-  AssertDimension (n_components, n_fe_components);
-                                // only need one component, but to avoid
-                                // compiler warnings, use n_components copies
-                                // here (but these will not be used)
-  const VectorType * src_data[n_components];
-  for (unsigned int d=0; d<n_components; ++d)
-    src_data[d] = &src;
-  read_dof_values (src_data);
-}
+    template <typename VectorType>
+    void process_constraint (const unsigned int index,
+                             const Number       weight,
+                             VectorType        &vec,
+                             Number            &res) const
+    {
+      res += weight * vector_access (vec, index);
+    }
 
+    void post_constraints (const Number &sum,
+                           Number       &write_pos) const
+    {
+      write_pos = sum;
+    }
 
+    void process_empty (Number &res) const
+    {
+      res = Number();
+    }
+  };
 
-template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
-template<typename VectorType>
-inline
-void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values (const std::vector<VectorType> &src,
-                const unsigned int             first_index)
-{
-  AssertIndexRange (first_index, src.size());
-  Assert (n_fe_components == 1, ExcNotImplemented());
-  Assert ((n_fe_components == 1 ?
-           (first_index+n_components <= src.size()) : true),
-          ExcIndexRange (first_index + n_components, 0, src.size()));
-  const VectorType * src_data [n_components];
-  for (unsigned int comp=0; comp<n_components; ++comp)
-    src_data[comp] = &src[comp+first_index];
-  read_dof_values (src_data);
-}
+                                // A class to use the same code to read from
+                                // and write to vector
+  template <typename Number>
+  struct VectorDistributorLocalToGlobal
+  {
+    template <typename VectorType>
+    void process_dof (const unsigned int  index,
+                      VectorType         &vec,
+                      Number             &res) const
+    {
+      vector_access (vec, index) += res;
+    }
 
+    void pre_constraints (const Number &input,
+                          Number       &res) const
+    {
+      res = input;
+    }
 
+    template <typename VectorType>
+    void process_constraint (const unsigned int index,
+                             const Number       weight,
+                             VectorType        &vec,
+                             Number            &res) const
+    {
+      vector_access (vec, index) += weight * res;
+    }
+
+    void post_constraints (const Number &,
+                           Number       &) const
+    {
+    }
+
+    void process_empty (Number &) const
+    {
+    }
+  };
+
+
+                                // A class to use the same code to read from
+                                // and write to vector
+  template <typename Number>
+  struct VectorSetter
+  {
+    template <typename VectorType>
+    void process_dof (const unsigned int  index,
+                      VectorType         &vec,
+                      Number             &res) const
+    {
+      vector_access (vec, index) = res;
+    }
+
+    void pre_constraints (const Number &,
+                          Number &) const
+    {
+    }
+
+    template <typename VectorType>
+    void process_constraint (const unsigned int,
+                             const Number,
+                             VectorType&,
+                             Number &) const
+    {
+    }
+
+    void post_constraints (const Number &,
+                           Number       &) const
+    {
+    }
+
+    void process_empty (Number &) const
+    {
+    }
+  };
 
-template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
-template<typename VectorType>
-inline
-void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values (const std::vector<VectorType*> &src,
-                const unsigned int              first_index)
-{
-  AssertIndexRange (first_index, src.size());
-  Assert (n_fe_components == 1, ExcNotImplemented());
-  Assert ((n_fe_components == 1 ?
-           (first_index+n_components <= src.size()) : true),
-          ExcIndexRange (first_index + n_components, 0, src.size()));
-  const VectorType * src_data [n_components];
-  for (unsigned int comp=0; comp<n_components; ++comp)
-    src_data[comp] = src[comp+first_index];
-  read_dof_values (src_data);
 }
 
 
 
+
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
-template<typename VectorType>
+          int n_components_, typename Number>
+template<typename VectorType, typename VectorOperation>
 inline
 void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values (const VectorType * src[])
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_write_operation (const VectorOperation &operation,
+                        VectorType            *src[]) const
 {
+                                // This functions processes all the functions
+                                // read_dof_values,
+                                // distribute_local_to_global, and
+                                // set_dof_values with the same code. The
+                                // distinction between these three cases is
+                                // made by the input VectorOperation that
+                                // either reads values from a vector and puts
+                                // the data into the local data field or write
+                                // local data into the vector. Certain
+                                // operations are no-ops for the given use
+                                // case.
+
   Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
 
                                 // loop over all local dofs. ind_local holds
@@ -1938,9 +2173,10 @@ read_dof_values (const VectorType * src[])
         VectorizedArray<Number>::n_array_elements * dofs_per_cell;
       for (unsigned int comp=0; comp<n_components; ++comp)
         internal::check_vector_compatibility (*src[comp], dof_info);
-      Number * local_src_number [n_components];
+      Number * local_data [n_components];
       for (unsigned int comp=0; comp<n_components; ++comp)
-        local_src_number[comp] = reinterpret_cast<Number*>(values_dofs[comp]);
+        local_data[comp] =
+          const_cast<Number*>(reinterpret_cast<const Number*>(values_dofs[comp]));
 
                                 // standard case where there are sufficiently
                                 // many cells to fill all vectors
@@ -1955,8 +2191,9 @@ read_dof_values (const VectorType * src[])
                                 // run through values up to next constraint
                   for (unsigned int j=0; j<indicators->first; ++j)
                     for (unsigned int comp=0; comp<n_components; ++comp)
-                      local_src_number[comp][ind_local+j] =
-                        internal::vector_access (*src[comp], dof_indices[j]);
+                      operation.process_dof (dof_indices[j], *src[comp],
+                                             local_data[comp][ind_local+j]);
+
                   ind_local += indicators->first;
                   dof_indices   += indicators->first;
 
@@ -1965,19 +2202,22 @@ read_dof_values (const VectorType * src[])
                                 // according to constraints
                   Number value [n_components];
                   for (unsigned int comp=0; comp<n_components; ++comp)
-                    value[comp] = 0;
+                    operation.pre_constraints (local_data[comp][ind_local],
+                                               value[comp]);
+
                   const Number * data_val =
                     matrix_info.constraint_pool_begin(indicators->second);
                   const Number * end_pool =
                     matrix_info.constraint_pool_end(indicators->second);
                   for ( ; data_val != end_pool; ++data_val, ++dof_indices)
                     for (unsigned int comp=0; comp<n_components; ++comp)
-                      value[comp] +=
-                        (internal::vector_access (*src[comp], *dof_indices) *
-                         (*data_val));
+                      operation.process_constraint (*dof_indices, *data_val,
+                                                    *src[comp], value[comp]);
 
                   for (unsigned int comp=0; comp<n_components; ++comp)
-                    local_src_number[comp][ind_local] = value[comp];
+                    operation.post_constraints (value[comp],
+                                                local_data[comp][ind_local]);
+
                   ind_local++;
                 }
 
@@ -1986,8 +2226,8 @@ read_dof_values (const VectorType * src[])
               for(; ind_local < n_local_dofs; ++dof_indices, ++ind_local)
                 {
                   for (unsigned int comp=0; comp<n_components; ++comp)
-                    local_src_number[comp][ind_local] =
-                      internal::vector_access (*src[comp], *dof_indices);
+                      operation.process_dof (*dof_indices, *src[comp],
+                                             local_data[comp][ind_local]);
                 }
             }
           else
@@ -1998,8 +2238,8 @@ read_dof_values (const VectorType * src[])
                                static_cast<int>(n_local_dofs));
               for (unsigned int j=0; j<n_local_dofs; ++j)
                 for (unsigned int comp=0; comp<n_components; ++comp)
-                  local_src_number[comp][j] =
-                    internal::vector_access (*src[comp], dof_indices[j]);
+                  operation.process_dof (dof_indices[j], *src[comp],
+                                         local_data[comp][j]);
             }
         }
 
@@ -2018,8 +2258,8 @@ read_dof_values (const VectorType * src[])
                                 // the global vector, src, to the local one,
                                 // local_src.
                   for (unsigned int comp=0; comp<n_components; ++comp)
-                    local_src_number[comp][ind_local] =
-                      internal::vector_access (*src[comp], dof_indices[j]);
+                    operation.process_dof (dof_indices[j], *src[comp],
+                                           local_data[comp][ind_local]);
 
                                 // here we jump over all the components that
                                 // are artificial
@@ -2028,7 +2268,7 @@ read_dof_values (const VectorType * src[])
                          >= n_irreg_components_filled)
                     {
                       for (unsigned int comp=0; comp<n_components; ++comp)
-                        local_src_number[comp][ind_local] = 0.;
+                        operation.process_empty (local_data[comp][ind_local]);
                       ++ind_local;
                     }
                 }
@@ -2039,23 +2279,28 @@ read_dof_values (const VectorType * src[])
                                 // according to constraint
               Number value [n_components];
               for (unsigned int comp=0; comp<n_components; ++comp)
-                value[comp] = 0.;
+                operation.pre_constraints (local_data[comp][ind_local],
+                                           value[comp]);
+
               const Number * data_val =
                 matrix_info.constraint_pool_begin(indicators->second);
               const Number * end_pool =
                 matrix_info.constraint_pool_end(indicators->second);
+
               for ( ; data_val != end_pool; ++data_val, ++dof_indices)
                 for (unsigned int comp=0; comp<n_components; ++comp)
-                  value[comp] +=
-                    internal::vector_access (*src[comp], *dof_indices) * (*data_val);
+                  operation.process_constraint (*dof_indices, *data_val,
+                                                *src[comp], value[comp]);
+
               for (unsigned int comp=0; comp<n_components; ++comp)
-                local_src_number[comp][ind_local] = value[comp];
+                operation.post_constraints (value[comp],
+                                            local_data[comp][ind_local]);
               ind_local++;
               while (ind_local % VectorizedArray<Number>::n_array_elements
                      >= n_irreg_components_filled)
                 {
                   for (unsigned int comp=0; comp<n_components; ++comp)
-                    local_src_number[comp][ind_local] = 0.;
+                    operation.process_empty (local_data[comp][ind_local]);
                   ++ind_local;
                 }
             }
@@ -2068,14 +2313,14 @@ read_dof_values (const VectorType * src[])
                                 // the global vector, src, to the local one,
                                 // local_dst.
               for (unsigned int comp=0; comp<n_components; ++comp)
-                local_src_number[comp][ind_local] =
-                  internal::vector_access (*src[comp], *dof_indices);
+                operation.process_dof (*dof_indices, *src[comp],
+                                       local_data[comp][ind_local]);
               ++ind_local;
               while (ind_local % VectorizedArray<Number>::n_array_elements
                      >= n_irreg_components_filled)
                 {
                   for (unsigned int comp=0; comp<n_components; ++comp)
-                    local_src_number[comp][ind_local] = 0.;
+                    operation.process_empty(local_data[comp][ind_local]);
                   ++ind_local;
                 }
             }
@@ -2091,10 +2336,11 @@ read_dof_values (const VectorType * src[])
                                 // out the indices.
     {
       internal::check_vector_compatibility (*src[0], dof_info);
-      Assert (n_fe_components == n_components, ExcNotImplemented());
+      Assert (n_fe_components == n_components_, ExcNotImplemented());
       const unsigned int n_local_dofs =
         dofs_per_cell*VectorizedArray<Number>::n_array_elements * n_components;
-      Number * local_src_number = reinterpret_cast<Number*>(values_dofs[0]);
+      Number *  local_data =
+          const_cast<Number*>(reinterpret_cast<const Number*>(values_dofs[0]));
       if (at_irregular_cell == false)
         {
                                 // check whether there is any constraint on
@@ -2105,33 +2351,35 @@ read_dof_values (const VectorType * src[])
                 {
                                 // run through values up to next constraint
                   for (unsigned int j=0; j<indicators->first; ++j)
-                      local_src_number[ind_local+j] =
-                        internal::vector_access (*src[0], dof_indices[j]);
+                    operation.process_dof (dof_indices[j], *src[0],
+                                           local_data[ind_local+j]);
                   ind_local += indicators->first;
                   dof_indices   += indicators->first;
 
                                 // constrained case: build the local value as
                                 // a linear combination of the global value
                                 // according to constraints
-                  Number value = 0;
+                  Number value;
+                  operation.pre_constraints (local_data[ind_local], value);
+
                   const Number * data_val =
                     matrix_info.constraint_pool_begin(indicators->second);
                   const Number * end_pool =
                     matrix_info.constraint_pool_end(indicators->second);
+
                   for ( ; data_val != end_pool; ++data_val, ++dof_indices)
-                    value +=
-                      (internal::vector_access (*src[0], *dof_indices) *
-                       (*data_val));
+                    operation.process_constraint (*dof_indices, *data_val,
+                                                  *src[0], value);
 
-                  local_src_number[ind_local] = value;
+                  operation.post_constraints (value, local_data[ind_local]);
                   ind_local++;
                 }
 
                                 // get the dof values past the last
                                 // constraint
               for(; ind_local<n_local_dofs; ++dof_indices, ++ind_local)
-                local_src_number[ind_local] =
-                  internal::vector_access (*src[0], *dof_indices);
+                operation.process_dof (*dof_indices, *src[0],
+                                       local_data[ind_local]);
               Assert (dof_indices == dof_info.end_indices(cell),
                       ExcInternalError());
             }
@@ -2142,8 +2390,8 @@ read_dof_values (const VectorType * src[])
               AssertDimension (dof_info.end_indices(cell)-dof_indices,
                                static_cast<int>(n_local_dofs));
               for (unsigned int j=0; j<n_local_dofs; ++j)
-                local_src_number[j] =
-                  internal::vector_access (*src[0], dof_indices[j]);
+                operation.process_dof (dof_indices[j], *src[0],
+                                       local_data[j]);
             }
         }
 
@@ -2161,8 +2409,8 @@ read_dof_values (const VectorType * src[])
                                 // non-constrained case: copy the data from
                                 // the global vector, src, to the local one,
                                 // local_src.
-                  local_src_number[ind_local] =
-                    internal::vector_access (*src[0], dof_indices[j]);
+                  operation.process_dof (dof_indices[j], *src[0],
+                                         local_data[ind_local]);
 
                                 // here we jump over all the components that
                                 // are artificial
@@ -2170,7 +2418,7 @@ read_dof_values (const VectorType * src[])
                   while (ind_local % VectorizedArray<Number>::n_array_elements
                          >= n_irreg_components_filled)
                     {
-                      local_src_number[ind_local] = 0.;
+                      operation.process_empty (local_data[ind_local]);
                       ++ind_local;
                     }
                 }
@@ -2179,20 +2427,24 @@ read_dof_values (const VectorType * src[])
                                 // constrained case: build the local value as
                                 // a linear combination of the global value
                                 // according to constraint
-              Number value = 0;
+              Number value;
+              operation.pre_constraints (local_data[ind_local], value);
+
               const Number * data_val =
                 matrix_info.constraint_pool_begin(indicators->second);
               const Number * end_pool =
                 matrix_info.constraint_pool_end(indicators->second);
+
               for ( ; data_val != end_pool; ++data_val, ++dof_indices)
-                value +=
-                  internal::vector_access (*src[0], *dof_indices) * (*data_val);
-              local_src_number[ind_local] = value;
+                operation.process_constraint (*dof_indices, *data_val,
+                                              *src[0], value);
+
+              operation.post_constraints (value, local_data[ind_local]);
               ind_local++;
               while (ind_local % VectorizedArray<Number>::n_array_elements
                      >= n_irreg_components_filled)
                 {
-                  local_src_number[ind_local] = 0.;
+                  operation.process_empty (local_data[ind_local]);
                   ++ind_local;
                 }
             }
@@ -2204,18 +2456,41 @@ read_dof_values (const VectorType * src[])
                                 // non-constrained case: copy the data from
                                 // the global vector, src, to the local one,
                                 // local_dst.
-              local_src_number[ind_local] =
-                internal::vector_access (*src[0], *dof_indices);
+              operation.process_dof (*dof_indices, *src[0],
+                                     local_data[ind_local]);
               ++ind_local;
               while (ind_local % VectorizedArray<Number>::n_array_elements
                      >= n_irreg_components_filled)
                 {
-                  local_src_number[ind_local] = 0.;
+                  operation.process_empty (local_data[ind_local]);
                   ++ind_local;
                 }
             }
         }
     }
+}
+
+
+
+template <int dim, int dofs_per_cell_, int n_q_points_,
+          int n_components_, typename Number>
+template<typename VectorType>
+inline
+void
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_dof_values (const VectorType &src)
+{
+  AssertDimension (n_components_, n_fe_components);
+
+                                // only need one component, but to silent
+                                // compiler warnings, use n_components copies
+                                // here (but these will not be used)
+  VectorType * src_data[n_components];
+  for (unsigned int d=0; d<n_components; ++d)
+    src_data[d] = const_cast<VectorType*>(&src);
+
+  internal::VectorReader<Number> reader;
+  read_write_operation (reader, src_data);
 
 #ifdef DEBUG
   dof_values_initialized = true;
@@ -2224,18 +2499,73 @@ read_dof_values (const VectorType * src[])
 
 
 
+template <int dim, int dofs_per_cell_, int n_q_points_,
+          int n_components_, typename Number>
+template<typename VectorType>
+inline
+void
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_dof_values (const std::vector<VectorType> &src,
+                   const unsigned int             first_index)
+{
+  AssertIndexRange (first_index, src.size());
+  Assert (n_fe_components == 1, ExcNotImplemented());
+  Assert ((n_fe_components == 1 ?
+           (first_index+n_components <= src.size()) : true),
+          ExcIndexRange (first_index + n_components_, 0, src.size()));
+
+  VectorType * src_data [n_components];
+  for (unsigned int comp=0; comp<n_components; ++comp)
+    src_data[comp] = const_cast<VectorType*>(&src[comp+first_index]);
+
+  internal::VectorReader<Number> reader;
+  read_write_operation (reader, src_data);
+
+#ifdef DEBUG
+  dof_values_initialized = true;
+#endif
+}
 
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 template<typename VectorType>
 inline
 void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values_plain (const VectorType &src)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_dof_values (const std::vector<VectorType*> &src,
+                   const unsigned int              first_index)
+{
+  AssertIndexRange (first_index, src.size());
+  Assert (n_fe_components == 1, ExcNotImplemented());
+  Assert ((n_fe_components == 1 ?
+           (first_index+n_components <= src.size()) : true),
+          ExcIndexRange (first_index + n_components_, 0, src.size()));
+
+  const VectorType * src_data [n_components];
+  for (unsigned int comp=0; comp<n_components; ++comp)
+    src_data[comp] = const_cast<VectorType*>(src[comp+first_index]);
+
+  internal::VectorReader<Number> reader;
+  read_write_operation (reader, src_data);
+
+#ifdef DEBUG
+  dof_values_initialized = true;
+#endif
+}
+
+
+
+template <int dim, int dofs_per_cell_, int n_q_points_,
+          int n_components_, typename Number>
+template<typename VectorType>
+inline
+void
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_dof_values_plain (const VectorType &src)
 {
-  AssertDimension (n_components, n_fe_components);
+  AssertDimension (n_components_, n_fe_components);
                                 // only need one component, but to avoid
                                 // compiler warnings, use n_components copies
                                 // here (but these will not be used)
@@ -2248,19 +2578,19 @@ read_dof_values_plain (const VectorType &src)
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 template<typename VectorType>
 inline
 void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values_plain (const std::vector<VectorType> &src,
-                      const unsigned int             first_index)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_dof_values_plain (const std::vector<VectorType> &src,
+                         const unsigned int             first_index)
 {
   AssertIndexRange (first_index, src.size());
   Assert (n_fe_components == 1, ExcNotImplemented());
   Assert ((n_fe_components == 1 ?
            (first_index+n_components <= src.size()) : true),
-          ExcIndexRange (first_index + n_components, 0, src.size()));
+          ExcIndexRange (first_index + n_components_, 0, src.size()));
   const VectorType * src_data [n_components];
   for (unsigned int comp=0; comp<n_components; ++comp)
     src_data[comp] = &src[comp+first_index];
@@ -2270,19 +2600,19 @@ read_dof_values_plain (const std::vector<VectorType> &src,
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 template<typename VectorType>
 inline
 void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values_plain (const std::vector<VectorType*> &src,
-                      const unsigned int              first_index)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_dof_values_plain (const std::vector<VectorType*> &src,
+                         const unsigned int              first_index)
 {
   AssertIndexRange (first_index, src.size());
   Assert (n_fe_components == 1, ExcNotImplemented());
   Assert ((n_fe_components == 1 ?
            (first_index+n_components <= src.size()) : true),
-          ExcIndexRange (first_index + n_components, 0, src.size()));
+          ExcIndexRange (first_index + n_components_, 0, src.size()));
   const VectorType * src_data [n_components];
   for (unsigned int comp=0; comp<n_components; ++comp)
     src_data[comp] = src[comp+first_index];
@@ -2292,213 +2622,186 @@ read_dof_values_plain (const std::vector<VectorType*> &src,
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 template<typename VectorType>
 inline
 void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values_plain (const VectorType * src[])
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::distribute_local_to_global (VectorType &dst) const
 {
-  Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
-  Assert (dof_info.store_plain_indices == true, ExcNotInitialized());
+  AssertDimension (n_components_, n_fe_components);
+  Assert (dof_values_initialized==true,
+          internal::ExcAccessToUninitializedField());
 
-                                // loop over all local dofs. ind_local holds
-                                // local number on cell, index iterates over
-                                // the elements of index_local_to_global and
-                                // dof_indices points to the global indices
-                                // stored in index_local_to_global
-  const unsigned int * dof_indices = dof_info.begin_indices_plain(cell);
+                                // only need one component, but to avoid
+                                // compiler warnings, use n_components copies
+                                // here (but these will not be used)
+  VectorType * dst_data [n_components];
+  for (unsigned int d=0; d<n_components; ++d)
+    dst_data[d] = &dst;
 
-                                // scalar case (or case when all components
-                                // have the same degrees of freedom and sit on
-                                // a different vector each)
-  if (n_fe_components == 1)
-    {
-      const unsigned int n_local_dofs =
-        VectorizedArray<Number>::n_array_elements * dofs_per_cell;
-      for (unsigned int comp=0; comp<n_components; ++comp)
-        internal::check_vector_compatibility (*src[comp], dof_info);
-      Number * local_src_number [n_components];
-      for (unsigned int comp=0; comp<n_components; ++comp)
-        local_src_number[comp] = reinterpret_cast<Number*>(values_dofs[comp]);
+  internal::VectorDistributorLocalToGlobal<Number> distributor;
+  read_write_operation (distributor, dst_data);
+}
 
-                                // standard case where there are sufficiently
-                                // many cells to fill all vectors
-      if (at_irregular_cell == false)
-        {
-          for (unsigned int j=0; j<n_local_dofs; ++j)
-            for (unsigned int comp=0; comp<n_components; ++comp)
-              local_src_number[comp][j] =
-                internal::vector_access (*src[comp], dof_indices[j]);
-        }
 
-                                // non-standard case: need to fill in zeros
-                                // for those components that are not present
-                                // (a bit more expensive), but there is not
-                                // more than one such cell
-      else
-        {
-          Assert (n_irreg_components_filled > 0, ExcInternalError());
-          for(unsigned int ind_local=0; ind_local<n_local_dofs;
-              ++dof_indices)
-            {
-                                // non-constrained case: copy the data from
-                                // the global vector, src, to the local one,
-                                // local_dst.
-              for (unsigned int comp=0; comp<n_components; ++comp)
-                local_src_number[comp][ind_local] =
-                  internal::vector_access (*src[comp], *dof_indices);
-              ++ind_local;
-              while (ind_local % VectorizedArray<Number>::n_array_elements >= n_irreg_components_filled)
-                {
-                  for (unsigned int comp=0; comp<n_components; ++comp)
-                    local_src_number[comp][ind_local] = 0.;
-                  ++ind_local;
-                }
-            }
-        }
-    }
-  else
-                                // case with vector-valued finite elements
-                                // where all components are included in one
-                                // single vector. Assumption: first come all
-                                // entries to the first component, then all
-                                // entries to the second one, and so on. This
-                                // is ensured by the way MatrixFree reads
-                                // out the indices.
-    {
-      internal::check_vector_compatibility (*src[0], dof_info);
-      Assert (n_fe_components == n_components, ExcNotImplemented());
-      const unsigned int n_local_dofs =
-        dofs_per_cell * VectorizedArray<Number>::n_array_elements * n_components;
-      Number * local_src_number = reinterpret_cast<Number*>(values_dofs[0]);
-      if (at_irregular_cell == false)
-        {
-          for (unsigned int j=0; j<n_local_dofs; ++j)
-            local_src_number[j] =
-              internal::vector_access (*src[0], dof_indices[j]);
-        }
 
-                                // non-standard case: need to fill in zeros
-                                // for those components that are not present
-                                // (a bit more expensive), but there is not
-                                // more than one such cell
-      else
-        {
-          Assert (n_irreg_components_filled > 0, ExcInternalError());
-          for(unsigned int ind_local=0; ind_local<n_local_dofs; ++dof_indices)
-            {
-                                // non-constrained case: copy the data from
-                                // the global vector, src, to the local one,
-                                // local_dst.
-              local_src_number[ind_local] =
-                internal::vector_access (*src[0], *dof_indices);
-              ++ind_local;
-              while (ind_local % VectorizedArray<Number>::n_array_elements >= n_irreg_components_filled)
-                {
-                  local_src_number[ind_local] = 0.;
-                  ++ind_local;
-                }
-            }
-        }
-    }
+template <int dim, int dofs_per_cell_, int n_q_points_,
+          int n_components_, typename Number>
+template<typename VectorType>
+inline
+void
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::distribute_local_to_global (std::vector<VectorType>  &dst,
+                              const unsigned int        first_index) const
+{
+  AssertIndexRange (first_index, dst.size());
+  Assert (n_fe_components == 1, ExcNotImplemented());
+  Assert ((n_fe_components == 1 ?
+           (first_index+n_components <= dst.size()) : true),
+          ExcIndexRange (first_index + n_components_, 0, dst.size()));
+  Assert (dof_values_initialized==true,
+          internal::ExcAccessToUninitializedField());
 
-#ifdef DEBUG
-  dof_values_initialized = true;
-#endif
+  VectorType * dst_data [n_components];
+  for (unsigned int comp=0; comp<n_components; ++comp)
+    dst_data[comp] = &dst[comp+first_index];
+
+  internal::VectorDistributorLocalToGlobal<Number> distributor;
+  read_write_operation (distributor, dst_data);
 }
 
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 template<typename VectorType>
 inline
 void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-distribute_local_to_global (VectorType &dst) const
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::distribute_local_to_global (std::vector<VectorType*>  &dst,
+                              const unsigned int         first_index) const
 {
-  AssertDimension (n_components, n_fe_components);
+  AssertIndexRange (first_index, dst.size());
+  Assert (n_fe_components == 1, ExcNotImplemented());
+  Assert ((n_fe_components == 1 ?
+           (first_index+n_components <= dst.size()) : true),
+          ExcIndexRange (first_index + n_components_, 0, dst.size()));
+  Assert (dof_values_initialized==true,
+          internal::ExcAccessToUninitializedField());
+
+  VectorType * dst_data [n_components];
+  for (unsigned int comp=0; comp<n_components; ++comp)
+    dst_data[comp] = dst[comp+first_index];
+
+  internal::VectorDistributorLocalToGlobal<Number> distributor;
+  read_write_operation (distributor, dst_data);
+}
+
+
+
+template <int dim, int dofs_per_cell_, int n_q_points_,
+          int n_components_, typename Number>
+template<typename VectorType>
+inline
+void
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::set_dof_values (VectorType &dst) const
+{
+  AssertDimension (n_components_, n_fe_components);
+  Assert (dof_values_initialized==true,
+          internal::ExcAccessToUninitializedField());
+
                                 // only need one component, but to avoid
                                 // compiler warnings, use n_components copies
                                 // here (but these will not be used)
   VectorType * dst_data [n_components];
   for (unsigned int d=0; d<n_components; ++d)
     dst_data[d] = &dst;
-  distribute_local_to_global (dst_data);
+
+  internal::VectorSetter<Number> setter;
+  read_write_operation (setter, dst_data);
 }
 
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 template<typename VectorType>
 inline
 void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-distribute_local_to_global (std::vector<VectorType>  &dst,
-                            const unsigned int        first_index) const
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::set_dof_values (std::vector<VectorType>  &dst,
+                  const unsigned int        first_index) const
 {
   AssertIndexRange (first_index, dst.size());
   Assert (n_fe_components == 1, ExcNotImplemented());
   Assert ((n_fe_components == 1 ?
            (first_index+n_components <= dst.size()) : true),
-          ExcIndexRange (first_index + n_components, 0, dst.size()));
+          ExcIndexRange (first_index + n_components_, 0, dst.size()));
+
+  Assert (dof_values_initialized==true,
+          internal::ExcAccessToUninitializedField());
 
   VectorType * dst_data [n_components];
   for (unsigned int comp=0; comp<n_components; ++comp)
     dst_data[comp] = &dst[comp+first_index];
-  distribute_local_to_global (dst_data);
+
+  internal::VectorSetter<Number> setter;
+  read_write_operation (setter, dst_data);
 }
 
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 template<typename VectorType>
 inline
 void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-distribute_local_to_global (std::vector<VectorType*>  &dst,
-                            const unsigned int         first_index) const
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::set_dof_values (std::vector<VectorType*>  &dst,
+                  const unsigned int         first_index) const
 {
   AssertIndexRange (first_index, dst.size());
   Assert (n_fe_components == 1, ExcNotImplemented());
   Assert ((n_fe_components == 1 ?
            (first_index+n_components <= dst.size()) : true),
-          ExcIndexRange (first_index + n_components, 0, dst.size()));
+          ExcIndexRange (first_index + n_components_, 0, dst.size()));
+
+  Assert (dof_values_initialized==true,
+          internal::ExcAccessToUninitializedField());
 
   VectorType * dst_data [n_components];
   for (unsigned int comp=0; comp<n_components; ++comp)
     dst_data[comp] = dst[comp+first_index];
-  distribute_local_to_global (dst_data);
+
+  internal::VectorSetter<Number> setter;
+  read_write_operation (setter, dst_data);
 }
 
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 template<typename VectorType>
 inline
 void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-distribute_local_to_global (VectorType * dst[]) const
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_dof_values_plain (const VectorType * src[])
 {
+                                // this is different from the other three
+                                // operations because we do not use
+                                // constraints here, so this is a separate
+                                // function.
   Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
-  Assert (dof_values_initialized==true,
-          internal::ExcAccessToUninitializedField());
+  Assert (dof_info.store_plain_indices == true, ExcNotInitialized());
 
                                 // loop over all local dofs. ind_local holds
                                 // local number on cell, index iterates over
                                 // the elements of index_local_to_global and
                                 // dof_indices points to the global indices
                                 // stored in index_local_to_global
-  const unsigned int * dof_indices = dof_info.begin_indices(cell);
-  const std::pair<unsigned short,unsigned short> * indicators =
-    dof_info.begin_indicators(cell);
-  const std::pair<unsigned short,unsigned short> * indicators_end =
-    dof_info.end_indicators(cell);
-  unsigned int ind_local = 0;
+  const unsigned int * dof_indices = dof_info.begin_indices_plain(cell);
 
                                 // scalar case (or case when all components
                                 // have the same degrees of freedom and sit on
@@ -2508,111 +2811,45 @@ distribute_local_to_global (VectorType * dst[]) const
       const unsigned int n_local_dofs =
         VectorizedArray<Number>::n_array_elements * dofs_per_cell;
       for (unsigned int comp=0; comp<n_components; ++comp)
-        internal::check_vector_compatibility (*dst[comp], dof_info);
-
-      const Number * local_dst_number [n_components];
+        internal::check_vector_compatibility (*src[comp], dof_info);
+      Number * local_src_number [n_components];
       for (unsigned int comp=0; comp<n_components; ++comp)
-        local_dst_number[comp] =
-          reinterpret_cast<const Number*>(values_dofs[comp]);
+        local_src_number[comp] = reinterpret_cast<Number*>(values_dofs[comp]);
+
+                                // standard case where there are sufficiently
+                                // many cells to fill all vectors
       if (at_irregular_cell == false)
         {
-                                // check whether there is no constraint at all
-          if (indicators != indicators_end)
-            {
-                                // run from one constraint to the next
-              for ( ; indicators != indicators_end; ++indicators)
-                {
-                                // distribute values up to the constraint
-                                // (values not constrained)
-                  for (unsigned int j=0; j<indicators->first; ++j)
-                    for (unsigned int comp=0; comp<n_components; ++comp)
-                      internal::vector_access (*dst[comp], dof_indices[j])
-                        += local_dst_number[comp][ind_local+j];
-                  dof_indices += indicators->first;
-                  ind_local   += indicators->first;
-
-                                // constrained case: build the local value as
-                                // a linear combination of the global value
-                                // according to constraint
-                  const Number * data_val =
-                    matrix_info.constraint_pool_begin(indicators->second);
-                  const Number * end_pool =
-                    matrix_info.constraint_pool_end(indicators->second);
-                  for ( ; data_val != end_pool; ++data_val, ++dof_indices)
-                    for (unsigned int comp=0; comp<n_components; ++comp)
-                      internal::vector_access (*dst[comp], *dof_indices)
-                        += local_dst_number[comp][ind_local] * (*data_val);
-                  ++ind_local;
-                }
-                                // distribute values after the last constraint
-                                // (values not constrained)
-              for(; ind_local<n_local_dofs; ++dof_indices, ++ind_local)
-                for (unsigned int comp=0; comp<n_components; ++comp)
-                  internal::vector_access (*dst[comp], *dof_indices)
-                    += local_dst_number[comp][ind_local];
-            }
-                                // no constraint at all: loop bounds are
-                                // known, compiler can unroll without checks
-          else
-            {
-              AssertDimension (dof_info.end_indices(cell)-dof_indices,
-                               static_cast<int>(n_local_dofs));
-              for (unsigned int j=0; j<n_local_dofs; ++j)
-                for (unsigned int comp=0; comp<n_components; ++comp)
-                  internal::vector_access (*dst[comp], dof_indices[j])
-                    += local_dst_number[comp][j];
-            }
-          return;
+          for (unsigned int j=0; j<n_local_dofs; ++j)
+            for (unsigned int comp=0; comp<n_components; ++comp)
+              local_src_number[comp][j] =
+                internal::vector_access (*src[comp], dof_indices[j]);
         }
 
-                                // irregular case
-      Assert (n_irreg_components_filled > 0, ExcInternalError());
-      for ( ; indicators != indicators_end; ++indicators)
+                                // non-standard case: need to fill in zeros
+                                // for those components that are not present
+                                // (a bit more expensive), but there is not
+                                // more than one such cell
+      else
         {
-          for(unsigned int j=0; j<indicators->first; ++j)
+          Assert (n_irreg_components_filled > 0, ExcInternalError());
+          for(unsigned int ind_local=0; ind_local<n_local_dofs;
+              ++dof_indices)
             {
                                 // non-constrained case: copy the data from
-                                // the local vector to the global one.
+                                // the global vector, src, to the local one,
+                                // local_dst.
               for (unsigned int comp=0; comp<n_components; ++comp)
-                internal::vector_access (*dst[comp], dof_indices[j])
-                  += local_dst_number[comp][ind_local];
+                local_src_number[comp][ind_local] =
+                  internal::vector_access (*src[comp], *dof_indices);
               ++ind_local;
-              if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
-                ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
-            }
-          dof_indices += indicators->first;
-
-                                // constrained case: distribute according to
-                                // the constraint
-          const Number * data_val =
-            matrix_info.constraint_pool_begin(indicators->second);
-          const Number * end_pool =
-            matrix_info.constraint_pool_end(indicators->second);
-          for ( ; data_val != end_pool; ++data_val, ++dof_indices)
-            {
-              for (unsigned int comp=0; comp<n_components; ++comp)
-                internal::vector_access (*dst[comp], *dof_indices)
-                  += local_dst_number[comp][ind_local] * (*data_val);
+              while (ind_local % VectorizedArray<Number>::n_array_elements >= n_irreg_components_filled)
+                {
+                  for (unsigned int comp=0; comp<n_components; ++comp)
+                    local_src_number[comp][ind_local] = 0.;
+                  ++ind_local;
+                }
             }
-          ++ind_local;
-          if (ind_local % VectorizedArray<Number>::n_array_elements ==
-              n_irreg_components_filled)
-            ind_local += VectorizedArray<Number>::n_array_elements-
-              n_irreg_components_filled;
-        }
-      for(; ind_local<n_local_dofs; ++dof_indices)
-        {
-          Assert (dof_indices != dof_info.end_indices(cell),
-                  ExcInternalError());
-
-                                // non-constrained case
-          for (unsigned int comp=0; comp<n_components; ++comp)
-            internal::vector_access (*dst[comp], *dof_indices)
-              += local_dst_number[comp][ind_local];
-          ++ind_local;
-          if (ind_local % VectorizedArray<Number>::n_array_elements ==
-              n_irreg_components_filled)
-            ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
         }
     }
   else
@@ -2624,408 +2861,144 @@ distribute_local_to_global (VectorType * dst[]) const
                                 // is ensured by the way MatrixFree reads
                                 // out the indices.
     {
-      internal::check_vector_compatibility (*dst[0], dof_info);
-      Assert (n_fe_components == n_components, ExcNotImplemented());
+      internal::check_vector_compatibility (*src[0], dof_info);
+      Assert (n_fe_components == n_components_, ExcNotImplemented());
       const unsigned int n_local_dofs =
         dofs_per_cell * VectorizedArray<Number>::n_array_elements * n_components;
-      const Number * local_dst_number =
-        reinterpret_cast<const Number*>(values_dofs[0]);
+      Number * local_src_number = reinterpret_cast<Number*>(values_dofs[0]);
       if (at_irregular_cell == false)
         {
-                                // check whether there is no constraint at all
-          if (indicators != indicators_end)
-            {
-                                // run from one constraint to the next
-              for ( ; indicators != indicators_end; ++indicators)
-                {
-                                // distribute values up to the constraint
-                                // (values not constrained)
-                  for (unsigned int j=0; j<indicators->first; ++j)
-                    internal::vector_access (*dst[0], dof_indices[j])
-                      += local_dst_number[ind_local+j];
-                  dof_indices += indicators->first;
-                  ind_local   += indicators->first;
-
-                                // constrained case: build the local value as
-                                // a linear combination of the global value
-                                // according to constraint
-                  const Number * data_val =
-                    matrix_info.constraint_pool_begin(indicators->second);
-                  const Number * end_pool =
-                    matrix_info.constraint_pool_end(indicators->second);
-                  for ( ; data_val != end_pool; ++data_val, ++dof_indices)
-                    internal::vector_access (*dst[0], *dof_indices)
-                      += local_dst_number[ind_local] * (*data_val);
-                  ++ind_local;
-                }
-                                // distribute values after the last constraint
-                                // (values not constrained)
-              for(; ind_local<n_local_dofs; ++dof_indices, ++ind_local)
-                internal::vector_access (*dst[0], *dof_indices)
-                  += local_dst_number[ind_local];
-            }
-                                // no constraint at all: loop bounds are
-                                // known, compiler can unroll without checks
-          else
-            {
-              AssertDimension (dof_info.end_indices(cell)-dof_indices,
-                               static_cast<int>(n_local_dofs));
-              for (unsigned int j=0; j<n_local_dofs; ++j)
-                internal::vector_access (*dst[0], dof_indices[j])
-                  += local_dst_number[j];
-            }
-          return;
+          for (unsigned int j=0; j<n_local_dofs; ++j)
+            local_src_number[j] =
+              internal::vector_access (*src[0], dof_indices[j]);
         }
 
-                                // irregular case
-      Assert (n_irreg_components_filled > 0, ExcInternalError());
-      for ( ; indicators != indicators_end; ++indicators)
+                                // non-standard case: need to fill in zeros
+                                // for those components that are not present
+                                // (a bit more expensive), but there is not
+                                // more than one such cell
+      else
         {
-          for(unsigned int j=0; j<indicators->first; ++j)
+          Assert (n_irreg_components_filled > 0, ExcInternalError());
+          for(unsigned int ind_local=0; ind_local<n_local_dofs; ++dof_indices)
             {
                                 // non-constrained case: copy the data from
-                                // the local vector to the global one.
-              internal::vector_access (*dst[0], dof_indices[j])
-                += local_dst_number[ind_local];
+                                // the global vector, src, to the local one,
+                                // local_dst.
+              local_src_number[ind_local] =
+                internal::vector_access (*src[0], *dof_indices);
               ++ind_local;
-              if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
-                ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
-            }
-          dof_indices += indicators->first;
-
-                                // constrained case: distribute according to
-                                // the constraint
-          const Number * data_val =
-            matrix_info.constraint_pool_begin(indicators->second);
-          const Number * end_pool =
-            matrix_info.constraint_pool_end(indicators->second);
-          for ( ; data_val != end_pool; ++data_val, ++dof_indices)
-            {
-              internal::vector_access (*dst[0], *dof_indices)
-                += local_dst_number[ind_local] * (*data_val);
+              while (ind_local % VectorizedArray<Number>::n_array_elements >= n_irreg_components_filled)
+                {
+                  local_src_number[ind_local] = 0.;
+                  ++ind_local;
+                }
             }
-          ++ind_local;
-          if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
-            ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
         }
-      for(; ind_local<n_local_dofs; ++dof_indices)
-        {
-          Assert (dof_indices != dof_info.end_indices(cell),
-                  ExcInternalError());
-
-                                // non-constrained case
-          internal::vector_access (*dst[0], *dof_indices)
-              += local_dst_number[ind_local];
-          ++ind_local;
-          if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
-            ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
-        }
-      Assert (dof_indices == dof_info.end_indices(cell),
-              ExcInternalError());
     }
+
+#ifdef DEBUG
+  dof_values_initialized = true;
+#endif
 }
 
 
 
+
+/*------------------------------ access to data fields ----------------------*/
+
 template <int dim, int dofs_per_cell_, int n_q_points_,
           int n_components, typename Number>
-template<typename VectorType>
 inline
-void
+const VectorizedArray<Number> *
 FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-set_dof_values (VectorType &dst) const
+begin_values () const
 {
-  AssertDimension (n_components, n_fe_components);
-                                // only need one component, but to avoid
-                                // compiler warnings, use n_components copies
-                                // here (but these will not be used)
-  VectorType * dst_data [n_components];
-  for (unsigned int d=0; d<n_components; ++d)
-    dst_data[d] = &dst;
-  set_dof_values (dst_data);
+  Assert (values_quad_initialized || values_quad_submitted,
+          ExcNotInitialized());
+  return &values_quad[0][0];
 }
 
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
           int n_components, typename Number>
-template<typename VectorType>
 inline
-void
+VectorizedArray<Number> *
 FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-set_dof_values (std::vector<VectorType>  &dst,
-                const unsigned int        first_index) const
+begin_values ()
 {
-  AssertIndexRange (first_index, dst.size());
-  Assert (n_fe_components == 1, ExcNotImplemented());
-  Assert ((n_fe_components == 1 ?
-           (first_index+n_components <= dst.size()) : true),
-          ExcIndexRange (first_index + n_components, 0, dst.size()));
-
-  VectorType * dst_data [n_components];
-  for (unsigned int comp=0; comp<n_components; ++comp)
-    dst_data[comp] = &dst[comp+first_index];
-  set_dof_values (dst_data);
+#ifdef DEBUG
+  values_quad_submitted = true;
+#endif
+  return &values_quad[0][0];
 }
 
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
           int n_components, typename Number>
-template<typename VectorType>
 inline
-void
+const VectorizedArray<Number> *
 FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-set_dof_values (std::vector<VectorType*>  &dst,
-                const unsigned int         first_index) const
+begin_gradients () const
 {
-  AssertIndexRange (first_index, dst.size());
-  Assert (n_fe_components == 1, ExcNotImplemented());
-  Assert ((n_fe_components == 1 ?
-           (first_index+n_components <= dst.size()) : true),
-          ExcIndexRange (first_index + n_components, 0, dst.size()));
-
-  VectorType * dst_data [n_components];
-  for (unsigned int comp=0; comp<n_components; ++comp)
-    dst_data[comp] = dst[comp+first_index];
-  set_dof_values (dst_data);
+  Assert (gradients_quad_initialized || gradients_quad_submitted,
+          ExcNotInitialized());
+  return &gradients_quad[0][0][0];
 }
 
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
           int n_components, typename Number>
-template<typename VectorType>
 inline
-void
+VectorizedArray<Number> *
 FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-set_dof_values (VectorType * dst[]) const
+begin_gradients ()
 {
-  Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
-  Assert (dof_values_initialized==true,
-          internal::ExcAccessToUninitializedField());
+#ifdef DEBUG
+  gradients_quad_submitted = true;
+#endif
+  return &gradients_quad[0][0][0];
+}
 
-                                // loop over all local dofs. ind_local holds
-                                // local number on cell, index iterates over
-                                // the elements of index_local_to_global and
-                                // glob_indices points to the global indices
-                                // stored in index_local_to_global
-  const unsigned int * dof_indices = dof_info.begin_indices(cell);
-  const std::pair<unsigned short,unsigned short> * indicators =
-    dof_info.begin_indicators(cell);
-  const std::pair<unsigned short,unsigned short> * indicators_end =
-    dof_info.end_indicators(cell);
-  unsigned int ind_local = 0;
 
-  if (n_fe_components == 1)
-    {
-      const unsigned int n_local_dofs =
-        VectorizedArray<Number>::n_array_elements * dofs_per_cell;
-      for (unsigned int comp=0; comp<n_components; ++comp)
-        AssertDimension (dst[comp]->size(),
-                         dof_info.vector_partitioner->size());
 
-      const Number * local_dst_number [n_components];
-      for (unsigned int comp=0; comp<n_components; ++comp)
-        local_dst_number[comp] =
-          reinterpret_cast<const Number*>(values_dofs[comp]);
-      if (at_irregular_cell == false)
-        {
-                                // check whether there is no constraint at all
-          if (indicators != indicators_end)
-            {
-                                // run from one constraint to the next
-              for ( ; indicators != indicators_end; ++indicators)
-                {
-                                // distribute values up to the constraint
-                                // (values not constrained)
-                  for (unsigned int j=0; j<indicators->first; ++j)
-                    for (unsigned int comp=0; comp<n_components; ++comp)
-                      internal::vector_access (*dst[comp], dof_indices[j])
-                        = local_dst_number[comp][ind_local+j];
-                  dof_indices += indicators->first;
-                  ind_local   += indicators->first;
-
-                                // jump over constraints
-                  const unsigned int row_length =
-                    matrix_info.constraint_pool_end(indicators->second)-
-                    matrix_info.constraint_pool_begin(indicators->second);
-                  dof_indices += row_length;
-                  ++ind_local;
-                }
-                                // distribute values after the last constraint
-                                // (values not constrained)
-              for(; ind_local<n_local_dofs; ++dof_indices, ++ind_local)
-                for (unsigned int comp=0; comp<n_components; ++comp)
-                  internal::vector_access (*dst[comp], *dof_indices)
-                    = local_dst_number[comp][ind_local];
-            }
-                                // no constraint at all: loop bounds are
-                                // known, compiler can unroll without checks
-          else
-            {
-              AssertDimension (dof_info.end_indices(cell)-dof_indices,
-                               n_local_dofs);
-              for (unsigned int j=0; j<n_local_dofs; ++j)
-                for (unsigned int comp=0; comp<n_components; ++comp)
-                  internal::vector_access (*dst[comp], dof_indices[j])
-                    = local_dst_number[comp][j];
-            }
-          return;
-        }
+template <int dim, int dofs_per_cell_, int n_q_points_,
+          int n_components, typename Number>
+inline
+const VectorizedArray<Number> *
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
+begin_hessians () const
+{
+  Assert (hessians_quad_initialized, ExcNotInitialized());
+  return &hessians_quad[0][0][0];
+}
 
-                                // irregular case
-      Assert (n_irreg_components_filled > 0, ExcInternalError());
-      for ( ; indicators != indicators_end; ++indicators)
-        {
-          for(unsigned int j=0; j<indicators->first; ++j)
-            {
-                                // non-constrained case
-              for (unsigned int comp=0; comp<n_components; ++comp)
-                internal::vector_access (*dst[comp], dof_indices[j])
-                  = local_dst_number[comp][ind_local];
-              ++ind_local;
-              if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
-                ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
-            }
-          dof_indices += indicators->first;
-
-                                // jump over constraint
-          const unsigned int row_length =
-            matrix_info.constraint_pool_end(indicators->second)-
-            matrix_info.constraint_pool_begin(indicators->second);
-          dof_indices += row_length;
-          ++ind_local;
-          if (ind_local % VectorizedArray<Number>::n_array_elements ==
-              n_irreg_components_filled)
-            ind_local += VectorizedArray<Number>::n_array_elements -
-              n_irreg_components_filled;
-        }
-      for(; ind_local<n_local_dofs; ++dof_indices)
-        {
-          Assert (dof_indices != dof_info.end_indices(cell),
-                  ExcInternalError());
-                                // non-constrained case
-          for (unsigned int comp=0; comp<n_components; ++comp)
-            internal::vector_access (*dst[comp], *dof_indices)
-              = local_dst_number[comp][ind_local];
-          ++ind_local;
-          if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
-            ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
-        }
-    }
-  else
-                                // case with vector-valued finite elements
-                                // where all components are included in one
-                                // single vector. Assumption: first come all
-                                // entries to the first component, then all
-                                // entries to the second one, and so on. This
-                                // is ensured by the way MatrixFree reads
-                                // out the indices.
-    {
-      AssertDimension (dst[0]->size(),
-                       dof_info.vector_partitioner->size());
-      Assert (n_fe_components == n_components, ExcNotImplemented());
-      const unsigned int n_local_dofs =
-        dofs_per_cell * VectorizedArray<Number>::n_array_elements * n_components;
-      const Number * local_dst_number =
-        reinterpret_cast<const Number*>(values_dofs[0]);
 
-      if (at_irregular_cell == false)
-        {
-                                // check whether there is no constraint at all
-          if (indicators != indicators_end)
-            {
-                                // run from one constraint to the next
-              for ( ; indicators != indicators_end; ++indicators)
-                {
-                                // distribute values up to the constraint
-                                // (values not constrained)
-                  for (unsigned int j=0; j<indicators->first; ++j)
-                    internal::vector_access (*dst[0], dof_indices[j])
-                      = local_dst_number[ind_local+j];
-                  dof_indices += indicators->first;
-                  ind_local   += indicators->first;
-
-                                // jump over constraints
-                  const unsigned int row_length =
-                    matrix_info.constraint_pool_end(indicators->second) -
-                    matrix_info.constraint_pool_begin(indicators->second);
-                  dof_indices += row_length;
-                  ++ind_local;
-                }
-                                // distribute values after the last constraint
-                                // (values not constrained)
-              for(; ind_local<n_local_dofs; ++dof_indices, ++ind_local)
-                internal::vector_access (*dst[0], *dof_indices)
-                  = local_dst_number[ind_local];
-            }
-                                // no constraint at all: loop bounds are
-                                // known, compiler can unroll without checks
-          else
-            {
-              AssertDimension (dof_info.end_indices(cell)-dof_indices,
-                               n_local_dofs);
-              for (unsigned int j=0; j<n_local_dofs; ++j)
-                internal::vector_access (*dst[0], dof_indices[j])
-                  = local_dst_number[j];
-            }
-          return;
-        }
 
-                                // irregular case
-      Assert (n_irreg_components_filled > 0, ExcInternalError());
-      for ( ; indicators != indicators_end; ++indicators)
-        {
-          for(unsigned int j=0; j<indicators->first; ++j)
-            {
-                                // non-constrained case
-              internal::vector_access (*dst[0], dof_indices[j])
-                = local_dst_number[ind_local];
-              ++ind_local;
-              if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
-                ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
-            }
-          dof_indices += indicators->first;
-
-                                // jump over constraint
-          const unsigned int row_length =
-            matrix_info.constraint_pool_end(indicators->second)-
-            matrix_info.constraint_pool_begin(indicators->second);
-          dof_indices += row_length;
-          ++ind_local;
-          if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
-            ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
-        }
-      for(; ind_local<n_local_dofs; ++dof_indices)
-        {
-          Assert (dof_indices != dof_info.end_indices(cell),
-                  ExcInternalError());
-                                // non-constrained case
-          internal::vector_access (*dst[0], *dof_indices)
-            = local_dst_number[ind_local];
-          ++ind_local;
-          if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
-            ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
-        }
-      Assert (dof_indices == dof_info.end_indices (cell),
-              ExcInternalError());
-    }
+template <int dim, int dofs_per_cell_, int n_q_points_,
+          int n_components, typename Number>
+inline
+VectorizedArray<Number> *
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
+begin_hessians ()
+{
+  return &hessians_quad[0][0][0];
 }
 
 
 
-// ------------------------------ access to data fields ---------------------
-
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 inline
-Tensor<1,n_components,VectorizedArray<Number> >
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_dof_value (const unsigned int dof) const
+Tensor<1,n_components_,VectorizedArray<Number> >
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_dof_value (const unsigned int dof) const
 {
   AssertIndexRange (dof, dofs_per_cell);
-  Tensor<1,n_components,VectorizedArray<Number> > return_value (false);
+  Tensor<1,n_components_,VectorizedArray<Number> > return_value (false);
   for(unsigned int comp=0;comp<n_components;comp++)
     return_value[comp] = this->values_dofs[comp][dof];
   return return_value;
@@ -3034,16 +3007,16 @@ get_dof_value (const unsigned int dof) const
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 inline
-Tensor<1,n_components,VectorizedArray<Number> >
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_value (const unsigned int q_point) const
+Tensor<1,n_components_,VectorizedArray<Number> >
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_value (const unsigned int q_point) const
 {
   Assert (this->values_quad_initialized==true,
           internal::ExcAccessToUninitializedField());
   AssertIndexRange (q_point, n_q_points);
-  Tensor<1,n_components,VectorizedArray<Number> > return_value (false);
+  Tensor<1,n_components_,VectorizedArray<Number> > return_value (false);
   for(unsigned int comp=0;comp<n_components;comp++)
     return_value[comp] = this->values_quad[comp][q_point];
   return return_value;
@@ -3052,17 +3025,17 @@ get_value (const unsigned int q_point) const
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 inline
-Tensor<1,n_components,Tensor<1,dim,VectorizedArray<Number> > >
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_gradient (const unsigned int q_point) const
+Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > >
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_gradient (const unsigned int q_point) const
 {
   Assert (this->gradients_quad_initialized==true,
           internal::ExcAccessToUninitializedField());
   AssertIndexRange (q_point, n_q_points);
 
-  Tensor<1,n_components,Tensor<1,dim,VectorizedArray<Number> > > grad_out (false);
+  Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > > grad_out (false);
 
                                 // Cartesian cell
   if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
@@ -3072,33 +3045,20 @@ get_gradient (const unsigned int q_point) const
           grad_out[comp][d] = (this->gradients_quad[comp][d][q_point] *
                                cartesian_data[0][d]);
     }
-                                // cell with general Jacobian
-  else if (this->cell_type == internal::MatrixFreeFunctions::general)
-    {
-      for(unsigned int comp=0;comp<n_components;comp++)
-        {
-          for (unsigned int d=0; d<dim; ++d)
-            {
-              grad_out[comp][d] = (jacobian[q_point][d][0] *
-                                   this->gradients_quad[comp][0][q_point]);
-              for (unsigned e=1; e<dim; ++e)
-                grad_out[comp][d] += (jacobian[q_point][d][e] *
-                                      this->gradients_quad[comp][e][q_point]);
-            }
-        }
-    }
-                                // cell with general Jacobian, but constant
-                                // within the cell
-  else // if (this->cell_type == internal::MatrixFreeFunctions::affine)
+                                // cell with general/affine Jacobian
+  else
     {
+      const Tensor<2,dim,VectorizedArray<Number> > & jac =
+        this->cell_type == internal::MatrixFreeFunctions::general ?
+        jacobian[q_point] : jacobian[0];
       for(unsigned int comp=0;comp<n_components;comp++)
         {
           for (unsigned int d=0; d<dim; ++d)
             {
-              grad_out[comp][d] = (jacobian[0][d][0] *
+              grad_out[comp][d] = (jac[d][0] *
                                    this->gradients_quad[comp][0][q_point]);
               for (unsigned e=1; e<dim; ++e)
-                grad_out[comp][d] += (jacobian[0][d][e] *
+                grad_out[comp][d] += (jac[d][e] *
                                       this->gradients_quad[comp][e][q_point]);
             }
         }
@@ -3108,12 +3068,57 @@ get_gradient (const unsigned int q_point) const
 
 
 
+namespace internal
+{
+                                // compute tmp = hess_unit(u) * J^T. do this
+                                // manually because we do not store the lower
+                                // diagonal because of symmetry
+  template <int dim, int n_q_points, typename Number>
+  inline
+  void
+  hessian_unit_times_jac (const Tensor<2,dim,VectorizedArray<Number> > &jac,
+                          const VectorizedArray<Number> hessians_quad[][n_q_points],
+                          const unsigned int            q_point,
+                          VectorizedArray<Number>       tmp[dim][dim])
+  {
+    for (unsigned int d=0; d<dim; ++d)
+      {
+        switch (dim)
+          {
+          case 1:
+            tmp[0][0] = jac[0][0] * hessians_quad[0][q_point];
+            break;
+          case 2:
+            tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] +
+                         jac[d][1] * hessians_quad[2][q_point]);
+            tmp[1][d] = (jac[d][0] * hessians_quad[2][q_point] +
+                         jac[d][1] * hessians_quad[1][q_point]);
+            break;
+          case 3:
+            tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] +
+                         jac[d][1] * hessians_quad[3][q_point] +
+                         jac[d][2] * hessians_quad[4][q_point]);
+            tmp[1][d] = (jac[d][0] * hessians_quad[3][q_point] +
+                         jac[d][1] * hessians_quad[1][q_point] +
+                         jac[d][2] * hessians_quad[5][q_point]);
+            tmp[2][d] = (jac[d][0] * hessians_quad[4][q_point] +
+                         jac[d][1] * hessians_quad[5][q_point] +
+                         jac[d][2] * hessians_quad[2][q_point]);
+            break;
+          default: Assert (false, ExcNotImplemented());
+          }
+      }
+  }
+}
+
+
+
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 inline
-Tensor<1,n_components,Tensor<2,dim,VectorizedArray<Number> > >
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_hessian (const unsigned int q_point) const
+Tensor<1,n_components_,Tensor<2,dim,VectorizedArray<Number> > >
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_hessian (const unsigned int q_point) const
 {
   Assert (this->hessians_quad_initialized==true,
           internal::ExcAccessToUninitializedField());
@@ -3167,37 +3172,9 @@ get_hessian (const unsigned int q_point) const
                                 // because it needs to access unscaled
                                 // gradient data
           VectorizedArray<Number> tmp[dim][dim];
+          internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
+                                            q_point, tmp);
 
-                                // compute tmp = hess_unit(u) * J^T. do this
-                                // manually because we do not store the lower
-                                // diagonal because of symmetry
-          for (unsigned int d=0; d<dim; ++d)
-            {
-              switch (dim)
-                {
-                case 1:
-                  tmp[0][0] = jac[0][0] * this->hessians_quad[comp][0][q_point];
-                  break;
-                case 2:
-                  tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][2][q_point]);
-                  tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][2][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][1][q_point]);
-                  break;
-                case 3:
-                  tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][3][q_point] +
-                               jac[d][2] * this->hessians_quad[comp][4][q_point]);
-                  tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][3][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][1][q_point] +
-                               jac[d][2] * this->hessians_quad[comp][5][q_point]);
-                  tmp[2][d] = (jac[d][0] * this->hessians_quad[comp][4][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][5][q_point] +
-                               jac[d][2] * this->hessians_quad[comp][2][q_point]);
-                  break;
-                default: Assert (false, ExcNotImplemented());
-                }
-            }
                                 // compute first part of hessian,
                                 // J * tmp = J * hess_unit(u) * J^T
           for (unsigned int d=0; d<dim; ++d)
@@ -3238,37 +3215,9 @@ get_hessian (const unsigned int q_point) const
                                 // because it needs to access unscaled
                                 // gradient data
           VectorizedArray<Number> tmp[dim][dim];
+          internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
+                                            q_point, tmp);
 
-                                // compute tmp = hess_unit(u) * J^T. do this
-                                // manually because we do not store the lower
-                                // diagonal because of symmetry
-          for (unsigned int d=0; d<dim; ++d)
-            {
-              switch (dim)
-                {
-                case 1:
-                  tmp[0][0] = jac[0][0] * this->hessians_quad[comp][0][q_point];
-                  break;
-                case 2:
-                  tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][2][q_point]);
-                  tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][2][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][1][q_point]);
-                  break;
-                case 3:
-                  tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][3][q_point] +
-                               jac[d][2] * this->hessians_quad[comp][4][q_point]);
-                  tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][3][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][1][q_point] +
-                               jac[d][2] * this->hessians_quad[comp][5][q_point]);
-                  tmp[2][d] = (jac[d][0] * this->hessians_quad[comp][4][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][5][q_point] +
-                               jac[d][2] * this->hessians_quad[comp][2][q_point]);
-                  break;
-                default: Assert (false, ExcNotImplemented());
-                }
-            }
                                 // compute first part of hessian,
                                 // J * tmp = J * hess_unit(u) * J^T
           for (unsigned int d=0; d<dim; ++d)
@@ -3289,23 +3238,23 @@ get_hessian (const unsigned int q_point) const
               hessian_out[comp][e][d] = hessian_out[comp][d][e];
         }
     }
-  return Tensor<1,n_components,Tensor<2,dim,VectorizedArray<Number> > >(hessian_out);
+  return Tensor<1,n_components_,Tensor<2,dim,VectorizedArray<Number> > >(hessian_out);
 }
 
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 inline
-Tensor<1,n_components,Tensor<1,dim,VectorizedArray<Number> > >
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_hessian_diagonal (const unsigned int q_point) const
+Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > >
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_hessian_diagonal (const unsigned int q_point) const
 {
   Assert (this->hessians_quad_initialized==true,
           internal::ExcAccessToUninitializedField());
   AssertIndexRange (q_point, n_q_points);
 
-  Tensor<1,n_components,Tensor<1,dim,VectorizedArray<Number> > > hessian_out (false);
+  Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > > hessian_out (false);
 
                                 // Cartesian cell
   if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
@@ -3329,37 +3278,9 @@ get_hessian_diagonal (const unsigned int q_point) const
                                 // because it needs to access unscaled
                                 // gradient data
           VectorizedArray<Number> tmp[dim][dim];
+          internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
+                                            q_point, tmp);
 
-                                // compute tmp = hess_unit(u) * J^T. do this
-                                // manually because we do not store the lower
-                                // diagonal because of symmetry
-          for (unsigned int d=0; d<dim; ++d)
-            {
-              switch (dim)
-                {
-                case 1:
-                  tmp[0][0] = jac[0][0] * this->hessians_quad[comp][0][q_point];
-                  break;
-                case 2:
-                  tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][2][q_point]);
-                  tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][2][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][1][q_point]);
-                  break;
-                case 3:
-                  tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][3][q_point] +
-                               jac[d][2] * this->hessians_quad[comp][4][q_point]);
-                  tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][3][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][1][q_point] +
-                               jac[d][2] * this->hessians_quad[comp][5][q_point]);
-                  tmp[2][d] = (jac[d][0] * this->hessians_quad[comp][4][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][5][q_point] +
-                               jac[d][2] * this->hessians_quad[comp][2][q_point]);
-                  break;
-                default: Assert (false, ExcNotImplemented());
-                }
-            }
                                 // compute only the trace part of hessian,
                                 // J * tmp = J * hess_unit(u) * J^T
           for (unsigned int d=0; d<dim; ++d)
@@ -3386,37 +3307,9 @@ get_hessian_diagonal (const unsigned int q_point) const
                                 // because it needs to access unscaled
                                 // gradient data
           VectorizedArray<Number> tmp[dim][dim];
+          internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
+                                            q_point, tmp);
 
-                                // compute tmp = hess_unit(u) * J^T. do this
-                                // manually because we do not store the lower
-                                // diagonal because of symmetry
-          for (unsigned int d=0; d<dim; ++d)
-            {
-              switch (dim)
-                {
-                case 1:
-                  tmp[0][0] = jac[0][0] * this->hessians_quad[comp][0][q_point];
-                  break;
-                case 2:
-                  tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][2][q_point]);
-                  tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][2][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][1][q_point]);
-                  break;
-                case 3:
-                  tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][3][q_point] +
-                               jac[d][2] * this->hessians_quad[comp][4][q_point]);
-                  tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][3][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][1][q_point] +
-                               jac[d][2] * this->hessians_quad[comp][5][q_point]);
-                  tmp[2][d] = (jac[d][0] * this->hessians_quad[comp][4][q_point] +
-                               jac[d][1] * this->hessians_quad[comp][5][q_point] +
-                               jac[d][2] * this->hessians_quad[comp][2][q_point]);
-                  break;
-                default: Assert (false, ExcNotImplemented());
-                }
-            }
                                 // compute only the trace part of hessian,
                                 // J * tmp = J * hess_unit(u) * J^T
           for (unsigned int d=0; d<dim; ++d)
@@ -3433,17 +3326,17 @@ get_hessian_diagonal (const unsigned int q_point) const
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 inline
-Tensor<1,n_components,VectorizedArray<Number> >
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_laplacian (const unsigned int q_point) const
+Tensor<1,n_components_,VectorizedArray<Number> >
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_laplacian (const unsigned int q_point) const
 {
   Assert (this->hessians_quad_initialized==true,
           internal::ExcAccessToUninitializedField());
   AssertIndexRange (q_point, n_q_points);
-  Tensor<1,n_components,VectorizedArray<Number> > laplacian_out (false);
-  const Tensor<1,n_components,Tensor<1,dim,VectorizedArray<Number> > > hess_diag
+  Tensor<1,n_components_,VectorizedArray<Number> > laplacian_out (false);
+  const Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > > hess_diag
     = get_hessian_diagonal(q_point);
   for (unsigned int comp=0; comp<n_components; ++comp)
     {
@@ -3457,12 +3350,12 @@ get_laplacian (const unsigned int q_point) const
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 inline
 void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-submit_dof_value (const Tensor<1,n_components,VectorizedArray<Number> > val_in,
-                  const unsigned int dof)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::submit_dof_value (const Tensor<1,n_components_,VectorizedArray<Number> > val_in,
+                    const unsigned int dof)
 {
 #ifdef DEBUG
   this->dof_values_initialized = true;
@@ -3475,12 +3368,12 @@ submit_dof_value (const Tensor<1,n_components,VectorizedArray<Number> > val_in,
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 inline
 void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-submit_value (const Tensor<1,n_components,VectorizedArray<Number> > val_in,
-              const unsigned int q_point)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::submit_value (const Tensor<1,n_components_,VectorizedArray<Number> > val_in,
+                const unsigned int q_point)
 {
 #ifdef DEBUG
   Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
@@ -3504,13 +3397,13 @@ submit_value (const Tensor<1,n_components,VectorizedArray<Number> > val_in,
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 inline
 void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-submit_gradient (const Tensor<1,n_components,
-                              Tensor<1,dim,VectorizedArray<Number> > > grad_in,
-                 const unsigned int q_point)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::submit_gradient (const Tensor<1,n_components_,
+                                Tensor<1,dim,VectorizedArray<Number> > >grad_in,
+                   const unsigned int q_point)
 {
 #ifdef DEBUG
   Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
@@ -3525,26 +3418,20 @@ submit_gradient (const Tensor<1,n_components,
           this->gradients_quad[comp][d][q_point] = (grad_in[comp][d] *
                                                     cartesian_data[0][d] * JxW);
     }
-  else if (this->cell_type == internal::MatrixFreeFunctions::general)
-    {
-      for (unsigned int comp=0; comp<n_components; ++comp)
-        for (unsigned int d=0; d<dim; ++d)
-          {
-            VectorizedArray<Number> new_val = jacobian[q_point][0][d] * grad_in[comp][0];
-            for (unsigned e=1; e<dim; ++e)
-              new_val += jacobian[q_point][e][d] * grad_in[comp][e];
-            this->gradients_quad[comp][d][q_point] = new_val * J_value[q_point];
-          }
-    }
-  else //if (this->cell_type == internal::MatrixFreeFunctions::affine)
+  else
     {
-      const VectorizedArray<Number> JxW = J_value[0] * quadrature_weights[q_point];
+      const Tensor<2,dim,VectorizedArray<Number> > &jac =
+        this->cell_type == internal::MatrixFreeFunctions::general ?
+        jacobian[q_point] : jacobian[0];
+      const VectorizedArray<Number> JxW =
+        this->cell_type == internal::MatrixFreeFunctions::general ?
+        J_value[q_point] : J_value[0] * quadrature_weights[q_point];
       for (unsigned int comp=0; comp<n_components; ++comp)
         for (unsigned int d=0; d<dim; ++d)
           {
-            VectorizedArray<Number> new_val = jacobian[0][0][d] * grad_in[comp][0];
+            VectorizedArray<Number> new_val = jac[0][d] * grad_in[comp][0];
             for (unsigned e=1; e<dim; ++e)
-              new_val += jacobian[0][e][d] * grad_in[comp][e];
+              new_val += (jac[e][d] * grad_in[comp][e]);
             this->gradients_quad[comp][d][q_point] = new_val * JxW;
           }
     }
@@ -3553,21 +3440,21 @@ submit_gradient (const Tensor<1,n_components,
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 inline
-Tensor<1,n_components,VectorizedArray<Number> >
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-integrate_value () const
+Tensor<1,n_components_,VectorizedArray<Number> >
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::integrate_value () const
 {
 #ifdef DEBUG
   Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
   Assert (this->values_quad_submitted == true,
           internal::ExcAccessToUninitializedField());
 #endif
-  Tensor<1,n_components,VectorizedArray<Number> > return_value (false);
+  Tensor<1,n_components_,VectorizedArray<Number> > return_value (false);
   for (unsigned int comp=0; comp<n_components; ++comp)
     return_value[comp] = this->values_quad[comp][0];
-  for (unsigned int q=0; q<n_q_points; ++q)
+  for (unsigned int q=1; q<n_q_points; ++q)
     for (unsigned int comp=0; comp<n_components; ++comp)
       return_value[comp] += this->values_quad[comp][q];
   return (return_value);
@@ -3575,33 +3462,33 @@ integrate_value () const
 
 
 
-/*----------------------- FEEvaluationAccess -------------------------------*/
+/*----------------------- FEEvaluationAccess --------------------------------*/
 
 
 template <int dim, int dofs_per_cell_, int n_q_points_,
-          int n_components, typename Number>
+          int n_components_, typename Number>
 inline
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
-                    const unsigned int fe_no,
-                    const unsigned int quad_no_in)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
+                      const unsigned int fe_no,
+                      const unsigned int quad_no_in)
   :
-  FEEvaluationBase <dim,dofs_per_cell_,n_q_points_,n_components,Number>
+  FEEvaluationBase <dim,dofs_per_cell_,n_q_points_,n_components_,Number>
   (data_in, fe_no, quad_no_in)
 {}
 
 
 
 
-/*-------------------- FEEvaluationAccess scalar --------------------------*/
+/*-------------------- FEEvaluationAccess scalar ----------------------------*/
 
 
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
-                    const unsigned int fe_no,
-                    const unsigned int quad_no_in)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
+                      const unsigned int fe_no,
+                      const unsigned int quad_no_in)
   :
   FEEvaluationBase <dim,dofs_per_cell_,n_q_points_,1,Number>
   (data_in, fe_no, quad_no_in)
@@ -3612,8 +3499,8 @@ FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 VectorizedArray<Number>
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-get_dof_value (const unsigned int dof) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::get_dof_value (const unsigned int dof) const
 {
   AssertIndexRange (dof, dofs_per_cell);
   return this->values_dofs[0][dof];
@@ -3624,8 +3511,8 @@ get_dof_value (const unsigned int dof) const
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 VectorizedArray<Number>
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-get_value (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::get_value (const unsigned int q_point) const
 {
   Assert (this->values_quad_initialized==true,
           internal::ExcAccessToUninitializedField());
@@ -3638,8 +3525,8 @@ get_value (const unsigned int q_point) const
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 Tensor<1,dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-get_gradient (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::get_gradient (const unsigned int q_point) const
 {
                                 // could use the base class gradient, but that
                                 // involves too many inefficient
@@ -3658,29 +3545,17 @@ get_gradient (const unsigned int q_point) const
         grad_out[d] = (this->gradients_quad[0][d][q_point] *
                        this->cartesian_data[0][d]);
     }
-                                // cell with general Jacobian
-  else if (this->cell_type == internal::MatrixFreeFunctions::general)
-    {
-      for (unsigned int d=0; d<dim; ++d)
-        {
-          grad_out[d] = (this->jacobian[q_point][d][0] *
-                         this->gradients_quad[0][0][q_point]);
-          for (unsigned e=1; e<dim; ++e)
-            grad_out[d] += (this->jacobian[q_point][d][e] *
-                            this->gradients_quad[0][e][q_point]);
-        }
-    }
-                                // cell with general Jacobian, but constant
-                                // within the cell
-  else // if (this->cell_type == internal::MatrixFreeFunctions::affine)
+                                // cell with general/constant Jacobian
+  else
     {
+      const Tensor<2,dim,VectorizedArray<Number> > &jac =
+        this->cell_type == internal::MatrixFreeFunctions::general ?
+        this->jacobian[q_point] : this->jacobian[0];
       for (unsigned int d=0; d<dim; ++d)
         {
-          grad_out[d] = (this->jacobian[0][d][0] *
-                         this->gradients_quad[0][0][q_point]);
+          grad_out[d] = (jac[d][0] * this->gradients_quad[0][0][q_point]);
           for (unsigned e=1; e<dim; ++e)
-            grad_out[d] += (this->jacobian[0][d][e] *
-                            this->gradients_quad[0][e][q_point]);
+            grad_out[d] += (jac[d][e] * this->gradients_quad[0][e][q_point]);
         }
     }
   return grad_out;
@@ -3691,8 +3566,8 @@ get_gradient (const unsigned int q_point) const
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 Tensor<2,dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-get_hessian (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::get_hessian (const unsigned int q_point) const
 {
   return BaseClass::get_hessian(q_point)[0];
 }
@@ -3702,8 +3577,8 @@ get_hessian (const unsigned int q_point) const
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 Tensor<1,dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-get_hessian_diagonal (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::get_hessian_diagonal (const unsigned int q_point) const
 {
   return BaseClass::get_hessian_diagonal(q_point)[0];
 }
@@ -3713,8 +3588,8 @@ get_hessian_diagonal (const unsigned int q_point) const
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 VectorizedArray<Number>
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-get_laplacian (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::get_laplacian (const unsigned int q_point) const
 {
   return BaseClass::get_laplacian(q_point)[0];
 }
@@ -3724,8 +3599,8 @@ get_laplacian (const unsigned int q_point) const
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 void
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-submit_dof_value (const VectorizedArray<Number> val_in,
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::submit_dof_value (const VectorizedArray<Number> val_in,
                   const unsigned int dof)
 {
 #ifdef DEBUG
@@ -3740,9 +3615,9 @@ submit_dof_value (const VectorizedArray<Number> val_in,
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 void
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-submit_value (const VectorizedArray<Number> val_in,
-              const unsigned int q_point)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::submit_value (const VectorizedArray<Number> val_in,
+                const unsigned int q_point)
 {
 #ifdef DEBUG
   Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
@@ -3766,9 +3641,9 @@ submit_value (const VectorizedArray<Number> val_in,
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 void
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-submit_gradient (const Tensor<1,dim,VectorizedArray<Number> > grad_in,
-                 const unsigned int q_point)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::submit_gradient (const Tensor<1,dim,VectorizedArray<Number> > grad_in,
+                   const unsigned int q_point)
 {
 #ifdef DEBUG
   Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
@@ -3783,24 +3658,20 @@ submit_gradient (const Tensor<1,dim,VectorizedArray<Number> > grad_in,
                                                this->cartesian_data[0][d] *
                                                JxW);
     }
-  else if (this->cell_type == internal::MatrixFreeFunctions::general)
-    {
-      for (unsigned int d=0; d<dim; ++d)
-        {
-          VectorizedArray<Number> new_val = this->jacobian[q_point][0][d] * grad_in[0];
-          for (unsigned e=1; e<dim; ++e)
-            new_val += this->jacobian[q_point][e][d] * grad_in[e];
-          this->gradients_quad[0][d][q_point] = new_val * this->J_value[q_point];
-        }
-    }
-  else //if (this->cell_type == internal::MatrixFreeFunctions::affine)
+                                // general/affine cell type
+  else
     {
-      const VectorizedArray<Number> JxW = this->J_value[0] * this->quadrature_weights[q_point];
+      const Tensor<2,dim,VectorizedArray<Number> > & jac =
+        this->cell_type == internal::MatrixFreeFunctions::general ?
+        this->jacobian[q_point] : this->jacobian[0];
+      const VectorizedArray<Number> JxW =
+        this->cell_type == internal::MatrixFreeFunctions::general ?
+        this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point];
       for (unsigned int d=0; d<dim; ++d)
         {
-          VectorizedArray<Number> new_val = this->jacobian[0][0][d] * grad_in[0];
+          VectorizedArray<Number> new_val = jac[0][d] * grad_in[0];
           for (unsigned e=1; e<dim; ++e)
-            new_val += this->jacobian[0][e][d] * grad_in[e];
+            new_val += jac[e][d] * grad_in[e];
           this->gradients_quad[0][d][q_point] = new_val * JxW;
         }
     }
@@ -3811,8 +3682,8 @@ submit_gradient (const Tensor<1,dim,VectorizedArray<Number> > grad_in,
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 VectorizedArray<Number>
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-integrate_value () const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::integrate_value () const
 {
   return BaseClass::integrate_value()[0];
 }
@@ -3820,15 +3691,15 @@ integrate_value () const
 
 
 
-/*----------------- FEEvaluationAccess vector-valued ----------------------*/
+/*----------------- FEEvaluationAccess vector-valued ------------------------*/
 
 
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
-                    const unsigned int fe_no,
-                    const unsigned int quad_no_in)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
+                      const unsigned int fe_no,
+                      const unsigned int quad_no_in)
   :
   FEEvaluationBase <dim,dofs_per_cell_,n_q_points_,dim,Number>
   (data_in, fe_no, quad_no_in)
@@ -3839,8 +3710,8 @@ FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 Tensor<2,dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-get_gradient (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::get_gradient (const unsigned int q_point) const
 {
   return BaseClass::get_gradient (q_point);
 }
@@ -3850,8 +3721,8 @@ get_gradient (const unsigned int q_point) const
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 VectorizedArray<Number>
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-get_divergence (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::get_divergence (const unsigned int q_point) const
 {
   Assert (this->gradients_quad_initialized==true,
           internal::ExcAccessToUninitializedField());
@@ -3868,32 +3739,18 @@ get_divergence (const unsigned int q_point) const
         divergence += (this->gradients_quad[d][d][q_point] *
                        this->cartesian_data[0][d]);
     }
-                                // cell with general Jacobian
-  else if (this->cell_type == internal::MatrixFreeFunctions::general)
-    {
-      divergence = (this->jacobian[q_point][0][0] *
-                    this->gradients_quad[0][0][q_point]);
-      for (unsigned e=1; e<dim; ++e)
-        divergence += (this->jacobian[q_point][0][e] *
-                       this->gradients_quad[0][e][q_point]);
-      for (unsigned int d=1; d<dim; ++d)
-        for (unsigned e=0; e<dim; ++e)
-          divergence += (this->jacobian[q_point][d][e] *
-                         this->gradients_quad[d][e][q_point]);
-    }
-                                // cell with general Jacobian, but constant
-                                // within the cell
-  else // if (this->cell_type == internal::MatrixFreeFunctions::affine)
+                                // cell with general/constant Jacobian
+  else
     {
-      divergence = (this->jacobian[0][0][0] *
-                    this->gradients_quad[0][0][q_point]);
+      const Tensor<2,dim,VectorizedArray<Number> > &jac =
+        this->cell_type == internal::MatrixFreeFunctions::general ?
+        this->jacobian[q_point] : this->jacobian[0];
+      divergence = (jac[0][0] * this->gradients_quad[0][0][q_point]);
       for (unsigned e=1; e<dim; ++e)
-        divergence += (this->jacobian[0][0][e] *
-                       this->gradients_quad[0][e][q_point]);
+        divergence += (jac[0][e] * this->gradients_quad[0][e][q_point]);
       for (unsigned int d=1; d<dim; ++d)
         for (unsigned e=0; e<dim; ++e)
-          divergence += (this->jacobian[0][d][e] *
-                         this->gradients_quad[d][e][q_point]);
+          divergence += (jac[d][e] * this->gradients_quad[d][e][q_point]);
     }
   return divergence;
 }
@@ -3903,8 +3760,8 @@ get_divergence (const unsigned int q_point) const
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 SymmetricTensor<2,dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-get_symmetric_gradient (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::get_symmetric_gradient (const unsigned int q_point) const
 {
                                 // copy from generic function into
                                 // dim-specialization function
@@ -3940,8 +3797,8 @@ get_symmetric_gradient (const unsigned int q_point) const
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 Tensor<1,dim==2?1:dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-get_curl (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::get_curl (const unsigned int q_point) const
 {
                                 // copy from generic function into
                                 // dim-specialization function
@@ -3972,8 +3829,8 @@ get_curl (const unsigned int q_point) const
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 Tensor<2,dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-get_hessian_diagonal (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::get_hessian_diagonal (const unsigned int q_point) const
 {
   Assert (this->hessians_quad_initialized==true,
           internal::ExcAccessToUninitializedField());
@@ -3987,8 +3844,8 @@ get_hessian_diagonal (const unsigned int q_point) const
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 Tensor<3,dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-get_hessian (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::get_hessian (const unsigned int q_point) const
 {
   Assert (this->hessians_quad_initialized==true,
           internal::ExcAccessToUninitializedField());
@@ -4001,9 +3858,9 @@ get_hessian (const unsigned int q_point) const
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 void
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-submit_gradient (const Tensor<2,dim,VectorizedArray<Number> > grad_in,
-                 const unsigned int q_point)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::submit_gradient (const Tensor<2,dim,VectorizedArray<Number> > grad_in,
+                   const unsigned int q_point)
 {
   BaseClass::submit_gradient (grad_in, q_point);
 }
@@ -4013,9 +3870,10 @@ submit_gradient (const Tensor<2,dim,VectorizedArray<Number> > grad_in,
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 void
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-submit_gradient (const Tensor<1,dim,Tensor<1,dim,VectorizedArray<Number> > > grad_in,
-                 const unsigned int q_point)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::submit_gradient (const Tensor<1,dim,Tensor<1,dim,VectorizedArray<Number> > >
+                                      grad_in,
+                   const unsigned int q_point)
 {
   BaseClass::submit_gradient(grad_in, q_point);
 }
@@ -4025,10 +3883,10 @@ submit_gradient (const Tensor<1,dim,Tensor<1,dim,VectorizedArray<Number> > > gra
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 void
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-submit_symmetric_gradient (const SymmetricTensor<2,dim,VectorizedArray<Number> >
- sym_grad,
-                           const unsigned int q_point)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::submit_symmetric_gradient(const SymmetricTensor<2,dim,VectorizedArray<Number> >
                                              sym_grad,
+                            const unsigned int q_point)
 {
                                 // could have used base class operator, but
                                 // that involves some overhead which is
@@ -4057,52 +3915,31 @@ submit_symmetric_gradient (const SymmetricTensor<2,dim,VectorizedArray<Number> >
                                                    this->cartesian_data[0][e]);
           }
     }
-  else if (this->cell_type == internal::MatrixFreeFunctions::general)
+                                // general/affine cell type
+  else
     {
+      const VectorizedArray<Number> JxW =
+        this->cell_type == internal::MatrixFreeFunctions::general ?
+        this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point];
+      const Tensor<2,dim,VectorizedArray<Number> > &jac =
+        this->cell_type == internal::MatrixFreeFunctions::general ?
+        this->jacobian[q_point] : this->jacobian[0];
       VectorizedArray<Number> weighted [dim][dim];
-      {
-        const VectorizedArray<Number> JxW = this->J_value[q_point];
-        for (unsigned int i=0; i<dim; ++i)
-          weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
-        for (unsigned int i=0, counter=dim; i<dim; ++i)
-          for (unsigned int j=i+1; j<dim; ++j, ++counter)
-            {
-              const VectorizedArray<Number> value = sym_grad.access_raw_entry(counter) * JxW;
-              weighted[i][j] = value;
-              weighted[j][i] = value;
-            }
-      }
-      for (unsigned int comp=0; comp<dim; ++comp)
-        for (unsigned int d=0; d<dim; ++d)
+      for (unsigned int i=0; i<dim; ++i)
+        weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
+      for (unsigned int i=0, counter=dim; i<dim; ++i)
+        for (unsigned int j=i+1; j<dim; ++j, ++counter)
           {
-            VectorizedArray<Number> new_val = this->jacobian[q_point][0][d] * weighted[comp][0];
-            for (unsigned e=1; e<dim; ++e)
-              new_val += this->jacobian[q_point][e][d] * weighted[comp][e];
-            this->gradients_quad[comp][d][q_point] = new_val;
+            const VectorizedArray<Number> value = sym_grad.access_raw_entry(counter) * JxW;
+            weighted[i][j] = value;
+            weighted[j][i] = value;
           }
-    }
-  else //if (this->cell_type == internal::MatrixFreeFunctions::affine)
-    {
-      VectorizedArray<Number> weighted [dim][dim];
-      {
-        const VectorizedArray<Number> JxW = (this->J_value[0] *
-                              this->quadrature_weights[q_point]);
-        for (unsigned int i=0; i<dim; ++i)
-          weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
-        for (unsigned int i=0, counter=dim; i<dim; ++i)
-          for (unsigned int j=i+1; j<dim; ++j, ++counter)
-            {
-              const VectorizedArray<Number> value = sym_grad.access_raw_entry(counter) * JxW;
-              weighted[i][j] = value;
-              weighted[j][i] = value;
-            }
-      }
       for (unsigned int comp=0; comp<dim; ++comp)
         for (unsigned int d=0; d<dim; ++d)
           {
-            VectorizedArray<Number> new_val = this->jacobian[q_point][0][d] * weighted[comp][0];
+            VectorizedArray<Number> new_val = jac[0][d] * weighted[comp][0];
             for (unsigned e=1; e<dim; ++e)
-              new_val += this->jacobian[q_point][e][d] * weighted[comp][e];
+              new_val += jac[e][d] * weighted[comp][e];
             this->gradients_quad[comp][d][q_point] = new_val;
           }
     }
@@ -4113,9 +3950,9 @@ submit_symmetric_gradient (const SymmetricTensor<2,dim,VectorizedArray<Number> >
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
 inline
 void
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-submit_curl (const Tensor<1,dim==2?1:dim,VectorizedArray<Number> > curl,
-             const unsigned int q_point)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::submit_curl (const Tensor<1,dim==2?1:dim,VectorizedArray<Number> > curl,
+               const unsigned int q_point)
 {
   Tensor<2,dim,VectorizedArray<Number> > grad;
   switch (dim)
@@ -4146,13 +3983,13 @@ submit_curl (const Tensor<1,dim==2?1:dim,VectorizedArray<Number> > curl,
 
 /*----------------------- FEEvaluationGeneral -------------------------------*/
 
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components,
+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 MatrixFree<dim,Number> &data_in,
-                     const unsigned int fe_no,
-                     const unsigned int quad_no_in)
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluationGeneral (const MatrixFree<dim,Number> &data_in,
+                       const unsigned int fe_no,
+                       const unsigned int quad_no_in)
   :
   BaseClass (data_in, fe_no, quad_no_in)
 {
@@ -4245,259 +4082,440 @@ FEEvaluationGeneral (const MatrixFree<dim,Number> &data_in,
 
 
 
-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>::
-evaluate (bool evaluate_val, bool evaluate_grad, bool evaluate_lapl)
+namespace internal
 {
-  Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
-  Assert (this->dof_values_initialized == true,
-          internal::ExcAccessToUninitializedField());
+                                // evaluates the given shape data in 1d-3d
+                                // using the tensor product form. does not use
+                                // the tensor product form and corresponds to
+                                // a usual matrix-matrix product
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number,
+            int direction, bool dof_to_quad, bool add>
+  inline
+  void
+  apply_tensor_product (const VectorizedArray<Number>*shape_data,
+                        const VectorizedArray<Number> in [],
+                        VectorizedArray<Number>       out [])
+  {
+    AssertIndexRange (direction, dim);
+    const int mm     = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
+              nn     = dof_to_quad ? n_q_points_1d : (fe_degree+1);
 
-  const VectorizedArray<Number> * val  = this->data.shape_values.begin();
-  const VectorizedArray<Number> * grad = this->data.shape_gradients.begin();
-  const VectorizedArray<Number> * hess = this->data.shape_hessians.begin();
+    const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
+    const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
+    const int stride    = Utilities::fixed_int_power<nn,direction>::value;
 
-  for(unsigned int comp=0;comp<n_components;comp++)
+    for (int i2=0; i2<n_blocks2; ++i2)
+      {
+        for (int i1=0; i1<n_blocks1; ++i1)
+          {
+            for (int col=0; col<nn; ++col)
+              {
+                VectorizedArray<Number> val0;
+                if (dof_to_quad == true)
+                  val0 = shape_data[col];
+                else
+                  val0 = shape_data[col*n_q_points_1d];
+                VectorizedArray<Number> res0 = val0 * in[0];
+                for (int ind=1; ind<mm; ++ind)
+                  {
+                    if (dof_to_quad == true)
+                      val0 = shape_data[ind*n_q_points_1d+col];
+                    else
+                val0 = shape_data[col*n_q_points_1d+ind];
+                    res0 += val0 * in[stride*ind];
+                  }
+                if (add == false)
+                  out[stride*col]         = res0;
+                else
+                  out[stride*col]        += res0;
+              }
+
+                                // increment: in regular case, just go to the
+                                // next point in x-direction. If we are at the
+                                // end of one chunk in x-dir, need to jump
+                                // over to the next layer in z-direction
+            switch (direction)
+              {
+              case 0:
+                in += mm;
+                out += nn;
+                break;
+              case 1:
+              case 2:
+                ++in;
+                ++out;
+                break;
+              default:
+                Assert (false, ExcNotImplemented());
+              }
+          }
+        if (direction == 1)
+          {
+            in += nn*(mm-1);
+            out += nn*(nn-1);
+          }
+      }
+  }
+
+
+
+                                // This performs the evaluation of function
+                                // values, gradients and Hessians for
+                                // tensor-product finite elements. The
+                                // operation is used for both
+                                // FEEvaluationGeneral and FEEvaluation, which
+                                // provide different functions apply_values,
+                                // apply_gradients in the individual
+                                // coordinate directions
+  template <typename FEEval>
+  inline
+  void
+  do_evaluate (FEEval    &fe_eval,
+               const bool evaluate_val,
+               const bool evaluate_grad,
+               const bool evaluate_lapl)
+  {
+    Assert (fe_eval.cell != numbers::invalid_unsigned_int,
+            ExcNotInitialized());
+    Assert (fe_eval.dof_values_initialized == true,
+            internal::ExcAccessToUninitializedField());
+
+    const unsigned int temp_size = fe_eval.dofs_per_cell > fe_eval.n_q_points ?
+      fe_eval.dofs_per_cell : fe_eval.n_q_points;
+    const unsigned int n_components = fe_eval.n_components;
+    const unsigned int dim = fe_eval.dimension;
+
+    for(unsigned int c=0; c<n_components; c++)
     {
-      VectorizedArray<Number> temp1[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
-      VectorizedArray<Number> temp2[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
+      VectorizedArray<typename FEEval::number_type> temp1[temp_size];
+      VectorizedArray<typename FEEval::number_type> temp2[temp_size];
 
-      if (dim == 3)
-        {
-          if (evaluate_grad == true)
-            {
-              // grad x
-              apply_tensor_prod<0,true,false> (grad, this->values_dofs[comp], temp1);
-              apply_tensor_prod<1,true,false> (val, temp1, temp2);
-              apply_tensor_prod<2,true,false> (val, temp2, this->gradients_quad[comp][0]);
-            }
+      switch (dim)
+      {
+      case 3:
 
-          if (evaluate_lapl == true)
-            {
-              // grad xz
-              if (evaluate_grad == false)
-                {
-                  apply_tensor_prod<0,true,false> (grad, this->values_dofs[comp], temp1);
-                  apply_tensor_prod<1,true,false> (val, temp1, temp2);
-                }
-              apply_tensor_prod<2,true,false>(grad, temp2, this->hessians_quad[comp][4]);
+        if (evaluate_grad == true)
+          {
+            // grad x
+            fe_eval.template apply_gradients<0,true,false>
+              (fe_eval.values_dofs[c], temp1);
+            fe_eval.template apply_values<1,true,false>
+              (temp1, temp2);
+            fe_eval.template apply_values<2,true,false>
+              (temp2, fe_eval.gradients_quad[c][0]);
+          }
 
-                  // grad xy
-              apply_tensor_prod<1,true,false>(grad, temp1, temp2);
-              apply_tensor_prod<2,true,false> (val, temp2, this->hessians_quad[comp][3]);
+        if (evaluate_lapl == true)
+          {
+            // grad xz
+            if (evaluate_grad == false)
+              {
+                fe_eval.template apply_gradients<0,true,false>
+                  (fe_eval.values_dofs[c], temp1);
+                fe_eval.template apply_values<1,true,false>
+                  (temp1, temp2);
+              }
+            fe_eval.template apply_gradients<2,true,false>
+              (temp2, fe_eval.hessians_quad[c][4]);
+
+            // grad xy
+            fe_eval.template apply_gradients<1,true,false>
+              (temp1, temp2);
+            fe_eval.template apply_values<2,true,false>
+              (temp2, fe_eval.hessians_quad[c][3]);
+
+            // grad xx
+            fe_eval.template apply_hessians<0,true,false>
+              (fe_eval.values_dofs[c], temp1);
+            fe_eval.template apply_values<1,true,false>
+              (temp1, temp2);
+            fe_eval.template apply_values<2,true,false>
+              (temp2, fe_eval.hessians_quad[c][0]);
+          }
 
-              // grad xx
-              apply_tensor_prod<0,true,false> (hess, this->values_dofs[comp], temp1);
-              apply_tensor_prod<1,true,false> (val, temp1, temp2);
-              apply_tensor_prod<2,true,false> (val, temp2, this->hessians_quad[comp][0]);
-            }
+        // grad y
+        fe_eval.template apply_values<0,true,false>
+          (fe_eval.values_dofs[c], temp1);
+        if (evaluate_grad == true)
+          {
+            fe_eval.template apply_gradients<1,true,false>
+              (temp1, temp2);
+            fe_eval.template apply_values<2,true,false>
+              (temp2, fe_eval.gradients_quad[c][1]);
+          }
 
-          // grad y
-          apply_tensor_prod<0,true,false> (val, this->values_dofs[comp], temp1);
-          if (evaluate_grad == true)
-            {
-              apply_tensor_prod<1,true,false> (grad, temp1, temp2);
-              apply_tensor_prod<2,true,false> (val, temp2, this->gradients_quad[comp][1]);
-            }
+        if (evaluate_lapl == true)
+          {
+            // grad yz
+            if (evaluate_grad == false)
+              fe_eval.template apply_gradients<1,true,false>
+                (temp1, temp2);
+            fe_eval.template apply_gradients<2,true,false>
+              (temp2, fe_eval.hessians_quad[c][5]);
+
+            // grad yy
+            fe_eval.template apply_hessians<1,true,false>
+              (temp1, temp2);
+            fe_eval.template apply_values<2,true,false>
+              (temp2, fe_eval.hessians_quad[c][1]);
+          }
 
-          if (evaluate_lapl == true)
-            {
-              // grad yz
-              if (evaluate_grad == false)
-                apply_tensor_prod<1,true,false> (grad, temp1, temp2);
-              apply_tensor_prod<2,true,false> (grad, temp2, this->hessians_quad[comp][5]);
+        // grad z: can use the values applied in x direction stored in temp1
+        fe_eval.template apply_values<1,true,false>
+          (temp1, temp2);
+        if (evaluate_grad == true)
+          fe_eval.template apply_gradients<2,true,false>
+            (temp2, fe_eval.gradients_quad[c][2]);
 
-              // grad yy
-              apply_tensor_prod<1,true,false> (hess, temp1, temp2);
-              apply_tensor_prod<2,true,false> (val, temp2, this->hessians_quad[comp][1]);
-            }
+        // grad zz: can use the values applied in x and y direction stored
+        // in temp2
+        if (evaluate_lapl == true)
+          fe_eval.template apply_hessians<2,true,false>
+            (temp2, fe_eval.hessians_quad[c][2]);
 
-          // grad z: can use the values applied in x direction stored in temp1
-          apply_tensor_prod<1,true,false> (val, temp1, temp2);
-          if (evaluate_grad == true)
-            apply_tensor_prod<2,true,false> (grad, temp2, this->gradients_quad[comp][2]);
+        // val: can use the values applied in x & y direction stored in temp2
+        if (evaluate_val == true)
+          fe_eval.template apply_values<2,true,false>
+            (temp2, fe_eval.values_quad[c]);
 
-          // grad zz: can use the values applied in x and y direction stored in temp2
-          if (evaluate_lapl == true)
-            apply_tensor_prod<2,true,false> (hess, temp2, this->hessians_quad[comp][2]);
+        break;
 
-          // val: can use the values applied in x & y direction stored in temp2
-          if (evaluate_val == true)
-            apply_tensor_prod<2,true,false> (val, temp2, this->values_quad[comp]);
-        }
-      else if (dim == 2)
-        {
-          // grad x
-          if (evaluate_grad == true)
-            {
-              apply_tensor_prod<0,true,false> (grad, this->values_dofs[comp], temp1);
-              apply_tensor_prod<1,true,false> (val, temp1, this->gradients_quad[comp][0]);
-            }
-          if (evaluate_lapl == true)
-            {
-              // grad xy
-              if (evaluate_grad == false)
-                apply_tensor_prod<0,true,false> (grad, this->values_dofs[comp], temp1);
-              apply_tensor_prod<1,true,false> (grad, temp1, this->hessians_quad[comp][2]);
+      case 2:
+
+        // grad x
+        if (evaluate_grad == true)
+          {
+            fe_eval.template apply_gradients<0,true,false>
+              (fe_eval.values_dofs[c], temp1);
+            fe_eval.template apply_values<1,true,false>
+              (temp1, fe_eval.gradients_quad[c][0]);
+          }
+        if (evaluate_lapl == true)
+          {
+            // grad xy
+            if (evaluate_grad == false)
+              fe_eval.template apply_gradients<0,true,false>
+                (fe_eval.values_dofs[c], temp1);
+            fe_eval.template apply_gradients<1,true,false>
+              (temp1, fe_eval.hessians_quad[c][2]);
+
+            // grad xx
+            fe_eval.template apply_hessians<0,true,false>
+              (fe_eval.values_dofs[c], temp1);
+            fe_eval.template apply_values<1,true,false>
+              (temp1, fe_eval.hessians_quad[c][0]);
+          }
+
+        // grad y
+        fe_eval.template apply_values<0,true,false>
+          (fe_eval.values_dofs[c], temp1);
+        if (evaluate_grad == true)
+          fe_eval.template apply_gradients<1,true,false>
+            (temp1, fe_eval.gradients_quad[c][1]);
+
+        // grad yy
+        if (evaluate_lapl == true)
+          fe_eval.template apply_hessians<1,true,false>
+            (temp1, fe_eval.hessians_quad[c][1]);
+
+        // val: can use values applied in x
+        if (evaluate_val == true)
+          fe_eval.template apply_values<1,true,false>
+            (temp1, fe_eval.values_quad[c]);
+
+        break;
+
+      case 1:
+        if (evaluate_val == true)
+          fe_eval.template apply_values<0,true,false>
+            (fe_eval.values_dofs[c], fe_eval.values_quad[c]);
+        if (evaluate_grad == true)
+          fe_eval.template apply_gradients<0,true,false>
+            (fe_eval.values_dofs[c], fe_eval.gradients_quad[c][0]);
+        if (evaluate_lapl == true)
+          fe_eval.template apply_hessians<0,true,false>
+            (fe_eval.values_dofs[c], fe_eval.hessians_quad[c][0]);
+        break;
+
+      default:
+        Assert (false, ExcNotImplemented());
+      }
+    }
+
+#ifdef DEBUG
+    if (evaluate_val == true)
+      fe_eval.values_quad_initialized = true;
+    if (evaluate_grad == true)
+      fe_eval.gradients_quad_initialized = true;
+    if (evaluate_lapl == true)
+      fe_eval.hessians_quad_initialized  = true;
+#endif
+  }
+
+
+
+  template <typename FEEval>
+  inline
+  void
+  do_integrate (FEEval    &fe_eval,
+                const bool integrate_val,
+                const bool integrate_grad)
+  {
+    Assert (fe_eval.cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+    if (integrate_val == true)
+      Assert (fe_eval.values_quad_submitted == true,
+              ExcAccessToUninitializedField());
+    if (integrate_grad == true)
+      Assert (fe_eval.gradients_quad_submitted == true,
+              ExcAccessToUninitializedField());
+
+    const unsigned int temp_size = fe_eval.dofs_per_cell > fe_eval.n_q_points ?
+      fe_eval.dofs_per_cell : fe_eval.n_q_points;
+    const unsigned int n_components = fe_eval.n_components;
+    const unsigned int dim = fe_eval.dimension;
+
+
+    for(unsigned int c=0; c<n_components; c++)
+    {
+      VectorizedArray<typename FEEval::number_type> temp1[temp_size];
+      VectorizedArray<typename FEEval::number_type> temp2[temp_size];
+
+      switch (dim)
+      {
+      case 3:
+
+        if (integrate_val == true)
+          {
+            // val
+            fe_eval.template apply_values<0,false,false>
+              (fe_eval.values_quad[c], temp1);
+          }
+        if (integrate_grad == true)
+          {
+            // grad x: can sum to temporary value in temp1
+            if (integrate_val == true)
+              fe_eval.template apply_gradients<0,false,true>
+                (fe_eval.gradients_quad[c][0], temp1);
+            else
+              fe_eval.template apply_gradients<0,false,false>
+                (fe_eval.gradients_quad[c][0], temp1);
+          }
+        fe_eval.template apply_values<1,false,false>
+          (temp1, temp2);
+        if (integrate_grad == true)
+          {
+            // grad y: can sum to temporary x value in temp2
+            fe_eval.template apply_values<0,false,false>
+              (fe_eval.gradients_quad[c][1], temp1);
+            fe_eval.template apply_gradients<1,false,true>
+              (temp1, temp2);
+          }
+        fe_eval.template apply_values<2,false,false>
+          (temp2, fe_eval.values_dofs[c]);
+        if (integrate_grad == true)
+          {
+            // grad z: can sum to temporary x and y value in output
+            fe_eval.template apply_values<0,false,false>
+              (fe_eval.gradients_quad[c][2], temp1);
+            fe_eval.template apply_values<1,false,false>
+              (temp1, temp2);
+            fe_eval.template apply_gradients<2,false,true>
+              (temp2, fe_eval.values_dofs[c]);
+          }
+
+        break;
+
+      case 2:
+
+        // val
+        if (integrate_val == true)
+          fe_eval.template apply_values<0,false,false>
+            (fe_eval.values_quad[c], temp1);
+        if (integrate_grad == true)
+          {
+            //grad x
+            if (integrate_val == true)
+              fe_eval.template apply_gradients<0,false,true>
+                (fe_eval.gradients_quad[c][0], temp1);
+            else
+              fe_eval.template apply_gradients<0,false,false>
+                (fe_eval.gradients_quad[c][0], temp1);
+          }
+        fe_eval.template apply_values<1,false,false>
+          (temp1, fe_eval.values_dofs[c]);
+        if (integrate_grad == true)
+          {
+            // grad y
+            fe_eval.template apply_values<0,false,false>
+              (fe_eval.gradients_quad[c][1], temp1);
+            fe_eval.template apply_gradients<1,false,true>
+              (temp1, fe_eval.values_dofs[c]);
+          }
 
-              // grad xx
-              apply_tensor_prod<0,true,false> (hess, this->values_dofs[comp], temp1);
-              apply_tensor_prod<1,true,false> (val, temp1, this->hessians_quad[comp][0]);
-            }
+        break;
 
-          // grad y
-          apply_tensor_prod<0,true,false> (val, this->values_dofs[comp], temp1);
-          if (evaluate_grad == true)
-            apply_tensor_prod<1,true,false> (grad, temp1, this->gradients_quad[comp][1]);
+      case 1:
 
-          // grad yy
-          if (evaluate_lapl == true)
-            apply_tensor_prod<1,true,false> (hess, temp1, this->hessians_quad[comp][1]);
+        if (integrate_grad == true)
+          fe_eval.template apply_gradients<0,false,false>
+            (fe_eval.gradients_quad[c][0], fe_eval.values_dofs[c]);
+        if (integrate_val == true)
+          {
+            if (integrate_grad == true)
+              fe_eval.template apply_values<0,false,true>
+                (fe_eval.values_quad[c], fe_eval.values_dofs[c]);
+            else
+              fe_eval.template apply_values<0,false,false>
+                (fe_eval.values_quad[c], fe_eval.values_dofs[c]);
+          }
+        break;
 
-          // val: can use values applied in x
-          if (evaluate_val == true)
-            apply_tensor_prod<1,true,false> (val, temp1, this->values_quad[comp]);
-        }
-      else if (dim == 1)
-        {
-          if (evaluate_val == true)
-            apply_tensor_prod<0,true,false> (val, this->values_dofs[comp],
-                                             this->values_quad[comp]);
-          if (evaluate_grad == true)
-            apply_tensor_prod<0,true,false> (grad, this->values_dofs[comp],
-                                             this->gradients_quad[comp][0]);
-          if (evaluate_lapl == true)
-            apply_tensor_prod<0,true,false> (hess, this->values_dofs[comp],
-                                             this->hessians_quad[comp][0]);
-        }
+      default:
+        Assert (false, ExcNotImplemented());
+      }
     }
 
 #ifdef DEBUG
-  if (evaluate_val == true)
-    this->values_quad_initialized = true;
-  if (evaluate_grad == true)
-    this->gradients_quad_initialized = true;
-  if (evaluate_lapl == true)
-    this->hessians_quad_initialized  = true;
+  fe_eval.dof_values_initialized = true;
 #endif
+  }
 }
 
 
 
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components,
+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>::
-integrate (bool integrate_val,bool integrate_grad)
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::evaluate (const bool evaluate_val,
+            const bool evaluate_grad,
+            const bool evaluate_lapl)
 {
-#ifdef DEBUG
-  Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
-  if (integrate_val == true)
-    Assert (this->values_quad_submitted == true,
-            internal::ExcAccessToUninitializedField());
-  if (integrate_grad == true)
-    Assert (this->gradients_quad_submitted == true,
-            internal::ExcAccessToUninitializedField());
-#endif
-
-  const VectorizedArray<Number> * val  = this->data.shape_values.begin();
-  const VectorizedArray<Number> * grad = this->data.shape_gradients.begin();
+  internal::do_evaluate (*this, evaluate_val, evaluate_grad, evaluate_lapl);
+}
 
-  for(unsigned int comp=0;comp<n_components;comp++)
-    {
-      VectorizedArray<Number> temp1[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
-      VectorizedArray<Number> temp2[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
 
-      if (dim == 3)
-        {
-          if (integrate_val == true)
-            {
-              // val
-              apply_tensor_prod<0,false,false> (val, this->values_quad[comp], temp1);
-            }
-          if (integrate_grad == true)
-            {
-              // grad x: can sum to temporary value in temp1
-              if (integrate_val == true)
-                apply_tensor_prod<0,false,true>
-                  (grad, this->gradients_quad[comp][0],temp1);
-              else
-                apply_tensor_prod<0,false,false>
-                  (grad, this->gradients_quad[comp][0],temp1);
-            }
-          apply_tensor_prod<1,false,false> (val, temp1, temp2);
-          if (integrate_grad == true)
-            {
-              // grad y: can sum to temporary x value in temp2
-              apply_tensor_prod<0,false,false> (val, this->gradients_quad[comp][1], temp1);
-              apply_tensor_prod<1,false,true> (grad, temp1, temp2);
-            }
-          apply_tensor_prod<2,false,false> (val, temp2, this->values_dofs[comp]);
-          if (integrate_grad == true)
-            {
-              // grad z: can sum to temporary x and y value in output
-              apply_tensor_prod<0,false,false> (val, this->gradients_quad[comp][2], temp1);
-              apply_tensor_prod<1,false,false> (val, temp1, temp2);
-              apply_tensor_prod<2,false,true> (grad, temp2, this->values_dofs[comp]);
-            }
-        }
-      else if (dim == 2)
-        {
-          // val
-          if (integrate_val == true)
-            apply_tensor_prod<0,false,false> (val, this->values_quad[comp], temp1);
-          if (integrate_grad == true)
-            {
-              //grad x
-              if (integrate_val == true)
-                apply_tensor_prod<0,false,true>
-                  (grad, this->gradients_quad[comp][0],temp1);
-              else
-                apply_tensor_prod<0,false,false>
-                  (grad, this->gradients_quad[comp][0],temp1);
-            }
-          apply_tensor_prod<1,false,false> (val, temp1, this->values_dofs[comp]);
-          if (integrate_grad == true)
-            {
-              // grad y
-              apply_tensor_prod<0,false,false> (grad, this->gradients_quad[comp][1], temp1);
-              apply_tensor_prod<1,false,true> (val, temp1, this->values_dofs[comp]);
-            }
-        }
-      else if (dim == 1)
-        {
-          if (integrate_grad == true)
-            apply_tensor_prod<0,false,false> (grad, this->gradients_quad[comp][0],
-                                              this->values_dofs[comp]);
-          if (integrate_val == true)
-            {
-              if (integrate_grad == true)
-                apply_tensor_prod<0,false,true> (val, this->values_quad[comp],
-                                                 this->values_dofs[comp]);
-              else
-                apply_tensor_prod<0,false,false> (val, this->values_quad[comp],
-                                                  this->values_dofs[comp]);
-            }
-        }
-    }
 
-#ifdef DEBUG
-  this->dof_values_initialized = true;
-#endif
+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>
+::integrate (const bool integrate_val,
+             const bool integrate_grad)
+{
+  internal::do_integrate (*this, integrate_val, integrate_grad);
 }
 
 
 
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components,
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
           typename Number>
 inline
 Point<dim,VectorizedArray<Number> >
-FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components,Number>::
-quadrature_point (const unsigned int q) const
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::quadrature_point (const unsigned int q) const
 {
   Assert (this->mapping_info.quadrature_points_initialized == true,
           ExcNotInitialized());
@@ -4535,95 +4553,64 @@ quadrature_point (const unsigned int q) const
 }
 
 
-                                // General tensor product application for up
-                                // to three spatial dimensions. Does not
-                                // assume any symmetry in the shape values
-                                // field
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components,
+
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
           typename Number>
 template <int direction, bool dof_to_quad, bool add>
 inline
 void
-FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components,Number>::
-apply_tensor_prod (const VectorizedArray<Number>*shape_data,
-                   const VectorizedArray<Number> input [],
-                   VectorizedArray<Number>       output [])
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::apply_values(const VectorizedArray<Number> in [],
+               VectorizedArray<Number>       out [])
 {
-  AssertIndexRange (direction, dim);
-  const int mm     = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
-            nn     = dof_to_quad ? n_q_points_1d : (fe_degree+1);
+  internal::apply_tensor_product<dim,fe_degree,n_q_points_1d,Number,
+                                 direction, dof_to_quad, add>
+    (this->data.shape_values.begin(), in, out);
+}
 
-  const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
-  const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
-  const int stride    = ((direction > 0 ? nn : 1 ) *
-                         (direction > 1 ? nn : 1));
 
-  const VectorizedArray<Number> * in = &input[0];
-  VectorizedArray<Number> * out = &output[0];
-  for (int i2=0; i2<n_blocks2; ++i2)
-  {
-    for (int i1=0; i1<n_blocks1; ++i1)
-    {
-      for (int col=0; col<nn; ++col)
-        {
-          VectorizedArray<Number> val0;
-          if (dof_to_quad == true)
-            val0 = shape_data[col];
-          else
-            val0 = shape_data[col*n_q_points_1d];
-          VectorizedArray<Number> res0 = val0 * in[0];
-          for (int ind=1; ind<mm; ++ind)
-            {
-              if (dof_to_quad == true)
-                val0 = shape_data[ind*n_q_points_1d+col];
-              else
-                val0 = shape_data[col*n_q_points_1d+ind];
-              res0 += val0 * in[stride*ind];
-            }
-          if (add == false)
-            out[stride*col]         = res0;
-          else
-            out[stride*col]        += res0;
-        }
 
-                                // increment: in regular case, just go to the
-                                // next point in x-direction. If we are at the
-                                // end of one chunk in x-dir, need to jump
-                                // over to the next layer in z-direction
-      switch (direction)
-        {
-        case 0:
-          in += mm;
-          out += nn;
-          break;
-        case 1:
-        case 2:
-          ++in;
-          ++out;
-          break;
-        default:
-          Assert (false, ExcNotImplemented());
-        }
-    }
-    if (direction == 1)
-      {
-        in += nn*(mm-1);
-        out += nn*(nn-1);
-      }
-  }
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
+          typename Number>
+template <int direction, bool dof_to_quad, bool add>
+inline
+void
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::apply_gradients(const VectorizedArray<Number> in [],
+                  VectorizedArray<Number>       out [])
+{
+  internal::apply_tensor_product<dim,fe_degree,n_q_points_1d,Number,
+                                 direction, dof_to_quad, add>
+    (this->data.shape_gradients.begin(), in, out);
+}
+
+
+
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
+          typename Number>
+template <int direction, bool dof_to_quad, bool add>
+inline
+void
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::apply_hessians(const VectorizedArray<Number> in [],
+                 VectorizedArray<Number>       out [])
+{
+  internal::apply_tensor_product<dim,fe_degree,n_q_points_1d,Number,
+                                 direction, dof_to_quad, add>
+    (this->data.shape_hessians.begin(), in, out);
 }
 
 
-/*----------------------- FEEvaluation -------------------------------*/
+/*-------------------------- FEEvaluation -----------------------------------*/
 
 
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components,
+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 MatrixFree<dim,Number> &data_in,
-              const unsigned int fe_no,
-              const unsigned int quad_no)
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluation (const MatrixFree<dim,Number> &data_in,
+                const unsigned int fe_no,
+                const unsigned int quad_no)
   :
   BaseClass (data_in, fe_no, quad_no)
 {
@@ -4685,257 +4672,42 @@ FEEvaluation (const MatrixFree<dim,Number> &data_in,
 
 
 
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components,
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
           typename Number>
 inline
 void
-FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,Number>::
-evaluate (bool evaluate_val, bool evaluate_grad, bool evaluate_lapl)
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::evaluate (const bool evaluate_val,
+            const bool evaluate_grad,
+            const bool evaluate_lapl)
 {
-  Assert (this->cell != numbers::invalid_unsigned_int,
-          ExcNotInitialized());
-  Assert (this->dof_values_initialized == true,
-          internal::ExcAccessToUninitializedField());
-
-  for(unsigned int comp=0;comp<n_components;comp++)
-    {
-      VectorizedArray<Number> temp1[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
-      VectorizedArray<Number> temp2[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
-
-      if (dim == 3)
-        {
-          if (evaluate_grad == true)
-            {
-              // grad x
-              apply_gradients<0,true,false> (this->values_dofs[comp], temp1);
-              apply_values<1,true,false> (temp1, temp2);
-              apply_values<2,true,false> (temp2, this->gradients_quad[comp][0]);
-            }
-
-          if (evaluate_lapl == true)
-            {
-              // grad xz
-              if (evaluate_grad == false)
-                {
-                  apply_gradients<0,true,false> (this->values_dofs[comp], temp1);
-                  apply_values<1,true,false> (temp1, temp2);
-                }
-              apply_gradients<2,true,false>(temp2, this->hessians_quad[comp][4]);
-
-                  // grad xy
-              apply_gradients<1,true,false>(temp1, temp2);
-              apply_values<2,true,false> (temp2, this->hessians_quad[comp][3]);
-
-              // grad xx
-              apply_hessians<0,true,false> (this->values_dofs[comp], temp1);
-              apply_values<1,true,false> (temp1, temp2);
-              apply_values<2,true,false> (temp2, this->hessians_quad[comp][0]);
-            }
-
-          // grad y
-          apply_values<0,true,false> (this->values_dofs[comp], temp1);
-          if (evaluate_grad == true)
-            {
-              apply_gradients<1,true,false> (temp1, temp2);
-              apply_values<2,true,false> (temp2, this->gradients_quad[comp][1]);
-            }
-
-          if (evaluate_lapl == true)
-            {
-              // grad yz
-              if (evaluate_grad == false)
-                apply_gradients<1,true,false> (temp1, temp2);
-              apply_gradients<2,true,false> (temp2, this->hessians_quad[comp][5]);
-
-              // grad yy
-              apply_hessians<1,true,false> (temp1, temp2);
-              apply_values<2,true,false> (temp2, this->hessians_quad[comp][1]);
-            }
-
-          // grad z: can use the values applied in x direction stored in temp1
-          apply_values<1,true,false> (temp1, temp2);
-          if (evaluate_grad == true)
-            apply_gradients<2,true,false> (temp2, this->gradients_quad[comp][2]);
-
-          // grad zz: can use the values applied in x and y direction stored in temp2
-          if (evaluate_lapl == true)
-            apply_hessians<2,true,false> (temp2, this->hessians_quad[comp][2]);
-
-          // val: can use the values applied in x & y direction stored in temp2
-          if (evaluate_val == true)
-            apply_values<2,true,false> (temp2, this->values_quad[comp]);
-        }
-      else if (dim == 2)
-        {
-          // grad x
-          if (evaluate_grad == true)
-            {
-              apply_gradients<0,true,false> (this->values_dofs[comp], temp1);
-              apply_values<1,true,false> (temp1, this->gradients_quad[comp][0]);
-            }
-          if (evaluate_lapl == true)
-            {
-              // grad xy
-              if (evaluate_grad == false)
-                apply_gradients<0,true,false> (this->values_dofs[comp], temp1);
-              apply_gradients<1,true,false> (temp1, this->hessians_quad[comp][2]);
-
-              // grad xx
-              apply_hessians<0,true,false> (this->values_dofs[comp], temp1);
-              apply_values<1,true,false> (temp1, this->hessians_quad[comp][0]);
-            }
-
-          // grad y
-          apply_values<0,true,false> (this->values_dofs[comp], temp1);
-          if (evaluate_grad == true)
-            apply_gradients<1,true,false> (temp1, this->gradients_quad[comp][1]);
-
-          // grad yy
-          if (evaluate_lapl == true)
-            apply_hessians<1,true,false> (temp1, this->hessians_quad[comp][1]);
-
-          // val: can use values applied in x
-          if (evaluate_val == true)
-            apply_values<1,true,false> (temp1, this->values_quad[comp]);
-        }
-      else if (dim == 1)
-        {
-          if (evaluate_val == true)
-            apply_values<0,true,false> (this->values_dofs[comp],
-                                        this->values_quad[comp]);
-          if (evaluate_grad == true)
-            apply_gradients<0,true,false> (this->values_dofs[comp],
-                                           this->gradients_quad[comp][0]);
-          if (evaluate_lapl == true)
-            apply_hessians<0,true,false> (this->values_dofs[comp],
-                                            this->hessians_quad[comp][0]);
-        }
-    }
-
-#ifdef DEBUG
-  if (evaluate_val == true)
-    this->values_quad_initialized = true;
-  if (evaluate_grad == true)
-    this->gradients_quad_initialized = true;
-  if (evaluate_lapl == true)
-    this->hessians_quad_initialized  = true;
-#endif
+  internal::do_evaluate (*this, evaluate_val, evaluate_grad, evaluate_lapl);
 }
 
 
 
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components,
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
           typename Number>
 inline
 void
-FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,Number>::
-integrate (bool integrate_val,bool integrate_grad)
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::integrate (bool integrate_val,bool integrate_grad)
 {
-#ifdef DEBUG
-  Assert (this->cell != numbers::invalid_unsigned_int,
-          ExcNotInitialized());
-  if (integrate_val == true)
-    Assert (this->values_quad_submitted == true,
-            internal::ExcAccessToUninitializedField());
-  if (integrate_grad == true)
-    Assert (this->gradients_quad_submitted == true,
-            internal::ExcAccessToUninitializedField());
-#endif
-
-  for(unsigned int comp=0;comp<n_components;comp++)
-    {
-      VectorizedArray<Number> temp1[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
-      VectorizedArray<Number> temp2[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
-
-      if (dim == 3)
-        {
-          if (integrate_val == true)
-            {
-              // val
-              apply_values<0,false,false> (this->values_quad[comp], temp1);
-            }
-          if (integrate_grad == true)
-            {
-              // grad x: can sum to temporary value in temp1
-              if (integrate_val == true)
-                apply_gradients<0,false,true> (this->gradients_quad[comp][0],
-                                               temp1);
-              else
-                apply_gradients<0,false,false> (this->gradients_quad[comp][0],
-                                                temp1);
-            }
-          apply_values<1,false,false> (temp1, temp2);
-          if (integrate_grad == true)
-            {
-              // grad y: can sum to temporary x value in temp2
-              apply_values<0,false,false> (this->gradients_quad[comp][1], temp1);
-              apply_gradients<1,false,true> (temp1, temp2);
-            }
-          apply_values<2,false,false> (temp2, this->values_dofs[comp]);
-          if (integrate_grad == true)
-            {
-              // grad z: can sum to temporary x and y value in output
-              apply_values<0,false,false> (this->gradients_quad[comp][2], temp1);
-              apply_values<1,false,false> (temp1, temp2);
-              apply_gradients<2,false,true> (temp2, this->values_dofs[comp]);
-            }
-        }
-      else if (dim == 2)
-        {
-          // val
-          if (integrate_val == true)
-            apply_values<0,false,false> (this->values_quad[comp], temp1);
-          if (integrate_grad == true)
-            {
-              //grad x
-              if (integrate_val == true)
-                apply_gradients<0,false,true> (this->gradients_quad[comp][0],
-                                               temp1);
-              else
-                apply_gradients<0,false,false> (this->gradients_quad[comp][0],
-                                               temp1);
-            }
-          apply_values<1,false,false> (temp1, this->values_dofs[comp]);
-          if (integrate_grad == true)
-            {
-              // grad y
-              apply_values<0,false,false> (this->gradients_quad[comp][1], temp1);
-              apply_gradients<1,false,true> (temp1, this->values_dofs[comp]);
-            }
-        }
-      else if (dim == 1)
-        {
-          if (integrate_grad == true)
-            apply_gradients<0,false,false> (this->gradients_quad[comp][0],
-                                            this->values_dofs[comp]);
-          if (integrate_val == true)
-            {
-              if (integrate_grad == true)
-                apply_values<0,false,true> (this->values_quad[comp],
-                                            this->values_dofs[comp]);
-              else
-                apply_values<0,false,false> (this->values_quad[comp],
-                                             this->values_dofs[comp]);
-            }
-        }
-    }
-#ifdef DEBUG
-  this->dof_values_initialized = true;
-#endif
+  internal::do_integrate (*this, integrate_val, integrate_grad);
 }
 
 
 
-// ----------------- optimized implementation tensor product symmetric case
+/*----------------- optimized implementation tensor product symmetric case --*/
 
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components,
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
           typename Number>
 template <int direction, bool dof_to_quad, bool add>
 inline
 void
-FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,Number>::
-apply_values (const VectorizedArray<Number> input [],
-              VectorizedArray<Number>       output [])
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::apply_values (const VectorizedArray<Number> in [],
+                VectorizedArray<Number>       out [])
 {
   AssertIndexRange (direction, dim);
   const int mm     = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
@@ -4945,50 +4717,78 @@ apply_values (const VectorizedArray<Number> input [],
 
   const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
   const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
-  const int stride    = ((direction > 0 ? nn : 1 ) *
-                         (direction > 1 ? nn : 1));
-
-  const VectorizedArray<Number> * in = &input[0];
-  VectorizedArray<Number> * out = &output[0];
+  const int stride    = Utilities::fixed_int_power<nn,direction>::value;
+
+                                // This loop specializes the general
+                                // application of tensor-product based
+                                // elements for "symmetric" finite elements,
+                                // i.e., when the shape functions are
+                                // symmetric about 0.5 and the quadrature
+                                // points are, too. In that case, the 1D shape
+                                // values read (sorted lexicographically, rows
+                                // run over 1D dofs, columns over quadrature
+                                // points):
+                                // Q2 --> [ 0.687  0 -0.087 ]
+                                //        [ 0.4    1  0.4   ]
+                                //        [-0.087  0  0.687 ]
+                                // Q3 --> [ 0.66   0.003  0.002  0.049 ]
+                                //        [ 0.521  1.005 -0.01  -0.230 ]
+                                //        [-0.230 -0.01   1.005  0.521 ]
+                                //        [ 0.049  0.002  0.003  0.66  ]
+                                // Q4 --> [ 0.658  0.022  0 -0.007 -0.032 ]
+                                //        [ 0.608  1.059  0  0.039  0.176 ]
+                                //        [-0.409 -0.113  1 -0.113 -0.409 ]
+                                //        [ 0.176  0.039  0  1.059  0.608 ]
+                                //        [-0.032 -0.007  0  0.022  0.658 ]
+                                //
+                                // In these matrices, we want to use avoid
+                                // computations involving zeros and ones and
+                                // in addition use the symmetry in entries to
+                                // reduce the number of read operations.
+  const VectorizedArray<Number> * shape_values = this->data.shape_values.begin();
   for (int i2=0; i2<n_blocks2; ++i2)
   {
     for (int i1=0; i1<n_blocks1; ++i1)
     {
       for (int col=0; col<n_cols; ++col)
         {
-          VectorizedArray<Number> val0, val1, res0, res1;
+          VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
           if (dof_to_quad == true)
             {
-              val0 = this->data.shape_values[col];
-              val1 = this->data.shape_values[nn-1-col];
+              val0 = shape_values[col];
+              val1 = shape_values[nn-1-col];
             }
           else
             {
-              val0 = this->data.shape_values[col*n_q_points_1d];
-              val1 = this->data.shape_values[(col+1)*n_q_points_1d-1];
+              val0 = shape_values[col*n_q_points_1d];
+              val1 = shape_values[(col+1)*n_q_points_1d-1];
             }
           if (mid > 0)
             {
-              res0 = val0 * in[0];
-              res1 = val1 * in[0];
-              res0 += val1 * in[stride*(mm-1)];
-              res1 += val0 * in[stride*(mm-1)];
+              in0 = in[0];
+              in1 = in[stride*(mm-1)];
+              res0 = val0 * in0;
+              res1 = val1 * in0;
+              res0 += val1 * in1;
+              res1 += val0 * in1;
               for (int ind=1; ind<mid; ++ind)
                 {
                   if (dof_to_quad == true)
                     {
-                      val0 = this->data.shape_values[ind*n_q_points_1d+col];
-                      val1 = this->data.shape_values[ind*n_q_points_1d+nn-1-col];
+                      val0 = shape_values[ind*n_q_points_1d+col];
+                      val1 = shape_values[ind*n_q_points_1d+nn-1-col];
                     }
                   else
                     {
-                      val0 = this->data.shape_values[col*n_q_points_1d+ind];
-                      val1 = this->data.shape_values[(col+1)*n_q_points_1d-1-ind];
+                      val0 = shape_values[col*n_q_points_1d+ind];
+                      val1 = shape_values[(col+1)*n_q_points_1d-1-ind];
                     }
-                  res0 += val0 * in[stride*ind];
-                  res1 += val1 * in[stride*ind];
-                  res0 += val1 * in[stride*(mm-1-ind)];
-                  res1 += val0 * in[stride*(mm-1-ind)];
+                  in0 = in[stride*ind];
+                  in1 = in[stride*(mm-1-ind)];
+                  res0 += val0 * in0;
+                  res1 += val1 * in0;
+                  res0 += val1 * in1;
+                  res1 += val0 * in1;
                 }
             }
           else
@@ -4997,7 +4797,7 @@ apply_values (const VectorizedArray<Number> input [],
             {
               if (mm % 2 == 1)
                 {
-                  val0 = this->data.shape_values[mid*n_q_points_1d+col];
+                  val0 = shape_values[mid*n_q_points_1d+col];
                   val1 = val0 * in[stride*mid];
                   res0 += val1;
                   res1 += val1;
@@ -5007,7 +4807,7 @@ apply_values (const VectorizedArray<Number> input [],
             {
               if (mm % 2 == 1 && nn % 2 == 0)
                 {
-                  val0 = this->data.shape_values[col*n_q_points_1d+mid];
+                  val0 = shape_values[col*n_q_points_1d+mid];
                   val1 = val0 * in[stride*mid];
                   res0 += val1;
                   res1 += val1;
@@ -5034,14 +4834,14 @@ apply_values (const VectorizedArray<Number> input [],
       else if (dof_to_quad == true && nn%2==1)
         {
           VectorizedArray<Number> res0;
-          VectorizedArray<Number> val0  = this->data.shape_values[n_cols];
+          VectorizedArray<Number> val0  = shape_values[n_cols];
           if (mid > 0)
             {
               res0  = in[0] + in[stride*(mm-1)];
               res0 *= val0;
               for (int ind=1; ind<mid; ++ind)
                 {
-                  val0  = this->data.shape_values[ind*n_q_points_1d+n_cols];
+                  val0  = shape_values[ind*n_q_points_1d+n_cols];
                   VectorizedArray<Number> val1  = in[stride*ind] + in[stride*(mm-1-ind)];
                   val1 *= val0;
                   res0 += val1;
@@ -5051,7 +4851,7 @@ apply_values (const VectorizedArray<Number> input [],
             res0 = VectorizedArray<Number>();
           if (mm % 2 == 1)
             {
-              val0  = this->data.shape_values[mid*n_q_points_1d+n_cols];
+              val0  = shape_values[mid*n_q_points_1d+n_cols];
               res0 += val0 * in[stride*mid];
             }
           if (add == false)
@@ -5064,12 +4864,12 @@ apply_values (const VectorizedArray<Number> input [],
           VectorizedArray<Number> res0;
           if (mid > 0)
             {
-              VectorizedArray<Number> val0 = this->data.shape_values[n_cols*n_q_points_1d];
+              VectorizedArray<Number> val0 = shape_values[n_cols*n_q_points_1d];
               res0 = in[0] + in[stride*(mm-1)];
               res0 *= val0;
               for (int ind=1; ind<mid; ++ind)
                 {
-                  val0  = this->data.shape_values[n_cols*n_q_points_1d+ind];
+                  val0  = shape_values[n_cols*n_q_points_1d+ind];
                   VectorizedArray<Number> val1 = in[stride*ind] + in[stride*(mm-1-ind)];
                   val1 *= val0;
                   res0 += val1;
@@ -5114,14 +4914,14 @@ apply_values (const VectorizedArray<Number> input [],
 
 
 
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components,
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
           typename Number>
 template <int direction, bool dof_to_quad, bool add>
 inline
 void
-FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,Number>::
-apply_gradients (const VectorizedArray<Number> input [],
-                 VectorizedArray<Number>       output [])
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::apply_gradients (const VectorizedArray<Number> in [],
+                   VectorizedArray<Number>       out [])
 {
   AssertIndexRange (direction, dim);
   const int mm     = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
@@ -5131,50 +4931,74 @@ apply_gradients (const VectorizedArray<Number> input [],
 
   const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
   const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
-  const int stride    = ((direction > 0 ? nn : 1 ) *
-                         (direction > 1 ? nn : 1));
+  const int stride    = Utilities::fixed_int_power<nn,direction>::value;
 
-  const VectorizedArray<Number> * in = &input[0];
-  VectorizedArray<Number> * out = &output[0];
+  const VectorizedArray<Number> * shape_gradients = this->data.shape_gradients.begin();
   for (int i2=0; i2<n_blocks2; ++i2)
   {
     for (int i1=0; i1<n_blocks1; ++i1)
     {
+                                // For the specialized loop used for the
+                                // gradient computation in here, the 1D shape
+                                // values read (sorted lexicographically, rows
+                                // run over 1D dofs, columns over quadrature
+                                // points):
+                                // Q2 --> [-2.549 -1  0.549 ]
+                                //        [ 3.098  0 -3.098 ]
+                                //        [-0.549  1  2.549 ]
+                                // Q3 --> [-4.315 -1.03  0.5  -0.44  ]
+                                //        [ 6.07  -1.44 -2.97  2.196 ]
+                                //        [-2.196  2.97  1.44 -6.07  ]
+                                //        [ 0.44  -0.5   1.03  4.315 ]
+                                // Q4 --> [-6.316 -1.3    0.333 -0.353  0.413 ]
+                                //        [10.111 -2.76  -2.667  2.066 -2.306 ]
+                                //        [-5.688  5.773  0     -5.773  5.688 ]
+                                //        [ 2.306 -2.066  2.667  2.76 -10.111 ]
+                                //        [-0.413  0.353 -0.333 -0.353  0.413 ]
+                                //
+                                // In these matrices, we want to use avoid
+                                // computations involving zeros and ones and
+                                // in addition use the symmetry in entries to
+                                // reduce the number of read operations.
       for (int col=0; col<n_cols; ++col)
         {
-          VectorizedArray<Number> val0, val1, res0, res1;
+          VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
           if (dof_to_quad == true)
             {
-              val0 = this->data.shape_gradients[col];
-              val1 = this->data.shape_gradients[nn-1-col];
+              val0 = shape_gradients[col];
+              val1 = shape_gradients[nn-1-col];
             }
           else
             {
-              val0 = this->data.shape_gradients[col*n_q_points_1d];
-              val1 = this->data.shape_gradients[(nn-col-1)*n_q_points_1d];
+              val0 = shape_gradients[col*n_q_points_1d];
+              val1 = shape_gradients[(nn-col-1)*n_q_points_1d];
             }
           if (mid > 0)
             {
-              res0 = val0 * in[0];
-              res1 = val1 * in[0];
-              res0 -= val1 * in[stride*(mm-1)];
-              res1 -= val0 * in[stride*(mm-1)];
+              in0 = in[0];
+              in1 = in[stride*(mm-1)];
+              res0 = val0 * in0;
+              res1 = val1 * in0;
+              res0 -= val1 * in1;
+              res1 -= val0 * in1;
               for (int ind=1; ind<mid; ++ind)
                 {
                   if (dof_to_quad == true)
                     {
-                      val0 = this->data.shape_gradients[ind*n_q_points_1d+col];
-                      val1 = this->data.shape_gradients[ind*n_q_points_1d+nn-1-col];
+                      val0 = shape_gradients[ind*n_q_points_1d+col];
+                      val1 = shape_gradients[ind*n_q_points_1d+nn-1-col];
                     }
                   else
                     {
-                      val0 = this->data.shape_gradients[col*n_q_points_1d+ind];
-                      val1 = this->data.shape_gradients[(nn-col-1)*n_q_points_1d+ind];
+                      val0 = shape_gradients[col*n_q_points_1d+ind];
+                      val1 = shape_gradients[(nn-col-1)*n_q_points_1d+ind];
                     }
-                  res0 += val0 * in[stride*ind];
-                  res1 += val1 * in[stride*ind];
-                  res0 -= val1 * in[stride*(mm-1-ind)];
-                  res1 -= val0 * in[stride*(mm-1-ind)];
+                  in0 = in[stride*ind];
+                  in1 = in[stride*(mm-1-ind)];
+                  res0 += val0 * in0;
+                  res1 += val1 * in0;
+                  res0 -= val1 * in1;
+                  res1 -= val0 * in1;
                 }
             }
           else
@@ -5182,9 +5006,9 @@ apply_gradients (const VectorizedArray<Number> input [],
           if (mm % 2 == 1)
             {
               if (dof_to_quad == true)
-                val0 = this->data.shape_gradients[mid*n_q_points_1d+col];
+                val0 = shape_gradients[mid*n_q_points_1d+col];
               else
-                val0 = this->data.shape_gradients[col*n_q_points_1d+mid];
+                val0 = shape_gradients[col*n_q_points_1d+mid];
               val1 = val0 * in[stride*mid];
               res0 += val1;
               res1 -= val1;
@@ -5204,17 +5028,17 @@ apply_gradients (const VectorizedArray<Number> input [],
         {
           VectorizedArray<Number> val0, res0;
           if (dof_to_quad == true)
-            val0 = this->data.shape_gradients[n_cols];
+            val0 = shape_gradients[n_cols];
           else
-            val0 = this->data.shape_gradients[n_cols*n_q_points_1d];
+            val0 = shape_gradients[n_cols*n_q_points_1d];
           res0  = in[0] - in[stride*(mm-1)];
           res0 *= val0;
           for (int ind=1; ind<mid; ++ind)
             {
               if (dof_to_quad == true)
-                val0 = this->data.shape_gradients[ind*n_q_points_1d+n_cols];
+                val0 = shape_gradients[ind*n_q_points_1d+n_cols];
               else
-                val0 = this->data.shape_gradients[n_cols*n_q_points_1d+ind];
+                val0 = shape_gradients[n_cols*n_q_points_1d+ind];
               VectorizedArray<Number> val1  = in[stride*ind] - in[stride*(mm-1-ind)];
               val1 *= val0;
               res0 += val1;
@@ -5261,14 +5085,14 @@ apply_gradients (const VectorizedArray<Number> input [],
                                 // same symmetry relations hold. However, it
                                 // is not possible to omit some values that
                                 // are zero for the values
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components,
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
           typename Number>
 template <int direction, bool dof_to_quad, bool add>
 inline
 void
-FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,Number>::
-apply_hessians (const VectorizedArray<Number> input [],
-                  VectorizedArray<Number>       output [])
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::apply_hessians (const VectorizedArray<Number> in [],
+                  VectorizedArray<Number>       out [])
 {
   AssertIndexRange (direction, dim);
   const int mm     = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
@@ -5278,18 +5102,15 @@ apply_hessians (const VectorizedArray<Number> input [],
 
   const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
   const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
-  const int stride    = ((direction > 0 ? nn : 1 ) *
-                         (direction > 1 ? nn : 1));
+  const int stride    = Utilities::fixed_int_power<nn,direction>::value;
 
-  const VectorizedArray<Number> * in = &input[0];
-  VectorizedArray<Number> * out = &output[0];
   for (int i2=0; i2<n_blocks2; ++i2)
   {
     for (int i1=0; i1<n_blocks1; ++i1)
     {
       for (int col=0; col<n_cols; ++col)
         {
-          VectorizedArray<Number> val0, val1, res0, res1;
+          VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
           if (dof_to_quad == true)
             {
               val0 = this->data.shape_hessians[col];
@@ -5302,10 +5123,12 @@ apply_hessians (const VectorizedArray<Number> input [],
             }
           if (mid > 0)
             {
-              res0 = val0 * in[0];
-              res1 = val1 * in[0];
-              res0 += val1 * in[stride*(mm-1)];
-              res1 += val0 * in[stride*(mm-1)];
+              in0 = in[0];
+              in1 = in[stride*(mm-1)];
+              res0 = val0 * in0;
+              res1 = val1 * in0;
+              res0 += val1 * in1;
+              res1 += val0 * in1;
               for (int ind=1; ind<mid; ++ind)
                 {
                   if (dof_to_quad == true)
@@ -5318,10 +5141,12 @@ apply_hessians (const VectorizedArray<Number> input [],
                       val0 = this->data.shape_hessians[col*n_q_points_1d+ind];
                       val1 = this->data.shape_hessians[(col+1)*n_q_points_1d-1-ind];
                     }
-                  res0 += val0 * in[stride*ind];
-                  res1 += val1 * in[stride*ind];
-                  res0 += val1 * in[stride*(mm-1-ind)];
-                  res1 += val0 * in[stride*(mm-1-ind)];
+                  in0 = in[stride*ind];
+                  in1 = in[stride*(mm-1-ind)];
+                  res0 += val0 * in0;
+                  res1 += val1 * in0;
+                  res0 += val1 * in1;
+                  res1 += val0 * in1;
                 }
             }
           else
@@ -5413,13 +5238,13 @@ apply_hessians (const VectorizedArray<Number> input [],
 }
 
 
-/*----------------------- FEEvaluationGL -------------------------------*/
 
+/*------------------------- FEEvaluationGL ----------------------------------*/
 
-template <int dim, int fe_degree, int n_components, typename Number>
+template <int dim, int fe_degree, int n_components_, typename Number>
 inline
-FEEvaluationGL<dim,fe_degree,n_components,Number>::
-FEEvaluationGL (const MatrixFree<dim,Number> &data_in,
+FEEvaluationGL<dim,fe_degree,n_components_,Number>
+::FEEvaluationGL (const MatrixFree<dim,Number> &data_in,
                   const unsigned int fe_no,
                   const unsigned int quad_no)
   :
@@ -5429,6 +5254,7 @@ FEEvaluationGL (const MatrixFree<dim,Number> &data_in,
   std::string error_mess = "FEEvaluationGL not appropriate. It assumes:\n";
   error_mess += "   - identity operation for shape values\n";
   error_mess += "   - zero diagonal at interior points for gradients\n";
+  error_mess += "   - gradient equal to unity at element boundary\n";
   error_mess += "Try FEEvaluation<...> instead!";
 
   const double zero_tol =
@@ -5451,16 +5277,21 @@ FEEvaluationGL (const MatrixFree<dim,Number> &data_in,
   for (unsigned int i=1; i<n_points_1d-1; ++i)
     Assert (std::fabs(this->data.shape_gradients[i*n_points_1d+i][0])<zero_tol,
             ExcMessage (error_mess.c_str()));
+  Assert (std::fabs(this->data.shape_gradients[n_points_1d-1][0]-
+                    (n_points_1d%2==0 ? -1. : 1.)) < zero_tol,
+          ExcMessage (error_mess.c_str()));
 #endif
 }
 
 
 
-template <int dim, int fe_degree, int n_components, typename Number>
+template <int dim, int fe_degree, int n_components_, typename Number>
 inline
 void
-FEEvaluationGL<dim,fe_degree,n_components,Number>::
-evaluate (bool evaluate_val,bool evaluate_grad,bool evaluate_lapl)
+FEEvaluationGL<dim,fe_degree,n_components_,Number>
+::evaluate (const bool evaluate_val,
+            const bool evaluate_grad,
+            const bool evaluate_lapl)
 {
   Assert (this->cell != numbers::invalid_unsigned_int,
           ExcNotInitialized());
@@ -5476,6 +5307,9 @@ evaluate (bool evaluate_val,bool evaluate_grad,bool evaluate_lapl)
       this->values_quad_initialized = true;
 #endif
     }
+                                // separate implementation here compared to
+                                // the general case because the values are an
+                                // identity operation
   if (evaluate_grad == true)
     {
       for(unsigned int comp=0;comp<n_components;comp++)
@@ -5560,11 +5394,11 @@ evaluate (bool evaluate_val,bool evaluate_grad,bool evaluate_lapl)
 
 
 
-template <int dim, int fe_degree, int n_components, typename Number>
+template <int dim, int fe_degree, int n_components_, typename Number>
 inline
 void
-FEEvaluationGL<dim,fe_degree,n_components,Number>::
-integrate (bool integrate_val, bool integrate_grad)
+FEEvaluationGL<dim,fe_degree,n_components_,Number>
+::integrate (const bool integrate_val, const bool integrate_grad)
 {
   Assert (this->cell != numbers::invalid_unsigned_int,
           ExcNotInitialized());
@@ -5634,13 +5468,13 @@ integrate (bool integrate_val, bool integrate_grad)
 
 
 
-template <int dim, int fe_degree, int n_components, typename Number>
+template <int dim, int fe_degree, int n_components_, typename Number>
 template <int direction, bool dof_to_quad, bool add>
 inline
 void
-FEEvaluationGL<dim,fe_degree,n_components,Number>::
-apply_gradients (const VectorizedArray<Number> input [],
-                 VectorizedArray<Number>       output [])
+FEEvaluationGL<dim,fe_degree,n_components_,Number>
+::apply_gradients (const VectorizedArray<Number> in [],
+                   VectorizedArray<Number>       out [])
 {
   AssertIndexRange (direction, dim);
   const int mm     = fe_degree+1;
@@ -5650,34 +5484,75 @@ apply_gradients (const VectorizedArray<Number> input [],
 
   const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
   const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
-  const int stride    = ((direction > 0 ? nn : 1 ) *
-                         (direction > 1 ? nn : 1));
-
-  const VectorizedArray<Number> * in = &input[0];
-  VectorizedArray<Number> * out = &output[0];
+  const int stride    = Utilities::fixed_int_power<nn,direction>::value;
+
+                                // This loop specializes the application of
+                                // the tensor product loop for Gauss-Lobatto
+                                // elements which are symmetric about 0.5 just
+                                // as the general class of elements treated by
+                                // FEEvaluation, have diagonal shape matrices
+                                // for the values and have the following
+                                // gradient matrices (notice the zeros on the
+                                // diagonal in the interior points, which is
+                                // due to the construction of Legendre
+                                // polynomials):
+                                // Q2 --> [-3 -1  1 ]
+                                //        [ 4  0 -4 ]
+                                //        [-1  1  3 ]
+                                // Q3 --> [-6    -1.618  0.618 -1    ]
+                                //        [ 8.09  0     -2.236  3.09 ]
+                                //        [-3.09  2.236  0     -8.09 ]
+                                //        [ 1    -0.618  1.618  6    ]
+                                // Q4 --> [-10    -2.482  0.75  -0.518  1     ]
+                                //        [ 13.51  0     -2.673  1.528 -2.82  ]
+                                //        [-5.333  3.491  0     -3.491  5.333 ]
+                                //        [ 2.82  -1.528  2.673  0    -13.51  ]
+                                //        [-1      0.518 -0.75   2.482 10     ]
   for (int i2=0; i2<n_blocks2; ++i2)
   {
     for (int i1=0; i1<n_blocks1; ++i1)
     {
       for (int col=0; col<n_cols; ++col)
         {
-          VectorizedArray<Number> val0, val1, res0, res1;
-          if (dof_to_quad == true)
-            {
-              val0 = this->data.shape_gradients[col];
-              val1 = this->data.shape_gradients[nn-1-col];
-            }
-          else
-            {
-              val0 = this->data.shape_gradients[col*mm];
-              val1 = this->data.shape_gradients[(nn-col-1)*mm];
-            }
+          VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
           if (mid > 0)
             {
-              res0 = val0 * in[0];
-              res1 = val1 * in[0];
-              res0 -= val1 * in[stride*(mm-1)];
-              res1 -= val0 * in[stride*(mm-1)];
+              if (dof_to_quad == true)
+                {
+                  val0 = this->data.shape_gradients[col];
+                  val1 = this->data.shape_gradients[nn-1-col];
+                }
+              else
+                {
+                  val0 = this->data.shape_gradients[col*mm];
+                  val1 = this->data.shape_gradients[(nn-col-1)*mm];
+                }
+              in0 = in[0];
+              in1 = in[stride*(mm-1)];
+              if (col == 0)
+                {
+                  if ((mm+dof_to_quad)%2 == 1)
+                    {
+                      res0 = val0 * in0;
+                      res1 = -in0;
+                      res0 += in1;
+                      res1 -= val0 * in1;
+                    }
+                  else
+                    {
+                      res0 = val0 * in0;
+                      res0 -= in1;
+                      res1 = in0;
+                      res1 -= val0 * in1;
+                    }
+                }
+              else
+                {
+                  res0 = val0 * in0;
+                  res1 = val1 * in0;
+                  res0 -= val1 * in1;
+                  res1 -= val0 * in1;
+                }
               for (int ind=1; ind<mid; ++ind)
                 {
                   if (dof_to_quad == true)
@@ -5693,17 +5568,19 @@ apply_gradients (const VectorizedArray<Number> input [],
 
                                 // at inner points, the gradient is zero for
                                 // ind==col
+                  in0 = in[stride*ind];
+                  in1 = in[stride*(mm-1-ind)];
                   if (ind == col)
                     {
-                      res1 += val1 * in[stride*ind];
-                      res0 -= val1 * in[stride*(mm-1-ind)];
+                      res1 += val1 * in0;
+                      res0 -= val1 * in1;
                     }
                   else
                     {
-                      res0 += val0 * in[stride*ind];
-                      res1 += val1 * in[stride*ind];
-                      res0 -= val1 * in[stride*(mm-1-ind)];
-                      res1 -= val0 * in[stride*(mm-1-ind)];
+                      res0 += val0 * in0;
+                      res1 += val1 * in0;
+                      res0 -= val1 * in1;
+                      res1 -= val0 * in1;
                     }
                 }
             }
index 1848193b73b253b935dad0dec38f1fd15ee3b9df..d8913efc2b5fc799a01be9d902bb784d6cdae4a0 100644 (file)
@@ -261,14 +261,15 @@ namespace MatrixFreeFunctions
                                  * underlying the problem (constructed from a
                                  * 1D tensor product quadrature formula).
                                  */
-      dealii::hp::QCollection<dim>    quadrature_formula;
+      dealii::hp::QCollection<dim>    quadrature;
 
                                 /**
                                  * The (dim-1)-dimensional quadrature formula
-                                 * underlying the problem (constructed from a
-                                 * 1D tensor product quadrature formula).
+                                 * corresponding to face evaluation
+                                 * (constructed from a 1D tensor product
+                                 * quadrature formula).
                                  */
-      dealii::hp::QCollection<dim-1>  quadrature_formula_faces;
+      dealii::hp::QCollection<dim-1>  face_quadrature;
 
                                 /**
                                  * The number of quadrature points for the
index c607df65e6510a67ccc91978f99c705c7d1fec43..9d1626f45a9e08af1301172e24433ece5a51d064 100644 (file)
@@ -130,11 +130,11 @@ namespace MatrixFreeFunctions
   template <int dim, typename Number>
   void
   MappingInfo<dim,Number>::initialize
-  (const dealii::Triangulation<dim>                                 &tria,
+  (const dealii::Triangulation<dim>                         &tria,
    const std::vector<std::pair<unsigned int,unsigned int> > &cells,
    const std::vector<unsigned int>                          &active_fe_index,
    const Mapping<dim>                                       &mapping,
-   const std::vector<dealii::hp::QCollection<1> >                   &quad,
+   const std::vector<dealii::hp::QCollection<1> >           &quad,
    const UpdateFlags                                         update_flags_input)
   {
     clear();
@@ -204,16 +204,16 @@ namespace MatrixFreeFunctions
 
             current_data.n_q_points_face.push_back 
               (Utilities::fixed_power<dim-1>(n_q_points_1d[q]));
-            current_data.quadrature_formula.push_back
+            current_data.quadrature.push_back
               (Quadrature<dim>(quad[my_q][q]));
-            current_data.quadrature_formula_faces.push_back
+            current_data.face_quadrature.push_back
               (Quadrature<dim-1>(quad[my_q][q]));
 
                                 // set quadrature weights in vectorized form
             current_data.quadrature_weights[q].resize(n_q_points);
             for (unsigned int i=0; i<n_q_points; ++i)
               current_data.quadrature_weights[q][i] =
-                current_data.quadrature_formula[q].get_weights()[i];
+                current_data.quadrature[q].get_weights()[i];
 
             if (n_hp_quads > 1)
               current_data.quad_index_conversion[q] = n_q_points;
@@ -260,7 +260,7 @@ namespace MatrixFreeFunctions
                                 // finite element, so just hold a vector of
                                 // FEValues
         std::vector<std_cxx1x::shared_ptr<FEValues<dim> > >
-          fe_values (current_data.quadrature_formula.size());
+          fe_values (current_data.quadrature.size());
         UpdateFlags update_flags_feval =
           (update_flags & update_inverse_jacobians ? update_jacobians : update_default) |
           (update_flags & update_jacobian_grads ? update_jacobian_grads : update_default) |
@@ -307,7 +307,7 @@ namespace MatrixFreeFunctions
             if (fe_values[fe_index].get() == 0)
               fe_values[fe_index].reset 
                 (new FEValues<dim> (mapping, dummy_fe,
-                                    current_data.quadrature_formula[fe_index],
+                                    current_data.quadrature[fe_index],
                                     update_flags_feval));
             FEValues<dim> &fe_val = *fe_values[fe_index];
             data.resize (n_q_points);
@@ -912,8 +912,8 @@ namespace MatrixFreeFunctions
     memory += MemoryConsumption::memory_consumption (jacobians_grad_upper);
     memory += MemoryConsumption::memory_consumption (rowstart_q_points);
     memory += MemoryConsumption::memory_consumption (quadrature_points);
-    memory += MemoryConsumption::memory_consumption (quadrature_formula);
-    memory += MemoryConsumption::memory_consumption (quadrature_formula_faces);
+    memory += MemoryConsumption::memory_consumption (quadrature);
+    memory += MemoryConsumption::memory_consumption (face_quadrature);
     memory += MemoryConsumption::memory_consumption (quadrature_weights);
     memory += MemoryConsumption::memory_consumption (n_q_points);
     memory += MemoryConsumption::memory_consumption (n_q_points_face);
index 8c7a87b2e9108775941291af87e880426dfdf7fb..bacc960834a45ecdb75d41262e656a232231c2d3 100644 (file)
@@ -932,8 +932,16 @@ public:
                                  * given hp index.
                                  */
   const Quadrature<dim> &
-  get_quad (const unsigned int quad_index = 0,
-            const unsigned int hp_active_fe_index = 0) const;
+  get_quadrature (const unsigned int quad_index = 0,
+                  const unsigned int hp_active_fe_index = 0) const;
+
+                                /**
+                                 * Returns the quadrature rule for
+                                 * given hp index.
+                                 */
+  const Quadrature<dim-1> &
+  get_face_quadrature (const unsigned int quad_index = 0,
+                       const unsigned int hp_active_fe_index = 0) const;
 
                                 /**
                                  * Queries whether or not the
@@ -1352,8 +1360,15 @@ MatrixFree<dim,Number>::create_cell_subrange_hp
  const unsigned int degree,
  const unsigned int vector_component) const
 {
-  if (dof_info[vector_component].cell_active_fe_index.size() == 0)
-    return range;
+  AssertIndexRange (vector_component, dof_info.size());
+  if (dof_info[vector_component].cell_active_fe_index.empty())
+    {
+      AssertDimension (dof_info[vector_component].fe_index_conversion.size(),1);
+      if (dof_info[vector_component].fe_index_conversion[0].first == degree)
+        return range;
+      else
+        return std::pair<unsigned int,unsigned int> (range.second,range.second);
+    }
 
   const unsigned int fe_index = 
     dof_info[vector_component].fe_index_from_degree(degree);
@@ -1691,12 +1706,25 @@ MatrixFree<dim,Number>::get_shape_info (const unsigned int index_fe,
 template <int dim, typename Number>
 inline
 const Quadrature<dim> &
-MatrixFree<dim,Number>::get_quad (const unsigned int quad_index,
-                                  const unsigned int active_fe_index) const
+MatrixFree<dim,Number>::get_quadrature (const unsigned int quad_index,
+                                        const unsigned int active_fe_index) const
+{
+  AssertIndexRange (quad_index, mapping_info.mapping_data_gen.size());
+  return mapping_info.mapping_data_gen[quad_index].
+    quadrature[active_fe_index];
+}
+
+
+
+template <int dim, typename Number>
+inline
+const Quadrature<dim-1> &
+MatrixFree<dim,Number>::get_face_quadrature (const unsigned int quad_index,
+                                             const unsigned int active_fe_index) const
 {
   AssertIndexRange (quad_index, mapping_info.mapping_data_gen.size());
   return mapping_info.mapping_data_gen[quad_index].
-    quadrature_formula[active_fe_index];
+    face_quadrature[active_fe_index];
 }
 
 
index 675f5679aef61ba4934c2b6add53cf4a518e1fa4..1fc7d26320a73c3c42dab253205180f8636bc82e 100644 (file)
@@ -184,7 +184,12 @@ internal_reinit(const Mapping<dim>                         &mapping,
                                 // cell, Jacobian determinants, quadrature
                                 // points in real space, based on the ordering
                                 // of the cells determined in @p
-                                // extract_local_to_global_indices.
+                                // extract_local_to_global_indices. The
+                                // algorithm assumes that the active FE index
+                                // for the transformations is given the active
+                                // FE index in the zeroth DoFHandler. TODO:
+                                // how do things look like in the more general
+                                // case?
   if(additional_data.initialize_mapping == true)
     {
       mapping_info.initialize (dof_handler[0]->get_tria(), cell_level_index,
@@ -486,14 +491,15 @@ void MatrixFree<dim,Number>::initialize_indices
           for (unsigned int f=0; f<fe.size(); ++f)
             fes.push_back (&fe[f]);
 
-          dof_info[no].cell_active_fe_index.resize(n_active_cells,
-                                                   numbers::invalid_unsigned_int);
           dof_info[no].max_fe_index = fe.size();
           dof_info[no].fe_index_conversion.resize (fe.size());
           for (unsigned int ind=0; ind<hpdof->get_fe().size(); ++ind)
             dof_info[no].fe_index_conversion[ind] =
               std::pair<unsigned int,unsigned int>(fe[ind].degree,
                                                    fe[ind].dofs_per_cell);
+          if (fe.size() > 1)
+            dof_info[no].cell_active_fe_index.resize(n_active_cells,
+                                                     numbers::invalid_unsigned_int);
         }
       else
         {
@@ -502,6 +508,10 @@ void MatrixFree<dim,Number>::initialize_indices
             &*dof_handlers.dof_handler[no] : &*dof_handlers.mg_dof_handler[no];
           fes.push_back (&dofh->get_fe());
           dof_info[no].max_fe_index = 1;
+          dof_info[no].fe_index_conversion.resize (1);
+          dof_info[no].fe_index_conversion[0] =
+            std::pair<unsigned int,unsigned int>(fes.back()->degree,
+                                                 fes.back()->dofs_per_cell);
         }
 
       lexicographic_inv[no].resize (fes.size());
@@ -654,8 +664,9 @@ void MatrixFree<dim,Number>::initialize_indices
                          cell_level_index[counter].first,
                          cell_level_index[counter].second,
                          dofh);
-              dof_info[no].cell_active_fe_index[counter] =
-                cell_it->active_fe_index();
+              if (dofh->get_fe().size() > 1)
+                dof_info[no].cell_active_fe_index[counter] =
+                  cell_it->active_fe_index();
               local_dof_indices.resize (cell_it->get_fe().dofs_per_cell);
               cell_it->get_dof_indices(local_dof_indices);
               dof_info[no].read_dof_indices (local_dof_indices,
index 913f3220b71d9226fde9cba1f5cf037a847a6b66..5c9d3882e114f6cc5a9b81e2d6de5d7cd0ca9e48 100644 (file)
@@ -17,36 +17,36 @@ for (deal_II_dimension : DIMENSIONS)
   template class MatrixFree<deal_II_dimension,float>;
 
   // reinit for DoFHandler
-  template void MatrixFree<deal_II_dimension,double>::reinit
+  template void MatrixFree<deal_II_dimension,double>::internal_reinit
   (const Mapping<deal_II_dimension>                           &,
    const std::vector<const DoFHandler<deal_II_dimension>*>    &,
    const std::vector<const ConstraintMatrix*>                 &,
    const std::vector<IndexSet>                                &,
-   const std::vector<Quadrature<1> >                          &,
+   const std::vector<hp::QCollection<1> >                     &,
    const MatrixFree<deal_II_dimension,double>::AdditionalData);
-  template void MatrixFree<deal_II_dimension,float>::reinit
+  template void MatrixFree<deal_II_dimension,float>::internal_reinit
   (const Mapping<deal_II_dimension>                           &,
    const std::vector<const DoFHandler<deal_II_dimension>*>    &,
    const std::vector<const ConstraintMatrix*>                 &,
    const std::vector<IndexSet>                                &,
-   const std::vector<Quadrature<1> >                          &,
+   const std::vector<hp::QCollection<1> >                     &,
    const MatrixFree<deal_II_dimension,float>::AdditionalData);
 
 
   // reinit for MGDoFHandler
-  template void MatrixFree<deal_II_dimension,double>::reinit
+  template void MatrixFree<deal_II_dimension,double>::internal_reinit
   (const Mapping<deal_II_dimension>                           &,
    const std::vector<const MGDoFHandler<deal_II_dimension>*>  &,
    const std::vector<const ConstraintMatrix*>                 &,
    const std::vector<IndexSet>                                &,
-   const std::vector<Quadrature<1> >                          &,
+   const std::vector<hp::QCollection<1> >                     &,
    const MatrixFree<deal_II_dimension,double>::AdditionalData);
-  template void MatrixFree<deal_II_dimension,float>::reinit
+  template void MatrixFree<deal_II_dimension,float>::internal_reinit
   (const Mapping<deal_II_dimension>                           &,
    const std::vector<const MGDoFHandler<deal_II_dimension>*>  &,
    const std::vector<const ConstraintMatrix*>                 &,
    const std::vector<IndexSet>                                &,
-   const std::vector<Quadrature<1> >                          &,
+   const std::vector<hp::QCollection<1> >                     &,
    const MatrixFree<deal_II_dimension,float>::AdditionalData);
 
   template void MatrixFree<deal_II_dimension,double>::
index 122db076614a5c1232e58e1fa60f82630dac1c27..f942ffd0268e33a9e092fbe3ed4dc67b87134cc7 100644 (file)
@@ -42,12 +42,10 @@ template <int dim, int fe_degree, int n_q_points_1d=fe_degree+1, typename Number
 class MatrixFreeTest
 {
  public:
-  static const std::size_t n_vectors = VectorizedArray<Number>::n_vectors;
-
   MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
     data   (data_in),
     fe_val (data.get_dof_handler().get_fe(),
-            Quadrature<dim>(data.get_quad(0)),
+            Quadrature<dim>(data.get_quadrature(0)),
             update_values | update_gradients | update_hessians)
   {};
 
@@ -55,7 +53,7 @@ class MatrixFreeTest
              const Mapping<dim>               &mapping):
     data   (data_in),
     fe_val (mapping, data.get_dof_handler().get_fe(),
-            Quadrature<dim>(data.get_quad(0)),
+            Quadrature<dim>(data.get_quadrature(0)),
             update_values | update_gradients | update_hessians)
   {};
 
@@ -70,7 +68,6 @@ class MatrixFreeTest
                const Vector<Number> &src,
                const std::pair<unsigned int,unsigned int> &cell_range) const
   {
-    typedef VectorizedArray<Number> vector_t;
     FEEvaluation<dim,fe_degree,n_q_points_1d,1,Number> fe_eval (data);
 
     std::vector<double> reference_values (fe_eval.n_q_points);
@@ -198,9 +195,6 @@ template <int dim, int fe_degree,typename Number>
 class MatrixFreeTest<dim,fe_degree,0,Number>
 {
  public:
-  typedef VectorizedArray<Number> vector_t;
-  static const std::size_t n_vectors = VectorizedArray<Number>::n_vectors;
-
   MatrixFreeTest(const MatrixFree<dim,Number> &)
   {};
 
index 3d80aa918211e9dddf4ac7903eda51ff33c7fcb6..62a7d65f2cf53192a0ef3e7d79b7ff92e9d86d37 100644 (file)
@@ -23,8 +23,6 @@ template <int dim, int fe_degree, typename Number>
 class MatrixFreeTestGL : public MatrixFreeTest<dim, fe_degree, fe_degree+1, Number>
 {
  public:
-  typedef VectorizedArray<Number> vector_t;
-  static const std::size_t n_vectors = VectorizedArray<Number>::n_array_elements;
 
   MatrixFreeTestGL(const MatrixFree<dim,Number> &data,
                const Mapping<dim>               &mapping):
index 7ecb02f458ee4282a1830a9a1fe20586f53c0531..dce1bb24e1b01fd9b63096d4582c1384decf4b70 100644 (file)
@@ -38,17 +38,15 @@ template <int dim, int fe_degree, int n_q_points_1d=fe_degree+1, typename Number
 class MatrixFreeTest
 {
  public:
-  typedef VectorizedArray<Number> vector_t;
   typedef std::vector<Vector<Number> > VectorType;
-  static const std::size_t n_vectors = VectorizedArray<Number>::n_array_elements;
 
   MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
     data    (data_in),
     fe_val0 (data.get_dof_handler(0).get_fe(),
-             Quadrature<dim>(data.get_quad(0)),
+             Quadrature<dim>(data.get_quadrature(0)),
              update_values | update_gradients | update_hessians),
     fe_val1 (data.get_dof_handler(1).get_fe(),
-             Quadrature<dim>(data.get_quad(1)),
+             Quadrature<dim>(data.get_quadrature(1)),
              update_values | update_gradients | update_hessians)
   {};
 
index 61499d7c9bee50fce48a4a4e308126413daa01a0..8c26693f5f7d0b0c56a9bb9a835580777883bda8 100644 (file)
@@ -40,20 +40,18 @@ template <int dim, int fe_degree, int n_q_points_1d=fe_degree+1, typename Number
 class MatrixFreeTest
 {
  public:
-  typedef VectorizedArray<Number> vector_t;
   typedef std::vector<Vector<Number> > VectorType;
-  static const std::size_t n_vectors = VectorizedArray<Number>::n_array_elements;
 
   MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
     data    (data_in),
     fe_val0 (data.get_dof_handler(0).get_fe(),
-             Quadrature<dim>(data.get_quad(0)),
+             Quadrature<dim>(data.get_quadrature(0)),
              update_values | update_gradients | update_hessians),
     fe_val1 (data.get_dof_handler(1).get_fe(),
-             Quadrature<dim>(data.get_quad(1)),
+             Quadrature<dim>(data.get_quadrature(1)),
              update_values | update_gradients | update_hessians),
     fe_val2 (data.get_dof_handler(2).get_fe(),
-             Quadrature<dim>(data.get_quad(1)),
+             Quadrature<dim>(data.get_quadrature(1)),
              update_values | update_gradients | update_hessians)
   {};
 
index b1ca17475b2ca3d60478ef568a6a9d143b2ecbff..9e2042debfcaf18127d3ff338d02d001a7ad5d07 100644 (file)
@@ -25,9 +25,6 @@ template <int dim, int fe_degree, typename Number>
 class MatrixFreeTestGen : public MatrixFreeTest<dim, fe_degree, fe_degree+1, Number>
 {
  public:
-  typedef VectorizedArray<Number> vector_t;
-  static const std::size_t n_vectors = VectorizedArray<Number>::n_array_elements;
-
   MatrixFreeTestGen(const MatrixFree<dim,Number> &data,
                 const Mapping<dim>               &mapping):
     MatrixFreeTest<dim, fe_degree, fe_degree+1, Number>(data, mapping)
index 5cb26618e439b3c8e24e581e8b6cead4c4b72d1e..86b2ec2dc03be5e5813c3e2a3ef237cbf5f306e4 100644 (file)
@@ -51,7 +51,7 @@ void test ()
   tria.begin(tria.n_levels()-1)->set_refine_flag();
   tria.last()->set_refine_flag();
   tria.execute_coarsening_and_refinement();
-  tria.refine_global (4-dim);
+  tria.refine_global (1);
 
   FE_Q<dim> fe (fe_degree);
   DoFHandler<dim> dof (tria);
index 5dd4c9d91c581a0d05cdafb67ada47238d88a62a..68035a30b895c2e18aaa92719b8b5f091b469ba6 100644 (file)
@@ -37,9 +37,7 @@ template <int dim, int fe_degree, typename Number>
 class MatrixFreeTest
 {
  public:
-  typedef VectorizedArray<Number> vector_t;
   typedef Vector<Number> VectorType;
-  static const std::size_t n_vectors = VectorizedArray<Number>::n_array_elements;
 
   MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
     data (data_in)
@@ -121,7 +119,7 @@ operator () (const MatrixFree<dim,Number> &data,
                                 // FEEvaluations. Those are tested in other
                                 // functions and seen as reference here
       for (unsigned int q=0; q<fe_eval.n_q_points; ++q)
-        for (unsigned int j=0; j<n_vectors; ++j)
+        for (unsigned int j=0; j<VectorizedArray<Number>::n_array_elements; ++j)
           {
             errors[0] += std::fabs(fe_eval.get_value(q)[j]-
                                    fe_eval2.get_value(q)[j]);
index 709f985288ba16c0ba3cdfb5379614c4a54876c0..5c97f82bb2f05e4d97f56e61012a79ef7d7b94b2 100644 (file)
@@ -38,8 +38,6 @@ template <int dim, int fe_degree, int n_q_points_1d=fe_degree+1, typename Number
 class MatrixFreeTest
 {
  public:
-  static const std::size_t n_vectors = VectorizedArray<Number>::n_array_elements;
-
   MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
     data   (data_in)
   {};
@@ -52,8 +50,6 @@ class MatrixFreeTest
                const Vector<Number> &src,
                const std::pair<unsigned int,unsigned int> &cell_range) const
   {
-    typedef VectorizedArray<Number> vector_t;
-    const unsigned int n_vectors = sizeof(vector_t)/sizeof(Number);
     FEEvaluation<dim,fe_degree,n_q_points_1d,1,Number> fe_eval (data);
     FEEvaluation<dim,fe_degree,n_q_points_1d,1,Number> fe_eval_plain (data);
     for(unsigned int cell=cell_range.first;cell<cell_range.second;++cell)
@@ -65,7 +61,7 @@ class MatrixFreeTest
         fe_eval_plain.read_dof_values_plain(src);
 
         for (unsigned int i=0; i<fe_eval.dofs_per_cell; ++i)
-          for (unsigned int j=0; j<n_vectors; ++j)
+          for (unsigned int j=0; j<VectorizedArray<Number>::n_array_elements; ++j)
             {
               error += std::fabs(fe_eval.get_dof_value(i)[j]-
                                  fe_eval_plain.get_dof_value(i)[j]);
index 12dce33f76179d8802001b18f6a20ffa25c98c5a..73cf3a85da6f8f28eff3150f6a17172008b6647e 100644 (file)
@@ -36,14 +36,12 @@ template <int dim, int fe_degree, typename Number>
 class MatrixFreeTest
 {
  public:
-  typedef VectorizedArray<Number> vector_t;
   typedef std::vector<Vector<Number>*> VectorType;
-  static const std::size_t n_vectors = VectorizedArray<Number>::n_array_elements;
 
   MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
     data   (data_in),
     fe_val (data.get_dof_handler().get_fe(),
-            Quadrature<dim>(data.get_quad(0)),
+            Quadrature<dim>(data.get_quadrature(0)),
             update_values | update_gradients | update_JxW_values)
   {};
 
@@ -83,8 +81,8 @@ operator () (const MatrixFree<dim,Number> &data,
   FEEvaluation<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
   const unsigned int n_q_points = fe_eval.n_q_points;
   const unsigned int dofs_per_cell = fe_eval.dofs_per_cell;
-  AlignedVector<vector_t> values (n_q_points);
-  AlignedVector<vector_t> gradients (dim*n_q_points);
+  AlignedVector<VectorizedArray<Number> > values (n_q_points);
+  AlignedVector<VectorizedArray<Number> > gradients (dim*n_q_points);
   std::vector<unsigned int> dof_indices (dofs_per_cell);
   for(unsigned int cell=cell_range.first;cell<cell_range.second;++cell)
     {
@@ -121,7 +119,7 @@ operator () (const MatrixFree<dim,Number> &data,
       for (unsigned int q=0; q<n_q_points; ++q)
         {
           fe_eval.submit_value (values[q], q);
-          Tensor<1,dim,vector_t> submit (false);
+          Tensor<1,dim,VectorizedArray<Number> > submit (false);
           for (unsigned int d=0; d<dim; ++d)
             submit[d] = gradients[q*dim+d];
           fe_eval.submit_gradient (submit, q);
index 5796664c9b74565a50dd373878fe87e21295270d..d8788de6220fb434efc6ff669c39cb11054ba994 100644 (file)
@@ -37,20 +37,18 @@ template <int dim, int fe_degree, typename Number>
 class MatrixFreeTest
 {
  public:
-  typedef VectorizedArray<Number> vector_t;
   typedef std::vector<Vector<Number> > VectorType;
-  static const std::size_t n_vectors = VectorizedArray<Number>::n_array_elements;
 
   MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
     data   (data_in),
     fe_val0 (data.get_dof_handler(0).get_fe(),
-             Quadrature<dim>(data.get_quad(0)),
+             Quadrature<dim>(data.get_quadrature(0)),
              update_values | update_gradients | update_JxW_values),
     fe_val01 (data.get_dof_handler(0).get_fe(),
-              Quadrature<dim>(data.get_quad(1)),
+              Quadrature<dim>(data.get_quadrature(1)),
               update_values | update_gradients | update_JxW_values),
     fe_val1 (data.get_dof_handler(1).get_fe(),
-             Quadrature<dim>(data.get_quad(1)),
+             Quadrature<dim>(data.get_quadrature(1)),
              update_values | update_gradients | update_JxW_values)
   {};
 
@@ -92,10 +90,10 @@ operator () (const MatrixFree<dim,Number> &data,
   const unsigned int n_q_points1 = fe_eval1.n_q_points;
   const unsigned int dofs_per_cell0 = fe_eval0.dofs_per_cell;
   const unsigned int dofs_per_cell1 = fe_eval1.dofs_per_cell;
-  AlignedVector<vector_t> values0 (n_q_points0);
-  AlignedVector<vector_t> gradients0 (dim*n_q_points0);
-  AlignedVector<vector_t> values1 (n_q_points1);
-  AlignedVector<vector_t> gradients1 (dim*n_q_points1);
+  AlignedVector<VectorizedArray<Number> > values0 (n_q_points0);
+  AlignedVector<VectorizedArray<Number> > gradients0 (dim*n_q_points0);
+  AlignedVector<VectorizedArray<Number> > values1 (n_q_points1);
+  AlignedVector<VectorizedArray<Number> > gradients1 (dim*n_q_points1);
   std::vector<unsigned int> dof_indices0 (dofs_per_cell0);
   std::vector<unsigned int> dof_indices1 (dofs_per_cell1);
   for(unsigned int cell=cell_range.first;cell<cell_range.second;++cell)
@@ -177,7 +175,7 @@ operator () (const MatrixFree<dim,Number> &data,
       for (unsigned int q=0; q<n_q_points0; ++q)
         {
           fe_eval0.submit_value (values0[q], q);
-          Tensor<1,dim,vector_t> submit (false);
+          Tensor<1,dim,VectorizedArray<Number> > submit (false);
           for (unsigned int d=0; d<dim; ++d)
             submit[d] = gradients0[q*dim+d];
           fe_eval0.submit_gradient (submit, q);
@@ -189,7 +187,7 @@ operator () (const MatrixFree<dim,Number> &data,
       for (unsigned int q=0; q<n_q_points1; ++q)
         {
           fe_eval1.submit_value (values1[q], q);
-          Tensor<1,dim,vector_t> submit (false);
+          Tensor<1,dim,VectorizedArray<Number> > submit (false);
           for (unsigned int d=0; d<dim; ++d)
             submit[d] = gradients1[q*dim+d];
           fe_eval1.submit_gradient (submit, q);
@@ -201,7 +199,7 @@ operator () (const MatrixFree<dim,Number> &data,
       for (unsigned int q=0; q<n_q_points1; ++q)
         {
           fe_eval01.submit_value (values1[q], q);
-          Tensor<1,dim,vector_t> submit (false);
+          Tensor<1,dim,VectorizedArray<Number> > submit (false);
           for (unsigned int d=0; d<dim; ++d)
             submit[d] = gradients1[q*dim+d];
           fe_eval01.submit_gradient (submit, q);
index c94c0077b96ed47897b383a96018f86e70e7c37c..6af0e155aa8a72ebb7597fbc5366a2f8ef87e47a 100644 (file)
@@ -29,6 +29,8 @@
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/numerics/vectors.h>
 
+#include "create_mesh.h"
+
 #include <iostream>
 
 std::ofstream logfile("integrate_functions_multife2/output");
@@ -38,20 +40,18 @@ template <int dim, int fe_degree, typename Number>
 class MatrixFreeTest
 {
  public:
-  typedef VectorizedArray<Number> vector_t;
   typedef std::vector<Vector<Number> > VectorType;
-  static const std::size_t n_vectors = VectorizedArray<Number>::n_array_elements;
 
   MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
     data   (data_in),
     fe_val0 (data.get_dof_handler(0).get_fe(),
-             Quadrature<dim>(data.get_quad(0)),
+             Quadrature<dim>(data.get_quadrature(0)),
              update_values | update_gradients | update_JxW_values),
     fe_val01 (data.get_dof_handler(0).get_fe(),
-              Quadrature<dim>(data.get_quad(1)),
+              Quadrature<dim>(data.get_quadrature(1)),
               update_values | update_gradients | update_JxW_values),
     fe_val1 (data.get_dof_handler(1).get_fe(),
-             Quadrature<dim>(data.get_quad(1)),
+             Quadrature<dim>(data.get_quadrature(1)),
              update_values | update_gradients | update_JxW_values)
   {};
 
@@ -93,10 +93,10 @@ operator () (const MatrixFree<dim,Number> &data,
   const unsigned int n_q_points1 = fe_eval1.n_q_points;
   const unsigned int dofs_per_cell0 = fe_eval0.dofs_per_cell;
   const unsigned int dofs_per_cell1 = fe_eval1.dofs_per_cell;
-  AlignedVector<vector_t> values0 (n_q_points0);
-  AlignedVector<vector_t> gradients0 (dim*n_q_points0);
-  AlignedVector<vector_t> values1 (n_q_points1);
-  AlignedVector<vector_t> gradients1 (dim*n_q_points1);
+  AlignedVector<VectorizedArray<Number> > values0 (n_q_points0);
+  AlignedVector<VectorizedArray<Number> > gradients0 (dim*n_q_points0);
+  AlignedVector<VectorizedArray<Number> > values1 (n_q_points1);
+  AlignedVector<VectorizedArray<Number> > gradients1 (dim*n_q_points1);
   std::vector<unsigned int> dof_indices0 (dofs_per_cell0);
   std::vector<unsigned int> dof_indices1 (dofs_per_cell1);
   for(unsigned int cell=cell_range.first;cell<cell_range.second;++cell)
@@ -178,7 +178,7 @@ operator () (const MatrixFree<dim,Number> &data,
       for (unsigned int q=0; q<n_q_points0; ++q)
         {
           fe_eval0.submit_value (values0[q], q);
-          Tensor<1,dim,vector_t> submit (false);
+          Tensor<1,dim,VectorizedArray<Number> > submit (false);
           for (unsigned int d=0; d<dim; ++d)
             submit[d] = gradients0[q*dim+d];
           fe_eval0.submit_gradient (submit, q);
@@ -190,7 +190,7 @@ operator () (const MatrixFree<dim,Number> &data,
       for (unsigned int q=0; q<n_q_points1; ++q)
         {
           fe_eval1.submit_value (values1[q], q);
-          Tensor<1,dim,vector_t> submit (false);
+          Tensor<1,dim,VectorizedArray<Number> > submit (false);
           for (unsigned int d=0; d<dim; ++d)
             submit[d] = gradients1[q*dim+d];
           fe_eval1.submit_gradient (submit, q);
@@ -202,7 +202,7 @@ operator () (const MatrixFree<dim,Number> &data,
       for (unsigned int q=0; q<n_q_points1; ++q)
         {
           fe_eval01.submit_value (values1[q], q);
-          Tensor<1,dim,vector_t> submit (false);
+          Tensor<1,dim,VectorizedArray<Number> > submit (false);
           for (unsigned int d=0; d<dim; ++d)
             submit[d] = gradients1[q*dim+d];
           fe_eval01.submit_gradient (submit, q);
@@ -220,30 +220,25 @@ void test ()
                                 // create hyper ball geometry and refine some
                                 // cells
   Triangulation<dim> tria;
-  GridGenerator::hyper_ball (tria);
-  static const HyperBallBoundary<dim> boundary;
-  tria.set_boundary (0, boundary);
-  typename Triangulation<dim>::active_cell_iterator
-    cell = tria.begin_active (),
-    endc = tria.end();
-  for (; cell!=endc; ++cell)
-    if (cell->center().norm()<1e-8)
-      cell->set_refine_flag();
+  create_mesh (tria);
+  tria.begin_active ()->set_refine_flag();
   tria.execute_coarsening_and_refinement();
+  typename Triangulation<dim>::active_cell_iterator cell, endc;
   cell = tria.begin_active ();
+  endc = tria.end();
   for (; cell!=endc; ++cell)
-    if (cell->center().norm()<0.2)
+    if (cell->center().norm()<0.5)
       cell->set_refine_flag();
   tria.execute_coarsening_and_refinement();
-  if (dim < 3 || fe_degree < 2)
-    tria.refine_global(1);
   tria.begin(tria.n_levels()-1)->set_refine_flag();
   tria.last()->set_refine_flag();
   tria.execute_coarsening_and_refinement();
+  tria.refine_global(1);
   cell = tria.begin_active ();
-  for (unsigned int i=0; i<7-2*dim; ++i)
+  for (unsigned int i=0; i<10-3*dim; ++i)
     {
       cell = tria.begin_active ();
+      endc = tria.end();
       unsigned int counter = 0;
       for (; cell!=endc; ++cell, ++counter)
         if (counter % (7-i) == 0)
@@ -324,7 +319,7 @@ int main ()
   deallog << std::setprecision (3);
 
   {
-    deallog.threshold_double(1.e-12);
+    deallog.threshold_double(1.e-11);
     deallog.push("2d");
     test<2,1,double>();
     test<2,2,double>();
index fc0d926027e4b9daa29f55a622ec9482389456f1..1708f8725580d36108f6fc1076d18b3c5eb898b0 100644 (file)
@@ -22,6 +22,8 @@ std::ofstream logfile("matrix_vector_06/output");
 template <int dim, int fe_degree>
 void test ()
 {
+  if (fe_degree > 1)
+    return;
   Triangulation<dim> tria;
   create_mesh (tria);
   tria.begin_active ()->set_refine_flag();
@@ -33,8 +35,6 @@ void test ()
     if (cell->center().norm()<0.5)
       cell->set_refine_flag();
   tria.execute_coarsening_and_refinement();
-  if (dim < 3 || fe_degree < 2)
-    tria.refine_global(1);
   tria.begin(tria.n_levels()-1)->set_refine_flag();
   tria.last()->set_refine_flag();
   tria.execute_coarsening_and_refinement();
index c8d1cd1e9ff2e786ac141702eb0d510f5059127d..e9e176a4c1a834620441fbe7738f5b598c0455a6 100644 (file)
@@ -2,12 +2,6 @@
 DEAL:2d::Testing FE_Q<2>(1)
 DEAL:2d::Norm of difference: 0
 DEAL:2d::
-DEAL:2d::Testing FE_Q<2>(2)
-DEAL:2d::Norm of difference: 0
-DEAL:2d::
 DEAL:3d::Testing FE_Q<3>(1)
 DEAL:3d::Norm of difference: 0
 DEAL:3d::
-DEAL:3d::Testing FE_Q<3>(2)
-DEAL:3d::Norm of difference: 0
-DEAL:3d::
index 3511799f8748ce73b89856a3da1093250881acff..1f468c9af8e58b85be65553a81895d458c1e154b 100644 (file)
@@ -20,6 +20,9 @@ std::ofstream logfile("matrix_vector_08/output");
 template <int dim, int fe_degree>
 void test ()
 {
+  if (fe_degree > 1)
+    return;
+
   Triangulation<dim> tria;
   create_mesh (tria, 1e20);
   tria.begin_active ()->set_refine_flag();
@@ -31,8 +34,6 @@ void test ()
     if (cell->center().norm()<0.5*1e20)
       cell->set_refine_flag();
   tria.execute_coarsening_and_refinement();
-  if (dim < 3 || fe_degree < 2)
-    tria.refine_global(1);
   tria.begin(tria.n_levels()-1)->set_refine_flag();
   tria.last()->set_refine_flag();
   tria.execute_coarsening_and_refinement();
index c8d1cd1e9ff2e786ac141702eb0d510f5059127d..e9e176a4c1a834620441fbe7738f5b598c0455a6 100644 (file)
@@ -2,12 +2,6 @@
 DEAL:2d::Testing FE_Q<2>(1)
 DEAL:2d::Norm of difference: 0
 DEAL:2d::
-DEAL:2d::Testing FE_Q<2>(2)
-DEAL:2d::Norm of difference: 0
-DEAL:2d::
 DEAL:3d::Testing FE_Q<3>(1)
 DEAL:3d::Norm of difference: 0
 DEAL:3d::
-DEAL:3d::Testing FE_Q<3>(2)
-DEAL:3d::Norm of difference: 0
-DEAL:3d::
index 20ca29bab8db3846927c090f5ff736f748852cf6..556f6940fe97adaf7f6cfb8a7fa3f97520b13edc 100644 (file)
@@ -20,6 +20,9 @@ std::ofstream logfile("matrix_vector_09/output");
 template <int dim, int fe_degree>
 void test ()
 {
+  if (fe_degree > 1)
+    return;
+
   Triangulation<dim> tria;
   create_mesh (tria, 1e-20);
   tria.begin_active ()->set_refine_flag();
@@ -31,12 +34,9 @@ void test ()
     if (cell->center().norm()<0.5*1e-20)
       cell->set_refine_flag();
   tria.execute_coarsening_and_refinement();
-  if (dim < 3 || fe_degree < 2)
-    tria.refine_global(1);
   tria.begin(tria.n_levels()-1)->set_refine_flag();
   tria.last()->set_refine_flag();
   tria.execute_coarsening_and_refinement();
-  tria.refine_global(1);
   cell = tria.begin_active ();
   for (unsigned int i=0; i<10-3*dim; ++i)
     {
index c8d1cd1e9ff2e786ac141702eb0d510f5059127d..e9e176a4c1a834620441fbe7738f5b598c0455a6 100644 (file)
@@ -2,12 +2,6 @@
 DEAL:2d::Testing FE_Q<2>(1)
 DEAL:2d::Norm of difference: 0
 DEAL:2d::
-DEAL:2d::Testing FE_Q<2>(2)
-DEAL:2d::Norm of difference: 0
-DEAL:2d::
 DEAL:3d::Testing FE_Q<3>(1)
 DEAL:3d::Norm of difference: 0
 DEAL:3d::
-DEAL:3d::Testing FE_Q<3>(2)
-DEAL:3d::Norm of difference: 0
-DEAL:3d::
index 0577fe6a374e9cf52965482c1c0c2e19dd58af35..fa7e807efa6a257ee5ff413473b76183cdc3c1d0 100644 (file)
@@ -25,9 +25,6 @@ template <int dim, typename Number>
 class MatrixFreeTestHP
 {
  public:
-  typedef VectorizedArray<Number> vector_t;
-  static const std::size_t n_vectors = VectorizedArray<Number>::n_array_elements;
-
   MatrixFreeTestHP(const MatrixFree<dim,Number> &data_in):
     data (data_in)
   {};
index f1caacaf72a833e2b2e88f0f0784f0da5054d78a..5c81f3effa7735bb383f1a6d569bba97d3b4ce06 100644 (file)
@@ -7,7 +7,7 @@
 
 // this function tests the correctness of the implementation of matrix free
 // matrix-vector products by comparing with the result of deal.II sparse
-// matrix. The mesh uses a hypercube mesh with no hanging nodes and no other
+// matrix. No hanging nodes and no other
 // constraints for a vector-valued problem (stokes equations).
 
 #include "../tests.h"
@@ -34,6 +34,8 @@ std::ofstream logfile("matrix_vector_stokes/output");
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/numerics/vectors.h>
 
+#include "create_mesh.h"
+
 #include <fstream>
 #include <iostream>
 #include <complex>
@@ -113,23 +115,11 @@ template <int dim, int fe_degree>
 void test ()
 {
   Triangulation<dim>   triangulation;
-  {
-    std::vector<unsigned int> subdivisions (dim, 1);
-    subdivisions[0] = 4;
-
-    const Point<dim> bottom_left = (dim == 2 ?
-                                    Point<dim>(-2,-1) :
-                                    Point<dim>(-2,0,-1));
-    const Point<dim> top_right   = (dim == 2 ?
-                                    Point<dim>(2,0) :
-                                    Point<dim>(2,1,0));
-
-    GridGenerator::subdivided_hyper_rectangle (triangulation,
-                                               subdivisions,
-                                               bottom_left,
-                                               top_right);
-  }
-  triangulation.refine_global (4-dim);
+  create_mesh (triangulation);
+  if (fe_degree == 1)
+    triangulation.refine_global (4-dim);
+  else
+    triangulation.refine_global (3-dim);
 
   FE_Q<dim>            fe_u (fe_degree+1);
   FE_Q<dim>            fe_p (fe_degree);
index e9c8909c602b29e7ee7892751f773b0caa6cef0f..6c58266926395066d07c77a1c5393b71e4322b6d 100644 (file)
@@ -341,11 +341,9 @@ int main ()
     test<2,1>();
     test<2,2>();
     test<2,3>();
-    test<2,4>();
     deallog.pop();
     deallog.push("3d");
     test<3,1>();
-    test<3,2>();
     deallog.pop();
   }
 }
index afab4527f16345c41479727bb71110633f5c7a08..0434059632a5efc3a790c9fdb06154d7c5f608ec 100644 (file)
@@ -8,9 +8,5 @@ DEAL:2d::Verification fe degree 2: 0
 DEAL:2d::
 DEAL:2d::Verification fe degree 3: 0
 DEAL:2d::
-DEAL:2d::Verification fe degree 4: 0
-DEAL:2d::
 DEAL:3d::Verification fe degree 1: 0
 DEAL:3d::
-DEAL:3d::Verification fe degree 2: 0
-DEAL:3d::
index ba030b88e9a6255b57ecf42c26e6796834056907..99c3b98af88fe8bcd4a386720594f2c73be26aeb 100644 (file)
@@ -63,9 +63,9 @@ void test ()
   }
 
   double error_points = 0, abs_points = 0;
-  const unsigned int n_cells = mf_data.get_size_info().n_macro_cells;
+  const unsigned int n_cells = mf_data.n_macro_cells();
   FEEvaluation<dim,fe_degree> fe_eval (mf_data);
-  FEValues<dim> fe_values (mapping, fe, mf_data.get_quad(),
+  FEValues<dim> fe_values (mapping, fe, mf_data.get_quadrature(),
                            update_quadrature_points);
 
   typedef VectorizedArray<double> vector_t;
index 84b58e609add44ad8783b202b24dd88b9d93e7fa..9b4940e2ab9887d731fbda4bf86282e919526f34 100644 (file)
@@ -26,9 +26,6 @@ template <int dim, typename Number>
 class MatrixFreeTestHP
 {
  public:
-  typedef VectorizedArray<Number> vector_t;
-  static const std::size_t n_vectors = VectorizedArray<Number>::n_array_elements;
-
   MatrixFreeTestHP(const MatrixFree<dim,Number> &data_in):
     data (data_in)
   {};
@@ -132,11 +129,11 @@ void do_test (const unsigned int parallel_option)
   dof.distribute_dofs(fe_collection);
   ConstraintMatrix constraints;
   DoFTools::make_hanging_node_constraints (dof,
-                                             constraints);
+                                           constraints);
   VectorTools::interpolate_boundary_values (dof,
-                                              0,
-                                              ZeroFunction<dim>(),
-                                              constraints);
+                                            0,
+                                            ZeroFunction<dim>(),
+                                            constraints);
   constraints.close ();
 
   //std::cout << "Number of cells: " << dof.get_tria().n_active_cells() << std::endl;
@@ -151,48 +148,55 @@ void do_test (const unsigned int parallel_option)
   mf_data.reinit (dof, constraints, quadrature_collection_mf, data);
   MatrixFreeTestHP<dim,number> mf (mf_data);
 
-  MatrixFree<dim,number> mf_data_par;
-  if (parallel_option == 0)
-    {
-      data.tasks_parallel_scheme =
-        MatrixFree<dim,number>::AdditionalData::partition_partition;
-      deallog << "Parallel option partition/partition" << std::endl;
-    }
-  else
+                                // test different block sizes, starting from
+                                // auto setting (= 0)
+  for (unsigned int block_size = 0; block_size < 5; ++block_size)
     {
-      data.tasks_parallel_scheme =
-        MatrixFree<dim,number>::AdditionalData::partition_color;
-      deallog << "Parallel option partition/color" << std::endl;
-    }
-  data.tasks_block_size = 1;
-  mf_data_par.reinit (dof, constraints, quadrature_collection_mf, data);
-  MatrixFreeTestHP<dim,number> mf_par(mf_data_par);
+      deallog.push ("blk_" + Utilities::int_to_string(block_size,1));
+      MatrixFree<dim,number> mf_data_par;
+      if (parallel_option == 0)
+        {
+          data.tasks_parallel_scheme =
+            MatrixFree<dim,number>::AdditionalData::partition_partition;
+          deallog << "Parallel option partition/partition" << std::endl;
+        }
+      else
+        {
+          data.tasks_parallel_scheme =
+            MatrixFree<dim,number>::AdditionalData::partition_color;
+          deallog << "Parallel option partition/color" << std::endl;
+        }
+      data.tasks_block_size = 1;
+      mf_data_par.reinit (dof, constraints, quadrature_collection_mf, data);
+      MatrixFreeTestHP<dim,number> mf_par(mf_data_par);
 
                                 // fill a right hand side vector with random
                                 // numbers in unconstrained degrees of freedom
-  Vector<number> src (dof.n_dofs());
-  Vector<number> result_ref(src), result_mf (src);
+      Vector<number> src (dof.n_dofs());
+      Vector<number> result_ref(src), result_mf (src);
 
-  for (unsigned int i=0; i<dof.n_dofs(); ++i)
-    {
-      if (constraints.is_constrained(i) == false)
-        src(i) = (double)rand()/RAND_MAX;
-    }
+      for (unsigned int i=0; i<dof.n_dofs(); ++i)
+        {
+          if (constraints.is_constrained(i) == false)
+            src(i) = (double)rand()/RAND_MAX;
+        }
 
-                                // now perform 50 matrix-vector products in
+                                // now perform 30 matrix-vector products in
                                 // parallel and check their correctness (take
                                 // many of them to make sure that we hit an
                                 // error)
-  mf.vmult (result_ref, src);
-  deallog << "Norm of difference: ";
-  for (unsigned int i=0; i<50; ++i)
-    {
-      mf_par.vmult (result_mf, src);
-      result_mf -= result_ref;
-      double diff_norm = result_mf.linfty_norm()/result_ref.linfty_norm();
-      deallog << diff_norm << "  ";
+      mf.vmult (result_ref, src);
+      deallog << "Norm of difference: ";
+      for (unsigned int i=0; i<50; ++i)
+        {
+          mf_par.vmult (result_mf, src);
+          result_mf -= result_ref;
+          double diff_norm = result_mf.linfty_norm()/result_ref.linfty_norm();
+          deallog << diff_norm << "  ";
+        }
+      deallog << std::endl << std::endl;
+      deallog.pop();
     }
-  deallog << std::endl << std::endl;
 }
 
 
index 0bfece69138edf31c9972e546d34af3a2a6dcf88..4164878d98ad6cc7016a64368b75f20a449d7f10 100644 (file)
 
-DEAL:2d:double::Parallel option partition/partition
-DEAL:2d:double::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
-DEAL:2d:double::
-DEAL:2d:float::Parallel option partition/partition
-DEAL:2d:float::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
-DEAL:2d:float::
-DEAL:2d:double::Parallel option partition/color
-DEAL:2d:double::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
-DEAL:2d:double::
-DEAL:2d:float::Parallel option partition/color
-DEAL:2d:float::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
-DEAL:2d:float::
-DEAL:3d:double::Parallel option partition/partition
-DEAL:3d:double::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
-DEAL:3d:double::
-DEAL:3d:float::Parallel option partition/partition
-DEAL:3d:float::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
-DEAL:3d:float::
-DEAL:3d:double::Parallel option partition/color
-DEAL:3d:double::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
-DEAL:3d:double::
-DEAL:3d:float::Parallel option partition/color
-DEAL:3d:float::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
-DEAL:3d:float::
+DEAL:2d:double:blk_0::Parallel option partition/partition
+DEAL:2d:double:blk_0::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:double:blk_0::
+DEAL:2d:double:blk_1::Parallel option partition/partition
+DEAL:2d:double:blk_1::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:double:blk_1::
+DEAL:2d:double:blk_2::Parallel option partition/partition
+DEAL:2d:double:blk_2::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:double:blk_2::
+DEAL:2d:double:blk_3::Parallel option partition/partition
+DEAL:2d:double:blk_3::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:double:blk_3::
+DEAL:2d:double:blk_4::Parallel option partition/partition
+DEAL:2d:double:blk_4::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:double:blk_4::
+DEAL:2d:float:blk_0::Parallel option partition/partition
+DEAL:2d:float:blk_0::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:float:blk_0::
+DEAL:2d:float:blk_1::Parallel option partition/partition
+DEAL:2d:float:blk_1::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:float:blk_1::
+DEAL:2d:float:blk_2::Parallel option partition/partition
+DEAL:2d:float:blk_2::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:float:blk_2::
+DEAL:2d:float:blk_3::Parallel option partition/partition
+DEAL:2d:float:blk_3::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:float:blk_3::
+DEAL:2d:float:blk_4::Parallel option partition/partition
+DEAL:2d:float:blk_4::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:float:blk_4::
+DEAL:2d:double:blk_0::Parallel option partition/color
+DEAL:2d:double:blk_0::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:double:blk_0::
+DEAL:2d:double:blk_1::Parallel option partition/color
+DEAL:2d:double:blk_1::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:double:blk_1::
+DEAL:2d:double:blk_2::Parallel option partition/color
+DEAL:2d:double:blk_2::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:double:blk_2::
+DEAL:2d:double:blk_3::Parallel option partition/color
+DEAL:2d:double:blk_3::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:double:blk_3::
+DEAL:2d:double:blk_4::Parallel option partition/color
+DEAL:2d:double:blk_4::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:double:blk_4::
+DEAL:2d:float:blk_0::Parallel option partition/color
+DEAL:2d:float:blk_0::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:float:blk_0::
+DEAL:2d:float:blk_1::Parallel option partition/color
+DEAL:2d:float:blk_1::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:float:blk_1::
+DEAL:2d:float:blk_2::Parallel option partition/color
+DEAL:2d:float:blk_2::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:float:blk_2::
+DEAL:2d:float:blk_3::Parallel option partition/color
+DEAL:2d:float:blk_3::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:float:blk_3::
+DEAL:2d:float:blk_4::Parallel option partition/color
+DEAL:2d:float:blk_4::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:2d:float:blk_4::
+DEAL:3d:double:blk_0::Parallel option partition/partition
+DEAL:3d:double:blk_0::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:double:blk_0::
+DEAL:3d:double:blk_1::Parallel option partition/partition
+DEAL:3d:double:blk_1::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:double:blk_1::
+DEAL:3d:double:blk_2::Parallel option partition/partition
+DEAL:3d:double:blk_2::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:double:blk_2::
+DEAL:3d:double:blk_3::Parallel option partition/partition
+DEAL:3d:double:blk_3::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:double:blk_3::
+DEAL:3d:double:blk_4::Parallel option partition/partition
+DEAL:3d:double:blk_4::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:double:blk_4::
+DEAL:3d:float:blk_0::Parallel option partition/partition
+DEAL:3d:float:blk_0::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:float:blk_0::
+DEAL:3d:float:blk_1::Parallel option partition/partition
+DEAL:3d:float:blk_1::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:float:blk_1::
+DEAL:3d:float:blk_2::Parallel option partition/partition
+DEAL:3d:float:blk_2::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:float:blk_2::
+DEAL:3d:float:blk_3::Parallel option partition/partition
+DEAL:3d:float:blk_3::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:float:blk_3::
+DEAL:3d:float:blk_4::Parallel option partition/partition
+DEAL:3d:float:blk_4::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:float:blk_4::
+DEAL:3d:double:blk_0::Parallel option partition/color
+DEAL:3d:double:blk_0::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:double:blk_0::
+DEAL:3d:double:blk_1::Parallel option partition/color
+DEAL:3d:double:blk_1::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:double:blk_1::
+DEAL:3d:double:blk_2::Parallel option partition/color
+DEAL:3d:double:blk_2::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:double:blk_2::
+DEAL:3d:double:blk_3::Parallel option partition/color
+DEAL:3d:double:blk_3::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:double:blk_3::
+DEAL:3d:double:blk_4::Parallel option partition/color
+DEAL:3d:double:blk_4::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:double:blk_4::
+DEAL:3d:float:blk_0::Parallel option partition/color
+DEAL:3d:float:blk_0::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:float:blk_0::
+DEAL:3d:float:blk_1::Parallel option partition/color
+DEAL:3d:float:blk_1::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:float:blk_1::
+DEAL:3d:float:blk_2::Parallel option partition/color
+DEAL:3d:float:blk_2::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:float:blk_2::
+DEAL:3d:float:blk_3::Parallel option partition/color
+DEAL:3d:float:blk_3::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:float:blk_3::
+DEAL:3d:float:blk_4::Parallel option partition/color
+DEAL:3d:float:blk_4::Norm of difference: 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  
+DEAL:3d:float:blk_4::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.