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
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
+template <int dim,
+ int n_components_,
+ typename Number,
+ bool is_face,
+ typename VectorizedArrayType>
+inline typename FEEvaluationBase<dim,
+ n_components_,
+ Number,
+ is_face,
+ VectorizedArrayType>::value_type
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+ 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 <int dim,
int n_components_,
typename Number,
+template <int dim,
+ int n_components_,
+ typename Number,
+ bool is_face,
+ typename VectorizedArrayType>
+inline DEAL_II_ALWAYS_INLINE void
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+ submit_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 <int dim,
int n_components_,
typename Number,
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2021 - 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+// Test 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.
+
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q_generic.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+test_hessians(const unsigned int degree,
+ const dealii::FE_Poly<dim> &fe,
+ const dealii::EvaluationFlags::EvaluationFlags evaluation_flags)
+{
+ using VectorizedArrayType = VectorizedArray<double>;
+
+ Triangulation<dim> 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<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+ MappingQGeneric<dim> mapping(1);
+
+ AffineConstraints<double> constraints;
+ DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+ VectorTools::interpolate_boundary_values(
+ mapping, dof_handler, 0, Functions::ZeroFunction<dim>(), constraints);
+ constraints.close();
+
+ // FEFaceEvaluation
+ typename MatrixFree<dim, double, VectorizedArrayType>::AdditionalData
+ additional_data;
+ additional_data.mapping_update_flags_inner_faces =
+ update_values | update_gradients | update_hessians;
+
+ MatrixFree<dim, double, VectorizedArrayType> matrix_free;
+ matrix_free.reinit(
+ mapping, dof_handler, constraints, QGauss<1>(degree + 1), additional_data);
+
+ if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+ deallog << "Working with " << fe.get_name() << " and "
+ << dof_handler.n_dofs() << " dofs" << std::endl;
+
+ LinearAlgebra::distributed::Vector<double> src, dst, dst2;
+ matrix_free.initialize_dof_vector(src);
+ for (auto &v : src)
+ v = random_value<double>();
+
+ matrix_free.initialize_dof_vector(dst);
+
+ src.update_ghost_values();
+
+ matrix_free.template loop<LinearAlgebra::distributed::Vector<double>,
+ LinearAlgebra::distributed::Vector<double>>(
+ [&](
+ const auto &matrix_free, auto &dst, const auto &src, const auto &range) {
+ (void)matrix_free;
+ (void)dst;
+ (void)src;
+ (void)range;
+ },
+ [&](
+ const auto &matrix_free, auto &dst, const auto &src, const auto &range) {
+ FEFaceEvaluation<dim, -1, 0, 1, double, VectorizedArrayType> fe_eval_m(
+ matrix_free, true);
+ FEFaceEvaluation<dim, -1, 0, 1, double, VectorizedArrayType> fe_eval_p(
+ matrix_free, false);
+ for (unsigned int face = range.first; face < range.second; ++face)
+ {
+ // FEFaceEvaluation
+ fe_eval_m.reinit(face);
+ fe_eval_m.gather_evaluate(src, evaluation_flags);
+ for (unsigned int q = 0; q < fe_eval_m.n_q_points; ++q)
+ {
+ if (evaluation_flags & EvaluationFlags::hessians)
+ fe_eval_m.submit_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);
+ }
+ }
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff