From c852c94c18c36dbf7f3ad4dac26c52d8d4ff4ace Mon Sep 17 00:00:00 2001
From: Luca Heltai <luca.heltai@sissa.it>
Date: Mon, 28 Mar 2016 12:52:55 +0200
Subject: [PATCH] Fixed JxW values on faces.

---
 source/fe/mapping_manifold.cc       | 667 +---------------------------
 tests/fe/mapping_manifold_07.cc     |  54 ++-
 tests/fe/mapping_manifold_07.output | 214 +--------
 3 files changed, 51 insertions(+), 884 deletions(-)

diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc
index 0357c2e8b4..a921dce6d5 100644
--- a/source/fe/mapping_manifold.cc
+++ b/source/fe/mapping_manifold.cc
@@ -96,21 +96,6 @@ initialize (const UpdateFlags      update_flags,
 
   if (this->update_each & update_volume_elements)
     volume_elements.resize(n_original_q_points);
-
-  // if (this->update_each &
-  //     (update_jacobian_grads | update_jacobian_pushed_forward_grads) )
-  //   shape_second_derivatives.resize(n_shape_functions * n_q_points);
-
-  // if (this->update_each &
-  //     (update_jacobian_2nd_derivatives | update_jacobian_pushed_forward_2nd_derivatives) )
-  //   shape_third_derivatives.resize(n_shape_functions * n_q_points);
-
-  // if (this->update_each &
-  //     (update_jacobian_3rd_derivatives | update_jacobian_pushed_forward_3rd_derivatives) )
-  //   shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
-
-  // // now also fill the various fields with their correct values
-  // compute_shape_function_values (q.get_points());
 }
 
 
@@ -483,440 +468,11 @@ namespace internal
             }
         }
     }
-
-    /**
-     * Update the Hessian of the transformation from unit to real cell, the
-     * Jacobian gradients.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim>
-    void
-    maybe_update_jacobian_grads (const CellSimilarity::Similarity                                   cell_similarity,
-                                 const typename QProjector<dim>::DataSetDescriptor                  data_set,
-                                 const typename dealii::MappingManifold<dim,spacedim>::InternalData      &data,
-                                 std::vector<DerivativeForm<2,dim,spacedim> >                      &jacobian_grads)
-    {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_grads)
-        {
-          const unsigned int n_q_points = jacobian_grads.size();
-
-          if (cell_similarity != CellSimilarity::translation)
-            {
-              for (unsigned int point=0; point<n_q_points; ++point)
-                {
-                  const Tensor<2,dim> *second =
-                    &data.second_derivative(point+data_set, 0);
-                  double result [spacedim][dim][dim];
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        result[i][j][l] = (second[0][j][l] *
-                                           data.mapping_support_points[0][i]);
-                  for (unsigned int k=1; k<data.n_shape_functions; ++k)
-                    for (unsigned int i=0; i<spacedim; ++i)
-                      for (unsigned int j=0; j<dim; ++j)
-                        for (unsigned int l=0; l<dim; ++l)
-                          result[i][j][l]
-                          += (second[k][j][l]
-                              *
-                              data.mapping_support_points[k][i]);
-
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        jacobian_grads[point][i][j][l] = result[i][j][l];
-                }
-            }
-        }
-    }
-
-    /**
-     * Update the Hessian of the transformation from unit to real cell, the
-     * Jacobian gradients, pushed forward to the real cell coordinates.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim>
-    void
-    maybe_update_jacobian_pushed_forward_grads (const CellSimilarity::Similarity                                   cell_similarity,
-                                                const typename QProjector<dim>::DataSetDescriptor                  data_set,
-                                                const typename dealii::MappingManifold<dim,spacedim>::InternalData      &data,
-                                                std::vector<Tensor<3,spacedim> >                      &jacobian_pushed_forward_grads)
-    {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_pushed_forward_grads)
-        {
-          const unsigned int n_q_points = jacobian_pushed_forward_grads.size();
-
-          if (cell_similarity != CellSimilarity::translation)
-            {
-              double tmp[spacedim][spacedim][spacedim];
-              for (unsigned int point=0; point<n_q_points; ++point)
-                {
-                  const Tensor<2,dim> *second =
-                    &data.second_derivative(point+data_set, 0);
-                  double result [spacedim][dim][dim];
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        result[i][j][l] = (second[0][j][l] *
-                                           data.mapping_support_points[0][i]);
-                  for (unsigned int k=1; k<data.n_shape_functions; ++k)
-                    for (unsigned int i=0; i<spacedim; ++i)
-                      for (unsigned int j=0; j<dim; ++j)
-                        for (unsigned int l=0; l<dim; ++l)
-                          result[i][j][l]
-                          += (second[k][j][l]
-                              *
-                              data.mapping_support_points[k][i]);
-
-                  // first push forward the j-components
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<spacedim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        {
-                          tmp[i][j][l] = result[i][0][l] *
-                                         data.covariant[point][j][0];
-                          for (unsigned int jr=1; jr<dim; ++jr)
-                            {
-                              tmp[i][j][l] += result[i][jr][l] *
-                                              data.covariant[point][j][jr];
-                            }
-                        }
-
-                  // now, pushing forward the l-components
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<spacedim; ++j)
-                      for (unsigned int l=0; l<spacedim; ++l)
-                        {
-                          jacobian_pushed_forward_grads[point][i][j][l] = tmp[i][j][0] *
-                                                                          data.covariant[point][l][0];
-                          for (unsigned int lr=1; lr<dim; ++lr)
-                            {
-                              jacobian_pushed_forward_grads[point][i][j][l] += tmp[i][j][lr] *
-                                                                               data.covariant[point][l][lr];
-                            }
-
-                        }
-                }
-            }
-        }
-    }
-
-    /**
-     * Update the third derivatives of the transformation from unit to real cell, the
-     * Jacobian hessians.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim>
-    void
-    maybe_update_jacobian_2nd_derivatives (const CellSimilarity::Similarity                              cell_similarity,
-                                           const typename QProjector<dim>::DataSetDescriptor             data_set,
-                                           const typename dealii::MappingManifold<dim,spacedim>::InternalData &data,
-                                           std::vector<DerivativeForm<3,dim,spacedim> >                 &jacobian_2nd_derivatives)
-    {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_2nd_derivatives)
-        {
-          const unsigned int n_q_points = jacobian_2nd_derivatives.size();
-
-          if (cell_similarity != CellSimilarity::translation)
-            {
-              for (unsigned int point=0; point<n_q_points; ++point)
-                {
-                  const Tensor<3,dim> *third =
-                    &data.third_derivative(point+data_set, 0);
-                  double result [spacedim][dim][dim][dim];
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          result[i][j][l][m] = (third[0][j][l][m] *
-                                                data.mapping_support_points[0][i]);
-                  for (unsigned int k=1; k<data.n_shape_functions; ++k)
-                    for (unsigned int i=0; i<spacedim; ++i)
-                      for (unsigned int j=0; j<dim; ++j)
-                        for (unsigned int l=0; l<dim; ++l)
-                          for (unsigned int m=0; m<dim; ++m)
-                            result[i][j][l][m]
-                            += (third[k][j][l][m]
-                                *
-                                data.mapping_support_points[k][i]);
-
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          jacobian_2nd_derivatives[point][i][j][l][m] = result[i][j][l][m];
-                }
-            }
-        }
-    }
-
-    /**
-     * Update the Hessian of the Hessian of the transformation from unit
-     * to real cell, the Jacobian Hessian gradients, pushed forward to the
-     * real cell coordinates.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim>
-    void
-    maybe_update_jacobian_pushed_forward_2nd_derivatives (const CellSimilarity::Similarity                                   cell_similarity,
-                                                          const typename QProjector<dim>::DataSetDescriptor                  data_set,
-                                                          const typename dealii::MappingManifold<dim,spacedim>::InternalData      &data,
-                                                          std::vector<Tensor<4,spacedim> >                      &jacobian_pushed_forward_2nd_derivatives)
-    {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_pushed_forward_2nd_derivatives)
-        {
-          const unsigned int n_q_points = jacobian_pushed_forward_2nd_derivatives.size();
-
-          if (cell_similarity != CellSimilarity::translation)
-            {
-              double tmp[spacedim][spacedim][spacedim][spacedim];
-              for (unsigned int point=0; point<n_q_points; ++point)
-                {
-                  const Tensor<3,dim> *third =
-                    &data.third_derivative(point+data_set, 0);
-                  double result [spacedim][dim][dim][dim];
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          result[i][j][l][m] = (third[0][j][l][m] *
-                                                data.mapping_support_points[0][i]);
-                  for (unsigned int k=1; k<data.n_shape_functions; ++k)
-                    for (unsigned int i=0; i<spacedim; ++i)
-                      for (unsigned int j=0; j<dim; ++j)
-                        for (unsigned int l=0; l<dim; ++l)
-                          for (unsigned int m=0; m<dim; ++m)
-                            result[i][j][l][m]
-                            += (third[k][j][l][m]
-                                *
-                                data.mapping_support_points[k][i]);
-
-                  // push forward the j-coordinate
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<spacedim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          {
-                            jacobian_pushed_forward_2nd_derivatives[point][i][j][l][m]
-                              = result[i][0][l][m]*
-                                data.covariant[point][j][0];
-                            for (unsigned int jr=1; jr<dim; ++jr)
-                              jacobian_pushed_forward_2nd_derivatives[point][i][j][l][m]
-                              += result[i][jr][l][m]*
-                                 data.covariant[point][j][jr];
-                          }
-
-                  // push forward the l-coordinate
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<spacedim; ++j)
-                      for (unsigned int l=0; l<spacedim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          {
-                            tmp[i][j][l][m]
-                              = jacobian_pushed_forward_2nd_derivatives[point][i][j][0][m]*
-                                data.covariant[point][l][0];
-                            for (unsigned int lr=1; lr<dim; ++lr)
-                              tmp[i][j][l][m]
-                              += jacobian_pushed_forward_2nd_derivatives[point][i][j][lr][m]*
-                                 data.covariant[point][l][lr];
-                          }
-
-                  // push forward the m-coordinate
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<spacedim; ++j)
-                      for (unsigned int l=0; l<spacedim; ++l)
-                        for (unsigned int m=0; m<spacedim; ++m)
-                          {
-                            jacobian_pushed_forward_2nd_derivatives[point][i][j][l][m]
-                              = tmp[i][j][l][0]*
-                                data.covariant[point][m][0];
-                            for (unsigned int mr=1; mr<dim; ++mr)
-                              jacobian_pushed_forward_2nd_derivatives[point][i][j][l][m]
-                              += tmp[i][j][l][mr]*
-                                 data.covariant[point][m][mr];
-                          }
-                }
-            }
-        }
-    }
-
-    /**
-         * Update the fourth derivatives of the transformation from unit to real cell, the
-         * Jacobian hessian gradients.
-         *
-         * Skip the computation if possible as indicated by the first argument.
-         */
-    template <int dim, int spacedim>
-    void
-    maybe_update_jacobian_3rd_derivatives (const CellSimilarity::Similarity                              cell_similarity,
-                                           const typename QProjector<dim>::DataSetDescriptor             data_set,
-                                           const typename dealii::MappingManifold<dim,spacedim>::InternalData &data,
-                                           std::vector<DerivativeForm<4,dim,spacedim> >                 &jacobian_3rd_derivatives)
-    {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_3rd_derivatives)
-        {
-          const unsigned int n_q_points = jacobian_3rd_derivatives.size();
-
-          if (cell_similarity != CellSimilarity::translation)
-            {
-              for (unsigned int point=0; point<n_q_points; ++point)
-                {
-                  const Tensor<4,dim> *fourth =
-                    &data.fourth_derivative(point+data_set, 0);
-                  double result [spacedim][dim][dim][dim][dim];
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          for (unsigned int n=0; n<dim; ++n)
-                            result[i][j][l][m][n] = (fourth[0][j][l][m][n] *
-                                                     data.mapping_support_points[0][i]);
-                  for (unsigned int k=1; k<data.n_shape_functions; ++k)
-                    for (unsigned int i=0; i<spacedim; ++i)
-                      for (unsigned int j=0; j<dim; ++j)
-                        for (unsigned int l=0; l<dim; ++l)
-                          for (unsigned int m=0; m<dim; ++m)
-                            for (unsigned int n=0; n<dim; ++n)
-                              result[i][j][l][m][n]
-                              += (fourth[k][j][l][m][n]
-                                  *
-                                  data.mapping_support_points[k][i]);
-
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          for (unsigned int n=0; n<dim; ++n)
-                            jacobian_3rd_derivatives[point][i][j][l][m][n] = result[i][j][l][m][n];
-                }
-            }
-        }
-    }
-
-    /**
-     * Update the Hessian gradient of the transformation from unit to real cell, the
-     * Jacobian Hessians, pushed forward to the real cell coordinates.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim>
-    void
-    maybe_update_jacobian_pushed_forward_3rd_derivatives (const CellSimilarity::Similarity                                   cell_similarity,
-                                                          const typename QProjector<dim>::DataSetDescriptor                  data_set,
-                                                          const typename dealii::MappingManifold<dim,spacedim>::InternalData      &data,
-                                                          std::vector<Tensor<5,spacedim> >                      &jacobian_pushed_forward_3rd_derivatives)
-    {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_pushed_forward_3rd_derivatives)
-        {
-          const unsigned int n_q_points = jacobian_pushed_forward_3rd_derivatives.size();
-
-          if (cell_similarity != CellSimilarity::translation)
-            {
-              double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
-              for (unsigned int point=0; point<n_q_points; ++point)
-                {
-                  const Tensor<4,dim> *fourth =
-                    &data.fourth_derivative(point+data_set, 0);
-                  double result [spacedim][dim][dim][dim][dim];
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          for (unsigned int n=0; n<dim; ++n)
-                            result[i][j][l][m][n] = (fourth[0][j][l][m][n] *
-                                                     data.mapping_support_points[0][i]);
-                  for (unsigned int k=1; k<data.n_shape_functions; ++k)
-                    for (unsigned int i=0; i<spacedim; ++i)
-                      for (unsigned int j=0; j<dim; ++j)
-                        for (unsigned int l=0; l<dim; ++l)
-                          for (unsigned int m=0; m<dim; ++m)
-                            for (unsigned int n=0; n<dim; ++n)
-                              result[i][j][l][m][n]
-                              += (fourth[k][j][l][m][n]
-                                  *
-                                  data.mapping_support_points[k][i]);
-
-                  // push-forward the j-coordinate
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<spacedim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          for (unsigned int n=0; n<dim; ++n)
-                            {
-                              tmp[i][j][l][m][n] = result[i][0][l][m][n] *
-                                                   data.covariant[point][j][0];
-                              for (unsigned int jr=1; jr<dim; ++jr)
-                                tmp[i][j][l][m][n] += result[i][jr][l][m][n] *
-                                                      data.covariant[point][j][jr];
-                            }
-
-                  // push-forward the l-coordinate
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<spacedim; ++j)
-                      for (unsigned int l=0; l<spacedim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          for (unsigned int n=0; n<dim; ++n)
-                            {
-                              jacobian_pushed_forward_3rd_derivatives[point][i][j][l][m][n]
-                                = tmp[i][j][0][m][n] *
-                                  data.covariant[point][l][0];
-                              for (unsigned int lr=1; lr<dim; ++lr)
-                                jacobian_pushed_forward_3rd_derivatives[point][i][j][l][m][n]
-                                += tmp[i][j][lr][m][n] *
-                                   data.covariant[point][l][lr];
-                            }
-
-                  // push-forward the m-coordinate
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<spacedim; ++j)
-                      for (unsigned int l=0; l<spacedim; ++l)
-                        for (unsigned int m=0; m<spacedim; ++m)
-                          for (unsigned int n=0; n<dim; ++n)
-                            {
-                              tmp[i][j][l][m][n]
-                                = jacobian_pushed_forward_3rd_derivatives[point][i][j][l][0][n] *
-                                  data.covariant[point][m][0];
-                              for (unsigned int mr=1; mr<dim; ++mr)
-                                tmp[i][j][l][m][n]
-                                += jacobian_pushed_forward_3rd_derivatives[point][i][j][l][mr][n] *
-                                   data.covariant[point][m][mr];
-                            }
-
-                  // push-forward the n-coordinate
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<spacedim; ++j)
-                      for (unsigned int l=0; l<spacedim; ++l)
-                        for (unsigned int m=0; m<spacedim; ++m)
-                          for (unsigned int n=0; n<spacedim; ++n)
-                            {
-                              jacobian_pushed_forward_3rd_derivatives[point][i][j][l][m][n]
-                                = tmp[i][j][l][m][0] *
-                                  data.covariant[point][n][0];
-                              for (unsigned int nr=1; nr<dim; ++nr)
-                                jacobian_pushed_forward_3rd_derivatives[point][i][j][l][m][n]
-                                += tmp[i][j][l][m][nr] *
-                                   data.covariant[point][n][nr];
-                            }
-                }
-            }
-        }
-    }
   }
 }
 
 
 
-
 template<int dim, int spacedim>
 CellSimilarity::Similarity
 MappingManifold<dim,spacedim>::
@@ -1049,36 +605,6 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
           output_data.inverse_jacobians[point] = data.covariant[point].transpose();
     }
 
-//   internal::maybe_update_jacobian_grads<dim,spacedim> (cell_similarity,
-//                                                        QProjector<dim>::DataSetDescriptor::cell (),
-//                                                        data,
-//                                                        output_data.jacobian_grads);
-
-//   internal::maybe_update_jacobian_pushed_forward_grads<dim,spacedim> (cell_similarity,
-//       QProjector<dim>::DataSetDescriptor::cell (),
-//       data,
-//       output_data.jacobian_pushed_forward_grads);
-
-//   internal::maybe_update_jacobian_2nd_derivatives<dim,spacedim> (cell_similarity,
-//       QProjector<dim>::DataSetDescriptor::cell (),
-//       data,
-//       output_data.jacobian_2nd_derivatives);
-
-//   internal::maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,spacedim> (cell_similarity,
-//       QProjector<dim>::DataSetDescriptor::cell (),
-//       data,
-//       output_data.jacobian_pushed_forward_2nd_derivatives);
-
-//   internal::maybe_update_jacobian_3rd_derivatives<dim,spacedim> (cell_similarity,
-//       QProjector<dim>::DataSetDescriptor::cell (),
-//       data,
-//       output_data.jacobian_3rd_derivatives);
-
-//   internal::maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,spacedim> (cell_similarity,
-//       QProjector<dim>::DataSetDescriptor::cell (),
-//       data,
-//       output_data.jacobian_pushed_forward_3rd_derivatives);
-
   return cell_similarity;
 }
 
@@ -1262,31 +788,6 @@ namespace internal
       maybe_update_Jacobians<dim,spacedim> (data_set,
                                             data);
 
-      // maybe_update_jacobian_grads<dim,spacedim> (CellSimilarity::none,
-      //                                            data_set,
-      //                                            data,
-      //                                            output_data.jacobian_grads);
-      // maybe_update_jacobian_pushed_forward_grads<dim,spacedim> (CellSimilarity::none,
-      //                                                           data_set,
-      //                                                           data,
-      //                                                           output_data.jacobian_pushed_forward_grads);
-      // maybe_update_jacobian_2nd_derivatives<dim,spacedim> (CellSimilarity::none,
-      //                                                      data_set,
-      //                                                      data,
-      //                                                      output_data.jacobian_2nd_derivatives);
-      // maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,spacedim> (CellSimilarity::none,
-      //     data_set,
-      //     data,
-      //     output_data.jacobian_pushed_forward_2nd_derivatives);
-      // maybe_update_jacobian_3rd_derivatives<dim,spacedim> (CellSimilarity::none,
-      //                                                      data_set,
-      //                                                      data,
-      //                                                      output_data.jacobian_3rd_derivatives);
-      // maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,spacedim> (CellSimilarity::none,
-      //     data_set,
-      //     data,
-      //     output_data.jacobian_pushed_forward_3rd_derivatives);
-
       maybe_compute_face_data (mapping,
                                cell, face_no, subface_no, quadrature.size(),
                                quadrature.get_weights(), data,
@@ -1686,7 +1187,7 @@ transform (const ArrayView<const Tensor<1, dim> >                  &input,
            const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
            const ArrayView<Tensor<1, spacedim> >                   &output) const
 {
-  //  transform_fields(input, mapping_type, mapping_data, output);
+  transform_fields(input, mapping_type, mapping_data, output);
 }
 
 
@@ -1799,172 +1300,6 @@ transform (const ArrayView<const  Tensor<3,dim> >                  &input,
     }
 }
 
-
-
-namespace
-{
-  /**
-   * Ask the manifold descriptor to return intermediate points on lines or
-   * faces. The function needs to return one or multiple points (depending on
-   * the number of elements in the output vector @p points that lie inside a
-   * line, quad or hex). Whether it is a line, quad or hex doesn't really
-   * matter to this function but it can be inferred from the number of input
-   * points in the @p surrounding_points vector.
-   */
-  template<int dim, int spacedim>
-  void
-  get_intermediate_points (const Manifold<dim, spacedim> &manifold,
-                           const QGaussLobatto<1>        &line_support_points,
-                           const std::vector<Point<spacedim> > &surrounding_points,
-                           std::vector<Point<spacedim> > &points)
-  {
-    Assert(surrounding_points.size() >= 2, ExcMessage("At least 2 surrounding points are required"));
-    const unsigned int n=points.size();
-    Assert(n>0, ExcMessage("You can't ask for 0 intermediate points."));
-    std::vector<double> w(surrounding_points.size());
-
-    switch (surrounding_points.size())
-      {
-      case 2:
-      {
-        // If two points are passed, these are the two vertices, and
-        // we can only compute degree-1 intermediate points.
-        for (unsigned int i=0; i<n; ++i)
-          {
-            const double x = line_support_points.point(i+1)[0];
-            w[1] = x;
-            w[0] = (1-x);
-            Quadrature<spacedim> quadrature(surrounding_points, w);
-            points[i] = manifold.get_new_point(quadrature);
-          }
-        break;
-      }
-
-      case 4:
-      {
-        Assert(spacedim >= 2, ExcImpossibleInDim(spacedim));
-        const unsigned m=
-          static_cast<unsigned int>(std::sqrt(static_cast<double>(n)));
-        // is n a square number
-        Assert(m*m==n, ExcInternalError());
-
-        // If four points are passed, these are the two vertices, and
-        // we can only compute (degree-1)*(degree-1) intermediate
-        // points.
-        for (unsigned int i=0; i<m; ++i)
-          {
-            const double y=line_support_points.point(1+i)[0];
-            for (unsigned int j=0; j<m; ++j)
-              {
-                const double x=line_support_points.point(1+j)[0];
-
-                w[0] = (1-x)*(1-y);
-                w[1] =     x*(1-y);
-                w[2] = (1-x)*y    ;
-                w[3] =     x*y    ;
-                Quadrature<spacedim> quadrature(surrounding_points, w);
-                points[i*m+j]=manifold.get_new_point(quadrature);
-              }
-          }
-        break;
-      }
-
-      case 8:
-        Assert(false, ExcNotImplemented());
-        break;
-      default:
-        Assert(false, ExcInternalError());
-        break;
-      }
-  }
-
-
-
-
-  /**
-   * Ask the manifold descriptor to return intermediate points on the object
-   * pointed to by the TriaIterator @p iter. This function tries to be
-   * backward compatible with respect to the differences between
-   * Boundary<dim,spacedim> and Manifold<dim,spacedim>, querying the first
-   * whenever the passed @p manifold can be upgraded to a
-   * Boundary<dim,spacedim>.
-   */
-  template <int dim, int spacedim, class TriaIterator>
-  void get_intermediate_points_on_object(const Manifold<dim, spacedim> &manifold,
-                                         const QGaussLobatto<1>        &line_support_points,
-                                         const TriaIterator &iter,
-                                         std::vector<Point<spacedim> > &points)
-  {
-    const unsigned int structdim = TriaIterator::AccessorType::structure_dimension;
-
-    // Try backward compatibility option.
-    if (const Boundary<dim,spacedim> *boundary
-        = dynamic_cast<const Boundary<dim,spacedim> *>(&manifold))
-      // This is actually a boundary. Call old methods.
-      {
-        switch (structdim)
-          {
-          case 1:
-          {
-            const typename Triangulation<dim,spacedim>::line_iterator line = iter;
-            boundary->get_intermediate_points_on_line(line, points);
-            return;
-          }
-          case 2:
-          {
-            const typename Triangulation<dim,spacedim>::quad_iterator quad = iter;
-            boundary->get_intermediate_points_on_quad(quad, points);
-            return;
-          }
-          default:
-            Assert(false, ExcInternalError());
-            return;
-          }
-      }
-    else
-      {
-        std::vector<Point<spacedim> > sp(GeometryInfo<structdim>::vertices_per_cell);
-        for (unsigned int i=0; i<sp.size(); ++i)
-          sp[i] = iter->vertex(i);
-        get_intermediate_points(manifold, line_support_points, sp, points);
-      }
-  }
-
-
-  /**
-   * Take a <tt>support_point_weights_on_hex(quad)</tt> and apply it to the vector
-   * @p a to compute the inner support points as a linear combination of the
-   * exterior points.
-   *
-   * The vector @p a initially contains the locations of the @p n_outer
-   * points, the @p n_inner computed inner points are appended.
-   *
-   * See equation (7) of the `mapping' report.
-   */
-  template <int spacedim>
-  void add_weighted_interior_points(const Table<2,double>   &lvs,
-                                    std::vector<Point<spacedim> > &a)
-  {
-    const unsigned int n_inner_apply=lvs.n_rows();
-    const unsigned int n_outer_apply=lvs.n_cols();
-    Assert(a.size()==n_outer_apply,
-           ExcDimensionMismatch(a.size(), n_outer_apply));
-
-    // compute each inner point as linear combination of the outer points. the
-    // weights are given by the lvs entries, the outer points are the first
-    // (existing) elements of a
-    for (unsigned int unit_point=0; unit_point<n_inner_apply; ++unit_point)
-      {
-        Assert(lvs.n_cols()==n_outer_apply, ExcInternalError());
-        Point<spacedim> p;
-        for (unsigned int k=0; k<n_outer_apply; ++k)
-          p+=lvs[unit_point][k]*a[k];
-
-        a.push_back(p);
-      }
-  }
-}
-
 //--------------------------- Explicit instantiations -----------------------
 #include "mapping_manifold.inst"
 
