From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Tue, 14 Mar 2017 09:24:42 +0000 (+0100)
Subject: Use simpler code path for quadratic operator evaluation path.
X-Git-Tag: v9.0.0-rc1~1768^2~6
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=56090ff690503220c1cc617eca1a210725681681;p=dealii.git

Use simpler code path for quadratic operator evaluation path.
---

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index 15c0ac4b08..893bd2ac44 100644
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -198,7 +198,7 @@ namespace internal
         AssertDimension(count_q, Utilities::fixed_power<dim>(shape_info.fe_degree+1));
       }
 
-    // These avoid compiler errors; they are only used in sensible context but
+    // These avoid compiler warnings; 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;
@@ -207,6 +207,46 @@ namespace internal
     const unsigned int d4 = dim>2?4:0;
     const unsigned int d5 = dim>2?5:0;
 
+    // check if we can go through the spectral evaluation option which is
+    // faster than the standard one
+    if (fe_degree+1 == n_q_points_1d &&
+        (type == MatrixFreeFunctions::tensor_symmetric ||
+         type == MatrixFreeFunctions::tensor_general) &&
+        evaluate_lapl == false)
+      {
+        Eval eval_grad(shape_info.shape_values,
+                       variant == evaluate_evenodd ? shape_info.shape_grad_spectral_eo :
+                       shape_info.shape_grad_spectral,
+                       shape_info.shape_hessians,
+                       shape_info.fe_degree,
+                       shape_info.n_q_points_1d);
+        for (unsigned int c=0; c<n_components; c++)
+          {
+            if (dim == 1)
+              eval.template values<0,true,false>(values_dofs[c], values_quad[c]);
+            else if (dim == 2)
+              {
+                eval.template values<0,true,false>(values_dofs[c], gradients_quad[c][0]);
+                eval.template values<1,true,false>(gradients_quad[c][0], values_quad[c]);
+              }
+            else if (dim == 3)
+              {
+                eval.template values<0,true,false>(values_dofs[c], values_quad[c]);
+                eval.template values<1,true,false>(values_quad[c], gradients_quad[c][0]);
+                eval.template values<2,true,false>(gradients_quad[c][0], values_quad[c]);
+              }
+            if (evaluate_grad == true)
+              {
+                eval_grad.template gradients<0,true,false>(values_quad[c], gradients_quad[c][0]);
+                if (dim >= 2)
+                  eval_grad.template gradients<1,true,false>(values_quad[c], gradients_quad[c][d1]);
+                if (dim >= 3)
+                  eval_grad.template gradients<2,true,false>(values_quad[c], gradients_quad[c][d2]);
+              }
+          }
+        return;
+      }
+
     switch (dim)
       {
       case 1:
@@ -399,12 +439,54 @@ namespace internal
                                    c*Utilities::fixed_power<dim>(shape_info.fe_degree+1);
       }
 
-    // These avoid compiler errors; they are only used in sensible context but
+    // These avoid compiler warnings; 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;
 
+    // check if we can go through the spectral evaluation option which is
+    // faster than the standard one
+    if (fe_degree+1 == n_q_points_1d &&
+        (type == MatrixFreeFunctions::tensor_symmetric ||
+         type == MatrixFreeFunctions::tensor_general))
+      {
+        Eval eval_grad(shape_info.shape_values,
+                       variant == evaluate_evenodd ? shape_info.shape_grad_spectral_eo :
+                       shape_info.shape_grad_spectral,
+                       shape_info.shape_hessians,
+                       shape_info.fe_degree,
+                       shape_info.n_q_points_1d);
+        for (unsigned int c=0; c<n_components; c++)
+          {
+            if (integrate_grad == true)
+              {
+                if (integrate_val)
+                  eval_grad.template gradients<0,false,true>(gradients_quad[c][0], values_quad[c]);
+                else
+                  eval_grad.template gradients<0,false,false>(gradients_quad[c][0], values_quad[c]);
+                if (dim >= 2)
+                  eval_grad.template gradients<1,false,true>(gradients_quad[c][d1], values_quad[c]);
+                if (dim >= 3)
+                  eval_grad.template gradients<2,false,true>(gradients_quad[c][d2], values_quad[c]);
+              }
+            if (dim == 1)
+              eval.template values<0,false,false>(values_quad[c], values_dofs[c]);
+            else if (dim == 2)
+              {
+                eval.template values<0,false,false>(values_quad[c], gradients_quad[c][0]);
+                eval.template values<1,false,false>(gradients_quad[c][0], values_dofs[c]);
+              }
+            else if (dim == 3)
+              {
+                eval.template values<0,false,false>(values_quad[c], gradients_quad[c][0]);
+                eval.template values<1,false,false>(gradients_quad[c][0], values_quad[c]);
+                eval.template values<2,false,false>(values_quad[c], values_dofs[c]);
+              }
+          }
+        return;
+      }
+
     switch (dim)
       {
       case 1:
@@ -527,8 +609,9 @@ namespace internal
       }
   }
 
