]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Simplify FEEvaluation routines considerably by merging all functionality into FEEvalu...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 27 Jun 2014 20:31:28 +0000 (20:31 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 27 Jun 2014 20:31:28 +0000 (20:31 +0000)
git-svn-id: https://svn.dealii.org/trunk@33094 0785d39b-7218-0410-832d-ea1e28bc413d

19 files changed:
deal.II/doc/news/changes.h
deal.II/examples/step-48/doc/intro.dox
deal.II/examples/step-48/step-48.cc
deal.II/include/deal.II/matrix_free/fe_evaluation.h
deal.II/include/deal.II/matrix_free/operators.h
deal.II/include/deal.II/matrix_free/shape_info.h
deal.II/include/deal.II/matrix_free/shape_info.templates.h
tests/matrix_free/estimate_condition_number_mass.cc
tests/matrix_free/evaluate_1d_shape.cc
tests/matrix_free/evaluate_1d_shape_evenodd.cc
tests/matrix_free/get_functions_gl.cc
tests/matrix_free/get_functions_q_hierarchical.cc
tests/matrix_free/matrix_vector_14.cc
tests/matrix_free/matrix_vector_14b.cc
tests/matrix_free/matrix_vector_17.cc
tests/matrix_free/matrix_vector_stokes_qdg0.cc
tests/matrix_free/step-48.cc
tests/matrix_free/step-48b.cc
tests/matrix_free/step-48c.cc

index c22b5985375b7e6cc08017c8356ec5d53ed0fb5f..ba57ec5880d9cb2d750f2fefe8d414b91fe9d871 100644 (file)
@@ -148,6 +148,20 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Simplified interfaces for FEEvaluation: Previously, the user had to
+  select the appropriate kernel (FEEvaluation, FEEvaluationGeneral,
+  FEEvaluationDGP, FEEvaluationGL) for the matrix-free evaluation
+  routines. This made it difficult to write compact code that needs to select
+  between different elements. Therefore, all the functionality has been merged
+  into FEEvaluation and one should only use FEEvaluation, while the other
+  interfaces have been marked deprecated and will be removed in a future
+  version. The internal data structures for the various special cases have
+  been kept in order to provide for the most efficient routines, which is
+  selected at construction of FEEvaluation.
+  <br>
+  (Martin Kronbichler, 2014/06/27)
+  </li>
+  
   <li> Bugfix: TBB specific include directories have to be added to the
   list of user include directories because of direct inclusion of header
   files in base/thread_local_storage.h.
@@ -180,8 +194,8 @@ inconvenience this causes.
   (Denis Davydov, 2014/06/19)
   </li>
 
-  <li> New: There is now a class FEEvaluationQ_DG0 that does
-  optimized matrix-free evaluation for FE_Q_DG0 elements.
+  <li> New: The class FEEvaluation now provides optimized matrix-free
+  evaluation routines for FE_Q_DG0 elements.
   <br>
   (Martin Kronbichler, 2014/06/13)
   </li>
index 6a89c11ca8ef7289a2c17ca0dac4718db7fc7041..a4715e6f80cf123265e0c4964356b9b20214947b 100644 (file)
@@ -127,7 +127,7 @@ variable <code>constraint_indicator</code> at hand that holds, for
 each cell, the local indices of DoFs that are constrained as well as
 the identifier of the type of constraint. Actually, you will not see
 these data structures in the example program since the class
-<code>FEEvaluationGL</code> takes care of the constraints without user
+<code>FEEvaluation</code> takes care of the constraints without user
 interaction.
 
 
index 77e1f719703b4bdb3b93cec89b6ac6d8e5ac8952..1c2843423849389a1815e07afd9deba4e66232d7 100644 (file)
@@ -125,7 +125,7 @@ namespace Step48
 
     data.initialize_dof_vector (inv_mass_matrix);
 
-    FEEvaluationGL<dim,fe_degree> fe_eval(data);
+    FEEvaluation<dim,fe_degree>   fe_eval(data);
     const unsigned int            n_q_points = fe_eval.n_q_points;
 
     for (unsigned int cell=0; cell<data.n_macro_cells(); ++cell)
@@ -151,9 +151,14 @@ namespace Step48
 
   // This operator implements the core operation of the program, the
   // integration over a range of cells for the nonlinear operator of the
-  // Sine-Gordon problem. The implementation is based on the FEEvaluationGL
-  // class since we are using the cell-based implementation for Gauss-Lobatto
-  // elements.
+  // Sine-Gordon problem. The implementation is based on the FEEvaluation
+  // class as in step-37. Due to the special structure in Gauss-Lobatto
+  // elements, certain operations become simpler, in particular the evaluation
+  // of shape function values on quadrature points which is simply the
+  // injection of the values of cell degrees of freedom. The MatrixFree class
+  // detects possible structure of the finite element at quadrature points
+  // when initializing, which is then used by FEEvaluation for selecting the
+  // most appropriate numerical kernel.
 
   // The nonlinear function that we have to evaluate for the time stepping
   // routine includes the value of the function at the present time @p current
@@ -183,7 +188,7 @@ namespace Step48
                const std::pair<unsigned int,unsigned int> &cell_range) const
   {
     AssertDimension (src.size(), 2);
-    FEEvaluationGL<dim,fe_degree> current (data), old (data);
+    FEEvaluation<dim,fe_degree> current (data), old (data);
     for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
       {
         current.reinit (cell);
index c00b23b0b95d8b90fdb983a934031f2324be2c66..58df9f018e0f443475fe154ecefe203e63c5d4d1 100644 (file)
@@ -53,13 +53,12 @@ namespace internal
 /**
  * This is the base class for the FEEvaluation classes. This class is a base
  * class and needs usually not be called in user code. It does not have any
- * public constructor. Use one of the derived classes FEEvaluationGeneral,
- * FEEvaluation, FEEvaluationGL, or FEEvaluationDGP instead. It implements a
- * reinit method that is used to set pointers so that operations on quadrature
- * points can be performed quickly, access functions to vectors for the @p
- * read_dof_values, @p set_dof_values, and @p distributed_local_to_global
- * functions, as well as methods to access values and gradients of finite
- * element functions.
+ * public constructor. The usage is through the class FEEvaluation instead. It
+ * implements a reinit method that is used to set pointers so that operations
+ * on quadrature points can be performed quickly, access functions to vectors
+ * for the @p read_dof_values, @p set_dof_values, and @p
+ * distributed_local_to_global functions, as well as methods to access values
+ * and gradients of finite element functions.
  *
  * This class has three template arguments:
  *
@@ -555,7 +554,7 @@ protected:
   FEEvaluationBase (const MatrixFree<dim,Number> &matrix_free,
                     const unsigned int            fe_no,
                     const unsigned int            quad_no,
-                    const unsigned int            dofs_per_cell,
+                    const unsigned int            fe_degree,
                     const unsigned int            n_q_points);
 
   /**
@@ -863,7 +862,7 @@ protected:
   FEEvaluationAccess (const MatrixFree<dim,Number> &matrix_free,
                       const unsigned int            fe_no,
                       const unsigned int            quad_no,
-                      const unsigned int            dofs_per_cell,
+                      const unsigned int            fe_degree,
                       const unsigned int            n_q_points);
 
   /**
@@ -1020,7 +1019,7 @@ protected:
   FEEvaluationAccess (const MatrixFree<dim,Number> &matrix_free,
                       const unsigned int            fe_no,
                       const unsigned int            quad_no,
-                      const unsigned int            dofs_per_cell,
+                      const unsigned int            fe_degree,
                       const unsigned int            n_q_points);
 
   /**
@@ -1227,196 +1226,6 @@ protected:
 
 
 
-/**
- * The class that provides all functions necessary to evaluate functions at
- * quadrature points and cell integrations. In functionality, this class is
- * similar to FEValues<dim>, however, it includes a lot of specialized
- * functions that make it much faster (between 5 and 500 times as fast,
- * depending on the polynomial order). Access to the data fields is provided
- * through functionality in the class FEEvaluationAccess.
- *
- * This class is designed for general local finite element operations based on
- * tensor products of 1D polynomials and quadrature points. Often, there are
- * some symmetries or zeros in the unit cell data that allow for a more
- * efficient operator application. FEEvaluation is specialized to standard
- * FE_Q/FE_DGQ elements and quadrature points symmetric around 0.5 (like Gauss
- * quadrature), and hence the most common situation. FEEvaluationGL is a
- * specialization for elements where quadrature formula and support points are
- * chosen so that a orthogonal relation between the values holds. This is the
- * case for FE_Q based on Gauss-Lobatto support points with Gauss-Lobatto
- * quadrature formula of the same order.
- *
- * This class has five template arguments:
- *
- * @param dim Dimension in which this class is to be used
- *
- * @param fe_degree Degree of the tensor product finite element with
- *                  fe_degree+1 degrees of freedom per coordinate direction
- *
- * @param n_q_points_1d Number of points in the quadrature formula in 1D,
- *                  usually chosen as fe_degree+1
- *
- * @param n_components Number of vector components when solving a system of
- *                  PDEs. If the same operation is applied to several
- *                  components of a PDE (e.g. a vector Laplace equation), they
- *                  can be applied simultaneously with one call (and often
- *                  more efficiently)
- *
- * @param Number Number format, usually @p double or @p float
- *
- * @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 >
-class FEEvaluationGeneral : public FEEvaluationAccess<dim, n_components_,Number>
-{
-public:
-  typedef FEEvaluationAccess<dim, 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 = Utilities::fixed_int_power<fe_degree+1,dim>::value;
-  static const unsigned int n_q_points    = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
-
-  /**
-   * Constructor. Takes all data stored in MatrixFree. If applied to problems
-   * with more than one finite element or more than one quadrature formula
-   * selected during construction of @p matrix_free, @p fe_no and @p quad_no
-   * allow to select the appropriate components.
-   */
-  FEEvaluationGeneral (const MatrixFree<dim,Number> &matrix_free,
-                       const unsigned int            fe_no   = 0,
-                       const unsigned int            quad_no = 0);
-
-  /**
-   * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type MappingFEEvaluation
-   * and a DoFHandler in order to allow for reading out the finite element
-   * data. It uses the data provided by dof_handler.get_fe(). If the element
-   * is vector-valued, the optional argument allows to specify the index of
-   * the base element (as long as the element is primitive, non-primitive are
-   * not supported currently).
-   *
-   * With initialization from a FEValues object, no call to a reinit method of
-   * this class is necessary. Instead, it is enough if the geometry is
-   * initialized to a given cell iterator. It can also read from or write to
-   * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
-   * types (which is less efficient with MPI since index translation has to be
-   * done), but of course only for one cell at a time. Hence, a kernel using
-   * this method does not vectorize over several elements (which is most
-   * efficient for vector operations), but only possibly within the element if
-   * the evaluate/integrate routines are combined (e.g. for computing cell
-   * matrices).
-   * With this initialization, no call to a reinit method of this
-   * class. Instead, it is enough if the geometry is initialized to a given
-   * cell iterator. Moreover, beware that a kernel using this method does not
-   * vectorize over several elements (which is most efficient for vector
-   * operations), but only possibly within the element if the
-   * evaluate/integrate routines are combined (e.g. for matrix assembly).
-   */
-  FEEvaluationGeneral (const MappingFEEvaluation<dim,Number> &geometry,
-                       const DoFHandler<dim>                 &dof_handler,
-                       const unsigned int                     first_selected_component = 0);
-
-  /**
-   * Copy constructor
-   */
-  FEEvaluationGeneral (const FEEvaluationGeneral &other);
-
-  /**
-   * Evaluates the function values, the gradients, and the Laplacians of the
-   * FE function given at the DoF values in the input vector at the quadrature
-   * points.  The function arguments specify which parts shall actually be
-   * computed. Needs to be called before the functions @p get_value(), @p
-   * get_gradient() or @p get_laplacian return useful information.
-   */
-  void evaluate (const bool evaluate_val,
-                 const bool evaluate_grad,
-                 const bool evaluate_hess = false);
-
-  /**
-   * This function takes the values and/or gradients that are stored on
-   * quadrature points, tests them by all the basis functions/gradients on the
-   * cell and performs the cell integration. The two function arguments @p
-   * integrate_val and @p integrate_grad are used to enable/disable some of
-   * values or gradients.
-   */
-  void integrate (const bool integrate_val,
-                  const bool integrate_grad);
-
-  /**
-   * Returns the q-th quadrature point stored in MappingInfo.
-   */
-  Point<dim,VectorizedArray<Number> >
-  quadrature_point (const unsigned int q_point) const;
-
-  /**
-   * 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 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_hessians (const VectorizedArray<Number> in [],
-                       VectorizedArray<Number> out []);
-
-protected:
-
-  /**
-   * Internally stored variables for the different data fields.
-   */
-  VectorizedArray<Number> my_data_array[n_components*(dofs_per_cell+1+(dim*dim+2*dim+1)*n_q_points)];
-
-  /**
-   * Checks the number of cell dofs with the value that we expect
-   */
-  void check_dofs_per_cell (const unsigned int dofs_per_cell) const;
-
-private:
-  /**
-   * Sets the pointers from the data array to values_dof, etc.
-   */
-  void set_data_pointers();
-
-  /**
-   * Checks that the element number of degrees of freedom given by the
-   * template arguments (via fe_degree) coincides with the number of degrees
-   * of freedom in the stored shape values.
-   */
-  void check_template_arguments(const unsigned int fe_no);
-};
-
-
-
 /**
  * The class that provides all functions necessary to evaluate functions at
  * quadrature points and cell integrations. In functionality, this class is
@@ -1434,14 +1243,16 @@ private:
  * possibly be vectorized by combining several operations (e.g. when
  * performing matrix assembly).
  *
- * This class is a specialization of FEEvaluationGeneral designed for standard
- * FE_Q or FE_DGQ elements and quadrature points symmetric around 0.5 (like
- * Gauss quadrature), and hence the most common situation. Note that many of
- * the operations available through this class are inherited from the base
- * class FEEvaluationBase, in particular reading from and writing to
- * vectors. Also, the class inherits from FEEvaluationAccess that implements
- * access to values, gradients and Hessians of the finite element function on
- * quadrature points.
+ * This class contains specialized evaluation routines for several elements,
+ * including standard FE_Q or FE_DGQ elements and quadrature points symmetric
+ * around 0.5 (like Gauss quadrature), FE_DGP elements based on truncated
+ * tensor products as well as the faster case of Gauss-Lobatto elements with
+ * Gauss-Lobatto quadrature which give diagonal mass matrices and quicker
+ * evaluation internally. Note that many of the operations available through
+ * this class are inherited from the base class FEEvaluationBase, in
+ * particular reading from and writing to vectors. Also, the class inherits
+ * from FEEvaluationAccess that implements access to values, gradients and
+ * Hessians of the finite element function on quadrature points.
  *
  * This class assumes that shape functions of the FiniteElement under
  * consideration do <em>not</em> depend on the actual shape of the cells in
@@ -1470,18 +1281,17 @@ private:
  */
 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>
+class FEEvaluation : public FEEvaluationAccess<dim,n_components_,Number>
 {
 public:
-  typedef FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number> BaseClass;
+  typedef FEEvaluationAccess<dim,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;
+  static const unsigned int n_q_points    = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
+  static const unsigned int tensor_dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
 
   /**
    * Constructor. Takes all data stored in MatrixFree. If applied to problems
@@ -1545,92 +1355,104 @@ public:
                   const bool integrate_grad);
 
   /**
-   * 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.
+   * Returns the q-th quadrature point stored in MappingInfo.
    */
-  template <int direction, bool dof_to_quad, bool add>
-  void apply_values (const VectorizedArray<Number> in [],
-                     VectorizedArray<Number> out []);
+  Point<dim,VectorizedArray<Number> >
+  quadrature_point (const unsigned int q_point) const;
 
   /**
-   * 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 previous content in the data fields or not.
+   * The number of scalar degrees of freedom on the cell. Usually close to
+   * tensor_dofs_per_cell, but depends on the actual element selected and is
+   * thus not static.
    */
-  template <int direction, bool dof_to_quad, bool add>
-  void apply_gradients (const VectorizedArray<Number> in [],
-                        VectorizedArray<Number> out []);
+  const unsigned int dofs_per_cell;
 
+private:
   /**
-   * 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.
+   * Internally stored variables for the different data fields.
    */
-  template <int direction, bool dof_to_quad, bool add>
-  void apply_hessians (const VectorizedArray<Number> in [],
-                       VectorizedArray<Number> out []);
+  VectorizedArray<Number> my_data_array[n_components*(tensor_dofs_per_cell+1+(dim*dim+2*dim+1)*n_q_points)];
 
-protected:
-  VectorizedArray<Number> shape_val_evenodd[fe_degree+1][(n_q_points_1d+1)/2];
-  VectorizedArray<Number> shape_gra_evenodd[fe_degree+1][(n_q_points_1d+1)/2];
-  VectorizedArray<Number> shape_hes_evenodd[fe_degree+1][(n_q_points_1d+1)/2];
+  /**
+   * Checks if the template arguments regarding degree of the element
+   * corresponds to the actual element used at initialization.
+   */
+  void check_template_arguments(const unsigned int fe_no);
+
+  /**
+   * Sets the pointers of the base class to my_data_array.
+   */
+  void set_data_pointers();
+
+  /**
+   * Function pointer for the evaluate function
+   */
+  void (*evaluate_funct) (const internal::MatrixFreeFunctions::ShapeInfo<Number> &,
+                          VectorizedArray<Number> *values_dofs_actual[],
+                          VectorizedArray<Number> *values_quad[],
+                          VectorizedArray<Number> *gradients_quad[][dim],
+                          VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+                          const bool               evaluate_val,
+                          const bool               evaluate_grad,
+                          const bool               evaluate_lapl);
 
-private:
   /**
-   * Fills the fields shape_???_evenodd, called in the constructor.
+   * Function pointer for the integrate function
    */
-  void compute_even_odd_factors();
+  void (*integrate_funct)(const internal::MatrixFreeFunctions::ShapeInfo<Number> &,
+                          VectorizedArray<Number> *values_dofs_actual[],
+                          VectorizedArray<Number> *values_quad[],
+                          VectorizedArray<Number> *gradients_quad[][dim],
+                          const bool               evaluate_val,
+                          const bool               evaluate_grad);
 };
 
 
 
 /**
- * The class that provides all functions necessary to evaluate functions at
- * quadrature points and cell integrations. In functionality, this class is
- * similar to FEValues<dim>, however, it includes a lot of specialized
- * functions that make it much faster (between 5 and 500, depending on the
- * polynomial order).
- *
- * This class is a specialization of FEEvaluation for elements where
- * quadrature formula and support points are chosen so that a orthonormal
- * relation between the values holds. This is the case for FE_Q based on
- * Gauss-Lobatto support points with Gauss-Lobatto quadrature formula of the
- * same order (QGaussLobatto). In that case, application of values is trivial
- * (as the transformation is the identity matrix), and application of
- * gradients is considerably simpler (since all value applications in
- * directions other than the gradient direction are again identity
- * operations).
- *
- * This class has four template arguments:
- *
- * @param dim Dimension in which this class is to be used
- *
- * @param fe_degree Degree of the tensor product finite element with
- *                  fe_degree+1 degrees of freedom per coordinate
- *                  direction. The quadrature formula is tied to the choice of
- *                  the element by setting n_q_points_1d = fe_degree+1, which
- *                  gives a diagonal mass matrix
- *
- * @param n_components Number of vector components when solving a system of
- *                  PDEs. If the same operation is applied to several
- *                  components of a PDE (e.g. a vector Laplace equation), they
- *                  can be applied simultaneously with one call (and often
- *                  more efficiently)
- *
- * @param Number Number format, usually @p double or @p float
- *
- * @author Katharina Kormann and Martin Kronbichler, 2010, 2011
+ * Deprecated. Functionality has been merged into FEEvaluation. Use class
+ * FEEvaluation instead.
+ */
+template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
+          int n_components_ = 1, typename Number = double >
+class FEEvaluationGeneral : public FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_,Number>
+{
+public:
+  typedef FEEvaluation<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 = Utilities::fixed_int_power<fe_degree+1,dim>::value;
+  static const unsigned int n_q_points    = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
+
+  /**
+   * Constructor.
+   */
+  FEEvaluationGeneral (const MatrixFree<dim,Number> &matrix_free,
+                       const unsigned int            fe_no   = 0,
+                       const unsigned int            quad_no = 0) DEAL_II_DEPRECATED
+    :
+    BaseClass (matrix_free, fe_no, quad_no)
+  {}
+
+  /**
+   * Constructor.
+   */
+  FEEvaluationGeneral (const MappingFEEvaluation<dim,Number> &geometry,
+                       const DoFHandler<dim>                 &dof_handler,
+                       const unsigned int                     first_selected_component = 0) DEAL_II_DEPRECATED
+    :
+    BaseClass (geometry, dof_handler, first_selected_component)
+  {}
+};
+
+
+
+/**
+ * Deprecated. Functionality has been merged into FEEvaluation. Use class
+ * FEEvaluation instead.
  */
 template <int dim, int fe_degree, int n_components_ = 1, typename Number = double >
 class FEEvaluationGL :
@@ -1643,81 +1465,28 @@ public:
   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 dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
   static const unsigned int n_q_points    = BaseClass::n_q_points;
 
   /**
-   * Constructor. Takes all data stored in MatrixFree. If applied to problems
-   * with more than one finite element or more than one quadrature formula
-   * selected during construction of @p matrix_free, @p fe_no and @p quad_no
-   * allow to select the appropriate components.
+   * Constructor.
    */
   FEEvaluationGL (const MatrixFree<dim,Number> &matrix_free,
                   const unsigned int          fe_no   = 0,
-                  const unsigned int          quad_no = 0);
+                  const unsigned int          quad_no = 0) DEAL_II_DEPRECATED
+    :
+    BaseClass (matrix_free, fe_no, quad_no)
+  {}
 
   /**
-   * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type MappingFEEvaluation
-   * and a DoFHandler in order to allow for reading out the finite element
-   * data. It uses the data provided by dof_handler.get_fe(). If the element
-   * is vector-valued, the optional argument allows to specify the index of
-   * the base element (as long as the element is primitive, non-primitive are
-   * not supported currently).
-   *
-   * With initialization from a FEValues object, no call to a reinit method of
-   * this class is necessary. Instead, it is enough if the geometry is
-   * initialized to a given cell iterator. It can also read from or write to
-   * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
-   * types (which is less efficient with MPI since index translation has to be
-   * done), but of course only for one cell at a time. Hence, a kernel using
-   * this method does not vectorize over several elements (which is most
-   * efficient for vector operations), but only possibly within the element if
-   * the evaluate/integrate routines are combined (e.g. for computing cell
-   * matrices).
+   * Constructor.
    */
   FEEvaluationGL (const MappingFEEvaluation<dim,Number> &geometry,
                   const DoFHandler<dim>                 &dof_handler,
-                  const unsigned int                     first_selected_component = 0);
-
-  /**
-   * Copy constructor
-   */
-  FEEvaluationGL (const FEEvaluationGL &other);
-
-  /**
-   * Evaluates the function values, the gradients, and the Hessians of the FE
-   * function given at the DoF values in the input vector at the quadrature
-   * points of the unit cell. The function arguments specify which parts shall
-   * actually be computed. Needs to be called before the functions @p
-   * get_value(), @p get_gradient() or @p get_laplacian give useful
-   * information (unless these values have been set manually).
-   */
-  void evaluate (const bool evaluate_val,
-                 const bool evaluate_grad,
-                 const bool evaluate_lapl = false);
-
-  /**
-   * This function takes the values and/or gradients that are stored on
-   * quadrature points, tests them by all the basis functions/gradients on the
-   * cell and performs the cell integration. The two function arguments @p
-   * integrate_val and @p integrate_grad are used to enable/disable some of
-   * values or gradients.
-   */
-  void integrate (const bool integrate_val,
-                  const bool integrate_grad);
-
-  /**
-   * 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.
-   */
-  template <int direction, bool dof_to_quad, bool add>
-  void apply_gradients (const VectorizedArray<Number> in [],
-                        VectorizedArray<Number> out []);
+                  const unsigned int                     first_selected_component = 0) DEAL_II_DEPRECATED
+    :
+    BaseClass (geometry, dof_handler, first_selected_component)
+  {}
 };
 
 