diff --git a/tests/fe/mapping_manifold_07.cc b/tests/fe/mapping_manifold_07.cc
index c8bd6af93d..5730688398 100644
--- a/tests/fe/mapping_manifold_07.cc
+++ b/tests/fe/mapping_manifold_07.cc
@@ -13,9 +13,7 @@
 //
 // ---------------------------------------------------------------------
 
-// Check MappingManifold on faces. Create a sphere and make sure all
-// points returned by FEValues.get_quadrature_points are actually on
-// the circle.
+// Try to compute the area of a circle/sphere using JxW values.
 
 #include "../tests.h"
 
@@ -48,33 +46,33 @@ void test()
   triangulation.set_all_manifold_ids_on_boundary(0);
   triangulation.set_manifold (0, manifold);
 
-  MappingManifold<dim,spacedim> map_manifold;
-  FE_Q<dim,spacedim> fe(1);
-  const QGauss<dim-1> quad(5);
-
-  FEFaceValues<dim,spacedim> fe_v(map_manifold, fe, quad,
-                                  update_quadrature_points);
-
-  out << "set size ratio -1" << std::endl
-      << "set terminal aqua " << dim+spacedim << std::endl
-      << (spacedim == 3 ? "s" : "") << "plot '-' w l " << std::endl;
-
-  for (typename Triangulation<dim,spacedim>::active_cell_iterator cell =
-         triangulation.begin_active(); cell != triangulation.end(); ++cell)
-    for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-      if (cell->face(f)->at_boundary())
-        {
-          fe_v.reinit(cell, f);
-          const std::vector<Point<spacedim> > &qps = fe_v.get_quadrature_points();
-          for (unsigned int i=0; i<qps.size(); ++i)
+  for (unsigned int cycle = 0; cycle<4; ++cycle)
+    {
+      MappingManifold<dim,spacedim> map_manifold;
+      FE_Q<dim,spacedim> fe(1);
+      const QGauss<dim-1> quad(3);
+
+      FEFaceValues<dim,spacedim> fe_v(map_manifold,
+                                      fe, quad,
+                                      update_JxW_values);
+      double area = 0;
+
+      for (typename Triangulation<dim,spacedim>::active_cell_iterator cell =
+             triangulation.begin_active(); cell != triangulation.end(); ++cell)
+        for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+          if (cell->face(f)->at_boundary())
             {
-              out << qps[i] << std::endl;
-              if (std::abs(qps[i].distance(center) -radius) > 1e-10)
-                out << "# Error! This should be on the sphere, but it's not!" << std::endl;
+              fe_v.reinit(cell, f);
+              for (unsigned int i=0; i<quad.size(); ++i)
+                area += fe_v.JxW(i);
             }
-          out << std::endl;
-        }
-  out << "e" << std::endl << std::endl;
+      deallog << "Cycle	      : " << cycle << std::endl;
+      deallog << "Surface Area  : " << area << std::endl;
+      deallog << "Error         : " << (area-(dim-1)*2*numbers::PI) << std::endl;
+
+      triangulation.refine_global(1);
+    }
+
 }
 
 
diff --git a/tests/fe/mapping_manifold_07.output b/tests/fe/mapping_manifold_07.output
index 97d058a327..1465b404ad 100644
--- a/tests/fe/mapping_manifold_07.output
+++ b/tests/fe/mapping_manifold_07.output
@@ -1,193 +1,27 @@
 
 # dim=2, spacedim=2
-set size ratio -1
-set terminal aqua 4
-plot '-' w l 
-0.846869 1.74276
-1.08958 1.58810
-1.50000 1.50000
-1.91042 1.58810
-2.15313 1.74276
-
-0.742755 1.84687
-0.588103 2.08958
-0.500000 2.50000
-0.588103 2.91042
-0.742755 3.15313
-
-2.25724 1.84687
-2.41190 2.08958
-2.50000 2.50000
-2.41190 2.91042
-2.25724 3.15313
-
-0.846869 3.25724
-1.08958 3.41190
-1.50000 3.50000
-1.91042 3.41190
-2.15313 3.25724
-
-e
-
+DEAL::Cycle	      : 0
+DEAL::Surface Area  : 6.28319
+DEAL::Error         : -1.77636e-15
+DEAL::Cycle	      : 1
+DEAL::Surface Area  : 6.28319
+DEAL::Error         : 0.00000
+DEAL::Cycle	      : 2
+DEAL::Surface Area  : 6.28319
+DEAL::Error         : 1.77636e-15
+DEAL::Cycle	      : 3
+DEAL::Surface Area  : 6.28319
+DEAL::Error         : -8.88178e-15
 # dim=3, spacedim=3
