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)
*/
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
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
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);
}
[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);
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)
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)
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;
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)
{
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
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;
+ }
+ }
}
}
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,
else
internal::FEEvaluationFactory<dim, Number, VectorizedArrayType>::evaluate(
n_components,
- evaluation_flags,
+ evaluation_flag_actual,
*this->data,
const_cast<VectorizedArrayType *>(values_array),
this->values_quad,
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
}
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(),
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,
else
internal::FEEvaluationFactory<dim, Number, VectorizedArrayType>::integrate(
n_components,
- integration_flag,
+ integration_flag_actual,
*this->data,
values_array,
this->values_quad,
!(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)
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],
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->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
}
!(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>(
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->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,
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;
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(
- /* ------------------------- initialization of cells ------------------- */
-
// Copy a vectorized array of one type to another type
template <typename VectorizedArrayType1, typename VectorizedArrayType2>
inline DEAL_II_ALWAYS_INLINE void
+ // 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
}
}
- // 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.
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));
}
}
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(
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)
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>
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] =
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+ }
+}
--- /dev/null
+
+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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+ }
+}
--- /dev/null
+
+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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+ }
+ }
+}
--- /dev/null
+
+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::