Improve documentation, make names of data structures uniform

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index 337b5a9618d784821bd8b33513cff91e0956e10a..c5f5af51a3b00565a6176196f10e2f2c2b468f2f 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -70,7 +70,7 @@ namespace internal
   };
 
   template <bool is_long>
-  struct EvaluatorSelector<MatrixFreeFunctions::tensor_gausslobatto,is_long>
+  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_collocation,is_long>
   {
     static const EvaluatorVariant variant = evaluate_evenodd;
   };
@@ -86,11 +86,12 @@ namespace internal
    * the template classes EvaluatorTensorProduct which in turn are selected
    * from the MatrixFreeFunctions::ElementType template argument.
    *
-   * There are two specialized implementation classes FEEvaluationImplSpectral
-   * (for Gauss-Lobatto elements where the 'values' operation is identity) and
-   * FEEvaluationImplTransformSpectral (which can be transformed to a spectral
-   * evaluation and uses the identity in these contexts), which both allow for
-   * shorter code.
+   * There are two specialized implementation classes
+   * FEEvaluationImplCollocation (for Gauss-Lobatto elements where the nodal
+   * points and the quadrature points coincide and the 'values' operation is
+   * identity) and FEEvaluationImplTransformToCollocation (which can be
+   * transformed to a collocation space and can then use the identity in these
+   * spaces), which both allow for shorter code.
    *
    * @author Katharina Kormann, Martin Kronbichler, 2012, 2014, 2017
    */
@@ -105,9 +106,9 @@ namespace internal
                    VectorizedArray<Number> *gradients_quad[][dim],
                    VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
                    VectorizedArray<Number> *scratch_data,
-                   const bool               evaluate_val,
-                   const bool               evaluate_grad,
-                   const bool               evaluate_lapl);
+                   const bool               evaluate_values,
+                   const bool               evaluate_gradients,
+                   const bool               evaluate_hessians);
 
     static
     void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
@@ -115,8 +116,8 @@ namespace internal
                     VectorizedArray<Number> *values_quad[],
                     VectorizedArray<Number> *gradients_quad[][dim],
                     VectorizedArray<Number> *scratch_data,
-                    const bool               evaluate_val,
-                    const bool               evaluate_grad);
+                    const bool               evaluate_values,
+                    const bool               evaluate_gradients);
   };
 
 
@@ -131,22 +132,22 @@ namespace internal
               VectorizedArray<Number> *gradients_quad[][dim],
               VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
               VectorizedArray<Number> *scratch_data,
-              const bool               evaluate_val,
-              const bool               evaluate_grad,
-              const bool               evaluate_lapl)
+              const bool               evaluate_values,
+              const bool               evaluate_gradients,
+              const bool               evaluate_hessians)
   {
-    if (evaluate_val == false && evaluate_grad == false && evaluate_lapl == false)
+    if (evaluate_values == false && evaluate_gradients == false && evaluate_hessians == false)
       return;
 
     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 :
+    Eval eval (variant == evaluate_evenodd ? shape_info.shape_values_eo :
                shape_info.shape_values,
-               variant == evaluate_evenodd ? shape_info.shape_gra_evenodd :
+               variant == evaluate_evenodd ? shape_info.shape_gradients_eo :
                shape_info.shape_gradients,
-               variant == evaluate_evenodd ? shape_info.shape_hes_evenodd :
+               variant == evaluate_evenodd ? shape_info.shape_hessians_eo :
                shape_info.shape_hessians,
                shape_info.fe_degree,
                shape_info.n_q_points_1d);
@@ -218,11 +219,11 @@ namespace internal
       case 1:
         for (unsigned int c=0; c<n_components; c++)
           {
-            if (evaluate_val == true)
+            if (evaluate_values == true)
               eval.template values<0,true,false> (values_dofs[c], values_quad[c]);
-            if (evaluate_grad == true)
+            if (evaluate_gradients == true)
               eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]);
-            if (evaluate_lapl == true)
+            if (evaluate_hessians == true)
               eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]);
           }
         break;
