Implement normal hessian

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 78ddb7e2988a7ee3df386fb38e283d88d5868659..8a36d9e0dbd992b78281bb571ed38b8d02a8842b 100644 (file)
--- 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<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
 
 
+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,
@@ -5428,6 +5538,84 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
 
 
+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,

diff --git a/tests/matrix_free/matrix_vector_normal_hessians.cc b/tests/matrix_free/matrix_vector_normal_hessians.cc
new file mode 100644 (file)
index 0000000..4df6c53
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_normal_hessians.cc
@@ -0,0 +1,179 @@
+// ------------------------------------------------------------------------
+//
+// 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);
+      }
+  }
+}

diff --git a/tests/matrix_free/matrix_vector_normal_hessians.output b/tests/matrix_free/matrix_vector_normal_hessians.output
new file mode 100644 (file)
index 0000000..24105dd
--- /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