]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Hessains on general cells and faces 12744/head
authorMaximilian Bergbauer <bergbauer@lnm.mw.tum.de>
Tue, 24 Aug 2021 07:40:34 +0000 (09:40 +0200)
committerMaximilian Bergbauer <bergbauer@lnm.mw.tum.de>
Tue, 7 Sep 2021 09:08:58 +0000 (11:08 +0200)
doc/news/changes/major/20210906MaximilianBergbauer [moved from doc/news/changes/minor/20210728MaximilianBergbauer with 64% similarity]
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/mapping_info.templates.h
tests/matrix_free/integrate_functions_hessians.cc [new file with mode: 0644]
tests/matrix_free/integrate_functions_hessians.output [new file with mode: 0644]
tests/matrix_free/matrix_vector_hessians_general_cells.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_hessians_general_cells.output [new file with mode: 0644]
tests/matrix_free/matrix_vector_hessians_general_faces.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_hessians_general_faces.output [new file with mode: 0644]

similarity index 64%
rename from doc/news/changes/minor/20210728MaximilianBergbauer
rename to doc/news/changes/major/20210906MaximilianBergbauer
index 1ccf8c26f87549caa9b6b1d5b2f47671f2daec8e..22dbc6a9aa11c7b8d6f659a4eb231759c3fdeaf8 100644 (file)
@@ -1,4 +1,4 @@
 Added: Evaluation and integration of hessians both on cells and on faces is now
 available in the matrix-free framework.
 <br>
-(Maximilian Bergbauer, 2021/07/28)
+(Maximilian Bergbauer, Martin Kronbichler, Peter Munch, 2021/09/06)
index 50e61feb1f001ea094eca93d595dddf61141f13f..972e122f09204c64ba808a18fad2c64da4583e5f 100644 (file)
@@ -1315,6 +1315,21 @@ protected:
    */
   VectorizedArrayType *gradients_quad;
 
+  /**
+   * This field stores the gradients of the finite element function on
+   * quadrature points after applying unit cell transformations or before
+   * integrating. The methods get_hessian() and submit_hessian() (as well as
+   * some specializations like get_hessian_diagonal() or get_laplacian())
+   * access this field for general cell/face types.
+   *
+   * The values of this array are stored in the start section of
+   * @p scratch_data_array. Due to its access as a thread local memory, the
+   * memory can get reused between different calls. As opposed to requesting
+   * memory on the stack, this approach allows for very large polynomial
+   * degrees.
+   */
+  VectorizedArrayType *gradients_from_hessians_quad;
+
   /**
    * This field stores the Hessians of the finite element function on
    * quadrature points after applying unit cell transformations. The methods
@@ -4330,7 +4345,8 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
     2 * n_quadrature_points;
   const unsigned int allocated_size =
     shift + n_components_ * dofs_per_component +
-    (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
+    (n_components_ * ((dim * (dim + 1)) / 2 + 2 * dim + 1) *
+     n_quadrature_points);
   this->scratch_data_array->resize_fast(allocated_size);
 
   // set the pointers to the correct position in the data array
@@ -4343,12 +4359,16 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
     this->scratch_data_array->begin() + n_components * dofs_per_component;
   gradients_quad = this->scratch_data_array->begin() +
                    n_components * (dofs_per_component + n_quadrature_points);
-  hessians_quad =
+  gradients_from_hessians_quad =
     this->scratch_data_array->begin() +
     n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
-  this->scratch_data =
-    this->scratch_data_array->begin() + n_components_ * dofs_per_component +
-    (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
+  hessians_quad =
+    this->scratch_data_array->begin() +
+    n_components * (dofs_per_component + (2 * dim + 1) * n_quadrature_points);
+  this->scratch_data = this->scratch_data_array->begin() +
+                       n_components_ * dofs_per_component +
+                       (n_components_ * ((dim * (dim + 1)) / 2 + 2 * dim + 1) *
+                        n_quadrature_points);
 }
 
 
@@ -6019,8 +6039,6 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
           [this->mapping_data->data_index_offsets[this->cell] + q_point];
       for (unsigned int comp = 0; comp < n_components; ++comp)
         {
-          // compute laplacian before the gradient because it needs to access
-          // unscaled gradient data
           VectorizedArrayType tmp[dim][dim];
           internal::hessian_unit_times_jac(
             jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
@@ -6039,7 +6057,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
             for (unsigned int e = 0; e < dim; ++e)
               hessian_out[comp][d][d] +=
                 jac_grad[d][e] *
-                gradients_quad[(comp * dim + e) * nqp + q_point];
+                gradients_from_hessians_quad[(comp * dim + e) * nqp + q_point];
 
           // add off-diagonal part of J' * grad(u)
           for (unsigned int d = 0, count = dim; d < dim; ++d)
@@ -6047,7 +6065,8 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
               for (unsigned int f = 0; f < dim; ++f)
                 hessian_out[comp][d][e] +=
                   jac_grad[count][f] *
-                  gradients_quad[(comp * dim + f) * nqp + q_point];
+                  gradients_from_hessians_quad[(comp * dim + f) * nqp +
+                                               q_point];
 
           // take symmetric part
           for (unsigned int d = 0; d < dim; ++d)
@@ -6144,7 +6163,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
             for (unsigned int e = 0; e < dim; ++e)
               hessian_out[comp][d] +=
                 jac_grad[d][e] *
-                gradients_quad[(comp * dim + e) * nqp + q_point];
+                gradients_from_hessians_quad[(comp * dim + e) * nqp + q_point];
         }
     }
   return hessian_out;
@@ -6389,7 +6408,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
               hessian_in[comp][d][d] * factor;
         }
 
-      // off diagonal part (use symmetry)
+      // off diagonal part
       for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
         for (unsigned int e = 0; e < d; ++e, ++off_dia)
           {
@@ -6398,7 +6417,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
             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;
+                (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
           }
     }
   // cell with general Jacobian, but constant within the cell
@@ -6413,39 +6432,87 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
           VectorizedArrayType tmp[dim][dim];
           for (unsigned int i = 0; i < dim; ++i)
             for (unsigned int j = 0; j < dim; ++j)
-              tmp[i][j] = 0.;
+              {
+                tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
+                for (unsigned int k = 1; k < dim; ++k)
+                  tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
+              }
 
+          // 2. hessian_unit = J^T * tmp
+          VectorizedArrayType tmp2[dim][dim];
           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
+              {
+                tmp2[i][j] = jac[0][i] * tmp[0][j];
+                for (unsigned int k = 1; k < dim; ++k)
+                  tmp2[i][j] += jac[k][i] * tmp[k][j];
+              }
 
           // 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];
-            }
+            hessians_quad[(comp * hdim + d) * nqp + q_point] = tmp2[d][d] * JxW;
 
-          // 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];
-              }
+          // off diagonal part
+          for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
+            for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
+              hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
+                (tmp2[d][e] + tmp2[e][d]) * JxW;
         }
     }
   else
     {
-      // TODO: implement general path submit_hessian
-      AssertThrow(false, ExcNotImplemented());
+      const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[q_point];
+      const VectorizedArrayType                 JxW = this->J_value[q_point];
+      const auto &                              jac_grad =
+        this->mapping_data->jacobian_gradients
+          [1 - this->is_interior_face]
+          [this->mapping_data->data_index_offsets[this->cell] + 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] = hessian_in[comp][i][0] * jac[0][j];
+                for (unsigned int k = 1; k < dim; ++k)
+                  tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
+              }
+
+          // 2. hessian_unit = J^T * tmp
+          VectorizedArrayType tmp2[dim][dim];
+          for (unsigned int i = 0; i < dim; ++i)
+            for (unsigned int j = 0; j < dim; ++j)
+              {
+                tmp2[i][j] = jac[0][i] * tmp[0][j];
+                for (unsigned int k = 1; k < dim; ++k)
+                  tmp2[i][j] += jac[k][i] * tmp[k][j];
+              }
+
+          // diagonal part
+          for (unsigned int d = 0; d < dim; ++d)
+            hessians_quad[(comp * hdim + d) * nqp + q_point] = tmp2[d][d] * JxW;
+
+          // off diagonal part
+          for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
+            for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
+              hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
+                (tmp2[d][e] + tmp2[e][d]) * JxW;
+
+          // 3. gradient_unit = J' ** hessian_in
+          for (unsigned int d = 0; d < dim; ++d)
+            {
+              VectorizedArrayType sum = 0;
+              for (unsigned int e = 0; e < dim; ++e)
+                sum += hessian_in[comp][e][e] * jac_grad[e][d];
+              for (unsigned int e = 0, count = dim; e < dim; ++e)
+                for (unsigned int f = e + 1; f < dim; ++f, ++count)
+                  sum += (hessian_in[comp][e][f] + hessian_in[comp][f][e]) *
+                         jac_grad[count][d];
+              gradients_from_hessians_quad[(comp * dim + d) * nqp + q_point] =
+                sum * JxW;
+            }
+        }
     }
 }
 
@@ -8281,12 +8348,19 @@ FEEvaluation<dim,
              Number,
              VectorizedArrayType>::
   evaluate(const VectorizedArrayType *            values_array,
-           const EvaluationFlags::EvaluationFlags evaluation_flags)
+           const EvaluationFlags::EvaluationFlags evaluation_flag)
 {
+  const bool hessians_on_general_cells =
+    evaluation_flag & EvaluationFlags::hessians &&
+    (this->cell_type > internal::MatrixFreeFunctions::affine);
+  EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
+  if (hessians_on_general_cells)
+    evaluation_flag_actual |= EvaluationFlags::gradients;
+
   if (fe_degree > -1)
     SelectEvaluator<dim, fe_degree, n_q_points_1d, VectorizedArrayType>::
       evaluate(n_components,
-               evaluation_flags,
+               evaluation_flag_actual,
                *this->data,
                const_cast<VectorizedArrayType *>(values_array),
                this->values_quad,
@@ -8296,7 +8370,7 @@ FEEvaluation<dim,
   else
     internal::FEEvaluationFactory<dim, Number, VectorizedArrayType>::evaluate(
       n_components,
-      evaluation_flags,
+      evaluation_flag_actual,
       *this->data,
       const_cast<VectorizedArrayType *>(values_array),
       this->values_quad,
@@ -8304,12 +8378,16 @@ FEEvaluation<dim,
       this->hessians_quad,
       this->scratch_data);
 
+  if (hessians_on_general_cells)
+    for (unsigned int i = 0; i < n_components * dim * n_q_points; ++i)
+      this->gradients_from_hessians_quad[i] = this->gradients_quad[i];
+
 #  ifdef DEBUG
-  if (evaluation_flags & EvaluationFlags::values)
+  if (evaluation_flag_actual & EvaluationFlags::values)
     this->values_quad_initialized = true;
-  if (evaluation_flags & EvaluationFlags::gradients)
+  if (evaluation_flag_actual & EvaluationFlags::gradients)
     this->gradients_quad_initialized = true;
-  if (evaluation_flags & EvaluationFlags::hessians)
+  if (evaluation_flag_actual & EvaluationFlags::hessians)
     this->hessians_quad_initialized = true;
 #  endif
 }
@@ -8678,6 +8756,9 @@ FEEvaluation<dim,
   if (integration_flag & EvaluationFlags::gradients)
     Assert(this->gradients_quad_submitted == true,
            internal::ExcAccessToUninitializedField());
+  if (integration_flag & EvaluationFlags::hessians)
+    Assert(this->hessians_quad_submitted == true,
+           internal::ExcAccessToUninitializedField());
 #  endif
   Assert(this->matrix_info != nullptr ||
            this->mapped_geometry->is_initialized(),
@@ -8689,10 +8770,28 @@ FEEvaluation<dim,
     ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, and "
                "EvaluationFlags::hessians are supported."));
 
+  EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
+  if (integration_flag & EvaluationFlags::hessians &&
+      (this->cell_type > internal::MatrixFreeFunctions::affine))
+    {
+      unsigned int size = n_components * dim * n_q_points;
+      if (integration_flag & EvaluationFlags::gradients)
+        {
+          for (unsigned int i = 0; i < size; ++i)
+            this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
+        }
+      else
+        {
+          for (unsigned int i = 0; i < size; ++i)
+            this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
+          integration_flag_actual |= EvaluationFlags::gradients;
+        }
+    }
+
   if (fe_degree > -1)
     SelectEvaluator<dim, fe_degree, n_q_points_1d, VectorizedArrayType>::
       integrate(n_components,
-                integration_flag,
+                integration_flag_actual,
                 *this->data,
                 values_array,
                 this->values_quad,
@@ -8703,7 +8802,7 @@ FEEvaluation<dim,
   else
     internal::FEEvaluationFactory<dim, Number, VectorizedArrayType>::integrate(
       n_components,
-      integration_flag,
+      integration_flag_actual,
       *this->data,
       values_array,
       this->values_quad,
@@ -9112,6 +9211,13 @@ FEFaceEvaluation<dim,
       !(evaluation_flag & EvaluationFlags::hessians))
     return;
 
+  const bool hessians_on_general_cells =
+    evaluation_flag & EvaluationFlags::hessians &&
+    (this->cell_type > internal::MatrixFreeFunctions::affine);
+  EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
+  if (hessians_on_general_cells)
+    evaluation_flag_actual |= EvaluationFlags::gradients;
+
   if (this->dof_access_index ==
         internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
       this->is_interior_face == false)
@@ -9140,9 +9246,9 @@ FEFaceEvaluation<dim,
           this->begin_gradients(),
           this->begin_hessians(),
           this->scratch_data,
-          evaluation_flag & EvaluationFlags::values,
-          evaluation_flag & EvaluationFlags::gradients,
-          evaluation_flag & EvaluationFlags::hessians,
+          evaluation_flag_actual & EvaluationFlags::values,
+          evaluation_flag_actual & EvaluationFlags::gradients,
+          evaluation_flag_actual & EvaluationFlags::hessians,
           face_nos[0],
           this->subface_index,
           face_orientations[0],
@@ -9161,9 +9267,9 @@ FEFaceEvaluation<dim,
             this->begin_gradients(),
             this->begin_hessians(),
             this->scratch_data,
-            evaluation_flag & EvaluationFlags::values,
-            evaluation_flag & EvaluationFlags::gradients,
-            evaluation_flag & EvaluationFlags::hessians,
+            evaluation_flag_actual & EvaluationFlags::values,
+            evaluation_flag_actual & EvaluationFlags::gradients,
+            evaluation_flag_actual & EvaluationFlags::hessians,
             this->face_no,
             this->subface_index,
             this->face_orientation,
@@ -9177,21 +9283,28 @@ FEFaceEvaluation<dim,
                    this->begin_gradients(),
                    this->begin_hessians(),
                    this->scratch_data,
-                   evaluation_flag & EvaluationFlags::values,
-                   evaluation_flag & EvaluationFlags::gradients,
-                   evaluation_flag & EvaluationFlags::hessians,
+                   evaluation_flag_actual & EvaluationFlags::values,
+                   evaluation_flag_actual & EvaluationFlags::gradients,
+                   evaluation_flag_actual & EvaluationFlags::hessians,
                    this->face_no,
                    this->subface_index,
                    this->face_orientation,
                    this->descriptor->face_orientations);
     }
 
+  if (hessians_on_general_cells)
+    {
+      unsigned int size = n_components * dim * n_q_points;
+      for (unsigned int i = 0; i < size; ++i)
+        this->gradients_from_hessians_quad[i] = this->gradients_quad[i];
+    }
+
 #  ifdef DEBUG
-  if (evaluation_flag & EvaluationFlags::values)
+  if (evaluation_flag_actual & EvaluationFlags::values)
     this->values_quad_initialized = true;
-  if (evaluation_flag & EvaluationFlags::gradients)
+  if (evaluation_flag_actual & EvaluationFlags::gradients)
     this->gradients_quad_initialized = true;
-  if (evaluation_flag & EvaluationFlags::hessians)
+  if (evaluation_flag_actual & EvaluationFlags::hessians)
     this->hessians_quad_initialized = true;
 #  endif
 }
@@ -9300,6 +9413,24 @@ FEFaceEvaluation<dim,
       !(evaluation_flag & EvaluationFlags::hessians))
     return;
 
+  EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
+  if (evaluation_flag & EvaluationFlags::hessians &&
+      (this->cell_type > internal::MatrixFreeFunctions::affine))
+    {
+      unsigned int size = n_components * dim * n_q_points;
+      if (evaluation_flag & EvaluationFlags::gradients)
+        {
+          for (unsigned int i = 0; i < size; ++i)
+            this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
+        }
+      else
+        {
+          for (unsigned int i = 0; i < size; ++i)
+            this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
+          evaluation_flag_actual |= EvaluationFlags::gradients;
+        }
+    }
+
   if (fe_degree > -1)
     internal::FEFaceEvaluationImplIntegrateSelector<dim, VectorizedArrayType>::
       template run<fe_degree, n_q_points_1d>(
@@ -9310,9 +9441,9 @@ FEFaceEvaluation<dim,
         this->begin_gradients(),
         this->begin_hessians(),
         this->scratch_data,
-        evaluation_flag & EvaluationFlags::values,
-        evaluation_flag & EvaluationFlags::gradients,
-        evaluation_flag & EvaluationFlags::hessians,
+        evaluation_flag_actual & EvaluationFlags::values,
+        evaluation_flag_actual & EvaluationFlags::gradients,
+        evaluation_flag_actual & EvaluationFlags::hessians,
         this->face_no,
         this->subface_index,
         this->face_orientation,
@@ -9326,9 +9457,9 @@ FEFaceEvaluation<dim,
                 this->begin_gradients(),
                 this->begin_hessians(),
                 this->scratch_data,
-                evaluation_flag & EvaluationFlags::values,
-                evaluation_flag & EvaluationFlags::gradients,
-                evaluation_flag & EvaluationFlags::hessians,
+                evaluation_flag_actual & EvaluationFlags::values,
+                evaluation_flag_actual & EvaluationFlags::gradients,
+                evaluation_flag_actual & EvaluationFlags::hessians,
                 this->face_no,
                 this->subface_index,
                 this->face_orientation,
index 2e3e1ecf174554de06c0afbe33fd4c2d7dbc6ec9..2c3072008f27b35e8fe281a7423c8fcb7e22e23d 100644 (file)
@@ -358,6 +358,10 @@ namespace internal
              update_quadrature_points ?
            update_quadrature_points :
            update_default) |
+        ((update_flags_inner_faces | update_flags_boundary_faces) &
+             (update_jacobian_grads | update_hessians) ?
+           update_jacobian_grads :
+           update_default) |
         update_normal_vectors | update_JxW_values | update_jacobians;
       this->update_flags_inner_faces    = this->update_flags_boundary_faces;
       this->update_flags_faces_by_cells = update_flags_faces_by_cells;
@@ -428,7 +432,7 @@ namespace internal
 
                   const auto quad_face = quad[my_q][hpq].get_tensor_basis()[0];
                   face_data[my_q].descriptor[hpq * scale].initialize(
-                    quad_face, update_flags_boundary_faces);
+                    quad_face, this->update_flags_boundary_faces);
                   face_data[my_q].q_collection[hpq] =
                     dealii::hp::QCollection<dim - 1>(quad_face);
                   face_data_by_cells[my_q].descriptor[hpq * scale].initialize(
@@ -496,8 +500,6 @@ namespace internal
 
 
 
-    /* ------------------------- initialization of cells ------------------- */
-
     // Copy a vectorized array of one type to another type
     template <typename VectorizedArrayType1, typename VectorizedArrayType2>
     inline DEAL_II_ALWAYS_INLINE void
