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