]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Put new code into separate function.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 28 Mar 2017 13:24:45 +0000 (15:24 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 28 Mar 2017 13:38:34 +0000 (15:38 +0200)
Do not use function pointers.

include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/shape_info.h
include/deal.II/matrix_free/shape_info.templates.h

index 893bd2ac443c9994f7586ac5f9217c46c33d7ee7..337b5a9618d784821bd8b33513cff91e0956e10a 100644 (file)
@@ -77,17 +77,23 @@ namespace internal
 
 
 
-  // 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.
+  /**
+   * 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 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.
+   *
+   * @author Katharina Kormann, Martin Kronbichler, 2012, 2014, 2017
+   */
   template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
             int n_q_points_1d, int n_components, typename Number>
   struct FEEvaluationImpl
@@ -207,46 +213,6 @@ 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:
@@ -445,48 +411,6 @@ namespace internal
     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:
@@ -609,12 +533,20 @@ namespace internal
       }
   }
 
-  // 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>
+
+
+  /**
+   * 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.
+   *
+   * @author Katharina Kormann, Martin Kronbichler, 2017
+   */
+  template <int dim, int fe_degree, int n_components, typename Number>
+  struct FEEvaluationImplTransformSpectral
   {
     static
     void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
@@ -637,148 +569,200 @@ namespace internal
                     const bool               integrate_grad);
   };
 
-  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
+  template <int dim, int fe_degree, int n_components, typename Number>
   inline
   void
-  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],
-                               VectorizedArray<Number> *scratch_data,
-                               const bool               evaluate_val,
-                               const bool               evaluate_grad,
-                               const bool               evaluate_lapl)
+  FEEvaluationImplTransformSpectral<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               ,
+              const bool               evaluate_grad,
+              const bool               evaluate_lapl)
   {
     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_val (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(shape_info.shape_values,
+              shape_info.shape_grad_spectral_eo,
+              shape_info.shape_hessian_spectral_eo,
+              shape_info.fe_degree,
+              shape_info.n_q_points_1d);
 
-    // 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;
-    const unsigned int d3 = dim>2?3:0;
-    const unsigned int d4 = dim>2?4:0;
-    const unsigned int d5 = dim>2?5:0;
+    const unsigned int d2 = dim>2?2:d1;
+    const unsigned int d3 = d1+d2;
+    const unsigned int d4 = dim>2?4:d3;
+    const unsigned int d5 = dim>2?5:d4;
 
-    switch (dim)
+    for (unsigned int c=0; c<n_components; c++)
       {
-      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++)
+        // transform to spectral coordinates
+        if (dim == 1)
+          eval_val.template values<0,true,false>(values_dofs[c], values_quad[c]);
+        else if (dim == 2)
           {
-            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]);
+            eval_val.template values<0,true,false>(values_dofs[c], gradients_quad[c][0]);
+            eval_val.template values<1,true,false>(gradients_quad[c][0], values_quad[c]);
           }
-        break;
-
-      case 2:
-        if (evaluate_val == true)
+        else if (dim == 3)
           {
-            std::memcpy (values_quad[0], values_dofs[0],
-                         Eval::dofs_per_cell * n_components *
-                         sizeof (values_dofs[0][0]));
+            eval_val.template values<0,true,false>(values_dofs[c], values_quad[c]);
+            eval_val.template values<1,true,false>(values_quad[c], gradients_quad[c][0]);
+            eval_val.template values<2,true,false>(gradients_quad[c][0], values_quad[c]);
           }
+
+        // apply derivatives in spectral space
         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]);
-            }
+          {
+            eval.template gradients<0,true,false>(values_quad[c], gradients_quad[c][0]);
+            if (dim >= 2)
+              eval.template gradients<1,true,false>(values_quad[c], gradients_quad[c][d1]);
+            if (dim >= 3)
+              eval.template gradients<2,true,false>(values_quad[c], gradients_quad[c][d2]);
+          }
         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]);
