Implement normal_derivative for FEPointEvaluation

diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h
index ef56931f7e4b1d41494b72c7b2a54f1fcc017793..b8c28ad587cd16b21f352edda4bded1775a08cdc 100644 (file)
--- a/include/deal.II/matrix_free/fe_point_evaluation.h
+++ b/include/deal.II/matrix_free/fe_point_evaluation.h
@@ -1349,6 +1349,22 @@ public:
   Tensor<1, spacedim, Number>
   normal_vector(const unsigned int point_index) const;
 
+  /**
+   * Return the normal gradient in real coordinates at the point with index
+   * `point_index` after a call to FEPointEvaluation::evaluate() with
+   * EvaluationFlags::gradients set.
+   */
+  const value_type
+  get_normal_derivative(const unsigned int point_index) const;
+
+  /**
+   * Write a contribution that is tested by the normal gradient to the field
+   * containing the values on points with the given `point_index`. Access to
+   * the same field as through set_gradient()/get_gradient.
+   */
+  void
+  submit_normal_derivative(const value_type &, const unsigned int point_index);
+
 private:
   static constexpr std::size_t n_lanes_user_interface =
     internal::VectorizedArrayTrait<Number>::width();
@@ -1707,6 +1723,22 @@ public:
   Tensor<1, spacedim, Number>
   normal_vector(const unsigned int point_index) const;
 
+  /**
+   * Return the normal gradient in real coordinates at the point with index
+   * `point_index` after a call to FEFacePointEvaluation::evaluate() with
+   * EvaluationFlags::gradients set.
+   */
+  const value_type
+  get_normal_derivative(const unsigned int point_index) const;
+
+  /**
+   * Write a contribution that is tested by the normal gradient to the field
+   * containing the values on points with the given `point_index`. Access to
+   * the same field as through set_gradient()/get_gradient.
+   */
+  void
+  submit_normal_derivative(const value_type &, const unsigned int point_index);
+
 private:
   static constexpr std::size_t n_lanes_user_interface =
     internal::VectorizedArrayTrait<Number>::width();
@@ -3086,6 +3118,36 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::normal_vector(
 
 
 
+template <int n_components_, int dim, int spacedim, typename Number>
+inline const typename FEPointEvaluation<n_components_,
+                                        dim,
+                                        spacedim,
+                                        Number>::value_type
+FEPointEvaluation<n_components_, dim, spacedim, Number>::get_normal_derivative(
+  const unsigned int point_index) const
+{
+  AssertIndexRange(point_index, this->gradients.size());
+  return this->gradients[point_index] * normal_vector(point_index);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::
+  submit_normal_derivative(const value_type  &value,
+                           const unsigned int point_index)
+{
+  AssertIndexRange(point_index, this->gradients.size());
+  if constexpr (n_components == 1)
+    this->gradients[point_index] = value * normal_vector(point_index);
+  else
+    this->gradients[point_index] =
+      outer_product(value, normal_vector(point_index));
+}
+
+
+
 template <int n_components_, int dim, int spacedim, typename Number>
 FEFacePointEvaluation<n_components_, dim, spacedim, Number>::
   FEFacePointEvaluation(
@@ -3817,6 +3879,36 @@ FEFacePointEvaluation<n_components_, dim, spacedim, Number>::normal_vector(
     }
 }
 
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline const typename FEFacePointEvaluation<n_components_,
+                                            dim,
+                                            spacedim,
+                                            Number>::value_type
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::
+  get_normal_derivative(const unsigned int point_index) const
+{
+  AssertIndexRange(point_index, this->gradients.size());
+  return this->gradients[point_index] * normal_vector(point_index);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline void
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::
+  submit_normal_derivative(const value_type  &value,
+                           const unsigned int point_index)
+{
+  AssertIndexRange(point_index, this->gradients.size());
+  if constexpr (n_components == 1)
+    this->gradients[point_index] = value * normal_vector(point_index);
+  else
+    this->gradients[point_index] =
+      outer_product(value, normal_vector(point_index));
+}
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif