From e78adaa20c8b7dc447b32c7addc38c44ebab0c11 Mon Sep 17 00:00:00 2001
From: Maien Hamed <tomaien@hotmail.com>
Date: Thu, 13 Aug 2015 12:02:29 +0200
Subject: [PATCH] Adding transform of Tensor<3> and DerivativeForm<2> to
 mapping classes

---
 include/deal.II/fe/mapping.h           |  92 ++++++++++-
 include/deal.II/fe/mapping_cartesian.h |  15 +-
 include/deal.II/fe/mapping_fe_field.h  |  16 +-
 include/deal.II/fe/mapping_q1.h        |  23 ++-
 source/fe/fe_poly_tensor.cc            |   4 +-
 source/fe/mapping_cartesian.cc         | 123 ++++++++++++++
 source/fe/mapping_fe_field.cc          |  73 ++++++++-
 source/fe/mapping_q1.cc                | 212 ++++++++++++++++++++++---
 8 files changed, 526 insertions(+), 32 deletions(-)

diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h
index 1778ed3022..e4e00988ca 100644
--- a/include/deal.II/fe/mapping.h
+++ b/include/deal.II/fe/mapping.h
@@ -106,7 +106,27 @@ enum MappingType
   /**
    * The mapping used for BDM elements.
    */
-  mapping_bdm = mapping_raviart_thomas
+  mapping_bdm = mapping_raviart_thomas,
+
+  /**
+   * The mappings for 2-forms and third order tensors.
+   *
+   * These are mappings typpically applied to hessians transformed to the
+   * reference cell.
+   *
+   * Mapping of the hessian of a covariant vector field (see Mapping::transform() for details).
+   */
+  mapping_covariant_hessian,
+
+  /**
+  * Mapping of the hessian of a contravariant vector field (see Mapping::transform() for details).
+  */
+  mapping_contravariant_hessian,
+
+  /**
+  * Mapping of the hessian of a piola vector field (see Mapping::transform() for details).
+  */
+  mapping_piola_hessian
 };
 
 
@@ -859,6 +879,76 @@ public:
              const InternalDataBase &internal,
              const MappingType type) const = 0;
 
+  /**
+   * Transform a tensor field from the reference cell to the physical cell.
+   * This tensors are most of times the hessians in the reference cell of
+   * vector fields that have been pulled back from the physical cell.
+   *
+   * The mapping types currently implemented by derived classes are:
+   * <ul>
+   * <li> @p mapping_covariant_gradient: maps a field of forms on the reference cell to
+   * a field of forms on the physical cell. Mathematically, it is the pull
+   * back of the differential form
+   * @f[
+   * \mathbf T_{ijk}(\mathbf x) = \hat{\mathbf  T}_{iJK}(\hat{\mathbf  x}) J_{jJ}^{\dagger} J_{kK}^{\dagger}@f],
+   *
+   * where @f[ J^{\dagger} = J(\hat{\mathbf  x})(J(\hat{\mathbf  x})^{T} J(\hat{\mathbf  x}))^{-1}.
+   * @f]
+   *
+   * Hessians of spacedim-vector valued differentiable functions are
+   * transformed this way (After subtraction of the product of the
+   * derivative with the Jacobian gradient).
+   *
+   * In the case when dim=spacedim the previous formula reduces to
+   * @f[J^{\dagger} = J^{-1}@f]
+   */
+  virtual
+  void
+  transform (const VectorSlice<const std::vector< DerivativeForm<2, dim, spacedim> > > input,
+             VectorSlice<std::vector<Tensor<3,spacedim> > >             output,
+             const InternalDataBase &internal,
+             const MappingType type) const = 0;
+
+  /**
+   * Transform a field of 3-differential forms from the reference cell to the
+   * physical cell.  It is useful to think of $\mathbf{T}_{ijk} = D^2_{jk} \mathbf u_i$ and
+   * $\mathbf{\hat T}_{IJK} = \hat D^2_{JK} \mathbf{\hat u}_I$, with $\mathbf u_i$ a vector
+   * field.
+   *
+   * The mapping types currently implemented by derived classes are:
+   * <ul>
+   * <li> @p mapping_contravariant_hessian: it assumes $\mathbf u_i(\mathbf x)
+   * = J_{iI} \hat{\mathbf  u}_I$ so that
+   * @f[
+   * \mathbf T_{ijk}(\mathbf x) =
+   * J_{iI}(\hat{\mathbf  x}) \hat{\mathbf  T}_{IJK}(\hat{\mathbf  x})
+   * J_{jJ}(\hat{\mathbf  x})^{-1} J_{kK}(\hat{\mathbf  x})^{-1}.
+   * @f]
+   * <li> @p mapping_covariant_hessian: it assumes $\mathbf u_i(\mathbf x) =
+   * J_{iI}^{-T} \hat{\mathbf  u}_I$ so that
+   * @f[
+   * \mathbf T_{ijk}(\mathbf x) =
+   * J_iI(\hat{\mathbf  x})^{-1} \hat{\mathbf  T}_{IJK}(\hat{\mathbf  x})
+   * J_{jJ}(\hat{\mathbf  x})^{-1} J_{kK}(\hat{\mathbf  x})^{-1}.
+   * @f]
+   * <li> @p mapping_piola_hessian: it assumes $\mathbf u_i(\mathbf x) =
+   * \frac{1}{\text{det}\;J(\mathbf x)} J_{iI}(\mathbf x) \hat{\mathbf  u}(\mathbf x)$
+   * so that
+   * @f[
+   * \mathbf T_{ijk}(\mathbf x) =
+   * \frac{1}{\text{det}\;J(\mathbf x)}
+   * J_{iI}(\hat{\mathbf  x}) \hat{\mathbf  T}_{IJK}(\hat{\mathbf  x})
+   * J_{jJ}(\hat{\mathbf  x})^{-1} J_{kK}(\hat{\mathbf  x})^{-1}.
+   * @f]
+   * </ul>
+   */
+  virtual
+  void
+  transform (const VectorSlice<const std::vector<Tensor<3, dim> > >     input,
+             VectorSlice<std::vector<Tensor<3,spacedim> > >             output,
+             const InternalDataBase &internal,
+             const MappingType type) const = 0;
+
   /**
    * @}
    */