-
-              // grad x grad y
-              eval.template gradients<0,true,false> (values_dofs[comp], scratch_data);
-              eval.template gradients<1,true,false> (scratch_data, hessians_quad[comp][d1+d1]);
-            }
-        break;
+          {
+            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]);
+                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]);
 
-      case 3:
-        if (evaluate_val == true)
+                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]);
+                eval.template hessians<2,true,false> (values_quad[c], hessians_quad[c][d2]);
+              }
+          }
+      }
+  }
+
+  template <int dim, int fe_degree, int n_components, typename Number>
+  inline
+  void
+  FEEvaluationImplTransformSpectral<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)
+  {
+    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,
+                   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,
+              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
+    // 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++)
+      {
+        // apply derivatives in spectral space
+        if (integrate_grad == true)
           {
-            std::memcpy (values_quad[0], values_dofs[0],
-                         Eval::dofs_per_cell * n_components *
-                         sizeof (values_dofs[0][0]));
+            if (integrate_val)
+              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)
+              eval.template gradients<1,false,true>(gradients_quad[c][d1], values_quad[c]);
+            if (dim >= 3)
+              eval.template gradients<2,false,true>(gradients_quad[c][d2], values_quad[c]);
+          }
+
+        // transform back to the usual space
+        if (dim == 1)
+          eval_val.template values<0,false,false>(values_quad[c], values_dofs[c]);
+        else if (dim == 2)
+          {
+            eval_val.template values<0,false,false>(values_quad[c], gradients_quad[c][0]);
+            eval_val.template values<1,false,false>(gradients_quad[c][0], values_dofs[c]);
+          }
+        else if (dim == 3)
+          {
+            eval_val.template values<0,false,false>(values_quad[c], gradients_quad[c][0]);
+            eval_val.template values<1,false,false>(gradients_quad[c][0], values_quad[c]);
+            eval_val.template values<2,false,false>(values_quad[c], values_dofs[c]);
           }
-        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 = scratch_data;
-              // 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 <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
+
+
+  /**
+   * 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.
+   *
+   * @author Katharina Kormann, 2012
+  */
+  template <int dim, int fe_degree, int n_components, typename Number>
+  struct FEEvaluationImplSpectral
+  {
+    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],
+                   VectorizedArray<Number> *scratch_data,
+                   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],
+                    VectorizedArray<Number> *scratch_data,
+                    const bool               integrate_val,
+                    const bool               integrate_grad);
+  };
+
+  template <int dim, int fe_degree, int n_components, typename Number>
   inline
   void
-  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],
-                                VectorizedArray<Number> *,
-                                const bool               integrate_val,
-                                const bool               integrate_grad)
+  FEEvaluationImplSpectral<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)
   {
     typedef EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree, fe_degree+1,
             VectorizedArray<Number> > Eval;
@@ -788,77 +772,92 @@ namespace internal
                shape_info.fe_degree,
                shape_info.n_q_points_1d);
 
-    // These avoid compiler errors; they are only used in sensible context but
-    // compilers typically cannot detect when we access something like
+    // 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;
+    const unsigned int d2 = dim>2?2:d1;
+    const unsigned int d3 = d1+d2;
+    const unsigned int d4 = dim>2?4:d3;
+    const unsigned int d5 = dim>2?5:d4;
 
-    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)
+    for (unsigned int c=0; c<n_components; c++)
       {
-      case 1:
-        for (unsigned int c=0; c<n_components; c++)
+        if (evaluate_val == 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 (integrate_grad == true)
+            eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]);
+            if (dim >= 2)
+              eval.template gradients<1,true,false>(values_dofs[c], gradients_quad[c][d1]);
+            if (dim >= 3)
+              eval.template gradients<2,true,false>(values_dofs[c], gradients_quad[c][d2]);
+          }
+        if (evaluate_lapl == true)
+          {
+            eval.template hessians<0,true,false> (values_quad[c], hessians_quad[c][0]);
+            if (dim > 1)
               {
-                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]);
+                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]);
+                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]);
+                eval.template hessians<2,true,false> (values_dofs[c], hessians_quad[c][d2]);
               }
           }