@@ -231,15 +232,15 @@ namespace internal
         for (unsigned int c=0; c<n_components; c++)
           {
             // grad x
-            if (evaluate_grad == true)
+            if (evaluate_gradients == true)
               {
                 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)
+            if (evaluate_hessians == true)
               {
                 // grad xy
-                if (evaluate_grad == false)
+                if (evaluate_gradients == false)
                   eval.template gradients<0,true,false>(values_dofs[c], temp1);
                 eval.template gradients<1,true,false>  (temp1, hessians_quad[c][d1+d1]);
 
@@ -250,15 +251,15 @@ namespace internal
 
             // grad y
             eval.template values<0,true,false> (values_dofs[c], temp1);
-            if (evaluate_grad == true)
+            if (evaluate_gradients == true)
               eval.template gradients<1,true,false> (temp1, gradients_quad[c][d1]);
 
             // grad yy
-            if (evaluate_lapl == true)
+            if (evaluate_hessians == true)
               eval.template hessians<1,true,false> (temp1, hessians_quad[c][d1]);
 
             // val: can use values applied in x
-            if (evaluate_val == true)
+            if (evaluate_values == true)
               eval.template values<1,true,false> (temp1, values_quad[c]);
           }
         break;
@@ -266,7 +267,7 @@ namespace internal
       case 3:
         for (unsigned int c=0; c<n_components; c++)
           {
-            if (evaluate_grad == true)
+            if (evaluate_gradients == true)
               {
                 // grad x
                 eval.template gradients<0,true,false> (values_dofs[c], temp1);
@@ -274,10 +275,10 @@ namespace internal
                 eval.template values<2,true,false> (temp2, gradients_quad[c][0]);
               }
 
-            if (evaluate_lapl == true)
+            if (evaluate_hessians == true)
               {
                 // grad xz
-                if (evaluate_grad == false)
+                if (evaluate_gradients == false)
                   {
                     eval.template gradients<0,true,false> (values_dofs[c], temp1);
                     eval.template values<1,true,false> (temp1, temp2);
@@ -296,16 +297,16 @@ namespace internal
 
             // grad y
             eval.template values<0,true,false> (values_dofs[c], temp1);
-            if (evaluate_grad == true)
+            if (evaluate_gradients == true)
               {
                 eval.template gradients<1,true,false>(temp1, temp2);
                 eval.template values<2,true,false>   (temp2, gradients_quad[c][d1]);
               }
 
-            if (evaluate_lapl == true)
+            if (evaluate_hessians == true)
               {
                 // grad yz
-                if (evaluate_grad == false)
+                if (evaluate_gradients == false)
                   eval.template gradients<1,true,false>(temp1, temp2);
                 eval.template gradients<2,true,false>  (temp2, hessians_quad[c][d5]);
 
@@ -316,16 +317,16 @@ namespace internal
 
             // 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)
+            if (evaluate_gradients == 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)
+            if (evaluate_hessians == true)
               eval.template hessians<2,true,false>(temp2, hessians_quad[c][d2]);
 
             // val: can use the values applied in x & y direction stored in temp2
-            if (evaluate_val == true)
+            if (evaluate_values == true)
               eval.template values<2,true,false> (temp2, values_quad[c]);
           }
         break;
@@ -336,7 +337,7 @@ namespace internal
 
     // 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)
+    if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 && evaluate_values)
       for (unsigned int c=0; c<n_components; ++c)
         for (unsigned int q=0; q<shape_info.n_q_points; ++q)
           values_quad[c][q] += values_dofs[c][shape_info.dofs_per_cell-1];
@@ -354,18 +355,18 @@ namespace internal
                VectorizedArray<Number> *values_quad[],
                VectorizedArray<Number> *gradients_quad[][dim],
                VectorizedArray<Number> *scratch_data,
-               const bool               integrate_val,
-               const bool               integrate_grad)
+               const bool               integrate_values,
+               const bool               integrate_gradients)
   {
     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 :
+    Eval eval (variant == evaluate_evenodd ? shape_info.shape_values_eo :
                shape_info.shape_values,
-               variant == evaluate_evenodd ? shape_info.shape_gra_evenodd :
+               variant == evaluate_evenodd ? shape_info.shape_gradients_eo :
                shape_info.shape_gradients,
-               variant == evaluate_evenodd ? shape_info.shape_hes_evenodd :
+               variant == evaluate_evenodd ? shape_info.shape_hessians_eo :
                shape_info.shape_hessians,
                shape_info.fe_degree,
                shape_info.n_q_points_1d);
@@ -416,11 +417,11 @@ namespace internal
       case 1:
         for (unsigned int c=0; c<n_components; c++)
           {
-            if (integrate_val == true)
+            if (integrate_values == true)
               eval.template values<0,false,false> (values_quad[c], values_dofs[c]);
-            if (integrate_grad == true)
+            if (integrate_gradients == true)
               {
-                if (integrate_val == true)
+                if (integrate_values == 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]);
@@ -431,20 +432,20 @@ namespace internal
       case 2:
         for (unsigned int c=0; c<n_components; c++)
           {
-            if (integrate_val == true)
+            if (integrate_values == true)
               {
                 // val
                 eval.template values<0,false,false> (values_quad[c], temp1);
                 //grad x
-                if (integrate_grad == true)
+                if (integrate_gradients == 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)
+            if (integrate_gradients == true)
               {
                 // grad y
                 eval.template values<0,false,false>  (gradients_quad[c][d1], temp1);
-                if (integrate_val == false)
+                if (integrate_values == false)
                   {
                     eval.template gradients<1,false,false>(temp1, values_dofs[c]);
                     //grad x
@@ -460,22 +461,22 @@ namespace internal
       case 3:
         for (unsigned int c=0; c<n_components; c++)
           {
-            if (integrate_val == true)
+            if (integrate_values == 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)
+                if (integrate_gradients == true)
                   eval.template gradients<0,false,true> (gradients_quad[c][0], temp1);
                 eval.template values<1,false,false>(temp1, temp2);
-                if (integrate_grad == true)
+                if (integrate_gradients == true)
                   {
                     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]);
               }
-            else if (integrate_grad == true)
+            else if (integrate_gradients == true)
               {
                 eval.template gradients<0,false,false>(gradients_quad[c][0], temp1);
                 eval.template values<1,false,false> (temp1, temp2);
@@ -483,7 +484,7 @@ namespace internal
                 eval.template gradients<1,false,true>(temp1, temp2);
                 eval.template values<2,false,false> (temp2, values_dofs[c]);
               }
-            if (integrate_grad == true)
+            if (integrate_gradients == true)
               {
                 // grad z: can sum to temporary x and y value in output
                 eval.template values<0,false,false> (gradients_quad[c][d2], temp1);
@@ -500,7 +501,7 @@ namespace internal
     // case FE_Q_DG0: add values, gradients and second derivatives are zero
     if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0)
       {
-        if (integrate_val)
+        if (integrate_values)
           for (unsigned int c=0; c<n_components; ++c)
             {
               values_dofs[c][shape_info.dofs_per_cell-1] = values_quad[c][0];
@@ -537,16 +538,18 @@ namespace internal
 
   /**
    * This struct performs the evaluation of function values, gradients and
-   * Hessians for tensor-product finite elements.  This a specialization for
-   * symmetric basis functions with the same number of quadrature points as
-   * degrees of freedom. In that case, we can first transform to a spectral
-   * basis and then perform the evaluation of the first and second derivatives
-   * in spectral space, using the identity operation for the shape values.
+   * Hessians for tensor-product finite elements. This a specialization for
+   * symmetric basis functions about the mid point 0.5 of the unit interval
+   * with the same number of quadrature points as degrees of freedom. In that
+   * case, we can first transform the basis to one that has the nodal points
+   * in the quadrature points (i.e., the collocation space) and then perform
+   * the evaluation of the first and second derivatives in this transformed
+   * space, using the identity operation for the shape values.
    *
    * @author Katharina Kormann, Martin Kronbichler, 2017
    */
   template <int dim, int fe_degree, int n_components, typename Number>
-  struct FEEvaluationImplTransformSpectral
+  struct FEEvaluationImplTransformToCollocation
   {
     static
     void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
@@ -555,9 +558,9 @@ namespace internal
                    VectorizedArray<Number> *gradients_quad[][dim],
                    VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
                    VectorizedArray<Number> *scratch_data,
-                   const bool               evaluate_val,
-                   const bool               evaluate_grad,
-                   const bool               evaluate_lapl);
+                   const bool               evaluate_values,
+                   const bool               evaluate_gradients,
+                   const bool               evaluate_hessians);
 
     static
     void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
@@ -565,14 +568,14 @@ namespace internal
                     VectorizedArray<Number> *values_quad[],
                     VectorizedArray<Number> *gradients_quad[][dim],
                     VectorizedArray<Number> *scratch_data,
-                    const bool               integrate_val,
-                    const bool               integrate_grad);
+                    const bool               integrate_values,
+                    const bool               integrate_gradients);
   };
 
   template <int dim, int fe_degree, int n_components, typename Number>
   inline
   void
-  FEEvaluationImplTransformSpectral<dim, fe_degree, n_components, Number>
+  FEEvaluationImplTransformToCollocation<dim, fe_degree, n_components, Number>
   ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
               VectorizedArray<Number> *values_dofs[],
               VectorizedArray<Number> *values_quad[],
@@ -580,19 +583,19 @@ namespace internal
               VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
               VectorizedArray<Number> *,
               const bool               ,
-              const bool               evaluate_grad,
-              const bool               evaluate_lapl)
+              const bool               evaluate_gradients,
+              const bool               evaluate_hessians)
   {
     typedef EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree, fe_degree+1,
             VectorizedArray<Number> > Eval;
-    Eval eval_val (shape_info.shape_val_evenodd,
-                   shape_info.shape_gra_evenodd,
-                   shape_info.shape_hes_evenodd,
+    Eval eval_val (shape_info.shape_values_eo,
+                   AlignedVector<VectorizedArray<Number> >(),
+                   AlignedVector<VectorizedArray<Number> >(),
                    shape_info.fe_degree,
                    shape_info.n_q_points_1d);
-    Eval eval(shape_info.shape_values,
-              shape_info.shape_grad_spectral_eo,
-              shape_info.shape_hessian_spectral_eo,
+    Eval eval(AlignedVector<VectorizedArray<Number> >(),
+              shape_info.shape_gradients_collocation_eo,
+              shape_info.shape_hessians_collocation_eo,
               shape_info.fe_degree,
               shape_info.n_q_points_1d);
 
@@ -607,7 +610,9 @@ namespace internal
 
     for (unsigned int c=0; c<n_components; c++)
       {
-        // transform to spectral coordinates
+        // transform to the basis functions of the collocation space. use
+        // gradients_quad[c][0] as a temporary array (it gets overwritten by
+        // the gradient contributions later)
         if (dim == 1)
           eval_val.template values<0,true,false>(values_dofs[c], values_quad[c]);
         else if (dim == 2)
@@ -622,31 +627,29 @@ namespace internal
             eval_val.template values<2,true,false>(gradients_quad[c][0], values_quad[c]);
           }
 
-        // apply derivatives in spectral space
-        if (evaluate_grad == true)
+        // apply derivatives in the collocation space
+        if (evaluate_gradients == true || evaluate_hessians == true)
           {
             eval.template gradients<0,true,false>(values_quad[c], gradients_quad[c][0]);
-            if (dim >= 2)
+            if (dim > 1)
               eval.template gradients<1,true,false>(values_quad[c], gradients_quad[c][d1]);
-            if (dim >= 3)
+            if (dim > 2)
               eval.template gradients<2,true,false>(values_quad[c], gradients_quad[c][d2]);
           }
-        if (evaluate_lapl == true)
+        if (evaluate_hessians == true)
           {
             eval.template hessians<0,true,false> (values_quad[c], hessians_quad[c][0]);
             if (dim > 1)
               {
-                eval.template gradients<0,true,false> (values_quad[c], hessians_quad[c][d2]);
-                eval.template gradients<1,true,false> (hessians_quad[c][d2], hessians_quad[c][d3]);
+                // re-use grad_x already in gradients
+                eval.template gradients<1,true,false> (gradients_quad[c][0], hessians_quad[c][d3]);
                 eval.template hessians<1,true,false> (values_quad[c], hessians_quad[c][d1]);
               }
             if (dim > 2)
               {
-                // note that grad x is already in hessians_quad[c][d2]
-                eval.template gradients<2,true,false> (hessians_quad[c][d2], hessians_quad[c][d4]);
-
-                eval.template gradients<1,true,false> (values_quad[c], hessians_quad[c][d2]);
-                eval.template gradients<2,true,false> (hessians_quad[c][d2], hessians_quad[c][d5]);
+                // re-use grad_x and grad_y already in gradients
+                eval.template gradients<2,true,false> (gradients_quad[c][0], hessians_quad[c][d4]);
+                eval.template gradients<2,true,false> (gradients_quad[c][d1], hessians_quad[c][d5]);
                 eval.template hessians<2,true,false> (values_quad[c], hessians_quad[c][d2]);
               }
           }
@@ -656,25 +659,25 @@ namespace internal
   template <int dim, int fe_degree, int n_components, typename Number>
   inline
   void
-  FEEvaluationImplTransformSpectral<dim, fe_degree, n_components, Number>
+  FEEvaluationImplTransformToCollocation<dim, fe_degree, n_components, Number>
   ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
                VectorizedArray<Number> *values_dofs[],
                VectorizedArray<Number> *values_quad[],
                VectorizedArray<Number> *gradients_quad[][dim],
                VectorizedArray<Number> *,
-               const bool               integrate_val,
-               const bool               integrate_grad)
+               const bool               integrate_values,
+               const bool               integrate_gradients)
   {
     typedef EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree, fe_degree+1,
             VectorizedArray<Number> > Eval;
-    Eval eval_val (shape_info.shape_val_evenodd,
-                   shape_info.shape_gra_evenodd,
-                   shape_info.shape_hes_evenodd,
+    Eval eval_val (shape_info.shape_values_eo,
+                   AlignedVector<VectorizedArray<Number> >(),
+                   AlignedVector<VectorizedArray<Number> >(),
                    shape_info.fe_degree,
                    shape_info.n_q_points_1d);
-    Eval eval(shape_info.shape_values,
-              shape_info.shape_grad_spectral_eo,
-              shape_info.shape_hessian_spectral_eo,
+    Eval eval(AlignedVector<VectorizedArray<Number> >(),
+              shape_info.shape_gradients_collocation_eo,
+              shape_info.shape_hessians_collocation_eo,
               shape_info.fe_degree,
               shape_info.n_q_points_1d);
 
@@ -686,20 +689,20 @@ namespace internal
 
     for (unsigned int c=0; c<n_components; c++)
       {
-        // apply derivatives in spectral space
-        if (integrate_grad == true)
+        // apply derivatives in collocation space
+        if (integrate_gradients == true)
           {
-            if (integrate_val)
+            if (integrate_values)
               eval.template gradients<0,false,true>(gradients_quad[c][0], values_quad[c]);
             else
               eval.template gradients<0,false,false>(gradients_quad[c][0], values_quad[c]);
-            if (dim >= 2)
+            if (dim > 1)
               eval.template gradients<1,false,true>(gradients_quad[c][d1], values_quad[c]);
-            if (dim >= 3)
+            if (dim > 2)
               eval.template gradients<2,false,true>(gradients_quad[c][d2], values_quad[c]);
           }
 
-        // transform back to the usual space
+        // transform back to the original space
         if (dim == 1)
           eval_val.template values<0,false,false>(values_quad[c], values_dofs[c]);
         else if (dim == 2)
@@ -720,14 +723,20 @@ namespace internal
 
   /**
    * This struct performs the evaluation of function values, gradients and
-   * Hessians for tensor-product finite elements.  This a specialization for
-   * "spectral" elements like Gauss-Lobatto elements where the 'values'
-   * operation is identity, which allows us to write shorter code.
+   * Hessians for tensor-product finite elements. This a specialization for
+   * elements where the nodal points coincide with the quadrature points like
+   * FE_Q shape functions on Gauss-Lobatto elements integrated with
+   * Gauss-Lobatto quadrature. The assumption of this class is that the shape
+   * 'values' operation is identity, which allows us to write shorter code.
+   *
+   * In literature, this form of evaluation is often called spectral
+   * evaluation, spectral collocation or simply collocation, meaning the same
+   * location for shape functions and evaluation space (quadrature points).
    *
    * @author Katharina Kormann, 2012
   */
   template <int dim, int fe_degree, int n_components, typename Number>
-  struct FEEvaluationImplSpectral
+  struct FEEvaluationImplCollocation
   {
     static
     void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
@@ -736,9 +745,9 @@ namespace internal
                    VectorizedArray<Number> *gradients_quad[][dim],
                    VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
                    VectorizedArray<Number> *scratch_data,
-                   const bool               evaluate_val,
-                   const bool               evaluate_grad,
-                   const bool               evaluate_lapl);
+                   const bool               evaluate_values,
+                   const bool               evaluate_gradients,
+                   const bool               evaluate_hessians);
 
     static
     void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
@@ -746,31 +755,31 @@ namespace internal
                     VectorizedArray<Number> *values_quad[],
                     VectorizedArray<Number> *gradients_quad[][dim],
                     VectorizedArray<Number> *scratch_data,
-                    const bool               integrate_val,
-                    const bool               integrate_grad);
+                    const bool               integrate_values,
+                    const bool               integrate_gradients);
   };
 
   template <int dim, int fe_degree, int n_components, typename Number>
   inline
   void
-  FEEvaluationImplSpectral<dim, fe_degree, n_components, Number>
+  FEEvaluationImplCollocation<dim, fe_degree, 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],
               VectorizedArray<Number> *,
-              const bool               evaluate_val,
-              const bool               evaluate_grad,
-              const bool               evaluate_lapl)
+              const bool               evaluate_values,
+              const bool               evaluate_gradients,
+              const bool               evaluate_hessians)
   {
     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,
-               shape_info.fe_degree,
-               shape_info.n_q_points_1d);
+    Eval eval(AlignedVector<VectorizedArray<Number> >(),
+              shape_info.shape_gradients_eo,
+              shape_info.shape_hessians_eo,
+              shape_info.fe_degree,
+              shape_info.n_q_points_1d);
 
     // These avoid compiler warnings; they are only used in sensible context
     // but compilers typically cannot detect when we access something like