-set size ratio -1
-set terminal aqua 6
-splot '-' w l 
-0.942530 1.94253 -0.615187
-1.12940 1.87632 -0.688247
-1.50000 1.82851 -0.741013
-1.87060 1.87632 -0.688247
-2.05747 1.94253 -0.615187
-0.876324 2.12940 -0.688247
-1.07160 2.07160 -0.795583
-1.50000 2.02589 -0.880468
-1.92840 2.07160 -0.795583
-2.12368 2.12940 -0.688247
-0.828509 2.50000 -0.741013
-1.02589 2.50000 -0.880468
-1.50000 2.50000 -1.00000
-1.97411 2.50000 -0.880468
-2.17149 2.50000 -0.741013
-0.876324 2.87060 -0.688247
-1.07160 2.92840 -0.795583
-1.50000 2.97411 -0.880468
-1.92840 2.92840 -0.795583
-2.12368 2.87060 -0.688247
-0.942530 3.05747 -0.615187
-1.12940 3.12368 -0.688247
-1.50000 3.17149 -0.741013
-1.87060 3.12368 -0.688247
-2.05747 3.05747 -0.615187
-
-2.11519 1.94253 -0.557470
-2.18825 1.87632 -0.370600
-2.24101 1.82851 -8.90587e-18
-2.18825 1.87632 0.370600
-2.11519 1.94253 0.557470
-2.18825 2.12940 -0.623676
-2.29558 2.07160 -0.428397
-2.38047 2.02589 -2.11638e-17
-2.29558 2.07160 0.428397
-2.18825 2.12940 0.623676
-2.24101 2.50000 -0.671491
-2.38047 2.50000 -0.474105
-2.50000 2.50000 0.00000
-2.38047 2.50000 0.474105
-2.24101 2.50000 0.671491
-2.18825 2.87060 -0.623676
-2.29558 2.92840 -0.428397
-2.38047 2.97411 0.00000
-2.29558 2.92840 0.428397
-2.18825 2.87060 0.623676
-2.11519 3.05747 -0.557470
-2.18825 3.12368 -0.370600
-2.24101 3.17149 0.00000
-2.18825 3.12368 0.370600
-2.11519 3.05747 0.557470
-
-0.942530 1.94253 0.615187
-0.876324 2.12940 0.688247
-0.828509 2.50000 0.741013
-0.876324 2.87060 0.688247
-0.942530 3.05747 0.615187
-1.12940 1.87632 0.688247
-1.07160 2.07160 0.795583
-1.02589 2.50000 0.880468
-1.07160 2.92840 0.795583
-1.12940 3.12368 0.688247
-1.50000 1.82851 0.741013
-1.50000 2.02589 0.880468
-1.50000 2.50000 1.00000
-1.50000 2.97411 0.880468
-1.50000 3.17149 0.741013
-1.87060 1.87632 0.688247
-1.92840 2.07160 0.795583
-1.97411 2.50000 0.880468
-1.92840 2.92840 0.795583
-1.87060 3.12368 0.688247
-2.05747 1.94253 0.615187
-2.12368 2.12940 0.688247
-2.17149 2.50000 0.741013
-2.12368 2.87060 0.688247
-2.05747 3.05747 0.615187
-
-0.884813 1.94253 -0.557470
-0.811753 2.12940 -0.623676
-0.758987 2.50000 -0.671491
-0.811753 2.87060 -0.623676
-0.884813 3.05747 -0.557470
-0.811753 1.87632 -0.370600
-0.704417 2.07160 -0.428397
-0.619532 2.50000 -0.474105
-0.704417 2.92840 -0.428397
-0.811753 3.12368 -0.370600
-0.758987 1.82851 -8.90587e-18
-0.619532 2.02589 -2.11638e-17
-0.500000 2.50000 0.00000
-0.619532 2.97411 0.00000
-0.758987 3.17149 0.00000
-0.811753 1.87632 0.370600
-0.704417 2.07160 0.428397
-0.619532 2.50000 0.474105
-0.704417 2.92840 0.428397
-0.811753 3.12368 0.370600
-0.884813 1.94253 0.557470
-0.811753 2.12940 0.623676
-0.758987 2.50000 0.671491
-0.811753 2.87060 0.623676
-0.884813 3.05747 0.557470
-
-0.942530 1.88481 -0.557470
-0.876324 1.81175 -0.370600
-0.828509 1.75899 -8.90587e-18
-0.876324 1.81175 0.370600
-0.942530 1.88481 0.557470
-1.12940 1.81175 -0.623676
-1.07160 1.70442 -0.428397
-1.02589 1.61953 -2.11638e-17
-1.07160 1.70442 0.428397
-1.12940 1.81175 0.623676
-1.50000 1.75899 -0.671491
-1.50000 1.61953 -0.474105
-1.50000 1.50000 0.00000
-1.50000 1.61953 0.474105
-1.50000 1.75899 0.671491
-1.87060 1.81175 -0.623676
-1.92840 1.70442 -0.428397
-1.97411 1.61953 0.00000
-1.92840 1.70442 0.428397
-1.87060 1.81175 0.623676
-2.05747 1.88481 -0.557470
-2.12368 1.81175 -0.370600
-2.17149 1.75899 0.00000
-2.12368 1.81175 0.370600
-2.05747 1.88481 0.557470
-
-0.942530 3.11519 -0.557470
-1.12940 3.18825 -0.623676
-1.50000 3.24101 -0.671491
-1.87060 3.18825 -0.623676
-2.05747 3.11519 -0.557470
-0.876324 3.18825 -0.370600
-1.07160 3.29558 -0.428397
-1.50000 3.38047 -0.474105
-1.92840 3.29558 -0.428397
-2.12368 3.18825 -0.370600
-0.828509 3.24101 -8.90587e-18
-1.02589 3.38047 -2.11638e-17
-1.50000 3.50000 0.00000
-1.97411 3.38047 0.00000
-2.17149 3.24101 0.00000
-0.876324 3.18825 0.370600
-1.07160 3.29558 0.428397
-1.50000 3.38047 0.474105
-1.92840 3.29558 0.428397
-2.12368 3.18825 0.370600
-0.942530 3.11519 0.557470
-1.12940 3.18825 0.623676
-1.50000 3.24101 0.671491
-1.87060 3.18825 0.623676
-2.05747 3.11519 0.557470
-
-e
-
+DEAL::Cycle	      : 0
+DEAL::Surface Area  : 9.63008
+DEAL::Error         : -2.93629
+DEAL::Cycle	      : 1
+DEAL::Surface Area  : 12.2948
+DEAL::Error         : -0.271577
+DEAL::Cycle	      : 2
+DEAL::Surface Area  : 12.4659
+DEAL::Error         : -0.100470
+DEAL::Cycle	      : 3
+DEAL::Surface Area  : 12.5217
+DEAL::Error         : -0.0446494
-- 
2.39.5