@@ -1748,27 +1517,16 @@ namespace internal
 
 
 /**
- * The class that provides all functions necessary to evaluate functions at
- * quadrature points and cell integrations. In functionality, this class is
- * similar to FEValues<dim>, however, it includes a lot of specialized
- * functions that make it much faster (between 5 and 500 times as fast,
- * depending on the polynomial order). Access to the data fields is provided
- * through functionality in the class FEEvaluationAccess.
- *
- * This class is an extension of FEEvaluationGeneral to work with elements of
- * complete polynomial degree p, FE_DGP. In this case, the polynomial basis is
- * a truncated tensor product, so the evaluate and integrate routines use a
- * truncation.
- *
- * @author Martin Kronbichler, 2014
+ * Functionality has been merged into FEEvaluation. Use class FEEvaluation
+ * instead.
  */
 template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
           int n_components_ = 1, typename Number = double >
 class FEEvaluationDGP :
-  public FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+  public FEEvaluation<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 FEEvaluation<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;
@@ -1778,160 +1536,28 @@ public:
   static const unsigned int n_q_points    = BaseClass::n_q_points;
 
   /**
-   * Constructor. Takes all data stored in MatrixFree. If applied to problems
-   * with more than one finite element or more than one quadrature formula
-   * selected during construction of @p matrix_free, @p fe_no and @p quad_no
-   * allow to select the appropriate components.
+   * Constructor.
    */
   FEEvaluationDGP (const MatrixFree<dim,Number> &matrix_free,
                    const unsigned int            fe_no   = 0,
-                   const unsigned int            quad_no = 0);
+                   const unsigned int            quad_no = 0) DEAL_II_DEPRECATED
+    :
+    BaseClass (matrix_free, fe_no, quad_no)
+  {}
 
   /**
-   * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type MappingFEEvaluation
-   * and a DoFHandler in order to allow for reading out the finite element
-   * data. It uses the data provided by dof_handler.get_fe(). If the element
-   * is vector-valued, the optional argument allows to specify the index of
-   * the base element (as long as the element is primitive, non-primitive are
-   * not supported currently).
-   *
-   * With initialization from a FEValues object, no call to a reinit method of
-   * this class is necessary. Instead, it is enough if the geometry is
-   * initialized to a given cell iterator. It can also read from or write to
-   * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
-   * types (which is less efficient with MPI since index translation has to be
-   * done), but of course only for one cell at a time. Hence, a kernel using
-   * this method does not vectorize over several elements (which is most
-   * efficient for vector operations), but only possibly within the element if
-   * the evaluate/integrate routines are combined (e.g. for computing cell
-   * matrices).
+   * Constructor.
    */
   FEEvaluationDGP (const MappingFEEvaluation<dim,Number> &geometry,
                    const DoFHandler<dim>                 &dof_handler,
-                   const unsigned int                     first_selected_component = 0);
-
-  /**
-   * Copy constructor
-   */
-  FEEvaluationDGP (const FEEvaluationDGP &other);
-
-  /**
-   * Evaluates the function values, the gradients, and the Hessians of the FE
-   * function given at the DoF values in the input vector at the quadrature
-   * points of the unit cell. The function arguments specify which parts shall
-   * actually be computed. Needs to be called before the functions @p
-   * get_value(), @p get_gradient() or @p get_laplacian give useful
-   * information (unless these values have been set manually).
-   */
-  void evaluate (const bool evaluate_val,
-                 const bool evaluate_grad,
-                 const bool evaluate_lapl = false);
-
-  /**
-   * This function takes the values and/or gradients that are stored on
-   * quadrature points, tests them by all the basis functions/gradients on the
-   * cell and performs the cell integration. The two function arguments @p
-   * integrate_val and @p integrate_grad are used to enable/disable some of
-   * values or gradients.
-   */
-  void integrate (const bool integrate_val,
-                  const bool integrate_grad);
+                   const unsigned int                     first_selected_component = 0) DEAL_II_DEPRECATED
+    :
+    BaseClass (geometry, dof_handler, first_selected_component)
+  {}
 };
 
 
 
-/**
- * The class that provides all functions necessary to evaluate functions at
- * quadrature points and cell integrations. In functionality, this class is
- * similar to FEValues<dim>, however, it includes a lot of specialized
- * functions that make it much faster (between 5 and 500 times as fast,
- * depending on the polynomial order). Access to the data fields is provided
- * through functionality in the class FEEvaluationAccess.
- *
- * This class is an extension of FEEvaluation to work with continuous elements
- * supplemented with a single discontinuous degree of freedom, i.e., FE_Q_DG0.
- *
- * @author Martin Kronbichler, 2014
- */
-template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
-          int n_components_ = 1, typename Number = double >
-class FEEvaluationQ_DG0 :
-  public FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
-{
-public:
-  typedef FEEvaluation<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 + 1;
-  static const unsigned int n_q_points    = BaseClass::n_q_points;
-
-  /**
-   * Constructor. Takes all data stored in MatrixFree. If applied to problems
-   * with more than one finite element or more than one quadrature formula
-   * selected during construction of @p matrix_free, @p fe_no and @p quad_no
-   * allow to select the appropriate components.
-   */
-  FEEvaluationQ_DG0 (const MatrixFree<dim,Number> &matrix_free,
-                     const unsigned int            fe_no   = 0,
-                     const unsigned int            quad_no = 0);
-
-  /**
-   * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type MappingFEEvaluation
-   * and a DoFHandler in order to allow for reading out the finite element
-   * data. It uses the data provided by dof_handler.get_fe(). If the element
-   * is vector-valued, the optional argument allows to specify the index of
-   * the base element (as long as the element is primitive, non-primitive are
-   * not supported currently).
-   *
-   * With initialization from a FEValues object, no call to a reinit method of
-   * this class is necessary. Instead, it is enough if the geometry is
-   * initialized to a given cell iterator. It can also read from or write to
-   * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
-   * types (which is less efficient with MPI since index translation has to be
-   * done), but of course only for one cell at a time. Hence, a kernel using
-   * this method does not vectorize over several elements (which is most
-   * efficient for vector operations), but only possibly within the element if
-   * the evaluate/integrate routines are combined (e.g. for computing cell
-   * matrices).
-   */
-  FEEvaluationQ_DG0 (const MappingFEEvaluation<dim,Number> &geometry,
-                     const DoFHandler<dim>                 &dof_handler,
-                     const unsigned int                     first_selected_component = 0);
-
-  /**
-   * Copy constructor
-   */
-  FEEvaluationQ_DG0 (const FEEvaluationQ_DG0 &other);
-
-  /**
-   * Evaluates the function values, the gradients, and the Hessians of the FE
-   * function given at the DoF values in the input vector at the quadrature
-   * points of the unit cell. The function arguments specify which parts shall
-   * actually be computed. Needs to be called before the functions @p
-   * get_value(), @p get_gradient() or @p get_laplacian give useful
-   * information (unless these values have been set manually).
-   */
-  void evaluate (const bool evaluate_val,
-                 const bool evaluate_grad,
-                 const bool evaluate_lapl = false);
-
-  /**
-   * This function takes the values and/or gradients that are stored on
-   * quadrature points, tests them by all the basis functions/gradients on the
-   * cell and performs the cell integration. The two function arguments @p
-   * integrate_val and @p integrate_grad are used to enable/disable some of
-   * values or gradients.
-   */
-  void integrate (const bool integrate_val,
-                  const bool integrate_grad);
-};
-
-
 /*----------------------- Inline functions ----------------------------------*/
 
 #ifndef DOXYGEN
@@ -1946,13 +1572,13 @@ FEEvaluationBase<dim,n_components_,Number>
 ::FEEvaluationBase (const MatrixFree<dim,Number> &data_in,
                     const unsigned int fe_no_in,
                     const unsigned int quad_no_in,
-                    const unsigned int dofs_per_cell,
+                    const unsigned int fe_degree,
                     const unsigned int n_q_points)
   :
   quad_no            (quad_no_in),
   n_fe_components    (data_in.get_dof_info(fe_no_in).n_components),
-  active_fe_index    (data_in.get_dof_info(fe_no_in).fe_index_from_dofs_per_cell
-                      (dofs_per_cell * n_fe_components)),
+  active_fe_index    (data_in.get_dof_info(fe_no_in).fe_index_from_degree
+                      (fe_degree)),
   active_quad_index  (data_in.get_mapping_info().
                       mapping_data_gen[quad_no_in].
                       quad_index_from_n_q_points(n_q_points)),
@@ -2065,7 +1691,7 @@ FEEvaluationBase<dim,n_components_,Number>
   active_quad_index  (other.active_quad_index),
   matrix_info        (other.matrix_info),
   dof_info           (other.dof_info),
-  mapping_info       (other.mapping_info), 
+  mapping_info       (other.mapping_info),
   stored_shape_info  (other.stored_shape_info),
   data               (other.data),
   cartesian_data     (other.cartesian_data),
@@ -2601,14 +2227,15 @@ FEEvaluationBase<dim,n_components_,Number>
             }
           else
             {
-              // no constraint at all: loop bounds are known, compiler can
-              // unroll without checks
+              // no constraint at all: compiler can unroll at least the
+              // vectorization loop
               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)
-                  operation.process_dof (dof_indices[j], *src[comp],
-                                         local_data[comp][j]);
+              for (unsigned int j=0; j<n_local_dofs; j+=VectorizedArray<Number>::n_array_elements)
+                for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+                  for (unsigned int comp=0; comp<n_components; ++comp)
+                    operation.process_dof (dof_indices[j+v], *src[comp],
+                                           local_data[comp][j+v]);
             }
         }
 
@@ -2743,13 +2370,14 @@ FEEvaluationBase<dim,n_components_,Number>
             }
           else
             {
-              // no constraint at all: loop bounds are known, compiler can
-              // unroll without checks
+              // no constraint at all: compiler can unroll at least the
+              // vectorization loop
               AssertDimension (dof_info->end_indices(cell)-dof_indices,
                                static_cast<int>(n_local_dofs));
-              for (unsigned int j=0; j<n_local_dofs; ++j)
-                operation.process_dof (dof_indices[j], *src[0],
-                                       local_data[j]);
+              for (unsigned int j=0; j<n_local_dofs; j+=VectorizedArray<Number>::n_array_elements)
+                for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+                  operation.process_dof (dof_indices[j+v], *src[0],
+                                         local_data[j+v]);
             }
         }
 
@@ -3826,11 +3454,11 @@ FEEvaluationAccess<dim,n_components_,Number>
 ::FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
                       const unsigned int fe_no,
                       const unsigned int quad_no_in,
-                      const unsigned int dofs_per_cell,
+                      const unsigned int fe_degree,
                       const unsigned int n_q_points)
   :
   FEEvaluationBase <dim,n_components_,Number>
-  (data_in, fe_no, quad_no_in, dofs_per_cell, n_q_points)
+  (data_in, fe_no, quad_no_in, fe_degree, n_q_points)
 {}
 
 
@@ -3867,11 +3495,11 @@ FEEvaluationAccess<dim,1,Number>
 ::FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
                       const unsigned int fe_no,
                       const unsigned int quad_no_in,
-                      const unsigned int dofs_per_cell,
+                      const unsigned int fe_degree,
                       const unsigned int n_q_points)
   :
   FEEvaluationBase <dim,1,Number>
-  (data_in, fe_no, quad_no_in, dofs_per_cell, n_q_points)
+  (data_in, fe_no, quad_no_in, fe_degree, n_q_points)
 {}
 
 
@@ -4101,11 +3729,11 @@ FEEvaluationAccess<dim,dim,Number>
 ::FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
                       const unsigned int fe_no,
                       const unsigned int quad_no_in,
-                      const unsigned int dofs_per_cell,
+                      const unsigned int fe_degree,
                       const unsigned int n_q_points)
   :
   FEEvaluationBase <dim,dim,Number>
-  (data_in, fe_no, quad_no_in, dofs_per_cell, n_q_points)
+  (data_in, fe_no, quad_no_in, fe_degree, n_q_points)
 {}
 
 
@@ -4448,278 +4076,140 @@ FEEvaluationAccess<dim,dim,Number>
 
 
 
-/*----------------------- FEEvaluationGeneral -------------------------------*/
-
-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)
-  :
-  BaseClass (data_in, fe_no, quad_no, dofs_per_cell, n_q_points)
+namespace internal
 {
-  check_template_arguments(fe_no);
-  set_data_pointers();
-}
+  /**
+   * In this namespace, the evaluator routines that evaluate the tensor
+   * products are implemented.
+   */
+  enum EvaluatorVariant
+    {
+      evaluate_general,
+      evaluate_symmetric,
+      evaluate_evenodd
+    };
 
+  /**
+   * Generic evaluator framework
+   */
+  template <EvaluatorVariant variant, int dim, int fe_degree, int n_q_points_1d,
+            typename Number>
+  struct EvaluatorTensorProduct
+  {};
 
+  /**
+   * Internal evaluator for 1d-3d shape function using the tensor product form
+   * of the basis functions
+   */
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+  struct EvaluatorTensorProduct<evaluate_general,dim,fe_degree,n_q_points_1d,Number>
+  {
+    static const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
+    static const unsigned int n_q_points = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
+
+    /**
+     * Empty constructor. Does nothing. Be careful when using 'values' and
+     * related methods because they need to be filled with the other pointer
+     */
+    EvaluatorTensorProduct ()
+      :
+      shape_values (0),
+      shape_gradients (0),
+      shape_hessians (0)
+    {}
+
+    /**
+     * Constructor, taking the data from ShapeInfo
+     */
+    EvaluatorTensorProduct (const AlignedVector<Number> &shape_values,
+                            const AlignedVector<Number> &shape_gradients,
+                            const AlignedVector<Number> &shape_hessians)
+      :
+      shape_values (shape_values.begin()),
+      shape_gradients (shape_gradients.begin()),
+      shape_hessians (shape_hessians.begin())
+    {}
+
+    template <int direction, bool dof_to_quad, bool add>
+    void
+    values (const Number in [],
+            Number       out[]) const
+    {
+      apply<direction,dof_to_quad,add>(shape_values, in, out);
+    }
 
-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 MappingFEEvaluation<dim,Number> &geometry,
-                       const DoFHandler<dim>                 &dof_handler,
-                       const unsigned int                     first_selected_component)
-  :
-  BaseClass (geometry, dof_handler, first_selected_component)
-{
-  check_template_arguments(numbers::invalid_unsigned_int);
-  set_data_pointers();
-}
+    template <int direction, bool dof_to_quad, bool add>
+    void
+    gradients (const Number in [],
+               Number       out[]) const
+    {
+      apply<direction,dof_to_quad,add>(shape_gradients, in, out);
+    }
 
+    template <int direction, bool dof_to_quad, bool add>
+    void
+    hessians (const Number in [],
+              Number       out[]) const
+    {
+      apply<direction,dof_to_quad,add>(shape_hessians, in, out);
+    }
 
+    template <int direction, bool dof_to_quad, bool add>
+    static void apply (const Number *shape_data,
+                       const Number in [],
+                       Number       out []);
 
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
-          typename Number>
-inline
-FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluationGeneral (const FEEvaluationGeneral &other)
-  :
-  BaseClass (other)
-{
-  set_data_pointers();
-}
+    const Number * shape_values;
+    const Number * shape_gradients;
+    const Number * shape_hessians;
+  };
 
+  // evaluates the given shape data in 1d-3d using the tensor product
+  // form. does not use a particular layout of entries in the matrices
+  // like the functions below and corresponds to a usual matrix-matrix
+  // product
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+  template <int direction, bool dof_to_quad, bool add>
+  inline
+  void
+  EvaluatorTensorProduct<evaluate_general,dim,fe_degree,n_q_points_1d,Number>
+  ::apply (const Number *shape_data,
+           const Number in [],
+           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 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;
 
-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>
-::check_template_arguments(const unsigned int fe_no)
-{
-#ifdef DEBUG
-  // print error message when the dimensions do not match. Propose a possible
-  // fix
-  if (fe_degree != this->data->fe_degree
-      ||
-      n_q_points != this->data->n_q_points)
-    {
-      std::string message =
-        "-------------------------------------------------------\n";
-      message += "Illegal arguments in constructor/wrong template arguments!\n";
-      message += "    Called -->   FEEvaluation<dim,";
-      message += Utilities::int_to_string(fe_degree) + ",";
-      message += Utilities::int_to_string(n_q_points_1d);
-      message += "," + Utilities::int_to_string(n_components);
-      message += ",Number>(data, ";
-      message += Utilities::int_to_string(fe_no) + ", ";
-      message += Utilities::int_to_string(this->quad_no) + ")\n";
-
-      // check whether some other vector component has the correct number of
-      // points
-      unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
-                   proposed_quad_comp = numbers::invalid_unsigned_int;
-      if (fe_no != numbers::invalid_unsigned_int)
-        {
-          if (fe_degree == this->data->fe_degree)
-            proposed_dof_comp = fe_no;
-          else
-            for (unsigned int no=0; no<this->matrix_info->n_components(); ++no)
-              if (this->matrix_info->get_shape_info(no,0,this->active_fe_index,0).fe_degree
-                  == fe_degree)
-                {
-                  proposed_dof_comp = no;
-                  break;
-                }
-          if (n_q_points ==
-              this->mapping_info->mapping_data_gen[this->quad_no].n_q_points[this->active_quad_index])
-            proposed_quad_comp = this->quad_no;
-          else
-            for (unsigned int no=0; no<this->mapping_info->mapping_data_gen.size(); ++no)
-              if (this->mapping_info->mapping_data_gen[no].n_q_points[this->active_quad_index]
-                  == n_q_points)
-                {
-                  proposed_quad_comp = no;
-                  break;
-                }
-        }
-      if (proposed_dof_comp  != numbers::invalid_unsigned_int &&
-          proposed_quad_comp != numbers::invalid_unsigned_int)
-        {
-          if (proposed_dof_comp != fe_no)
-            message += "Wrong vector component selection:\n";
-          else
-            message += "Wrong quadrature formula selection:\n";
-          message += "    Did you mean FEEvaluation<dim,";
-          message += Utilities::int_to_string(fe_degree) + ",";
-          message += Utilities::int_to_string(n_q_points_1d);
-          message += "," + Utilities::int_to_string(n_components);
-          message += ",Number>(data, ";
-          message += Utilities::int_to_string(proposed_dof_comp) + ", ";
-          message += Utilities::int_to_string(proposed_quad_comp) + ")?\n";
-          std::string correct_pos;
-          if (proposed_dof_comp != fe_no)
-            correct_pos = " ^ ";
-          else
-            correct_pos = "   ";
-          if (proposed_quad_comp != this->quad_no)
-            correct_pos += " ^\n";
-          else
-            correct_pos += "  \n";
-          message += "                                                     " + correct_pos;
-        }
-      // ok, did not find the numbers specified by the template arguments in
-      // the given list. Suggest correct template arguments
-      const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(std::pow(1.001*this->data->n_q_points,1./dim));
-      message += "Wrong template arguments:\n";
-      message += "    Did you mean FEEvaluation<dim,";
-      message += Utilities::int_to_string(this->data->fe_degree) + ",";
-      message += Utilities::int_to_string(proposed_n_q_points_1d);
-      message += "," + Utilities::int_to_string(n_components);
-      message += ",Number>(data, ";
-      message += Utilities::int_to_string(fe_no) + ", ";
-      message += Utilities::int_to_string(this->quad_no) + ")?\n";
-      std::string correct_pos;
-      if (this->data->fe_degree != fe_degree)
-        correct_pos = " ^";
-      else
-        correct_pos = "  ";
-      if (proposed_n_q_points_1d != n_q_points_1d)
-        correct_pos += " ^\n";
-      else
-        correct_pos += "  \n";
-      message += "                                 " + correct_pos;
-
-      Assert (fe_degree == this->data->fe_degree &&
-              n_q_points == this->data->n_q_points,
-              ExcMessage(message));
-    }
-  if (fe_no != numbers::invalid_unsigned_int)
-    {
-      AssertDimension (n_q_points,
-                       this->mapping_info->mapping_data_gen[this->quad_no].
-                       n_q_points[this->active_quad_index]);
-      AssertDimension (this->data->dofs_per_cell * this->n_fe_components,
-                       this->dof_info->dofs_per_cell[this->active_fe_index]);
-    }
-#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>
-::set_data_pointers()
-{
-  AssertIndexRange(this->data->dofs_per_cell, dofs_per_cell+2);
-  const unsigned int desired_dofs_per_cell = this->data->dofs_per_cell;
-
-  // set the pointers to the correct position in the data array
-  for (unsigned int c=0; c<n_components_; ++c)
-    {
-      this->values_dofs[c] = &my_data_array[c*desired_dofs_per_cell];
-      this->values_quad[c] = &my_data_array[n_components*desired_dofs_per_cell+c*n_q_points];
-      for (unsigned int d=0; d<dim; ++d)
-        this->gradients_quad[c][d] = &my_data_array[n_components*(desired_dofs_per_cell+
-                                                                  n_q_points)
-                                                    +
-                                                    (c*dim+d)*n_q_points];
-      for (unsigned int d=0; d<(dim*dim+dim)/2; ++d)
-        this->hessians_quad[c][d] = &my_data_array[n_components*((dim+1)*n_q_points+
-                                                                 desired_dofs_per_cell)
-                                                   +
-                                                   (c*(dim*dim+dim)+d)*n_q_points];
-    }
-}
-
-
-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>
-::check_dofs_per_cell (const unsigned int given_dofs_per_cell) const
-{
-#ifdef DEBUG
-  if (given_dofs_per_cell != this->data->dofs_per_cell)
-    {
-      std::ostringstream str;
-      str << "Dofs per cell in FEEvaluation* class does not match the element ("
-          << given_dofs_per_cell << " != " << this->data->dofs_per_cell << "). "
-          << "Suggestion: ";
-      if (this->data->dofs_per_cell == dofs_per_cell)
-        str << "FEEvaluation/FEEvaluationGeneral";
-      else if (this->data->dofs_per_cell ==
-               internal::MatrixFreeFunctions::DGP_dofs_per_cell<dim,fe_degree>::value)
-        str << "FEEvaluationDGP";
-      else if (this->data->dofs_per_cell == dofs_per_cell+1)
-        str << "FEEvaluationQ_DG0";
-      else
-        str << "No matching suggestion found";
-
-      Assert(given_dofs_per_cell == this->data->dofs_per_cell,
-             ExcMessage(str.str().c_str()));
-    }
-#endif         
-}
-
-
-
-namespace internal
-{
-  // evaluates the given shape data in 1d-3d using the tensor product
-  // form. does not use a particular layout of entries in the matrices
-  // like the functions below 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 Number *shape_data,
-                        const Number in [],
-                        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 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 (int i2=0; i2<n_blocks2; ++i2)
-      {
-        for (int i1=0; i1<n_blocks1; ++i1)
-          {
-            for (int col=0; col<nn; ++col)
-              {
-                Number val0;
-                if (dof_to_quad == true)
-                  val0 = shape_data[col];
-                else
-                  val0 = shape_data[col*n_q_points_1d];
-                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;
-              }
+    for (int i2=0; i2<n_blocks2; ++i2)
+      {
+        for (int i1=0; i1<n_blocks1; ++i1)
+          {
+            for (int col=0; col<nn; ++col)
+              {
+                Number val0;
+                if (dof_to_quad == true)
+                  val0 = shape_data[col];
+                else
+                  val0 = shape_data[col*n_q_points_1d];
+                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
@@ -4840,11 +4330,51 @@ namespace internal
 
 
 
-  // This method specializes the general application of tensor-product based
+  // This class 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):
+  // are symmetric about 0.5 and the quadrature points are, too.
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+  struct EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
+  {
+    static const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
+    static const unsigned int n_q_points = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
+
+    /**
+     * Constructor, taking the data from ShapeInfo
+     */
+    EvaluatorTensorProduct (const AlignedVector<Number> &shape_values,
+                            const AlignedVector<Number> &shape_gradients,
+                            const AlignedVector<Number> &shape_hessians)
+      :
+      shape_values (shape_values.begin()),
+      shape_gradients (shape_gradients.begin()),
+      shape_hessians (shape_hessians.begin())
+    {}
+
+    template <int direction, bool dof_to_quad, bool add>
+    void
+    values (const Number in [],
+            Number       out[]) const;
+
+    template <int direction, bool dof_to_quad, bool add>
+    void
+    gradients (const Number in [],
+               Number       out[]) const;
+
+    template <int direction, bool dof_to_quad, bool add>
+    void
+    hessians (const Number in [],
+              Number       out[]) const;
+
+    const Number * shape_values;
+    const Number * shape_gradients;
+    const Number * shape_hessians;
+  };
+
+
+
+  // In this 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 ]
@@ -4861,13 +4391,13 @@ namespace internal
   // 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.
-  template <int dim, int fe_degree, int n_q_points_1d, typename Number,
-            int direction, bool dof_to_quad, bool add>
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+  template <int direction, bool dof_to_quad, bool add>
   inline
   void
-  apply_tensor_product_values (const Number *shape_values,
-                               const Number in [],
-                               Number       out [])
+  EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
+  ::values (const Number in [],
+            Number       out []) const
   {
     AssertIndexRange (direction, dim);
     const int mm     = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
@@ -5041,10 +4571,6 @@ namespace internal
 
 
 
-  // evaluates the given shape data in 1d-3d using the tensor product
-  // form assuming the symmetries of unit cell shape gradients for
-  // finite elements in FEEvaluation
-
   // 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):
@@ -5064,13 +4590,13 @@ namespace internal
   // 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.
-  template <int dim, int fe_degree, int n_q_points_1d, typename Number,
-            int direction, bool dof_to_quad, bool add>
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+  template <int direction, bool dof_to_quad, bool add>
   inline
   void
-  apply_tensor_product_gradients (const Number *shape_gradients,
-                                  const Number in [],
-                                  Number       out [])
+  EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
+  ::gradients (const Number in [],
+               Number       out []) const
   {
     AssertIndexRange (direction, dim);
     const int mm     = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
@@ -5208,13 +4734,13 @@ namespace internal
   // evaluates the given shape data in 1d-3d using the tensor product
   // form assuming the symmetries of unit cell shape hessians for
   // finite elements in FEEvaluation
-  template <int dim, int fe_degree, int n_q_points_1d, typename Number,
-            int direction, bool dof_to_quad, bool add>
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+  template <int direction, bool dof_to_quad, bool add>
   inline
   void
-  apply_tensor_product_hessians (const Number *shape_hessians,
-                                 const Number in [],
-                                 Number       out [])
+  EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
+  ::hessians (const Number in [],
+              Number       out []) const
   {
     AssertIndexRange (direction, dim);
     const int mm     = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
@@ -5360,7 +4886,7 @@ namespace internal
 
 
 
-  // This method implements a different approach to the symmetric case for
+  // This class implements a different approach to the symmetric case for
   // values, gradients, and Hessians also treated with the above functions: It
   // is possible to reduce the cost per dimension from N^2 to N^2/2, where N
   // is the number of 1D dofs (there are only N^2/2 different entries in the
@@ -5373,13 +4899,80 @@ namespace internal
   // experiments in the book say that the method is not efficient for N<20, it
   // is more efficient in the context where the loop bounds are compile-time
   // constants (templates).
-  template <int dim, int fe_degree, int n_q_points_1d, typename Number,
-            int direction, bool dof_to_quad, bool add, int type>
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+  struct EvaluatorTensorProduct<evaluate_evenodd,dim,fe_degree,n_q_points_1d,Number>
+  {
+    static const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
+    static const unsigned int n_q_points = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
+
+    /**
+     * Empty constructor. Does nothing. Be careful when using 'values' and
+     * related methods because they need to be filled with the other pointer
+     */
+    EvaluatorTensorProduct ()
+      :
+      shape_values (0),
+      shape_gradients (0),
+      shape_hessians (0)
+    {}
+
+    /**
+     * Constructor, taking the data from ShapeInfo (using the even-odd
+     * variants stored there)
+     */
+    EvaluatorTensorProduct (const AlignedVector<Number> &shape_values,
+                            const AlignedVector<Number> &shape_gradients,
+                            const AlignedVector<Number> &shape_hessians)
+      :
+      shape_values (shape_values.begin()),
+      shape_gradients (shape_gradients.begin()),
+      shape_hessians (shape_hessians.begin())
+    {}
+
+    template <int direction, bool dof_to_quad, bool add>
+    void
+    values (const Number in [],
+            Number       out[]) const
+    {
+      apply<direction,dof_to_quad,add,0>(shape_values, in, out);
+    }
+
+    template <int direction, bool dof_to_quad, bool add>
+    void
+    gradients (const Number in [],
+               Number       out[]) const
+    {
+      apply<direction,dof_to_quad,add,1>(shape_gradients, in, out);
+    }
+
+    template <int direction, bool dof_to_quad, bool add>
+    void
+    hessians (const Number in [],
+              Number       out[]) const
+    {
+      apply<direction,dof_to_quad,add,2>(shape_hessians, in, out);
+    }
+
+    template <int direction, bool dof_to_quad, bool add, int type>
+    static void apply (const Number *shape_data,
+                       const Number  in [],
+                       Number        out []);
+
+    const Number * shape_values;
+    const Number * shape_gradients;
+    const Number * shape_hessians;
+  };
+
+
+
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+  template <int direction, bool dof_to_quad, bool add, int type>
   inline
   void
-  apply_tensor_product_evenodd (const Number shapes [][(n_q_points_1d+1)/2],
-                                const Number in [],
-                                Number       out [])
+  EvaluatorTensorProduct<evaluate_evenodd,dim,fe_degree,n_q_points_1d,Number>
+  ::apply (const Number *shapes,
+           const Number  in [],
+           Number        out [])
   {
     AssertIndexRange (type, 3);
     AssertIndexRange (direction, dim);
@@ -5392,6 +4985,8 @@ namespace internal
     const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
     const int stride    = Utilities::fixed_int_power<nn,direction>::value;
 
+    const int offset = (n_q_points_1d+1)/2;
+
     // this code may look very inefficient at first sight due to the many
     // different cases with if's at the innermost loop part, but all of the
     // conditionals can be evaluated at compile time because they are
@@ -5421,25 +5016,25 @@ namespace internal
                   {
                     if (dof_to_quad == true)
                       {
-                        r0 = shapes[0][col]         * xp[0];
-                        r1 = shapes[fe_degree][col] * xm[0];
+                        r0 = shapes[col]                    * xp[0];
+                        r1 = shapes[fe_degree*offset + col] * xm[0];
                       }
                     else
                       {
-                        r0 = shapes[col][0]           * xp[0];
-                        r1 = shapes[fe_degree-col][0] * xm[0];
+                        r0 = shapes[col*offset]             * xp[0];
+                        r1 = shapes[(fe_degree-col)*offset] * xm[0];
                       }
                     for (int ind=1; ind<mid; ++ind)
                       {
                         if (dof_to_quad == true)
                           {
-                            r0 += shapes[ind][col]           * xp[ind];
-                            r1 += shapes[fe_degree-ind][col] * xm[ind];
+                            r0 += shapes[ind*offset+col]             * xp[ind];
+                            r1 += shapes[(fe_degree-ind)*offset+col] * xm[ind];
                           }
                         else
                           {
-                            r0 += shapes[col][ind]           * xp[ind];
-                            r1 += shapes[fe_degree-col][ind] * xm[ind];
+                            r0 += shapes[col*offset+ind]             * xp[ind];
+                            r1 += shapes[(fe_degree-col)*offset+ind] * xm[ind];
                           }
                       }
                   }
@@ -5448,12 +5043,12 @@ namespace internal
                 if (mm % 2 == 1 && dof_to_quad == true)
                   {
                     if (type == 1)
-                      r1 += shapes[mid][col] * in[stride*mid];
+                      r1 += shapes[mid*offset+col] * in[stride*mid];
                     else
-                      r0 += shapes[mid][col] * in[stride*mid];
+                      r0 += shapes[mid*offset+col] * in[stride*mid];
                   }
                 else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0))
-                  r0 += shapes[col][mid] * in[stride*mid];
+                  r0 += shapes[col*offset+mid] * in[stride*mid];
 
                 if (add == false)
                   {
@@ -5484,14 +5079,14 @@ namespace internal
                 Number r0;
                 if (mid > 0)
                   {
-                    r0  = shapes[0][n_cols] * xp[0];
+                    r0  = shapes[n_cols] * xp[0];
                     for (int ind=1; ind<mid; ++ind)
-                      r0 += shapes[ind][n_cols] * xp[ind];
+                      r0 += shapes[ind*offset+n_cols] * xp[ind];
                   }
                 else
                   r0 = Number();
                 if (type != 1 && mm % 2 == 1)
-                  r0 += shapes[mid][n_cols] * in[stride*mid];
+                  r0 += shapes[mid*offset+n_cols] * in[stride*mid];
 
                 if (add == false)
                   out[stride*n_cols]  = r0;
@@ -5505,15 +5100,15 @@ namespace internal
                   {
                     if (type == 1)
                       {
-                        r0 = shapes[n_cols][0] * xm[0];
+                        r0 = shapes[n_cols*offset] * xm[0];
                         for (int ind=1; ind<mid; ++ind)
-                          r0 += shapes[n_cols][ind] * xm[ind];
+                          r0 += shapes[n_cols*offset+ind] * xm[ind];
                       }
                     else
                       {
-                        r0 = shapes[n_cols][0] * xp[0];
+                        r0 = shapes[n_cols*offset] * xp[0];
                         for (int ind=1; ind<mid; ++ind)
-                          r0 += shapes[n_cols][ind] * xp[ind];
+                          r0 += shapes[n_cols*offset+ind] * xp[ind];
                       }
                   }
                 else
@@ -5522,7 +5117,7 @@ namespace internal
                 if (type == 0 && mm % 2 == 1)
                   r0 += in[stride*mid];
                 else if (type == 2 && mm % 2 == 1)
-                  r0 += shapes[n_cols][mid] * in[stride*mid];
+                  r0 += shapes[n_cols*offset+mid] * in[stride*mid];
 
                 if (add == false)
                   out[stride*n_cols]  = r0;
@@ -5558,1284 +5153,751 @@ namespace internal
 
 
 
-  // evaluates the given shape data in 1d-3d using the tensor product
-  // form assuming the symmetries of unit cell shape gradients for
-  // finite elements in FEEvaluationGL
-
-  // This function 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     ]
-  template <int dim, int fe_degree, typename Number,
-            int direction, bool dof_to_quad, bool add>
+  // Select evaluator type from element shape function type
+  template <MatrixFreeFunctions::ElementType element, bool is_long>
+  struct EvaluatorSelector {};
+
+  template <bool is_long>
+  struct EvaluatorSelector<MatrixFreeFunctions::tensor_general,is_long>
+  {
+    static const EvaluatorVariant variant = evaluate_general;
+  };
+
+  template <>
+  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric,false>
+  {
+    static const EvaluatorVariant variant = evaluate_symmetric;
+  };
+
+  template <> struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric,true>
+  {
+    static const EvaluatorVariant variant = evaluate_evenodd;
+  };
+
+  template <bool is_long>
+  struct EvaluatorSelector<MatrixFreeFunctions::truncated_tensor,is_long>
+  {
+    static const EvaluatorVariant variant = evaluate_general;
+  };
+
+  template <>
+  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,false>
+  {
+    static const EvaluatorVariant variant = evaluate_symmetric;
+  };
+
+  template <>
+  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,true>
+  {
+    static const EvaluatorVariant variant = evaluate_evenodd;
+  };
+
+  template <bool is_long>
+  struct EvaluatorSelector<MatrixFreeFunctions::tensor_gausslobatto,is_long>
+  {
+    static const EvaluatorVariant variant = evaluate_evenodd;
+  };
+
+
+
+  // This struct performs the evaluation of function values, gradients and
+  // Hessians for tensor-product finite elements. The operation is used for
+  // both the symmetric and non-symmetric case, which use different apply
+  // functions 'values', 'gradients' in the individual coordinate
+  // directions. The apply functions for values are provided through one of
+  // the template classes EvaluatorTensorProduct which in turn are selected
+  // from the MatrixFreeFunctions::ElementType template argument.
+  //
+  // There is a specialization made for Gauss-Lobatto elements further down
+  // where the 'values' operation is identity, which allows us to write
+  // shorter code.
+  template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
+            int n_q_points_1d, int n_components, typename Number>
+  struct FEEvaluationImpl
+  {
+    static
+    void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+                   VectorizedArray<Number> *values_dofs_actual[],
+                   VectorizedArray<Number> *values_quad[],
+                   VectorizedArray<Number> *gradients_quad[][dim],
+                   VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+                   const bool               evaluate_val,
+                   const bool               evaluate_grad,
+                   const bool               evaluate_lapl);
+
+    static
+    void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+                    VectorizedArray<Number> *values_dofs_actual[],
+                    VectorizedArray<Number> *values_quad[],
+                    VectorizedArray<Number> *gradients_quad[][dim],
+                    const bool               evaluate_val,
+                    const bool               evaluate_grad);
+  };
+
+
+  template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
+            int n_q_points_1d, int n_components, typename Number>
   inline
   void
-  apply_tensor_product_gradients_gl (const Number *shape_gradients,
-                                     const Number in [],
-                                     Number       out [])
+  FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
+  ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+              VectorizedArray<Number> *values_dofs_actual[],
+              VectorizedArray<Number> *values_quad[],
+              VectorizedArray<Number> *gradients_quad[][dim],
+              VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+              const bool               evaluate_val,
+              const bool               evaluate_grad,
+              const bool               evaluate_lapl)
   {
-    AssertIndexRange (direction, dim);
-    const int mm     = fe_degree+1;
-    const int nn     = fe_degree+1;
-    const int n_cols = nn / 2;
-    const int mid    = mm / 2;
+    const EvaluatorVariant variant =
+      EvaluatorSelector<type,(fe_degree+n_q_points_1d>4)>::variant;
+    typedef EvaluatorTensorProduct<variant, dim, fe_degree, n_q_points_1d,
+                                   VectorizedArray<Number> > Eval;
+    Eval eval (variant == evaluate_evenodd ? shape_info.shape_val_evenodd :
+               shape_info.shape_values,
+               variant == evaluate_evenodd ? shape_info.shape_gra_evenodd :
+               shape_info.shape_gradients,
+               variant == evaluate_evenodd ? shape_info.shape_hes_evenodd :
+               shape_info.shape_hessians);
+
+    const unsigned int temp_size = Eval::dofs_per_cell > Eval::n_q_points ?
+                                   Eval::dofs_per_cell : Eval::n_q_points;
+
+    VectorizedArray<Number> **values_dofs = values_dofs_actual;
+    VectorizedArray<Number> data_array[type!=MatrixFreeFunctions::truncated_tensor ? 1 :
+                                       n_components*Eval::dofs_per_cell];
+    VectorizedArray<Number> *expanded_dof_values[n_components];
+    if (type == MatrixFreeFunctions::truncated_tensor)
+      {
+        for (unsigned int c=0; c<n_components; ++c)
+          expanded_dof_values[c] = &data_array[c*Eval::dofs_per_cell];
+        values_dofs = expanded_dof_values;
 
-    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;
+        unsigned int count_p = 0, count_q = 0;
+        for (unsigned int i=0; i<(dim>2?fe_degree+1:1); ++i)
+          {
+            for (unsigned int j=0; j<(dim>1?fe_degree+1-i:1); ++j)
+              {
+                for (unsigned int k=0; k<fe_degree+1-j-i; ++k, ++count_p, ++count_q)
+                  for (unsigned int c=0; c<n_components; ++c)
+                    expanded_dof_values[c][count_q] = values_dofs_actual[c][count_p];
+                for (unsigned int k=fe_degree+1-j-i; k<fe_degree+1; ++k, ++count_q)
+                  for (unsigned int c=0; c<n_components; ++c)
+                    expanded_dof_values[c][count_q] = VectorizedArray<Number>();
+              }
+            for (unsigned int j=fe_degree+1-i; j<fe_degree+1; ++j)
+              for (unsigned int k=0; k<fe_degree+1; ++k, ++count_q)
+                for (unsigned int c=0; c<n_components; ++c)
+                  expanded_dof_values[c][count_q] = VectorizedArray<Number>();
+          }
+        AssertDimension(count_q, Eval::dofs_per_cell);
+      }
 
-    for (int i2=0; i2<n_blocks2; ++i2)
+    // These avoid compiler errors; they are only used in sensible context but
+    // compilers typically cannot detect when we access something like
+    // gradients_quad[2] only for dim==3.
+    const unsigned int d1 = dim>1?1:0;
+    const unsigned int d2 = dim>2?2:0;
+    const unsigned int d3 = dim>2?3:0;
+    const unsigned int d4 = dim>2?4:0;
+    const unsigned int d5 = dim>2?5:0;
+
+    switch (dim)
       {
-        for (int i1=0; i1<n_blocks1; ++i1)
+      case 1:
+        for (unsigned int c=0; c<n_components; c++)
           {
-            for (int col=0; col<n_cols; ++col)
+            if (evaluate_val == true)
+              eval.template values<0,true,false> (values_dofs[c], values_quad[c]);
+            if (evaluate_grad == true)
+              eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]);
+            if (evaluate_lapl == true)
+              eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]);
+          }
+        break;
+
+      case 2:
+        for (unsigned int c=0; c<n_components; c++)
+          {
+            VectorizedArray<Number> temp1[temp_size];
+            VectorizedArray<Number> temp2[temp_size];
+
+            // grad x
+            if (evaluate_grad == true)
               {
-                Number val0, val1, in0, in1, res0, res1;
-                if (mid > 0)
+                eval.template gradients<0,true,false> (values_dofs[c], temp1);
+                eval.template values<1,true,false> (temp1, gradients_quad[c][0]);
+              }
+            if (evaluate_lapl == true)
+              {
+                // grad xy
+                if (evaluate_grad == false)
+                  eval.template gradients<0,true,false>(values_dofs[c], temp1);
+                eval.template gradients<1,true,false>  (temp1, hessians_quad[c][d1+d1]);
+
+                // grad xx
+                eval.template hessians<0,true,false>(values_dofs[c], temp1);
+                eval.template values<1,true,false>  (temp1, hessians_quad[c][0]);
+              }
+
+            // grad y
+            eval.template values<0,true,false> (values_dofs[c], temp1);
+            if (evaluate_grad == true)
+              eval.template gradients<1,true,false> (temp1, gradients_quad[c][d1]);
+
+            // grad yy
+            if (evaluate_lapl == true)
+              eval.template hessians<1,true,false> (temp1, hessians_quad[c][d1]);
+
+            // val: can use values applied in x
+            if (evaluate_val == true)
+              eval.template values<1,true,false> (temp1, values_quad[c]);
+          }
+        break;
+
+      case 3:
+        for (unsigned int c=0; c<n_components; c++)
+          {
+            VectorizedArray<Number> temp1[temp_size];
+            VectorizedArray<Number> temp2[temp_size];
+
+            if (evaluate_grad == true)
+              {
+                // grad x
+                eval.template gradients<0,true,false> (values_dofs[c], temp1);
+                eval.template values<1,true,false> (temp1, temp2);
+                eval.template values<2,true,false> (temp2, gradients_quad[c][0]);
+              }
+
+            if (evaluate_lapl == true)
+              {
+                // grad xz
+                if (evaluate_grad == false)
                   {
-                    if (dof_to_quad == true)
-                      {
-                        val0 = shape_gradients[col];
-                        val1 = shape_gradients[nn-1-col];
-                      }
-                    else
-                      {
-                        val0 = shape_gradients[col*mm];
-                        val1 = 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)
-                          {
-                            val0 = shape_gradients[ind*mm+col];
-                            val1 = shape_gradients[ind*mm+nn-1-col];
-                          }
-                        else
-                          {
-                            val0 = shape_gradients[col*mm+ind];
-                            val1 = shape_gradients[(nn-col-1)*mm+ind];
-                          }
-
-                        // 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 * in0;
-                            res0 -= val1 * in1;
-                          }
-                        else
-                          {
-                            res0 += val0 * in0;
-                            res1 += val1 * in0;
-                            res0 -= val1 * in1;
-                            res1 -= val0 * in1;
-                          }
-                      }
-                  }
-                else
-                  res0 = res1 = Number();
-                if (mm % 2 == 1)
-                  {
-                    if (dof_to_quad == true)
-                      val0 = shape_gradients[mid*mm+col];
-                    else
-                      val0 = shape_gradients[col*mm+mid];
-                    val1 = val0 * in[stride*mid];
-                    res0 += val1;
-                    res1 -= val1;
-                  }
-                if (add == false)
-                  {
-                    out[stride*col]         = res0;
-                    out[stride*(nn-1-col)]  = res1;
-                  }
-                else
-                  {
-                    out[stride*col]        += res0;
-                    out[stride*(nn-1-col)] += res1;
+                    eval.template gradients<0,true,false> (values_dofs[c], temp1);
+                    eval.template values<1,true,false> (temp1, temp2);
                   }
+                eval.template gradients<2,true,false> (temp2, hessians_quad[c][d4]);
+
+                // grad xy
+                eval.template gradients<1,true,false> (temp1, temp2);
+                eval.template values<2,true,false> (temp2, hessians_quad[c][d3]);
+
+                // grad xx
+                eval.template hessians<0,true,false>(values_dofs[c], temp1);
+                eval.template values<1,true,false>  (temp1, temp2);
+                eval.template values<2,true,false>  (temp2, hessians_quad[c][0]);
               }
-            if ( nn%2 == 1 )
+
+            // grad y
+            eval.template values<0,true,false> (values_dofs[c], temp1);
+            if (evaluate_grad == true)
               {
-                Number val0, res0;
-                if (dof_to_quad == true)
-                  val0 = shape_gradients[n_cols];
-                else
-                  val0 = shape_gradients[n_cols*mm];
-                if (mid > 0)
-                  {
-                    res0  = in[0] - in[stride*(mm-1)];
-                    res0 *= val0;
-                    for (int ind=1; ind<mid; ++ind)
-                      {
-                        if (dof_to_quad == true)
-                          val0 = shape_gradients[ind*mm+n_cols];
-                        else
-                          val0 = shape_gradients[n_cols*mm+ind];
-                        Number val1  = in[stride*ind] - in[stride*(mm-1-ind)];
-                        val1 *= val0;
-                        res0 += val1;
-                      }
-                  }
-                else
-                  res0 = Number();
-                if (add == false)
-                  out[stride*n_cols]  = res0;
-                else
-                  out[stride*n_cols] += res0;
+                eval.template gradients<1,true,false>(temp1, temp2);
+                eval.template values<2,true,false>   (temp2, gradients_quad[c][d1]);
               }
 
-            // increment: in regular case, just go to the next point in
-            // x-direction. for y-part in 3D and 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)
+            if (evaluate_lapl == true)
               {
-              case 0:
-                in += mm;
-                out += nn;
-                break;
-              case 1:
-              case 2:
-                ++in;
-                ++out;
-                break;
-              default:
-                Assert (false, ExcNotImplemented());
+                // grad yz
+                if (evaluate_grad == false)
+                  eval.template gradients<1,true,false>(temp1, temp2);
+                eval.template gradients<2,true,false>  (temp2, hessians_quad[c][d5]);
+
+                // grad yy
+                eval.template hessians<1,true,false> (temp1, temp2);
+                eval.template values<2,true,false> (temp2, hessians_quad[c][d1]);
               }
-          }
-
-        if (direction == 1)
-          {
-            in  += nn * (mm-1);
-            out += nn * (nn-1);
-          }
-      }
-  }
 
+            // grad z: can use the values applied in x direction stored in temp1
+            eval.template values<1,true,false> (temp1, temp2);
+            if (evaluate_grad == true)
+              eval.template gradients<2,true,false> (temp2, gradients_quad[c][d2]);
 
+            // grad zz: can use the values applied in x and y direction stored
+            // in temp2
+            if (evaluate_lapl == true)
+              eval.template hessians<2,true,false>(temp2, hessians_quad[c][d2]);
 
-  // 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.
-  // use different versions for 1d, 2d, 3d to avoid nasty compiler warnings
-  template <typename FEEval>
-  inline
-  void
-  do_evaluate (FEEval    &fe_eval,
-               VectorizedArray<typename FEEval::number_type>* values_dofs[],
-               VectorizedArray<typename FEEval::number_type>* values_quad[],
-               VectorizedArray<typename FEEval::number_type>* gradients_quad[][1],
-               VectorizedArray<typename FEEval::number_type>* hessians_quad[][1],
-               const bool evaluate_val,
-               const bool evaluate_grad,
-               const bool evaluate_lapl,
-               internal::int2type<1>)
-  {
-    const unsigned int n_components = FEEval::n_components;
+            // val: can use the values applied in x & y direction stored in temp2
+            if (evaluate_val == true)
+              eval.template values<2,true,false> (temp2, values_quad[c]);
+          }
+        break;
 
-    for (unsigned int c=0; c<n_components; c++)
-      {
-        if (evaluate_val == true)
-          fe_eval.template apply_values<0,true,false>
-            (values_dofs[c], fe_eval.values_quad[c]);
-        if (evaluate_grad == true)
-          fe_eval.template apply_gradients<0,true,false>
-            (values_dofs[c], fe_eval.gradients_quad[c][0]);
-        if (evaluate_lapl == true)
-          fe_eval.template apply_hessians<0,true,false>
-            (values_dofs[c], fe_eval.hessians_quad[c][0]);
+      default:
+        AssertThrow(false, ExcNotImplemented());
       }
+
+    // case additional dof for FE_Q_DG0: add values; gradients and second
+    // derivatives evaluate to zero
+    if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 && evaluate_val)
+      for (unsigned int c=0; c<n_components; ++c)
+        for (unsigned int q=0; q<Eval::n_q_points; ++q)
+          values_quad[c][q] += values_dofs[c][Eval::dofs_per_cell];
   }
 
 