diff --git a/include/deal.II/fe/mapping_cartesian.h b/include/deal.II/fe/mapping_cartesian.h
index c003cf1193..1ec35ff102 100644
--- a/include/deal.II/fe/mapping_cartesian.h
+++ b/include/deal.II/fe/mapping_cartesian.h
@@ -119,7 +119,6 @@ public:
              const typename Mapping<dim,spacedim>::InternalDataBase &internal,
              const MappingType type) const;
 
-
   virtual
   void
   transform (const VectorSlice<const std::vector<Tensor<2, dim> > >     input,
@@ -127,6 +126,20 @@ public:
              const typename Mapping<dim,spacedim>::InternalDataBase &internal,
              const MappingType type) const;
 
+  virtual
+  void
+  transform (const VectorSlice<const std::vector< DerivativeForm<2, dim, spacedim> > > input,
+             VectorSlice<std::vector<Tensor<3,spacedim> > >             output,
+             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+             const MappingType type) const;
+
+  virtual
+  void
+  transform (const VectorSlice<const std::vector<Tensor<3, dim> > >     input,
+             VectorSlice<std::vector<Tensor<3,spacedim> > >             output,
+             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+             const MappingType type) const;
+
   virtual Point<spacedim>
   transform_unit_to_real_cell (
     const typename Triangulation<dim,spacedim>::cell_iterator &cell,
diff --git a/include/deal.II/fe/mapping_fe_field.h b/include/deal.II/fe/mapping_fe_field.h
index 341adbc4dc..2e64529b34 100644
--- a/include/deal.II/fe/mapping_fe_field.h
+++ b/include/deal.II/fe/mapping_fe_field.h
@@ -170,6 +170,20 @@ public:
              const typename Mapping<dim,spacedim>::InternalDataBase &internal,
              const MappingType type) const;
 