@@ -783,33 +792,31 @@ namespace internal
 
     for (unsigned int c=0; c<n_components; c++)
       {
-        if (evaluate_val == true)
+        if (evaluate_values == true)
           for (unsigned int i=0; i<Eval::dofs_per_cell; ++i)
             values_quad[c][i] = values_dofs[c][i];
-        if (evaluate_grad == true)
+        if (evaluate_gradients == true || evaluate_hessians == true)
           {
             eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]);
-            if (dim >= 2)
+            if (dim > 1)
               eval.template gradients<1,true,false>(values_dofs[c], gradients_quad[c][d1]);
-            if (dim >= 3)
+            if (dim > 2)
               eval.template gradients<2,true,false>(values_dofs[c], gradients_quad[c][d2]);
           }
-        if (evaluate_lapl == true)
+        if (evaluate_hessians == true)
           {
             eval.template hessians<0,true,false> (values_quad[c], hessians_quad[c][0]);
             if (dim > 1)
               {
-                eval.template gradients<0,true,false> (values_dofs[c], hessians_quad[c][d2]);
-                eval.template gradients<1,true,false> (hessians_quad[c][d2], hessians_quad[c][d3]);
+                // re-use grad_x already in gradients
+                eval.template gradients<1,true,false> (gradients_quad[c][0], hessians_quad[c][d3]);
                 eval.template hessians<1,true,false> (values_dofs[c], hessians_quad[c][d1]);
               }
             if (dim > 2)
               {
-                // note that grad x is already in hessians_quad[c][d2]
-                eval.template gradients<2,true,false> (hessians_quad[c][d2], hessians_quad[c][d4]);
-
-                eval.template gradients<1,true,false> (values_dofs[c], hessians_quad[c][d2]);
-                eval.template gradients<2,true,false> (hessians_quad[c][d2], hessians_quad[c][d5]);
+                // re-use grad_x already in gradients
+                eval.template gradients<2,true,false> (gradients_quad[c][0], hessians_quad[c][d4]);
+                eval.template gradients<2,true,false> (gradients_quad[c][d1], hessians_quad[c][d5]);
                 eval.template hessians<2,true,false> (values_dofs[c], hessians_quad[c][d2]);
               }
           }
@@ -819,22 +826,22 @@ namespace internal
   template <int dim, int fe_degree, int n_components, typename Number>
   inline
   void
-  FEEvaluationImplSpectral<dim, fe_degree, n_components, Number>
+  FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>
   ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
                VectorizedArray<Number> *values_dofs[],
                VectorizedArray<Number> *values_quad[],
                VectorizedArray<Number> *gradients_quad[][dim],
                VectorizedArray<Number> *,
-               const bool               integrate_val,
-               const bool               integrate_grad)
+               const bool               integrate_values,
+               const bool               integrate_gradients)
   {
     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,
-               shape_info.fe_degree,
-               shape_info.n_q_points_1d);
+    Eval eval(AlignedVector<VectorizedArray<Number> >(),
+              shape_info.shape_gradients_eo,
+              shape_info.shape_hessians_eo,
+              shape_info.fe_degree,
+              shape_info.n_q_points_1d);
 
     // These avoid compiler warnings; they are only used in sensible context
     // but compilers typically cannot detect when we access something like
@@ -844,18 +851,18 @@ namespace internal
 
     for (unsigned int c=0; c<n_components; c++)
       {
-        if (integrate_val == true)
+        if (integrate_values == true)
           for (unsigned int i=0; i<Eval::dofs_per_cell; ++i)
             values_dofs[c][i] = values_quad[c][i];
-        if (integrate_grad == true)
+        if (integrate_gradients == true)
           {
-            if (integrate_val == true)
+            if (integrate_values == 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]);
-            if (dim >= 2)
+            if (dim > 1)
               eval.template gradients<1,false,true>(gradients_quad[c][d1], values_dofs[c]);
-            if (dim >= 3)
+            if (dim > 2)
               eval.template gradients<2,false,true>(gradients_quad[c][d2], values_dofs[c]);
           }
       }

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 906fe539116b597441a6428caf245aca2cc656d4..3eb419d96716ff817ed291df096fd23058b3a242 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -2741,6 +2741,8 @@ namespace internal
       res = vector_access (const_cast<const VectorType &>(vec), index);
     }
 