-  template <typename FEEval>
+
+  template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
+            int n_q_points_1d, int n_components, typename Number>
   inline
   void
-  do_evaluate (FEEval    &fe_eval,
-               VectorizedArray<typename FEEval::number_type>* values_dofs[],
-               VectorizedArray<typename FEEval::number_type>* values_quad[],
-               VectorizedArray<typename FEEval::number_type>* gradients_quad[][2],
-               VectorizedArray<typename FEEval::number_type>* hessians_quad[][3],
-               const bool evaluate_val,
-               const bool evaluate_grad,
-               const bool evaluate_lapl,
-               internal::int2type<2>)
+  FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
+  ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+               VectorizedArray<Number> *values_dofs_actual[],
+               VectorizedArray<Number> *values_quad[],
+               VectorizedArray<Number> *gradients_quad[][dim],
+               const bool               integrate_val,
+               const bool               integrate_grad)
   {
-    const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ?
-                                   FEEval::dofs_per_cell : FEEval::n_q_points;
-    const unsigned int n_components = FEEval::n_components;
-
-    for (unsigned int c=0; c<n_components; c++)
+    const EvaluatorVariant variant =
+      EvaluatorSelector<type,(fe_degree+n_q_points_1d>4)>::variant;
+    typedef EvaluatorTensorProduct<variant, dim, fe_degree, n_q_points_1d,
+                                   VectorizedArray<Number> > Eval;
+    Eval eval (variant == evaluate_evenodd ? shape_info.shape_val_evenodd :
+               shape_info.shape_values,
+               variant == evaluate_evenodd ? shape_info.shape_gra_evenodd :
+               shape_info.shape_gradients,
+               variant == evaluate_evenodd ? shape_info.shape_hes_evenodd :
+               shape_info.shape_hessians);
+
+    const unsigned int temp_size = Eval::dofs_per_cell > Eval::n_q_points ?
+                                   Eval::dofs_per_cell : Eval::n_q_points;
+    VectorizedArray<Number> temp1[temp_size];
+    VectorizedArray<Number> temp2[temp_size];
+
+    // expand dof_values to tensor product for truncated tensor products
+    VectorizedArray<Number> ** values_dofs = values_dofs_actual;
+    VectorizedArray<Number> data_array[type!=MatrixFreeFunctions::truncated_tensor ? 1 :
+                                       n_components*Eval::dofs_per_cell];
+    VectorizedArray<Number> *expanded_dof_values[n_components];
+    if (type == MatrixFreeFunctions::truncated_tensor)
       {
-        VectorizedArray<typename FEEval::number_type> temp1[temp_size];
-        VectorizedArray<typename FEEval::number_type> temp2[temp_size];
-
-        // grad x
-        if (evaluate_grad == true)
-          {
-            fe_eval.template apply_gradients<0,true,false>
-              (values_dofs[c], temp1);
-            fe_eval.template apply_values<1,true,false>
-              (temp1, gradients_quad[c][0]);
-          }
-        if (evaluate_lapl == true)
-          {
-            // grad xy
-            if (evaluate_grad == false)
-              fe_eval.template apply_gradients<0,true,false>
-                (values_dofs[c], temp1);
-            fe_eval.template apply_gradients<1,true,false>
-              (temp1, hessians_quad[c][2]);
-
-            // grad xx
-            fe_eval.template apply_hessians<0,true,false>
-              (values_dofs[c], temp1);
-            fe_eval.template apply_values<1,true,false>
-              (temp1, hessians_quad[c][0]);
-          }
-
-        // grad y
-        fe_eval.template apply_values<0,true,false>
-          (values_dofs[c], temp1);
-        if (evaluate_grad == true)
-          fe_eval.template apply_gradients<1,true,false>
-            (temp1, gradients_quad[c][1]);
-
-        // grad yy
-        if (evaluate_lapl == true)
-          fe_eval.template apply_hessians<1,true,false>
-            (temp1, hessians_quad[c][1]);
-
-        // val: can use values applied in x
-        if (evaluate_val == true)
-          fe_eval.template apply_values<1,true,false>
-            (temp1, values_quad[c]);
+        for (unsigned int c=0; c<n_components; ++c)
+          expanded_dof_values[c] = &data_array[c*Eval::dofs_per_cell];
+        values_dofs = expanded_dof_values;
       }
-  }
 
-  template <typename FEEval>
-  inline
-  void
-  do_evaluate (FEEval    &fe_eval,
-               VectorizedArray<typename FEEval::number_type>* values_dofs[],
-               VectorizedArray<typename FEEval::number_type>* values_quad[],
-               VectorizedArray<typename FEEval::number_type>* gradients_quad[][3],
-               VectorizedArray<typename FEEval::number_type>* hessians_quad[][6],
-               const bool evaluate_val,
-               const bool evaluate_grad,
-               const bool evaluate_lapl,
-               internal::int2type<3>)
-  {
-    const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ?
-                                   FEEval::dofs_per_cell : FEEval::n_q_points;
-    const unsigned int n_components = FEEval::n_components;
+    // These avoid compiler errors; they are only used in sensible context but
+    // compilers typically cannot detect when we access something like
+    // gradients_quad[2] only for dim==3.
+    const unsigned int d1 = dim>1?1:0;
+    const unsigned int d2 = dim>2?2:0;
 
-    for (unsigned int c=0; c<n_components; c++)
+    switch (dim)
       {
-        VectorizedArray<typename FEEval::number_type> temp1[temp_size];
-        VectorizedArray<typename FEEval::number_type> temp2[temp_size];
-
-        if (evaluate_grad == true)
-          {
-            // grad x
-            fe_eval.template apply_gradients<0,true,false>
-              (values_dofs[c], temp1);
-            fe_eval.template apply_values<1,true,false>
-              (temp1, temp2);
-            fe_eval.template apply_values<2,true,false>
-              (temp2, gradients_quad[c][0]);
-          }
-
-        if (evaluate_lapl == true)
+      case 1:
+        for (unsigned int c=0; c<n_components; c++)
           {
-            // grad xz
-            if (evaluate_grad == false)
+            if (integrate_val == true)
+              eval.template values<0,false,false> (values_quad[c], values_dofs[c]);
+            if (integrate_grad == true)
               {
-                fe_eval.template apply_gradients<0,true,false>
-                  (values_dofs[c], temp1);
-                fe_eval.template apply_values<1,true,false>
-                  (temp1, temp2);
+                if (integrate_val == true)
+                  eval.template gradients<0,false,true> (gradients_quad[c][0], values_dofs[c]);
+                else
+                  eval.template gradients<0,false,false> (gradients_quad[c][0], values_dofs[c]);
               }
-            fe_eval.template apply_gradients<2,true,false>
-              (temp2, 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, hessians_quad[c][3]);
-
-            // grad xx
-            fe_eval.template apply_hessians<0,true,false>
-              (values_dofs[c], temp1);
-            fe_eval.template apply_values<1,true,false>
-              (temp1, temp2);
-            fe_eval.template apply_values<2,true,false>
-              (temp2, hessians_quad[c][0]);
           }
+        break;
 
-        // grad y
-        fe_eval.template apply_values<0,true,false>
-          (values_dofs[c], temp1);
-        if (evaluate_grad == true)
+      case 2:
+        for (unsigned int c=0; c<n_components; c++)
           {
-            fe_eval.template apply_gradients<1,true,false>
-              (temp1, temp2);
-            fe_eval.template apply_values<2,true,false>
-              (temp2, gradients_quad[c][1]);
+            if (integrate_val == true)
+              {
+                // val
+                eval.template values<0,false,false> (values_quad[c], temp1);
+                //grad x
+                if (integrate_grad == true)
+                  eval.template gradients<0,false,true> (gradients_quad[c][0], temp1);
+                eval.template values<1,false,false>(temp1, values_dofs[c]);
+              }
+            if (integrate_grad == true)
+              {
+                // grad y
+                eval.template values<0,false,false>  (gradients_quad[c][d1], temp1);
+                if (integrate_val == false)
+                  {
+                    eval.template gradients<1,false,false>(temp1, values_dofs[c]);
+                    //grad x
+                    eval.template gradients<0,false,false> (gradients_quad[c][0], temp1);
+                    eval.template values<1,false,true> (temp1, values_dofs[c]);
+                  }
+                else
+                  eval.template gradients<1,false,true>(temp1, values_dofs[c]);
+              }
           }
+        break;
 
-        if (evaluate_lapl == true)
+      case 3:
+        for (unsigned int c=0; c<n_components; c++)
           {
-            // 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, 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, hessians_quad[c][1]);
+            if (integrate_val == true)
+              {
+                // val
+                eval.template values<0,false,false> (values_quad[c], temp1);
+                //grad x: can sum to temporary value in temp1
+                if (integrate_grad == true)
+                  eval.template gradients<0,false,true> (gradients_quad[c][0], temp1);
+                eval.template values<1,false,false>(temp1, temp2);
+                eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
+                if (integrate_grad == true)
+                  eval.template gradients<1,false,true>(temp1, temp2);
+                eval.template values<2,false,false> (temp2, values_dofs[c]);
+              }
+            else if (integrate_grad == true)
+              {
+                eval.template gradients<0,false,false>(gradients_quad[c][0], temp1);
+                eval.template values<1,false,false> (temp1, temp2);
+                eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
+                eval.template gradients<1,false,true>(temp1, temp2);
+                eval.template values<2,false,false> (temp2, values_dofs[c]);
+              }
+            if (integrate_grad == true)
+              {
+                // grad z: can sum to temporary x and y value in output
+                eval.template values<0,false,false> (gradients_quad[c][d2], temp1);
+                eval.template values<1,false,false> (temp1, temp2);
+                eval.template gradients<2,false,true> (temp2, values_dofs[c]);
+              }
           }
+        break;
 
-        // 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, gradients_quad[c][2]);
-
-        // 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, hessians_quad[c][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, values_quad[c]);
+      default:
+        AssertThrow(false, ExcNotImplemented());
       }
-  }
 
+    // case FE_Q_DG0: add values, gradients and second derivatives are zero
+    if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0)
+      {
+        if (integrate_val)
+          for (unsigned int c=0; c<n_components; ++c)
+            {
+              values_dofs[c][Eval::dofs_per_cell] = values_quad[c][0];
+              for (unsigned int q=1; q<Eval::n_q_points; ++q)
+                values_dofs[c][Eval::dofs_per_cell] += values_quad[c][q];
+            }
+        else
+          for (unsigned int c=0; c<n_components; ++c)
+            values_dofs[c][Eval::dofs_per_cell] = VectorizedArray<Number>();
+      }
 
-
-  template <typename FEEval>
-  inline
-  void
-  do_integrate (FEEval    &fe_eval,
-                VectorizedArray<typename FEEval::number_type>* values_dofs[],
-                VectorizedArray<typename FEEval::number_type>* values_quad[],
-                VectorizedArray<typename FEEval::number_type>* gradients_quad[][1],
-                const bool integrate_val,
-                const bool integrate_grad,
-                internal::int2type<1>)
-  {
-    const unsigned int n_components = FEEval::n_components;
-
-    for (unsigned int c=0; c<n_components; c++)
+    if (type == MatrixFreeFunctions::truncated_tensor)
       {
-        if (integrate_grad == true)
-          fe_eval.template apply_gradients<0,false,false>
-            (gradients_quad[c][0], fe_eval.values_dofs[c]);
-        if (integrate_val == true)
+        unsigned int count_p = 0, count_q = 0;
+        for (unsigned int i=0; i<(dim>2?fe_degree+1:1); ++i)
           {
-            if (integrate_grad == true)
-              fe_eval.template apply_values<0,false,true>
-                (values_quad[c], values_dofs[c]);
-            else
-              fe_eval.template apply_values<0,false,false>
-                (values_quad[c], values_dofs[c]);
+            for (unsigned int j=0; j<(dim>1?fe_degree+1-i:1); ++j)
+              {
+                for (unsigned int k=0; k<fe_degree+1-j-i; ++k, ++count_p, ++count_q)
+                  {
+                    for (unsigned int c=0; c<n_components; ++c)
+                      values_dofs_actual[c][count_p] = expanded_dof_values[c][count_q];
+                  }
+                count_q += j+i;
+              }
+            count_q += i*(fe_degree+1);
           }
+        AssertDimension(count_q, Eval::dofs_per_cell);
       }
   }
 
-  template <typename FEEval>
+  // This a specialization for Gauss-Lobatto elements where the 'values'
+  // operation is identity, which allows us to write shorter code.
+  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
+  struct FEEvaluationImpl<MatrixFreeFunctions::tensor_gausslobatto, dim,
+                          fe_degree, n_q_points_1d, n_components, Number>
+  {
+    static
+    void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+                   VectorizedArray<Number> *values_dofs[],
+                   VectorizedArray<Number> *values_quad[],
+                   VectorizedArray<Number> *gradients_quad[][dim],
+                   VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+                   const bool               evaluate_val,
+                   const bool               evaluate_grad,
+                   const bool               evaluate_lapl);
+
+    static
+    void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+                    VectorizedArray<Number> *values_dofs[],
+                    VectorizedArray<Number> *values_quad[],
+                    VectorizedArray<Number> *gradients_quad[][dim],
+                    const bool               integrate_val,
+                    const bool               integrate_grad);
+  };
+
+  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
   inline
   void
-  do_integrate (FEEval    &fe_eval,
-                VectorizedArray<typename FEEval::number_type>* values_dofs[],
-                VectorizedArray<typename FEEval::number_type>* values_quad[],
-                VectorizedArray<typename FEEval::number_type>* gradients_quad[][2],
-                const bool integrate_val,
-                const bool integrate_grad,
-                internal::int2type<2>)
+  FEEvaluationImpl<MatrixFreeFunctions::tensor_gausslobatto, dim,
+                   fe_degree, n_q_points_1d, n_components, Number>
+  ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+              VectorizedArray<Number> *values_dofs[],
+              VectorizedArray<Number> *values_quad[],
+              VectorizedArray<Number> *gradients_quad[][dim],
+              VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+              const bool               evaluate_val,
+              const bool               evaluate_grad,
+              const bool               evaluate_lapl)
   {
-    const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ?
-                                   FEEval::dofs_per_cell : FEEval::n_q_points;
-    const unsigned int n_components = FEEval::n_components;
-
-    for (unsigned int c=0; c<n_components; c++)
+    typedef EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree, fe_degree+1,
+                                   VectorizedArray<Number> > Eval;
+    Eval eval (shape_info.shape_val_evenodd, shape_info.shape_gra_evenodd,
+               shape_info.shape_hes_evenodd);
+
+    // These avoid compiler errors; they are only used in sensible context but
+    // compilers typically cannot detect when we access something like
+    // gradients_quad[2] only for dim==3.
+    const unsigned int d1 = dim>1?1:0;
+    const unsigned int d2 = dim>2?2:0;
+    const unsigned int d3 = dim>2?3:0;
+    const unsigned int d4 = dim>2?4:0;
+    const unsigned int d5 = dim>2?5:0;
+
+    switch (dim)
       {
-        VectorizedArray<typename FEEval::number_type> temp1[temp_size];
-        VectorizedArray<typename FEEval::number_type> temp2[temp_size];
+      case 1:
+        if (evaluate_val == true)
+          std::memcpy (values_quad[0], values_dofs[0],
+                       eval.dofs_per_cell * n_components *
+                       sizeof (values_dofs[0][0]));
+        for (unsigned int c=0; c<n_components; c++)
+          {
+            if (evaluate_grad == true)
+              eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]);
+            if (evaluate_lapl == true)
+              eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]);
+          }
+        break;
 
-        // val
-        if (integrate_val == true)
-          fe_eval.template apply_values<0,false,false>
-            (values_quad[c], temp1);
-        if (integrate_grad == true)
+      case 2:
+        if (evaluate_val == true)
           {
-            //grad x
-            if (integrate_val == true)
-              fe_eval.template apply_gradients<0,false,true>
-                (gradients_quad[c][0], temp1);
-            else
-              fe_eval.template apply_gradients<0,false,false>
-                (gradients_quad[c][0], temp1);
+            std::memcpy (values_quad[0], values_dofs[0],
+                         Eval::dofs_per_cell * n_components *
+                         sizeof (values_dofs[0][0]));
           }
-        if (integrate_val || integrate_grad)
-          fe_eval.template apply_values<1,false,false>
-            (temp1, values_dofs[c]);
-        if (integrate_grad == true)
+        if (evaluate_grad == true)
+          for (unsigned int comp=0; comp<n_components; comp++)
+            {
+              // grad x
+              eval.template gradients<0,true,false> (values_dofs[comp],
+                                                     gradients_quad[comp][0]);
+              // grad y
+              eval.template gradients<1,true,false> (values_dofs[comp],
+                                                     gradients_quad[comp][d1]);
+            }
+        if (evaluate_lapl == true)
+          for (unsigned int comp=0; comp<n_components; comp++)
+            {
+              // hess x
+              eval.template hessians<0,true,false> (values_dofs[comp],
+                                                    hessians_quad[comp][0]);
+              // hess y
+              eval.template hessians<1,true,false> (values_dofs[comp],
+                                                    hessians_quad[comp][d1]);
+
+              VectorizedArray<Number> temp1[Eval::dofs_per_cell];
+              // grad x grad y
+              eval.template gradients<0,true,false> (values_dofs[comp], temp1);
+              eval.template gradients<1,true,false> (temp1, hessians_quad[comp][d1+d1]);
+            }
+        break;
+
+      case 3:
+        if (evaluate_val == true)
           {
-            // grad y
-            fe_eval.template apply_values<0,false,false>
-              (gradients_quad[c][1], temp1);
-            fe_eval.template apply_gradients<1,false,true>
-              (temp1, values_dofs[c]);
+            std::memcpy (values_quad[0], values_dofs[0],
+                         Eval::dofs_per_cell * n_components *
+                         sizeof (values_dofs[0][0]));
           }
+        if (evaluate_grad == true)
+          for (unsigned int comp=0; comp<n_components; comp++)
+            {
+              // grad x
+              eval.template gradients<0,true,false> (values_dofs[comp],
+                                                     gradients_quad[comp][0]);
+              // grad y
+              eval.template gradients<1,true,false> (values_dofs[comp],
+                                                     gradients_quad[comp][d1]);
+              // grad y
+              eval.template gradients<2,true,false> (values_dofs[comp],
+                                                     gradients_quad[comp][d2]);
+            }
+        if (evaluate_lapl == true)
+          for (unsigned int comp=0; comp<n_components; comp++)
+            {
+              // grad x
+              eval.template hessians<0,true,false> (values_dofs[comp],
+                                                    hessians_quad[comp][0]);
+              // grad y
+              eval.template hessians<1,true,false> (values_dofs[comp],
+                                                    hessians_quad[comp][d1]);
+              // grad y
+              eval.template hessians<2,true,false> (values_dofs[comp],
+                                                    hessians_quad[comp][d2]);
+
+              VectorizedArray<Number> temp1[Eval::dofs_per_cell];
+              // grad xy
+              eval.template gradients<0,true,false> (values_dofs[comp], temp1);
+              eval.template gradients<1,true,false> (temp1, hessians_quad[comp][d3]);
+              // grad xz
+              eval.template gradients<2,true,false> (temp1, hessians_quad[comp][d4]);
+              // grad yz
+              eval.template gradients<1,true,false> (values_dofs[comp], temp1);
+              eval.template gradients<2,true,false> (temp1, hessians_quad[comp][d5]);
+            }
+        break;
+      default:
+        AssertThrow(false, ExcNotImplemented());
       }
   }
 