+      }
+  }
 
-        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]);
-
-              // grad y
-              eval.template gradients<1, false, true> (gradients_quad[comp][d1],
-                                                       values_dofs[comp]);
-            }
-        break;
+  template <int dim, int fe_degree, int n_components, typename Number>
+  inline
+  void
+  FEEvaluationImplSpectral<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)
+  {
+    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);
 
-      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]);
-
-              // 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;
+    // 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;
 
-      default:
-        AssertThrow(false, ExcNotImplemented());
+    for (unsigned int c=0; c<n_components; c++)
+      {
+        if (integrate_val == 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_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]);
+            if (dim >= 2)
+              eval.template gradients<1,false,true>(gradients_quad[c][d1], values_dofs[c]);
+            if (dim >= 3)
+              eval.template gradients<2,false,true>(gradients_quad[c][d2], values_dofs[c]);
+          }
       }
   }
 
index a8e8bef3e863d490a5f9872d80390e6d127c46bb..143002fc7c200ad9967a27d3dab82049ed23a443 100644 (file)
@@ -2086,30 +2086,6 @@ private:
    */
   void check_template_arguments(const unsigned int fe_no,
                                 const unsigned int first_selected_component);
-
-  /**
-   * 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],
-                          VectorizedArray<Number> *scratch_data,
-                          const bool               evaluate_val,
-                          const bool               evaluate_grad,
-                          const bool               evaluate_lapl);
-
-  /**
-   * Function pointer for the integrate function
-   */
-  void (*integrate_funct)(const internal::MatrixFreeFunctions::ShapeInfo<Number> &,
-                          VectorizedArray<Number> *values_dofs_actual[],
-                          VectorizedArray<Number> *values_quad[],
-                          VectorizedArray<Number> *gradients_quad[][dim],
-                          VectorizedArray<Number> *scratch_data,
-                          const bool               evaluate_val,
-                          const bool               evaluate_grad);
 };
 
 
