From 0c5d0bcdc47ec4ebca85e98571695cb67c78cf25 Mon Sep 17 00:00:00 2001 From: Maximilian Bergbauer Date: Fri, 19 Jul 2024 17:53:32 +0200 Subject: [PATCH] Implement normal hessian --- include/deal.II/matrix_free/fe_evaluation.h | 188 ++++++++++++++++++ .../matrix_vector_normal_hessians.cc | 179 +++++++++++++++++ .../matrix_vector_normal_hessians.output | 26 +++ 3 files changed, 393 insertions(+) create mode 100644 tests/matrix_free/matrix_vector_normal_hessians.cc create mode 100644 tests/matrix_free/matrix_vector_normal_hessians.output diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 78ddb7e298..8a36d9e0db 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -469,6 +469,18 @@ public: value_type get_laplacian(const unsigned int q_point) const; + /** + * Return the second derivative along the normal direction $\partial_{n}^2 + * u_h$ (i.e., the Hessian of the function $u_h$ contracted twice with the + * direction of the normal vector) of the finite element + * function interpolated to the quadrature point with index + * @p q_point after a call to @p evaluate(EvaluationFlags::hessians). + * Compared to the case when computing the full Hessian, some operations can + * be saved when only the normal Hessian is requested. + */ + value_type + get_normal_hessian(const unsigned int q_point) const; + /** * Write a contribution that gets multiplied by the Hessian of the test * function to the field containing the Hessians at quadrature points with @@ -486,6 +498,25 @@ public: void submit_hessian(const hessian_type hessian_in, const unsigned int q_point); + /** + * Write a contribution that gets multiplied by the Hessian of the test + * function times the normal projector to the field containing the Hessians at + * quadrature points with index @p q_point. When this function has queued + * information for all quadrature points and followed by a call to the + * function @p integrate(EvaluationFlags::hessians), the result is an + * array of entries, each representing the result of the integral of product + * of the test function Hessian multiplied by the values times the normal + * projector passed to this function. + * + * @note This function accesses the same field as through get_hessian() and + * get_normal_hessian() so make sure to not call it after calling + * submit_hessian() or submit_normal_hessian() for a specific quadrature point + * index. + */ + void + submit_normal_hessian(const value_type normal_hessian_in, + const unsigned int q_point); + /** * Return the divergence of a vector-valued finite element at quadrature * point number @p q_point after a call to @@ -4777,6 +4808,85 @@ FEEvaluationBase:: +template +inline typename FEEvaluationBase::value_type +FEEvaluationBase:: + get_normal_hessian(const unsigned int q_point) const +{ +# ifdef DEBUG + Assert(this->hessians_quad_initialized == true, + internal::ExcAccessToUninitializedField()); +# endif + AssertIndexRange(q_point, this->n_quadrature_points); + + Assert(this->normal_x_jacobian != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_hessians")); + + Tensor<1, n_components, VectorizedArrayType> hessian_out; + + const std::size_t nqp = this->n_quadrature_points; + constexpr unsigned int hdim = (dim * (dim + 1)) / 2; + + if (this->cell_type <= internal::MatrixFreeFunctions::affine) + { + const auto nxj = this->normal_x_jacobian[0]; + + for (unsigned int comp = 0; comp < n_components; ++comp) + { + for (unsigned int d = 0; d < dim; ++d) + hessian_out[comp] += + this->hessians_quad[(comp * hdim + d) * nqp + q_point] * + (nxj[d]) * (nxj[d]); + + switch (dim) + { + case 1: + break; + case 2: + hessian_out[comp] += + this->hessians_quad[(comp * hdim + 2) * nqp + q_point] * + (nxj[0] * nxj[1]); + break; + case 3: + hessian_out[comp] += + 2. * this->hessians_quad[(comp * hdim + 3) * nqp + q_point] * + (nxj[0] * nxj[1]); + hessian_out[comp] += + 2. * this->hessians_quad[(comp * hdim + 4) * nqp + q_point] * + (nxj[0] * nxj[2]); + hessian_out[comp] += + 2. * this->hessians_quad[(comp * hdim + 5) * nqp + q_point] * + (nxj[1] * nxj[2]); + break; + default: + DEAL_II_NOT_IMPLEMENTED(); + } + } + + if constexpr (n_components == 1) + return hessian_out[0]; + else + return hessian_out; + } + // cell with general Jacobian + else + { + const auto normal = this->normal_vector(q_point); + return get_hessian(q_point) * normal * normal; + } +} + + + template :: +template +inline DEAL_II_ALWAYS_INLINE void +FEEvaluationBase:: + submit_normal_hessian(const value_type normal_hessian_in, + const unsigned int q_point) +{ +# ifdef DEBUG + Assert(this->is_reinitialized, ExcNotInitialized()); +# endif + 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 (this->cell_type <= internal::MatrixFreeFunctions::affine) + { + const VectorizedArrayType JxW = + this->J_value[0] * this->quadrature_weights[q_point]; + + const auto nxj = this->normal_x_jacobian[0]; + + // diagonal part + for (unsigned int d = 0; d < dim; ++d) + { + const auto nxj_d = nxj[d]; + const VectorizedArrayType factor = nxj_d * nxj_d * JxW; + for (unsigned int comp = 0; comp < n_components; ++comp) + if constexpr (n_components == 1) + this->hessians_quad[d * nqp + q_point] = + normal_hessian_in * factor; + else + this->hessians_quad[(comp * hdim + d) * nqp + q_point] = + normal_hessian_in[comp] * factor; + } + + // off diagonal part + for (unsigned int d = 1, off_dia = dim; d < dim; ++d) + for (unsigned int e = 0; e < d; ++e, ++off_dia) + { + const auto jac_d = nxj[d]; + const auto jac_e = nxj[e]; + const VectorizedArrayType factor = jac_d * jac_e * JxW; + for (unsigned int comp = 0; comp < n_components; ++comp) + if constexpr (n_components == 1) + this->hessians_quad[off_dia * nqp + q_point] = + 2. * normal_hessian_in * factor; + else + this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] = + 2. * normal_hessian_in[comp] * factor; + } + } + else + { + const auto normal = this->normal_vector(q_point); + const auto normal_projector = outer_product(normal, normal); + if constexpr (n_components == 1) + submit_hessian(normal_hessian_in * normal_projector, q_point); + else + submit_hessian(outer_product(normal_hessian_in, normal_projector), + q_point); + } +} + + + template + +#include + +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + + +template +void +test_hessians(const unsigned int degree, + const dealii::FE_Poly &fe, + const dealii::EvaluationFlags::EvaluationFlags evaluation_flags) +{ + using VectorizedArrayType = VectorizedArray; + + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(4 - dim); + + 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; + DoFTools::make_hanging_node_constraints(dof_handler, 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); + + src.update_ghost_values(); + + 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_normal_hessian(fe_eval_m.get_normal_hessian(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_normal_hessian(fe_eval_p.get_normal_hessian(q), + q); + } + fe_eval_p.integrate(evaluation_flags); + fe_eval_p.distribute_local_to_global(dst); + } + }, + [&]( + 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; +} + +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 < 3; ++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_normal_hessians.output b/tests/matrix_free/matrix_vector_normal_hessians.output new file mode 100644 index 0000000000..24105ddbeb --- /dev/null +++ b/tests/matrix_free/matrix_vector_normal_hessians.output @@ -0,0 +1,26 @@ + +DEAL::test_hessians_only +DEAL::Working with FE_Q<1>(1) and 12 dofs +DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 +DEAL::Working with FE_Q<2>(1) and 51 dofs +DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 +DEAL::Working with FE_Q<3>(1) and 81 dofs +DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 +DEAL::Working with FE_DGQ<1>(1) and 22 dofs +DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 +DEAL::Working with FE_DGQ<2>(1) and 136 dofs +DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 +DEAL::Working with FE_DGQ<3>(1) and 232 dofs +DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 +DEAL::Working with FE_Q<1>(2) and 23 dofs +DEAL::dst FEE: 1465345.061 40571.23729 11202.30012 186184.9221 -92344.77470 1048204.969 49574.60308 10790.27125 666909.3801 +DEAL::Working with FE_Q<2>(2) and 181 dofs +DEAL::dst FEE: 0.000000000 0.000000000 -11330.88975 -5767.500760 -10780.28907 -1865.442471 0.000000000 -1626.733233 3933.131615 +DEAL::Working with FE_Q<3>(2) and 446 dofs +DEAL::dst FEE: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 25.74078069 0.000000000 0.000000000 +DEAL::Working with FE_DGQ<1>(2) and 33 dofs +DEAL::dst FEE: 9211.388432 -18422.77686 9211.388432 49499.63824 -98999.27647 49499.63824 -11803.61329 23607.22657 -11803.61329 +DEAL::Working with FE_DGQ<2>(2) and 306 dofs +DEAL::dst FEE: -87.24207970 115.4371186 -231.0527514 560.7971199 -1003.500158 848.4184633 -161.0790796 263.1111184 -304.8897513 +DEAL::Working with FE_DGQ<3>(2) and 783 dofs +DEAL::dst FEE: 2.138422207 4.212276390 -2.105358266 -0.9089585142 -15.50152582 0.1505972516 -0.1426991866 0.9662025725 0.2225716328 -- 2.39.5