-  template <typename FEEval>
+  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
   inline
   void
-  do_integrate (FEEval    &fe_eval,
-                VectorizedArray<typename FEEval::number_type>* values_dofs[],
-                VectorizedArray<typename FEEval::number_type>* values_quad[],
-                VectorizedArray<typename FEEval::number_type>* gradients_quad[][3],
-                const bool integrate_val,
-                const bool integrate_grad,
-                internal::int2type<3>)
+  FEEvaluationImpl<MatrixFreeFunctions::tensor_gausslobatto, dim,
+                   fe_degree, n_q_points_1d, n_components, Number>
+  ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+               VectorizedArray<Number> *values_dofs[],
+               VectorizedArray<Number> *values_quad[],
+               VectorizedArray<Number> *gradients_quad[][dim],
+               const bool               integrate_val,
+               const bool               integrate_grad)
   {
-    const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ?
-                                   FEEval::dofs_per_cell : FEEval::n_q_points;
-    const unsigned int n_components = FEEval::n_components;
-
-    for (unsigned int c=0; c<n_components; c++)
+    typedef EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree, fe_degree+1,
+                                   VectorizedArray<Number> > Eval;
+    Eval eval (shape_info.shape_val_evenodd, shape_info.shape_gra_evenodd,
+               shape_info.shape_hes_evenodd);
+
+    // These avoid compiler errors; they are only used in sensible context but
+    // compilers typically cannot detect when we access something like
+    // gradients_quad[2] only for dim==3.
+    const unsigned int d1 = dim>1?1:0;
+    const unsigned int d2 = dim>2?2:0;
+
+    if (integrate_val == true)
+      std::memcpy (values_dofs[0], values_quad[0],
+                   Eval::dofs_per_cell * n_components *
+                   sizeof (values_dofs[0][0]));
+    switch (dim)
       {
-        VectorizedArray<typename FEEval::number_type> temp1[temp_size];
-        VectorizedArray<typename FEEval::number_type> temp2[temp_size];
-
-        if (integrate_val == true)
-          {
-            // val
-            fe_eval.template apply_values<0,false,false>
-              (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>
-                (gradients_quad[c][0], temp1);
-            else
-              fe_eval.template apply_gradients<0,false,false>
-                (gradients_quad[c][0], temp1);
-          }
-        if (integrate_val || integrate_grad)
-          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>
-              (gradients_quad[c][1], temp1);
-            fe_eval.template apply_gradients<1,false,true>
-              (temp1, temp2);
-          }
-        if (integrate_val || integrate_grad)
-          fe_eval.template apply_values<2,false,false>
-            (temp2, values_dofs[c]);
-        if (integrate_grad == true)
+      case 1:
+        for (unsigned int c=0; c<n_components; c++)
           {
-            // grad z: can sum to temporary x and y value in output
-            fe_eval.template apply_values<0,false,false>
-              (gradients_quad[c][2], temp1);
-            fe_eval.template apply_values<1,false,false>
-              (temp1, temp2);
-            fe_eval.template apply_gradients<2,false,true>
-              (temp2, values_dofs[c]);
+            if (integrate_grad == true)
+              {
+                if (integrate_val == true)
+                  eval.template gradients<0,false,true> (gradients_quad[c][0],
+                                                         values_dofs[c]);
+                else
+                  eval.template gradients<0,false,false> (gradients_quad[c][0],
+                                                          values_dofs[c]);
+              }
           }
-      }
-  }
-
-} // end of namespace internal
-
 
+        break;
+      case 2:
+        if (integrate_grad == true)
+          for (unsigned int comp=0; comp<n_components; comp++)
+            {
+              // grad x: If integrate_val == true we have to add to the
+              // previous output
+              if (integrate_val == true)
+                eval.template gradients<0, false, true> (gradients_quad[comp][0],
+                                                         values_dofs[comp]);
+              else
+                eval.template gradients<0, false, false> (gradients_quad[comp][0],
+                                                          values_dofs[comp]);
 
-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 (const bool evaluate_val,
-            const bool evaluate_grad,
-            const bool evaluate_lapl)
-{
-  this->check_dofs_per_cell(dofs_per_cell);
-  Assert (this->dof_values_initialized == true,
-          internal::ExcAccessToUninitializedField());
-  internal::do_evaluate (*this, this->values_dofs, this->values_quad,
-                         this->gradients_quad, this->hessians_quad,
-                         evaluate_val, evaluate_grad, evaluate_lapl,
-                         internal::int2type<dim>());
+              // grad y
+              eval.template gradients<1, false, true> (gradients_quad[comp][d1],
+                                                       values_dofs[comp]);
+            }
+        break;
 
-#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
-}
-
-
-
-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)
-{
-  this->check_dofs_per_cell(dofs_per_cell);
-  if (integrate_val == true)
-    Assert (this->values_quad_submitted == true,
-            internal::ExcAccessToUninitializedField());
-  if (integrate_grad == true)
-    Assert (this->gradients_quad_submitted == true,
-            internal::ExcAccessToUninitializedField());
-
-  internal::do_integrate (*this, this->values_dofs, this->values_quad,
-                          this->gradients_quad, integrate_val, integrate_grad,
-                          internal::int2type<dim>());
-
-#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
-Point<dim,VectorizedArray<Number> >
-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());
-  AssertIndexRange (q, n_q_points);
-
-  // Cartesian mesh: not all quadrature points are stored, only the
-  // diagonal. Hence, need to find the tensor product index and retrieve the
-  // value from that
-  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
-    {
-      Point<dim,VectorizedArray<Number> > point (false);
-      switch (dim)
-        {
-        case 1:
-          return this->quadrature_points[q];
-        case 2:
-          point[0] = this->quadrature_points[q%n_q_points_1d][0];
-          point[1] = this->quadrature_points[q/n_q_points_1d][1];
-          return point;
-        case 3:
-          point[0] = this->quadrature_points[q%n_q_points_1d][0];
-          point[1] = this->quadrature_points[(q/n_q_points_1d)%n_q_points_1d][1];
-          point[2] = this->quadrature_points[q/(n_q_points_1d*n_q_points_1d)][2];
-          return point;
-        default:
-          Assert (false, ExcNotImplemented());
-          return point;
-        }
-    }
-  // all other cases: just return the respective data as it is fully stored
-  else
-    return this->quadrature_points[q];
-}
-
-
-
-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_values(const VectorizedArray<Number> in [],
-               VectorizedArray<Number>       out [])
-{
-  internal::apply_tensor_product<dim,fe_degree,n_q_points_1d,
-           VectorizedArray<Number>, direction, dof_to_quad, add>
-           (this->data->shape_values.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_gradients(const VectorizedArray<Number> in [],
-                  VectorizedArray<Number>       out [])
-{
-  internal::apply_tensor_product<dim,fe_degree,n_q_points_1d,
-           VectorizedArray<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,
-           VectorizedArray<Number>, direction, dof_to_quad, add>
-           (this->data->shape_hessians.begin(), in, out);
-}
-
-
-/*-------------------------- FEEvaluation -----------------------------------*/
-
-
-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)
-  :
-  BaseClass (data_in, fe_no, quad_no)
-{
-  compute_even_odd_factors();
-}
-
-
-
-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 MappingFEEvaluation<dim,Number> &geometry,
-                const DoFHandler<dim>               &dof_handler,
-                const unsigned int                   first_selected_component)
-  :
-  BaseClass (geometry, dof_handler, first_selected_component)
-{
-  compute_even_odd_factors();
-}
-
-
-
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
-          typename Number>
-inline
-FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluation (const FEEvaluation &other)
-  :
-  BaseClass (other)
-{
-  compute_even_odd_factors();
-}
-
-
-
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
-          typename Number>
-inline
-void
-FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::compute_even_odd_factors()
-{
-  // check whether element is appropriate
-#ifdef DEBUG
-  const double zero_tol =
-    types_are_equal<Number,double>::value==true?1e-8:1e-7;
-  std::string error_message = "FEEvaluation not appropriate.\n";
-  error_message += "  It assumes symmetry of quadrature points w.r.t. 0.5 \n";
-  error_message += " and the basis functions starting from left and right.\n";
-  error_message += "Try FEEvaluationGeneral<...> instead!";
-
-  // symmetry for values
-  const unsigned int n_dofs_1d = fe_degree + 1;
-  for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
-    for (unsigned int j=0; j<n_q_points_1d; ++j)
-      Assert (std::fabs(this->data->shape_values[i*n_q_points_1d+j][0] -
-                        this->data->shape_values[(n_dofs_1d-i)*n_q_points_1d
-                                                -j-1][0]) < zero_tol,
-              ExcMessage(error_message));
-
-  // shape values should be zero at for all basis functions except for one
-  // where they are one in the middle
-  if (n_q_points_1d%2 == 1 && n_dofs_1d%2 == 1)
-    {
-      for (int i=0; i<static_cast<int>(n_dofs_1d/2); ++i)
-        Assert (std::fabs(this->data->shape_values[i*n_q_points_1d+
-                                                  n_q_points_1d/2][0]) < zero_tol,
-                ExcMessage(error_message));
-      Assert (std::fabs(this->data->shape_values[(n_dofs_1d/2)*n_q_points_1d+
-                                                n_q_points_1d/2][0]-1.)< zero_tol,
-              ExcMessage(error_message));
-    }
-
-  // skew-symmetry for gradient, zero of middle basis function in middle
-  // quadrature point
-  for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
-    for (unsigned int j=0; j<n_q_points_1d; ++j)
-      Assert (std::fabs(this->data->shape_gradients[i*n_q_points_1d+j][0] +
-                        this->data->shape_gradients[(n_dofs_1d-i)*n_q_points_1d-
-                                                   j-1][0]) < zero_tol,
-              ExcMessage(error_message));
-  if (n_dofs_1d%2 == 1 && n_q_points_1d%2 == 1)
-    Assert (std::fabs(this->data->shape_gradients[(n_dofs_1d/2)*n_q_points_1d+
-                                                 (n_q_points_1d/2)][0]) < zero_tol,
-            ExcMessage(error_message));
-
-
-  // symmetry for Laplacian
-  for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
-    for (unsigned int j=0; j<n_q_points_1d; ++j)
-      Assert (std::fabs(this->data->shape_hessians[i*n_q_points_1d+j][0] -
-                        this->data->shape_hessians[(n_dofs_1d-i)*n_q_points_1d-
-                                                  j-1][0]) < zero_tol,
-              ExcMessage(error_message));
-#endif
-
-  // Compute symmetric and skew-symmetric part of shape values for even-odd
-  // decomposition
-  for (int i=0; i<(fe_degree+1)/2; ++i)
-    for (unsigned int q=0; q<(n_q_points_1d+1)/2; ++q)
-      {
-        shape_val_evenodd[i][q] =
-          0.5 * (this->data->shape_values[i*n_q_points_1d+q] +
-                 this->data->shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
-        shape_val_evenodd[fe_degree-i][q] =
-          0.5 * (this->data->shape_values[i*n_q_points_1d+q] -
-                 this->data->shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
-
-        shape_gra_evenodd[i][q] =
-          0.5 * (this->data->shape_gradients[i*n_q_points_1d+q] +
-                 this->data->shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
-        shape_gra_evenodd[fe_degree-i][q] =
-          0.5 * (this->data->shape_gradients[i*n_q_points_1d+q] -
-                 this->data->shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
-
-        shape_hes_evenodd[i][q] =
-          0.5 * (this->data->shape_hessians[i*n_q_points_1d+q] +
-                 this->data->shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
-        shape_hes_evenodd[fe_degree-i][q] =
-          0.5 * (this->data->shape_hessians[i*n_q_points_1d+q] -
-                 this->data->shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
-      }
-  if (fe_degree % 2 == 0)
-    for (unsigned int q=0; q<(n_q_points_1d+1)/2; ++q)
-      {
-        shape_val_evenodd[fe_degree/2][q] =
-          this->data->shape_values[(fe_degree/2)*n_q_points_1d+q];
-        shape_gra_evenodd[fe_degree/2][q] =
-          this->data->shape_gradients[(fe_degree/2)*n_q_points_1d+q];
-        shape_hes_evenodd[fe_degree/2][q] =
-          this->data->shape_hessians[(fe_degree/2)*n_q_points_1d+q];
-      }
-}
-
-
-
-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 (const bool evaluate_val,
-            const bool evaluate_grad,
-            const bool evaluate_lapl)
-{
-  this->check_dofs_per_cell(dofs_per_cell);
-  Assert (this->dof_values_initialized == true,
-          internal::ExcAccessToUninitializedField());
-  internal::do_evaluate (*this, this->values_dofs, this->values_quad,
-                         this->gradients_quad, this->hessians_quad,
-                         evaluate_val, evaluate_grad, evaluate_lapl,
-                         internal::int2type<dim>());
-
-#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
-}
-
-
-
-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)
-{
-  this->check_dofs_per_cell(dofs_per_cell);
-  if (integrate_val == true)
-    Assert (this->values_quad_submitted == true,
-            internal::ExcAccessToUninitializedField());
-  if (integrate_grad == true)
-    Assert (this->gradients_quad_submitted == true,
-            internal::ExcAccessToUninitializedField());
-
-  internal::do_integrate (*this, this->values_dofs, this->values_quad,
-                          this->gradients_quad, integrate_val, integrate_grad,
-                          internal::int2type<dim>());
-
-#ifdef DEBUG
-  this->dof_values_initialized = true;
-#endif
-}
-
-
-
-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> in [],
-                VectorizedArray<Number>       out [])
-{
-  // for linear elements, the even-odd decomposition is slower than the plain
-  // evaluation (additions to create the symmetric and anti-symmetric part),
-  // for all other orders, we choose even-odd
-  if (fe_degree > 1 || n_q_points_1d > 3)
-    internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
-             VectorizedArray<Number>, direction, dof_to_quad, add, 0>
-             (&shape_val_evenodd[0], in, out);
-  else
-    internal::apply_tensor_product_values<dim,fe_degree,n_q_points_1d,
-             VectorizedArray<Number>, direction, dof_to_quad, add>
-             (this->data->shape_values.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
-FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::apply_gradients (const VectorizedArray<Number> in [],
-                   VectorizedArray<Number>       out [])
-{
-  if (fe_degree > 1 || n_q_points_1d > 3)
-    internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
-             VectorizedArray<Number>, direction, dof_to_quad, add, 1>
-             (&shape_gra_evenodd[0], in, out);
-  else
-    internal::apply_tensor_product_gradients<dim,fe_degree,n_q_points_1d,
-             VectorizedArray<Number>, direction, dof_to_quad, add>
-             (this->data->shape_gradients.begin(), in, out);
-}
-
-
-
-// Laplacian operator application. Very similar to value application because
-// the 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_,
-          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> in [],
-                  VectorizedArray<Number>       out [])
-{
-  if (fe_degree > 1 || n_q_points_1d > 3)
-    internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
-             VectorizedArray<Number>, direction, dof_to_quad, add, 2>
-             (&shape_hes_evenodd[0], in, out);
-  else
-    internal::apply_tensor_product_hessians<dim,fe_degree,n_q_points_1d,
-             VectorizedArray<Number>, direction, dof_to_quad, add>
-             (this->data->shape_hessians.begin(), in, out);
-}
-
-
-
-/*------------------------- FEEvaluationGL ----------------------------------*/
-
-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,
-                  const unsigned int fe_no,
-                  const unsigned int quad_no)
-  :
-  BaseClass (data_in, fe_no, quad_no)
-{
-#ifdef DEBUG
-  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 =
-    types_are_equal<Number,double>::value==true?1e-9:1e-7;
-
-  const unsigned int n_points_1d = fe_degree+1;
-  for (unsigned int i=0; i<n_points_1d; ++i)
-    for (unsigned int j=0; j<n_points_1d; ++j)
-      if (i!=j)
-        {
-          Assert (std::fabs(this->data->shape_values[i*n_points_1d+j][0])<zero_tol,
-                  ExcMessage (error_mess.c_str()));
-        }
-      else
-        {
-          Assert (std::fabs(this->data->shape_values[i*n_points_1d+
-                                                    j][0]-1.)<zero_tol,
-                  ExcMessage (error_mess.c_str()));
-        }
-  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>
-inline
-FEEvaluationGL<dim,fe_degree,n_components_,Number>
-::FEEvaluationGL (const MappingFEEvaluation<dim,Number> &geometry,
-                  const DoFHandler<dim>               &dof_handler,
-                  const unsigned int                   first_selected_component)
-  :
-  BaseClass (geometry, dof_handler, first_selected_component)
-{}
-
-
-
-template <int dim, int fe_degree, int n_components_, typename Number>
-inline
-FEEvaluationGL<dim,fe_degree,n_components_,Number>
-::FEEvaluationGL (const FEEvaluationGL &other)
-  :
-  BaseClass (other)
-{}
-
-
-
-template <int dim, int fe_degree, int n_components_, typename Number>
-inline
-void
-FEEvaluationGL<dim,fe_degree,n_components_,Number>
-::evaluate (const bool evaluate_val,
-            const bool evaluate_grad,
-            const bool evaluate_lapl)
-{
-  this->check_dofs_per_cell(dofs_per_cell);
-  Assert (this->cell != numbers::invalid_unsigned_int,
-          ExcNotInitialized());
-  Assert (this->dof_values_initialized == true,
-          internal::ExcAccessToUninitializedField());
-
-  if (evaluate_val == true)
-    {
-      std::memcpy (&this->values_quad[0][0], &this->values_dofs[0][0],
-                   dofs_per_cell * n_components *
-                   sizeof (this->values_dofs[0][0]));
-#ifdef DEBUG
-      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++)
-        {
-          if (dim == 3)
-            {
-              // grad x
-              apply_gradients<0,true,false> (this->values_dofs[comp],
-                                             this->gradients_quad[comp][0]);
-              // grad y
-              apply_gradients<1,true,false> (this->values_dofs[comp],
-                                             this->gradients_quad[comp][1]);
-              // grad y
-              apply_gradients<2,true,false> (this->values_dofs[comp],
-                                             this->gradients_quad[comp][2]);
-            }
-          else if (dim == 2)
-            {
-              // grad x
-              apply_gradients<0,true,false> (this->values_dofs[comp],
-                                             this->gradients_quad[comp][0]);
-              // grad y
-              apply_gradients<1,true,false> (this->values_dofs[comp],
-                                             this->gradients_quad[comp][1]);
-            }
-          else if (dim == 1)
-            apply_gradients<0,true,false> (this->values_dofs[comp],
-                                           this->gradients_quad[comp][0]);
-        }
-#ifdef DEBUG
-      this->gradients_quad_initialized = true;
-#endif
-    }
-  if (evaluate_lapl == true)
-    {
-      for (unsigned int comp=0; comp<n_components; comp++)
-        {
-          if (dim == 3)
-            {
-              // grad x
-              this->template apply_hessians<0,true,false> (this->values_dofs[comp],
-                                                           this->hessians_quad[comp][0]);
-              // grad y
-              this->template apply_hessians<1,true,false> (this->values_dofs[comp],
-                                                           this->hessians_quad[comp][1]);
-              // grad y
-              this->template apply_hessians<2,true,false> (this->values_dofs[comp],
-                                                           this->hessians_quad[comp][2]);
-
-              VectorizedArray<Number> temp1[n_q_points];
-              // grad xy
-              apply_gradients<0,true,false> (this->values_dofs[comp], temp1);
-              apply_gradients<1,true,false> (temp1, this->hessians_quad[comp][3]);
-              // grad xz
-              apply_gradients<2,true,false> (temp1, this->hessians_quad[comp][4]);
-              // grad yz
-              apply_gradients<1,true,false> (this->values_dofs[comp], temp1);
-              apply_gradients<2,true,false> (temp1, this->hessians_quad[comp][5]);
-            }
-          else if (dim == 2)
-            {
-              // grad x
-              this->template apply_hessians<0,true,false> (this->values_dofs[comp],
-                                                           this->hessians_quad[comp][0]);
-              // grad y
-              this->template apply_hessians<1,true,false> (this->values_dofs[comp],
-                                                           this->hessians_quad[comp][1]);
-              VectorizedArray<Number> temp1[n_q_points];
-              // grad xy
-              apply_gradients<0,true,false> (this->values_dofs[comp], temp1);
-              apply_gradients<1,true,false> (temp1, this->hessians_quad[comp][2]);
-            }
-          else if (dim == 1)
-            this->template apply_hessians<0,true,false> (this->values_dofs[comp],
-                                                       this->hessians_quad[comp][0]);
-        }
-#ifdef DEBUG
-      this->hessians_quad_initialized = true;
-#endif
-    }
-}
-
-
-
-template <int dim, int fe_degree, int n_components_, typename Number>
-inline
-void
-FEEvaluationGL<dim,fe_degree,n_components_,Number>
-::integrate (const bool integrate_val, const bool integrate_grad)
-{
-  this->check_dofs_per_cell(dofs_per_cell);
-  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());
-  if (integrate_val == true)
-    std::memcpy (&this->values_dofs[0][0], &this->values_quad[0][0],
-                 dofs_per_cell * n_components *
-                 sizeof (this->values_dofs[0][0]));
-  if (integrate_grad == true)
-    {
-      for (unsigned int comp=0; comp<n_components; comp++)
-        {
-          if (dim == 3)
-            {
-              // grad x: If integrate_val == true we have to add to the previous output
-              if (integrate_val == true)
-                apply_gradients<0, false, true> (this->gradients_quad[comp][0],
-                                                 this->values_dofs[comp]);
-              else
-                apply_gradients<0, false, false> (this->gradients_quad[comp][0],
-                                                  this->values_dofs[comp]);
-
-              // grad y: can sum to temporary x value in temp2
-              apply_gradients<1, false, true> (this->gradients_quad[comp][1],
-                                               this->values_dofs[comp]);
-
-              // grad z: can sum to temporary x and y value in output
-              apply_gradients<2, false, true> (this->gradients_quad[comp][2],
-                                               this->values_dofs[comp]);
-            }
-          else if (dim == 2)
-            {
-              // grad x: If integrate_val == true we have to add to the previous output
-              if (integrate_val == true)
-                apply_gradients<0, false, true> (this->gradients_quad[comp][0],
-                                                 this->values_dofs[comp]);
-              else
-                apply_gradients<0, false, false> (this->gradients_quad[comp][0],
-                                                  this->values_dofs[comp]);
-
-              // grad y: can sum to temporary x value in temp2
-              apply_gradients<1, false, true> (this->gradients_quad[comp][1],
-                                               this->values_dofs[comp]);
-            }
-          else if (dim == 1)
-            {
-              if (integrate_val == true)
-                apply_gradients<0, false, true> (this->gradients_quad[comp][0],
-                                                 this->values_dofs[comp]);
-              else
-                apply_gradients<0, false, false> (this->gradients_quad[comp][0],
-                                                  this->values_dofs[comp]);
-
-            }
-        }
-    }
+      case 3:
+        if (integrate_grad == true)
+          for (unsigned int comp=0; comp<n_components; comp++)
+            {
+              // grad x: If integrate_val == true we have to add to the
+              // previous output
+              if (integrate_val == true)
+                eval.template gradients<0, false, true> (gradients_quad[comp][0],
+                                                         values_dofs[comp]);
+              else
+                eval.template gradients<0, false, false> (gradients_quad[comp][0],
+                                                          values_dofs[comp]);
 
-#ifdef DEBUG
-  this->dof_values_initialized = true;
-#endif
-}
+              // grad y
+              eval.template gradients<1, false, true> (gradients_quad[comp][d1],
+                                                       values_dofs[comp]);
 
+              // grad z
+              eval.template gradients<2, false, true> (gradients_quad[comp][d2],
+                                                       values_dofs[comp]);
+            }
+        break;
+        
+      default:
+        AssertThrow(false, ExcNotImplemented());
+      }
+  }
 
+} // end of namespace internal
 
-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> in [],
-                   VectorizedArray<Number>       out [])
-{
-  internal::apply_tensor_product_gradients_gl<dim,fe_degree,
-           VectorizedArray<Number>, direction, dof_to_quad, add>
-           (this->data->shape_gradients.begin(), in, out);
-}
 
 
+/*-------------------------- FEEvaluation -----------------------------------*/
 
-/*------------------------- FEEvaluationDGP ---------------------------------*/
 
-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
-FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluationDGP (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)
+  BaseClass (data_in, fe_no, quad_no, fe_degree, n_q_points),
+  dofs_per_cell (this->data->dofs_per_cell)
 {
-  AssertDimension(static_cast<unsigned int>(this->values_quad[0]-this->values_dofs[0]),
-                  n_components * dofs_per_cell);
+  check_template_arguments(fe_no);
+  set_data_pointers();
 }
 
 
 
-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
-FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluationDGP (const MappingFEEvaluation<dim,Number> &geometry,
-                const DoFHandler<dim>                  &dof_handler,
-                const unsigned int                      first_selected_component)
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluation (const MappingFEEvaluation<dim,Number> &geometry,
+                const DoFHandler<dim>               &dof_handler,
+                const unsigned int                   first_selected_component)
   :
-  BaseClass (geometry, dof_handler, first_selected_component)
+  BaseClass (geometry, dof_handler, first_selected_component),
+  dofs_per_cell (this->data->dofs_per_cell)
 {
-  AssertDimension(static_cast<unsigned int>(this->values_quad[0]-this->values_dofs[0]),
-                  n_components * dofs_per_cell);
+  check_template_arguments(numbers::invalid_unsigned_int);
+  set_data_pointers();
 }
 
 
 
-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
-FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluationDGP (const FEEvaluationDGP &other)
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluation (const FEEvaluation &other)
   :
-  BaseClass (other)
+  BaseClass (other),
+  dofs_per_cell (this->data->dofs_per_cell)
 {
+  set_data_pointers();
 }
 
 
@@ -6844,148 +5906,249 @@ template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
           typename Number>
 inline
 void
-FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::evaluate (const bool evaluate_val,
-            const bool evaluate_grad,
-            const bool evaluate_lapl)
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::set_data_pointers()
 {
-  this->check_dofs_per_cell(dofs_per_cell);
-  Assert (this->dof_values_initialized == true,
-          internal::ExcAccessToUninitializedField());
-
-  // expand dof_values to tensor product
-  VectorizedArray<Number> data_array[n_components*BaseClass::dofs_per_cell];
-  VectorizedArray<Number> *expanded_dof_values[n_components];
-  for (unsigned int c=0; c<n_components; ++c)
-    expanded_dof_values[c] = &data_array[c*BaseClass::dofs_per_cell];
+  AssertIndexRange(this->data->dofs_per_cell, tensor_dofs_per_cell+2);
+  const unsigned int desired_dofs_per_cell = this->data->dofs_per_cell;
 
-  unsigned int count_p = 0, count_q = 0;
-  for (unsigned int i=0; i<(dim>2?fe_degree+1:1); ++i)
+  // set the pointers to the correct position in the data array
+  for (unsigned int c=0; c<n_components_; ++c)
     {
-      for (unsigned int j=0; j<(dim>1?fe_degree+1-i:1); ++j)
-        {
-          for (unsigned int k=0; k<fe_degree+1-j-i; ++k, ++count_p, ++count_q)
-            for (unsigned int c=0; c<n_components; ++c)
-              expanded_dof_values[c][count_q] = this->values_dofs[c][count_p];
-          for (unsigned int k=fe_degree+1-j-i; k<fe_degree+1; ++k, ++count_q)
-            for (unsigned int c=0; c<n_components; ++c)
-              expanded_dof_values[c][count_q] = VectorizedArray<Number>();
-        }
-      for (unsigned int j=fe_degree+1-i; j<fe_degree+1; ++j)
-        for (unsigned int k=0; k<fe_degree+1; ++k, ++count_q)
-          for (unsigned int c=0; c<n_components; ++c)
-            expanded_dof_values[c][count_q] = VectorizedArray<Number>();
+      this->values_dofs[c] = &my_data_array[c*desired_dofs_per_cell];
+      this->values_quad[c] = &my_data_array[n_components*desired_dofs_per_cell+c*n_q_points];
+      for (unsigned int d=0; d<dim; ++d)
+        this->gradients_quad[c][d] = &my_data_array[n_components*(desired_dofs_per_cell+
+                                                                  n_q_points)
+                                                    +
+                                                    (c*dim+d)*n_q_points];
+      for (unsigned int d=0; d<(dim*dim+dim)/2; ++d)
+        this->hessians_quad[c][d] = &my_data_array[n_components*((dim+1)*n_q_points+
+                                                                 desired_dofs_per_cell)
+                                                   +
+                                                   (c*(dim*dim+dim)+d)*n_q_points];
     }
-  AssertDimension(count_q, BaseClass::dofs_per_cell);
-  AssertDimension(count_p, dofs_per_cell);
 
-  internal::do_evaluate (*this, expanded_dof_values, this->values_quad,
-                         this->gradients_quad, this->hessians_quad,
-                         evaluate_val, evaluate_grad, evaluate_lapl,
-                         internal::int2type<dim>());
+  switch (this->data->element_type)
+    {
+    case internal::MatrixFreeFunctions::tensor_symmetric:
+      evaluate_funct =
+        internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric,
+                                   dim, fe_degree, n_q_points_1d, n_components_,
+                                   Number>::evaluate;
+      integrate_funct =
+        internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric,
+                                   dim, fe_degree, n_q_points_1d, n_components_,
+                                   Number>::integrate;
+      break;
 
-#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
-}
+    case internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0:
+      evaluate_funct =
+        internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
+                                   dim, fe_degree, n_q_points_1d, n_components_,
+                                   Number>::evaluate;
+      integrate_funct =
+        internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
+                                   dim, fe_degree, n_q_points_1d, n_components_,
+                                   Number>::integrate;
+      break;
 
+    case internal::MatrixFreeFunctions::tensor_general:
+      evaluate_funct =
+        internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
+                                   dim, fe_degree, n_q_points_1d, n_components_,
+                                   Number>::evaluate;
+      integrate_funct =
+        internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
+                                   dim, fe_degree, n_q_points_1d, n_components_,
+                                   Number>::integrate;
+      break;
 
+    case internal::MatrixFreeFunctions::tensor_gausslobatto:
+      evaluate_funct =
+        internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_gausslobatto,
+                                   dim, fe_degree, n_q_points_1d, n_components_,
+                                   Number>::evaluate;
+      integrate_funct =
+        internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_gausslobatto,
+                                   dim, fe_degree, n_q_points_1d, n_components_,
+                                   Number>::integrate;
+      break;
 
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
-          typename Number>
-inline
-void
-FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::integrate (bool integrate_val,bool integrate_grad)
-{
-  this->check_dofs_per_cell(dofs_per_cell);
-  if (integrate_val == true)
-    Assert (this->values_quad_submitted == true,
-            internal::ExcAccessToUninitializedField());
-  if (integrate_grad == true)
-    Assert (this->gradients_quad_submitted == true,
-            internal::ExcAccessToUninitializedField());
+    case internal::MatrixFreeFunctions::truncated_tensor:
+      evaluate_funct =
+        internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
+                                   dim, fe_degree, n_q_points_1d, n_components_,
+                                   Number>::evaluate;
+      integrate_funct =
+        internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
+                                   dim, fe_degree, n_q_points_1d, n_components_,
+                                   Number>::integrate;
+      break;
 
-  VectorizedArray<Number> data_array[n_components*BaseClass::dofs_per_cell];
-  VectorizedArray<Number> *expanded_dof_values[n_components];
-  for (unsigned int c=0; c<n_components; ++c)
-    expanded_dof_values[c] = &data_array[c*BaseClass::dofs_per_cell];
-  internal::do_integrate (*this, expanded_dof_values, this->values_quad,
-                          this->gradients_quad, integrate_val, integrate_grad,
-                          internal::int2type<dim>());
-
-  // truncate tensor product
-  unsigned int count_p = 0, count_q = 0;
-  for (unsigned int i=0; i<(dim>2?fe_degree+1:1); ++i)
-    {
-      for (unsigned int j=0; j<(dim>1?fe_degree+1-i:1); ++j)
-        {
-          for (unsigned int k=0; k<fe_degree+1-j-i; ++k, ++count_p, ++count_q)
-            {
-              for (unsigned int c=0; c<n_components; ++c)
-                this->values_dofs[c][count_p] = expanded_dof_values[c][count_q];
-            }
-          count_q += j+i;
-        }
-      count_q += i*(fe_degree+1);
+    default:
+      AssertThrow(false, ExcNotImplemented());
     }
-  AssertDimension(count_q, BaseClass::dofs_per_cell);
-  AssertDimension(count_p, dofs_per_cell);
 
-#ifdef DEBUG
-  this->dof_values_initialized = true;
-#endif
 }
 
 
 
-/*------------------------- FEEvaluationQ_DG0 -------------------------------*/
-
-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
-FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluationQ_DG0 (const MatrixFree<dim,Number> &data_in,
-                     const unsigned int fe_no,
-                     const unsigned int quad_no)
-  :
-  BaseClass (data_in, fe_no, quad_no)
+void
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::check_template_arguments(const unsigned int fe_no)
 {
-  AssertDimension(static_cast<unsigned int>(this->values_quad[0]-this->values_dofs[0]),
-                  n_components * dofs_per_cell);
-}
-
+#ifdef DEBUG
+  // print error message when the dimensions do not match. Propose a possible
+  // fix
+  if (fe_degree != this->data->fe_degree
+      ||
+      n_q_points != this->data->n_q_points)
+    {
+      std::string message =
+        "-------------------------------------------------------\n";
+      message += "Illegal arguments in constructor/wrong template arguments!\n";
+      message += "    Called -->   FEEvaluation<dim,";
+      message += Utilities::int_to_string(fe_degree) + ",";
+      message += Utilities::int_to_string(n_q_points_1d);
+      message += "," + Utilities::int_to_string(n_components);
+      message += ",Number>(data, ";
+      message += Utilities::int_to_string(fe_no) + ", ";
+      message += Utilities::int_to_string(this->quad_no) + ")\n";
 
+      // check whether some other vector component has the correct number of
+      // points
+      unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
+                   proposed_quad_comp = numbers::invalid_unsigned_int;
+      if (fe_no != numbers::invalid_unsigned_int)
+        {
+          if (fe_degree == this->data->fe_degree)
+            proposed_dof_comp = fe_no;
+          else
+            for (unsigned int no=0; no<this->matrix_info->n_components(); ++no)
+              if (this->matrix_info->get_shape_info(no,0,this->active_fe_index,0).fe_degree
+                  == fe_degree)
+                {
+                  proposed_dof_comp = no;
+                  break;
+                }
+          if (n_q_points ==
+              this->mapping_info->mapping_data_gen[this->quad_no].n_q_points[this->active_quad_index])
+            proposed_quad_comp = this->quad_no;
+          else
+            for (unsigned int no=0; no<this->mapping_info->mapping_data_gen.size(); ++no)
+              if (this->mapping_info->mapping_data_gen[no].n_q_points[this->active_quad_index]
+                  == n_q_points)
+                {
+                  proposed_quad_comp = no;
+                  break;
+                }
+        }
+      if (proposed_dof_comp  != numbers::invalid_unsigned_int &&
+          proposed_quad_comp != numbers::invalid_unsigned_int)
+        {
+          if (proposed_dof_comp != fe_no)
+            message += "Wrong vector component selection:\n";
+          else
+            message += "Wrong quadrature formula selection:\n";
+          message += "    Did you mean FEEvaluation<dim,";
+          message += Utilities::int_to_string(fe_degree) + ",";
+          message += Utilities::int_to_string(n_q_points_1d);
+          message += "," + Utilities::int_to_string(n_components);
+          message += ",Number>(data, ";
+          message += Utilities::int_to_string(proposed_dof_comp) + ", ";
+          message += Utilities::int_to_string(proposed_quad_comp) + ")?\n";
+          std::string correct_pos;
+          if (proposed_dof_comp != fe_no)
+            correct_pos = " ^ ";
+          else
+            correct_pos = "   ";
+          if (proposed_quad_comp != this->quad_no)
+            correct_pos += " ^\n";
+          else
+            correct_pos += "  \n";
+          message += "                                                     " + correct_pos;
+        }
+      // ok, did not find the numbers specified by the template arguments in
+      // the given list. Suggest correct template arguments
+      const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(std::pow(1.001*this->data->n_q_points,1./dim));
+      message += "Wrong template arguments:\n";
+      message += "    Did you mean FEEvaluation<dim,";
+      message += Utilities::int_to_string(this->data->fe_degree) + ",";
+      message += Utilities::int_to_string(proposed_n_q_points_1d);
+      message += "," + Utilities::int_to_string(n_components);
+      message += ",Number>(data, ";
+      message += Utilities::int_to_string(fe_no) + ", ";
+      message += Utilities::int_to_string(this->quad_no) + ")?\n";
+      std::string correct_pos;
+      if (this->data->fe_degree != fe_degree)
+        correct_pos = " ^";
+      else
+        correct_pos = "  ";
+      if (proposed_n_q_points_1d != n_q_points_1d)
+        correct_pos += " ^\n";
+      else
+        correct_pos += "  \n";
+      message += "                                 " + correct_pos;
 
-template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
-          typename Number>
-inline
-FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluationQ_DG0 (const MappingFEEvaluation<dim,Number> &geometry,
-                     const DoFHandler<dim>                  &dof_handler,
-                     const unsigned int                      first_selected_component)
-  :
-  BaseClass (geometry, dof_handler, first_selected_component)
-{
-  AssertDimension(static_cast<unsigned int>(this->values_quad[0]-this->values_dofs[0]),
-                  n_components * dofs_per_cell);
+      Assert (fe_degree == this->data->fe_degree &&
+              n_q_points == this->data->n_q_points,
+              ExcMessage(message));
+    }
+  if (fe_no != numbers::invalid_unsigned_int)
+    {
+      AssertDimension (n_q_points,
+                       this->mapping_info->mapping_data_gen[this->quad_no].
+                       n_q_points[this->active_quad_index]);
+      AssertDimension (this->data->dofs_per_cell * this->n_fe_components,
+                       this->dof_info->dofs_per_cell[this->active_fe_index]);
+    }
+#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
-FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluationQ_DG0 (const FEEvaluationQ_DG0 &other)
-  :
-  BaseClass (other)
+Point<dim,VectorizedArray<Number> >
+FEEvaluation<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());
+  AssertIndexRange (q, n_q_points);
+
+  // Cartesian mesh: not all quadrature points are stored, only the
+  // diagonal. Hence, need to find the tensor product index and retrieve the
+  // value from that
+  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
+    {
+      Point<dim,VectorizedArray<Number> > point (false);
+      switch (dim)
+        {
+        case 1:
+          return this->quadrature_points[q];
+        case 2:
+          point[0] = this->quadrature_points[q%n_q_points_1d][0];
+          point[1] = this->quadrature_points[q/n_q_points_1d][1];
+          return point;
+        case 3:
+          point[0] = this->quadrature_points[q%n_q_points_1d][0];
+          point[1] = this->quadrature_points[(q/n_q_points_1d)%n_q_points_1d][1];
+          point[2] = this->quadrature_points[q/(n_q_points_1d*n_q_points_1d)][2];
+          return point;
+        default:
+          Assert (false, ExcNotImplemented());
+          return point;
+        }
+    }
+  // all other cases: just return the respective data as it is fully stored
+  else
+    return this->quadrature_points[q];
 }
 
 
@@ -6994,22 +6157,19 @@ template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
           typename Number>
 inline
 void
-FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 ::evaluate (const bool evaluate_val,
             const bool evaluate_grad,
             const bool evaluate_lapl)
 {
-  this->check_dofs_per_cell(dofs_per_cell);
   Assert (this->dof_values_initialized == true,
           internal::ExcAccessToUninitializedField());
-  internal::do_evaluate (*this, this->values_dofs, this->values_quad,
-                         this->gradients_quad, this->hessians_quad,
-                         evaluate_val, evaluate_grad, evaluate_lapl,
-                         internal::int2type<dim>());
-  if (evaluate_val)
-    for (unsigned int c=0; c<n_components; ++c)
-      for (unsigned int q=0; q<n_q_points; ++q)
-        this->values_quad[c][q] += this->values_dofs[c][dofs_per_cell-1];
+
+  // Select algorithm matching the element type at run time (the function
+  // pointer is easy to predict, so negligible in cost)
+  evaluate_funct (*this->data, &this->values_dofs[0],
+                  this->values_quad, this->gradients_quad, this->hessians_quad,
+                  evaluate_val, evaluate_grad, evaluate_lapl);
 
 #ifdef DEBUG
   if (evaluate_val == true)
@@ -7027,10 +6187,9 @@ template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
           typename Number>
 inline
 void
-FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 ::integrate (bool integrate_val,bool integrate_grad)
 {
-  this->check_dofs_per_cell(dofs_per_cell);
   if (integrate_val == true)
     Assert (this->values_quad_submitted == true,
             internal::ExcAccessToUninitializedField());
@@ -7038,20 +6197,10 @@ FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
     Assert (this->gradients_quad_submitted == true,
             internal::ExcAccessToUninitializedField());
 
-  if (integrate_val)
-    for (unsigned int c=0; c<n_components; ++c)
-      {
-        this->values_dofs[c][dofs_per_cell-1] = this->values_quad[c][0];
-        for (unsigned int q=1; q<n_q_points; ++q)
-          this->values_dofs[c][dofs_per_cell-1] += this->values_quad[c][q];
-      }
-  else
-    for (unsigned int c=0; c<n_components; ++c)
-      this->values_dofs[c][dofs_per_cell-1] = VectorizedArray<Number>();
-
-  internal::do_integrate (*this, this->values_dofs, this->values_quad,
-                          this->gradients_quad, integrate_val, integrate_grad,
-                          internal::int2type<dim>());
+  // Select algorithm matching the element type at run time (the function
+  // pointer is easy to predict, so negligible in cost)
+  integrate_funct (*this->data, this->values_dofs, this->values_quad,
+                   this->gradients_quad, integrate_val, integrate_grad);
 
 #ifdef DEBUG
   this->dof_values_initialized = true;
index ac5670f928c4f72ad91b6754804abc04bb58a344..d246ae4cdb5981aea91ea802e790d70ebbad8e28 100644 (file)
@@ -86,14 +86,9 @@ namespace MatrixFreeOperators
     const FEEvaluationBase<dim,n_components,Number> &fe_eval;
 
     /**
-     * A structure to hold inverse shape functions, put into an aligned vector
-     * in order to make sure it can be allocated.
+     * A structure to hold inverse shape functions
      */
-    struct Inverse1DShape
-    {
-      VectorizedArray<Number> evenodd[fe_degree+1][fe_degree/2+1];
-    };
-    AlignedVector<Inverse1DShape> inverse_shape;
+    AlignedVector<VectorizedArray<Number> > inverse_shape;
   };
 
 
@@ -112,18 +107,19 @@ namespace MatrixFreeOperators
       for (unsigned int j=0; j<shapes_1d.n(); ++j, ++c)
         shapes_1d(i,j) = fe_eval.get_shape_info().shape_values_number[c];
     shapes_1d.gauss_jordan();
-    inverse_shape.resize(1);
-    for (int i=0; i<(fe_degree+1)/2; ++i)
+    const unsigned int stride = (fe_degree+2)/2;
+    inverse_shape.resize(stride*(fe_degree+1));
+    for (unsigned int i=0; i<stride; ++i)
       for (unsigned int q=0; q<(fe_degree+2)/2; ++q)
         {
-          inverse_shape[0].evenodd[i][q] =
+          inverse_shape[i*stride+q] =
             0.5 * (shapes_1d(i,q) + shapes_1d(i,fe_degree-q));
-          inverse_shape[0].evenodd[fe_degree-i][q] =
+          inverse_shape[(fe_degree-i)*stride+q] =
             0.5 * (shapes_1d(i,q) - shapes_1d(i,fe_degree-q));
         }
     if (fe_degree % 2 == 0)
       for (unsigned int q=0; q<(fe_degree+2)/2; ++q)
-        inverse_shape[0].evenodd[fe_degree/2][q] = shapes_1d(fe_degree/2,q);
+        inverse_shape[fe_degree/2*stride+q] = shapes_1d(fe_degree/2,q);
   }
 
 
@@ -175,6 +171,10 @@ namespace MatrixFreeOperators
 
     Assert(dim == 2 || dim == 3, ExcNotImplemented());
 
+    internal::EvaluatorTensorProduct<internal::evaluate_evenodd,dim,fe_degree,
+                                     fe_degree+1, VectorizedArray<Number> >
+      evaluator(inverse_shape, inverse_shape, inverse_shape);
+
     const unsigned int shift_coefficient =
       inverse_coefficients.size() > dofs_per_cell ? dofs_per_cell : 0;
     const VectorizedArray<Number> *inv_coefficient = &inverse_coefficients[0];
@@ -183,34 +183,24 @@ namespace MatrixFreeOperators
       {
         const VectorizedArray<Number>* in = in_array+d*dofs_per_cell;
         VectorizedArray<Number>* out = out_array+d*dofs_per_cell;
-        internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
-                                               VectorizedArray<Number>,0,false,false,2>
-          (inverse_shape[0].evenodd, in, temp_data_field);
-        internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
-                                               VectorizedArray<Number>,1,false,false,2>
-          (inverse_shape[0].evenodd, temp_data_field, out);
+        // Need to select 'apply' method with hessian slot because values
+        // assume symmetries that do not exist in the inverse shapes
+        evaluator.template hessians<0,false,false> (in, temp_data_field);
+        evaluator.template hessians<1,false,false> (temp_data_field, out);
 
         if (dim == 3)
           {
-            internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
-                                                   VectorizedArray<Number>,2,false,false,2>
-              (inverse_shape[0].evenodd, out, temp_data_field);
+            evaluator.template hessians<2,false,false> (out, temp_data_field);
             for (unsigned int q=0; q<dofs_per_cell; ++q)
               temp_data_field[q] *= inv_coefficient[q];
-            internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
-                                                   VectorizedArray<Number>,2,true,false,2>
-              (inverse_shape[0].evenodd, temp_data_field, out);
+            evaluator.template hessians<2,true,false> (temp_data_field, out);
           }
         else if (dim == 2)
           for (unsigned int q=0; q<dofs_per_cell; ++q)
             out[q] *= inv_coefficient[q];
 
-        internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
-                                               VectorizedArray<Number>,1,true,false,2>
-          (inverse_shape[0].evenodd, out, temp_data_field);
-        internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
-                                               VectorizedArray<Number>,0,true,false,2>
-          (inverse_shape[0].evenodd, temp_data_field, out);
+        evaluator.template hessians<1,true,false>(out, temp_data_field);
+        evaluator.template hessians<0,true,false>(temp_data_field, out);
 
         inv_coefficient += shift_coefficient;
       }