@@ -5214,94 +5190,6 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 ::check_template_arguments(const unsigned int fe_no,
                            const unsigned int first_selected_component)
 {
-  const unsigned int fe_degree_templ = fe_degree != -1 ? fe_degree : 0;
-  const unsigned int n_q_points_1d_templ = n_q_points_1d > 0 ? n_q_points_1d : 1;
-  if (fe_degree == -1)
-    {
-      if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
-        {
-          evaluate_funct = internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
-          dim, -1, 0, n_components_, Number>::evaluate;
-          integrate_funct = internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
-          dim, -1, 0, n_components_, Number>::integrate;
-        }
-      else if (this->data->element_type == internal::MatrixFreeFunctions::truncated_tensor)
-        {
-          evaluate_funct = internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
-          dim, -1, 0, n_components_, Number>::evaluate;
-          integrate_funct = internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
-          dim, -1, 0, n_components_, Number>::integrate;
-        }
-      else
-        {
-          evaluate_funct = internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
-          dim, -1, 0, n_components_, Number>::evaluate;
-          integrate_funct = internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
-          dim, -1, 0, n_components_, Number>::integrate;
-        }
-    }
-  else
-    switch (this->data->element_type)
-      {
-      case internal::MatrixFreeFunctions::tensor_symmetric:
-        evaluate_funct =
-          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric,
-          dim, fe_degree_templ, n_q_points_1d_templ, n_components_,
-          Number>::evaluate;
-        integrate_funct =
-          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric,
-          dim, fe_degree_templ, n_q_points_1d_templ, n_components_,
-          Number>::integrate;
-        break;
-
-      case internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0:
-        evaluate_funct =
-          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
-          dim, fe_degree_templ, n_q_points_1d_templ, n_components_,
-          Number>::evaluate;
-        integrate_funct =
-          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
-          dim, fe_degree_templ, n_q_points_1d_templ, n_components_,
-          Number>::integrate;
-        break;
-
-      case internal::MatrixFreeFunctions::tensor_general:
-        evaluate_funct =
-          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
-          dim, fe_degree_templ, n_q_points_1d_templ, n_components_,
-          Number>::evaluate;
-        integrate_funct =
-          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
-          dim, fe_degree_templ, n_q_points_1d_templ, n_components_,
-          Number>::integrate;
-        break;
-
-      case internal::MatrixFreeFunctions::tensor_gausslobatto:
-        evaluate_funct =
-          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_gausslobatto,
-          dim, fe_degree_templ, n_q_points_1d_templ, n_components_,
-          Number>::evaluate;
-        integrate_funct =
-          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_gausslobatto,
-          dim, fe_degree_templ, n_q_points_1d_templ, n_components_,
-          Number>::integrate;
-        break;
-
-      case internal::MatrixFreeFunctions::truncated_tensor:
-        evaluate_funct =
-          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
-          dim, fe_degree_templ, n_q_points_1d_templ, n_components_,
-          Number>::evaluate;
-        integrate_funct =
-          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
-          dim, fe_degree_templ, n_q_points_1d_templ, n_components_,
-          Number>::integrate;
-        break;
-
-      default:
-        AssertThrow(false, ExcNotImplemented());
-      }
-
   (void)fe_no;
   (void)first_selected_component;
 
@@ -5484,11 +5372,89 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
   Assert(this->matrix_info != 0 ||
          this->mapped_geometry->is_initialized(), ExcNotInitialized());
 
-  // 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, this->scratch_data,
-                  evaluate_val, evaluate_grad, evaluate_lapl);
+  // Select algorithm matching the element type at run time
+  const unsigned int fe_degree_templ = fe_degree != -1 ? fe_degree : 0;
+  const unsigned int n_q_points_1d_templ = n_q_points_1d > 0 ? n_q_points_1d : 1;
+  if (fe_degree == -1)
+    {
+      if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
+        {
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
+                   dim, -1, 0, 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);
+        }
+      else if (this->data->element_type == internal::MatrixFreeFunctions::truncated_tensor)
+        {
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
+                   dim, -1, 0, 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);
+        }
+      else
+        {
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
+                   dim, -1, 0, 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);
+        }
+    }
+  else
+    {
+      if (fe_degree+1 == n_q_points_1d &&
+          this->data->element_type == internal::MatrixFreeFunctions::tensor_gausslobatto)
+        {
+          internal::FEEvaluationImplSpectral<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);
+        }
+      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>
+          ::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);
+        }
+      else if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric)
+        {
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric,
+                   dim, fe_degree_templ, n_q_points_1d_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);
+        }
+      else if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
+        {
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
+                   dim, fe_degree_templ, n_q_points_1d_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);
+        }
+      else if (this->data->element_type == internal::MatrixFreeFunctions::tensor_general)
+        {
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
+                   dim, fe_degree_templ, n_q_points_1d_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);
+        }
+      else if (this->data->element_type == internal::MatrixFreeFunctions::truncated_tensor)
+        {
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
+                   dim, fe_degree_templ, n_q_points_1d_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);
+        }
+      else
+        AssertThrow(false, ExcNotImplemented());
+    }
 
 #ifdef DEBUG
   if (evaluate_val == true)
@@ -5518,11 +5484,89 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
   Assert(this->matrix_info != 0 ||
          this->mapped_geometry->is_initialized(), ExcNotInitialized());
 