-  // This a specialization for Gauss-Lobatto elements where the 'values'
-  // operation is identity, which allows us to write shorter code.
+  // This a specialization for "spectral" elements like 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>
diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h
index cc4962f206..2cda877895 100644
--- a/include/deal.II/matrix_free/shape_info.h
+++ b/include/deal.II/matrix_free/shape_info.h
@@ -144,6 +144,18 @@ namespace internal
        */
       AlignedVector<VectorizedArray<Number> > shape_hes_evenodd;
 
+      /**
+       * Stores the shape gradients of the spectral element space
+       * FE_DGQ<1>(Quadrature<1>) for faster evaluation in even-odd format.
+       */
+      AlignedVector<VectorizedArray<Number> > shape_grad_spectral_eo;
+
+      /**
+       * Stores the shape gradients of the spectral element space
+       * FE_DGQ<1>(Quadrature<1>) for faster evaluation in standard format.
+       */
+      AlignedVector<VectorizedArray<Number> > shape_grad_spectral;
+
       /**
        * 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.
diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index 88b302494a..eb6a2e96b6 100644
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -24,6 +24,7 @@
 #include <deal.II/base/polynomials_piecewise.h>
 #include <deal.II/fe/fe_poly.h>
 #include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_dgq.h>
 #include <deal.II/fe/fe_q_dg0.h>
 
 #include <deal.II/matrix_free/shape_info.h>
@@ -208,6 +209,16 @@ namespace internal
           this->face_gradient[1][i] = fe->shape_grad(my_i,q_point)[0];
         }
 
+      // get spectral evaluation points
+      if (fe_degree+1 == n_q_points_1d)
+        {
+          FE_DGQArbitraryNodes<1> fe(quad.get_points());
+          shape_grad_spectral.resize(n_dofs_1d*n_dofs_1d);
+          for (unsigned int i=0; i<n_dofs_1d; ++i)
+            for (unsigned int q=0; q<n_dofs_1d; ++q)
+              shape_grad_spectral[i*n_q_points_1d+q] = fe.shape_grad(i, quad.get_points()[q])[0];
+        }
+
       if (element_type == tensor_general &&
           check_1d_shapes_symmetric(n_q_points_1d))
         {
@@ -320,6 +331,7 @@ namespace internal
       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_grad_spectral_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)
           {
@@ -343,6 +355,16 @@ namespace internal
             shape_hes_evenodd[(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+1 == n_q_points_1d)
+              {
+                shape_grad_spectral_eo[i*stride+q] =
+                  Number(0.5) * (shape_grad_spectral[i*n_q_points_1d+q] +
+                                 shape_grad_spectral[i*n_q_points_1d+n_q_points_1d-1-q]);
+                shape_grad_spectral_eo[(fe_degree-i)*stride+q] =
+                  Number(0.5) * (shape_grad_spectral[i*n_q_points_1d+q] -
+                                 shape_grad_spectral[i*n_q_points_1d+n_q_points_1d-1-q]);
+              }
           }
       if (fe_degree % 2 == 0)
         for (unsigned int q=0; q<stride; ++q)
@@ -353,6 +375,9 @@ namespace internal
               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];
+            if (fe_degree+1 == n_q_points_1d)
+              shape_grad_spectral_eo[fe_degree/2*stride+q] =
+                shape_grad_spectral[(fe_degree/2)*n_q_points_1d+q];
           }
 
       return true;
@@ -409,6 +434,8 @@ namespace internal
       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);
+      memory += MemoryConsumption::memory_consumption(shape_grad_spectral_eo);
       memory += face_indices.memory_consumption();
       for (unsigned int i=0; i<2; ++i)
         {