index f9b51e21f6c6f1787e3d337d5e2ed6e972fc8605..01be04f08b15fc608a53b53c38512929187dd20c 100644 (file)
@@ -34,6 +34,20 @@ namespace internal
 {
   namespace MatrixFreeFunctions
   {
+    /**
+     * An enum that encodes the type of element detected during
+     * initialization. FEEvaluation will select the most efficient algorithm
+     * based on the given element type.
+     */
+    enum ElementType
+      {
+        tensor_general,
+        tensor_symmetric,
+        truncated_tensor,
+        tensor_symmetric_plus_dg0,
+        tensor_gausslobatto
+      };
+
     /**
      * The class that stores the shape functions, gradients and Hessians
      * evaluated for a tensor product finite element and tensor product
@@ -76,6 +90,13 @@ namespace internal
        */
       std::size_t memory_consumption () const;
 
+      /**
+       * Encodes the type of element detected at construction. FEEvaluation
+       * will select the most efficient algorithm based on the given element
+       * type.
+       */
+      ElementType element_type;
+          
       /**
        * Stores the shape values of the 1D finite element evaluated on all 1D
        * quadrature points in vectorized format, i.e., as an array of
@@ -103,6 +124,27 @@ namespace internal
        */
       AlignedVector<VectorizedArray<Number> > shape_hessians;
 
+      /**
+       * Stores the shape values in a different format, namely the so-called
+       * even-odd scheme where the symmetries in shape_values are used for
+       * faster evaluation.
+       */
+      AlignedVector<VectorizedArray<Number> > shape_val_evenodd;
+
+      /**
+       * Stores the shape gradients in a different format, namely the
+       * so-called even-odd scheme where the symmetries in shape_gradients are
+       * used for faster evaluation.
+       */
+      AlignedVector<VectorizedArray<Number> > shape_gra_evenodd;
+
+      /**
+       * Stores the shape second derivatives in a different format, namely the
+       * so-called even-odd scheme where the symmetries in shape_hessians are
+       * used for faster evaluation.
+       */
+      AlignedVector<VectorizedArray<Number> > shape_hes_evenodd;
+
       /**
        * Stores the indices from cell DoFs to face DoFs. The rows go through
        * the <tt>2*dim</tt> faces, and the columns the DoFs on the faces.
@@ -173,6 +215,19 @@ namespace internal
        * Stores the number of DoFs per face in @p dim dimensions.
        */
       unsigned int dofs_per_face;
+
+      /**
+       * Checks whether we have symmetries in the shape values. In that case,
+       * also fill the shape_???_evenodd fields.
+       */
+      bool check_1d_shapes_symmetric(const unsigned int n_q_points_1d);
+
+      /**
+       * Checks whether symmetric 1D basis functions are such that the shape
+       * values form a diagonal matrix, which allows to use specialized
+       * algorithms that save some operations.
+       */
+      bool check_1d_shapes_gausslobatto();
     };
 
 
index a38ea4fe72b60ad99ed2867417819e6ad8433bfc..82e6566cf9719f5f345466888e9090ae44b03331 100644 (file)
@@ -40,6 +40,7 @@ namespace internal
     template <typename Number>
     ShapeInfo<Number>::ShapeInfo ()
       :
+      element_type (tensor_general),
       n_q_points (0),
       dofs_per_cell (0)
     {}
@@ -85,6 +86,7 @@ namespace internal
 
         const FE_Q_DG0<dim> *fe_q_dg0 = dynamic_cast<const FE_Q_DG0<dim>*>(fe);
 
+        element_type = tensor_general;
         if (fe_poly != 0)
           scalar_lexicographic = fe_poly->get_poly_space_numbering_inverse();
         else if (fe_poly_piece != 0)
@@ -94,9 +96,13 @@ namespace internal
             scalar_lexicographic.resize(fe_dgp->dofs_per_cell);
             for (unsigned int i=0; i<fe_dgp->dofs_per_cell; ++i)
               scalar_lexicographic[i] = i;
+            element_type = truncated_tensor;
           }
         else if (fe_q_dg0 != 0)
-          scalar_lexicographic = fe_q_dg0->get_poly_space_numbering_inverse();
+          {
+            scalar_lexicographic = fe_q_dg0->get_poly_space_numbering_inverse();
+            element_type = tensor_symmetric_plus_dg0;
+          }
         else
           Assert(false, ExcNotImplemented());
 
@@ -184,6 +190,17 @@ namespace internal
           this->face_gradient[1][i] = fe->shape_grad(my_i,q_point)[0];
         }
 
+      if (element_type == tensor_general &&
+          check_1d_shapes_symmetric(n_q_points_1d))
+        {
+          if (check_1d_shapes_gausslobatto())
+            element_type = tensor_gausslobatto;
+          else
+            element_type = tensor_symmetric;
+        }
+      else if (element_type == tensor_symmetric_plus_dg0)
+        check_1d_shapes_symmetric(n_q_points_1d);
+
       // face information
       unsigned int n_faces = GeometryInfo<dim>::faces_per_cell;
       this->face_indices.reinit(n_faces, this->dofs_per_face);
@@ -229,6 +246,137 @@ namespace internal
 
 
 
+    template <typename Number>
+    bool
+    ShapeInfo<Number>::check_1d_shapes_symmetric(const unsigned int n_q_points_1d)
+    {
+      const double zero_tol =
+        types_are_equal<Number,double>::value==true?1e-10:1e-7;
+      // symmetry for values
+      const unsigned int n_dofs_1d = fe_degree + 1;
+      for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
+        for (unsigned int j=0; j<n_q_points_1d; ++j)
+          if (std::fabs(shape_values[i*n_q_points_1d+j][0] -
+                        shape_values[(n_dofs_1d-i)*n_q_points_1d
+                                                 -j-1][0]) > zero_tol)
+            return false;
+
+      // shape values should be zero at x=0.5 for all basis functions except
+      // for one which is one
+      if (n_q_points_1d%2 == 1 && n_dofs_1d%2 == 1)
+        {
+          for (unsigned int i=0; i<n_dofs_1d/2; ++i)
+            if (std::fabs(shape_values[i*n_q_points_1d+
+                                                   n_q_points_1d/2][0]) > zero_tol)
+              return false;
+          if (std::fabs(shape_values[(n_dofs_1d/2)*n_q_points_1d+
+                                                 n_q_points_1d/2][0]-1.)> zero_tol)
+            return false;
+        }
+
+      // skew-symmetry for gradient, zero of middle basis function in middle
+      // quadrature point
+      for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
+        for (unsigned int j=0; j<n_q_points_1d; ++j)
+          if (std::fabs(shape_gradients[i*n_q_points_1d+j][0] +
+                        shape_gradients[(n_dofs_1d-i)*n_q_points_1d-
+                                                    j-1][0]) > zero_tol)
+            return false;
+      if (n_dofs_1d%2 == 1 && n_q_points_1d%2 == 1)
+        if (std::fabs(shape_gradients[(n_dofs_1d/2)*n_q_points_1d+
+                                                  (n_q_points_1d/2)][0]) > zero_tol)
+          return false;
+
+      // symmetry for Laplacian
+      for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
+        for (unsigned int j=0; j<n_q_points_1d; ++j)
+          if (std::fabs(shape_hessians[i*n_q_points_1d+j][0] -
+                        shape_hessians[(n_dofs_1d-i)*n_q_points_1d-
+                                                   j-1][0]) > zero_tol)
+            return false;
+
+      const unsigned int stride = (n_q_points_1d+1)/2;
+      shape_val_evenodd.resize((fe_degree+1)*stride);
+      shape_gra_evenodd.resize((fe_degree+1)*stride);
+      shape_hes_evenodd.resize((fe_degree+1)*stride);
+      for (unsigned int i=0; i<(fe_degree+1)/2; ++i)
+        for (unsigned int q=0; q<stride; ++q)
+          {
+            shape_val_evenodd[i*stride+q] =
+              0.5 * (shape_values[i*n_q_points_1d+q] +
+                     shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
+            shape_val_evenodd[(fe_degree-i)*stride+q] =
+              0.5 * (shape_values[i*n_q_points_1d+q] -
+                     shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
+
+            shape_gra_evenodd[i*stride+q] =
+              0.5 * (shape_gradients[i*n_q_points_1d+q] +
+                     shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
+            shape_gra_evenodd[(fe_degree-i)*stride+q] =
+              0.5 * (shape_gradients[i*n_q_points_1d+q] -
+                     shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
+
+            shape_hes_evenodd[i*stride+q] =
+              0.5 * (shape_hessians[i*n_q_points_1d+q] +
+                     shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
+            shape_hes_evenodd[(fe_degree-i)*stride+q] =
+              0.5 * (shape_hessians[i*n_q_points_1d+q] -
+                     shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
+          }
+      if (fe_degree % 2 == 0)
+        for (unsigned int q=0; q<stride; ++q)
+          {
+            shape_val_evenodd[fe_degree/2*stride+q] =
+              shape_values[(fe_degree/2)*n_q_points_1d+q];
+            shape_gra_evenodd[fe_degree/2*stride+q] =
+              shape_gradients[(fe_degree/2)*n_q_points_1d+q];
+            shape_hes_evenodd[fe_degree/2*stride+q] =
+              shape_hessians[(fe_degree/2)*n_q_points_1d+q];
+          }
+
+      return true;
+    }
+
+
+
+    template <typename Number>
+    bool
+    ShapeInfo<Number>::check_1d_shapes_gausslobatto()
+    {
+      if (dofs_per_cell != n_q_points)
+        return false;
+
+      const double zero_tol =
+        types_are_equal<Number,double>::value==true?1e-10:1e-7;
+      // check: - identity operation for shape values
+      //        - zero diagonal at interior points for gradients
+      //        - gradient equal to unity at element boundary
+      const unsigned int n_points_1d = fe_degree+1;
+      for (unsigned int i=0; i<n_points_1d; ++i)
+        for (unsigned int j=0; j<n_points_1d; ++j)
+          if (i!=j)
+            {
+              if (std::fabs(shape_values[i*n_points_1d+j][0])>zero_tol)
+                return false;
+            }
+          else
+            {
+              if (std::fabs(shape_values[i*n_points_1d+
+                                                     j][0]-1.)>zero_tol)
+                return false;
+            }
+      for (unsigned int i=1; i<n_points_1d-1; ++i)
+        if (std::fabs(shape_gradients[i*n_points_1d+i][0])>zero_tol)
+          return false;
+      if (std::fabs(shape_gradients[n_points_1d-1][0]-
+                    (n_points_1d%2==0 ? -1. : 1.)) > zero_tol)
+        return false;
+
+      return true;
+    }
+
+
+
     template <typename Number>
     std::size_t
     ShapeInfo<Number>::memory_consumption () const
@@ -236,7 +384,10 @@ namespace internal
       std::size_t memory = sizeof(*this);
       memory += MemoryConsumption::memory_consumption(shape_values);
       memory += MemoryConsumption::memory_consumption(shape_gradients);
-      memory += MemoryConsumption::memory_consumption(shape_hessians);
+      memory += MemoryConsumption::memory_consumption(shape_hessians); 
+      memory += MemoryConsumption::memory_consumption(shape_val_evenodd);
+      memory += MemoryConsumption::memory_consumption(shape_gra_evenodd);
+      memory += MemoryConsumption::memory_consumption(shape_hes_evenodd);
       memory += face_indices.memory_consumption();
       for (unsigned int i=0; i<2; ++i)
         {
index dfd39b0910d9c574dd5504634b72c5b72b5ed79d..bea4f8ac5929f6cb27de624b2118967abb554523 100644 (file)
@@ -45,7 +45,7 @@ mass_operator (const MatrixFree<dim,Number>  &data,
                const Vector<Number> &src,
                const std::pair<unsigned int,unsigned int> &cell_range)
 {
-  FEEvaluationGeneral<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
+  FEEvaluation<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
   const unsigned int n_q_points = fe_eval.n_q_points;
 
   for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
index 44cc266fe552fdd28b0e7a914b27aeb07e12f0ab..56c0787dbfaee05e8e7f058a1dfec1610e67d373 100644 (file)
@@ -32,24 +32,24 @@ template <int M, int N, int type, bool add>
 void test()
 {
   deallog << "Test " << M << " x " << N << std::endl;
-  double shape[M][N];
+  AlignedVector<double> shape(M*N);
   for (unsigned int i=0; i<(M+1)/2; ++i)
     for (unsigned int j=0; j<N; ++j)
       {
-        shape[i][j] = -1. + 2. * (double)Testing::rand()/RAND_MAX;
+        shape[i*N+j] = -1. + 2. * (double)Testing::rand()/RAND_MAX;
         if (type == 1)
-          shape[M-1-i][N-1-j] = -shape[i][j];
+          shape[(M-1-i)*N+N-1-j] = -shape[i*N+j];
         else
-          shape[M-1-i][N-1-j] = shape[i][j];
+          shape[(M-1-i)*N+N-1-j] = shape[i*N+j];
       }
   if (type == 0 && M%2 == 1 && N%2 == 1)
     {
       for (unsigned int i=0; i<M; ++i)
-        shape[i][N/2] = 0.;
-      shape[M/2][N/2] = 1;
+        shape[i*N+N/2] = 0.;
+      shape[M/2*N+N/2] = 1;
     }
   if (type == 1 && M%2 == 1 && N%2 == 1)
-    shape[M/2][N/2] = 0.;
+    shape[M/2*N+N/2] = 0.;
 
   double x[N], x_ref[N], y[M], y_ref[M];
   for (unsigned int i=0; i<N; ++i)
@@ -61,19 +61,17 @@ void test()
       y[i] = 1.;
       y_ref[i] = add ? y[i] : 0.;
       for (unsigned int j=0; j<N; ++j)
-        y_ref[i] += shape[i][j] * x[j];
+        y_ref[i] += shape[i*N+j] * x[j];
     }
 
   // apply function for tensor product
+  internal::EvaluatorTensorProduct<internal::evaluate_symmetric,1,M-1,N,double> evaluator(shape, shape, shape);
   if (type == 0)
-    internal::apply_tensor_product_values<1,M-1,N,double,0,false,add>
-      (&shape[0][0],x,y);
+    evaluator.template values<0,false,add> (x,y);
   if (type == 1)
-    internal::apply_tensor_product_gradients<1,M-1,N,double,0,false,add>
-      (&shape[0][0],x,y);
+    evaluator.template gradients<0,false,add> (x,y);
   if (type == 2)
-    internal::apply_tensor_product_hessians<1,M-1,N,double,0,false,add>
-      (&shape[0][0],x,y);
+    evaluator.template hessians<0,false,add> (x,y);
 
   deallog << "Errors no transpose: ";
   for (unsigned int i=0; i<M; ++i)
@@ -90,19 +88,16 @@ void test()
       x[i] = 2.;
       x_ref[i] = add ? x[i] : 0.;
       for (unsigned int j=0; j<M; ++j)
-        x_ref[i] += shape[j][i] * y[j];
+        x_ref[i] += shape[j*N+i] * y[j];
     }
 
   // apply function for tensor product
   if (type == 0)
-    internal::apply_tensor_product_values<1,M-1,N,double,0,true,add>
-      (&shape[0][0],y,x);
+    evaluator.template values<0,true,add> (y,x);
   if (type == 1)
-    internal::apply_tensor_product_gradients<1,M-1,N,double,0,true,add>
-      (&shape[0][0],y,x);
+    evaluator.template gradients<0,true,add> (y,x);
   if (type == 2)
-    internal::apply_tensor_product_hessians<1,M-1,N,double,0,true,add>
-      (&shape[0][0],y,x);
+    evaluator.template hessians<0,true,add> (y,x);
 
   deallog << "Errors transpose:    ";
   for (unsigned int i=0; i<N; ++i)
index 4db1704541c6507004f1d37d50d804b15b6185ff..a0eebb553a3e8ec795384e2edcddada1d1f95445 100644 (file)
@@ -54,16 +54,16 @@ void test()
 
   // create symmetrized shape array exactly as expected by the evenodd
   // function
-  double shape_sym[M][(N+1)/2];
+  AlignedVector<double> shape_sym(M*((N+1)/2));
   for (unsigned int i=0; i<M/2; ++i)
     for (unsigned int q=0; q<(N+1)/2; ++q)
       {
-        shape_sym[i][q] = 0.5 * (shape[i][q] + shape[i][N-1-q]);
-        shape_sym[M-1-i][q] = 0.5 * (shape[i][q] - shape[i][N-1-q]);
+        shape_sym[i*((N+1)/2)+q] = 0.5 * (shape[i][q] + shape[i][N-1-q]);
+        shape_sym[(M-1-i)*((N+1)/2)+q] = 0.5 * (shape[i][q] - shape[i][N-1-q]);
       }
   if (M % 2 == 1)
     for (unsigned int q=0; q<(N+1)/2; ++q)
-      shape_sym[(M-1)/2][q] = shape[(M-1)/2][q];
+      shape_sym[(M-1)/2*((N+1)/2)+q] = shape[(M-1)/2][q];
 
   double x[N], x_ref[N], y[M], y_ref[M];
   for (unsigned int i=0; i<N; ++i)
@@ -77,9 +77,16 @@ void test()
       for (unsigned int j=0; j<N; ++j)
         y_ref[i] += shape[i][j] * x[j];
     }
-
+                
   // apply function for tensor product
-  internal::apply_tensor_product_evenodd<1,M-1,N,double,0,false,add,type>(shape_sym,x,y);
+  internal::EvaluatorTensorProduct<internal::evaluate_evenodd,1,M-1,N,double> evaluator(shape_sym, shape_sym, shape_sym);
+  if (type == 0)
+    evaluator.template values<0,false,add> (x,y);
+  if (type == 1)
+    evaluator.template gradients<0,false,add> (x,y);
+  if (type == 2)
+    evaluator.template hessians<0,false,add> (x,y);
+
 
   deallog << "Errors no transpose: ";
   for (unsigned int i=0; i<M; ++i)
@@ -100,7 +107,12 @@ void test()
     }
 
   // apply function for tensor product
-  internal::apply_tensor_product_evenodd<1,M-1,N,double,0,true,add,type>(shape_sym,y,x);
+  if (type == 0)
+    evaluator.template values<0,true,add> (y,x);
+  if (type == 1)
+    evaluator.template gradients<0,true,add> (y,x);
+  if (type == 2)
+    evaluator.template hessians<0,true,add> (y,x);
 
   deallog << "Errors transpose:    ";
   for (unsigned int i=0; i<N; ++i)
index c585c8e235ab7509df5e1a8478acc2674ad08f48..34a3799c49dde1c4b407c1ba215f1ba6f84f510b 100644 (file)
@@ -30,60 +30,6 @@ std::ofstream logfile("output");
 #include "get_functions_common.h"
 
 
-template <int dim, int fe_degree, typename Number>
-class MatrixFreeTestGL : public MatrixFreeTest<dim, fe_degree, fe_degree+1, Number>
-{
-public:
-
-  MatrixFreeTestGL(const MatrixFree<dim,Number> &data,
-                   const Mapping<dim>               &mapping):
-    MatrixFreeTest<dim, fe_degree, fe_degree+1, Number>(data, mapping)
-  {};
-
-  void operator() (const MatrixFree<dim,Number> &data,
-                   Vector<Number> &,
-                   const Vector<Number> &src,
-                   const std::pair<unsigned int,unsigned int> &cell_range) const
-  {
-    FEEvaluationGL<dim,fe_degree,1,Number> fe_eval (this->data);
-    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
-      {
-        fe_eval.reinit (cell);
-        std::vector<double> reference_values (fe_eval.n_q_points);
-        std::vector<Tensor<1,dim> > reference_grads (fe_eval.n_q_points);
-        std::vector<Tensor<2,dim> > reference_hess (fe_eval.n_q_points);
-        fe_eval.read_dof_values(src);
-        fe_eval.evaluate (true,true,true);
-
-        // compare values with the ones the FEValues
-        // gives us. Those are seen as reference
-        for (unsigned int j=0; j<data.n_components_filled(cell); ++j)
-          {
-            this->fe_val.reinit (data.get_cell_iterator(cell,j));
-            this->fe_val.get_function_values(src, reference_values);
-            this->fe_val.get_function_gradients(src, reference_grads);
-            this->fe_val.get_function_hessians(src, reference_hess);
-
-            for (int q=0; q<(int)fe_eval.n_q_points; q++)
-              {
-                this->errors[0] += std::fabs(fe_eval.get_value(q)[j]-
-                                             reference_values[q]);
-                for (unsigned int d=0; d<dim; ++d)
-                  this->errors[1] += std::fabs(fe_eval.get_gradient(q)[d][j]-
-                                               reference_grads[q][d]);
-                this->errors[2] += std::fabs(fe_eval.get_laplacian(q)[j]-
-                                             trace(reference_hess[q]));
-                this->total[0] += std::fabs(reference_values[q]);
-                for (unsigned int d=0; d<dim; ++d)
-                  this->total[1] += std::fabs(reference_grads[q][d]);
-                this->total[2] += std::fabs(fe_eval.get_laplacian(q)[j]);
-              }
-          }
-      }
-  }
-};
-
-
 
 template <int dim, int fe_degree>
 void test ()
@@ -134,7 +80,7 @@ void test ()
   data.tasks_parallel_scheme = MatrixFree<dim,number>::AdditionalData::none;
   data.mapping_update_flags = update_gradients | update_second_derivatives;
   mf_data.reinit (mapping, dof, constraints, quad, data);
-  MatrixFreeTestGL<dim,fe_degree,number> mf (mf_data, mapping);
+  MatrixFreeTest<dim,fe_degree,fe_degree+1,number> mf (mf_data, mapping);
   mf.test_functions (solution);
 }
 
index ce3c796e2093548d8cf5e3693170c95915dff11f..6e430c9dc7e0ea4ef522c83d134fca7319bc9d6d 100644 (file)
@@ -32,61 +32,6 @@ std::ofstream logfile("output");
 #include "get_functions_common.h"
 
 
-template <int dim, int fe_degree, typename Number>
-class MatrixFreeTestGen : public MatrixFreeTest<dim, fe_degree, fe_degree+1, Number>
-{
-public:
-  MatrixFreeTestGen(const MatrixFree<dim,Number> &data,
-                    const Mapping<dim>               &mapping):
-    MatrixFreeTest<dim, fe_degree, fe_degree+1, Number>(data, mapping)
-  {};
-
-  void operator () (const MatrixFree<dim,Number> &data,
-                    Vector<Number> &,
-                    const Vector<Number> &src,
-                    const std::pair<unsigned int,unsigned int> &cell_range) const
-  {
-    FEEvaluationGeneral<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
-    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
-      {
-        fe_eval.reinit (cell);
-        std::vector<double> reference_values (fe_eval.n_q_points);
-        std::vector<Tensor<1,dim> > reference_grads (fe_eval.n_q_points);
-        std::vector<Tensor<2,dim> > reference_hess (fe_eval.n_q_points);
-        fe_eval.read_dof_values(src);
-        fe_eval.evaluate (true,true,true);
-
-        // compare values with the ones the FEValues
-        // gives us. Those are seen as reference
-        for (unsigned int j=0; j<data.n_components_filled(cell); ++j)
-          {
-            this->fe_val.reinit (data.get_cell_iterator(cell,j));
-            this->fe_val.get_function_values(src, reference_values);
-            this->fe_val.get_function_gradients(src, reference_grads);
-            this->fe_val.get_function_hessians(src, reference_hess);
-
-            for (int q=0; q<(int)fe_eval.n_q_points; q++)
-              {
-                this->errors[0] += std::fabs(fe_eval.get_value(q)[j]-
-                                             reference_values[q]);
-                for (unsigned int d=0; d<dim; ++d)
-                  {
-                    this->errors[1] += std::fabs(fe_eval.get_gradient(q)[d][j]-
-                                                 reference_grads[q][d]);
-                  }
-                this->errors[2] += std::fabs(fe_eval.get_laplacian(q)[j]-
-                                             trace(reference_hess[q]));
-                this->total[0] += std::fabs(reference_values[q]);
-                for (unsigned int d=0; d<dim; ++d)
-                  this->total[1] += std::fabs(reference_grads[q][d]);
-                this->total[2] += std::fabs(fe_eval.get_laplacian(q)[j]);
-              }
-          }
-      }
-  }
-};
-
-
 
 template <int dim, int fe_degree>
 void test ()
@@ -137,7 +82,7 @@ void test ()
   data.tasks_parallel_scheme = MatrixFree<dim,number>::AdditionalData::none;
   data.mapping_update_flags = update_gradients | update_second_derivatives;
   mf_data.reinit (mapping, dof, constraints, quad, data);
-  MatrixFreeTestGen<dim,fe_degree,number> mf (mf_data, mapping);
+  MatrixFreeTest<dim,fe_degree,fe_degree+1,number> mf (mf_data, mapping);
   mf.test_functions (solution);
 }
 
index 80126fc3ea8315a3e038629304acc3c6bb285390..6d18bd6199ee8f7d1b3ae74a180b1d0783929539 100644 (file)
 // approximation of the Laplacian. It only contains cell terms.
 
 #include "../tests.h"
-#include <deal.II/base/function.h>
 #include <deal.II/fe/fe_dgp.h>
-#include <deal.II/matrix_free/matrix_free.h>
-#include <deal.II/matrix_free/fe_evaluation.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/lac/constraint_matrix.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/numerics/vector_tools.h>
-
-#include <deal.II/lac/vector.h>
 
 
 std::ofstream logfile("output");
 
-
-
-template <int dim, int fe_degree, typename VECTOR>
-void
-helmholtz_operator_dgp (const MatrixFree<dim,typename VECTOR::value_type>  &data,
-                        VECTOR       &dst,
-                        const VECTOR &src,
-                        const std::pair<unsigned int,unsigned int> &cell_range)
-{
-  typedef typename VECTOR::value_type Number;
-  FEEvaluationDGP<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
-  const unsigned int n_q_points = fe_eval.n_q_points;
-
-  for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
-    {
-      fe_eval.reinit (cell);
-      fe_eval.read_dof_values (src);
-      fe_eval.evaluate (true, true, false);
-      for (unsigned int q=0; q<n_q_points; ++q)
-        {
-          fe_eval.submit_value (Number(10)*fe_eval.get_value(q),q);
-          fe_eval.submit_gradient (fe_eval.get_gradient(q),q);
-        }
-      fe_eval.integrate (true,true);
-      fe_eval.distribute_local_to_global (dst);
-    }
-}
-
-
-
-template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
-class MatrixFreeTest
-{
-public:
-  typedef VectorizedArray<Number> vector_t;
-
-  MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
-    data (data_in)
-  {};
-
-  void vmult (VECTOR       &dst,
-              const VECTOR &src) const
-  {
-    dst = 0;
-    const std_cxx1x::function<void(const MatrixFree<dim,typename VECTOR::value_type> &,
-                                   VECTOR &,
-                                   const VECTOR &,
-                                   const std::pair<unsigned int,unsigned int> &)>
-      wrap = helmholtz_operator_dgp<dim,fe_degree,VECTOR>;
-    data.cell_loop (wrap, dst, src);
-  };
-
-private:
-  const MatrixFree<dim,Number> &data;
-};
-
-
-
-template <int dim, int fe_degree, typename number>
-void do_test (const DoFHandler<dim> &dof,
-              const ConstraintMatrix &constraints)
-{
-
-  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
-
-  MatrixFree<dim,number> mf_data;
-  {
-    const QGauss<1> quad (fe_degree+1);
-    typename MatrixFree<dim,number>::AdditionalData data;
-    data.tasks_parallel_scheme =
-      MatrixFree<dim,number>::AdditionalData::partition_color;
-    data.tasks_block_size = 3;
-
-    mf_data.reinit (dof, constraints, quad, data);
-  }
-
-  MatrixFreeTest<dim,fe_degree,number> mf (mf_data);
-  Vector<number> in (dof.n_dofs()), out (dof.n_dofs());
-  Vector<number> in_dist (dof.n_dofs());
-  Vector<number> out_dist (in_dist);
-
-  for (unsigned int i=0; i<dof.n_dofs(); ++i)
-    {
-      if (constraints.is_constrained(i))
-        continue;
-      const double entry = Testing::rand()/(double)RAND_MAX;
-      in(i) = entry;
-      in_dist(i) = entry;
-    }
-
-  mf.vmult (out_dist, in_dist);
-
-
-  // assemble sparse matrix with (\nabla v, \nabla u) + (v, 10 * u)
-  SparsityPattern sparsity;
-  {
-    CompressedSimpleSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
-    DoFTools::make_sparsity_pattern (dof, csp, constraints, true);
-    sparsity.copy_from(csp);
-  }
-  SparseMatrix<double> sparse_matrix (sparsity);
-  {
-    QGauss<dim>  quadrature_formula(fe_degree+1);
-
-    FEValues<dim> fe_values (dof.get_fe(), quadrature_formula,
-                             update_values    |  update_gradients |
-                             update_JxW_values);
-
-    const unsigned int   dofs_per_cell = dof.get_fe().dofs_per_cell;
-    const unsigned int   n_q_points    = quadrature_formula.size();
-
-    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
-    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
-
-    typename DoFHandler<dim>::active_cell_iterator
-    cell = dof.begin_active(),
-    endc = dof.end();
-    for (; cell!=endc; ++cell)
-      {
-        cell_matrix = 0;
-        fe_values.reinit (cell);
-
-        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-          for (unsigned int i=0; i<dofs_per_cell; ++i)
-            {
-              for (unsigned int j=0; j<dofs_per_cell; ++j)
-                cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point) *
-                                      fe_values.shape_grad(j,q_point)
-                                      +
-                                      10. *
-                                      fe_values.shape_value(i,q_point) *
-                                      fe_values.shape_value(j,q_point)) *
-                                     fe_values.JxW(q_point));
-            }
-
-        cell->get_dof_indices(local_dof_indices);
-        constraints.distribute_local_to_global (cell_matrix,
-                                                local_dof_indices,
-                                                sparse_matrix);
-      }
-  }
-
-  sparse_matrix.vmult (out, in);
-  out -= out_dist;
-  const double diff_norm = out.linfty_norm() / out_dist.linfty_norm();
-
-  deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
-}
+#include "matrix_vector_common.h"
 
 
 