-  // 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, this->scratch_data,
-                   integrate_val, integrate_grad);
+  // Select algorithm matching the element type at run time
+  const unsigned int fe_degree_templ = fe_degree != -1 ? fe_degree : 0;
+  const unsigned int n_q_points_1d_templ = n_q_points_1d > 0 ? n_q_points_1d : 1;
+  if (fe_degree == -1)
+    {
+      if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
+        {
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
+                   dim, -1, 0, n_components_, Number>
+                   ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
+                               this->gradients_quad, this->scratch_data,
+                               integrate_val, integrate_grad);
+        }
+      else if (this->data->element_type == internal::MatrixFreeFunctions::truncated_tensor)
+        {
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
+                   dim, -1, 0, n_components_, Number>
+                   ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
+                               this->gradients_quad, this->scratch_data,
+                               integrate_val, integrate_grad);
+        }
+      else
+        {
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
+                   dim, -1, 0, n_components_, Number>
+                   ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
+                               this->gradients_quad, this->scratch_data,
+                               integrate_val, integrate_grad);
+        }
+    }
+  else
+    {
+      if (fe_degree+1 == n_q_points_1d &&
+          this->data->element_type == internal::MatrixFreeFunctions::tensor_gausslobatto)
+        {
+          internal::FEEvaluationImplSpectral<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);
+        }
+      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>
+          ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
+                      this->gradients_quad, this->scratch_data,
+                      integrate_val, integrate_grad);
+        }
+      else if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric)
+        {
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric,
+                   dim, fe_degree_templ, n_q_points_1d_templ, n_components_, Number>
+                   ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
+                               this->gradients_quad, this->scratch_data,
+                               integrate_val, integrate_grad);
+        }
+      else if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
+        {
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
+                   dim, fe_degree_templ, n_q_points_1d_templ, n_components_, Number>
+                   ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
+                               this->gradients_quad, this->scratch_data,
+                               integrate_val, integrate_grad);
+        }
+      else if (this->data->element_type == internal::MatrixFreeFunctions::tensor_general)
+        {
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
+                   dim, fe_degree_templ, n_q_points_1d_templ, n_components_, Number>
+                   ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
+                               this->gradients_quad, this->scratch_data,
+                               integrate_val, integrate_grad);
+        }
+      else if (this->data->element_type == internal::MatrixFreeFunctions::truncated_tensor)
+        {
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
+                   dim, fe_degree_templ, n_q_points_1d_templ, n_components_, Number>
+                   ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
+                               this->gradients_quad, this->scratch_data,
+                               integrate_val, integrate_grad);
+        }
+      else
+        AssertThrow(false, ExcNotImplemented());
+    }
 
 #ifdef DEBUG
   this->dof_values_initialized = true;
index 2cda877895c6cc8c7aaa3ca624a3f2c75abad95c..95b981f43baddace09399f1964516f140b9c4dfe 100644 (file)
@@ -151,10 +151,10 @@ namespace internal
       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.
+       * Stores the shape hessians of the spectral element space
+       * FE_DGQ<1>(Quadrature<1>) for faster evaluation in even-odd format.
        */
-      AlignedVector<VectorizedArray<Number> > shape_grad_spectral;
+      AlignedVector<VectorizedArray<Number> > shape_hessian_spectral_eo;
 
       /**
        * Stores the indices from cell DoFs to face DoFs. The rows go through
index eb6a2e96b653b2fe1b1ccbf152cb75ce86e17ca8..9d7b32818b53a80320622c4e0146888c2e3aa579 100644 (file)
@@ -209,23 +209,47 @@ 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))
         {
           if (check_1d_shapes_gausslobatto())
             element_type = tensor_gausslobatto;
           else
-            element_type = tensor_symmetric;
+            {
+              element_type = tensor_symmetric;
+              // get spectral evaluation points
+              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);
+                  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] =
+                          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] =
+                          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] =
+                          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] =
+                          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] =
+                          fe.shape_grad(fe_degree/2, quad.get_points()[q])[0];
+                        shape_hessian_spectral_eo[fe_degree/2*stride+q] =
+                          fe.shape_grad_grad(fe_degree/2, quad.get_points()[q])[0][0];
+                      }
+                }
+            }
         }
       else if (element_type == tensor_symmetric_plus_dg0)
         check_1d_shapes_symmetric(n_q_points_1d);
@@ -331,7 +355,6 @@ 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)
           {
@@ -355,16 +378,6 @@ 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)
@@ -375,9 +388,6 @@ 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;
@@ -434,8 +444,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 += MemoryConsumption::memory_consumption(shape_hessian_spectral_eo);
       memory += face_indices.memory_consumption();
       for (unsigned int i=0; i<2; ++i)
         {

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.