+    // variant where VectorType::value_type is the same as Number -> can call
+    // gather
     template <typename VectorType>
     void process_dof_gather (const unsigned int      *indices,
                              VectorType              &vec,
@@ -2750,6 +2752,8 @@ namespace internal
       res.gather(vec.begin(), indices);
     }
 
+    // variant where VectorType::value_type is not the same as Number -> must
+    // manually load the data
     template <typename VectorType>
     void process_dof_gather (const unsigned int      *indices,
                              VectorType              &vec,
@@ -2807,6 +2811,8 @@ namespace internal
       vector_access (vec, index) += res;
     }
 
+    // variant where VectorType::value_type is the same as Number -> can call
+    // scatter
     template <typename VectorType>
     void process_dof_gather (const unsigned int      *indices,
                              VectorType              &vec,
@@ -2814,7 +2820,6 @@ namespace internal
                              internal::bool2type<true>) const
     {
       // TODO: enable scatter path when indices are fixed
-
       //#if DEAL_II_COMPILER_VECTORIZATION_LEVEL < 3
 #if 1
       for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
@@ -2828,6 +2833,8 @@ namespace internal
 #endif
     }
 
+    // variant where VectorType::value_type is not the same as Number -> must
+    // manually append all data
     template <typename VectorType>
     void process_dof_gather (const unsigned int      *indices,
                              VectorType              &vec,
@@ -5408,9 +5415,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
   else
     {
       if (fe_degree+1 == n_q_points_1d &&
-          this->data->element_type == internal::MatrixFreeFunctions::tensor_gausslobatto)
+          this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation)
         {
-          internal::FEEvaluationImplSpectral<dim, fe_degree_templ, n_components_, Number>
+          internal::FEEvaluationImplCollocation<dim, fe_degree_templ, n_components_, Number>
           ::evaluate(*this->data, &this->values_dofs[0], this->values_quad,
                      this->gradients_quad, this->hessians_quad, this->scratch_data,
                      evaluate_val, evaluate_grad, evaluate_lapl);
@@ -5418,7 +5425,7 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
       else if (fe_degree+1 == n_q_points_1d &&
                this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric)
         {
-          internal::FEEvaluationImplTransformSpectral<dim, fe_degree_templ, n_components_, Number>
+          internal::FEEvaluationImplTransformToCollocation<dim, fe_degree_templ, n_components_, Number>
           ::evaluate(*this->data, &this->values_dofs[0], this->values_quad,
                      this->gradients_quad, this->hessians_quad, this->scratch_data,
                      evaluate_val, evaluate_grad, evaluate_lapl);
@@ -5520,9 +5527,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
   else
     {
       if (fe_degree+1 == n_q_points_1d &&
-          this->data->element_type == internal::MatrixFreeFunctions::tensor_gausslobatto)
+          this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation)
         {
-          internal::FEEvaluationImplSpectral<dim, fe_degree_templ, n_components_, Number>
+          internal::FEEvaluationImplCollocation<dim, fe_degree_templ, n_components_, Number>
           ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
                       this->gradients_quad, this->scratch_data,
                       integrate_val, integrate_grad);
@@ -5530,7 +5537,7 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
       else if (fe_degree+1 == n_q_points_1d &&
                this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric)
         {
-          internal::FEEvaluationImplTransformSpectral<dim, fe_degree_templ, n_components_, Number>
+          internal::FEEvaluationImplTransformToCollocation<dim, fe_degree_templ, n_components_, Number>
           ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
                       this->gradients_quad, this->scratch_data,
                       integrate_val, integrate_grad);

diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h
index 95b981f43baddace09399f1964516f140b9c4dfe..eafeda2f72f4c88e207c8b631a365cf1343f2290 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.h
+++ b/include/deal.II/matrix_free/shape_info.h
@@ -40,11 +40,42 @@ namespace internal
      */
     enum ElementType
     {
-      tensor_general,
-      tensor_symmetric,
-      truncated_tensor,
-      tensor_symmetric_plus_dg0,
-      tensor_gausslobatto
+      /**
+       * Tensor product shape function where the shape value evaluation in the
+       * quadrature point corresponds to the identity operation and no
+       * interpolation needs to be performed (collocation approach, also
+       * called spectral evaluation). This is for example the case for an
+       * element with nodes in the Gauss-Lobatto support points and
+       * integration in the Gauss-Lobatto quadrature points of the same order.
+       */
+      tensor_symmetric_collocation = 0,
+      /**
+       * Symmetric tensor product shape functions fulfilling a Hermite
+       * identity with values and first derivatives zero at the element end
+       * points in 1D.
+       */
+      tensor_symmetric_hermite = 1,
+      /**
+       * Usual tensor product shape functions whose shape values and
+       * quadrature points are symmetric about the midpoint of the unit
+       * interval 0.5
+       */
+      tensor_symmetric = 2,
+      /**
+       * Tensor product shape functions without further particular properties
+       */
+      tensor_general = 3,
+      /**
+       * Polynomials of complete degree rather than tensor degree which can be
+       * described by a truncated tensor product
+       */
+      truncated_tensor = 4,
+      /**
+       * Tensor product shape functions that are symmetric about the midpoint
+       * of the unit interval 0.5 that additionally add a constant shape
+       * function according to FE_Q_DG0.
+       */
+      tensor_symmetric_plus_dg0 = 5
     };
 
     /**
@@ -128,33 +159,35 @@ namespace internal
        * even-odd scheme where the symmetries in shape_values are used for
        * faster evaluation.
        */
-      AlignedVector<VectorizedArray<Number> > shape_val_evenodd;
+      AlignedVector<VectorizedArray<Number> > shape_values_eo;
 
       /**
        * 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;
+      AlignedVector<VectorizedArray<Number> > shape_gradients_eo;
 
       /**
        * 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;
+      AlignedVector<VectorizedArray<Number> > shape_hessians_eo;
 
       /**
-       * Stores the shape gradients of the spectral element space
-       * FE_DGQ<1>(Quadrature<1>) for faster evaluation in even-odd format.
+       * Stores the shape gradients of the shape function space associated to
+       * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). For
+       * faster evaluation only the even-odd format is necessary.
        */
-      AlignedVector<VectorizedArray<Number> > shape_grad_spectral_eo;
+      AlignedVector<VectorizedArray<Number> > shape_gradients_collocation_eo;
 
       /**
-       * Stores the shape hessians of the spectral element space
-       * FE_DGQ<1>(Quadrature<1>) for faster evaluation in even-odd format.
+       * Stores the shape hessians of the shape function space associated to
+       * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). For
+       * faster evaluation only the even-odd format is necessary.
        */
-      AlignedVector<VectorizedArray<Number> > shape_hessian_spectral_eo;
+      AlignedVector<VectorizedArray<Number> > shape_hessians_collocation_eo;
 
       /**
        * Stores the indices from cell DoFs to face DoFs. The rows go through
@@ -234,16 +267,17 @@ namespace internal
 
       /**
        * Check whether we have symmetries in the shape values. In that case,
-       * also fill the shape_???_evenodd fields.
+       * also fill the shape_???_eo fields.
        */
       bool check_1d_shapes_symmetric(const unsigned int n_q_points_1d);
 
       /**
        * Check 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.
+       * values form a diagonal matrix, i.e., the nodal points are collocated
+       * with the quadrature points. This allows for specialized algorithms
+       * that save some operations in the evaluation.
        */
-      bool check_1d_shapes_gausslobatto();
+      bool check_1d_shapes_collocation();
     };
 
 
@@ -257,6 +291,7 @@ namespace internal
                                   const FiniteElement<dim> &fe_in,
                                   const unsigned int base_element_number)
       :
+      element_type(tensor_general),
       fe_degree (0),
       n_q_points_1d (0),
       n_q_points (0),
@@ -267,9 +302,8 @@ namespace internal
       reinit (quad, fe_in, base_element_number);
     }
 
