]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Evaluate and integrate hessians on faces 12593/head
authorMaximilian Bergbauer <bergbauer@lnm.mw.tum.de>
Thu, 17 Jun 2021 11:50:11 +0000 (13:50 +0200)
committerMaximilian Bergbauer <bergbauer@lnm.mw.tum.de>
Mon, 9 Aug 2021 07:26:37 +0000 (09:26 +0200)
16 files changed:
doc/news/changes/minor/20210728MaximilianBergbauer [new file with mode: 0644]
examples/step-76/step-76.cc
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/evaluation_selector.h
include/deal.II/matrix_free/evaluation_template_factory.h
include/deal.II/matrix_free/evaluation_template_factory.templates.h
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/mapping_info.templates.h
tests/matrix_free/interpolate_constant.output
tests/matrix_free/interpolate_cubic.output
tests/matrix_free/interpolate_functions_common.h
tests/matrix_free/interpolate_linear.output
tests/matrix_free/interpolate_quadratic.output
tests/matrix_free/interpolate_quadrature_01.cc
tests/matrix_free/matrix_vector_hessians_faces.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_hessians_faces.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20210728MaximilianBergbauer b/doc/news/changes/minor/20210728MaximilianBergbauer
new file mode 100644 (file)
index 0000000..91300d9
--- /dev/null
@@ -0,0 +1,4 @@
+Added: Evaluation and integration of hessians on faces is now available
+ in the matrix free framework.
+<br>
+(Maximilian Bergbauer, 2021/07/28)
index fb93b5ea495411c7882f7c3b1c6143c3030287f9..c7b2a8f987ec7b1ee2c93f37fe0bf8e416487d00 100644 (file)
@@ -794,6 +794,7 @@ namespace Euler_DG
                     buffer.data(),
                     phi_m.begin_values(),
                     false,
+                    false,
                     face);
 
                 // Check if the face is an internal or a boundary face and
@@ -896,6 +897,7 @@ namespace Euler_DG
                     phi_m.begin_values(),
                     phi.begin_values(),
                     false,
+                    false,
                     face);
               }
 
index 8a978f403d48521eaa73b7474fcf0e7e00d0fbe9..1438db4da614d556e55fccecf3f902759b673ba0 100644 (file)
@@ -131,6 +131,7 @@ namespace internal
               Number *                                      values_dofs_actual,
               Number *                                      values_quad,
               Number *                                      gradients_quad,
+              Number *                                      hessians_quad,
               Number *                                      scratch_data,
               const bool add_into_values_array);
   };
@@ -165,6 +166,7 @@ namespace internal
               Number *                                      values_dofs_actual,
               Number *                                      values_quad,
               Number *                                      gradients_quad,
+              Number *                                      hessians_quad,
               Number *                                      scratch_data,
               const bool add_into_values_array);
   };
@@ -459,9 +461,15 @@ namespace internal
     Number *                                      values_dofs_actual,
     Number *                                      values_quad,
     Number *                                      gradients_quad,
+    Number *                                      hessians_quad,
     Number *                                      scratch_data,
     const bool                                    add_into_values_array)
   {
+    // TODO: implement hessians
+    (void)hessians_quad;
+    AssertThrow(!(integration_flag & EvaluationFlags::hessians),
+                ExcNotImplemented());
+
     const EvaluatorVariant variant =
       EvaluatorSelector<type, (fe_degree + n_q_points_1d > 4)>::variant;
     using Eval = EvaluatorTensorProduct<variant,
@@ -764,9 +772,14 @@ namespace internal
                        Number *   values_dofs_actual,
                        Number *   values_quad,
                        Number *   gradients_quad,
+                       Number *   hessians_quad,
                        Number *   scratch_data,
                        const bool add_into_values_array)
   {
+    // TODO: implement hessians
+    (void)hessians_quad;
+    AssertThrow(!(integration_flag & EvaluationFlags::hessians),
+                ExcNotImplemented());
     (void)scratch_data;
 
     const unsigned int n_dofs     = shape_info.dofs_per_component_on_cell;
@@ -1263,6 +1276,7 @@ namespace internal
               Number *                                      values_dofs,
               Number *                                      values_quad,
               Number *                                      gradients_quad,
+              Number *                                      hessians_quad,
               Number *                                      scratch_data,
               const bool add_into_values_array);
   };
@@ -1356,9 +1370,16 @@ namespace internal
     Number *                                      values_dofs,
     Number *                                      values_quad,
     Number *                                      gradients_quad,
+    Number *                                      hessians_quad,
     Number *,
     const bool add_into_values_array)
   {
+    // TODO: implement hessians
+    (void)hessians_quad;
+    AssertThrow(!(integration_flag & EvaluationFlags::hessians),
+                ExcNotImplemented());
+
+
     AssertDimension(
       shape_info.data.front().shape_gradients_collocation_eo.size(),
       (fe_degree + 2) / 2 * (fe_degree + 1));
@@ -1440,6 +1461,7 @@ namespace internal
               Number *                                      values_dofs,
               Number *                                      values_quad,
               Number *                                      gradients_quad,
+              Number *                                      hessians_quad,
               Number *                                      scratch_data,
               const bool add_into_values_array);
   };
@@ -1517,6 +1539,7 @@ namespace internal
                        Number *values_dofs,
                        Number *values_quad,
                        Number *gradients_quad,
+                       Number *hessians_quad,
                        Number *,
                        const bool add_into_values_array)
   {
@@ -1536,11 +1559,13 @@ namespace internal
         if (integration_flag & EvaluationFlags::gradients)
           FEEvaluationImplCollocation<dim, n_q_points_1d - 1, Number>::
             integrate(1,
-                      integration_flag & EvaluationFlags::gradients,
+                      integration_flag & (EvaluationFlags::gradients |
+                                          EvaluationFlags::hessians),
                       shape_info,
                       values_quad,
                       nullptr,
                       gradients_quad,
+                      hessians_quad,
                       nullptr,
                       /*add_into_values_array=*/integration_flag &
                         EvaluationFlags::values);
@@ -1754,6 +1779,7 @@ namespace internal
         Number *   values_dofs_actual,
         Number *   values_quad,
         Number *   gradients_quad,
+        Number *   hessians_quad,
         Number *   scratch_data,
         const bool sum_into_values_array)
     {
@@ -1777,6 +1803,7 @@ namespace internal
                       values_dofs_actual,
                       values_quad,
                       gradients_quad,
+                      hessians_quad,
                       scratch_data,
                       sum_into_values_array);
         }
@@ -1796,6 +1823,7 @@ namespace internal
                                values_dofs_actual,
                                values_quad,
                                gradients_quad,
+                               hessians_quad,
                                scratch_data,
                                sum_into_values_array);
         }
@@ -1814,6 +1842,7 @@ namespace internal
                                values_dofs_actual,
                                values_quad,
                                gradients_quad,
+                               hessians_quad,
                                scratch_data,
                                sum_into_values_array);
         }
@@ -1831,6 +1860,7 @@ namespace internal
                                values_dofs_actual,
                                values_quad,
                                gradients_quad,
+                               hessians_quad,
                                scratch_data,
                                sum_into_values_array);
         }
@@ -1848,6 +1878,7 @@ namespace internal
                                values_dofs_actual,
                                values_quad,
                                gradients_quad,
+                               hessians_quad,
                                scratch_data,
                                sum_into_values_array);
         }
@@ -1864,6 +1895,7 @@ namespace internal
                                                         values_dofs_actual,
                                                         values_quad,
                                                         gradients_quad,
+                                                        hessians_quad,
                                                         scratch_data,
                                                         sum_into_values_array);
         }
@@ -1880,6 +1912,7 @@ namespace internal
                                values_dofs_actual,
                                values_quad,
                                gradients_quad,
+                               hessians_quad,
                                scratch_data,
                                sum_into_values_array);
         }
@@ -1913,37 +1946,52 @@ namespace internal
                      Number *                                      values_dofs,
                      Number *                                      values_quad,
                      Number *           gradients_quad,
+                     Number *           hessians_quad,
                      Number *           scratch_data,
                      const bool         evaluate_val,
                      const bool         evaluate_grad,
+                     const bool         evaluate_hessian,
                      const unsigned int subface_index)
     {
-      const AlignedVector<Number> &val1 =
+      const AlignedVector<Number> &val0 =
         symmetric_evaluate ?
           data.data.front().shape_values_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
              data.data.front().shape_values :
              data.data.front().values_within_subface[subface_index % 2]);
-      const AlignedVector<Number> &val2 =
+      const AlignedVector<Number> &val1 =
         symmetric_evaluate ?
           data.data.front().shape_values_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
              data.data.front().shape_values :
              data.data.front().values_within_subface[subface_index / 2]);
 
-      const AlignedVector<Number> &grad1 =
+      const AlignedVector<Number> &grad0 =
         symmetric_evaluate ?
           data.data.front().shape_gradients_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
              data.data.front().shape_gradients :
              data.data.front().gradients_within_subface[subface_index % 2]);
-      const AlignedVector<Number> &grad2 =
+      const AlignedVector<Number> &grad1 =
         symmetric_evaluate ?
           data.data.front().shape_gradients_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
              data.data.front().shape_gradients :
              data.data.front().gradients_within_subface[subface_index / 2]);
 
+      const AlignedVector<Number> &hessian0 =
+        symmetric_evaluate ?
+          data.data.front().shape_hessians_eo :
+          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
+             data.data.front().shape_hessians :
+             data.data.front().hessians_within_subface[subface_index % 2]);
+      const AlignedVector<Number> &hessian1 =
+        symmetric_evaluate ?
+          data.data.front().shape_hessians_eo :
+          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
+             data.data.front().shape_hessians :
+             data.data.front().hessians_within_subface[subface_index / 2]);
+
       using Eval =
         internal::EvaluatorTensorProduct<symmetric_evaluate ?
                                            internal::evaluate_evenodd :
@@ -1952,14 +2000,14 @@ namespace internal
                                          fe_degree + 1,
                                          n_q_points_1d,
                                          Number>;
