From: Simon Sticko <simon@sticko.se>
Date: Fri, 14 Feb 2020 13:09:21 +0000 (+0100)
Subject: Add function FE_Poly::correct_hessians.
X-Git-Tag: v9.2.0-rc1~401^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=dc9ad3c3ebe36af31d6ee53731fc115a232567a2;p=dealii.git

Add function FE_Poly::correct_hessians.

In the same way as FE_Poly::correct_third_derivatives, collect the
terms involing the derivative of the Jacobian in a function
correct_hessians.
---

diff --git a/include/deal.II/fe/fe_poly.h b/include/deal.II/fe/fe_poly.h
index b8f30add18..31628d934a 100644
--- a/include/deal.II/fe/fe_poly.h
+++ b/include/deal.II/fe/fe_poly.h
@@ -454,6 +454,30 @@ protected:
     Table<2, Tensor<3, dim>> shape_3rd_derivatives;
   };
 
+  /**
+   * Correct the shape Hessians by subtracting the terms corresponding to the
+   * Jacobian pushed forward gradient.
+   *
+   * Before the correction, the Hessians would be given by
+   * @f[
+   * D_{ijk} = \frac{d^2\phi_i}{d \hat x_J d \hat x_K} (J_{jJ})^{-1}
+   * (J_{kK})^{-1},
+   * @f]
+   * where $J_{iI}=\frac{d x_i}{d \hat x_I}$. After the correction, the
+   * correct Hessians would be given by
+   * @f[
+   * \frac{d^2 \phi_i}{d x_j d x_k} = D_{ijk} - H_{mjk} \frac{d \phi_i}{d x_m},
+   * @f]
+   * where $H_{ijk}$ is the Jacobian pushed-forward derivative.
+   */
+  void
+  correct_hessians(
+    internal::FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>
+      &output_data,
+    const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+      &                mapping_data,
+    const unsigned int n_q_points) const;
+
   /**
    * Correct the shape third derivatives by subtracting the terms
    * corresponding to the Jacobian pushed forward gradient and second
@@ -485,6 +509,7 @@ protected:
     const unsigned int n_q_points,
     const unsigned int dof) const;
 
+
   /**
    * The polynomial space. Its type is given by the template parameter
    * PolynomialType.
diff --git a/include/deal.II/fe/fe_poly.templates.h b/include/deal.II/fe/fe_poly.templates.h
index 55d395d98c..890a6e7a2d 100644
--- a/include/deal.II/fe/fe_poly.templates.h
+++ b/include/deal.II/fe/fe_poly.templates.h
@@ -268,11 +268,7 @@ FE_Poly<PolynomialType, dim, spacedim>::fill_fe_values(
                           mapping_internal,
                           make_array_view(output_data.shape_hessians, k));
 
-      for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
-        for (unsigned int i = 0; i < quadrature.size(); ++i)
-          output_data.shape_hessians[k][i] -=
-            output_data.shape_gradients[k][i] *
-            mapping_data.jacobian_pushed_forward_grads[i];
+      correct_hessians(output_data, mapping_data, quadrature.size());
     }
 
   if (flags & update_3rd_derivatives &&
@@ -357,12 +353,7 @@ FE_Poly<PolynomialType, dim, spacedim>::fill_fe_face_values(
           mapping_internal,
           make_array_view(output_data.shape_hessians, k));
 
-      for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
-        for (unsigned int i = 0; i < quadrature.size(); ++i)
-          for (unsigned int j = 0; j < spacedim; ++j)
-            output_data.shape_hessians[k][i] -=
-              mapping_data.jacobian_pushed_forward_grads[i][j] *
-              output_data.shape_gradients[k][i][j];
+      correct_hessians(output_data, mapping_data, quadrature.size());
     }
 
   if (flags & update_3rd_derivatives)
@@ -452,12 +443,7 @@ FE_Poly<PolynomialType, dim, spacedim>::fill_fe_subface_values(
           mapping_internal,
           make_array_view(output_data.shape_hessians, k));
 
-      for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
-        for (unsigned int i = 0; i < quadrature.size(); ++i)
-          for (unsigned int j = 0; j < spacedim; ++j)
-            output_data.shape_hessians[k][i] -=
-              mapping_data.jacobian_pushed_forward_grads[i][j] *
-              output_data.shape_gradients[k][i][j];
+      correct_hessians(output_data, mapping_data, quadrature.size());
     }
 
   if (flags & update_3rd_derivatives)
@@ -482,6 +468,25 @@ FE_Poly<PolynomialType, dim, spacedim>::fill_fe_subface_values(
 
 
 
+template <class PolynomialType, int dim, int spacedim>
+inline void
+FE_Poly<PolynomialType, dim, spacedim>::correct_hessians(
+  internal::FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>
+    &output_data,
+  const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+    &                mapping_data,
+  const unsigned int n_q_points) const
+{
+  for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
+    for (unsigned int i = 0; i < n_q_points; ++i)
+      for (unsigned int j = 0; j < spacedim; ++j)
+        output_data.shape_hessians[dof][i] -=
+          mapping_data.jacobian_pushed_forward_grads[i][j] *
+          output_data.shape_gradients[dof][i][j];
+}
+
+
+
 template <class PolynomialType, int dim, int spacedim>
 inline void
 FE_Poly<PolynomialType, dim, spacedim>::correct_third_derivatives(
diff --git a/source/fe/fe_poly.cc b/source/fe/fe_poly.cc
index d845fd05e0..016cd82715 100644
--- a/source/fe/fe_poly.cc
+++ b/source/fe/fe_poly.cc
@@ -70,12 +70,7 @@ FE_Poly<TensorProductPolynomials<1>, 1, 2>::fill_fe_values(
                           mapping_internal,
                           make_array_view(output_data.shape_hessians, k));
 
-      for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
-        for (unsigned int i = 0; i < quadrature.size(); ++i)
-          for (unsigned int j = 0; j < 2; ++j)
-            output_data.shape_hessians[k][i] -=
-              mapping_data.jacobian_pushed_forward_grads[i][j] *
-              output_data.shape_gradients[k][i][j];
+      correct_hessians(output_data, mapping_data, quadrature.size());
     }
 
   if (fe_data.update_each & update_3rd_derivatives &&
@@ -139,12 +134,7 @@ FE_Poly<TensorProductPolynomials<2>, 2, 3>::fill_fe_values(
                           mapping_internal,
                           make_array_view(output_data.shape_hessians, k));
 
-      for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
-        for (unsigned int i = 0; i < quadrature.size(); ++i)
-          for (unsigned int j = 0; j < 3; ++j)
-            output_data.shape_hessians[k][i] -=
-              mapping_data.jacobian_pushed_forward_grads[i][j] *
-              output_data.shape_gradients[k][i][j];
+      correct_hessians(output_data, mapping_data, quadrature.size());
     }
 
   if (fe_data.update_each & update_3rd_derivatives &&