@@ -516,6 +518,76 @@ namespace internal
 
 
 
+    // For second derivatives on the real cell, we need the gradient of the
+    // inverse Jacobian J. This involves some calculus and is done
+    // vectorized. If L is the gradient of the jacobian on the unit cell,
+    // the gradient of the inverse is given by (multidimensional calculus) -
+    // J * (J * L) * J (the third J is because we need to transform the
+    // gradient L from the unit to the real cell, and then apply the inverse
+    // Jacobian). Compare this with 1D with j(x) = 1/k(phi(x)), where j =
+    // phi' is the inverse of the jacobian and k is the derivative of the
+    // jacobian on the unit cell. Then j' = phi' k'/k^2 = j k' j^2.
+    template <int dim, typename Number>
+    Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, Number>>
+    process_jacobian_gradient(const Tensor<2, dim, Number> &inv_jac_permut,
+                              const Tensor<2, dim, Number> &inv_jac,
+                              const Tensor<3, dim, Number> &jac_grad)
+    {
+      Number inv_jac_grad[dim][dim][dim];
+
+      // compute: inv_jac_grad = inv_jac_permut * grad_unit(jac)
+      for (unsigned int d = 0; d < dim; ++d)
+        for (unsigned int e = 0; e < dim; ++e)
+          for (unsigned int f = 0; f < dim; ++f)
+            {
+              inv_jac_grad[f][e][d] =
+                (inv_jac_permut[f][0] * jac_grad[d][e][0]);
+              for (unsigned int g = 1; g < dim; ++g)
+                inv_jac_grad[f][e][d] +=
+                  (inv_jac_permut[f][g] * jac_grad[d][e][g]);
+            }
+
+      // compute: transpose (-inv_jac_permut * inv_jac_grad[d] * inv_jac)
+      Number tmp[dim];
+      Number grad_jac_inv[dim][dim][dim];
+      for (unsigned int d = 0; d < dim; ++d)
+        for (unsigned int e = 0; e < dim; ++e)
+          {
+            for (unsigned int f = 0; f < dim; ++f)
+              {
+                tmp[f] = -inv_jac_grad[d][f][0] * inv_jac[0][e];
+                for (unsigned int g = 1; g < dim; ++g)
+                  tmp[f] -= inv_jac_grad[d][f][g] * inv_jac[g][e];
+              }
+
+            // needed for non-diagonal part of Jacobian grad
+            for (unsigned int f = 0; f < dim; ++f)
+              {
+                grad_jac_inv[f][d][e] = inv_jac_permut[f][0] * tmp[0];
+                for (unsigned int g = 1; g < dim; ++g)
+                  grad_jac_inv[f][d][e] += inv_jac_permut[f][g] * tmp[g];
+              }
+          }
+
+      Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, Number>> result;
+
+      // the diagonal part of Jacobian gradient comes first
+      for (unsigned int d = 0; d < dim; ++d)
+        for (unsigned int f = 0; f < dim; ++f)
+          result[d][f] = grad_jac_inv[d][d][f];
+
+      // then the upper-diagonal part
+      for (unsigned int d = 0, count = 0; d < dim; ++d)
+        for (unsigned int e = d + 1; e < dim; ++e, ++count)
+          for (unsigned int f = 0; f < dim; ++f)
+            result[dim + count][f] = grad_jac_inv[d][e][f];
+      return result;
+    }
+
+
+
+    /* ------------------------- initialization of cells ------------------- */
+
     // Namespace with implementation of extraction of values on cell
     // range
     namespace ExtractCellHelper
@@ -589,68 +661,7 @@ namespace internal
           }
       }
 