@@ -223,26 +60,3 @@ void test ()
   do_test<dim, fe_degree, double> (dof, constraints);
 }
 
-
-
-int main ()
-{
-  deallog.attach(logfile);
-  deallog.depth_console(0);
-
-  deallog << std::setprecision (3);
-
-  {
-    deallog.threshold_double(5.e-11);
-    deallog.push("2d");
-    test<2,1>();
-    test<2,2>();
-    deallog.pop();
-    deallog.push("3d");
-    test<3,1>();
-    test<3,2>();
-    deallog.pop();
-  }
-}
-
-
index 04e33dd37f4418cee0f2484dd76a59c4939246b7..65c44ec4b7d7e1fe5406c1d99536ae9e86a400f0 100644 (file)
 // degree would suggest
 
 #include "../tests.h"
-#include <deal.II/base/function.h>
 #include <deal.II/fe/fe_dgp.h>
-#include <deal.II/matrix_free/matrix_free.h>
-#include <deal.II/matrix_free/fe_evaluation.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/lac/constraint_matrix.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/numerics/vector_tools.h>
-
-#include <deal.II/lac/vector.h>
-
 
 std::ofstream logfile("output");
 
-
-
-template <int dim, int fe_degree, typename VECTOR>
-void
-helmholtz_operator_dgp (const MatrixFree<dim,typename VECTOR::value_type>  &data,
-                        VECTOR       &dst,
-                        const VECTOR &src,
-                        const std::pair<unsigned int,unsigned int> &cell_range)
-{
-  typedef typename VECTOR::value_type Number;
-  FEEvaluationDGP<dim,fe_degree,fe_degree+2,1,Number> fe_eval (data);
-  const unsigned int n_q_points = fe_eval.n_q_points;
-
-  for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
-    {
-      fe_eval.reinit (cell);
-      fe_eval.read_dof_values (src);
-      fe_eval.evaluate (true, true, false);
-      for (unsigned int q=0; q<n_q_points; ++q)
-        {
-          fe_eval.submit_value (Number(10)*fe_eval.get_value(q),q);
-          fe_eval.submit_gradient (fe_eval.get_gradient(q),q);
-        }
-      fe_eval.integrate (true,true);
-      fe_eval.distribute_local_to_global (dst);
-    }
-}
-
-
-
-template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
-class MatrixFreeTest
-{
-public:
-  typedef VectorizedArray<Number> vector_t;
-
-  MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
-    data (data_in)
-  {};
-
-  void vmult (VECTOR       &dst,
-              const VECTOR &src) const
-  {
-    dst = 0;
-    const std_cxx1x::function<void(const MatrixFree<dim,typename VECTOR::value_type> &,
-                                   VECTOR &,
-                                   const VECTOR &,
-                                   const std::pair<unsigned int,unsigned int> &)>
-      wrap = helmholtz_operator_dgp<dim,fe_degree,VECTOR>;
-    data.cell_loop (wrap, dst, src);
-  };
-
-private:
-  const MatrixFree<dim,Number> &data;
-};
-
-
-
-template <int dim, int fe_degree, typename number>
-void do_test (const DoFHandler<dim> &dof,
-              const ConstraintMatrix &constraints)
-{
-
-  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
-
-  MatrixFree<dim,number> mf_data;
-  {
-    const QGauss<1> quad (fe_degree+2);
-    typename MatrixFree<dim,number>::AdditionalData data;
-    data.tasks_parallel_scheme =
-      MatrixFree<dim,number>::AdditionalData::partition_color;
-    data.tasks_block_size = 3;
-
-    mf_data.reinit (dof, constraints, quad, data);
-  }
-
-  MatrixFreeTest<dim,fe_degree,number> mf (mf_data);
-  Vector<number> in (dof.n_dofs()), out (dof.n_dofs());
-  Vector<number> in_dist (dof.n_dofs());
-  Vector<number> out_dist (in_dist);
-
-  for (unsigned int i=0; i<dof.n_dofs(); ++i)
-    {
-      if (constraints.is_constrained(i))
-        continue;
-      const double entry = Testing::rand()/(double)RAND_MAX;
-      in(i) = entry;
-      in_dist(i) = entry;
-    }
-
-  mf.vmult (out_dist, in_dist);
-
-
-  // assemble sparse matrix with (\nabla v, \nabla u) + (v, 10 * u)
-  SparsityPattern sparsity;
-  {
-    CompressedSimpleSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
-    DoFTools::make_sparsity_pattern (dof, csp, constraints, true);
-    sparsity.copy_from(csp);
-  }
-  SparseMatrix<double> sparse_matrix (sparsity);
-  {
-    QGauss<dim>  quadrature_formula(fe_degree+2);
-
-    FEValues<dim> fe_values (dof.get_fe(), quadrature_formula,
-                             update_values    |  update_gradients |
-                             update_JxW_values);
-
-    const unsigned int   dofs_per_cell = dof.get_fe().dofs_per_cell;
-    const unsigned int   n_q_points    = quadrature_formula.size();
-
-    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
-    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
-
-    typename DoFHandler<dim>::active_cell_iterator
-    cell = dof.begin_active(),
-    endc = dof.end();
-    for (; cell!=endc; ++cell)
-      {
-        cell_matrix = 0;
-        fe_values.reinit (cell);
-
-        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-          for (unsigned int i=0; i<dofs_per_cell; ++i)
-            {
-              for (unsigned int j=0; j<dofs_per_cell; ++j)
-                cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point) *
-                                      fe_values.shape_grad(j,q_point)
-                                      +
-                                      10. *
-                                      fe_values.shape_value(i,q_point) *
-                                      fe_values.shape_value(j,q_point)) *
-                                     fe_values.JxW(q_point));
-            }
-
-        cell->get_dof_indices(local_dof_indices);
-        constraints.distribute_local_to_global (cell_matrix,
-                                                local_dof_indices,
-                                                sparse_matrix);
-      }
-  }
-
-  sparse_matrix.vmult (out, in);
-  out -= out_dist;
-  const double diff_norm = out.linfty_norm() / out_dist.linfty_norm();
-
-  deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
-}
+#include "matrix_vector_common.h"
 
 
 
@@ -220,26 +56,3 @@ void test ()
   do_test<dim, fe_degree, double> (dof, constraints);
 }
 
-
-
-int main ()
-{
-  deallog.attach(logfile);
-  deallog.depth_console(0);
-
-  deallog << std::setprecision (3);
-
-  {
-    deallog.threshold_double(5.e-11);
-    deallog.push("2d");
-    test<2,1>();
-    test<2,2>();
-    deallog.pop();
-    deallog.push("3d");
-    test<3,1>();
-    test<3,2>();
-    deallog.pop();
-  }
-}
-
-
index 4ab84fd05673a61ab0dfd4508cbde5543ba641b5..3fd4d4bf4fc1bf83c1e6a89e947305b0838c5527 100644 (file)
 // in an approximation of the Laplacian. It only contains cell terms.
 
 #include "../tests.h"
-#include <deal.II/base/function.h>
 #include <deal.II/fe/fe_q_dg0.h>
-#include <deal.II/matrix_free/matrix_free.h>
-#include <deal.II/matrix_free/fe_evaluation.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/lac/constraint_matrix.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/numerics/vector_tools.h>
-
-#include <deal.II/lac/vector.h>
-
 
 std::ofstream logfile("output");
 
-
-
-template <int dim, int fe_degree, typename VECTOR>
-void
-helmholtz_operator_qdg0 (const MatrixFree<dim,typename VECTOR::value_type>  &data,
-                         VECTOR       &dst,
-                         const VECTOR &src,
-                         const std::pair<unsigned int,unsigned int> &cell_range)
-{
-  typedef typename VECTOR::value_type Number;
-  FEEvaluationQ_DG0<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
-  const unsigned int n_q_points = fe_eval.n_q_points;
-
-  for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
-    {
-      fe_eval.reinit (cell);
-      fe_eval.read_dof_values (src);
-      fe_eval.evaluate (true, true, false);
-      for (unsigned int q=0; q<n_q_points; ++q)
-        {
-          fe_eval.submit_value (Number(10)*fe_eval.get_value(q),q);
-          fe_eval.submit_gradient (fe_eval.get_gradient(q),q);
-        }
-      fe_eval.integrate (true,true);
-      fe_eval.distribute_local_to_global (dst);
-    }
-}
-
-
-
-template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
-class MatrixFreeTest
-{
-public:
-  typedef VectorizedArray<Number> vector_t;
-
-  MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
-    data (data_in)
-  {};
-
-  void vmult (VECTOR       &dst,
-              const VECTOR &src) const
-  {
-    dst = 0;
-    const std_cxx1x::function<void(const MatrixFree<dim,typename VECTOR::value_type> &,
-                                   VECTOR &,
-                                   const VECTOR &,
-                                   const std::pair<unsigned int,unsigned int> &)>
-      wrap = helmholtz_operator_qdg0<dim,fe_degree,VECTOR>;
-    data.cell_loop (wrap, dst, src);
-  };
-
-private:
-  const MatrixFree<dim,Number> &data;
-};
-
-
-
-template <int dim, int fe_degree, typename number>
-void do_test (const DoFHandler<dim> &dof,
-              const ConstraintMatrix &constraints)
-{
-
-  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
-
-  MatrixFree<dim,number> mf_data;
-  {
-    const QGauss<1> quad (fe_degree+1);
-    typename MatrixFree<dim,number>::AdditionalData data;
-    data.tasks_parallel_scheme =
-      MatrixFree<dim,number>::AdditionalData::partition_color;
-    data.tasks_block_size = 3;
-
-    mf_data.reinit (dof, constraints, quad, data);
-  }
-
-  MatrixFreeTest<dim,fe_degree,number> mf (mf_data);
-  Vector<number> in (dof.n_dofs()), out (dof.n_dofs());
-  Vector<number> in_dist (dof.n_dofs());
-  Vector<number> out_dist (in_dist);
-
-  for (unsigned int i=0; i<dof.n_dofs(); ++i)
-    {
-      if (constraints.is_constrained(i))
-        continue;
-      const double entry = Testing::rand()/(double)RAND_MAX;
-      in(i) = entry;
-      in_dist(i) = entry;
-    }
-
-  mf.vmult (out_dist, in_dist);
-
-
-  // assemble sparse matrix with (\nabla v, \nabla u) + (v, 10 * u)
-  SparsityPattern sparsity;
-  {
-    CompressedSimpleSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
-    DoFTools::make_sparsity_pattern (dof, csp, constraints, true);
-    sparsity.copy_from(csp);
-  }
-  SparseMatrix<double> sparse_matrix (sparsity);
-  {
-    QGauss<dim>  quadrature_formula(fe_degree+1);
-
-    FEValues<dim> fe_values (dof.get_fe(), quadrature_formula,
-                             update_values    |  update_gradients |
-                             update_JxW_values);
-
-    const unsigned int   dofs_per_cell = dof.get_fe().dofs_per_cell;
-    const unsigned int   n_q_points    = quadrature_formula.size();
-
-    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
-    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
-
-    typename DoFHandler<dim>::active_cell_iterator
-    cell = dof.begin_active(),
-    endc = dof.end();
-    for (; cell!=endc; ++cell)
-      {
-        cell_matrix = 0;
-        fe_values.reinit (cell);
-
-        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-          for (unsigned int i=0; i<dofs_per_cell; ++i)
-            {
-              for (unsigned int j=0; j<dofs_per_cell; ++j)
-                cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point) *
-                                      fe_values.shape_grad(j,q_point)
-                                      +
-                                      10. *
-                                      fe_values.shape_value(i,q_point) *
-                                      fe_values.shape_value(j,q_point)) *
-                                     fe_values.JxW(q_point));
-            }
-
-        cell->get_dof_indices(local_dof_indices);
-        constraints.distribute_local_to_global (cell_matrix,
-                                                local_dof_indices,
-                                                sparse_matrix);
-      }
-  }
-
-  sparse_matrix.vmult (out, in);
-  out -= out_dist;
-  const double diff_norm = out.linfty_norm() / out_dist.linfty_norm();
-
-  deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
-}
+#include "matrix_vector_common.h"
 
 
 
@@ -219,30 +55,13 @@ void test ()
   DoFHandler<dim> dof (tria);
   dof.distribute_dofs(fe);
   ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints(dof, constraints);
+  VectorTools::interpolate_boundary_values (dof, 0, ZeroFunction<dim>(),
+                                            constraints);
+  constraints.close();
 
   do_test<dim, fe_degree, double> (dof, constraints);
 }
 
 
 
-int main ()
-{
-  deallog.attach(logfile);
-  deallog.depth_console(0);
-
-  deallog << std::setprecision (3);
-
-  {
-    deallog.threshold_double(5.e-11);
-    deallog.push("2d");
-    test<2,1>();
-    test<2,2>();
-    deallog.pop();
-    deallog.push("3d");
-    test<3,1>();
-    test<3,2>();
-    deallog.pop();
-  }
-}
-
-
index c746da6c998f4b6168b38b9b687464526598d0fc..3ed93664409925403fd74585c693ba480d7d538c 100644 (file)
@@ -72,8 +72,8 @@ public:
                const std::pair<unsigned int,unsigned int> &cell_range) const
   {
     typedef VectorizedArray<Number> vector_t;
-    FEEvaluation     <dim,degree_p+1,degree_p+2,dim,Number> velocity (data, 0);
-    FEEvaluationQ_DG0<dim,degree_p  ,degree_p+2,1,  Number> pressure (data, 1);
+    FEEvaluation<dim,degree_p+1,degree_p+2,dim,Number> velocity (data, 0);
+    FEEvaluation<dim,degree_p  ,degree_p+2,1,  Number> pressure (data, 1);
 
     for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
       {
index b920f8b52c49d19621ff524a07dfc629805d532f..af1625e49ad4d9fb536a379158a88b6f6cd0f0d4 100644 (file)
@@ -92,8 +92,8 @@ namespace Step48
 
     data.initialize_dof_vector (inv_mass_matrix);
 
-    FEEvaluationGL<dim,fe_degree> fe_eval(data);
-    const unsigned int            n_q_points = fe_eval.n_q_points;
+    FEEvaluation<dim,fe_degree> fe_eval(data);
+    const unsigned int          n_q_points = fe_eval.n_q_points;
 
     for (unsigned int cell=0; cell<data.get_size_info().n_macro_cells; ++cell)
       {
@@ -124,7 +124,7 @@ namespace Step48
                const std::pair<unsigned int,unsigned int> &cell_range) const
   {
     AssertDimension (src.size(), 2);
-    FEEvaluationGL<dim,fe_degree> current (data), old (data);
+    FEEvaluation<dim,fe_degree> current (data), old (data);
     for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
       {
         current.reinit (cell);
index 0d90bc53048a57e7780dcb877ab9c82a25a5b7e9..7ee790599e4221484fd942d147ca1697a5478903 100644 (file)
@@ -87,8 +87,8 @@ namespace Step48
 
     data.initialize_dof_vector (inv_mass_matrix);
 
-    FEEvaluationGL<dim,fe_degree> fe_eval(data);
-    const unsigned int            n_q_points = fe_eval.n_q_points;
+    FEEvaluation<dim,fe_degree> fe_eval(data);
+    const unsigned int          n_q_points = fe_eval.n_q_points;
 
     for (unsigned int cell=0; cell<data.get_size_info().n_macro_cells; ++cell)
       {
@@ -119,7 +119,7 @@ namespace Step48
                const std::pair<unsigned int,unsigned int> &cell_range) const
   {
     AssertDimension (src.size(), 2);
-    FEEvaluationGL<dim,fe_degree> current (data), old (data);
+    FEEvaluation<dim,fe_degree> current (data), old (data);
     deallog << "submit / sine values: ";
     for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
       {
index c14b9ebce509c9f50d4e2f8b0f8fac1758f3eceb..66cb180d0b4f849b0f5bc33bf13c5716551ae698 100644 (file)
@@ -92,8 +92,8 @@ namespace Step48
 
     data.initialize_dof_vector (inv_mass_matrix);
 
-    FEEvaluationGL<dim,fe_degree> fe_eval(data);
-    const unsigned int            n_q_points = fe_eval.n_q_points;
+    FEEvaluation<dim,fe_degree> fe_eval(data);
+    const unsigned int          n_q_points = fe_eval.n_q_points;
 
     for (unsigned int cell=0; cell<data.get_size_info().n_macro_cells; ++cell)
       {
@@ -124,7 +124,7 @@ namespace Step48
                const std::pair<unsigned int,unsigned int> &cell_range) const
   {
     AssertDimension (src.size(), 2);
-    FEEvaluationGL<dim,fe_degree> current (data), old (data);
+    FEEvaluation<dim,fe_degree> current (data), old (data);
     for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
       {
         current.reinit (cell);

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.