-      Eval eval1(val1,
-                 grad1,
-                 AlignedVector<Number>(),
+      Eval eval0(val0,
+                 grad0,
+                 hessian0,
                  data.data.front().fe_degree + 1,
                  data.data.front().n_q_points_1d);
-      Eval eval2(val2,
-                 grad2,
-                 AlignedVector<Number>(),
+      Eval eval1(val1,
+                 grad1,
+                 hessian1,
                  data.data.front().fe_degree + 1,
                  data.data.front().n_q_points_1d);
 
@@ -1974,19 +2022,22 @@ namespace internal
                                         Utilities::pow(n_q_points_1d, dim - 1) :
                                         data.n_q_points_face;
 
-      if (evaluate_grad == false)
+      // keep a copy of the original pointer for the case of the Hessians
+      Number *values_dofs_ptr = values_dofs;
+
+      if (evaluate_val && !evaluate_grad)
         for (unsigned int c = 0; c < n_components; ++c)
           {
             switch (dim)
               {
                 case 3:
-                  eval1.template values<0, true, false>(values_dofs,
+                  eval0.template values<0, true, false>(values_dofs,
                                                         values_quad);
-                  eval2.template values<1, true, false>(values_quad,
+                  eval1.template values<1, true, false>(values_quad,
                                                         values_quad);
                   break;
                 case 2:
-                  eval1.template values<0, true, false>(values_dofs,
+                  eval0.template values<0, true, false>(values_dofs,
                                                         values_quad);
                   break;
                 case 1:
@@ -1995,10 +2046,10 @@ namespace internal
                 default:
                   Assert(false, ExcNotImplemented());
               }
-            values_dofs += 2 * size_deg;
+            values_dofs += 3 * size_deg;
             values_quad += n_q_points;
           }
-      else
+      else if (evaluate_grad)
         for (unsigned int c = 0; c < n_components; ++c)
           {
             switch (dim)
@@ -2006,9 +2057,9 @@ namespace internal
                 case 3:
                   if (use_collocation)
                     {
-                      eval1.template values<0, true, false>(values_dofs,
+                      eval0.template values<0, true, false>(values_dofs,
                                                             values_quad);
-                      eval1.template values<1, true, false>(values_quad,
+                      eval0.template values<1, true, false>(values_quad,
                                                             values_quad);
                       internal::EvaluatorTensorProduct<
                         internal::evaluate_evenodd,
@@ -2027,36 +2078,39 @@ namespace internal
                     }
                   else
                     {
-                      eval1.template gradients<0, true, false>(values_dofs,
+                      // grad x
+                      eval0.template gradients<0, true, false>(values_dofs,
                                                                scratch_data);
-                      eval2.template values<1, true, false>(scratch_data,
+                      eval1.template values<1, true, false>(scratch_data,
                                                             gradients_quad);
 
-                      eval1.template values<0, true, false>(values_dofs,
+                      // grad y
+                      eval0.template values<0, true, false>(values_dofs,
                                                             scratch_data);
-                      eval2.template gradients<1, true, false>(scratch_data,
+                      eval1.template gradients<1, true, false>(scratch_data,
                                                                gradients_quad +
                                                                  n_q_points);
 
                       if (evaluate_val == true)
-                        eval2.template values<1, true, false>(scratch_data,
+                        eval1.template values<1, true, false>(scratch_data,
                                                               values_quad);
                     }
-                  eval1.template values<0, true, false>(values_dofs + size_deg,
+                  // grad z
+                  eval0.template values<0, true, false>(values_dofs + size_deg,
                                                         scratch_data);
-                  eval2.template values<1, true, false>(
+                  eval1.template values<1, true, false>(
                     scratch_data, gradients_quad + (dim - 1) * n_q_points);
 
                   break;
                 case 2:
-                  eval1.template values<0, true, false>(values_dofs + size_deg,
+                  eval0.template values<0, true, false>(values_dofs + size_deg,
                                                         gradients_quad +
                                                           (dim - 1) *
                                                             n_q_points);
-                  eval1.template gradients<0, true, false>(values_dofs,
+                  eval0.template gradients<0, true, false>(values_dofs,
                                                            gradients_quad);
                   if (evaluate_val == true)
-                    eval1.template values<0, true, false>(values_dofs,
+                    eval0.template values<0, true, false>(values_dofs,
                                                           values_quad);
                   break;
                 case 1:
@@ -2066,10 +2120,86 @@ namespace internal
                 default:
                   AssertThrow(false, ExcNotImplemented());
               }
-            values_dofs += 2 * size_deg;
+            values_dofs += 3 * size_deg;
             values_quad += n_q_points;
             gradients_quad += dim * n_q_points;
           }
+
+      if (evaluate_hessian)
+        {
+          values_dofs                 = values_dofs_ptr;
+          constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
+          for (unsigned int c = 0; c < n_components; ++c)
+            {
+              switch (dim)
+                {
+                  case 3:
+                    // grad xx
+                    eval0.template hessians<0, true, false>(values_dofs,
+                                                            scratch_data);
+                    eval1.template values<1, true, false>(scratch_data,
+                                                          hessians_quad);
+
+                    // grad yy
+                    eval0.template values<0, true, false>(values_dofs,
+                                                          scratch_data);
+                    eval1.template hessians<1, true, false>(scratch_data,
+                                                            hessians_quad +
+                                                              n_q_points);
+
+                    // grad zz
+                    eval0.template values<0, true, false>(values_dofs +
+                                                            2 * size_deg,
+                                                          scratch_data);
+                    eval1.template values<1, true, false>(scratch_data,
+                                                          hessians_quad +
+                                                            2 * n_q_points);
+
+                    // grad xy
+                    eval0.template gradients<0, true, false>(values_dofs,
+                                                             scratch_data);
+                    eval1.template gradients<1, true, false>(scratch_data,
+                                                             hessians_quad +
+                                                               3 * n_q_points);
+
+                    // grad xz
+                    eval0.template gradients<0, true, false>(values_dofs +
+                                                               size_deg,
+                                                             scratch_data);
+                    eval1.template values<1, true, false>(scratch_data,
+                                                          hessians_quad +
+                                                            4 * n_q_points);
+
+                    // grad yz
+                    eval0.template values<0, true, false>(values_dofs +
+                                                            size_deg,
+                                                          scratch_data);
+                    eval1.template gradients<1, true, false>(scratch_data,
+                                                             hessians_quad +
+                                                               5 * n_q_points);
+
+                    break;
+                  case 2:
+                    // grad xx
+                    eval0.template hessians<0, true, false>(values_dofs,
+                                                            hessians_quad);
+                    // grad yy
+                    eval0.template values<0, true, false>(
+                      values_dofs + 2 * size_deg, hessians_quad + n_q_points);
+                    // grad xy
+                    eval0.template gradients<0, true, false>(
+                      values_dofs + size_deg, hessians_quad + 2 * n_q_points);
+                    break;
+                  case 1:
+                    hessians_quad[0] = values_dofs[2];
+                    break;
+                  default:
+                    AssertThrow(false, ExcNotImplemented());
+                }
+              values_dofs += 3 * size_deg;
+              hessians_quad += hdim * n_q_points;
+            }
+        }
     }
 
     static void
@@ -2078,37 +2208,52 @@ namespace internal
                       Number *                                      values_dofs,
                       Number *                                      values_quad,
                       Number *           gradients_quad,
+                      Number *           hessians_quad,
                       Number *           scratch_data,
                       const bool         integrate_val,
                       const bool         integrate_grad,
+                      const bool         integrate_hessian,
                       const unsigned int subface_index)
     {
-      const AlignedVector<Number> &val1 =
+      const AlignedVector<Number> &val0 =
         symmetric_evaluate ?
           data.data.front().shape_values_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
              data.data.front().shape_values :
              data.data.front().values_within_subface[subface_index % 2]);
-      const AlignedVector<Number> &val2 =
+      const AlignedVector<Number> &val1 =
         symmetric_evaluate ?
           data.data.front().shape_values_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
              data.data.front().shape_values :
              data.data.front().values_within_subface[subface_index / 2]);
 
-      const AlignedVector<Number> &grad1 =
+      const AlignedVector<Number> &grad0 =
         symmetric_evaluate ?
           data.data.front().shape_gradients_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
              data.data.front().shape_gradients :
              data.data.front().gradients_within_subface[subface_index % 2]);
-      const AlignedVector<Number> &grad2 =
+      const AlignedVector<Number> &grad1 =
         symmetric_evaluate ?
           data.data.front().shape_gradients_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
              data.data.front().shape_gradients :
              data.data.front().gradients_within_subface[subface_index / 2]);
 
+      const AlignedVector<Number> &hessian0 =
+        symmetric_evaluate ?
+          data.data.front().shape_hessians_eo :
+          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
+             data.data.front().shape_hessians :
+             data.data.front().hessians_within_subface[subface_index % 2]);
+      const AlignedVector<Number> &hessian1 =
+        symmetric_evaluate ?
+          data.data.front().shape_hessians_eo :
+          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
+             data.data.front().shape_hessians :
+             data.data.front().hessians_within_subface[subface_index / 2]);
+
       using Eval =
         internal::EvaluatorTensorProduct<symmetric_evaluate ?
                                            internal::evaluate_evenodd :
@@ -2117,14 +2262,14 @@ namespace internal
                                          fe_degree + 1,
                                          n_q_points_1d,
                                          Number>;
-      Eval eval1(val1,
-                 grad1,
-                 val1,
+      Eval eval0(val0,
+                 grad0,
+                 hessian0,
                  data.data.front().fe_degree + 1,
                  data.data.front().n_q_points_1d);
-      Eval eval2(val2,
-                 grad2,
-                 val1,
+      Eval eval1(val1,
+                 grad1,
+                 hessian1,
                  data.data.front().fe_degree + 1,
                  data.data.front().n_q_points_1d);
 
@@ -2139,19 +2284,22 @@ namespace internal
                                         Utilities::pow(n_q_points_1d, dim - 1) :
                                         data.n_q_points_face;
 
-      if (integrate_grad == false)
+      // keep a copy of the original pointer for the case of the Hessians
+      Number *values_dofs_ptr = values_dofs;
+
+      if (integrate_val && !integrate_grad)
         for (unsigned int c = 0; c < n_components; ++c)
           {
             switch (dim)
               {
                 case 3:
-                  eval2.template values<1, false, false>(values_quad,
+                  eval1.template values<1, false, false>(values_quad,
                                                          values_quad);
-                  eval1.template values<0, false, false>(values_quad,
+                  eval0.template values<0, false, false>(values_quad,
                                                          values_dofs);
                   break;
                 case 2:
-                  eval1.template values<0, false, false>(values_quad,
+                  eval0.template values<0, false, false>(values_quad,
                                                          values_dofs);
                   break;
                 case 1:
@@ -2160,20 +2308,21 @@ namespace internal
                 default:
                   Assert(false, ExcNotImplemented());
               }
-            values_dofs += 2 * size_deg;
+            values_dofs += 3 * size_deg;
             values_quad += n_q_points;
           }
-      else
+      else if (integrate_grad)
         for (unsigned int c = 0; c < n_components; ++c)
           {
             switch (dim)
               {
                 case 3:
-                  eval2.template values<1, false, false>(gradients_quad +
+                  // grad z
+                  eval1.template values<1, false, false>(gradients_quad +
                                                            2 * n_q_points,
                                                          gradients_quad +
                                                            2 * n_q_points);
-                  eval1.template values<0, false, false>(
+                  eval0.template values<0, false, false>(
                     gradients_quad + 2 * n_q_points, values_dofs + size_deg);
                   if (use_collocation)
                     {
@@ -2195,39 +2344,42 @@ namespace internal
                           gradients_quad + n_q_points, values_quad);
                       eval_grad.template gradients<0, false, true>(
                         gradients_quad, values_quad);
-                      eval1.template values<1, false, false>(values_quad,
+                      eval0.template values<1, false, false>(values_quad,
                                                              values_quad);
-                      eval1.template values<0, false, false>(values_quad,
+                      eval0.template values<0, false, false>(values_quad,
                                                              values_dofs);
                     }
                   else
                     {
                       if (integrate_val)
                         {
-                          eval2.template values<1, false, false>(values_quad,
+                          eval1.template values<1, false, false>(values_quad,
                                                                  scratch_data);
-                          eval2.template gradients<1, false, true>(
+                          eval1.template gradients<1, false, true>(
                             gradients_quad + n_q_points, scratch_data);
                         }
                       else
-                        eval2.template gradients<1, false, false>(
+                        eval1.template gradients<1, false, false>(
                           gradients_quad + n_q_points, scratch_data);
 
-                      eval1.template values<0, false, false>(scratch_data,
+                      // grad y
+                      eval0.template values<0, false, false>(scratch_data,
                                                              values_dofs);
-                      eval2.template values<1, false, false>(gradients_quad,
+
+                      // grad x
+                      eval1.template values<1, false, false>(gradients_quad,
                                                              scratch_data);
-                      eval1.template gradients<0, false, true>(scratch_data,
+                      eval0.template gradients<0, false, true>(scratch_data,
                                                                values_dofs);
                     }
                   break;
                 case 2:
-                  eval1.template values<0, false, false>(
+                  eval0.template values<0, false, false>(
                     gradients_quad + n_q_points, values_dofs + size_deg);
-                  eval1.template gradients<0, false, false>(gradients_quad,
+                  eval0.template gradients<0, false, false>(gradients_quad,
                                                             values_dofs);
                   if (integrate_val == true)
-                    eval1.template values<0, false, true>(values_quad,
+                    eval0.template values<0, false, true>(values_quad,
                                                           values_dofs);
                   break;
                 case 1:
@@ -2237,10 +2389,108 @@ namespace internal
                 default:
                   AssertThrow(false, ExcNotImplemented());
               }
-            values_dofs += 2 * size_deg;
+            values_dofs += 3 * size_deg;
             values_quad += n_q_points;
             gradients_quad += dim * n_q_points;
           }
+
+      if (integrate_hessian)
+        {
+          values_dofs                 = values_dofs_ptr;
+          constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
+          for (unsigned int c = 0; c < n_components; ++c)
+            {
+              switch (dim)
+                {
+                  case 3:
+                    // grad xx
+                    eval1.template values<1, false, false>(hessians_quad,
+                                                           scratch_data);
+                    if (integrate_val || integrate_grad)
+                      eval0.template hessians<0, false, true>(scratch_data,
+                                                              values_dofs);
+                    else
+                      eval0.template hessians<0, false, false>(scratch_data,
+                                                               values_dofs);
+
+                    // grad yy
+                    eval1.template hessians<1, false, false>(hessians_quad +
+                                                               n_q_points,
+                                                             scratch_data);
+                    eval0.template values<0, false, true>(scratch_data,
+                                                          values_dofs);
+
+                    // grad zz
+                    eval1.template values<1, false, false>(hessians_quad +
+                                                             2 * n_q_points,
+                                                           scratch_data);
+                    eval0.template values<0, false, false>(scratch_data,
+                                                           values_dofs +
+                                                             2 * size_deg);
+
+                    // grad xy
+                    eval1.template gradients<1, false, false>(hessians_quad +
+                                                                3 * n_q_points,
+                                                              scratch_data);
+                    eval0.template gradients<0, false, true>(scratch_data,
+                                                             values_dofs);
+
+                    // grad xz
+                    eval1.template values<1, false, false>(hessians_quad +
+                                                             4 * n_q_points,
+                                                           scratch_data);
+                    if (integrate_grad)
+                      eval0.template gradients<0, false, true>(scratch_data,
+                                                               values_dofs +
+                                                                 size_deg);
+                    else
+                      eval0.template gradients<0, false, false>(scratch_data,
+                                                                values_dofs +
+                                                                  size_deg);
+
+                    // grad yz
+                    eval1.template gradients<1, false, false>(hessians_quad +
+                                                                5 * n_q_points,
+                                                              scratch_data);
+                    eval0.template values<0, false, true>(scratch_data,
+                                                          values_dofs +
+                                                            size_deg);
+
+                    break;
+                  case 2:
+                    // grad xx
+                    if (integrate_val || integrate_grad)
+                      eval0.template hessians<0, false, true>(hessians_quad,
+                                                              values_dofs);
+                    else
+                      eval0.template hessians<0, false, false>(hessians_quad,
+                                                               values_dofs);
+
+                    // grad yy
+                    eval0.template values<0, false, false>(
+                      hessians_quad + n_q_points, values_dofs + 2 * size_deg);
+                    // grad xy
+                    if (integrate_grad)
+                      eval0.template gradients<0, false, true>(
+                        hessians_quad + 2 * n_q_points, values_dofs + size_deg);
+                    else
+                      eval0.template gradients<0, false, false>(
+                        hessians_quad + 2 * n_q_points, values_dofs + size_deg);
+                    break;
+                  case 1:
+                    values_dofs[2] = hessians_quad[0];
+                    if (!integrate_val)
+                      values_dofs[0] = 0;
+                    if (!integrate_grad)
+                      values_dofs[1] = 0;
+                    break;
+                  default:
+                    AssertThrow(false, ExcNotImplemented());
+                }
+              values_dofs += 3 * size_deg;
+              hessians_quad += hdim * n_q_points;
+            }
+        }
     }
   };
 
@@ -2256,6 +2506,7 @@ namespace internal
                 const Number *                                input,
                 Number *                                      output,
                 const bool                                    do_gradients,
+                const bool                                    do_hessians,
                 const unsigned int                            face_no)
     {
       Assert(static_cast<unsigned int>(fe_degree) ==
@@ -2268,11 +2519,12 @@ namespace internal
         input,
         output,
         do_gradients,
+        do_hessians,
         face_no,
         data.data.front().fe_degree + 1,
         data.data.front().shape_data_on_face,
         data.dofs_per_component_on_cell,
-        2 * data.dofs_per_component_on_face);
+        3 * data.dofs_per_component_on_face);
     }
 
     /**
@@ -2285,6 +2537,7 @@ namespace internal
                            const Number *                                input,
                            Number *                                      output,
                            const bool         do_gradients,
+                           const bool         do_hessians,
                            const unsigned int face_no)
     {
       Assert(static_cast<unsigned int>(fe_degree + 1) ==
@@ -2297,6 +2550,7 @@ namespace internal
         input,
         output,
         do_gradients,
+        do_hessians,
         face_no,
         data.data.front().quadrature.size(),
         data.data.front().quadrature_data_on_face,
@@ -2311,6 +2565,7 @@ namespace internal
                         const Number *     input,
                         Number *           output,
                         const bool         do_gradients,
+                        const bool         do_hessians,
                         const unsigned int face_no,
                         const unsigned int n_points_1d,
                         const std::array<AlignedVector<Number>, 2> &shape_data,
@@ -2339,7 +2594,13 @@ namespace internal
 
           for (unsigned int c = 0; c < n_components; c++)
             {
-              if (do_gradients)
+              if (do_hessians)
+                evalf.template apply_face<face_direction,
+                                          do_evaluate,
+                                          add_into_output,
+                                          2,
+                                          lex_faces>(input, output);
+              else if (!do_hessians && do_gradients)
                 evalf.template apply_face<face_direction,
                                           do_evaluate,
                                           add_into_output,
@@ -2364,6 +2625,7 @@ namespace internal
             input,
             output,
             do_gradients,
+            do_hessians,
             face_no,
             n_points_1d,
             shape_data,
@@ -2489,10 +2751,12 @@ namespace internal
                               const bool                    integrate,
                               const bool                    values,
                               const bool                    gradients,
+                              const bool                    hessians,
                               const unsigned int            n_q_points,
                               Number *                      tmp_values,
                               Number *                      values_quad,
-                              Number *                      gradients_quad)
+                              Number *                      gradients_quad,
+                              Number *                      hessians_quad)
   {
     Assert(face_orientation, ExcInternalError());
     const unsigned int *orientation = &orientation_map[face_orientation][0];
@@ -2523,6 +2787,24 @@ namespace internal
               for (unsigned int q = 0; q < n_q_points; ++q)
                 gradients_quad[(c * dim + d) * n_q_points + q] = tmp_values[q];
             }
+        if (hessians == true)
+          {
+            const unsigned int hdim = (dim * (dim + 1)) / 2;
+            for (unsigned int d = 0; d < hdim; ++d)
+              {
+                if (integrate)
+                  for (unsigned int q = 0; q < n_q_points; ++q)
+                    tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
+                                                  orientation[q]];
+                else
+                  for (unsigned int q = 0; q < n_q_points; ++q)
+                    tmp_values[orientation[q]] =
+                      hessians_quad[(c * hdim + d) * n_q_points + q];
+                for (unsigned int q = 0; q < n_q_points; ++q)
+                  hessians_quad[(c * hdim + d) * n_q_points + q] =
+                    tmp_values[q];
+              }
+          }
       }
   }
 
@@ -2538,9 +2820,11 @@ namespace internal
         const VectorizedArrayType *                                values_array,
         VectorizedArrayType *                                      values_quad,
         VectorizedArrayType *         gradients_quad,
+        VectorizedArrayType *         hessians_quad,
         VectorizedArrayType *         scratch_data,
         const bool                    evaluate_values,
         const bool                    evaluate_gradients,
+        const bool                    evaluate_hessians,
         const unsigned int            face_no,
         const unsigned int            subface_index,
         const unsigned int            face_orientation,
@@ -2607,6 +2891,10 @@ namespace internal
                 }
             }
 
+          if (evaluate_hessians)
+            {
+              AssertThrow(false, ExcNotImplemented());
+            }
 
           return true;
         }
@@ -2622,8 +2910,13 @@ namespace internal
       VectorizedArrayType *temp1 = scratch_data;
 
       FEFaceNormalEvaluationImpl<dim, fe_degree, VectorizedArrayType>::
-        template interpolate<true, false>(
-          n_components, data, values_array, temp1, evaluate_gradients, face_no);
+        template interpolate<true, false>(n_components,
+                                          data,
+                                          values_array,
+                                          temp1,
+                                          evaluate_gradients,
+                                          evaluate_hessians,
+                                          face_no);
 
       const unsigned int n_q_points_1d_actual =
         fe_degree > -1 ? n_q_points_1d : 0;
@@ -2640,11 +2933,13 @@ namespace internal
                                                  temp1,
                                                  values_quad,
                                                  gradients_quad,
-                                                 scratch_data + 2 *
+                                                 hessians_quad,
+                                                 scratch_data + 3 *
                                                                   n_components *
                                                                   dofs_per_face,
                                                  evaluate_values,
                                                  evaluate_gradients,
+                                                 evaluate_hessians,
                                                  subface_index);
       else
         FEFaceEvaluationImpl<
@@ -2657,11 +2952,13 @@ namespace internal
                                                  temp1,
                                                  values_quad,
                                                  gradients_quad,
-                                                 scratch_data + 2 *
+                                                 hessians_quad,
+                                                 scratch_data + 3 *
                                                                   n_components *
                                                                   dofs_per_face,
                                                  evaluate_values,
                                                  evaluate_gradients,
+                                                 evaluate_hessians,
                                                  subface_index);
 
       if (face_orientation)
@@ -2672,10 +2969,12 @@ namespace internal
                                     false,
                                     evaluate_values,
                                     evaluate_gradients,
+                                    evaluate_hessians,
                                     data.n_q_points_face,
                                     scratch_data,
                                     values_quad,
-                                    gradients_quad);
+                                    gradients_quad,
+                                    hessians_quad);
 
       return false;
     }
@@ -2693,9 +2992,11 @@ namespace internal
         VectorizedArrayType *                                      values_array,
         VectorizedArrayType *                                      values_quad,
         VectorizedArrayType *         gradients_quad,
+        VectorizedArrayType *         hessians_quad,
         VectorizedArrayType *         scratch_data,
         const bool                    integrate_values,
         const bool                    integrate_gradients,
+        const bool                    integrate_hessians,
         const unsigned int            face_no,
         const unsigned int            subface_index,
         const unsigned int            face_orientation,
@@ -2766,6 +3067,10 @@ namespace internal
                 }
             }
 
+          if (integrate_hessians)
+            {
+              AssertThrow(false, ExcNotImplemented());
+            }
 
           return true;
         }
@@ -2778,10 +3083,12 @@ namespace internal
                                     true,
                                     integrate_values,
                                     integrate_gradients,
+                                    integrate_hessians,
                                     data.n_q_points_face,
                                     scratch_data,
                                     values_quad,
-                                    gradients_quad);
+                                    gradients_quad,
+                                    hessians_quad);
 
       constexpr unsigned int static_dofs_per_face =
         fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
@@ -2808,11 +3115,13 @@ namespace internal
                                                   temp1,
                                                   values_quad,
                                                   gradients_quad,
+                                                  hessians_quad,
                                                   scratch_data +
-                                                    2 * n_components *
+                                                    3 * n_components *
                                                       dofs_per_face,
                                                   integrate_values,
                                                   integrate_gradients,
+                                                  integrate_hessians,
                                                   subface_index);
       else
         FEFaceEvaluationImpl<
@@ -2825,11 +3134,13 @@ namespace internal
                                                   temp1,
                                                   values_quad,
                                                   gradients_quad,
+                                                  hessians_quad,
                                                   scratch_data +
-                                                    2 * n_components *
+                                                    3 * n_components *
                                                       dofs_per_face,
                                                   integrate_values,
                                                   integrate_gradients,
+                                                  integrate_hessians,
                                                   subface_index);
 
       FEFaceNormalEvaluationImpl<dim, fe_degree, VectorizedArrayType>::
@@ -2838,6 +3149,7 @@ namespace internal
                                            temp1,
                                            values_array,
                                            integrate_gradients,
+                                           integrate_hessians,
                                            face_no);
       return false;
     }
@@ -2857,9 +3169,11 @@ namespace internal
     auto &dof_info                 = proc.dof_info;
     auto  values_quad              = proc.values_quad;
     auto  gradients_quad           = proc.gradients_quad;
+    auto  hessians_quad            = proc.hessians_quad;
     auto  scratch_data             = proc.scratch_data;
     auto  do_values                = proc.do_values;
     auto  do_gradients             = proc.do_gradients;
+    auto  do_hessians              = proc.do_hessians;
     auto  active_fe_index          = proc.active_fe_index;
     auto  first_selected_component = proc.first_selected_component;
     auto  cells                    = proc.cells;
@@ -2915,10 +3229,12 @@ namespace internal
                                     true,
                                     do_values,
                                     do_gradients,
+                                    do_hessians,
                                     data.n_q_points_face,
                                     scratch_data,
                                     values_quad,
-                                    gradients_quad);
+                                    gradients_quad,
+                                    hessians_quad);
       }
 
     // we know that the gradient weights for the Hermite case on the
@@ -3337,10 +3653,11 @@ namespace internal
               }
 
             // case 4: contiguous indices without interleaving
-            else if (n_face_orientations > 1 ||
-                     dof_info.index_storage_variants[dof_access_index][cell] ==
-                       MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-                         contiguous)
+            else if ((n_face_orientations > 1 ||
+                      dof_info.index_storage_variants[dof_access_index][cell] ==
+                        MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                          contiguous) &&
+                     !do_hessians)
               {
                 const unsigned int *indices =
                   &dof_info.dof_indices_contiguous[dof_access_index]
@@ -3562,10 +3879,12 @@ namespace internal
                                     false,
                                     do_values,
                                     do_gradients,
+                                    do_hessians,
                                     data.n_q_points_face,
                                     scratch_data,
                                     values_quad,
-                                    gradients_quad);
+                                    gradients_quad,
+                                    hessians_quad);
       }
 
     return accesses_global_vector;
@@ -3589,9 +3908,11 @@ namespace internal
         const MatrixFreeFunctions::DoFInfo &                       dof_info,
         VectorizedArrayType *                                      values_quad,
         VectorizedArrayType *gradients_quad,
+        VectorizedArrayType *hessians_quad,
         VectorizedArrayType *scratch_data,
         const bool           evaluate_values,
         const bool           evaluate_gradients,
+        const bool           evaluate_hessians,
         const unsigned int   active_fe_index,
         const unsigned int   first_selected_component,
         const std::array<unsigned int, VectorizedArrayType::size()> cells,
@@ -3618,9 +3939,11 @@ namespace internal
                                             dof_info,
                                             values_quad,
                                             gradients_quad,
+                                            hessians_quad,
                                             scratch_data,
                                             evaluate_values,
                                             evaluate_gradients,
+                                            evaluate_hessians,
                                             active_fe_index,
                                             first_selected_component,
                                             cells,
@@ -3657,9 +3980,11 @@ namespace internal
         const MatrixFreeFunctions::DoFInfo &                       dof_info,
         VectorizedArrayType *                                      values_quad,
         VectorizedArrayType *gradients_quad,
+        VectorizedArrayType *hessians_quad,
         VectorizedArrayType *scratch_data,
         const bool           do_values,
         const bool           do_gradients,
+        const bool           do_hessians,
         const unsigned int   active_fe_index,
         const unsigned int   first_selected_component,
         const std::array<unsigned int, VectorizedArrayType::size()> cells,
@@ -3677,9 +4002,11 @@ namespace internal
         , dof_info(dof_info)
         , values_quad(values_quad)
         , gradients_quad(gradients_quad)
+        , hessians_quad(hessians_quad)
         , scratch_data(scratch_data)
         , do_values(do_values)
         , do_gradients(do_gradients)
+        , do_hessians(do_hessians)
         , active_fe_index(active_fe_index)
         , first_selected_component(first_selected_component)
         , cells(cells)
@@ -3771,6 +4098,7 @@ namespace internal
         const unsigned int n_q_points =
           fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
                            data.n_q_points_face;
+        constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
         if (fe_degree > -1 &&
             subface_index >= GeometryInfo<dim>::max_children_per_cell &&
             data.element_type <= MatrixFreeFunctions::tensor_symmetric)
@@ -3784,9 +4112,11 @@ namespace internal
                              temp1,
                              values_quad + comp * n_q_points,
                              gradients_quad + comp * dim * n_q_points,
-                             scratch_data + 2 * dofs_per_face,
+                             hessians_quad + comp * hdim * n_q_points,
+                             scratch_data + 3 * dofs_per_face,
                              do_values,
                              do_gradients,
+                             do_hessians,
                              subface_index);
         else
           FEFaceEvaluationImpl<false,
@@ -3799,9 +4129,11 @@ namespace internal
                              temp1,
                              values_quad + comp * n_q_points,
                              gradients_quad + comp * dim * n_q_points,
-                             scratch_data + 2 * dofs_per_face,
+                             hessians_quad + comp * hdim * n_q_points,
+                             scratch_data + 3 * dofs_per_face,
                              do_values,
                              do_gradients,
+                             do_hessians,
                              subface_index);
       }
 
@@ -3813,9 +4145,11 @@ namespace internal
       const MatrixFreeFunctions::DoFInfo &                       dof_info;
       VectorizedArrayType *                                      values_quad;
       VectorizedArrayType *                                      gradients_quad;
+      VectorizedArrayType *                                      hessians_quad;
       VectorizedArrayType *                                      scratch_data;
       const bool                                                 do_values;
       const bool                                                 do_gradients;
+      const bool                                                 do_hessians;
       const unsigned int active_fe_index;
       const unsigned int first_selected_component;
       const std::array<unsigned int, VectorizedArrayType::size()> cells;
@@ -3845,9 +4179,11 @@ namespace internal
         VectorizedArrayType *                                      values_array,
         VectorizedArrayType *                                      values_quad,
         VectorizedArrayType *gradients_quad,
+        VectorizedArrayType *hessians_quad,
         VectorizedArrayType *scratch_data,
         const bool           integrate_values,
         const bool           integrate_gradients,
+        const bool           integrate_hessians,
         const unsigned int   active_fe_index,
         const unsigned int   first_selected_component,
         const std::array<unsigned int, VectorizedArrayType::size()> cells,
@@ -3872,9 +4208,11 @@ namespace internal
                                                    values_array,
                                                    values_quad,
                                                    gradients_quad,
+                                                   hessians_quad,
                                                    scratch_data,
                                                    integrate_values,
                                                    integrate_gradients,
+                                                   integrate_hessians,
                                                    face_nos[0],
                                                    subface_index,
                                                    face_orientations[0],
@@ -3894,9 +4232,11 @@ namespace internal
                                             dof_info,
                                             values_quad,
                                             gradients_quad,
+                                            hessians_quad,
                                             scratch_data,
                                             integrate_values,
                                             integrate_gradients,
+                                            integrate_hessians,
                                             active_fe_index,
                                             first_selected_component,
                                             cells,
@@ -3935,9 +4275,11 @@ namespace internal
         const MatrixFreeFunctions::DoFInfo &                       dof_info,
         VectorizedArrayType *                                      values_quad,
         VectorizedArrayType *gradients_quad,
+        VectorizedArrayType *hessians_quad,
         VectorizedArrayType *scratch_data,
         const bool           do_values,
         const bool           do_gradients,
+        const bool           do_hessians,
         const unsigned int   active_fe_index,
         const unsigned int   first_selected_component,
         const std::array<unsigned int, VectorizedArrayType::size()> cells,
@@ -3956,9 +4298,11 @@ namespace internal
         , dof_info(dof_info)
         , values_quad(values_quad)
         , gradients_quad(gradients_quad)
+        , hessians_quad(hessians_quad)
         , scratch_data(scratch_data)
         , do_values(do_values)
         , do_gradients(do_gradients)
+        , do_hessians(do_hessians)
         , active_fe_index(active_fe_index)
         , first_selected_component(first_selected_component)
         , cells(cells)
@@ -4054,6 +4398,7 @@ namespace internal
             temp1,
             values_array + comp * data.dofs_per_component_on_cell,
             do_gradients,
+            do_hessians,
             face_nos[0]);
       }
 
@@ -4068,6 +4413,7 @@ namespace internal
         const unsigned int n_q_points =
           fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
                            data.n_q_points_face;
+        constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
         if (fe_degree > -1 &&
             subface_index >= GeometryInfo<dim>::max_children_per_cell &&
             data.element_type <=
@@ -4082,9 +4428,11 @@ namespace internal
                               temp1,
                               values_quad + comp * n_q_points,
                               gradients_quad + dim * comp * n_q_points,
-                              scratch_data + 2 * dofs_per_face,
+                              hessians_quad + hdim * comp * n_q_points,
+                              scratch_data + 3 * dofs_per_face,
                               do_values,
                               do_gradients,
+                              do_hessians,
                               subface_index);
         else
           internal::FEFaceEvaluationImpl<false,
@@ -4097,9 +4445,11 @@ namespace internal
                               temp1,
                               values_quad + comp * n_q_points,
                               gradients_quad + dim * comp * n_q_points,
-                              scratch_data + 2 * dofs_per_face,
+                              hessians_quad + hdim * comp * n_q_points,
+                              scratch_data + 3 * dofs_per_face,
                               do_values,
                               do_gradients,
+                              do_hessians,
                               subface_index);
       }
 
@@ -4114,9 +4464,11 @@ namespace internal
       const MatrixFreeFunctions::DoFInfo &                       dof_info;
       VectorizedArrayType *                                      values_quad;
       VectorizedArrayType *                                      gradients_quad;
+      VectorizedArrayType *                                      hessians_quad;
       VectorizedArrayType *                                      scratch_data;
       const bool                                                 do_values;
       const bool                                                 do_gradients;
+      const bool                                                 do_hessians;
       const unsigned int active_fe_index;
       const unsigned int first_selected_component;
       const std::array<unsigned int, VectorizedArrayType::size()> cells;
index 990263531ba60000f3201b85eac6a989ae41e845..62fc91bcd0d80b10afad5ced814689f0c04404b5 100644 (file)
@@ -64,6 +64,7 @@ struct SelectEvaluator
             Number *   values_dofs_actual,
             Number *   values_quad,
             Number *   gradients_quad,
+            Number *   hessians_quad,
             Number *   scratch_data,
             const bool sum_into_values_array = false);
 };
@@ -107,6 +108,7 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, Number>::integrate(
   Number *                                                values_dofs_actual,
   Number *                                                values_quad,
   Number *                                                gradients_quad,
+  Number *                                                hessians_quad,
   Number *                                                scratch_data,
   const bool                                              sum_into_values_array)
 {
@@ -119,6 +121,7 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, Number>::integrate(
                                            values_dofs_actual,
                                            values_quad,
                                            gradients_quad,
+                                           hessians_quad,
                                            scratch_data,
                                            sum_into_values_array);
 }
index 9299bd4921ec6dacf53fd3cf810ca4ab290f56e4..c496ad4edb36459b79c2d7075569d108d6b7692e 100644 (file)
@@ -58,6 +58,7 @@ namespace internal
       VectorizedArrayType *values_dofs_actual,
       VectorizedArrayType *values_quad,
       VectorizedArrayType *gradients_quad,
+      VectorizedArrayType *hessians_quad,
       VectorizedArrayType *scratch_data,
       const bool           sum_into_values_array);
 
@@ -79,9 +80,11 @@ namespace internal
              const VectorizedArrayType *   values_array,
              VectorizedArrayType *         values_quad,
              VectorizedArrayType *         gradients_quad,
+             VectorizedArrayType *         hessians_quad,
              VectorizedArrayType *         scratch_data,
              const bool                    evaluate_values,
              const bool                    evaluate_gradients,
+             const bool                    evaluate_hessians,
              const unsigned int            face_no,
              const unsigned int            subface_index,
              const unsigned int            face_orientation,
@@ -93,9 +96,11 @@ namespace internal
               VectorizedArrayType *         values_array,
               VectorizedArrayType *         values_quad,
               VectorizedArrayType *         gradients_quad,
+              VectorizedArrayType *         hessians_quad,
               VectorizedArrayType *         scratch_data,
               const bool                    integrate_values,
               const bool                    integrate_gradients,
+              const bool                    integrate_hessians,
               const unsigned int            face_no,
               const unsigned int            subface_index,
               const unsigned int            face_orientation,
@@ -111,9 +116,11 @@ namespace internal
       const MatrixFreeFunctions::DoFInfo &                       dof_info,
       VectorizedArrayType *                                      values_quad,
       VectorizedArrayType *                                      gradients_quad,
+      VectorizedArrayType *                                      hessians_quad,
       VectorizedArrayType *                                      scratch_data,
       const bool         evaluate_values,
       const bool         evaluate_gradients,
+      const bool         evaluate_hessians,
       const unsigned int active_fe_index,
       const unsigned int first_selected_component,
       const std::array<unsigned int, VectorizedArrayType::size()> cells,
@@ -135,9 +142,11 @@ namespace internal
       VectorizedArrayType *                                      values_array,
       VectorizedArrayType *                                      values_quad,
       VectorizedArrayType *                                      gradients_quad,
+      VectorizedArrayType *                                      hessians_quad,
       VectorizedArrayType *                                      scratch_data,
       const bool         integrate_values,
       const bool         integrate_gradients,
+      const bool         integrate_hessians,
       const unsigned int active_fe_index,
       const unsigned int first_selected_component,
       const std::array<unsigned int, VectorizedArrayType::size()> cells,
index 027aa939997c18e7b82217575e4a5ed007a31f93..4d26a68279497bf136ce2dae87ed70f795cac278 100644 (file)
@@ -113,6 +113,7 @@ namespace internal
     VectorizedArrayType *values_dofs_actual,
     VectorizedArrayType *values_quad,
     VectorizedArrayType *gradients_quad,
+    VectorizedArrayType *hessians_quad,
     VectorizedArrayType *scratch_data,
     const bool           sum_into_values_array)
   {
@@ -127,6 +128,7 @@ namespace internal
       values_dofs_actual,
       values_quad,
       gradients_quad,
+      hessians_quad,
       scratch_data,
       sum_into_values_array);
   }
@@ -153,9 +155,11 @@ namespace internal
     const VectorizedArrayType *                                values_array,
     VectorizedArrayType *                                      values_quad,
     VectorizedArrayType *                                      gradients_quad,
+    VectorizedArrayType *                                      hessians_quad,
     VectorizedArrayType *                                      scratch_data,
     const bool                                                 evaluate_values,
     const bool                    evaluate_gradients,
+    const bool                    evaluate_hessians,
     const unsigned int            face_no,
     const unsigned int            subface_index,
     const unsigned int            face_orientation,
@@ -171,9 +175,11 @@ namespace internal
       values_array,
       values_quad,
       gradients_quad,
+      hessians_quad,
       scratch_data,
       evaluate_values,
       evaluate_gradients,
+      evaluate_hessians,
       face_no,
       subface_index,
       face_orientation,
@@ -190,9 +196,11 @@ namespace internal
     VectorizedArrayType *                                      values_array,
     VectorizedArrayType *                                      values_quad,
     VectorizedArrayType *                                      gradients_quad,
+    VectorizedArrayType *                                      hessians_quad,
     VectorizedArrayType *                                      scratch_data,
     const bool                                                 integrate_values,
     const bool                    integrate_gradients,
+    const bool                    integrate_hessians,
     const unsigned int            face_no,
     const unsigned int            subface_index,
     const unsigned int            face_orientation,
@@ -208,9 +216,11 @@ namespace internal
       values_array,
       values_quad,
       gradients_quad,
+      hessians_quad,
       scratch_data,
       integrate_values,
       integrate_gradients,
+      integrate_hessians,
       face_no,
       subface_index,
       face_orientation,
@@ -230,9 +240,11 @@ namespace internal
     const MatrixFreeFunctions::DoFInfo &                       dof_info,
     VectorizedArrayType *                                      values_quad,
     VectorizedArrayType *                                      gradients_quad,
+    VectorizedArrayType *                                      hessians_quad,
     VectorizedArrayType *                                      scratch_data,
     const bool                                                 evaluate_values,
     const bool         evaluate_gradients,
+    const bool         evaluate_hessians,
     const unsigned int active_fe_index,
     const unsigned int first_selected_component,
     const std::array<unsigned int, VectorizedArrayType::size()> cells,
@@ -258,9 +270,11 @@ namespace internal
       dof_info,
       values_quad,
       gradients_quad,
+      hessians_quad,
       scratch_data,
       evaluate_values,
       evaluate_gradients,
+      evaluate_hessians,
       active_fe_index,
       first_selected_component,
       cells,
@@ -285,9 +299,11 @@ namespace internal
     VectorizedArrayType *                                      values_array,
     VectorizedArrayType *                                      values_quad,
     VectorizedArrayType *                                      gradients_quad,
+    VectorizedArrayType *                                      hessians_quad,
     VectorizedArrayType *                                      scratch_data,
     const bool                                                 integrate_values,
     const bool         integrate_gradients,
+    const bool         integrate_hessians,
     const unsigned int active_fe_index,
     const unsigned int first_selected_component,
     const std::array<unsigned int, VectorizedArrayType::size()> cells,
@@ -314,9 +330,11 @@ namespace internal
       values_array,
       values_quad,
       gradients_quad,
+      hessians_quad,
       scratch_data,
       integrate_values,
       integrate_gradients,
+      integrate_hessians,
       active_fe_index,
       first_selected_component,
       cells,
index c1eb1afd98fa577c7d83a7409e9ec76807de48ab..66d7a750dc840a40356ce80a6bb9d1fac0aebd5b 100644 (file)
@@ -549,6 +549,8 @@ public:
   using value_type  = Tensor<1, n_components_, VectorizedArrayType>;
   using gradient_type =
     Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>>;
+  using hessian_type =
+    Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>;
   static constexpr unsigned int dimension    = dim;
   static constexpr unsigned int n_components = n_components_;
 
@@ -855,6 +857,21 @@ public:
   submit_normal_derivative(const value_type   grad_in,
                            const unsigned int q_point);
 
+  /**
+   * Write a contribution that is tested by the Hessian to the field
+   * containing the values at quadrature points with component @p q_point.
+   * Access to the same field as through get_hessian(). If applied before the
+   * function FEEvaluation::integrate(EvaluationFlags::hessians) is called,
+   * this specifies what is tested by the Hessians of all basis functions on the
+   * current cell and integrated over.
+   *
+   * Note that the derived class FEEvaluationAccess overloads this operation
+   * with specializations for the scalar case (n_components == 1) and for the
+   * vector-valued case (n_components == dim).
+   */
+  void
+  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
+
   /**
    * Return the Hessian of a finite element function at quadrature point
    * number @p q_point after a call to
@@ -1359,6 +1376,13 @@ protected:
    */
   bool gradients_quad_submitted;
 
+  /**
+   * Debug information to track whether hessians on quadrature points have
+   * been submitted for integration before the integration is actually stared.
+   * Used to control exceptions when uninitialized data is used.
+   */
+  bool hessians_quad_submitted;
+
   /**
    * For a FiniteElement with more than one base element, select at which
    * component this data structure should start.
@@ -1481,6 +1505,7 @@ public:
   using number_type                       = Number;
   using value_type                        = VectorizedArrayType;
   using gradient_type                     = Tensor<1, dim, VectorizedArrayType>;
+  using hessian_type                      = Tensor<2, dim, VectorizedArrayType>;
   static constexpr unsigned int dimension = dim;
   using BaseClass =
     FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>;
@@ -1544,7 +1569,7 @@ public:
   /**
    * @copydoc FEEvaluationBase<dim,1,Number,is_face>::get_hessian()
    */
-  Tensor<2, dim, VectorizedArrayType>
+  hessian_type
   get_hessian(unsigned int q_point) const;
 
   /**
@@ -1553,6 +1578,12 @@ public:
   gradient_type
   get_hessian_diagonal(const unsigned int q_point) const;
 
+  /**
+   * @copydoc FEEvaluationBase<dim,1,Number,is_face>::submit_hessian()
+   */
+  void
+  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
+
   /**
    * @copydoc FEEvaluationBase<dim,1,Number,is_face>::get_laplacian()
    */
@@ -1797,6 +1828,7 @@ public:
   using number_type                       = Number;
   using value_type                        = VectorizedArrayType;
   using gradient_type                     = Tensor<1, 1, VectorizedArrayType>;
+  using hessian_type                      = Tensor<2, 1, VectorizedArrayType>;
   static constexpr unsigned int dimension = 1;
   using BaseClass =
     FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>;
@@ -1885,7 +1917,7 @@ public:
   /**
    * @copydoc FEEvaluationBase<1,1,Number,is_face>::get_hessian()
    */
-  Tensor<2, 1, VectorizedArrayType>
+  hessian_type
   get_hessian(unsigned int q_point) const;
 
   /**
@@ -1894,6 +1926,12 @@ public:
   gradient_type
   get_hessian_diagonal(const unsigned int q_point) const;
 
+  /**
+   * @copydoc FEEvaluationBase<1,1,Number,is_face>::submit_hessian()
+   */
+  void
+  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
+
   /**
    * @copydoc FEEvaluationBase<1,1,Number,is_face>::get_laplacian()
    */
@@ -5774,7 +5812,6 @@ inline Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>
 FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   get_hessian(const unsigned int q_point) const
 {
-  Assert(!is_face, ExcNotImplemented());
 #  ifdef DEBUG
   Assert(this->hessians_quad_initialized == true,
          internal::ExcAccessToUninitializedField());
@@ -5795,7 +5832,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
 
   // Cartesian cell
-  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
+  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
     {
       for (unsigned int comp = 0; comp < n_components; comp++)
         {
@@ -5832,7 +5869,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
         }
     }
   // cell with general Jacobian, but constant within the cell
-  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
+  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
     {
       for (unsigned int comp = 0; comp < n_components; comp++)
         {
@@ -6195,6 +6232,110 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
 
 
+template <int dim,
+          int n_components_,
+          typename Number,
+          bool is_face,
+          typename VectorizedArrayType>
+inline DEAL_II_ALWAYS_INLINE void
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+  submit_hessian(
+    const Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>
+                       hessian_in,
+    const unsigned int q_point)
+{
+  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  AssertIndexRange(q_point, this->n_quadrature_points);
+  Assert(this->J_value != nullptr,
+         internal::ExcMatrixFreeAccessToUninitializedMappingField(
+           "update_hessians"));
+  Assert(this->jacobian != nullptr,
+         internal::ExcMatrixFreeAccessToUninitializedMappingField(
+           "update_hessians"));
+#  ifdef DEBUG
+  this->hessians_quad_submitted = true;
+#  endif
+
+  // compute hessian_unit = J^T * hessian_in(u) * J
+  const std::size_t      nqp  = this->n_quadrature_points;
+  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
+  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
+    {
+      const VectorizedArrayType JxW =
+        this->J_value[0] * this->quadrature_weights[q_point];
+
+      // diagonal part
+      for (unsigned int d = 0; d < dim; ++d)
+        {
+          const auto                jac_d  = this->jacobian[0][d][d];
+          const VectorizedArrayType factor = jac_d * jac_d * JxW;
+          for (unsigned int comp = 0; comp < n_components; comp++)
+            hessians_quad[(comp * hdim + d) * nqp + q_point] =
+              hessian_in[comp][d][d] * factor;
+        }
+
+      // off diagonal part (use symmetry)
+      for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
+        for (unsigned int e = 0; e < d; ++e, ++off_dia)
+          {
+            const auto                jac_d  = this->jacobian[0][d][d];
+            const auto                jac_e  = this->jacobian[0][e][e];
+            const VectorizedArrayType factor = jac_d * jac_e * JxW;
+            for (unsigned int comp = 0; comp < n_components; comp++)
+              hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
+                2.0 * hessian_in[comp][d][e] * factor;
+          }
+    }
+  // cell with general Jacobian, but constant within the cell
+  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
+    {
+      const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
+      const VectorizedArrayType                 JxW =
+        this->J_value[0] * this->quadrature_weights[q_point];
+      for (unsigned int comp = 0; comp < n_components; ++comp)
+        {
+          // 1. tmp = hessian_in(u) * J
+          VectorizedArrayType tmp[dim][dim];
+          for (unsigned int i = 0; i < dim; ++i)
+            for (unsigned int j = 0; j < dim; ++j)
+              tmp[i][j] = 0.;
+
+          for (unsigned int i = 0; i < dim; ++i)
+            for (unsigned int j = 0; j < dim; ++j)
+              for (unsigned int k = 0; k < dim; ++k)
+                tmp[i][j] += hessian_in[comp][i][k] * jac[k][j] * JxW;
+
+          // 2. hessian_unit = J^T * tmp
+
+          // diagonal part
+          for (unsigned int d = 0; d < dim; ++d)
+            {
+              hessians_quad[(comp * hdim + d) * nqp + q_point] = 0;
+              for (unsigned int i = 0; i < dim; ++i)
+                hessians_quad[(comp * hdim + d) * nqp + q_point] +=
+                  jac[i][d] * tmp[i][d];
+            }
+
+          // off diagonal part (use symmetry)
+          for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
+            for (unsigned int e = 0; e < d; ++e, ++off_dia)
+              {
+                hessians_quad[(comp * hdim + off_dia) * nqp + q_point] = 0;
+                for (unsigned int i = 0; i < dim; ++i)
+                  hessians_quad[(comp * hdim + off_dia) * nqp + q_point] +=
+                    2.0 * jac[i][e] * tmp[i][d];
+              }
+        }
+    }
+  else
+    {
+      // TODO: implement general path submit_hessian
+      AssertThrow(false, ExcNotImplemented());
+    }
+}
+
+
+
 template <int dim,
           int n_components_,
           typename Number,
@@ -6640,6 +6781,19 @@ FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
 
 
 
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline DEAL_II_ALWAYS_INLINE void
+FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
+  submit_hessian(const Tensor<2, dim, VectorizedArrayType> hessian_in,
+                 const unsigned int                        q_point)
+{
+  Tensor<1, 1, Tensor<2, dim, VectorizedArrayType>> hessian;
+  hessian[0] = hessian_in;
+  BaseClass::submit_hessian(hessian, q_point);
+}
+
+
+
 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
 inline VectorizedArrayType
 FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
@@ -7359,6 +7513,17 @@ FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
 }
 
 
+template <typename Number, bool is_face, typename VectorizedArrayType>
+inline DEAL_II_ALWAYS_INLINE void
+FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::submit_hessian(
+  const Tensor<2, 1, VectorizedArrayType> hessian_in,
+  const unsigned int                      q_point)
+{
+  Tensor<1, 1, Tensor<2, 1, VectorizedArrayType>> hessian;
+  hessian[0] = hessian_in;
+  BaseClass::submit_hessian(hessian, q_point);
+}
+
 
 template <typename Number, bool is_face, typename VectorizedArrayType>
 inline VectorizedArrayType
@@ -8173,6 +8338,7 @@ namespace internal
     const internal::MatrixFreeFunctions::DoFInfo *dof_info,
     VectorizedArrayType *                         values_quad,
     VectorizedArrayType *                         gradients_quad,
+    VectorizedArrayType *                         hessians_quad,
     VectorizedArrayType *                         scratch_data,
     const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> *data,
     const EvaluationFlags::EvaluationFlags integration_flag)
@@ -8212,6 +8378,7 @@ namespace internal
                       vec_values,
                       values_quad,
                       gradients_quad,
+                      hessians_quad,
                       scratch_data,
                       true);
         else
@@ -8222,6 +8389,7 @@ namespace internal
             vec_values,
             values_quad,
             gradients_quad,
+            hessians_quad,
             scratch_data,
             true);
 
@@ -8256,6 +8424,7 @@ namespace internal
     const VectorizedArrayType *,
     const VectorizedArrayType *,
     const VectorizedArrayType *,
+    const VectorizedArrayType *,
     const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> *,
     const EvaluationFlags::EvaluationFlags)
   {
@@ -8400,10 +8569,10 @@ FEEvaluation<dim,
          ExcNotInitialized());
 
   Assert(
-    (integration_flag &
-     ~(EvaluationFlags::values | EvaluationFlags::gradients)) == 0,
-    ExcMessage(
-      "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
+    (integration_flag & ~(EvaluationFlags::values | EvaluationFlags::gradients |
+                          EvaluationFlags::hessians)) == 0,
+    ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, and "
+               "EvaluationFlags::hessians are supported."));
 
   if (fe_degree > -1)
     SelectEvaluator<dim, fe_degree, n_q_points_1d, VectorizedArrayType>::
@@ -8413,6 +8582,7 @@ FEEvaluation<dim,
                 values_array,
                 this->values_quad,
                 this->gradients_quad,
+                this->hessians_quad,
                 this->scratch_data,
                 false);
   else
@@ -8423,6 +8593,7 @@ FEEvaluation<dim,
       values_array,
       this->values_quad,
       this->gradients_quad,
+      this->hessians_quad,
       this->scratch_data,
       false);
 
@@ -8491,6 +8662,7 @@ FEEvaluation<dim,
         this->dof_info,
         this->values_quad,
         this->gradients_quad,
+        this->hessians_quad,
         this->scratch_data,
         this->data,
         integration_flag) == false)
@@ -8815,13 +8987,14 @@ FEFaceEvaluation<dim,
            const EvaluationFlags::EvaluationFlags evaluation_flag)
 {
   Assert(
-    (evaluation_flag &
-     ~(EvaluationFlags::values | EvaluationFlags::gradients)) == 0,
+    (evaluation_flag & ~(EvaluationFlags::values | EvaluationFlags::gradients |
+                         EvaluationFlags::hessians)) == 0,
     ExcMessage(
-      "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
+      "Only EvaluationFlags::values, EvaluationFlags::gradients, and EvaluationFlags::hessians are supported."));
 
   if (!(evaluation_flag & EvaluationFlags::values) &&
-      !(evaluation_flag & EvaluationFlags::gradients))
+      !(evaluation_flag & EvaluationFlags::gradients) &&
+      !(evaluation_flag & EvaluationFlags::hessians))
     return;
 
   if (this->dof_access_index ==
@@ -8850,9 +9023,11 @@ FEFaceEvaluation<dim,
           values_array,
           this->begin_values(),
           this->begin_gradients(),
+          this->begin_hessians(),
           this->scratch_data,
           evaluation_flag & EvaluationFlags::values,
           evaluation_flag & EvaluationFlags::gradients,
+          evaluation_flag & EvaluationFlags::hessians,
           face_nos[0],
           this->subface_index,
           face_orientations[0],
@@ -8869,9 +9044,11 @@ FEFaceEvaluation<dim,
             values_array,
             this->begin_values(),
             this->begin_gradients(),
+            this->begin_hessians(),
             this->scratch_data,
             evaluation_flag & EvaluationFlags::values,
             evaluation_flag & EvaluationFlags::gradients,
+            evaluation_flag & EvaluationFlags::hessians,
             this->face_no,
             this->subface_index,
             this->face_orientation,
@@ -8883,9 +9060,11 @@ FEFaceEvaluation<dim,
                    values_array,
                    this->begin_values(),
                    this->begin_gradients(),
+                   this->begin_hessians(),
                    this->scratch_data,
                    evaluation_flag & EvaluationFlags::values,
                    evaluation_flag & EvaluationFlags::gradients,
+                   evaluation_flag & EvaluationFlags::hessians,
                    this->face_no,
                    this->subface_index,
                    this->face_orientation,
@@ -8897,6 +9076,8 @@ FEFaceEvaluation<dim,
     this->values_quad_initialized = true;
   if (evaluation_flag & EvaluationFlags::gradients)
     this->gradients_quad_initialized = true;
+  if (evaluation_flag & EvaluationFlags::hessians)
+    this->hessians_quad_initialized = true;
 #  endif
 }
 
@@ -8994,13 +9175,14 @@ FEFaceEvaluation<dim,
             VectorizedArrayType *                  values_array)
 {
   Assert(
-    (evaluation_flag &
-     ~(EvaluationFlags::values | EvaluationFlags::gradients)) == 0,
+    (evaluation_flag & ~(EvaluationFlags::values | EvaluationFlags::gradients |
+                         EvaluationFlags::hessians)) == 0,
     ExcMessage(
-      "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
+      "Only EvaluationFlags::values, EvaluationFlags::gradients, and EvaluationFlags::hessians are supported."));
 
   if (!(evaluation_flag & EvaluationFlags::values) &&
-      !(evaluation_flag & EvaluationFlags::gradients))
+      !(evaluation_flag & EvaluationFlags::gradients) &&
+      !(evaluation_flag & EvaluationFlags::hessians))
     return;
 
   if (fe_degree > -1)
@@ -9011,9 +9193,11 @@ FEFaceEvaluation<dim,
         values_array,
         this->begin_values(),
         this->begin_gradients(),
+        this->begin_hessians(),
         this->scratch_data,
         evaluation_flag & EvaluationFlags::values,
         evaluation_flag & EvaluationFlags::gradients,
+        evaluation_flag & EvaluationFlags::hessians,
         this->face_no,
         this->subface_index,
         this->face_orientation,
@@ -9025,9 +9209,11 @@ FEFaceEvaluation<dim,
                 values_array,
                 this->begin_values(),
                 this->begin_gradients(),
+                this->begin_hessians(),
                 this->scratch_data,
                 evaluation_flag & EvaluationFlags::values,
                 evaluation_flag & EvaluationFlags::gradients,
+                evaluation_flag & EvaluationFlags::hessians,
                 this->face_no,
                 this->subface_index,
                 this->face_orientation,
@@ -9082,10 +9268,10 @@ FEFaceEvaluation<dim,
                   const EvaluationFlags::EvaluationFlags evaluation_flag)
 {
   Assert(
-    (evaluation_flag &
-     ~(EvaluationFlags::values | EvaluationFlags::gradients)) == 0,
+    (evaluation_flag & ~(EvaluationFlags::values | EvaluationFlags::gradients |
+                         EvaluationFlags::hessians)) == 0,
     ExcMessage(
-      "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
+      "Only EvaluationFlags::values, EvaluationFlags::gradients, and EvaluationFlags::hessians are supported."));
 
   const auto fu = [&]() {
     const auto shared_vector_data = internal::get_shared_vector_data(
@@ -9116,9 +9302,11 @@ FEFaceEvaluation<dim,
           *this->dof_info,
           this->begin_values(),
           this->begin_gradients(),
+          this->begin_hessians(),
           this->scratch_data,
           evaluation_flag & EvaluationFlags::values,
           evaluation_flag & EvaluationFlags::gradients,
+          evaluation_flag & EvaluationFlags::hessians,
           this->active_fe_index,
           this->first_selected_component,
           cells,
@@ -9156,9 +9344,11 @@ FEFaceEvaluation<dim,
               *this->dof_info,
               this->begin_values(),
               this->begin_gradients(),
+              this->begin_hessians(),
               this->scratch_data,
               evaluation_flag & EvaluationFlags::values,
               evaluation_flag & EvaluationFlags::gradients,
+              evaluation_flag & EvaluationFlags::hessians,
               this->active_fe_index,
               this->first_selected_component,
               cells_,
@@ -9179,9 +9369,11 @@ FEFaceEvaluation<dim,
                               *this->dof_info,
                               this->begin_values(),
                               this->begin_gradients(),
+                              this->begin_hessians(),
                               this->scratch_data,
                               evaluation_flag & EvaluationFlags::values,
                               evaluation_flag & EvaluationFlags::gradients,
+                              evaluation_flag & EvaluationFlags::hessians,
                               this->active_fe_index,
                               this->first_selected_component,
                               cells_,
@@ -9204,6 +9396,8 @@ FEFaceEvaluation<dim,
     this->values_quad_initialized = true;
   if (evaluation_flag & EvaluationFlags::gradients)
     this->gradients_quad_initialized = true;
+  if (evaluation_flag & EvaluationFlags::hessians)
+    this->hessians_quad_initialized = true;
 #  endif
 }
 
@@ -9292,9 +9486,11 @@ FEFaceEvaluation<dim,
             this->begin_dof_values(),
             this->begin_values(),
             this->begin_gradients(),
+            this->begin_hessians(),
             this->scratch_data,
             evaluation_flag & EvaluationFlags::values,
             evaluation_flag & EvaluationFlags::gradients,
+            evaluation_flag & EvaluationFlags::hessians,
             this->active_fe_index,
             this->first_selected_component,
             cells_,
@@ -9323,9 +9519,11 @@ FEFaceEvaluation<dim,
                               this->begin_dof_values(),
                               this->begin_values(),
                               this->begin_gradients(),
+                              this->begin_hessians(),
                               this->scratch_data,
                               evaluation_flag & EvaluationFlags::values,
                               evaluation_flag & EvaluationFlags::gradients,
+                              evaluation_flag & EvaluationFlags::hessians,
                               this->active_fe_index,
                               this->first_selected_component,
                               cells_,
index 474d02b171d3a47284ff5b4491aa9d81ee2dfdd1..2e3e1ecf174554de06c0afbe33fd4c2d7dbc6ec9 100644 (file)
@@ -2292,10 +2292,14 @@ namespace internal
         const unsigned int n_q_points = my_data.descriptor[0].n_q_points;
         const unsigned int n_mapping_points =
           shape_info.dofs_per_component_on_cell;
+        constexpr unsigned int hess_dim = dim * (dim + 1) / 2;
 
         AlignedVector<VectorizedDouble> cell_points(dim * n_mapping_points);
         AlignedVector<VectorizedDouble> face_quads(dim * n_q_points);
         AlignedVector<VectorizedDouble> face_grads(dim * dim * n_q_points);
+        AlignedVector<VectorizedDouble> face_grad_grads(dim * hess_dim *
+                                                        n_q_points);
+
         AlignedVector<VectorizedDouble> scratch_data(
           dim * (2 * n_q_points + 3 * n_mapping_points));
 
@@ -2326,9 +2330,11 @@ namespace internal
                 cell_points.data(),
                 face_quads.data(),
                 face_grads.data(),
+                face_grad_grads.data(),
                 scratch_data.data(),
                 true,
                 true,
+                update_flags_faces & update_jacobian_grads,
                 face_no,
                 GeometryInfo<dim>::max_children_per_cell,
                 faces[face].face_orientation > 8 ?
@@ -2377,6 +2383,12 @@ namespace internal
                           my_data.jacobians[0][offset + q][d][e]);
                     }
 
+                  if (update_flags_faces & update_jacobian_grads)
+                    {
+                      // TODO: implement jacobian grads for general path
+                      AssertThrow(false, ExcNotImplemented());
+                    }
+
                   std::array<Tensor<1, dim, VectorizedDouble>, dim - 1>
                     tangential_vectors;
                   for (unsigned int d = 0; d != dim - 1; ++d)
@@ -2443,9 +2455,11 @@ namespace internal
                              cell_points.data(),
                              face_quads.data(),
                              face_grads.data(),
+                             face_grad_grads.data(),
                              scratch_data.data(),
                              false,
                              true,
+                             update_flags_faces & update_jacobian_grads,
                              faces[face].exterior_face_no,
                              faces[face].subface_index,
                              faces[face].face_orientation < 8 ?
@@ -2477,6 +2491,13 @@ namespace internal
                               vv,
                               my_data.jacobians[1][offset + q][d][e]);
                         }
+
+                      if (update_flags_faces & update_jacobian_grads)
+                        {
+                          // TODO: implement jacobian grads for general path
+                          AssertThrow(false, ExcNotImplemented());
+                        }
+
                       my_data.normals_times_jacobians[1][offset + q] =
                         my_data.normal_vectors[offset + q] *
                         my_data.jacobians[1][offset + q];
index 2785518be34360910258c3bbd46c9003db86ca9c..8a50e884c41bb30a84ea7a86bbc57dc196f46215 100644 (file)
@@ -5,12 +5,15 @@ DEAL:1d::Error cell gradients:   0
 DEAL:1d::Error cell Hessians:    0
 DEAL:1d::Error face- values:     0
 DEAL:1d::Error face- gradients:  0
+DEAL:1d::Error face- Hessians:   0
 DEAL:1d::Error face- grad*normal:0
 DEAL:1d::Error face+ values:     0
 DEAL:1d::Error face+ gradients:  0
+DEAL:1d::Error face+ Hessians:   0
 DEAL:1d::Error face+ grad*normal:0
 DEAL:1d::Error face b values:    0
 DEAL:1d::Error face b gradients: 0
+DEAL:1d::Error face b Hessians:  0
 DEAL:1d::Error face b grad*norm: 0
 DEAL:1d::
 DEAL:1d::Test without templates on FEEvaluation
@@ -20,12 +23,15 @@ DEAL:1d::Error cell gradients:   0
 DEAL:1d::Error cell Hessians:    0
 DEAL:1d::Error face- values:     0
 DEAL:1d::Error face- gradients:  0
+DEAL:1d::Error face- Hessians:   0
 DEAL:1d::Error face- grad*normal:0
 DEAL:1d::Error face+ values:     0
 DEAL:1d::Error face+ gradients:  0
+DEAL:1d::Error face+ Hessians:   0
 DEAL:1d::Error face+ grad*normal:0
 DEAL:1d::Error face b values:    0
 DEAL:1d::Error face b gradients: 0
+DEAL:1d::Error face b Hessians:  0
 DEAL:1d::Error face b grad*norm: 0
 DEAL:1d::
 DEAL:1d::Testing FE_Q<1>(2)
@@ -34,12 +40,15 @@ DEAL:1d::Error cell gradients:   0
 DEAL:1d::Error cell Hessians:    0
 DEAL:1d::Error face- values:     0
 DEAL:1d::Error face- gradients:  0
+DEAL:1d::Error face- Hessians:   0
 DEAL:1d::Error face- grad*normal:0
 DEAL:1d::Error face+ values:     0
 DEAL:1d::Error face+ gradients:  0
+DEAL:1d::Error face+ Hessians:   0
 DEAL:1d::Error face+ grad*normal:0
 DEAL:1d::Error face b values:    0
 DEAL:1d::Error face b gradients: 0
+DEAL:1d::Error face b Hessians:  0
 DEAL:1d::Error face b grad*norm: 0
 DEAL:1d::
 DEAL:1d::Testing FE_DGQ<1>(2)
@@ -48,12 +57,15 @@ DEAL:1d::Error cell gradients:   0
 DEAL:1d::Error cell Hessians:    0
 DEAL:1d::Error face- values:     0
 DEAL:1d::Error face- gradients:  0
+DEAL:1d::Error face- Hessians:   0
 DEAL:1d::Error face- grad*normal:0
 DEAL:1d::Error face+ values:     0
 DEAL:1d::Error face+ gradients:  0
+DEAL:1d::Error face+ Hessians:   0
 DEAL:1d::Error face+ grad*normal:0
 DEAL:1d::Error face b values:    0
 DEAL:1d::Error face b gradients: 0
+DEAL:1d::Error face b Hessians:  0
 DEAL:1d::Error face b grad*norm: 0
 DEAL:1d::
 DEAL:1d::Test without templates on FEEvaluation
@@ -63,12 +75,15 @@ DEAL:1d::Error cell gradients:   0
 DEAL:1d::Error cell Hessians:    0
 DEAL:1d::Error face- values:     0
 DEAL:1d::Error face- gradients:  0
+DEAL:1d::Error face- Hessians:   0
 DEAL:1d::Error face- grad*normal:0
 DEAL:1d::Error face+ values:     0
 DEAL:1d::Error face+ gradients:  0
+DEAL:1d::Error face+ Hessians:   0
 DEAL:1d::Error face+ grad*normal:0
 DEAL:1d::Error face b values:    0
 DEAL:1d::Error face b gradients: 0
+DEAL:1d::Error face b Hessians:  0
 DEAL:1d::Error face b grad*norm: 0
 DEAL:1d::
 DEAL:2d::Testing FE_DGQ<2>(0)
@@ -77,12 +92,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Test without templates on FEEvaluation
@@ -92,12 +110,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_Q<2>(1)
@@ -106,12 +127,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_DGQ<2>(1)
@@ -120,12 +144,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Test without templates on FEEvaluation
@@ -135,12 +162,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_Q<2>(2)
@@ -149,12 +179,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_DGQ<2>(2)
@@ -163,12 +196,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Test without templates on FEEvaluation
@@ -178,12 +214,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_Q<2>(3)
@@ -192,12 +231,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_DGQ<2>(3)
@@ -206,12 +248,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Test without templates on FEEvaluation
@@ -221,12 +266,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_Q<2>(4)
@@ -235,12 +283,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_DGQ<2>(4)
@@ -249,12 +300,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Test without templates on FEEvaluation
@@ -264,12 +318,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:3d::Testing FE_Q<3>(1)
@@ -278,12 +335,15 @@ DEAL:3d::Error cell gradients:   0
 DEAL:3d::Error cell Hessians:    0
 DEAL:3d::Error face- values:     0
 DEAL:3d::Error face- gradients:  0
+DEAL:3d::Error face- Hessians:   0
 DEAL:3d::Error face- grad*normal:0
 DEAL:3d::Error face+ values:     0
 DEAL:3d::Error face+ gradients:  0
+DEAL:3d::Error face+ Hessians:   0
 DEAL:3d::Error face+ grad*normal:0
 DEAL:3d::Error face b values:    0
 DEAL:3d::Error face b gradients: 0
+DEAL:3d::Error face b Hessians:  0
 DEAL:3d::Error face b grad*norm: 0
 DEAL:3d::
 DEAL:3d::Testing FE_DGQ<3>(1)
@@ -292,12 +352,15 @@ DEAL:3d::Error cell gradients:   0
 DEAL:3d::Error cell Hessians:    0
 DEAL:3d::Error face- values:     0
 DEAL:3d::Error face- gradients:  0
+DEAL:3d::Error face- Hessians:   0
 DEAL:3d::Error face- grad*normal:0
 DEAL:3d::Error face+ values:     0
 DEAL:3d::Error face+ gradients:  0
+DEAL:3d::Error face+ Hessians:   0
 DEAL:3d::Error face+ grad*normal:0
 DEAL:3d::Error face b values:    0
 DEAL:3d::Error face b gradients: 0
+DEAL:3d::Error face b Hessians:  0
 DEAL:3d::Error face b grad*norm: 0
 DEAL:3d::
 DEAL:3d::Test without templates on FEEvaluation
@@ -307,12 +370,15 @@ DEAL:3d::Error cell gradients:   0
 DEAL:3d::Error cell Hessians:    0
 DEAL:3d::Error face- values:     0
 DEAL:3d::Error face- gradients:  0
+DEAL:3d::Error face- Hessians:   0
 DEAL:3d::Error face- grad*normal:0
 DEAL:3d::Error face+ values:     0
 DEAL:3d::Error face+ gradients:  0
+DEAL:3d::Error face+ Hessians:   0
 DEAL:3d::Error face+ grad*normal:0
 DEAL:3d::Error face b values:    0
 DEAL:3d::Error face b gradients: 0
+DEAL:3d::Error face b Hessians:  0
 DEAL:3d::Error face b grad*norm: 0
 DEAL:3d::
 DEAL:3d::Testing FE_Q<3>(2)
@@ -321,12 +387,15 @@ DEAL:3d::Error cell gradients:   0
 DEAL:3d::Error cell Hessians:    0
 DEAL:3d::Error face- values:     0
 DEAL:3d::Error face- gradients:  0
+DEAL:3d::Error face- Hessians:   0
 DEAL:3d::Error face- grad*normal:0
 DEAL:3d::Error face+ values:     0
 DEAL:3d::Error face+ gradients:  0
+DEAL:3d::Error face+ Hessians:   0
 DEAL:3d::Error face+ grad*normal:0
 DEAL:3d::Error face b values:    0
 DEAL:3d::Error face b gradients: 0
+DEAL:3d::Error face b Hessians:  0
 DEAL:3d::Error face b grad*norm: 0
 DEAL:3d::
 DEAL:3d::Testing FE_DGQ<3>(2)
@@ -335,12 +404,15 @@ DEAL:3d::Error cell gradients:   0
 DEAL:3d::Error cell Hessians:    0
 DEAL:3d::Error face- values:     0
 DEAL:3d::Error face- gradients:  0
+DEAL:3d::Error face- Hessians:   0
 DEAL:3d::Error face- grad*normal:0
 DEAL:3d::Error face+ values:     0
 DEAL:3d::Error face+ gradients:  0
+DEAL:3d::Error face+ Hessians:   0
 DEAL:3d::Error face+ grad*normal:0
 DEAL:3d::Error face b values:    0
 DEAL:3d::Error face b gradients: 0
+DEAL:3d::Error face b Hessians:  0
 DEAL:3d::Error face b grad*norm: 0
 DEAL:3d::
 DEAL:3d::Test without templates on FEEvaluation
@@ -350,11 +422,14 @@ DEAL:3d::Error cell gradients:   0
 DEAL:3d::Error cell Hessians:    0
 DEAL:3d::Error face- values:     0
 DEAL:3d::Error face- gradients:  0
+DEAL:3d::Error face- Hessians:   0
 DEAL:3d::Error face- grad*normal:0
 DEAL:3d::Error face+ values:     0
 DEAL:3d::Error face+ gradients:  0
+DEAL:3d::Error face+ Hessians:   0
 DEAL:3d::Error face+ grad*normal:0
 DEAL:3d::Error face b values:    0
 DEAL:3d::Error face b gradients: 0
+DEAL:3d::Error face b Hessians:  0
 DEAL:3d::Error face b grad*norm: 0
 DEAL:3d::
index b2cb4c4a1479343bdc46d676bf11ab6c332641a7..61a6acd6fd3c5e502e7e8a954603df17244de25d 100644 (file)
@@ -5,12 +5,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_DGQ<2>(3)
@@ -19,12 +22,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Test without templates on FEEvaluation
@@ -34,12 +40,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_Q<2>(4)
@@ -48,12 +57,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_DGQ<2>(4)
@@ -62,12 +74,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Test without templates on FEEvaluation
@@ -77,11 +92,14 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
index a5ec16e84236dfaef96c5e09ed17b037af488d0a..dbbf1588f7ed12d887d266566e9d5a711a3864fa 100644 (file)
@@ -126,10 +126,12 @@ public:
       {
         fe_evalm.reinit(face);
         fe_evalm.read_dof_values(src);
-        fe_evalm.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+        fe_evalm.evaluate(EvaluationFlags::values | EvaluationFlags::gradients |
+                          EvaluationFlags::hessians);
         fe_evalp.reinit(face);
         fe_evalp.read_dof_values(src);
-        fe_evalp.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+        fe_evalp.evaluate(EvaluationFlags::values | EvaluationFlags::gradients |
+                          EvaluationFlags::hessians);
 
         for (unsigned int j = 0; j < VectorizedArray<Number>::size(); ++j)
           {
@@ -142,37 +144,41 @@ public:
                 ++facem_times;
                 ++facep_times;
                 Point<dim> p;
+
+                // interior face
                 for (unsigned int d = 0; d < dim; ++d)
                   p[d] = fe_evalm.quadrature_point(q)[d][j];
                 facem_errors[0] +=
                   std::abs(fe_evalm.get_value(q)[j] - function.value(p, 0));
                 for (unsigned int d = 0; d < dim; ++d)
-                  {
-                    // std::cout << fe_evalm.get_gradient(q)[d][j] << " ";
-                    facem_errors[1] += std::abs(fe_evalm.get_gradient(q)[d][j] -
-                                                function.gradient(p, 0)[d]);
-                  }
+                  facem_errors[1] += std::abs(fe_evalm.get_gradient(q)[d][j] -
+                                              function.gradient(p, 0)[d]);
+                for (unsigned int d = 0; d < dim; ++d)
+                  for (unsigned int e = 0; e < dim; ++e)
+                    facem_errors[2] +=
+                      std::abs(fe_evalm.get_hessian(q)[d][e][j] -
+                               function.hessian(p, 0)[d][e]);
                 double normal_derivative = 0;
                 for (unsigned int d = 0; d < dim; ++d)
                   normal_derivative += function.gradient(p, 0)[d] *
                                        fe_evalm.get_normal_vector(q)[d][j];
                 facem_errors[3] += std::abs(
                   fe_evalm.get_normal_derivative(q)[j] - normal_derivative);
+
+                // exterior face
                 facep_errors[0] +=
                   std::abs(fe_evalp.get_value(q)[j] - function.value(p, 0));
                 for (unsigned int d = 0; d < dim; ++d)
                   facep_errors[1] += std::abs(fe_evalp.get_gradient(q)[d][j] -
                                               function.gradient(p, 0)[d]);
+                for (unsigned int d = 0; d < dim; ++d)
+                  for (unsigned int e = 0; e < dim; ++e)
+                    facep_errors[2] +=
+                      std::abs(fe_evalp.get_hessian(q)[d][e][j] -
+                               function.hessian(p, 0)[d][e]);
                 facep_errors[3] += std::abs(
                   fe_evalp.get_normal_derivative(q)[j] - normal_derivative);
-                // hessians not currently implemented in FEFaceEvaluation
-                // for (unsigned int d=0; d<dim; ++d)
-                //  for (unsigned int e=0; e<dim; ++e)
-                //    facem_errors[2] +=
-                //    std::abs(fe_evalm.get_hessian(q)[d][e][j]-
-                //                                function.hessian(q)[d][e]);
               }
-            // std::cout << std::endl;
           }
       }
   }
@@ -192,7 +198,8 @@ public:
       {
         fe_evalm.reinit(face);
         fe_evalm.read_dof_values(src);
-        fe_evalm.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+        fe_evalm.evaluate(EvaluationFlags::values | EvaluationFlags::gradients |
+                          EvaluationFlags::hessians);
 
         for (unsigned int j = 0; j < VectorizedArray<Number>::size(); ++j)
           {
@@ -212,18 +219,17 @@ public:
                   boundary_errors[1] +=
                     std::abs(fe_evalm.get_gradient(q)[d][j] -
                              function.gradient(p, 0)[d]);
+                for (unsigned int d = 0; d < dim; ++d)
+                  for (unsigned int e = 0; e < dim; ++e)
+                    boundary_errors[2] +=
+                      std::abs(fe_evalm.get_hessian(q)[d][e][j] -
+                               function.hessian(p, 0)[d][e]);
                 double normal_derivative = 0;
                 for (unsigned int d = 0; d < dim; ++d)
                   normal_derivative += function.gradient(p, 0)[d] *
                                        fe_evalm.get_normal_vector(q)[d][j];
                 boundary_errors[3] += std::abs(
                   fe_evalm.get_normal_derivative(q)[j] - normal_derivative);
-                // hessians not currently implemented in FEFaceEvaluation
-                // for (unsigned int d=0; d<dim; ++d)
-                //  for (unsigned int e=0; e<dim; ++e)
-                //    boundary_errors[2] +=
-                //    std::abs(fe_evalm.get_hessian(q)[d][e][j]-
-                //                                   function.hessian(q)[d][e]);
               }
           }
       }
@@ -273,24 +279,24 @@ public:
             << std::endl;
     deallog << "Error face- gradients:  " << facem_errors[1] / facem_times
             << std::endl;
-    // deallog << "Error face- Hessians:   " << facem_errors[2]/facem_times <<
-    // std::endl;
+    deallog << "Error face- Hessians:   " << facem_errors[2] / facem_times
+            << std::endl;
     deallog << "Error face- grad*normal:" << facem_errors[3] / facem_times
             << std::endl;
     deallog << "Error face+ values:     " << facep_errors[0] / facep_times
             << std::endl;
     deallog << "Error face+ gradients:  " << facep_errors[1] / facep_times
             << std::endl;
-    // deallog << "Error face+ Hessians:   " << facep_errors[2]/facep_times <<
-    // std::endl;
-    deallog << "Error face+ grad*normal:" << facem_errors[3] / facep_times
+    deallog << "Error face+ Hessians:   " << facep_errors[2] / facep_times
+            << std::endl;
+    deallog << "Error face+ grad*normal:" << facep_errors[3] / facep_times
             << std::endl;
     deallog << "Error face b values:    " << boundary_errors[0] / boundary_times
             << std::endl;
     deallog << "Error face b gradients: " << boundary_errors[1] / boundary_times
             << std::endl;
-    // deallog << "Error face b Hessians:  " <<
-    // boundary_errors[2]/boundary_times << std::endl;
+    deallog << "Error face b Hessians:  " << boundary_errors[2] / boundary_times
+            << std::endl;
     deallog << "Error face b grad*norm: " << boundary_errors[3] / boundary_times
             << std::endl;
     deallog << std::endl;
index f432deda3300b045017062a1d4e3d6be732b3d50..7d6076d705d2158d9c70da71cc848f12272571ce 100644 (file)
@@ -5,12 +5,15 @@ DEAL:1d::Error cell gradients:   0
 DEAL:1d::Error cell Hessians:    0
 DEAL:1d::Error face- values:     0
 DEAL:1d::Error face- gradients:  0
+DEAL:1d::Error face- Hessians:   0
 DEAL:1d::Error face- grad*normal:0
 DEAL:1d::Error face+ values:     0
 DEAL:1d::Error face+ gradients:  0
+DEAL:1d::Error face+ Hessians:   0
 DEAL:1d::Error face+ grad*normal:0
 DEAL:1d::Error face b values:    0
 DEAL:1d::Error face b gradients: 0
+DEAL:1d::Error face b Hessians:  0
 DEAL:1d::Error face b grad*norm: 0
 DEAL:1d::
 DEAL:1d::Testing FE_DGQ<1>(2)
@@ -19,12 +22,15 @@ DEAL:1d::Error cell gradients:   0
 DEAL:1d::Error cell Hessians:    0
 DEAL:1d::Error face- values:     0
 DEAL:1d::Error face- gradients:  0
+DEAL:1d::Error face- Hessians:   0
 DEAL:1d::Error face- grad*normal:0
 DEAL:1d::Error face+ values:     0
 DEAL:1d::Error face+ gradients:  0
+DEAL:1d::Error face+ Hessians:   0
 DEAL:1d::Error face+ grad*normal:0
 DEAL:1d::Error face b values:    0
 DEAL:1d::Error face b gradients: 0
+DEAL:1d::Error face b Hessians:  0
 DEAL:1d::Error face b grad*norm: 0
 DEAL:1d::
 DEAL:1d::Test without templates on FEEvaluation
@@ -34,12 +40,15 @@ DEAL:1d::Error cell gradients:   0
 DEAL:1d::Error cell Hessians:    0
 DEAL:1d::Error face- values:     0
 DEAL:1d::Error face- gradients:  0
+DEAL:1d::Error face- Hessians:   0
 DEAL:1d::Error face- grad*normal:0
 DEAL:1d::Error face+ values:     0
 DEAL:1d::Error face+ gradients:  0
+DEAL:1d::Error face+ Hessians:   0
 DEAL:1d::Error face+ grad*normal:0
 DEAL:1d::Error face b values:    0
 DEAL:1d::Error face b gradients: 0
+DEAL:1d::Error face b Hessians:  0
 DEAL:1d::Error face b grad*norm: 0
 DEAL:1d::
 DEAL:2d::Testing FE_Q<2>(1)
@@ -48,12 +57,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_DGQ<2>(1)
@@ -62,12 +74,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Test without templates on FEEvaluation
@@ -77,12 +92,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_Q<2>(2)
@@ -91,12 +109,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_DGQ<2>(2)
@@ -105,12 +126,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Test without templates on FEEvaluation
@@ -120,12 +144,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_Q<2>(3)
@@ -134,12 +161,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_DGQ<2>(3)
@@ -148,12 +178,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Test without templates on FEEvaluation
@@ -163,12 +196,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_Q<2>(4)
@@ -177,12 +213,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_DGQ<2>(4)
@@ -191,12 +230,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Test without templates on FEEvaluation
@@ -206,12 +248,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:3d::Testing FE_Q<3>(1)
@@ -220,12 +265,15 @@ DEAL:3d::Error cell gradients:   0
 DEAL:3d::Error cell Hessians:    0
 DEAL:3d::Error face- values:     0
 DEAL:3d::Error face- gradients:  0
+DEAL:3d::Error face- Hessians:   0
 DEAL:3d::Error face- grad*normal:0
 DEAL:3d::Error face+ values:     0
 DEAL:3d::Error face+ gradients:  0
+DEAL:3d::Error face+ Hessians:   0
 DEAL:3d::Error face+ grad*normal:0
 DEAL:3d::Error face b values:    0
 DEAL:3d::Error face b gradients: 0
+DEAL:3d::Error face b Hessians:  0
 DEAL:3d::Error face b grad*norm: 0
 DEAL:3d::
 DEAL:3d::Testing FE_DGQ<3>(1)
@@ -234,12 +282,15 @@ DEAL:3d::Error cell gradients:   0
 DEAL:3d::Error cell Hessians:    0
 DEAL:3d::Error face- values:     0
 DEAL:3d::Error face- gradients:  0
+DEAL:3d::Error face- Hessians:   0
 DEAL:3d::Error face- grad*normal:0
 DEAL:3d::Error face+ values:     0
 DEAL:3d::Error face+ gradients:  0
+DEAL:3d::Error face+ Hessians:   0
 DEAL:3d::Error face+ grad*normal:0
 DEAL:3d::Error face b values:    0
 DEAL:3d::Error face b gradients: 0
+DEAL:3d::Error face b Hessians:  0
 DEAL:3d::Error face b grad*norm: 0
 DEAL:3d::
 DEAL:3d::Test without templates on FEEvaluation
@@ -249,12 +300,15 @@ DEAL:3d::Error cell gradients:   0
 DEAL:3d::Error cell Hessians:    0
 DEAL:3d::Error face- values:     0
 DEAL:3d::Error face- gradients:  0
+DEAL:3d::Error face- Hessians:   0
 DEAL:3d::Error face- grad*normal:0
 DEAL:3d::Error face+ values:     0
 DEAL:3d::Error face+ gradients:  0
+DEAL:3d::Error face+ Hessians:   0
 DEAL:3d::Error face+ grad*normal:0
 DEAL:3d::Error face b values:    0
 DEAL:3d::Error face b gradients: 0
+DEAL:3d::Error face b Hessians:  0
 DEAL:3d::Error face b grad*norm: 0
 DEAL:3d::
 DEAL:3d::Testing FE_Q<3>(2)
@@ -263,12 +317,15 @@ DEAL:3d::Error cell gradients:   0
 DEAL:3d::Error cell Hessians:    0
 DEAL:3d::Error face- values:     0
 DEAL:3d::Error face- gradients:  0
+DEAL:3d::Error face- Hessians:   0
 DEAL:3d::Error face- grad*normal:0
 DEAL:3d::Error face+ values:     0
 DEAL:3d::Error face+ gradients:  0
+DEAL:3d::Error face+ Hessians:   0
 DEAL:3d::Error face+ grad*normal:0
 DEAL:3d::Error face b values:    0
 DEAL:3d::Error face b gradients: 0
+DEAL:3d::Error face b Hessians:  0
 DEAL:3d::Error face b grad*norm: 0
 DEAL:3d::
 DEAL:3d::Testing FE_DGQ<3>(2)
@@ -277,12 +334,15 @@ DEAL:3d::Error cell gradients:   0
 DEAL:3d::Error cell Hessians:    0
 DEAL:3d::Error face- values:     0
 DEAL:3d::Error face- gradients:  0
+DEAL:3d::Error face- Hessians:   0
 DEAL:3d::Error face- grad*normal:0
 DEAL:3d::Error face+ values:     0
 DEAL:3d::Error face+ gradients:  0
+DEAL:3d::Error face+ Hessians:   0
 DEAL:3d::Error face+ grad*normal:0
 DEAL:3d::Error face b values:    0
 DEAL:3d::Error face b gradients: 0
+DEAL:3d::Error face b Hessians:  0
 DEAL:3d::Error face b grad*norm: 0
 DEAL:3d::
 DEAL:3d::Test without templates on FEEvaluation
@@ -292,11 +352,14 @@ DEAL:3d::Error cell gradients:   0
 DEAL:3d::Error cell Hessians:    0
 DEAL:3d::Error face- values:     0
 DEAL:3d::Error face- gradients:  0
+DEAL:3d::Error face- Hessians:   0
 DEAL:3d::Error face- grad*normal:0
 DEAL:3d::Error face+ values:     0
 DEAL:3d::Error face+ gradients:  0
+DEAL:3d::Error face+ Hessians:   0
 DEAL:3d::Error face+ grad*normal:0
 DEAL:3d::Error face b values:    0
 DEAL:3d::Error face b gradients: 0
+DEAL:3d::Error face b Hessians:  0
 DEAL:3d::Error face b grad*norm: 0
 DEAL:3d::
index 8886b3804898f86286ce4d293fc9a66bc4924af5..fc2208e70abf7742ee00afa17d45a0cde5a70184 100644 (file)
@@ -5,12 +5,15 @@ DEAL:1d::Error cell gradients:   0
 DEAL:1d::Error cell Hessians:    0
 DEAL:1d::Error face- values:     0
 DEAL:1d::Error face- gradients:  0
+DEAL:1d::Error face- Hessians:   0
 DEAL:1d::Error face- grad*normal:0
 DEAL:1d::Error face+ values:     0
 DEAL:1d::Error face+ gradients:  0
+DEAL:1d::Error face+ Hessians:   0
 DEAL:1d::Error face+ grad*normal:0
 DEAL:1d::Error face b values:    0
 DEAL:1d::Error face b gradients: 0
+DEAL:1d::Error face b Hessians:  0
 DEAL:1d::Error face b grad*norm: 0
 DEAL:1d::
 DEAL:1d::Testing FE_DGQ<1>(2)
@@ -19,12 +22,15 @@ DEAL:1d::Error cell gradients:   0
 DEAL:1d::Error cell Hessians:    0
 DEAL:1d::Error face- values:     0
 DEAL:1d::Error face- gradients:  0
+DEAL:1d::Error face- Hessians:   0
 DEAL:1d::Error face- grad*normal:0
 DEAL:1d::Error face+ values:     0
 DEAL:1d::Error face+ gradients:  0
+DEAL:1d::Error face+ Hessians:   0
 DEAL:1d::Error face+ grad*normal:0
 DEAL:1d::Error face b values:    0
 DEAL:1d::Error face b gradients: 0
+DEAL:1d::Error face b Hessians:  0
 DEAL:1d::Error face b grad*norm: 0
 DEAL:1d::
 DEAL:1d::Test without templates on FEEvaluation
@@ -34,12 +40,15 @@ DEAL:1d::Error cell gradients:   0
 DEAL:1d::Error cell Hessians:    0
 DEAL:1d::Error face- values:     0
 DEAL:1d::Error face- gradients:  0
+DEAL:1d::Error face- Hessians:   0
 DEAL:1d::Error face- grad*normal:0
 DEAL:1d::Error face+ values:     0
 DEAL:1d::Error face+ gradients:  0
+DEAL:1d::Error face+ Hessians:   0
 DEAL:1d::Error face+ grad*normal:0
 DEAL:1d::Error face b values:    0
 DEAL:1d::Error face b gradients: 0
+DEAL:1d::Error face b Hessians:  0
 DEAL:1d::Error face b grad*norm: 0
 DEAL:1d::
 DEAL:2d::Testing FE_Q<2>(2)
@@ -48,12 +57,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_DGQ<2>(2)
@@ -62,12 +74,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Test without templates on FEEvaluation
@@ -77,12 +92,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_Q<2>(3)
@@ -91,12 +109,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_DGQ<2>(3)
@@ -105,12 +126,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Test without templates on FEEvaluation
@@ -120,12 +144,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_Q<2>(4)
@@ -134,12 +161,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Testing FE_DGQ<2>(4)
@@ -148,12 +178,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:2d::Test without templates on FEEvaluation
@@ -163,12 +196,15 @@ DEAL:2d::Error cell gradients:   0
 DEAL:2d::Error cell Hessians:    0
 DEAL:2d::Error face- values:     0
 DEAL:2d::Error face- gradients:  0
+DEAL:2d::Error face- Hessians:   0
 DEAL:2d::Error face- grad*normal:0
 DEAL:2d::Error face+ values:     0
 DEAL:2d::Error face+ gradients:  0
+DEAL:2d::Error face+ Hessians:   0
 DEAL:2d::Error face+ grad*normal:0
 DEAL:2d::Error face b values:    0
 DEAL:2d::Error face b gradients: 0
+DEAL:2d::Error face b Hessians:  0
 DEAL:2d::Error face b grad*norm: 0
 DEAL:2d::
 DEAL:3d::Testing FE_Q<3>(2)
@@ -177,12 +213,15 @@ DEAL:3d::Error cell gradients:   0
 DEAL:3d::Error cell Hessians:    0
 DEAL:3d::Error face- values:     0
 DEAL:3d::Error face- gradients:  0
+DEAL:3d::Error face- Hessians:   0
 DEAL:3d::Error face- grad*normal:0
 DEAL:3d::Error face+ values:     0
 DEAL:3d::Error face+ gradients:  0
+DEAL:3d::Error face+ Hessians:   0
 DEAL:3d::Error face+ grad*normal:0
 DEAL:3d::Error face b values:    0
 DEAL:3d::Error face b gradients: 0
+DEAL:3d::Error face b Hessians:  0
 DEAL:3d::Error face b grad*norm: 0
 DEAL:3d::
 DEAL:3d::Testing FE_DGQ<3>(2)
@@ -191,12 +230,15 @@ DEAL:3d::Error cell gradients:   0
 DEAL:3d::Error cell Hessians:    0
 DEAL:3d::Error face- values:     0
 DEAL:3d::Error face- gradients:  0
+DEAL:3d::Error face- Hessians:   0
 DEAL:3d::Error face- grad*normal:0
 DEAL:3d::Error face+ values:     0
 DEAL:3d::Error face+ gradients:  0
+DEAL:3d::Error face+ Hessians:   0
 DEAL:3d::Error face+ grad*normal:0
 DEAL:3d::Error face b values:    0
 DEAL:3d::Error face b gradients: 0
+DEAL:3d::Error face b Hessians:  0
 DEAL:3d::Error face b grad*norm: 0
 DEAL:3d::
 DEAL:3d::Test without templates on FEEvaluation
@@ -206,11 +248,14 @@ DEAL:3d::Error cell gradients:   0
 DEAL:3d::Error cell Hessians:    0
 DEAL:3d::Error face- values:     0
 DEAL:3d::Error face- gradients:  0
+DEAL:3d::Error face- Hessians:   0
 DEAL:3d::Error face- grad*normal:0
 DEAL:3d::Error face+ values:     0
 DEAL:3d::Error face+ gradients:  0
+DEAL:3d::Error face+ Hessians:   0
 DEAL:3d::Error face+ grad*normal:0
 DEAL:3d::Error face b values:    0
 DEAL:3d::Error face b gradients: 0
+DEAL:3d::Error face b Hessians:  0
 DEAL:3d::Error face b grad*norm: 0
 DEAL:3d::
index 4e0da38c5b590164d7e27a6328c6572154e90aa0..6706195df71f9f6b4f2d782623c3e2be075a45c0 100644 (file)
@@ -124,6 +124,7 @@ test(const unsigned int n_refinements = 1)
                   phi.begin_values(),
                   temp.data(),
                   false,
+                  false,
                   face);
 
 
diff --git a/tests/matrix_free/matrix_vector_hessians_faces.cc b/tests/matrix_free/matrix_vector_hessians_faces.cc
new file mode 100644 (file)
index 0000000..801271e
--- /dev/null
@@ -0,0 +1,396 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2013 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q_generic.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+// Tests if applying a matrix vector product with a matrix containing hessians
+// on faces produces the same result with FEFaceEvaluation and FEFaceValues.
+// This is checked for different combinations of EvaluationFlags, FE types and
+// polynomial degrees.
+
+template <int dim>
+void
+test_hessians(const unsigned int                             degree,
+              const dealii::FE_Poly<dim> &                   fe,
+              const dealii::EvaluationFlags::EvaluationFlags evaluation_flags)
+{
+  using namespace dealii;
+  using VectorizedArrayType = VectorizedArray<double>;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(3);
+
+  // TODO: make Hessians work for adaptive refinement
+
+  //  unsigned int counter = 0;
+  //  for (auto cell = tria.begin_active(); cell != tria.end(); ++cell,
+  //  ++counter)
+  //    if (cell->is_locally_owned() && counter % 3 == 0)
+  //      cell->set_refine_flag();
+  //  tria.execute_coarsening_and_refinement();
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+  MappingQGeneric<dim> mapping(1);
+
+  AffineConstraints<double> constraints;
+  VectorTools::interpolate_boundary_values(
+    mapping, dof_handler, 0, Functions::ZeroFunction<dim>(), constraints);
+  constraints.close();
+
+  // FEFaceEvaluation
+  typename MatrixFree<dim, double, VectorizedArrayType>::AdditionalData
+    additional_data;
+  additional_data.mapping_update_flags_inner_faces =
+    update_values | update_gradients | update_hessians;
+
+  MatrixFree<dim, double, VectorizedArrayType> matrix_free;
+  matrix_free.reinit(
+    mapping, dof_handler, constraints, QGauss<1>(degree + 1), additional_data);
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    deallog << "Working with " << fe.get_name() << " and "
+            << dof_handler.n_dofs() << " dofs" << std::endl;
+
+  LinearAlgebra::distributed::Vector<double> src, dst, dst2;
+  matrix_free.initialize_dof_vector(src);
+  for (auto &v : src)
+    v = random_value<double>();
+
+  matrix_free.initialize_dof_vector(dst);
+  matrix_free.initialize_dof_vector(dst2);
+
+  // Setup FEFaceValues
+  QGauss<dim - 1>                      quad(degree + 1);
+  FEFaceValues<dim>                    fe_face_values_m(mapping,
+                                     fe,
+                                     quad,
+                                     update_values | update_gradients |
+                                       update_hessians | update_JxW_values);
+  FEFaceValues<dim>                    fe_face_values_p(mapping,
+                                     fe,
+                                     quad,
+                                     update_values | update_gradients |
+                                       update_hessians | update_JxW_values);
+  Vector<double>                       solution_values_m(fe.dofs_per_cell);
+  Vector<double>                       solution_values_p(fe.dofs_per_cell);
+  std::vector<Tensor<2, dim>>          solution_hessians(quad.size());
+  std::vector<Tensor<1, dim>>          solution_gradients(quad.size());
+  std::vector<double>                  solution_values(quad.size());
+  std::vector<types::global_dof_index> dof_indices_m(fe.dofs_per_cell);
+  std::vector<types::global_dof_index> dof_indices_p(fe.dofs_per_cell);
+  src.update_ghost_values();
+  dst2 = 0;
+
+
+
+  matrix_free.template loop<LinearAlgebra::distributed::Vector<double>,
+                            LinearAlgebra::distributed::Vector<double>>(
+    [&](
+      const auto &matrix_free, auto &dst, const auto &src, const auto &range) {
+      (void)matrix_free;
+      (void)dst;
+      (void)src;
+      (void)range;
+    },
+    [&](
+      const auto &matrix_free, auto &dst, const auto &src, const auto &range) {
+      FEFaceEvaluation<dim, -1, 0, 1, double, VectorizedArrayType> fe_eval_m(
+        matrix_free, true);
+      FEFaceEvaluation<dim, -1, 0, 1, double, VectorizedArrayType> fe_eval_p(
+        matrix_free, false);
+      for (unsigned int face = range.first; face < range.second; ++face)
+        {
+          // FEFaceEvaluation
+          fe_eval_m.reinit(face);
+          fe_eval_m.gather_evaluate(src, evaluation_flags);
+          for (unsigned int q = 0; q < fe_eval_m.n_q_points; ++q)
+            {
+              if (evaluation_flags & EvaluationFlags::hessians)
+                fe_eval_m.submit_hessian(fe_eval_m.get_hessian(q), q);
+              if (evaluation_flags & EvaluationFlags::gradients)
+                fe_eval_m.submit_gradient(fe_eval_m.get_gradient(q), q);
+              if (evaluation_flags & EvaluationFlags::values)
+                fe_eval_m.submit_value(fe_eval_m.get_value(q), q);
+            }
+          fe_eval_m.integrate(evaluation_flags);
+          fe_eval_m.distribute_local_to_global(dst);
+
+          fe_eval_p.reinit(face);
+          fe_eval_p.gather_evaluate(src, evaluation_flags);
+          for (unsigned int q = 0; q < fe_eval_p.n_q_points; ++q)
+            {
+              if (evaluation_flags & EvaluationFlags::hessians)
+                fe_eval_p.submit_hessian(fe_eval_p.get_hessian(q), q);
+              if (evaluation_flags & EvaluationFlags::gradients)
+                fe_eval_p.submit_gradient(fe_eval_p.get_gradient(q), q);
+              if (evaluation_flags & EvaluationFlags::values)
+                fe_eval_p.submit_value(fe_eval_p.get_value(q), q);
+            }
+          fe_eval_p.integrate(evaluation_flags);
+          fe_eval_p.distribute_local_to_global(dst);
+
+          // FEFaceValues
+          for (unsigned int v = 0;
+               v < matrix_free.n_active_entries_per_face_batch(face);
+               ++v)
+            {
+              const auto face_accessor_inside =
+                matrix_free.get_face_iterator(face, v, true);
+              const auto face_accessor_outside =
+                matrix_free.get_face_iterator(face, v, false);
+
+              fe_face_values_m.reinit(face_accessor_inside.first,
+                                      face_accessor_inside.second);
+
+              fe_face_values_p.reinit(face_accessor_outside.first,
+                                      face_accessor_outside.second);
+
+              face_accessor_inside.first->get_dof_indices(dof_indices_m);
+              face_accessor_outside.first->get_dof_indices(dof_indices_p);
+              constraints.get_dof_values(src,
+                                         dof_indices_m.begin(),
+                                         solution_values_m.begin(),
+                                         solution_values_m.end());
+
+              constraints.get_dof_values(src,
+                                         dof_indices_p.begin(),
+                                         solution_values_p.begin(),
+                                         solution_values_p.end());
+
+              for (unsigned int q = 0; q < quad.size(); ++q)
+                {
+                  double         values = 0.;
+                  Tensor<1, dim> gradients;
+                  Tensor<2, dim> hessians;
+                  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
+                    {
+                      if (evaluation_flags & EvaluationFlags::hessians)
+                        hessians += solution_values_m(i) *
+                                    fe_face_values_m.shape_hessian(i, q);
+                      if (evaluation_flags & EvaluationFlags::gradients)
+                        gradients += solution_values_m(i) *
+                                     fe_face_values_m.shape_grad(i, q);
+                      if (evaluation_flags & EvaluationFlags::values)
+                        values += solution_values_m(i) *
+                                  fe_face_values_m.shape_value(i, q);
+                    }
+                  solution_hessians[q]  = hessians * fe_face_values_m.JxW(q);
+                  solution_gradients[q] = gradients * fe_face_values_m.JxW(q);
+                  solution_values[q]    = values * fe_face_values_m.JxW(q);
+                }
+              for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
+                {
+                  double sum_hessians  = 0.;
+                  double sum_gradients = 0.;
+                  double sum_values    = 0.;
+                  for (unsigned int q = 0; q < quad.size(); ++q)
+                    {
+                      if (evaluation_flags & EvaluationFlags::hessians)
+                        sum_hessians += double_contract<0, 0, 1, 1>(
+                          solution_hessians[q],
+                          fe_face_values_m.shape_hessian(i, q));
+                      if (evaluation_flags & EvaluationFlags::gradients)
+                        sum_gradients += solution_gradients[q] *
+                                         fe_face_values_m.shape_grad(i, q);
+                      if (evaluation_flags & EvaluationFlags::values)
+                        sum_values += solution_values[q] *
+                                      fe_face_values_m.shape_value(i, q);
+                    }
+                  solution_values_m(i) =
+                    sum_hessians + sum_gradients + sum_values;
+                }
+              constraints.distribute_local_to_global(solution_values_m,
+                                                     dof_indices_m,
+                                                     dst2);
+
+              for (unsigned int q = 0; q < quad.size(); ++q)
+                {
+                  double         values = 0.;
+                  Tensor<1, dim> gradients;
+                  Tensor<2, dim> hessians;
+                  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
+                    {
+                      if (evaluation_flags & EvaluationFlags::hessians)
+                        hessians += solution_values_p(i) *
+                                    fe_face_values_p.shape_hessian(i, q);
+                      if (evaluation_flags & EvaluationFlags::gradients)
+                        gradients += solution_values_p(i) *
+                                     fe_face_values_p.shape_grad(i, q);
+                      if (evaluation_flags & EvaluationFlags::values)
+                        values += solution_values_p(i) *
+                                  fe_face_values_p.shape_value(i, q);
+                    }
+                  solution_hessians[q]  = hessians * fe_face_values_p.JxW(q);
+                  solution_gradients[q] = gradients * fe_face_values_p.JxW(q);
+                  solution_values[q]    = values * fe_face_values_p.JxW(q);
+                }
+              for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
+                {
+                  double sum_hessians  = 0.;
+                  double sum_gradients = 0.;
+                  double sum_values    = 0.;
+                  for (unsigned int q = 0; q < quad.size(); ++q)
+                    {
+                      if (evaluation_flags & EvaluationFlags::hessians)
+                        sum_hessians += double_contract<0, 0, 1, 1>(
+                          solution_hessians[q],
+                          fe_face_values_p.shape_hessian(i, q));
+                      if (evaluation_flags & EvaluationFlags::gradients)
+                        sum_gradients += solution_gradients[q] *
+                                         fe_face_values_p.shape_grad(i, q);
+                      if (evaluation_flags & EvaluationFlags::values)
+                        sum_values += solution_values[q] *
+                                      fe_face_values_p.shape_value(i, q);
+                    }
+                  solution_values_p(i) =
+                    sum_hessians + sum_gradients + sum_values;
+                }
+              constraints.distribute_local_to_global(solution_values_p,
+                                                     dof_indices_p,
+                                                     dst2);
+            }
+        }
+    },
+    [&](
+      const auto &matrix_free, auto &dst, const auto &src, const auto &range) {
+      (void)matrix_free;
+      (void)dst;
+      (void)src;
+      (void)range;
+    },
+    dst,
+    src,
+    true);
+
+  const unsigned int end_of_print_dst =
+    dof_handler.n_dofs() > 9 ? 9 : dof_handler.n_dofs();
+
+  deallog << "dst FEE: ";
+  for (unsigned int i = 0; i < end_of_print_dst; ++i)
+    deallog << dst[i] << " ";
+  deallog << std::endl;
+
+  deallog << "dst FEV: ";
+  for (unsigned int i = 0; i < end_of_print_dst; ++i)
+    deallog << dst2[i] << " ";
+  deallog << std::endl;
+
+  // compare solutions of matrix vector product
+  {
+    dst2 -= dst;
+
+    double error = 0.;
+    if (dst.l2_norm() > 0)
+      error = dst2.l2_norm() / dst.l2_norm();
+    else
+      error = dst2.l2_norm();
+
+    if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+      deallog << "FEValues verification: " << error << std::endl << std::endl;
+  }
+}
+
+int
+main(int argc, char **argv)
+{
+  dealii::Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+
+  initlog();
+  deallog << std::setprecision(10);
+
+  {
+    dealii::EvaluationFlags::EvaluationFlags evaluation_flags =
+      EvaluationFlags::hessians;
+    deallog << "test_hessians_only" << std::endl;
+    for (unsigned int i = 1; i < 4; ++i)
+      {
+        test_hessians<1>(i, dealii::FE_Q<1>(i), evaluation_flags);
+        test_hessians<2>(i, dealii::FE_Q<2>(i), evaluation_flags);
+        test_hessians<3>(i, dealii::FE_Q<3>(i), evaluation_flags);
+        test_hessians<1>(i, dealii::FE_DGQ<1>(i), evaluation_flags);
+        test_hessians<2>(i, dealii::FE_DGQ<2>(i), evaluation_flags);
+        test_hessians<3>(i, dealii::FE_DGQ<3>(i), evaluation_flags);
+      }
+  }
+
+  {
+    dealii::EvaluationFlags::EvaluationFlags evaluation_flags =
+      EvaluationFlags::values | EvaluationFlags::hessians;
+    deallog << "test_hessians_with_values" << std::endl;
+    for (unsigned int i = 1; i < 4; ++i)
+      {
+        test_hessians<1>(i, dealii::FE_Q<1>(i), evaluation_flags);
+        test_hessians<2>(i, dealii::FE_Q<2>(i), evaluation_flags);
+        test_hessians<3>(i, dealii::FE_Q<3>(i), evaluation_flags);
+        test_hessians<1>(i, dealii::FE_DGQ<1>(i), evaluation_flags);
+        test_hessians<2>(i, dealii::FE_DGQ<2>(i), evaluation_flags);
+        test_hessians<3>(i, dealii::FE_DGQ<3>(i), evaluation_flags);
+      }
+  }
+
+  {
+    dealii::EvaluationFlags::EvaluationFlags evaluation_flags =
+      EvaluationFlags::gradients | EvaluationFlags::hessians;
+    deallog << "test_hessians_with_gradients" << std::endl;
+    for (unsigned int i = 1; i < 4; ++i)
+      {
+        test_hessians<1>(i, dealii::FE_Q<1>(i), evaluation_flags);
+        test_hessians<2>(i, dealii::FE_Q<2>(i), evaluation_flags);
+        test_hessians<3>(i, dealii::FE_Q<3>(i), evaluation_flags);
+        test_hessians<1>(i, dealii::FE_DGQ<1>(i), evaluation_flags);
+        test_hessians<2>(i, dealii::FE_DGQ<2>(i), evaluation_flags);
+        test_hessians<3>(i, dealii::FE_DGQ<3>(i), evaluation_flags);
+      }
+  }
+
+  {
+    dealii::EvaluationFlags::EvaluationFlags evaluation_flags =
+      EvaluationFlags::values | EvaluationFlags::gradients |
+      EvaluationFlags::hessians;
+    deallog << "test_hessians_with_gradients_and_values" << std::endl;
+    for (unsigned int i = 1; i < 4; ++i)
+      {
+        test_hessians<1>(i, dealii::FE_Q<1>(i), evaluation_flags);
+        test_hessians<2>(i, dealii::FE_Q<2>(i), evaluation_flags);
+        test_hessians<3>(i, dealii::FE_Q<3>(i), evaluation_flags);
+        test_hessians<1>(i, dealii::FE_DGQ<1>(i), evaluation_flags);
+        test_hessians<2>(i, dealii::FE_DGQ<2>(i), evaluation_flags);
+        test_hessians<3>(i, dealii::FE_DGQ<3>(i), evaluation_flags);
+      }
+  }
+}
diff --git a/tests/matrix_free/matrix_vector_hessians_faces.output b/tests/matrix_free/matrix_vector_hessians_faces.output
new file mode 100644 (file)
index 0000000..a3ad9ef
--- /dev/null
@@ -0,0 +1,365 @@
+
+DEAL::test_hessians_only
+DEAL::Working with FE_Q<1>(1) and 9 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<2>(1) and 81 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 -3961.844511 0.000000000 3975.516664 0.000000000 2180.676992 4279.873135 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 -3961.844511 0.000000000 3975.516664 0.000000000 2180.676992 4279.873135 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 729 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 2671.158976 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 2671.158976 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<1>(1) and 16 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 256 dofs
+DEAL::dst FEE: 186.7674628 -186.7674628 -186.7674628 186.7674628 1095.845132 -1095.845132 -1095.845132 1095.845132 -1770.726878 
+DEAL::dst FEV: 186.7674628 -186.7674628 -186.7674628 186.7674628 1095.845132 -1095.845132 -1095.845132 1095.845132 -1770.726878 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 4096 dofs
+DEAL::dst FEE: -41.82729979 116.7876889 -72.84484080 -2.115548314 24.66611142 -99.62650053 90.00602917 -15.04564006 -135.6029421 
+DEAL::dst FEV: -41.82729979 116.7876889 -72.84484080 -2.115548314 24.66611142 -99.62650053 90.00602917 -15.04564006 -135.6029421 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<1>(2) and 17 dofs
+DEAL::dst FEE: 0.000000000 67338.98941 16221.91627 174612.5300 -150899.8951 93741.61943 -198325.1649 100444.2774 10841.92603 
+DEAL::dst FEV: 0.000000000 67338.98941 16221.91627 174612.5300 -150899.8951 93741.61943 -198325.1649 100444.2774 10841.92603 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 289 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 -104920.8906 0.000000000 -49521.70596 0.000000000 150903.6170 6452.694464 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 -104920.8906 0.000000000 -49521.70596 0.000000000 150903.6170 6452.694464 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 4913 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 31698.05993 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 31698.05993 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<1>(2) and 24 dofs
+DEAL::dst FEE: -57370.44455 114740.8891 -57370.44455 -130179.7286 260359.4573 -130179.7286 -132118.6579 264237.3158 -132118.6579 
+DEAL::dst FEV: -57370.44455 114740.8891 -57370.44455 -130179.7286 260359.4573 -130179.7286 -132118.6579 264237.3158 -132118.6579 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 576 dofs
+DEAL::dst FEE: 2621.950491 -1575.677282 -3332.535102 -31645.71115 63719.74128 -27501.50635 24104.29342 -52305.12954 25914.57422 
+DEAL::dst FEV: 2621.950491 -1575.677282 -3332.535102 -31645.71115 63719.74128 -27501.50635 24104.29342 -52305.12954 25914.57422 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 13824 dofs
+DEAL::dst FEE: 50.88230764 -830.0219522 -818.2722402 -1019.544361 3222.190645 -744.7435660 811.4995692 1670.530726 -893.6889450 
+DEAL::dst FEV: 50.88230764 -830.0219522 -818.2722402 -1019.544361 3222.190645 -744.7435660 811.4995692 1670.530726 -893.6889450 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<1>(3) and 25 dofs
+DEAL::dst FEE: 0.000000000 -737280.1891 -22080.92721 29817.48585 -1201734.837 1302772.885 -799181.1739 -1342610.859 3103193.420 
+DEAL::dst FEV: 0.000000000 -737280.1891 -22080.92721 29817.48585 -1201734.837 1302772.885 -799181.1739 -1342610.859 3103193.420 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 625 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 -1547335.199 0.000000000 0.000000000 368454.5924 67296.69792 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 -1547335.199 0.000000000 0.000000000 368454.5924 67296.69792 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 15625 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 70210.21641 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 70210.21641 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<1>(3) and 32 dofs
+DEAL::dst FEE: -254241.5072 725631.1855 -979872.6927 508483.0143 369638.5526 -633813.2719 321979.8414 -57805.12216 -296675.7776 
+DEAL::dst FEV: -254241.5072 725631.1855 -979872.6927 508483.0143 369638.5526 -633813.2719 321979.8414 -57805.12216 -296675.7776 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 1024 dofs
+DEAL::dst FEE: 19404.30284 11401.17609 53539.65109 -156046.9045 -56073.49671 -37352.39201 -103984.8509 359118.5689 51864.00926 
+DEAL::dst FEV: 19404.30284 11401.17609 53539.65109 -156046.9045 -56073.49671 -37352.39201 -103984.8509 359118.5689 51864.00926 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 32768 dofs
+DEAL::dst FEE: 118.4885551 -4107.177606 1812.382519 -363.6809581 -2068.812140 9700.899062 -6925.304917 4256.899099 -2307.998977 
+DEAL::dst FEV: 118.4885551 -4107.177606 1812.382519 -363.6809581 -2068.812140 9700.899062 -6925.304917 4256.899099 -2307.998977 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::test_hessians_with_values
+DEAL::Working with FE_Q<1>(1) and 9 dofs
+DEAL::dst FEE: 0.000000000 1.485241358 0.4717695818 1.336048384 0.09186878712 0.5541558529 0.6274125085 0.1848661360 0.000000000 
+DEAL::dst FEV: 0.000000000 1.485241358 0.4717695818 1.336048384 0.09186878712 0.5541558529 0.6274125085 0.1848661360 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<2>(1) and 81 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 1070.969383 0.000000000 2878.096119 0.000000000 -3268.391675 -412.5247158 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 1070.969383 0.000000000 2878.096119 0.000000000 -3268.391675 -412.5247158 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 729 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 3500.106392 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 3500.106392 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<1>(1) and 16 dofs
+DEAL::dst FEE: 0.000000000 0.9178876015 0.3258156191 0.6197056345 0.1537631360 0.2945842553 0.3555839413 0.1329097907 0.7070026150 
+DEAL::dst FEV: 0.000000000 0.9178876015 0.3258156191 0.6197056345 0.1537631360 0.2945842553 0.3555839413 0.1329097907 0.7070026150 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 256 dofs
+DEAL::dst FEE: -569.5871632 569.6347435 569.6092785 -569.5287350 -1511.130148 1511.164480 1511.202041 -1511.108327 -2579.882528 
+DEAL::dst FEV: -569.5871632 569.6347435 569.6092785 -569.5287350 -1511.130148 1511.164480 1511.202041 -1511.108327 -2579.882528 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 4096 dofs
+DEAL::dst FEE: 7.424074176 241.6432657 -274.2813513 25.22374027 79.00171305 -328.0612417 187.8648553 61.21168368 78.30753564 
+DEAL::dst FEV: 7.424074176 241.6432657 -274.2813513 25.22374027 79.00171305 -328.0612417 187.8648553 61.21168368 78.30753564 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<1>(2) and 17 dofs
+DEAL::dst FEE: 0.000000000 -167135.6577 229538.9869 -97797.09182 104732.8852 -39369.24028 90862.39763 -109611.4362 -12122.77601 
+DEAL::dst FEV: 0.000000000 -167135.6577 229538.9869 -97797.09182 104732.8852 -39369.24028 90862.39763 -109611.4362 -12122.77601 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 289 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 60209.63375 0.000000000 69766.67093 0.000000000 -58791.04965 -39753.16339 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 60209.63375 0.000000000 69766.67093 0.000000000 -58791.04965 -39753.16339 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 4913 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 12306.73607 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 12306.73607 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<1>(2) and 24 dofs
+DEAL::dst FEE: 47813.68795 -95627.37590 47814.42813 -96189.12898 192379.2740 -96189.51817 55225.52064 -110450.7480 55226.32803 
+DEAL::dst FEV: 47813.68795 -95627.37590 47814.42813 -96189.12898 192379.2740 -96189.51817 55225.52064 -110450.7480 55226.32803 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 576 dofs
+DEAL::dst FEE: -5287.684869 6572.726856 -20184.16775 24229.94350 -48365.22252 61933.61384 -25037.39476 53982.83575 -47844.48685 
+DEAL::dst FEV: -5287.684869 6572.726856 -20184.16775 24229.94350 -48365.22252 61933.61384 -25037.39476 53982.83575 -47844.48685 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 13824 dofs
+DEAL::dst FEE: 370.0566403 -292.3081868 -720.4052507 -622.7558219 238.6508182 1862.634403 281.8888045 2091.765602 -1871.063632 
+DEAL::dst FEV: 370.0566403 -292.3081868 -720.4052507 -622.7558219 238.6508182 1862.634403 281.8888045 2091.765602 -1871.063632 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<1>(3) and 25 dofs
+DEAL::dst FEE: 0.000000000 299377.0807 289587.9028 -391051.6578 -3248169.708 -134872.5658 4149.691300 2108415.055 7239160.344 
+DEAL::dst FEV: 0.000000000 299377.0807 289587.9028 -391051.6578 -3248169.708 -134872.5658 4149.691300 2108415.055 7239160.344 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 625 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 -293073.3740 0.000000000 0.000000000 675561.0838 -427952.8498 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 -293073.3740 0.000000000 0.000000000 675561.0838 -427952.8498 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 15625 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 -85712.86767 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 -85712.86767 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<1>(3) and 32 dofs
+DEAL::dst FEE: -252775.6199 721447.3938 -974223.0137 505551.7303 1611958.374 -3761443.564 4015469.154 -1865982.968 -51596.77459 
+DEAL::dst FEV: -252775.6199 721447.3938 -974223.0137 505551.7303 1611958.374 -3761443.564 4015469.154 -1865982.968 -51596.77459 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 1024 dofs
+DEAL::dst FEE: -10040.20263 53844.53594 -182876.7764 97992.24016 38492.66295 -176621.9667 403072.9022 -164030.7460 150191.7481 
+DEAL::dst FEV: -10040.20263 53844.53594 -182876.7764 97992.24016 38492.66295 -176621.9667 403072.9022 -164030.7460 150191.7481 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 32768 dofs
+DEAL::dst FEE: -674.3169640 814.0460283 -4216.899024 5057.280398 2850.106346 -2615.447288 14600.65907 -8650.675187 -3053.718113 
+DEAL::dst FEV: -674.3169640 814.0460283 -4216.899024 5057.280398 2850.106346 -2615.447288 14600.65907 -8650.675187 -3053.718113 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::test_hessians_with_gradients
+DEAL::Working with FE_Q<1>(1) and 9 dofs
+DEAL::dst FEE: 0.000000000 64.15179569 36.90356960 -40.28644325 -68.37967620 39.09859355 46.88549312 -62.19046491 31.72474919 
+DEAL::dst FEV: 0.000000000 64.15179569 36.90356960 -40.28644325 -68.37967620 39.09859355 46.88549312 -62.19046491 31.72474919 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<2>(1) and 81 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 442.0107683 0.000000000 3542.943855 0.000000000 -9575.588132 7499.031588 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 442.0107683 0.000000000 3542.943855 0.000000000 -9575.588132 7499.031588 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 729 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 -2821.551619 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 -2821.551619 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<1>(1) and 16 dofs
+DEAL::dst FEE: 10.06237412 -10.06237412 -27.76560022 27.76560022 34.04355665 -34.04355665 -21.65914739 21.65914739 -7.379041318 
+DEAL::dst FEV: 10.06237412 -10.06237412 -27.76560022 27.76560022 34.04355665 -34.04355665 -21.65914739 21.65914739 -7.379041318 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 256 dofs
+DEAL::dst FEE: -1170.185020 1177.791704 1171.983635 -1179.590319 -370.5322257 366.0001104 361.6813287 -357.1492134 373.2952228 
+DEAL::dst FEV: -1170.185020 1177.791704 1171.983635 -1179.590319 -370.5322257 366.0001104 361.6813287 -357.1492134 373.2952228 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 4096 dofs
+DEAL::dst FEE: -21.58586802 226.0914318 181.3343315 -385.5024982 28.95060815 -234.0415895 -189.5067703 394.2603544 -219.6450550 
+DEAL::dst FEV: -21.58586802 226.0914318 181.3343315 -385.5024982 28.95060815 -234.0415895 -189.5067703 394.2603544 -219.6450550 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<1>(2) and 17 dofs
+DEAL::dst FEE: 0.000000000 22439.52331 189626.7649 11978.51334 -234916.2278 -265869.0386 211155.8797 -262006.1754 320433.1750 
+DEAL::dst FEV: 0.000000000 22439.52331 189626.7649 11978.51334 -234916.2278 -265869.0386 211155.8797 -262006.1754 320433.1750 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 289 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 -10856.65148 0.000000000 -25321.75472 0.000000000 95280.37568 4511.212957 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 -10856.65148 0.000000000 -25321.75472 0.000000000 95280.37568 4511.212957 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 4913 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 9080.423198 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 9080.423198 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<1>(2) and 24 dofs
+DEAL::dst FEE: 9696.728295 -19428.05746 9731.329162 -57259.37641 114443.5555 -57184.17906 159527.5899 -319143.7449 159616.1550 
+DEAL::dst FEV: 9696.728295 -19428.05746 9731.329162 -57259.37641 114443.5555 -57184.17906 159527.5899 -319143.7449 159616.1550 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 576 dofs
+DEAL::dst FEE: 10194.28422 -35453.44115 24293.82892 -45055.92628 130503.0063 -83489.42850 30711.39427 -86727.69612 55023.97831 
+DEAL::dst FEV: 10194.28422 -35453.44115 24293.82892 -45055.92628 130503.0063 -83489.42850 30711.39427 -86727.69612 55023.97831 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 13824 dofs
+DEAL::dst FEE: 448.1446413 -510.9189963 1634.494590 755.7615764 -6223.772600 800.6586105 529.4850620 1959.380310 711.9024984 
+DEAL::dst FEV: 448.1446413 -510.9189963 1634.494590 755.7615764 -6223.772600 800.6586105 529.4850620 1959.380310 711.9024984 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<1>(3) and 25 dofs
+DEAL::dst FEE: 0.000000000 2252292.265 1573930.091 -2126197.154 -174516.4842 -2501369.733 2393718.162 -640706.0941 -2054230.489 
+DEAL::dst FEV: 0.000000000 2252292.265 1573930.091 -2126197.154 -174516.4842 -2501369.733 2393718.162 -640706.0941 -2054230.489 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 625 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 -1317837.470 0.000000000 0.000000000 -264992.6854 1082118.987 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 -1317837.470 0.000000000 0.000000000 -264992.6854 1082118.987 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 15625 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 -65742.12001 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 -65742.12001 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<1>(3) and 32 dofs
+DEAL::dst FEE: 607362.1400 -1733480.886 2340968.548 -1214849.802 -1047359.729 2045725.963 -1566917.253 568551.0199 -3864497.016 
+DEAL::dst FEV: 607362.1400 -1733480.886 2340968.548 -1214849.802 -1047359.729 2045725.963 -1566917.253 568551.0199 -3864497.016 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 1024 dofs
+DEAL::dst FEE: -27040.61278 3573.267393 57805.44291 58432.19141 166809.4545 -257530.0470 225343.1275 -358312.5442 -334186.2219 
+DEAL::dst FEV: -27040.61278 3573.267393 57805.44291 58432.19141 166809.4545 -257530.0470 225343.1275 -358312.5442 -334186.2219 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 32768 dofs
+DEAL::dst FEE: -464.3635240 2028.142271 763.8205053 -549.6966376 327.0786350 -596.5268953 -20821.36159 10310.01212 -817.8506617 
+DEAL::dst FEV: -464.3635240 2028.142271 763.8205053 -549.6966376 327.0786350 -596.5268953 -20821.36159 10310.01212 -817.8506617 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::test_hessians_with_gradients_and_values
+DEAL::Working with FE_Q<1>(1) and 9 dofs
+DEAL::dst FEE: 0.000000000 -37.00478557 48.59625945 0.8581540048 -7.323365868 59.63530424 -26.02881961 -15.72262317 -0.3217377366 
+DEAL::dst FEV: 0.000000000 -37.00478557 48.59625945 0.8581540048 -7.323365868 59.63530424 -26.02881961 -15.72262317 -0.3217377366 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<2>(1) and 81 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 9507.556791 0.000000000 -7321.070916 0.000000000 -12431.59166 12502.87908 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 9507.556791 0.000000000 -7321.070916 0.000000000 -12431.59166 12502.87908 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 729 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 280.4491881 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 280.4491881 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<1>(1) and 16 dofs
+DEAL::dst FEE: 39.29021196 -39.01623275 -9.963182797 11.39805324 -15.84908361 16.43351670 80.59747247 -79.58483093 53.75457114 
+DEAL::dst FEV: 39.29021196 -39.01623275 -9.963182797 11.39805324 -15.84908361 16.43351670 80.59747247 -79.58483093 53.75457114 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 256 dofs
+DEAL::dst FEE: 1516.712105 -1527.764366 -1523.902913 1535.086538 -1666.971584 1682.840844 1665.583768 -1681.331729 -76.91218429 
+DEAL::dst FEV: 1516.712105 -1527.764366 -1523.902913 1535.086538 -1666.971584 1682.840844 1665.583768 -1681.331729 -76.91218429 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 4096 dofs
+DEAL::dst FEE: -1.189912404 234.8477478 180.7860225 -413.6249483 -121.5440256 -112.1298409 -57.51853159 290.3907440 289.0256656 
+DEAL::dst FEV: -1.189912404 234.8477478 180.7860225 -413.6249483 -121.5440256 -112.1298409 -57.51853159 290.3907440 289.0256656 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<1>(2) and 17 dofs
+DEAL::dst FEE: 0.000000000 29141.70579 65351.23242 -37541.71163 -123874.9385 15975.67678 199194.1660 179480.0205 -231343.3646 
+DEAL::dst FEV: 0.000000000 29141.70579 65351.23242 -37541.71163 -123874.9385 15975.67678 199194.1660 179480.0205 -231343.3646 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 289 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 5622.410746 0.000000000 -124785.4497 0.000000000 -28083.28907 140221.5951 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 5622.410746 0.000000000 -124785.4497 0.000000000 -28083.28907 140221.5951 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 4913 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 4030.831352 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 4030.831352 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<1>(2) and 24 dofs
+DEAL::dst FEE: 1260.253162 -2590.873853 1331.260025 -26516.06494 53162.96332 -26645.70214 192196.9013 -384333.9650 192138.6030 
+DEAL::dst FEV: 1260.253162 -2590.873853 1331.260025 -26516.06494 53162.96332 -26645.70214 192196.9013 -384333.9650 192138.6030 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 576 dofs
+DEAL::dst FEE: -5169.333357 29063.34663 -23105.45307 24704.96903 -97773.79532 71498.06891 -4767.981572 39154.65283 -33604.35418 
+DEAL::dst FEV: -5169.333357 29063.34663 -23105.45307 24704.96903 -97773.79532 71498.06891 -4767.981572 39154.65283 -33604.35418 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 13824 dofs
+DEAL::dst FEE: -95.58676656 1276.489009 -1334.575175 274.9405999 -4042.776356 2530.457071 1137.743263 1405.963407 -1676.022846 
+DEAL::dst FEV: -95.58676656 1276.489009 -1334.575175 274.9405999 -4042.776356 2530.457071 1137.743263 1405.963407 -1676.022846 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<1>(3) and 25 dofs
+DEAL::dst FEE: 0.000000000 -3677129.436 -674168.6624 910870.7211 167797.2252 7072904.503 -6925910.014 1212429.567 6407996.850 
+DEAL::dst FEV: 0.000000000 -3677129.436 -674168.6624 910870.7211 167797.2252 7072904.503 -6925910.014 1212429.567 6407996.850 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 625 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 901522.2209 0.000000000 0.000000000 497122.6409 -783742.6748 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 901522.2209 0.000000000 0.000000000 497122.6409 -783742.6748 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 15625 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 61931.57932 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 61931.57932 0.000000000 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<1>(3) and 32 dofs
+DEAL::dst FEE: -561890.2156 1603751.370 -2166648.089 1124787.675 -928689.6365 1915272.283 -1654708.155 668125.8408 -2711210.037 
+DEAL::dst FEV: -561890.2156 1603751.370 -2166648.089 1124787.675 -928689.6365 1915272.283 -1654708.155 668125.8408 -2711210.037 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 1024 dofs
+DEAL::dst FEE: -2674.684464 105342.7462 -38524.91232 -48556.52335 48969.79796 -415911.0458 256625.4211 38926.74483 -131829.3972 
+DEAL::dst FEV: -2674.684464 105342.7462 -38524.91232 -48556.52335 48969.79796 -415911.0458 256625.4211 38926.74483 -131829.3972 
+DEAL::FEValues verification: 0.000000000
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 32768 dofs
+DEAL::dst FEE: 692.2851052 972.7257414 820.4660691 1018.320908 -949.9325546 -10318.08673 14538.53642 -19371.16378 5741.320571 
+DEAL::dst FEV: 692.2851052 972.7257414 820.4660691 1018.320908 -949.9325546 -10318.08673 14538.53642 -19371.16378 5741.320571 
+DEAL::FEValues verification: 0.000000000
+DEAL::

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.