+  virtual
+  void
+  transform (const VectorSlice<const std::vector< DerivativeForm<2, dim, spacedim> > > input,
+             VectorSlice<std::vector<Tensor<3,spacedim> > >             output,
+             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+             const MappingType type) const;
+
+  virtual
+  void
+  transform (const VectorSlice<const std::vector<Tensor<3, dim> > >     input,
+             VectorSlice<std::vector<Tensor<3,spacedim> > >             output,
+             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+             const MappingType type) const;
+
 
 
   /**
@@ -462,7 +476,7 @@ protected:
   void
   transform_differential_forms(
     const VectorSlice<const std::vector<DerivativeForm<rank, dim,spacedim> > >    input,
-    VectorSlice<std::vector<DerivativeForm<rank, spacedim,spacedim> > > output,
+    VectorSlice<std::vector<Tensor<rank+1, spacedim> > > output,
     const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
     const MappingType mapping_type) const;
 
diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h
index 7f050b1ccd..0ac0ab8197 100644
--- a/include/deal.II/fe/mapping_q1.h
+++ b/include/deal.II/fe/mapping_q1.h
@@ -107,6 +107,20 @@ public:
              const typename Mapping<dim,spacedim>::InternalDataBase &internal,
              const MappingType type) const;
 
+  virtual
+  void
+  transform (const VectorSlice<const std::vector< DerivativeForm<2, dim, spacedim> > > input,
+             VectorSlice<std::vector<Tensor<3,spacedim> > >             output,
+             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+             const MappingType type) const;
+
+  virtual
+  void
+  transform (const VectorSlice<const std::vector<Tensor<3, dim> > >     input,
+             VectorSlice<std::vector<Tensor<3,spacedim> > >             output,
+             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+             const MappingType type) const;
+
 
 protected:
   /**
@@ -129,6 +143,13 @@ protected:
                       VectorSlice<      std::vector<Tensor<rank,spacedim> > > output,
                       const typename Mapping<dim,spacedim>::InternalDataBase &internal,
                       const MappingType type) const;
+
+  void
+  transform_hessians(const VectorSlice<const std::vector<Tensor<3,dim> > > input,
+                     VectorSlice<std::vector<Tensor<3,spacedim> > > output,
+                     const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+                     const MappingType mapping_type) const;
+
   /**
    * see doc in transform_fields
    */
@@ -136,7 +157,7 @@ protected:
   void
   transform_differential_forms(
     const VectorSlice<const std::vector<DerivativeForm<rank, dim, spacedim> > >    input,
-    VectorSlice<std::vector<DerivativeForm<rank, spacedim, spacedim> > > output,
+    VectorSlice<std::vector<Tensor<rank+1, spacedim> > > output,
     const typename Mapping<dim,spacedim>::InternalDataBase &internal,
     const MappingType type) const;
 
diff --git a/source/fe/fe_poly_tensor.cc b/source/fe/fe_poly_tensor.cc
index fc9e1e752a..529710f9e5 100644
--- a/source/fe/fe_poly_tensor.cc
+++ b/source/fe/fe_poly_tensor.cc
@@ -430,7 +430,7 @@ fill_fe_values (const Mapping<dim,spacedim>                      &mapping,
               break;
             }
             case mapping_raviart_thomas:
-            case mapping_piola_gradient:
+            case mapping_piola:
             {
               for (unsigned int k=0; k<n_q_points; ++k)
                 fe_data.untransformed_shape_grads[k] = fe_data.shape_grads[i][k];
@@ -639,7 +639,7 @@ fill_fe_face_values (const Mapping<dim,spacedim>                   &mapping,
               break;
             }
             case mapping_raviart_thomas:
-            case mapping_piola_gradient:
+            case mapping_piola:
             {
               for (unsigned int k=0; k<n_q_points; ++k)
                 fe_data.untransformed_shape_grads[k+offset] = fe_data.shape_grads[i][k+offset];
diff --git a/source/fe/mapping_cartesian.cc b/source/fe/mapping_cartesian.cc
index 8aca482fe4..cc0f59c8b2 100644
--- a/source/fe/mapping_cartesian.cc
+++ b/source/fe/mapping_cartesian.cc
@@ -779,6 +779,7 @@ MappingCartesian<dim,spacedim>::transform (
                                 / data.volume_element;
       return;
     }
+
     case mapping_piola_gradient:
     {
       Assert (data.update_each & update_contravariant_transformation,
@@ -801,6 +802,128 @@ MappingCartesian<dim,spacedim>::transform (
 }
 
 
+template<int dim, int spacedim>
+void
+MappingCartesian<dim,spacedim>::transform (
+  const VectorSlice<const std::vector< DerivativeForm<2, dim, spacedim> > > input,
+  VectorSlice<std::vector<Tensor<3,spacedim> > >             output,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+  const MappingType mapping_type) const
+{
+
+  AssertDimension (input.size(), output.size());
+  Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
+          ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(mapping_data);
+
+  switch (mapping_type)
+    {
+    case mapping_covariant_gradient:
+    {
+      Assert (data.update_each & update_contravariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+
+      for (unsigned int q=0; q<output.size(); ++q)
+        for (unsigned int i=0; i<spacedim; ++i)
+          for (unsigned int j=0; j<spacedim; ++j)
+            for (unsigned int k=0; k<spacedim; ++k)
+              {
+
+                output[q][i][j][k] = input[q][i][j][k] / data.length[j] / data.length[k];
+
+              }
+    }
+    default:
+      Assert(false, ExcNotImplemented());
+    }
+
+}
+
+template<int dim, int spacedim>
+void
+MappingCartesian<dim,spacedim>::transform (
+  const VectorSlice<const std::vector< Tensor<3,dim> > > input,
+  VectorSlice<std::vector<Tensor<3,spacedim> > >             output,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+  const MappingType mapping_type) const
+{
+
+  AssertDimension (input.size(), output.size());
+  Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
+          ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(mapping_data);
+
+  switch (mapping_type)
+    {
+    case mapping_contravariant_hessian:
+    {
+      Assert (data.update_each & update_covariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+      Assert (data.update_each & update_contravariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+
+      for (unsigned int q=0; q<output.size(); ++q)
+        for (unsigned int i=0; i<spacedim; ++i)
+          for (unsigned int j=0; j<spacedim; ++j)
+            for (unsigned int k=0; k<spacedim; ++k)
+              {
+                output[q][i][j][k] =    input[q][i][j][k]
+                                        * data.length[i]
+                                        / data.length[j]
+                                        / data.length[k];
+              }
+      return;
+    }
+
+    case mapping_covariant_hessian:
+    {
+      Assert (data.update_each & update_covariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+
+      for (unsigned int q=0; q<output.size(); ++q)
+        for (unsigned int i=0; i<spacedim; ++i)
+          for (unsigned int j=0; j<spacedim; ++j)
+            for (unsigned int k=0; k<spacedim; ++k)
+              {
+                output[q][i][j][k] =    input[q][i][j][k]
+                                        / data.length[i]
+                                        / data.length[j]
+                                        / data.length[k];
+              }
+
+      return;
+    }
+
+    case mapping_piola_hessian:
+    {
+      Assert (data.update_each & update_covariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+      Assert (data.update_each & update_contravariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+      Assert (data.update_each & update_volume_elements,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_volume_elements"));
+
+      for (unsigned int q=0; q<output.size(); ++q)
+        for (unsigned int i=0; i<spacedim; ++i)
+          for (unsigned int j=0; j<spacedim; ++j)
+            for (unsigned int k=0; k<spacedim; ++k)
+              {
+                output[q][i][j][k] =    input[q][i][j][k]
+                                        * data.length[i]
+                                        / data.volume_element
+                                        / data.length[j]
+                                        / data.length[k];
+              }
+
+      return;
+    }
+
+    default:
+      Assert(false, ExcNotImplemented());
+    }
+}
+
+
 template <int dim, int spacedim>
 Point<spacedim>
 MappingCartesian<dim, spacedim>::transform_unit_to_real_cell (
diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc
index a136a0be41..61c514ecb7 100644
--- a/source/fe/mapping_fe_field.cc
+++ b/source/fe/mapping_fe_field.cc
@@ -911,13 +911,8 @@ MappingFEField<dim,spacedim,VECTOR,DH>::transform (
 {
   AssertDimension (input.size(), output.size());
 
-  std::vector<DerivativeForm<1, spacedim,spacedim> > aux_output1(output.size());
-  VectorSlice< std::vector<DerivativeForm<1, spacedim,spacedim> > >  aux_output( aux_output1);
+  transform_differential_forms(input, output, mapping_data, mapping_type);
 
-  transform_differential_forms(input, aux_output, mapping_data, mapping_type);
-
-  for (unsigned int i=0; i<output.size(); i++)
-    output[i] = aux_output[i];
 }
 
 
@@ -937,6 +932,70 @@ void MappingFEField<dim,spacedim,VECTOR,DH>::transform
 }
 
 
+
+template<int dim, int spacedim, class VECTOR, class DH>
+void MappingFEField<dim,spacedim,VECTOR,DH>::transform (
+  const VectorSlice<const std::vector< DerivativeForm<2, dim, spacedim> > > input,
+  VectorSlice<std::vector<Tensor<3,spacedim> > >             output,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+  const MappingType mapping_type) const
+{
+
+  AssertDimension (input.size(), output.size());
+  Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
+          ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(mapping_data);
+
+  switch (mapping_type)
+    {
+    case mapping_covariant_gradient:
+    {
+      Assert (data.update_each & update_contravariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+
+      for (unsigned int q=0; q<output.size(); ++q)
+        for (unsigned int i=0; i<spacedim; ++i)
+          for (unsigned int j=0; j<spacedim; ++j)
+            for (unsigned int k=0; k<spacedim; ++k)
+              {
+                output[q][i][j][k] = data.covariant[q][j][0]
+                                     * data.covariant[q][k][0]
+                                     * input[q][i][0][0];
+                for (unsigned int J=0; J<dim; ++J)
+                  {
+                    const unsigned int K0 = (0==J)? 1 : 0;
+                    for (unsigned int K=K0; K<dim; ++K)
+                      output[q][i][j][k] += data.covariant[q][j][J]
+                                            * data.covariant[q][k][K]
+                                            * input[q][i][J][K];
+                  }
+
+              }
+    }
+    default:
+      Assert(false, ExcNotImplemented());
+    }
+
+}
+
+template<int dim, int spacedim, class VECTOR, class DH>
+void MappingFEField<dim,spacedim,VECTOR,DH>::transform (
+  const VectorSlice<const std::vector< Tensor<3,dim> > > input,
+  VectorSlice<std::vector<Tensor<3,spacedim> > >             output,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+  const MappingType /*mapping_type*/) const
+{
+
+  (void)input;
+  (void)output;
+  (void)mapping_data;
+  AssertDimension (input.size(), output.size());
+
+  AssertThrow(false, ExcNotImplemented());
+
+}
+
+
 template<int dim, int spacedim, class VECTOR, class DH>
 template < int rank >
 void MappingFEField<dim,spacedim,VECTOR,DH>::transform_fields(
@@ -1003,7 +1062,7 @@ template<int dim, int spacedim, class VECTOR, class DH>
 template < int rank >
 void MappingFEField<dim,spacedim,VECTOR,DH>::transform_differential_forms(
   const VectorSlice<const std::vector<DerivativeForm<rank, dim,spacedim> > >    input,
-  VectorSlice<std::vector<DerivativeForm<rank, spacedim,spacedim> > > output,
+  VectorSlice<std::vector<Tensor<rank+1, spacedim> > > output,
   const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
   const MappingType mapping_type) const
 {
diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc
index fb913e9dc9..52eeb241d6 100644
--- a/source/fe/mapping_q1.cc
+++ b/source/fe/mapping_q1.cc
@@ -1318,13 +1318,8 @@ MappingQ1<dim,spacedim>::transform (
   const MappingType mapping_type) const
 {
 
-  std::vector<DerivativeForm<1, spacedim,spacedim> > aux_output1(output.size());
-  VectorSlice< std::vector<DerivativeForm<1, spacedim,spacedim> > >  aux_output( aux_output1);
+  transform_differential_forms(input, output, mapping_data, mapping_type);
 
-  transform_differential_forms(input, aux_output, mapping_data, mapping_type);
-
-  for (unsigned int i=0; i<output.size(); i++)
-    output[i] = aux_output[i];
 }
 
 
@@ -1337,6 +1332,7 @@ MappingQ1<dim,spacedim>::transform (
   const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
   const MappingType mapping_type) const
 {
+
   switch (mapping_type)
     {
     case mapping_contravariant:
@@ -1351,9 +1347,78 @@ MappingQ1<dim,spacedim>::transform (
     default:
       Assert(false, ExcNotImplemented());
     }
+
+}
+
+
+
+template<int dim, int spacedim>
+void
+MappingQ1<dim,spacedim>::transform (
+  const VectorSlice<const std::vector< DerivativeForm<2, dim, spacedim> > > input,
+  VectorSlice<std::vector<Tensor<3,spacedim> > >             output,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+  const MappingType mapping_type) const
+{
+
+  AssertDimension (input.size(), output.size());
+  Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
+          ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(mapping_data);
+
+  switch (mapping_type)
+    {
+    case mapping_covariant_gradient:
+    {
+      Assert (data.update_each & update_contravariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+
+      for (unsigned int q=0; q<output.size(); ++q)
+        for (unsigned int i=0; i<spacedim; ++i)
+          for (unsigned int j=0; j<spacedim; ++j)
+            for (unsigned int k=0; k<spacedim; ++k)
+              {
+                output[q][i][j][k] = data.covariant[q][j][0]
+                                     * data.covariant[q][k][0]
+                                     * input[q][i][0][0];
+                for (unsigned int J=0; J<dim; ++J)
+                  {
+                    const unsigned int K0 = (0==J)? 1 : 0;
+                    for (unsigned int K=K0; K<dim; ++K)
+                      output[q][i][j][k] += data.covariant[q][j][J]
+                                            * data.covariant[q][k][K]
+                                            * input[q][i][J][K];
+                  }
+
+              }
+    }
+    default:
+      Assert(false, ExcNotImplemented());
+    }
+
 }
 
+template<int dim, int spacedim>
+void
+MappingQ1<dim,spacedim>::transform (
+  const VectorSlice<const std::vector< Tensor<3,dim> > > input,
+  VectorSlice<std::vector<Tensor<3,spacedim> > >             output,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+  const MappingType mapping_type) const
+{
+
+  switch (mapping_type)
+    {
+    case mapping_piola_hessian:
+    case mapping_contravariant_hessian:
+    case mapping_covariant_hessian:
+      transform_hessians(input, output, mapping_data, mapping_type);
+      return;
+    default:
+      Assert(false, ExcNotImplemented());
+    }
 
+}
 
 template<int dim, int spacedim>
 template < int rank >
@@ -1398,15 +1463,13 @@ void MappingQ1<dim,spacedim>::transform_fields(
         }
       return;
     }
-
-
     //We still allow this operation as in the
     //reference cell Derivatives are Tensor
     //rather than DerivativeForm
     case mapping_covariant:
     {
       Assert (data.update_each & update_contravariant_transformation,
-              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
 
       for (unsigned int i=0; i<output.size(); ++i)
         output[i] = apply_transformation(data.covariant[i], input[i]);
@@ -1457,8 +1520,6 @@ void MappingQ1<dim,spacedim>::transform_gradients(
     {
       Assert (data.update_each & update_covariant_transformation,
               typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
-      Assert (data.update_each & update_contravariant_transformation,
-              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
       Assert (rank==2, ExcMessage("Only for rank 2"));
 
       for (unsigned int i=0; i<output.size(); ++i)
@@ -1466,15 +1527,11 @@ void MappingQ1<dim,spacedim>::transform_gradients(
           DerivativeForm<1,spacedim,dim> A =
             apply_transformation(data.covariant[i], transpose(input[i]) );
           output[i] = apply_transformation(data.covariant[i], A.transpose() );
-
         }
 
       return;
     }
 
-
-
-
     case mapping_piola_gradient:
     {
       Assert (data.update_each & update_covariant_transformation,
@@ -1507,11 +1564,130 @@ void MappingQ1<dim,spacedim>::transform_gradients(
 
 
 
+template<int dim, int spacedim>
+void MappingQ1<dim,spacedim>::transform_hessians(
+  const VectorSlice<const std::vector<Tensor<3,dim> > > input,
+  VectorSlice<std::vector<Tensor<3,spacedim> > > output,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+  const MappingType mapping_type) const
+{
+  AssertDimension (input.size(), output.size());
+  Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
+          ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(mapping_data);
+
+  switch (mapping_type)
+    {
+    case mapping_contravariant_hessian:
+    {
+      Assert (data.update_each & update_covariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+      Assert (data.update_each & update_contravariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+
+      for (unsigned int q=0; q<output.size(); ++q)
+        for (unsigned int i=0; i<spacedim; ++i)
+          for (unsigned int j=0; j<spacedim; ++j)
+            for (unsigned int k=0; k<spacedim; ++k)
+              {
+                output[q][i][j][k] =    data.contravariant[q][i][0]
+                                        * data.covariant[q][j][0]
+                                        * data.covariant[q][k][0]
+                                        * input[q][0][0][0];
+                for (unsigned int I=0; I<dim; ++I)
+                  for (unsigned int J=0; J<dim; ++J)
+                    {
+                      const unsigned int K0 = (0==(I+J))? 1 : 0;
+                      for (unsigned int K=K0; K<dim; ++K)
+                        output[q][i][j][k] +=    data.contravariant[q][i][I]
+                                                 * data.covariant[q][j][J]
+                                                 * data.covariant[q][k][K]
+                                                 * input[q][I][J][K];
+                    }
+
+              }
+      return;
+    }
+
+    case mapping_covariant_hessian:
+    {
+      Assert (data.update_each & update_covariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+
+      for (unsigned int q=0; q<output.size(); ++q)
+        for (unsigned int i=0; i<spacedim; ++i)
+          for (unsigned int j=0; j<spacedim; ++j)
+            for (unsigned int k=0; k<spacedim; ++k)
+              {
+                output[q][i][j][k] =    data.covariant[q][i][0]
+                                        * data.covariant[q][j][0]
+                                        * data.covariant[q][k][0]
+                                        * input[q][0][0][0];
+                for (unsigned int I=0; I<dim; ++I)
+                  for (unsigned int J=0; J<dim; ++J)
+                    {
+                      const unsigned int K0 = (0==(I+J))? 1 : 0;
+                      for (unsigned int K=K0; K<dim; ++K)
+                        output[q][i][j][k] +=   data.covariant[q][i][I]
+                                                * data.covariant[q][j][J]
+                                                * data.covariant[q][k][K]
+                                                * input[q][I][J][K];
+                    }
+
+              }
+
+      return;
+    }
+
+    case mapping_piola_hessian:
+    {
+      Assert (data.update_each & update_covariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+      Assert (data.update_each & update_contravariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+      Assert (data.update_each & update_volume_elements,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_volume_elements"));
+
+      for (unsigned int q=0; q<output.size(); ++q)
+        for (unsigned int i=0; i<spacedim; ++i)
+          for (unsigned int j=0; j<spacedim; ++j)
+            for (unsigned int k=0; k<spacedim; ++k)
+              {
+                output[q][i][j][k] =    data.contravariant[q][i][0]
+                                        / data.volume_elements[q]
+                                        * data.covariant[q][j][0]
+                                        * data.covariant[q][k][0]
+                                        * input[q][0][0][0];
+                for (unsigned int I=0; I<dim; ++I)
+                  for (unsigned int J=0; J<dim; ++J)
+                    {
+                      const unsigned int K0 = (0==(I+J))? 1 : 0;
+                      for (unsigned int K=K0; K<dim; ++K)
+                        output[q][i][j][k] +=    data.contravariant[q][i][I]
+                                                 / data.volume_elements[q]
+                                                 * data.covariant[q][j][J]
+                                                 * data.covariant[q][k][K]
+                                                 * input[q][I][J][K];
+                    }
+
+              }
+
+      return;
+    }
+
+    default:
+      Assert(false, ExcNotImplemented());
+    }
+}
+
+
+
+
 template<int dim, int spacedim>
 template < int rank >
 void MappingQ1<dim,spacedim>::transform_differential_forms(
   const VectorSlice<const std::vector<DerivativeForm<rank, dim,spacedim> > >    input,
-  VectorSlice<std::vector<DerivativeForm<rank, spacedim,spacedim> > > output,
+  VectorSlice<std::vector<Tensor<rank+1, spacedim> > > output,
   const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
   const MappingType mapping_type) const
 {
@@ -1526,7 +1702,7 @@ void MappingQ1<dim,spacedim>::transform_differential_forms(
     case mapping_covariant:
     {
       Assert (data.update_each & update_contravariant_transformation,
-              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
 
       for (unsigned int i=0; i<output.size(); ++i)
         output[i] = apply_transformation(data.covariant[i], input[i]);
@@ -1542,8 +1718,6 @@ void MappingQ1<dim,spacedim>::transform_differential_forms(
 
 
 
-
-
 template<int dim, int spacedim>
 Point<spacedim>
 MappingQ1<dim,spacedim>::transform_unit_to_real_cell (
-- 
2.39.5