-      // For second derivatives on the real cell, we need the gradient of the
-      // inverse Jacobian J. This involves some calculus and is done
-      // vectorized. If L is the gradient of the jacobian on the unit cell,
-      // the gradient of the inverse is given by (multidimensional calculus) -
-      // J * (J * L) * J (the third J is because we need to transform the
-      // gradient L from the unit to the real cell, and then apply the inverse
-      // Jacobian). Compare this with 1D with j(x) = 1/k(phi(x)), where j =
-      // phi' is the inverse of the jacobian and k is the derivative of the
-      // jacobian on the unit cell. Then j' = phi' k'/k^2 = j k' j^2.
-      template <int dim, typename Number>
-      Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, Number>>
-      process_jacobian_gradient(const Tensor<2, dim, Number> &inv_jac,
-                                const Tensor<3, dim, Number> &jac_grad)
-      {
-        Number inv_jac_grad[dim][dim][dim];
-
-        // compute: inv_jac_grad = J*grad_unit(J^-1)
-        for (unsigned int d = 0; d < dim; ++d)
-          for (unsigned int e = 0; e < dim; ++e)
-            for (unsigned int f = 0; f < dim; ++f)
-              {
-                inv_jac_grad[f][e][d] = (inv_jac[f][0] * jac_grad[d][e][0]);
-                for (unsigned int g = 1; g < dim; ++g)
-                  inv_jac_grad[f][e][d] += (inv_jac[f][g] * jac_grad[d][e][g]);
-              }
-
-        // compute: transpose (-jac * jac_grad[d] * jac)
-        Number tmp[dim];
-        Number grad_jac_inv[dim][dim][dim];
-        for (unsigned int d = 0; d < dim; ++d)
-          for (unsigned int e = 0; e < dim; ++e)
-            {
-              for (unsigned int f = 0; f < dim; ++f)
-                {
-                  tmp[f] = Number();
-                  for (unsigned int g = 0; g < dim; ++g)
-                    tmp[f] -= inv_jac_grad[d][f][g] * inv_jac[g][e];
-                }
 
-              // needed for non-diagonal part of Jacobian grad
-              for (unsigned int f = 0; f < dim; ++f)
-                {
-                  grad_jac_inv[f][d][e] = inv_jac[f][0] * tmp[0];
-                  for (unsigned int g = 1; g < dim; ++g)
-                    grad_jac_inv[f][d][e] += inv_jac[f][g] * tmp[g];
-                }
-            }
-
-        Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, Number>> result;
-
-        // the diagonal part of Jacobian gradient comes first
-        for (unsigned int d = 0; d < dim; ++d)
-          for (unsigned int e = 0; e < dim; ++e)
-            result[d][e] = grad_jac_inv[d][d][e];
-
-        // then the upper-diagonal part
-        for (unsigned int d = 0, count = 0; d < dim; ++d)
-          for (unsigned int e = d + 1; e < dim; ++e, ++count)
-            for (unsigned int f = 0; f < dim; ++f)
-              result[dim + count][f] = grad_jac_inv[d][e][f];
-        return result;
-      }
 
       /**
        * Helper function called internally during the initialize function.
@@ -1049,7 +1060,9 @@ namespace internal
 
                       if (update_flags & update_jacobian_grads)
                         data.first[my_q].jacobian_gradients[0].push_back(
-                          process_jacobian_gradient(inv_jac, jacobian_grad));
+                          process_jacobian_gradient(inv_jac,
+                                                    inv_jac,
+                                                    jacobian_grad));
                     }
                 }
 
@@ -1504,7 +1517,9 @@ namespace internal
                                   cell_grad_grads[q + (d * hess_dim + c) *
                                                         n_q_points];
                             const auto inv_jac_grad =
-                              process_jacobian_gradient(inv_jac, jac_grad);
+                              process_jacobian_gradient(inv_jac,
+                                                        inv_jac,
+                                                        jac_grad);
                             for (unsigned int d = 0; d < hess_dim; ++d)
                               for (unsigned int e = 0; e < dim; ++e)
                                 store_vectorized_array(
@@ -2355,8 +2370,56 @@ namespace internal
               if (process_face[face] == false)
                 continue;
 
-              // go through the faces and fill the result
               const unsigned int offset = my_data.data_index_offsets[face];
+
+              const auto compute_jacobian_grad =
+                [&](unsigned int                     face_no,
+                    int                              is_exterior,
+                    unsigned int                     q,
+                    Tensor<2, dim, VectorizedDouble> inv_jac) {
+                  Tensor<2, dim, VectorizedDouble> inv_transp_jac_permut;
+                  for (unsigned int d = 0; d < dim; ++d)
+                    for (unsigned int e = 0; e < dim; ++e)
+                      {
+                        const unsigned int ee =
+                          ExtractFaceHelper::reorder_face_derivative_indices<
+                            dim>(face_no, e);
+                        inv_transp_jac_permut[d][e] = inv_jac[ee][d];
+                      }
+                  Tensor<2, dim, VectorizedDouble> jacobi;
+                  for (unsigned int e = 0; e < dim; ++e)
+                    for (unsigned int d = 0; d < dim; ++d)
+                      jacobi[d][e] = face_grads[(d * dim + e) * n_q_points + q];
+                  Tensor<2, dim, VectorizedDouble> inv_transp_jac =
+                    transpose(invert(jacobi));
+                  Tensor<3, dim, VectorizedDouble> jac_grad;
+                  for (unsigned int d = 0; d < dim; ++d)
+                    {
+                      for (unsigned int e = 0; e < dim; ++e)
+                        jac_grad[d][e][e] =
+                          face_grad_grads[q + (d * hess_dim + e) * n_q_points];
+                      for (unsigned int c = dim, e = 0; e < dim; ++e)
+                        for (unsigned int f = e + 1; f < dim; ++f, ++c)
+                          jac_grad[d][e][f] = jac_grad[d][f][e] =
+                            face_grad_grads[q +
+                                            (d * hess_dim + c) * n_q_points];
+                      const auto inv_jac_grad =
+                        process_jacobian_gradient(inv_transp_jac_permut,
+                                                  inv_transp_jac,
+                                                  jac_grad);
+                      for (unsigned int e = 0; e < dim; ++e)
+                        {
+                          for (unsigned int d = 0; d < hess_dim; ++d)
+                            store_vectorized_array(
+                              inv_jac_grad[d][e],
+                              vv,
+                              my_data.jacobian_gradients[is_exterior]
+                                                        [offset + q][d][e]);
+                        }
+                    }
+                };
+
+              // go through the faces and fill the result
               const unsigned int n_points_compute =
                 face_type[face] <= affine ? 1 : n_q_points;
               for (unsigned int q = 0; q < n_points_compute; ++q)
@@ -2385,8 +2448,7 @@ namespace internal
 
                   if (update_flags_faces & update_jacobian_grads)
                     {
-                      // TODO: implement jacobian grads for general path
-                      AssertThrow(false, ExcNotImplemented());
+                      compute_jacobian_grad(face_no, 0, q, inv_jac);
                     }
 
                   std::array<Tensor<1, dim, VectorizedDouble>, dim - 1>
@@ -2494,8 +2556,10 @@ namespace internal
 
                       if (update_flags_faces & update_jacobian_grads)
                         {
-                          // TODO: implement jacobian grads for general path
-                          AssertThrow(false, ExcNotImplemented());
+                          compute_jacobian_grad(faces[face].exterior_face_no,
+                                                1,
+                                                q,
+                                                inv_jac);
                         }
 
                       my_data.normals_times_jacobians[1][offset + q] =
diff --git a/tests/matrix_free/integrate_functions_hessians.cc b/tests/matrix_free/integrate_functions_hessians.cc
new file mode 100644 (file)
index 0000000..11363a5
--- /dev/null
@@ -0,0 +1,263 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// this function tests the correctness of the implementation of matrix free
+// operations in integrating functions, gradients and hessians on a hyperball
+// mesh with adaptive refinement.
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+
+template <int dim, int fe_degree, typename Number>
+class MatrixFreeTest
+{
+public:
+  using VectorType = std::vector<Vector<Number> *>;
+
+  MatrixFreeTest(const MatrixFree<dim, Number> &data_in)
+    : data(data_in)
+    , fe_val(data.get_dof_handler().get_fe(),
+             Quadrature<dim>(data.get_quadrature(0)),
+             update_values | update_gradients | update_hessians |
+               update_JxW_values){};
+
+  void
+  operator()(const MatrixFree<dim, Number> &              data,
+             VectorType &                                 dst,
+             const VectorType &                           src,
+             const std::pair<unsigned int, unsigned int> &cell_range) const;
+
+  void
+  test_functions(Vector<Number> &dst, Vector<Number> &dst_deal) const
+  {
+    dst      = 0;
+    dst_deal = 0;
+    VectorType dst_data(2);
+    dst_data[0] = &dst;
+    dst_data[1] = &dst_deal;
+    VectorType src_dummy;
+    data.cell_loop(&MatrixFreeTest<dim, fe_degree, Number>::operator(),
+                   this,
+                   dst_data,
+                   src_dummy);
+  };
+
+private:
+  const MatrixFree<dim, Number> &data;
+  mutable FEValues<dim>          fe_val;
+};
+
+
+
+template <int dim, int fe_degree, typename Number>
+void
+MatrixFreeTest<dim, fe_degree, Number>::operator()(
+  const MatrixFree<dim, Number> &data,
+  std::vector<Vector<Number> *> &dst,
+  const std::vector<Vector<Number> *> &,
+  const std::pair<unsigned int, unsigned int> &cell_range) const
+{
+  FEEvaluation<dim, fe_degree, fe_degree + 1, 1, Number> fe_eval(data);
+  const unsigned int                     n_q_points    = fe_eval.n_q_points;
+  const unsigned int                     dofs_per_cell = fe_eval.dofs_per_cell;
+  AlignedVector<VectorizedArray<Number>> values(n_q_points);
+  AlignedVector<VectorizedArray<Number>> gradients(dim * n_q_points);
+  AlignedVector<VectorizedArray<Number>> hessians(dim * dim * n_q_points);
+  std::vector<types::global_dof_index>   dof_indices(dofs_per_cell);
+  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+    {
+      fe_eval.reinit(cell);
+      // compare values with the ones the FEValues
+      // gives us. Those are seen as reference
+      for (unsigned int j = 0; j < data.n_components_filled(cell); ++j)
+        {
+          // generate random numbers at quadrature
+          // points and test them with basis functions
+          // and their gradients
+          for (unsigned int q = 0; q < n_q_points; ++q)
+            {
+              values[q][j] = random_value<double>();
+              for (unsigned int d = 0; d < dim; ++d)
+                gradients[q * dim + d][j] = -1. + 2. * (random_value<double>());
+              for (unsigned int d = 0; d < dim; ++d)
+                for (unsigned int e = 0; e < dim; ++e)
+                  hessians[q * dim * dim + d * dim + e][j] =
+                    -1. + 2. * (random_value<double>());
+            }
+          fe_val.reinit(data.get_cell_iterator(cell, j));
+          data.get_cell_iterator(cell, j)->get_dof_indices(dof_indices);
+
+          for (unsigned int i = 0; i < dofs_per_cell; ++i)
+            {
+              double sum = 0.;
+              for (unsigned int q = 0; q < n_q_points; ++q)
+                {
+                  sum +=
+                    values[q][j] * fe_val.shape_value(i, q) * fe_val.JxW(q);
+                  for (unsigned int d = 0; d < dim; ++d)
+                    sum += (gradients[q * dim + d][j] *
+                            fe_val.shape_grad(i, q)[d] * fe_val.JxW(q));
+                  for (unsigned int d = 0; d < dim; ++d)
+                    for (unsigned int e = 0; e < dim; ++e)
+                      sum += (hessians[q * dim * dim + d * dim + e][j] *
+                              fe_val.shape_hessian(i, q)[d][e] * fe_val.JxW(q));
+                }
+              (*dst[1])(dof_indices[i]) += sum;
+            }
+        }
+      for (unsigned int q = 0; q < n_q_points; ++q)
+        {
+          fe_eval.submit_value(values[q], q);
+          Tensor<1, dim, VectorizedArray<Number>> submit;
+          for (unsigned int d = 0; d < dim; ++d)
+            submit[d] = gradients[q * dim + d];
+          fe_eval.submit_gradient(submit, q);
+          Tensor<2, dim, VectorizedArray<Number>> submit2;
+          for (unsigned int d = 0; d < dim; ++d)
+            for (unsigned int e = 0; e < dim; ++e)
+              submit2[d][e] = hessians[q * dim * dim + d * dim + e];
+          fe_eval.submit_hessian(submit2, q);
+        }
+      fe_eval.integrate(EvaluationFlags::values | EvaluationFlags::gradients |
+                        EvaluationFlags::hessians);
+      fe_eval.distribute_local_to_global(*dst[0]);
+    }
+}
+
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  using number = double;
+  const SphericalManifold<dim> manifold;
+  Triangulation<dim>           tria;
+  GridGenerator::hyper_ball(tria);
+  typename Triangulation<dim>::active_cell_iterator cell = tria.begin_active(),
+                                                    endc = tria.end();
+  for (; cell != endc; ++cell)
+    for (const unsigned int f : GeometryInfo<dim>::face_indices())
+      if (cell->at_boundary(f))
+        cell->face(f)->set_all_manifold_ids(0);
+  tria.set_manifold(0, manifold);
+
+  cell = tria.begin_active();
+  for (; cell != endc; ++cell)
+    if (cell->center().norm() < 1e-8)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  cell = tria.begin_active();
+  for (; cell != endc; ++cell)
+    if (cell->center().norm() < 0.2)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  if (dim < 3 || fe_degree < 2)
+    tria.refine_global(1);
+  tria.begin(tria.n_levels() - 1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  cell = tria.begin_active();
+  for (unsigned int i = 0; i < 7 - 2 * dim; ++i)
+    {
+      cell                 = tria.begin_active();
+      unsigned int counter = 0;
+      for (; cell != endc; ++cell, ++counter)
+        if (counter % (7 - i) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+    }
+
+  FE_Q<dim>       fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+  deallog << "Testing " << fe.get_name() << std::endl;
+  // std::cout << "Number of cells: " << tria.n_active_cells() << std::endl;
+  // std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
+
+  AffineConstraints<double> constraints;
+  DoFTools::make_hanging_node_constraints(dof, constraints);
+  constraints.close();
+
+  MatrixFree<dim, number> mf_data;
+  {
+    typename MatrixFree<dim, number>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+    data.mapping_update_flags |= update_hessians;
+    const QGauss<1> quad(fe_degree + 1);
+    mf_data.reinit(dof, constraints, quad, data);
+  }
+
+  MatrixFreeTest<dim, fe_degree, number> mf(mf_data);
+  Vector<number>                         solution(dof.n_dofs());
+  Vector<number>                         solution_dist(dof.n_dofs());
+
+  mf.test_functions(solution_dist, solution);
+
+  constraints.condense(solution);
+
+  Vector<number> compare(solution_dist);
+  compare -= solution;
+  const double diff_norm = compare.linfty_norm();
+
+  deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(3);
+
+  {
+    deallog.push("2d");
+    test<2, 1>();
+    test<2, 2>();
+    test<2, 3>();
+    test<2, 4>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3, 1>();
+    test<3, 2>();
+    test<3, 3>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/integrate_functions_hessians.output b/tests/matrix_free/integrate_functions_hessians.output
new file mode 100644 (file)
index 0000000..3829c8f
--- /dev/null
@@ -0,0 +1,22 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(2)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(4)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(3)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
diff --git a/tests/matrix_free/matrix_vector_hessians_general_cells.cc b/tests/matrix_free/matrix_vector_hessians_general_cells.cc
new file mode 100644 (file)
index 0000000..d9171bc
--- /dev/null
@@ -0,0 +1,310 @@
+// ---------------------------------------------------------------------
+//
+// 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 FEEvaluation and FEValues.
+// This is checked for different combinations of EvaluationFlags, FE types and
+// polynomial degrees.
+
+template <int dim>
+void
+test_hessians(const dealii::FE_Poly<dim> &                   fe,
+              const dealii::Quadrature<dim> &                quad,
+              const dealii::EvaluationFlags::EvaluationFlags evaluation_flags)
+{
+  using namespace dealii;
+  using VectorizedArrayType = VectorizedArray<double>;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_ball(tria);
+  tria.refine_global(1);
+
+  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();
+
+  // FEEvaluation
+  typename MatrixFree<dim, double, VectorizedArrayType>::AdditionalData
+    additional_data;
+  additional_data.mapping_update_flags =
+    update_values | update_gradients | update_hessians;
+
+  MatrixFree<dim, double, VectorizedArrayType> matrix_free;
+  matrix_free.reinit(mapping, dof_handler, constraints, quad, 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
+  FEValues<dim>                        fe_values(mapping,
+                          fe,
+                          quad,
+                          update_values | update_gradients | update_hessians |
+                            update_JxW_values);
+  Vector<double>                       solution_values_local(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(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) {
+      FEEvaluation<dim, -1, 0, 1, double, VectorizedArrayType> fe_eval(
+        matrix_free);
+      for (unsigned int cell = range.first; cell < range.second; ++cell)
+        {
+          // FEEvaluation
+          fe_eval.reinit(cell);
+          fe_eval.gather_evaluate(src, evaluation_flags);
+          for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
+            {
+              if (evaluation_flags & EvaluationFlags::hessians)
+                fe_eval.submit_hessian(fe_eval.get_hessian(q), q);
+              if (evaluation_flags & EvaluationFlags::gradients)
+                fe_eval.submit_gradient(fe_eval.get_gradient(q), q);
+              if (evaluation_flags & EvaluationFlags::values)
+                fe_eval.submit_value(fe_eval.get_value(q), q);
+            }
+          fe_eval.integrate(evaluation_flags);
+          fe_eval.distribute_local_to_global(dst);
+
+          // FEValues
+          for (unsigned int v = 0;
+               v < matrix_free.n_active_entries_per_cell_batch(cell);
+               ++v)
+            {
+              const auto cell_iterator = matrix_free.get_cell_iterator(cell, v);
+
+              fe_values.reinit(cell_iterator);
+
+              cell_iterator->get_dof_indices(dof_indices);
+              constraints.get_dof_values(src,
+                                         dof_indices.begin(),
+                                         solution_values_local.begin(),
+                                         solution_values_local.end());
+
+              std::vector<Tensor<2, dim>> hessians_function(quad.size());
+              fe_values.get_function_hessians(src, hessians_function);
+
+              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_local(i) *
+                                    fe_values.shape_hessian(i, q);
+                      if (evaluation_flags & EvaluationFlags::gradients)
+                        gradients +=
+                          solution_values_local(i) * fe_values.shape_grad(i, q);
+                      if (evaluation_flags & EvaluationFlags::values)
+                        values += solution_values_local(i) *
+                                  fe_values.shape_value(i, q);
+                    }
+                  solution_hessians[q]  = hessians * fe_values.JxW(q);
+                  solution_gradients[q] = gradients * fe_values.JxW(q);
+                  solution_values[q]    = values * fe_values.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_values.shape_hessian(i, q));
+                      if (evaluation_flags & EvaluationFlags::gradients)
+                        sum_gradients +=
+                          solution_gradients[q] * fe_values.shape_grad(i, q);
+                      if (evaluation_flags & EvaluationFlags::values)
+                        sum_values +=
+                          solution_values[q] * fe_values.shape_value(i, q);
+                    }
+                  solution_values_local(i) =
+                    sum_hessians + sum_gradients + sum_values;
+                }
+              constraints.distribute_local_to_global(solution_values_local,
+                                                     dof_indices,
+                                                     dst2);
+            }
+        }
+    },
+    [&](
+      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) {
+      (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;
+  }
+}
+
+void
+test_qiterated(dealii::EvaluationFlags::EvaluationFlags evaluation_flags)
+{
+  using namespace dealii;
+  deallog << "test_qiterated" << std::endl;
+  for (unsigned int i = 1; i < 4; ++i)
+    {
+      QIterated<2> quad_2(QGauss<1>(i + 1), 2);
+      QIterated<3> quad_3(QGauss<1>(i + 1), 2);
+      test_hessians<2>(FE_Q<2>(i), quad_2, evaluation_flags);
+      test_hessians<3>(FE_Q<3>(i), quad_3, evaluation_flags);
+      test_hessians<2>(FE_DGQ<2>(i), quad_2, evaluation_flags);
+      test_hessians<3>(FE_DGQ<3>(i), quad_3, evaluation_flags);
+    }
+}
+
+void
+test_qgauss(dealii::EvaluationFlags::EvaluationFlags evaluation_flags)
+{
+  using namespace dealii;
+  deallog << "test_qgauss" << std::endl;
+  for (unsigned int i = 1; i < 4; ++i)
+    {
+      QGauss<2> quad_2(i + 1);
+      QGauss<3> quad_3(i + 1);
+      test_hessians<2>(FE_Q<2>(i), quad_2, evaluation_flags);
+      test_hessians<3>(FE_Q<3>(i), quad_3, evaluation_flags);
+      test_hessians<2>(FE_DGQ<2>(i), quad_2, evaluation_flags);
+      test_hessians<3>(FE_DGQ<3>(i), quad_3, evaluation_flags);
+    }
+}
+
+int
+main(int argc, char **argv)
+{
+  using namespace dealii;
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+
+  initlog();
+  deallog << std::setprecision(10);
+
+  {
+    EvaluationFlags::EvaluationFlags evaluation_flags =
+      EvaluationFlags::hessians;
+    deallog << "test_hessians_only" << std::endl;
+    test_qgauss(evaluation_flags);
+    test_qiterated(evaluation_flags);
+  }
+
+  {
+    EvaluationFlags::EvaluationFlags evaluation_flags =
+      EvaluationFlags::values | EvaluationFlags::hessians;
+    deallog << "test_hessians_with_values" << std::endl;
+    test_qgauss(evaluation_flags);
+    test_qiterated(evaluation_flags);
+  }
+
+  {
+    EvaluationFlags::EvaluationFlags evaluation_flags =
+      EvaluationFlags::gradients | EvaluationFlags::hessians;
+    deallog << "test_hessians_with_gradients" << std::endl;
+    test_qgauss(evaluation_flags);
+    test_qiterated(evaluation_flags);
+  }
+
+  {
+    EvaluationFlags::EvaluationFlags evaluation_flags =
+      EvaluationFlags::values | EvaluationFlags::gradients |
+      EvaluationFlags::hessians;
+    deallog << "test_hessians_with_gradients_and_values" << std::endl;
+    test_qgauss(evaluation_flags);
+    test_qiterated(evaluation_flags);
+  }
+}
diff --git a/tests/matrix_free/matrix_vector_hessians_general_cells.output b/tests/matrix_free/matrix_vector_hessians_general_cells.output
new file mode 100644 (file)
index 0000000..412bae7
--- /dev/null
@@ -0,0 +1,493 @@
+
+DEAL::test_hessians_only
+DEAL::test_qgauss
+DEAL::Working with FE_Q<2>(1) and 25 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 45.22219818 -7.037353151 0.000000000 11.88687292 -61.49576085 37.57202054 -35.68314307 
+DEAL::dst FEV: 0.000000000 0.000000000 45.22219818 -7.037353151 0.000000000 11.88687292 -61.49576085 37.57202054 -35.68314307 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 79 dofs
+DEAL::dst FEE: 33.10696022 -22.71237546 -68.25761235 64.86438841 -52.49601275 16.86513406 30.83643936 -7.704940388 -34.30613074 
+DEAL::dst FEV: 33.10696022 -22.71237546 -68.25761235 64.86438841 -52.49601275 16.86513406 30.83643936 -7.704940388 -34.30613074 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 80 dofs
+DEAL::dst FEE: -10.10812224 7.655480511 14.29504356 -11.84240183 -2.477974473 3.271861099 3.833223715 -4.627110340 8.081683648 
+DEAL::dst FEV: -10.10812224 7.655480511 14.29504356 -11.84240183 -2.477974473 3.271861099 3.833223715 -4.627110340 8.081683648 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 448 dofs
+DEAL::dst FEE: 6.235165197 -7.858320343 0.1523344538 1.470820692 -4.214851267 5.838006413 -2.172648384 0.5494932382 3.923280641 
+DEAL::dst FEV: 6.235165197 -7.858320343 0.1523344538 1.470820692 -4.214851267 5.838006413 -2.172648384 0.5494932382 3.923280641 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 89 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 -207.3090085 -399.8368570 365.2209748 349.7535808 0.000000000 438.9124523 -467.2428521 
+DEAL::dst FEV: 0.000000000 0.000000000 -207.3090085 -399.8368570 365.2209748 349.7535808 0.000000000 438.9124523 -467.2428521 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 517 dofs
+DEAL::dst FEE: -89.44696014 33.15396135 43.60705670 122.9701105 202.3600767 30.62778421 36.41906605 26.08882904 -14.27490495 
+DEAL::dst FEV: -89.44696014 33.15396135 43.60705670 122.9701105 202.3600767 30.62778421 36.41906605 26.08882904 -14.27490495 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 180 dofs
+DEAL::dst FEE: 66.71712651 251.4306103 -12.27594710 -170.3467774 -483.4850601 35.48250021 117.8761478 214.8381899 -20.23679013 
+DEAL::dst FEV: 66.71712651 251.4306103 -12.27594710 -170.3467774 -483.4850601 35.48250021 117.8761478 214.8381899 -20.23679013 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 1512 dofs
+DEAL::dst FEE: -9.588563680 19.11134992 -2.739954593 46.74620058 -92.01359251 28.89086228 -5.744434228 -16.14481508 19.53643851 
+DEAL::dst FEV: -9.588563680 19.11134992 -2.739954593 46.74620058 -92.01359251 28.89086228 -5.744434228 -16.14481508 19.53643851 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 193 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 115.2792512 1514.207777 129.6997082 -385.2037949 562.9695734 -1541.447662 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 115.2792512 1514.207777 129.6997082 -385.2037949 562.9695734 -1541.447662 0.000000000 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 1651 dofs
+DEAL::dst FEE: 182.6670278 -173.2496744 42.37257073 233.4142098 294.1859475 286.3887194 -175.4075400 -146.4736019 -864.9517464 
+DEAL::dst FEV: 182.6670278 -173.2496744 42.37257073 233.4142098 294.1859475 286.3887194 -175.4075400 -146.4736019 -864.9517464 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 320 dofs
+DEAL::dst FEE: -19.97069668 -2263.033643 1394.341009 -180.1207746 205.5919186 4843.256322 -3031.471695 363.4727737 -584.4861111 
+DEAL::dst FEV: -19.97069668 -2263.033643 1394.341009 -180.1207746 205.5919186 4843.256322 -3031.471695 363.4727737 -584.4861111 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 3584 dofs
+DEAL::dst FEE: 5.120375759 -102.7942743 35.13091124 26.35855211 -3.524129155 270.9555244 -269.4815327 0.9612496484 122.2229695 
+DEAL::dst FEV: 5.120375759 -102.7942743 35.13091124 26.35855211 -3.524129155 270.9555244 -269.4815327 0.9612496484 122.2229695 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::test_qiterated
+DEAL::Working with FE_Q<2>(1) and 25 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 -78.31750939 8.028509910 0.000000000 31.94914913 112.2418132 -8.404529286 -57.17545763 
+DEAL::dst FEV: 0.000000000 0.000000000 -78.31750939 8.028509910 0.000000000 31.94914913 112.2418132 -8.404529286 -57.17545763 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 79 dofs
+DEAL::dst FEE: -28.62894665 12.83395912 -3.022244288 11.75298450 54.81383740 -40.08721361 -25.26188941 -8.528413431 28.07903454 
+DEAL::dst FEV: -28.62894665 12.83395912 -3.022244288 11.75298450 54.81383740 -40.08721361 -25.26188941 -8.528413431 28.07903454 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 80 dofs
+DEAL::dst FEE: 8.563495216 -6.485642852 -12.11061108 10.03275871 -0.4091739206 0.5402639324 0.6329585687 -0.7640485805 -0.7596303223 
+DEAL::dst FEV: 8.563495216 -6.485642852 -12.11061108 10.03275871 -0.4091739206 0.5402639324 0.6329585687 -0.7640485805 -0.7596303223 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 448 dofs
+DEAL::dst FEE: -0.04925434952 4.129829593 5.129818715 -9.210393958 -5.559736708 1.479161465 0.4791723432 3.601402900 -0.9781090982 
+DEAL::dst FEV: -0.04925434952 4.129829593 5.129818715 -9.210393958 -5.559736708 1.479161465 0.4791723432 3.601402900 -0.9781090982 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 89 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 -360.3224836 -234.0666285 579.7961020 132.2256446 0.000000000 313.0605472 -3.884425931 
+DEAL::dst FEV: 0.000000000 0.000000000 -360.3224836 -234.0666285 579.7961020 132.2256446 0.000000000 313.0605472 -3.884425931 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 517 dofs
+DEAL::dst FEE: 139.4934997 117.6177543 38.83313679 49.94720211 -56.48270646 -7.179565374 38.65311026 -21.47267519 -149.4955229 
+DEAL::dst FEV: 139.4934997 117.6177543 38.83313679 49.94720211 -56.48270646 -7.179565374 38.65311026 -21.47267519 -149.4955229 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 180 dofs
+DEAL::dst FEE: 26.94269068 24.79455169 -50.15544477 -9.626629392 -173.3905474 175.5729219 -18.80078412 158.8729847 -134.2097432 
+DEAL::dst FEV: 26.94269068 24.79455169 -50.15544477 -9.626629392 -173.3905474 175.5729219 -18.80078412 158.8729847 -134.2097432 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 1512 dofs
+DEAL::dst FEE: -22.35910006 48.84059752 -5.609505384 36.26733584 -116.6366681 34.78761915 -11.86986656 34.27145167 -13.50944896 
+DEAL::dst FEV: -22.35910006 48.84059752 -5.609505384 36.26733584 -116.6366681 34.78761915 -11.86986656 34.27145167 -13.50944896 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 193 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 2281.203835 53.14339166 799.5213508 -1660.550741 -61.70284234 218.9584292 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 2281.203835 53.14339166 799.5213508 -1660.550741 -61.70284234 218.9584292 0.000000000 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 1651 dofs
+DEAL::dst FEE: 411.0414196 19.69629401 -45.09261305 -160.0937108 157.9595318 -63.70495743 -285.2088913 -88.92577089 854.5629205 
+DEAL::dst FEV: 411.0414196 19.69629401 -45.09261305 -160.0937108 157.9595318 -63.70495743 -285.2088913 -88.92577089 854.5629205 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 320 dofs
+DEAL::dst FEE: -255.5427784 812.0180402 956.7833695 -196.0315102 455.8205352 -1387.185608 -2097.876988 411.3732530 -14.36300381 
+DEAL::dst FEV: -255.5427784 812.0180402 956.7833695 -196.0315102 455.8205352 -1387.185608 -2097.876988 411.3732530 -14.36300381 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 3584 dofs
+DEAL::dst FEE: 17.32821876 -95.95926190 145.4413654 -15.35003562 -43.17407779 208.8850676 -131.1611595 -40.83909124 3.768804606 
+DEAL::dst FEV: 17.32821876 -95.95926190 145.4413654 -15.35003562 -43.17407779 208.8850676 -131.1611595 -40.83909124 3.768804606 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::test_hessians_with_values
+DEAL::test_qgauss
+DEAL::Working with FE_Q<2>(1) and 25 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 -16.58940782 -4.782648286 0.000000000 -19.53954235 59.54922800 -66.05480814 64.73504286 
+DEAL::dst FEV: 0.000000000 0.000000000 -16.58940782 -4.782648286 0.000000000 -19.53954235 59.54922800 -66.05480814 64.73504286 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 79 dofs
+DEAL::dst FEE: 32.15282885 21.33825418 -46.88409574 -5.642916990 4.268878289 21.02923008 -6.835454147 1.347064398 -35.37416089 
+DEAL::dst FEV: 32.15282885 21.33825418 -46.88409574 -5.642916990 4.268878289 21.02923008 -6.835454147 1.347064398 -35.37416089 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 80 dofs
+DEAL::dst FEE: 10.79722351 -8.139267187 -15.21835329 12.64734630 -1.429251592 1.946800889 2.271259133 -2.693814615 18.10526726 
+DEAL::dst FEV: 10.79722351 -8.139267187 -15.21835329 12.64734630 -1.429251592 1.946800889 2.271259133 -2.693814615 18.10526726 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 448 dofs
+DEAL::dst FEE: -3.738241635 8.395409955 0.6254841474 -5.280834462 -0.9752822717 -3.679925410 4.089943235 0.5673075883 0.3087058789 
+DEAL::dst FEV: -3.738241635 8.395409955 0.6254841474 -5.280834462 -0.9752822717 -3.679925410 4.089943235 0.5673075883 0.3087058789 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 89 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 380.5559204 -105.0922063 -145.8738327 65.76715372 0.000000000 -153.4360487 450.3778700 
+DEAL::dst FEV: 0.000000000 0.000000000 380.5559204 -105.0922063 -145.8738327 65.76715372 0.000000000 -153.4360487 450.3778700 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 517 dofs
+DEAL::dst FEE: 2.937489230 132.6971713 50.52850669 -15.20168086 -139.5233581 6.756647264 39.08008100 26.94906549 304.6782079 
+DEAL::dst FEV: 2.937489230 132.6971713 50.52850669 -15.20168086 -139.5233581 6.756647264 39.08008100 26.94906549 304.6782079 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 180 dofs
+DEAL::dst FEE: -81.08771912 334.1346990 13.63891270 221.2008497 -725.0794917 -21.98122432 -146.3265523 390.8131029 14.77884105 
+DEAL::dst FEV: -81.08771912 334.1346990 13.63891270 221.2008497 -725.0794917 -21.98122432 -146.3265523 390.8131029 14.77884105 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 1512 dofs
+DEAL::dst FEE: -4.494825168 -37.24295230 3.658353095 46.00761389 17.77100136 -4.702653614 -19.70483155 37.19542383 -11.33928896 
+DEAL::dst FEV: -4.494825168 -37.24295230 3.658353095 46.00761389 17.77100136 -4.702653614 -19.70483155 37.19542383 -11.33928896 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 193 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 1015.004669 418.3204606 890.1073227 -1057.084816 752.3299712 -759.4644896 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 1015.004669 418.3204606 890.1073227 -1057.084816 752.3299712 -759.4644896 0.000000000 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 1651 dofs
+DEAL::dst FEE: -161.7750270 -180.4296561 -16.05309565 -18.79545541 34.84814822 -168.2536254 -18.17675410 -29.88875304 618.8759516 
+DEAL::dst FEV: -161.7750270 -180.4296561 -16.05309565 -18.79545541 34.84814822 -168.2536254 -18.17675410 -29.88875304 618.8759516 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 320 dofs
+DEAL::dst FEE: -129.9833408 232.5435605 -1068.363295 159.0051914 261.6113130 -257.6833985 2175.493847 -287.8736806 -200.4477312 
+DEAL::dst FEV: -129.9833408 232.5435605 -1068.363295 159.0051914 261.6113130 -257.6833985 2175.493847 -287.8736806 -200.4477312 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 3584 dofs
+DEAL::dst FEE: 8.573820572 -70.82155289 -103.2850267 17.61450370 -51.13396496 124.0614029 312.1554282 -84.46207607 67.21754006 
+DEAL::dst FEV: 8.573820572 -70.82155289 -103.2850267 17.61450370 -51.13396496 124.0614029 312.1554282 -84.46207607 67.21754006 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::test_qiterated
+DEAL::Working with FE_Q<2>(1) and 25 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 18.53584598 -36.66169783 0.000000000 84.22806902 -40.72609295 109.2107485 -146.1168399 
+DEAL::dst FEV: 0.000000000 0.000000000 18.53584598 -36.66169783 0.000000000 84.22806902 -40.72609295 109.2107485 -146.1168399 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 79 dofs
+DEAL::dst FEE: 75.52280024 -48.98798905 -17.99142928 17.81212936 -39.28533289 -43.32180002 2.565804398 42.49911425 43.87285076 
+DEAL::dst FEV: 75.52280024 -48.98798905 -17.99142928 17.81212936 -39.28533289 -43.32180002 2.565804398 42.49911425 43.87285076 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 80 dofs
+DEAL::dst FEE: 3.331492342 -2.476604953 -4.673138931 3.912736086 -1.153920507 1.551152426 1.809397922 -2.168905382 22.77362202 
+DEAL::dst FEV: 3.331492342 -2.476604953 -4.673138931 3.912736086 -1.153920507 1.551152426 1.809397922 -2.168905382 22.77362202 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 448 dofs
+DEAL::dst FEE: -1.084466294 -1.172066883 6.039579603 -3.780049703 -2.346661305 4.606206981 -2.605353929 0.3489988567 -0.4208673164 
+DEAL::dst FEV: -1.084466294 -1.172066883 6.039579603 -3.780049703 -2.346661305 4.606206981 -2.605353929 0.3489988567 -0.4208673164 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 89 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 -309.5249641 -65.11577579 378.5689837 81.96287877 0.000000000 -18.49549227 -95.28776221 
+DEAL::dst FEV: 0.000000000 0.000000000 -309.5249641 -65.11577579 378.5689837 81.96287877 0.000000000 -18.49549227 -95.28776221 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 517 dofs
+DEAL::dst FEE: 92.07855370 37.90950393 32.56690227 -101.3425164 -77.46971424 73.09994955 -38.69427253 1.202702857 -86.04224867 
+DEAL::dst FEV: 92.07855370 37.90950393 32.56690227 -101.3425164 -77.46971424 73.09994955 -38.69427253 1.202702857 -86.04224867 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 180 dofs
+DEAL::dst FEE: -47.82882231 -102.6420815 76.20983134 145.0526142 166.6021274 -155.7349457 -123.9589010 -22.78484976 65.20942197 
+DEAL::dst FEV: -47.82882231 -102.6420815 76.20983134 145.0526142 166.6021274 -155.7349457 -123.9589010 -22.78484976 65.20942197 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 1512 dofs
+DEAL::dst FEE: -15.97522003 7.821566253 17.51495995 45.64355144 -65.09372864 -6.523935492 3.012415787 -1.810943197 9.600318192 
+DEAL::dst FEV: -15.97522003 7.821566253 17.51495995 45.64355144 -65.09372864 -6.523935492 3.012415787 -1.810943197 9.600318192 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 193 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 -209.0416960 -346.9818449 564.0182222 -1099.889775 575.2820845 -212.5628742 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 -209.0416960 -346.9818449 564.0182222 -1099.889775 575.2820845 -212.5628742 0.000000000 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 1651 dofs
+DEAL::dst FEE: 456.2392295 -118.0476821 229.4245288 -73.52012482 282.7625847 70.52133886 -142.1717202 163.4847811 271.8552957 
+DEAL::dst FEV: 456.2392295 -118.0476821 229.4245288 -73.52012482 282.7625847 70.52133886 -142.1717202 163.4847811 271.8552957 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 320 dofs
+DEAL::dst FEE: 243.7688278 -2179.037732 70.00926913 93.64438209 -357.9827311 4397.708618 149.2573164 -206.5012884 -35.52626013 
+DEAL::dst FEV: 243.7688278 -2179.037732 70.00926913 93.64438209 -357.9827311 4397.708618 149.2573164 -206.5012884 -35.52626013 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 3584 dofs
+DEAL::dst FEE: 42.60953463 1.905386931 4.105747246 -32.68370149 -66.55720248 -119.7933249 135.5493554 -15.93352778 22.40456232 
+DEAL::dst FEV: 42.60953463 1.905386931 4.105747246 -32.68370149 -66.55720248 -119.7933249 135.5493554 -15.93352778 22.40456232 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::test_hessians_with_gradients
+DEAL::test_qgauss
+DEAL::Working with FE_Q<2>(1) and 25 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 7.311269713 -5.156295237 0.000000000 -19.59465764 65.07757006 -96.82039038 98.14454960 
+DEAL::dst FEV: 0.000000000 0.000000000 7.311269713 -5.156295237 0.000000000 -19.59465764 65.07757006 -96.82039038 98.14454960 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 79 dofs
+DEAL::dst FEE: -1.505294846 -24.57206121 52.11921332 39.11055377 -9.219780077 13.43074516 -61.01667614 9.190914949 -12.82165271 
+DEAL::dst FEV: -1.505294846 -24.57206121 52.11921332 39.11055377 -9.219780077 13.43074516 -61.01667614 9.190914949 -12.82165271 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 80 dofs
+DEAL::dst FEE: -6.317046236 4.302693679 8.925812636 -6.911460079 -0.3150298140 0.1815725597 0.4627341489 -0.3292768946 -1.150968081 
+DEAL::dst FEV: -6.317046236 4.302693679 8.925812636 -6.911460079 -0.3150298140 0.1815725597 0.4627341489 -0.3292768946 -1.150968081 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 448 dofs
+DEAL::dst FEE: 0.1232404125 0.01022319883 1.416411963 -1.568321407 2.617554109 -2.833395117 -4.137893678 4.372180519 -1.115248686 
+DEAL::dst FEV: 0.1232404125 0.01022319883 1.416411963 -1.568321407 2.617554109 -2.833395117 -4.137893678 4.372180519 -1.115248686 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 89 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 -205.1606811 -109.2816798 169.2814957 25.66427265 0.000000000 -115.4425259 402.3574467 
+DEAL::dst FEV: 0.000000000 0.000000000 -205.1606811 -109.2816798 169.2814957 25.66427265 0.000000000 -115.4425259 402.3574467 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 517 dofs
+DEAL::dst FEE: -229.0756197 -169.2001752 0.4466035661 -74.58939161 -142.9520136 -29.63421128 68.81675020 57.23240132 -54.66379827 
+DEAL::dst FEV: -229.0756197 -169.2001752 0.4466035661 -74.58939161 -142.9520136 -29.63421128 68.81675020 57.23240132 -54.66379827 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 180 dofs
+DEAL::dst FEE: 56.25534771 319.1648241 -91.58165669 -170.0871302 -603.7836113 196.7016090 151.6182738 225.7414787 -84.02913509 
+DEAL::dst FEV: 56.25534771 319.1648241 -91.58165669 -170.0871302 -603.7836113 196.7016090 151.6182738 225.7414787 -84.02913509 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 1512 dofs
+DEAL::dst FEE: -15.84166619 0.2353443988 -2.115338897 23.85811429 49.84048831 -7.585819571 -14.10981493 -20.17993682 15.77644128 
+DEAL::dst FEV: -15.84166619 0.2353443988 -2.115338897 23.85811429 49.84048831 -7.585819571 -14.10981493 -20.17993682 15.77644128 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 193 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 257.9275170 938.3914844 1444.117736 -1047.901673 1340.863608 -1542.857351 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 257.9275170 938.3914844 1444.117736 -1047.901673 1340.863608 -1542.857351 0.000000000 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 1651 dofs
+DEAL::dst FEE: -272.6279247 29.31795285 -149.4439137 17.98294068 82.28977435 -333.7417511 170.4438834 50.90367540 571.9652347 
+DEAL::dst FEV: -272.6279247 29.31795285 -149.4439137 17.98294068 82.28977435 -333.7417511 170.4438834 50.90367540 571.9652347 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 320 dofs
+DEAL::dst FEE: 167.8239843 937.4290293 373.7473509 -19.98309967 -453.9251491 -1886.074462 -704.8832500 -71.02701366 558.1818119 
+DEAL::dst FEV: 167.8239843 937.4290293 373.7473509 -19.98309967 -453.9251491 -1886.074462 -704.8832500 -71.02701366 558.1818119 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 3584 dofs
+DEAL::dst FEE: -37.20119468 142.7420109 -3.007994596 1.347129280 58.28259413 -229.2925423 20.73287757 31.08558645 -49.78341272 
+DEAL::dst FEV: -37.20119468 142.7420109 -3.007994596 1.347129280 58.28259413 -229.2925423 20.73287757 31.08558645 -49.78341272 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::test_qiterated
+DEAL::Working with FE_Q<2>(1) and 25 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 -39.34040403 53.15927296 0.000000000 -83.88672176 36.18288001 -85.82630612 133.8837425 
+DEAL::dst FEV: 0.000000000 0.000000000 -39.34040403 53.15927296 0.000000000 -83.88672176 36.18288001 -85.82630612 133.8837425 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 79 dofs
+DEAL::dst FEE: -6.763525875 -29.84983569 8.787229829 29.57510339 35.20802355 -6.515494560 3.199548103 -17.09671493 61.28472734 
+DEAL::dst FEV: -6.763525875 -29.84983569 8.787229829 29.57510339 35.20802355 -6.515494560 3.199548103 -17.09671493 61.28472734 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 80 dofs
+DEAL::dst FEE: -9.799798542 6.710971666 13.84374919 -10.75492232 -1.428274410 1.045904167 1.947944821 -1.565574578 5.993067307 
+DEAL::dst FEV: -9.799798542 6.710971666 13.84374919 -10.75492232 -1.428274410 1.045904167 1.947944821 -1.565574578 5.993067307 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 448 dofs
+DEAL::dst FEE: -0.8390714818 3.507305214 2.477527164 -5.202526796 3.371341804 -5.961446157 -5.043747453 7.690617706 -1.205573925 
+DEAL::dst FEV: -0.8390714818 3.507305214 2.477527164 -5.202526796 3.371341804 -5.961446157 -5.043747453 7.690617706 -1.205573925 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 89 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 76.24995757 -151.4652529 139.3454674 112.8200482 0.000000000 -86.44907039 268.3337342 
+DEAL::dst FEV: 0.000000000 0.000000000 76.24995757 -151.4652529 139.3454674 112.8200482 0.000000000 -86.44907039 268.3337342 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 517 dofs
+DEAL::dst FEE: -109.3389581 -55.45400574 -15.87490632 -43.83124134 -209.9211798 40.69143970 -52.74249772 -18.77328792 -148.3336084 
+DEAL::dst FEV: -109.3389581 -55.45400574 -15.87490632 -43.83124134 -209.9211798 40.69143970 -52.74249772 -18.77328792 -148.3336084 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 180 dofs
+DEAL::dst FEE: 59.12248949 145.3263730 20.39413228 -145.1793073 -240.0765308 -62.93699888 99.98252044 67.73515150 55.63217032 
+DEAL::dst FEV: 59.12248949 145.3263730 20.39413228 -145.1793073 -240.0765308 -62.93699888 99.98252044 67.73515150 55.63217032 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 1512 dofs
+DEAL::dst FEE: 0.1855879691 22.48428433 4.315303570 -54.04857952 108.0103722 -48.99698879 23.31684388 -30.68686048 11.60957481 
+DEAL::dst FEV: 0.1855879691 22.48428433 4.315303570 -54.04857952 108.0103722 -48.99698879 23.31684388 -30.68686048 11.60957481 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 193 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 158.5095907 -1232.437308 572.2153320 -244.6069909 -498.6954353 994.6613382 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 158.5095907 -1232.437308 572.2153320 -244.6069909 -498.6954353 994.6613382 0.000000000 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 1651 dofs
+DEAL::dst FEE: -43.97463300 -73.92611758 -241.0763625 16.24844181 -292.3228855 -131.0871170 71.84204246 -200.0407478 -502.1313285 
+DEAL::dst FEV: -43.97463300 -73.92611758 -241.0763625 16.24844181 -292.3228855 -131.0871170 71.84204246 -200.0407478 -502.1313285 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 320 dofs
+DEAL::dst FEE: -79.52218037 -253.6335126 -828.4329925 -48.67475436 315.0517120 165.4890443 1618.475137 295.2068982 -521.4704040 
+DEAL::dst FEV: -79.52218037 -253.6335126 -828.4329925 -48.67475436 315.0517120 165.4890443 1618.475137 295.2068982 -521.4704040 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 3584 dofs
+DEAL::dst FEE: 5.976291794 41.32606645 28.53127386 15.47828625 29.83484294 -321.4991292 -201.2675788 50.72178853 -45.30966870 
+DEAL::dst FEV: 5.976291794 41.32606645 28.53127386 15.47828625 29.83484294 -321.4991292 -201.2675788 50.72178853 -45.30966870 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::test_hessians_with_gradients_and_values
+DEAL::test_qgauss
+DEAL::Working with FE_Q<2>(1) and 25 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 10.65897096 -39.52678335 0.000000000 103.0744784 -17.64340333 94.38631066 -164.5032332 
+DEAL::dst FEV: 0.000000000 0.000000000 10.65897096 -39.52678335 0.000000000 103.0744784 -17.64340333 94.38631066 -164.5032332 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 79 dofs
+DEAL::dst FEE: 1.102627435 -23.74365223 12.79134451 15.14153610 20.35616429 11.08290824 -9.664669182 14.19970889 -61.95383325 
+DEAL::dst FEV: 1.102627435 -23.74365223 12.79134451 15.14153610 20.35616429 11.08290824 -9.664669182 14.19970889 -61.95383325 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 80 dofs
+DEAL::dst FEE: -5.658488958 4.454226542 8.341036788 -7.037923219 3.277233078 -3.216068598 -4.875453948 4.916510798 -38.29668774 
+DEAL::dst FEV: -5.658488958 4.454226542 8.341036788 -7.037923219 3.277233078 -3.216068598 -4.875453948 4.916510798 -38.29668774 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 448 dofs
+DEAL::dst FEE: 3.537497829 -7.088058162 3.891919744 -0.3425928421 -1.788199416 5.329676319 -5.720480567 2.184775324 -1.566814782 
+DEAL::dst FEV: 3.537497829 -7.088058162 3.891919744 -0.3425928421 -1.788199416 5.329676319 -5.720480567 2.184775324 -1.566814782 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 89 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 227.5797473 350.2646716 14.46622536 -299.5647869 0.000000000 -246.4604494 81.46044522 
+DEAL::dst FEV: 0.000000000 0.000000000 227.5797473 350.2646716 14.46622536 -299.5647869 0.000000000 -246.4604494 81.46044522 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 517 dofs
+DEAL::dst FEE: -179.0434278 23.06902656 82.29275118 15.02013198 71.55602698 -124.1388389 48.89168750 20.79615558 269.8927219 
+DEAL::dst FEV: -179.0434278 23.06902656 82.29275118 15.02013198 71.55602698 -124.1388389 48.89168750 20.79615558 269.8927219 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 180 dofs
+DEAL::dst FEE: 39.15418328 -29.11997591 67.98689134 -159.4974460 182.2488961 -183.0088824 126.0094564 -155.9890929 112.3426688 
+DEAL::dst FEV: 39.15418328 -29.11997591 67.98689134 -159.4974460 182.2488961 -183.0088824 126.0094564 -155.9890929 112.3426688 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 1512 dofs
+DEAL::dst FEE: -1.274756990 -43.23619053 2.299787190 -40.17458947 170.9257619 8.426704364 11.03438425 -43.08065331 0.6067802486 
+DEAL::dst FEV: -1.274756990 -43.23619053 2.299787190 -40.17458947 170.9257619 8.426704364 11.03438425 -43.08065331 0.6067802486 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 193 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 390.6713500 -209.1815116 453.4456649 -417.3813654 133.4270439 -135.8576152 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 390.6713500 -209.1815116 453.4456649 -417.3813654 133.4270439 -135.8576152 0.000000000 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 1651 dofs
+DEAL::dst FEE: 154.9190172 -360.4883318 -356.1558831 157.2464180 117.8532397 252.9031712 -66.01918249 96.95024383 287.9441329 
+DEAL::dst FEV: 154.9190172 -360.4883318 -356.1558831 157.2464180 117.8532397 252.9031712 -66.01918249 96.95024383 287.9441329 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 320 dofs
+DEAL::dst FEE: -89.45060520 259.6951529 848.1451912 -41.23094175 178.0046639 -291.4637914 -1901.434532 0.9182705078 -205.8252742 
+DEAL::dst FEV: -89.45060520 259.6951529 848.1451912 -41.23094175 178.0046639 -291.4637914 -1901.434532 0.9182705078 -205.8252742 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 3584 dofs
+DEAL::dst FEE: 16.88952936 -26.19076464 166.9424606 -5.798811776 44.16364356 -140.5368456 -273.2726993 49.83102002 -17.72949070 
+DEAL::dst FEV: 16.88952936 -26.19076464 166.9424606 -5.798811776 44.16364356 -140.5368456 -273.2726993 49.83102002 -17.72949070 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::test_qiterated
+DEAL::Working with FE_Q<2>(1) and 25 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 11.52339170 33.12164009 0.000000000 -31.44757979 31.42806080 -91.11960114 74.83230770 
+DEAL::dst FEV: 0.000000000 0.000000000 11.52339170 33.12164009 0.000000000 -31.44757979 31.42806080 -91.11960114 74.83230770 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 79 dofs
+DEAL::dst FEE: -26.72770100 16.57344973 38.86689629 -30.27824765 33.97005374 -24.12586571 -42.87121363 43.97939426 -35.86348976 
+DEAL::dst FEV: -26.72770100 16.57344973 38.86689629 -30.27824765 33.97005374 -24.12586571 -42.87121363 43.97939426 -35.86348976 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 80 dofs
+DEAL::dst FEE: 3.371302359 -2.385095954 -5.105225396 4.213858047 -1.101697459 1.730125280 1.611434656 -2.168999319 23.58480366 
+DEAL::dst FEV: 3.371302359 -2.385095954 -5.105225396 4.213858047 -1.101697459 1.730125280 1.611434656 -2.168999319 23.58480366 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 448 dofs
+DEAL::dst FEE: 1.207718911 -3.743044786 4.895166874 -2.351876459 -3.005052078 5.530366961 -3.058534796 0.5306904929 0.1734889179 
+DEAL::dst FEV: 1.207718911 -3.743044786 4.895166874 -2.351876459 -3.005052078 5.530366961 -3.058534796 0.5306904929 0.1734889179 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 89 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 304.2620687 -12.11556339 17.19536074 -102.9230813 0.000000000 -102.4623557 -65.22429151 
+DEAL::dst FEV: 0.000000000 0.000000000 304.2620687 -12.11556339 17.19536074 -102.9230813 0.000000000 -102.4623557 -65.22429151 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 517 dofs
+DEAL::dst FEE: 176.7401543 -83.32894758 71.39944018 88.63852540 -190.9847564 -135.9080833 -40.88743431 53.54002808 -70.69405833 
+DEAL::dst FEV: 176.7401543 -83.32894758 71.39944018 88.63852540 -190.9847564 -135.9080833 -40.88743431 53.54002808 -70.69405833 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 180 dofs
+DEAL::dst FEE: -88.91974786 39.17050152 34.41914498 239.8242491 -151.0085069 -51.10455761 -177.7689307 155.8642570 -0.3653264004 
+DEAL::dst FEV: -88.91974786 39.17050152 34.41914498 239.8242491 -151.0085069 -51.10455761 -177.7689307 155.8642570 -0.3653264004 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 1512 dofs
+DEAL::dst FEE: -16.00930361 29.32089522 -4.691924057 13.64430384 -67.03254929 35.16094568 -2.368318354 0.8062879409 -2.621673503 
+DEAL::dst FEV: -16.00930361 29.32089522 -4.691924057 13.64430384 -67.03254929 35.16094568 -2.368318354 0.8062879409 -2.621673503 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 193 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 602.3311513 349.8347114 865.7547210 -1258.503940 -599.7204024 190.5161394 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 602.3311513 349.8347114 865.7547210 -1258.503940 -599.7204024 190.5161394 0.000000000 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 1651 dofs
+DEAL::dst FEE: 891.3558785 -64.82585347 144.7351978 -144.2971717 348.2630789 -131.1675690 -350.8844751 -25.03595897 145.5246669 
+DEAL::dst FEV: 891.3558785 -64.82585347 144.7351978 -144.2971717 348.2630789 -131.1675690 -350.8844751 -25.03595897 145.5246669 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 320 dofs
+DEAL::dst FEE: 75.49607652 271.9395872 206.8530002 51.99120011 -332.7199127 -226.1663066 -472.3855211 -160.5810108 800.9167395 
+DEAL::dst FEV: 75.49607652 271.9395872 206.8530002 51.99120011 -332.7199127 -226.1663066 -472.3855211 -160.5810108 800.9167395 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 3584 dofs
+DEAL::dst FEE: 10.13419666 26.50225800 -89.47868398 15.45030973 44.40065804 -53.19843171 279.5017995 -63.22517727 -2.347590852 
+DEAL::dst FEV: 10.13419666 26.50225800 -89.47868398 15.45030973 44.40065804 -53.19843171 279.5017995 -63.22517727 -2.347590852 
+DEAL::FEValues verification: 0.0
+DEAL::
diff --git a/tests/matrix_free/matrix_vector_hessians_general_faces.cc b/tests/matrix_free/matrix_vector_hessians_general_faces.cc
new file mode 100644 (file)
index 0000000..d327a0f
--- /dev/null
@@ -0,0 +1,383 @@
+// ---------------------------------------------------------------------
+//
+// 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/grid/grid_out.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/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+#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_ball(tria, Point<dim>(), 1., true);
+  tria.refine_global(1);
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+  MappingQGeneric<dim> mapping(2);
+
+  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.read_dof_values(src);
+          fe_eval_m.evaluate(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)
+{
+  using namespace dealii;
+  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<2>(i, dealii::FE_Q<2>(i), evaluation_flags);
+        test_hessians<3>(i, dealii::FE_Q<3>(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<2>(i, dealii::FE_Q<2>(i), evaluation_flags);
+        test_hessians<3>(i, dealii::FE_Q<3>(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<2>(i, dealii::FE_Q<2>(i), evaluation_flags);
+        test_hessians<3>(i, dealii::FE_Q<3>(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<2>(i, dealii::FE_Q<2>(i), evaluation_flags);
+        test_hessians<3>(i, dealii::FE_Q<3>(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_general_faces.output b/tests/matrix_free/matrix_vector_hessians_general_faces.output
new file mode 100644 (file)
index 0000000..90ff70b
--- /dev/null
@@ -0,0 +1,245 @@
+
+DEAL::test_hessians_only
+DEAL::Working with FE_Q<2>(1) and 25 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 1023.282938 -241.4191116 0.000000000 404.0338925 -1610.938705 830.1574688 -815.6478478 
+DEAL::dst FEV: 0.000000000 0.000000000 1023.282938 -241.4191116 0.000000000 404.0338925 -1610.938705 830.1574688 -815.6478478 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 79 dofs
+DEAL::dst FEE: 2085.339062 -984.5969128 -3156.855020 2414.250954 -2457.242215 914.7438039 1549.567722 -598.1729479 -1625.597093 
+DEAL::dst FEV: 2085.339062 -984.5969128 -3156.855020 2414.250954 -2457.242215 914.7438039 1549.567722 -598.1729479 -1625.597093 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 80 dofs
+DEAL::dst FEE: -70.25497193 56.55633625 99.11541602 -85.41678034 -18.84008428 22.15098884 28.29206961 -31.60297417 224.7350853 
+DEAL::dst FEV: -70.25497193 56.55633625 99.11541602 -85.41678034 -18.84008428 22.15098884 28.29206961 -31.60297417 224.7350853 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 448 dofs
+DEAL::dst FEE: 230.8404882 -276.9256039 -49.48465092 95.56976663 -173.4791192 219.5642349 -7.876718073 -38.20839764 154.6373788 
+DEAL::dst FEV: 230.8404882 -276.9256039 -49.48465092 95.56976663 -173.4791192 219.5642349 -7.876718073 -38.20839764 154.6373788 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 89 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 1995.727331 -5872.174659 4254.122091 4893.906481 0.000000000 5978.307265 -5884.314781 
+DEAL::dst FEV: 0.000000000 0.000000000 1995.727331 -5872.174659 4254.122091 4893.906481 0.000000000 5978.307265 -5884.314781 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 517 dofs
+DEAL::dst FEE: -18232.67381 -63.84155745 8315.900756 10959.22675 18566.85721 1642.297186 5567.506956 3577.708281 3777.934353 
+DEAL::dst FEV: -18232.67381 -63.84155745 8315.900756 10959.22675 18566.85721 1642.297186 5567.506956 3577.708281 3777.934353 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 180 dofs
+DEAL::dst FEE: 781.6486510 699.7352769 -40.96207937 -1563.480011 -1942.133396 521.4080939 724.0336595 1483.639308 -663.8895026 
+DEAL::dst FEV: 781.6486510 699.7352769 -40.96207937 -1563.480011 -1942.133396 521.4080939 724.0336595 1483.639308 -663.8895026 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 1512 dofs
+DEAL::dst FEE: -1145.622366 1317.428972 -661.2759029 3000.099316 -4245.714206 1889.713020 -394.8445158 -1418.351983 1342.575424 
+DEAL::dst FEV: -1145.622366 1317.428972 -661.2759029 3000.099316 -4245.714206 1889.713020 -394.8445158 -1418.351983 1342.575424 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 193 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 -40598.10535 46880.86712 3405.824484 -6942.530759 19670.93356 -35093.42245 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 -40598.10535 46880.86712 3405.824484 -6942.530759 19670.93356 -35093.42245 0.000000000 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 1651 dofs
+DEAL::dst FEE: 76024.91722 -21521.48567 7429.238626 34892.37890 61041.37781 39385.92760 -23643.19008 -28773.40415 -133970.1244 
+DEAL::dst FEV: 76024.91722 -21521.48567 7429.238626 34892.37890 61041.37781 39385.92760 -23643.19008 -28773.40415 -133970.1244 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 320 dofs
+DEAL::dst FEE: 1588.303554 -14029.09500 10487.91417 -4673.201800 -2725.353375 37066.63526 -27593.29159 11153.52368 5951.323359 
+DEAL::dst FEV: 1588.303554 -14029.09500 10487.91417 -4673.201800 -2725.353375 37066.63526 -27593.29159 11153.52368 5951.323359 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 3584 dofs
+DEAL::dst FEE: 374.2238240 -10778.38656 5315.827310 3061.915154 -3526.094877 27036.97679 -25348.35949 3956.883602 14680.05201 
+DEAL::dst FEV: 374.2238240 -10778.38656 5315.827310 3061.915154 -3526.094877 27036.97679 -25348.35949 3956.883602 14680.05201 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::test_hessians_with_values
+DEAL::Working with FE_Q<2>(1) and 25 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 -1584.165210 65.43910208 0.000000000 872.7146628 2609.017889 -88.44154865 -1485.472063 
+DEAL::dst FEV: 0.000000000 0.000000000 -1584.165210 65.43910208 0.000000000 872.7146628 2609.017889 -88.44154865 -1485.472063 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 79 dofs
+DEAL::dst FEE: -1621.730760 643.4195686 293.0631459 279.0288608 2304.285927 -1409.468004 -1137.774174 -131.8455981 984.9798012 
+DEAL::dst FEV: -1621.730760 643.4195686 293.0631459 279.0288608 2304.285927 -1409.468004 -1137.774174 -131.8455981 984.9798012 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 80 dofs
+DEAL::dst FEE: 57.37233989 -52.08768998 -81.52364724 76.90652431 -1.864777112 3.864526751 3.668178727 -5.283659601 -41.91915588 
+DEAL::dst FEV: 57.37233989 -52.08768998 -81.52364724 76.90652431 -1.864777112 3.864526751 3.668178727 -5.283659601 -41.91915588 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 448 dofs
+DEAL::dst FEE: -35.93990419 151.8225756 180.2088970 -296.0414848 -123.2946863 7.466972238 -20.93161579 136.8238756 -43.26334548 
+DEAL::dst FEV: -35.93990419 151.8225756 180.2088970 -296.0414848 -123.2946863 7.466972238 -20.93161579 136.8238756 -43.26334548 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 89 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 880.6007238 -2516.066873 7387.748162 1913.628818 0.000000000 2864.036480 -2737.954009 
+DEAL::dst FEV: 0.000000000 0.000000000 880.6007238 -2516.066873 7387.748162 1913.628818 0.000000000 2864.036480 -2737.954009 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 517 dofs
+DEAL::dst FEE: 18596.12613 15528.57481 6884.355281 5777.026592 -2674.374463 -1007.329718 4432.749714 -2911.678883 -13831.23738 
+DEAL::dst FEV: 18596.12613 15528.57481 6884.355281 5777.026592 -2674.374463 -1007.329718 4432.749714 -2911.678883 -13831.23738 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 180 dofs
+DEAL::dst FEE: 103.2620949 532.2380623 -633.7928408 419.9454288 -2675.339147 2218.294506 -622.4386240 2469.622503 -1811.077935 
+DEAL::dst FEV: 103.2620949 532.2380623 -633.7928408 419.9454288 -2675.339147 2218.294506 -622.4386240 2469.622503 -1811.077935 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 1512 dofs
+DEAL::dst FEE: -2094.840483 3896.320685 -974.7402212 3024.687552 -6977.400539 2671.522120 -1333.566959 3011.785402 -1481.004655 
+DEAL::dst FEV: -2094.840483 3896.320685 -974.7402212 3024.687552 -6977.400539 2671.522120 -1333.566959 3011.785402 -1481.004655 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 193 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 126464.1026 -1087.849656 22905.39195 -37835.72538 -736.0791996 9064.904645 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 126464.1026 -1087.849656 22905.39195 -37835.72538 -736.0791996 9064.904645 0.000000000 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 1651 dofs
+DEAL::dst FEE: -11637.00801 -886.8694441 472.9210085 -28012.52627 36260.35838 -16255.65596 -35494.55425 -2416.417224 123004.4958 
+DEAL::dst FEV: -11637.00801 -886.8694441 472.9210085 -28012.52627 36260.35838 -16255.65596 -35494.55425 -2416.417224 123004.4958 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 320 dofs
+DEAL::dst FEE: -3827.887005 -560.1149871 6336.511188 -2424.587270 7291.112492 3998.764339 -19339.34113 6836.645423 613.0280050 
+DEAL::dst FEV: -3827.887005 -560.1149871 6336.511188 -2424.587270 7291.112492 3998.764339 -19339.34113 6836.645423 613.0280050 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 3584 dofs
+DEAL::dst FEE: 3267.074706 -11654.92306 15702.51697 -2120.107249 -3708.542275 13910.79441 -12561.65555 -4917.218117 -643.8540446 
+DEAL::dst FEV: 3267.074706 -11654.92306 15702.51697 -2120.107249 -3708.542275 13910.79441 -12561.65555 -4917.218117 -643.8540446 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::test_hessians_with_gradients
+DEAL::Working with FE_Q<2>(1) and 25 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 -621.4111919 114.1857906 0.000000000 -564.6157944 1395.687107 -1160.831541 1354.790422 
+DEAL::dst FEV: 0.000000000 0.000000000 -621.4111919 114.1857906 0.000000000 -564.6157944 1395.687107 -1160.831541 1354.790422 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 79 dofs
+DEAL::dst FEE: 1278.032634 776.1433800 -1937.270825 -166.3396816 164.8182756 646.8512138 73.41373899 27.59426419 -1571.717176 
+DEAL::dst FEV: 1278.032634 776.1433800 -1937.270825 -166.3396816 164.8182756 646.8512138 73.41373899 27.59426419 -1571.717176 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 80 dofs
+DEAL::dst FEE: 75.36615516 -64.43827144 -107.2029444 96.27506071 -10.15977218 14.65805710 15.33219542 -19.83048034 485.1973589 
+DEAL::dst FEV: 75.36615516 -64.43827144 -107.2029444 96.27506071 -10.15977218 14.65805710 15.33219542 -19.83048034 485.1973589 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 448 dofs
+DEAL::dst FEE: -147.6510186 280.3568712 58.16675185 -191.7314131 12.97777322 -145.4542511 76.29975751 57.03552894 21.57976763 
+DEAL::dst FEV: -147.6510186 280.3568712 58.16675185 -191.7314131 12.97777322 -145.4542511 76.29975751 57.03552894 21.57976763 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 89 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 22159.26719 2260.566599 -2765.305051 714.0319473 0.000000000 -6768.530153 2345.620782 
+DEAL::dst FEV: 0.000000000 0.000000000 22159.26719 2260.566599 -2765.305051 714.0319473 0.000000000 -6768.530153 2345.620782 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 517 dofs
+DEAL::dst FEE: -3899.489149 15520.20968 5914.037728 -1363.044631 -15353.90033 2244.139771 4155.919950 3503.216608 22003.29974 
+DEAL::dst FEV: -3899.489149 15520.20968 5914.037728 -1363.044631 -15353.90033 2244.139771 4155.919950 3503.216608 22003.29974 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 180 dofs
+DEAL::dst FEE: -1292.549827 1556.133239 6.674065417 3834.116898 -4965.432205 581.6592615 -2903.715710 4050.659028 -867.5447504 
+DEAL::dst FEV: -1292.549827 1556.133239 6.674065417 3834.116898 -4965.432205 581.6592615 -2903.715710 4050.659028 -867.5447504 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 1512 dofs
+DEAL::dst FEE: -169.5401972 -2238.676603 552.5699787 2580.643343 483.3065112 -585.8533834 -1693.719654 1992.806308 -614.3341662 
+DEAL::dst FEV: -169.5401972 -2238.676603 552.5699787 2580.643343 483.3065112 -585.8533834 -1693.719654 1992.806308 -614.3341662 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 193 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 56775.09128 12613.03416 12887.52510 -18193.82181 17692.30415 -22298.47264 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 56775.09128 12613.03416 12887.52510 -18193.82181 17692.30415 -22298.47264 0.000000000 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 1651 dofs
+DEAL::dst FEE: -42336.21123 -26120.26166 -2297.564253 795.3646393 195.9688350 -26011.43162 -6710.972876 -1437.849217 111350.4393 
+DEAL::dst FEV: -42336.21123 -26120.26166 -2297.564253 795.3646393 195.9688350 -26011.43162 -6710.972876 -1437.849217 111350.4393 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 320 dofs
+DEAL::dst FEE: -5791.865257 1722.236271 -6101.199466 3515.520543 12454.33999 -1254.417975 14202.82726 -7529.176527 -12023.46778 
+DEAL::dst FEV: -5791.865257 1722.236271 -6101.199466 3515.520543 12454.33999 -1254.417975 14202.82726 -7529.176527 -12023.46778 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 3584 dofs
+DEAL::dst FEE: 873.2174032 -6121.609196 -13551.64405 4309.630048 -4242.986432 9017.066410 28639.88440 -10457.78950 5235.775871 
+DEAL::dst FEV: 873.2174032 -6121.609196 -13551.64405 4309.630048 -4242.986432 9017.066410 28639.88440 -10457.78950 5235.775871 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::test_hessians_with_gradients_and_values
+DEAL::Working with FE_Q<2>(1) and 25 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 520.1325322 -690.2874624 0.000000000 1816.271572 -973.5840476 2143.792127 -3254.583688 
+DEAL::dst FEV: 0.000000000 0.000000000 520.1325322 -690.2874624 0.000000000 1816.271572 -973.5840476 2143.792127 -3254.583688 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(1) and 79 dofs
+DEAL::dst FEE: 3244.549979 -1911.946038 -996.0944826 666.3414459 -1630.451414 -1193.463312 174.5157649 1414.847115 1650.753381 
+DEAL::dst FEV: 3244.549979 -1911.946038 -996.0944826 666.3414459 -1630.451414 -1193.463312 174.5157649 1414.847115 1650.753381 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(1) and 80 dofs
+DEAL::dst FEE: 22.28659957 -15.47011943 -34.59406422 28.35461984 -4.581461453 13.39644349 8.875622598 -17.53301884 575.4344949 
+DEAL::dst FEV: 22.28659957 -15.47011943 -34.59406422 28.35461984 -4.581461453 13.39644349 8.875622598 -17.53301884 575.4344949 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(1) and 448 dofs
+DEAL::dst FEE: -44.84164597 -20.30585088 186.0485414 -121.5573779 -53.47487885 118.0871478 -87.61030991 23.83004676 -39.39444906 
+DEAL::dst FEV: -44.84164597 -20.30585088 186.0485414 -121.5573779 -53.47487885 118.0871478 -87.61030991 23.83004676 -39.39444906 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(2) and 89 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 -7209.671714 -3603.056480 4376.848314 994.6112263 0.000000000 2021.933328 -981.7698190 
+DEAL::dst FEV: 0.000000000 0.000000000 -7209.671714 -3603.056480 4376.848314 994.6112263 0.000000000 2021.933328 -981.7698190 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(2) and 517 dofs
+DEAL::dst FEE: 1367.885593 4424.107523 5513.989531 -7434.495343 -3908.076539 6480.955358 -330.9300163 -1136.937005 -6414.069171 
+DEAL::dst FEV: 1367.885593 4424.107523 5513.989531 -7434.495343 -3908.076539 6480.955358 -330.9300163 -1136.937005 -6414.069171 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(2) and 180 dofs
+DEAL::dst FEE: -437.8541486 -496.2779202 809.0422323 1332.294523 719.5400451 -1784.743121 -1085.021345 29.27757678 914.5814941 
+DEAL::dst FEV: -437.8541486 -496.2779202 809.0422323 1332.294523 719.5400451 -1784.743121 -1085.021345 29.27757678 914.5814941 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(2) and 1512 dofs
+DEAL::dst FEE: -1566.061648 375.3459780 1252.028307 3338.351514 -3056.749505 -440.7156811 -382.1943741 -181.2283054 682.1878856 
+DEAL::dst FEV: -1566.061648 375.3459780 1252.028307 3338.351514 -3056.749505 -440.7156811 -382.1943741 -181.2283054 682.1878856 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<2>(3) and 193 dofs
+DEAL::dst FEE: 0.000000000 0.000000000 -72513.09849 2440.045444 5936.767657 366.5967395 13588.30105 -15481.35114 0.000000000 
+DEAL::dst FEV: 0.000000000 0.000000000 -72513.09849 2440.045444 5936.767657 366.5967395 13588.30105 -15481.35114 0.000000000 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_Q<3>(3) and 1651 dofs
+DEAL::dst FEE: 33810.56033 -17394.35348 34025.74741 -12344.52484 48066.48418 4634.340963 -24021.08665 26673.43063 30975.17620 
+DEAL::dst FEV: 33810.56033 -17394.35348 34025.74741 -12344.52484 48066.48418 4634.340963 -24021.08665 26673.43063 30975.17620 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<2>(3) and 320 dofs
+DEAL::dst FEE: 5527.094242 -10178.50933 242.7247396 -696.3185404 -11066.69999 24621.98310 2200.859720 1533.485743 11062.88761 
+DEAL::dst FEV: 5527.094242 -10178.50933 242.7247396 -696.3185404 -11066.69999 24621.98310 2200.859720 1533.485743 11062.88761 
+DEAL::FEValues verification: 0.0
+DEAL::
+DEAL::Working with FE_DGQ<3>(3) and 3584 dofs
+DEAL::dst FEE: 7134.467624 -3170.261114 -745.5916238 -4621.556435 -6399.737030 -7674.914630 14942.24049 -3133.337521 5780.416539 
+DEAL::dst FEV: 7134.467624 -3170.261114 -745.5916238 -4621.556435 -6399.737030 -7674.914630 14942.24049 -3133.337521 5780.416539 
+DEAL::FEValues verification: 0.0
+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.