-
-
   } // end of namespace MatrixFreeFunctions
+
 } // end of namespace internal
 
 DEAL_II_NAMESPACE_CLOSE

diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index 9d7b32818b53a80320622c4e0146888c2e3aa579..45918788d90c9aa92b48c9808e3c35095a475e1b 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -212,40 +212,42 @@ namespace internal
       if (element_type == tensor_general &&
           check_1d_shapes_symmetric(n_q_points_1d))
         {
-          if (check_1d_shapes_gausslobatto())
-            element_type = tensor_gausslobatto;
+          if (check_1d_shapes_collocation())
+            element_type = tensor_symmetric_collocation;
           else
             {
               element_type = tensor_symmetric;
-              // get spectral evaluation points
+              // get gradient and Hessian transformation matrix for the
+              // polynomial space associated with the quadrature rule
+              // (collocation space)
               if (fe_degree+1 == n_q_points_1d)
                 {
                   const unsigned int stride = fe_degree/2+1;
-                  shape_grad_spectral_eo.resize((fe_degree+1)*stride);
-                  shape_hessian_spectral_eo.resize((fe_degree+1)*stride);
+                  shape_gradients_collocation_eo.resize((fe_degree+1)*stride);
+                  shape_hessians_collocation_eo.resize((fe_degree+1)*stride);
                   FE_DGQArbitraryNodes<1> fe(quad.get_points());
                   for (unsigned int i=0; i<(fe_degree+1)/2; ++i)
                     for (unsigned int q=0; q<stride; ++q)
                       {
-                        shape_grad_spectral_eo[i*stride+q] =
+                        shape_gradients_collocation_eo[i*stride+q] =
                           0.5* (fe.shape_grad(i, quad.get_points()[q])[0] +
                                 fe.shape_grad(i, quad.get_points()[n_q_points_1d-1-q])[0]);
-                        shape_grad_spectral_eo[(fe_degree-i)*stride+q] =
+                        shape_gradients_collocation_eo[(fe_degree-i)*stride+q] =
                           0.5* (fe.shape_grad(i, quad.get_points()[q])[0] -
                                 fe.shape_grad(i, quad.get_points()[n_q_points_1d-1-q])[0]);
-                        shape_hessian_spectral_eo[i*stride+q] =
+                        shape_hessians_collocation_eo[i*stride+q] =
                           0.5* (fe.shape_grad_grad(i, quad.get_points()[q])[0][0] +
                                 fe.shape_grad_grad(i, quad.get_points()[n_q_points_1d-1-q])[0][0]);
-                        shape_hessian_spectral_eo[(fe_degree-i)*stride+q] =
+                        shape_hessians_collocation_eo[(fe_degree-i)*stride+q] =
                           0.5* (fe.shape_grad_grad(i, quad.get_points()[q])[0][0] -
                                 fe.shape_grad_grad(i, quad.get_points()[n_q_points_1d-1-q])[0][0]);
                       }
                   if (fe_degree % 2 == 0)
                     for (unsigned int q=0; q<stride; ++q)
                       {
-                        shape_grad_spectral_eo[fe_degree/2*stride+q] =
+                        shape_gradients_collocation_eo[fe_degree/2*stride+q] =
                           fe.shape_grad(fe_degree/2, quad.get_points()[q])[0];
-                        shape_hessian_spectral_eo[fe_degree/2*stride+q] =
+                        shape_hessians_collocation_eo[fe_degree/2*stride+q] =
                           fe.shape_grad_grad(fe_degree/2, quad.get_points()[q])[0][0];
                       }
                 }
@@ -352,41 +354,41 @@ namespace internal
             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);
+      shape_values_eo.resize((fe_degree+1)*stride);
+      shape_gradients_eo.resize((fe_degree+1)*stride);
+      shape_hessians_eo.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] =
+            shape_values_eo[i*stride+q] =
               Number(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] =
+            shape_values_eo[(fe_degree-i)*stride+q] =
               Number(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] =
+            shape_gradients_eo[i*stride+q] =
               Number(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] =
+            shape_gradients_eo[(fe_degree-i)*stride+q] =
               Number(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] =
+            shape_hessians_eo[i*stride+q] =
               Number(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] =
+            shape_hessians_eo[(fe_degree-i)*stride+q] =
               Number(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_eo[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_eo[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_eo[fe_degree/2*stride+q] =
               shape_hessians[(fe_degree/2)*n_q_points_1d+q];
           }
 
@@ -397,7 +399,7 @@ namespace internal
 
     template <typename Number>
     bool
-    ShapeInfo<Number>::check_1d_shapes_gausslobatto()
+    ShapeInfo<Number>::check_1d_shapes_collocation()
     {
       if (dofs_per_cell != n_q_points)
         return false;
@@ -441,11 +443,11 @@ namespace internal
       memory += MemoryConsumption::memory_consumption(shape_values);
       memory += MemoryConsumption::memory_consumption(shape_gradients);
       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 += MemoryConsumption::memory_consumption(shape_grad_spectral_eo);
-      memory += MemoryConsumption::memory_consumption(shape_hessian_spectral_eo);
+      memory += MemoryConsumption::memory_consumption(shape_values_eo);
+      memory += MemoryConsumption::memory_consumption(shape_gradients_eo);
+      memory += MemoryConsumption::memory_consumption(shape_hessians_eo);
+      memory += MemoryConsumption::memory_consumption(shape_gradients_collocation_eo);
+      memory += MemoryConsumption::memory_consumption(shape_hessians_collocation_eo);
       memory += face_indices.memory_consumption();
       for (unsigned int i=0; i<2; ++i)
         {