mapping_q_generic.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2000 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------


#include <deal.II/base/array_view.h>
#include <deal.II/base/derivative_form.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/qprojector.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/table.h>
#include <deal.II/base/tensor_product_polynomials.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/std_cxx14/memory.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/tensor_product_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_boundary.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q_generic.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/matrix_free/tensor_product_kernels.h>
#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/matrix_free/evaluation_kernels.h>
#include <deal.II/matrix_free/evaluation_selector.h>

#include <cmath>
#include <algorithm>
#include <numeric>
#include <array>
#include <memory>


DEAL_II_NAMESPACE_OPEN


namespace internal
{
  namespace MappingQGenericImplementation
  {
    namespace
    {
      template <int dim>
      std::vector<unsigned int>
      get_dpo_vector (const unsigned int degree)
      {
        std::vector<unsigned int> dpo(dim+1, 1U);
        for (unsigned int i=1; i<dpo.size(); ++i)
          dpo[i]=dpo[i-1]*(degree-1);
        return dpo;
      }
    }
  }

  namespace MappingQ1
  {
    namespace
    {
      // These are left as templates on the spatial dimension (even though dim
      // == spacedim must be true for them to make sense) because templates are
      // expanded before the compiler eliminates code due to the 'if (dim ==
      // spacedim)' statement (see the body of the general
      // transform_real_to_unit_cell).
      template <int spacedim>
      Point<1>
      transform_real_to_unit_cell
      (const std::array<Point<spacedim>, GeometryInfo<1>::vertices_per_cell> &vertices,
       const Point<spacedim> &p)
      {
        Assert(spacedim == 1, ExcInternalError());
        return Point<1>((p[0] - vertices[0](0))/(vertices[1](0) - vertices[0](0)));
      }



      template <int spacedim>
      Point<2>
      transform_real_to_unit_cell
      (const std::array<Point<spacedim>, GeometryInfo<2>::vertices_per_cell> &vertices,
       const Point<spacedim> &p)
      {
        Assert(spacedim == 2, ExcInternalError());

        // For accuracy reasons, we do all arithmetics in extended precision
        // (long double). This has a noticeable effect on the hit rate for
        // borderline cases and thus makes the algorithm more robust.
        const long double x = p(0);
        const long double y = p(1);

        const long double x0 = vertices[0](0);
        const long double x1 = vertices[1](0);
        const long double x2 = vertices[2](0);
        const long double x3 = vertices[3](0);

        const long double y0 = vertices[0](1);
        const long double y1 = vertices[1](1);
        const long double y2 = vertices[2](1);
        const long double y3 = vertices[3](1);

        const long double a = (x1 - x3)*(y0 - y2) - (x0 - x2)*(y1 - y3);
        const long double b = -(x0 - x1 - x2 + x3)*y + (x - 2*x1 + x3)*y0 - (x - 2*x0 + x2)*y1
                              - (x - x1)*y2 + (x - x0)*y3;
        const long double c = (x0 - x1)*y - (x - x1)*y0 + (x - x0)*y1;

        const long double discriminant = b*b - 4*a*c;
        // exit if the point is not in the cell (this is the only case where the
        // discriminant is negative)
        AssertThrow (discriminant > 0.0,
                     (typename Mapping<spacedim,spacedim>::ExcTransformationFailed()));

        long double eta1;
        long double eta2;
        const long double sqrt_discriminant = std::sqrt(discriminant);
        // special case #1: if a is near-zero to make the discriminant exactly
        // equal b, then use the linear formula
        if (b != 0.0 && std::abs(b) == sqrt_discriminant)
          {
            eta1 = -c/b;
            eta2 = -c/b;
          }
        // special case #2: a is zero for parallelograms and very small for
        // near-parallelograms:
        else if (std::abs(a) < 1e-8*std::abs(b))
          {
            // if both a and c are very small then the root should be near
            // zero: this first case will capture that
            eta1 = 2*c / (-b - sqrt_discriminant);
            eta2 = 2*c / (-b + sqrt_discriminant);
          }
        // finally, use the plain version:
        else
          {
            eta1 = (-b - sqrt_discriminant) / (2*a);
            eta2 = (-b + sqrt_discriminant) / (2*a);
          }
        // pick the one closer to the center of the cell.
        const long double eta = (std::abs(eta1 - 0.5) < std::abs(eta2 - 0.5)) ? eta1 : eta2;

        /*
         * There are two ways to compute xi from eta, but either one may have a
         * zero denominator.
         */
        const long double subexpr0 = -eta*x2 + x0*(eta - 1);
        const long double xi_denominator0 = eta*x3 - x1*(eta - 1) + subexpr0;
        const long double max_x = std::max(std::max(std::abs(x0), std::abs(x1)),
                                           std::max(std::abs(x2), std::abs(x3)));

        if (std::abs(xi_denominator0) > 1e-10*max_x)
          {
            const double xi = (x + subexpr0)/xi_denominator0;
            return Point<2>(xi, eta);
          }
        else
          {
            const long double max_y = std::max(std::max(std::abs(y0), std::abs(y1)),
                                               std::max(std::abs(y2), std::abs(y3)));
            const long double subexpr1 = -eta*y2 + y0*(eta - 1);
            const long double xi_denominator1 = eta*y3 - y1*(eta - 1) + subexpr1;
            if (std::abs(xi_denominator1) > 1e-10*max_y)
              {
                const double xi = (subexpr1 + y)/xi_denominator1;
                return Point<2>(xi, eta);
              }
            else // give up and try Newton iteration
              {
                AssertThrow (false,
                             (typename Mapping<spacedim,spacedim>::ExcTransformationFailed()));
              }
          }
        // bogus return to placate compiler. It should not be possible to get
        // here.
        Assert(false, ExcInternalError());
        return Point<2>(std::numeric_limits<double>::quiet_NaN(),
                        std::numeric_limits<double>::quiet_NaN());
      }



      template <int spacedim>
      Point<3>
      transform_real_to_unit_cell
      (const std::array<Point<spacedim>, GeometryInfo<3>::vertices_per_cell> &/*vertices*/,
       const Point<spacedim> &/*p*/)
      {
        // It should not be possible to get here
        Assert(false, ExcInternalError());
        return Point<3>();
      }



      template <int dim, int spacedim>
      void compute_shape_function_values_general (const unsigned int            n_shape_functions,
                                                  const std::vector<Point<dim> > &unit_points,
                                                  typename ::MappingQGeneric<dim,spacedim>::InternalData &data)
      {
        const unsigned int n_points=unit_points.size();

        // Construct the tensor product polynomials used as shape functions for the
        // Qp mapping of cells at the boundary.
        const TensorProductPolynomials<dim>
        tensor_pols (Polynomials::generate_complete_Lagrange_basis(data.line_support_points.get_points()));
        Assert (n_shape_functions==tensor_pols.n(),
                ExcInternalError());

        // then also construct the mapping from lexicographic to the Qp shape function numbering
        const std::vector<unsigned int>
        renumber (FETools::
                  lexicographic_to_hierarchic_numbering
                  (FiniteElementData<dim> (internal::MappingQGenericImplementation::get_dpo_vector<dim>
                                           (data.polynomial_degree), 1, data.polynomial_degree)));

        std::vector<double> values;
        std::vector<Tensor<1,dim> > grads;
        if (data.shape_values.size()!=0)
          {
            Assert(data.shape_values.size()==n_shape_functions*n_points,
                   ExcInternalError());
            values.resize(n_shape_functions);
          }
        if (data.shape_derivatives.size()!=0)
          {
            Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
                   ExcInternalError());
            grads.resize(n_shape_functions);
          }

        std::vector<Tensor<2,dim> > grad2;
        if (data.shape_second_derivatives.size()!=0)
          {
            Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points,
                   ExcInternalError());
            grad2.resize(n_shape_functions);
          }

        std::vector<Tensor<3,dim> > grad3;
        if (data.shape_third_derivatives.size()!=0)
          {
            Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points,
                   ExcInternalError());
            grad3.resize(n_shape_functions);
          }

        std::vector<Tensor<4,dim> > grad4;
        if (data.shape_fourth_derivatives.size()!=0)
          {
            Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points,
                   ExcInternalError());
            grad4.resize(n_shape_functions);
          }


        if (data.shape_values.size()!=0 ||
            data.shape_derivatives.size()!=0 ||
            data.shape_second_derivatives.size()!=0 ||
            data.shape_third_derivatives.size()!=0 ||
            data.shape_fourth_derivatives.size()!=0 )
          for (unsigned int point=0; point<n_points; ++point)
            {
              tensor_pols.compute(unit_points[point], values, grads, grad2, grad3, grad4);

              if (data.shape_values.size()!=0)
                for (unsigned int i=0; i<n_shape_functions; ++i)
                  data.shape(point,renumber[i]) = values[i];

              if (data.shape_derivatives.size()!=0)
                for (unsigned int i=0; i<n_shape_functions; ++i)
                  data.derivative(point,renumber[i]) = grads[i];

              if (data.shape_second_derivatives.size()!=0)
                for (unsigned int i=0; i<n_shape_functions; ++i)
                  data.second_derivative(point,renumber[i]) = grad2[i];

              if (data.shape_third_derivatives.size()!=0)
                for (unsigned int i=0; i<n_shape_functions; ++i)
                  data.third_derivative(point,renumber[i]) = grad3[i];

              if (data.shape_fourth_derivatives.size()!=0)
                for (unsigned int i=0; i<n_shape_functions; ++i)
                  data.fourth_derivative(point,renumber[i]) = grad4[i];
            }
      }


      void
      compute_shape_function_values_hardcode (const unsigned int            n_shape_functions,
                                              const std::vector<Point<1> > &unit_points,
                                              ::MappingQGeneric<1,1>::InternalData &data)
      {
        (void)n_shape_functions;
        const unsigned int n_points=unit_points.size();
        for (unsigned int k = 0 ; k < n_points ; ++k)
          {
            double x = unit_points[k](0);

            if (data.shape_values.size()!=0)
              {
                Assert(data.shape_values.size()==n_shape_functions*n_points,
                       ExcInternalError());
                data.shape(k,0) = 1.-x;
                data.shape(k,1) = x;
              }
            if (data.shape_derivatives.size()!=0)
              {
                Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
                       ExcInternalError());
                data.derivative(k,0)[0] = -1.;
                data.derivative(k,1)[0] = 1.;
              }
            if (data.shape_second_derivatives.size()!=0)
              {
                Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points,
                       ExcInternalError());
                data.second_derivative(k,0)[0][0] = 0;
                data.second_derivative(k,1)[0][0] = 0;
              }
            if (data.shape_third_derivatives.size()!=0)
              {
                Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points,
                       ExcInternalError());

                Tensor<3,1> zero;
                data.third_derivative(k,0) = zero;
                data.third_derivative(k,1) = zero;
              }
            if (data.shape_fourth_derivatives.size()!=0)
              {
                Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points,
                       ExcInternalError());

                Tensor<4,1> zero;
                data.fourth_derivative(k,0) = zero;
                data.fourth_derivative(k,1) = zero;
              }
          }
      }


      void
      compute_shape_function_values_hardcode (const unsigned int            n_shape_functions,
                                              const std::vector<Point<2> > &unit_points,
                                              ::MappingQGeneric<2,2>::InternalData &data)
      {

        (void)n_shape_functions;
        const unsigned int n_points=unit_points.size();
        for (unsigned int k = 0 ; k < n_points ; ++k)
          {
            double x = unit_points[k](0);
            double y = unit_points[k](1);

            if (data.shape_values.size()!=0)
              {
                Assert(data.shape_values.size()==n_shape_functions*n_points,
                       ExcInternalError());
                data.shape(k,0) = (1.-x)*(1.-y);
                data.shape(k,1) = x*(1.-y);
                data.shape(k,2) = (1.-x)*y;
                data.shape(k,3) = x*y;
              }
            if (data.shape_derivatives.size()!=0)
              {
                Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
                       ExcInternalError());
                data.derivative(k,0)[0] = (y-1.);
                data.derivative(k,1)[0] = (1.-y);
                data.derivative(k,2)[0] = -y;
                data.derivative(k,3)[0] = y;
                data.derivative(k,0)[1] = (x-1.);
                data.derivative(k,1)[1] = -x;
                data.derivative(k,2)[1] = (1.-x);
                data.derivative(k,3)[1] = x;
              }
            if (data.shape_second_derivatives.size()!=0)
              {
                Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points,
                       ExcInternalError());
                data.second_derivative(k,0)[0][0] = 0;
                data.second_derivative(k,1)[0][0] = 0;
                data.second_derivative(k,2)[0][0] = 0;
                data.second_derivative(k,3)[0][0] = 0;
                data.second_derivative(k,0)[0][1] = 1.;
                data.second_derivative(k,1)[0][1] = -1.;
                data.second_derivative(k,2)[0][1] = -1.;
                data.second_derivative(k,3)[0][1] = 1.;
                data.second_derivative(k,0)[1][0] = 1.;
                data.second_derivative(k,1)[1][0] = -1.;
                data.second_derivative(k,2)[1][0] = -1.;
                data.second_derivative(k,3)[1][0] = 1.;
                data.second_derivative(k,0)[1][1] = 0;
                data.second_derivative(k,1)[1][1] = 0;
                data.second_derivative(k,2)[1][1] = 0;
                data.second_derivative(k,3)[1][1] = 0;
              }
            if (data.shape_third_derivatives.size()!=0)
              {
                Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points,
                       ExcInternalError());

                Tensor<3,2> zero;
                for (unsigned int i=0; i<4; ++i)
                  data.third_derivative(k,i) = zero;
              }
            if (data.shape_fourth_derivatives.size()!=0)
              {
                Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points,
                       ExcInternalError());
                Tensor<4,2> zero;
                for (unsigned int i=0; i<4; ++i)
                  data.fourth_derivative(k,i) = zero;
              }
          }
      }



      void
      compute_shape_function_values_hardcode (const unsigned int            n_shape_functions,
                                              const std::vector<Point<3> > &unit_points,
                                              ::MappingQGeneric<3,3>::InternalData &data)
      {
        (void)n_shape_functions;
        const unsigned int n_points=unit_points.size();
        for (unsigned int k = 0 ; k < n_points ; ++k)
          {
            double x = unit_points[k](0);
            double y = unit_points[k](1);
            double z = unit_points[k](2);

            if (data.shape_values.size()!=0)
              {
                Assert(data.shape_values.size()==n_shape_functions*n_points,
                       ExcInternalError());
                data.shape(k,0) = (1.-x)*(1.-y)*(1.-z);
                data.shape(k,1) = x*(1.-y)*(1.-z);
                data.shape(k,2) = (1.-x)*y*(1.-z);
                data.shape(k,3) = x*y*(1.-z);
                data.shape(k,4) = (1.-x)*(1.-y)*z;
                data.shape(k,5) = x*(1.-y)*z;
                data.shape(k,6) = (1.-x)*y*z;
                data.shape(k,7) = x*y*z;
              }
            if (data.shape_derivatives.size()!=0)
              {
                Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
                       ExcInternalError());
                data.derivative(k,0)[0] = (y-1.)*(1.-z);
                data.derivative(k,1)[0] = (1.-y)*(1.-z);
                data.derivative(k,2)[0] = -y*(1.-z);
                data.derivative(k,3)[0] = y*(1.-z);
                data.derivative(k,4)[0] = (y-1.)*z;
                data.derivative(k,5)[0] = (1.-y)*z;
                data.derivative(k,6)[0] = -y*z;
                data.derivative(k,7)[0] = y*z;
                data.derivative(k,0)[1] = (x-1.)*(1.-z);
                data.derivative(k,1)[1] = -x*(1.-z);
                data.derivative(k,2)[1] = (1.-x)*(1.-z);
                data.derivative(k,3)[1] = x*(1.-z);
                data.derivative(k,4)[1] = (x-1.)*z;
                data.derivative(k,5)[1] = -x*z;
                data.derivative(k,6)[1] = (1.-x)*z;
                data.derivative(k,7)[1] = x*z;
                data.derivative(k,0)[2] = (x-1)*(1.-y);
                data.derivative(k,1)[2] = x*(y-1.);
                data.derivative(k,2)[2] = (x-1.)*y;
                data.derivative(k,3)[2] = -x*y;
                data.derivative(k,4)[2] = (1.-x)*(1.-y);
                data.derivative(k,5)[2] = x*(1.-y);
                data.derivative(k,6)[2] = (1.-x)*y;
                data.derivative(k,7)[2] = x*y;
              }
            if (data.shape_second_derivatives.size()!=0)
              {
                Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points,
                       ExcInternalError());
                data.second_derivative(k,0)[0][0] = 0;
                data.second_derivative(k,1)[0][0] = 0;
                data.second_derivative(k,2)[0][0] = 0;
                data.second_derivative(k,3)[0][0] = 0;
                data.second_derivative(k,4)[0][0] = 0;
                data.second_derivative(k,5)[0][0] = 0;
                data.second_derivative(k,6)[0][0] = 0;
                data.second_derivative(k,7)[0][0] = 0;
                data.second_derivative(k,0)[1][1] = 0;
                data.second_derivative(k,1)[1][1] = 0;
                data.second_derivative(k,2)[1][1] = 0;
                data.second_derivative(k,3)[1][1] = 0;
                data.second_derivative(k,4)[1][1] = 0;
                data.second_derivative(k,5)[1][1] = 0;
                data.second_derivative(k,6)[1][1] = 0;
                data.second_derivative(k,7)[1][1] = 0;
                data.second_derivative(k,0)[2][2] = 0;
                data.second_derivative(k,1)[2][2] = 0;
                data.second_derivative(k,2)[2][2] = 0;
                data.second_derivative(k,3)[2][2] = 0;
                data.second_derivative(k,4)[2][2] = 0;
                data.second_derivative(k,5)[2][2] = 0;
                data.second_derivative(k,6)[2][2] = 0;
                data.second_derivative(k,7)[2][2] = 0;

                data.second_derivative(k,0)[0][1] = (1.-z);
                data.second_derivative(k,1)[0][1] = -(1.-z);
                data.second_derivative(k,2)[0][1] = -(1.-z);
                data.second_derivative(k,3)[0][1] = (1.-z);
                data.second_derivative(k,4)[0][1] = z;
                data.second_derivative(k,5)[0][1] = -z;
                data.second_derivative(k,6)[0][1] = -z;
                data.second_derivative(k,7)[0][1] = z;
                data.second_derivative(k,0)[1][0] = (1.-z);
                data.second_derivative(k,1)[1][0] = -(1.-z);
                data.second_derivative(k,2)[1][0] = -(1.-z);
                data.second_derivative(k,3)[1][0] = (1.-z);
                data.second_derivative(k,4)[1][0] = z;
                data.second_derivative(k,5)[1][0] = -z;
                data.second_derivative(k,6)[1][0] = -z;
                data.second_derivative(k,7)[1][0] = z;

                data.second_derivative(k,0)[0][2] = (1.-y);
                data.second_derivative(k,1)[0][2] = -(1.-y);
                data.second_derivative(k,2)[0][2] = y;
                data.second_derivative(k,3)[0][2] = -y;
                data.second_derivative(k,4)[0][2] = -(1.-y);
                data.second_derivative(k,5)[0][2] = (1.-y);
                data.second_derivative(k,6)[0][2] = -y;
                data.second_derivative(k,7)[0][2] = y;
                data.second_derivative(k,0)[2][0] = (1.-y);
                data.second_derivative(k,1)[2][0] = -(1.-y);
                data.second_derivative(k,2)[2][0] = y;
                data.second_derivative(k,3)[2][0] = -y;
                data.second_derivative(k,4)[2][0] = -(1.-y);
                data.second_derivative(k,5)[2][0] = (1.-y);
                data.second_derivative(k,6)[2][0] = -y;
                data.second_derivative(k,7)[2][0] = y;

                data.second_derivative(k,0)[1][2] = (1.-x);
                data.second_derivative(k,1)[1][2] = x;
                data.second_derivative(k,2)[1][2] = -(1.-x);
                data.second_derivative(k,3)[1][2] = -x;
                data.second_derivative(k,4)[1][2] = -(1.-x);
                data.second_derivative(k,5)[1][2] = -x;
                data.second_derivative(k,6)[1][2] = (1.-x);
                data.second_derivative(k,7)[1][2] = x;
                data.second_derivative(k,0)[2][1] = (1.-x);
                data.second_derivative(k,1)[2][1] = x;
                data.second_derivative(k,2)[2][1] = -(1.-x);
                data.second_derivative(k,3)[2][1] = -x;
                data.second_derivative(k,4)[2][1] = -(1.-x);
                data.second_derivative(k,5)[2][1] = -x;
                data.second_derivative(k,6)[2][1] = (1.-x);
                data.second_derivative(k,7)[2][1] = x;
              }
            if (data.shape_third_derivatives.size()!=0)
              {
                Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points,
                       ExcInternalError());

                for (unsigned int i=0; i<3; ++i)
                  for (unsigned int j=0; j<3; ++j)
                    for (unsigned int l=0; l<3; ++l)
                      if ((i==j)||(j==l)||(l==i))
                        {
                          for (unsigned int m=0; m<8; ++m)
                            data.third_derivative(k,m)[i][j][l] = 0;
                        }
                      else
                        {
                          data.third_derivative(k,0)[i][j][l] = -1.;
                          data.third_derivative(k,1)[i][j][l] = 1.;
                          data.third_derivative(k,2)[i][j][l] = 1.;
                          data.third_derivative(k,3)[i][j][l] = -1.;
                          data.third_derivative(k,4)[i][j][l] = 1.;
                          data.third_derivative(k,5)[i][j][l] = -1.;
                          data.third_derivative(k,6)[i][j][l] = -1.;
                          data.third_derivative(k,7)[i][j][l] = 1.;
                        }

              }
            if (data.shape_fourth_derivatives.size()!=0)
              {
                Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points,
                       ExcInternalError());
                Tensor<4,3> zero;
                for (unsigned int i=0; i<8; ++i)
                  data.fourth_derivative(k,i) = zero;
              }
          }
      }
    }
  }
}



template <int dim, int spacedim>
MappingQGeneric<dim,spacedim>::InternalData::InternalData (const unsigned int polynomial_degree)
  :
  polynomial_degree (polynomial_degree),
  n_shape_functions (Utilities::fixed_power<dim>(polynomial_degree+1)),
  line_support_points(QGaussLobatto<1>(polynomial_degree+1)),
  tensor_product_quadrature(false)
{}



template <int dim, int spacedim>
std::size_t
MappingQGeneric<dim,spacedim>::InternalData::memory_consumption () const
{
  return (Mapping<dim,spacedim>::InternalDataBase::memory_consumption() +
          MemoryConsumption::memory_consumption (shape_values) +
          MemoryConsumption::memory_consumption (shape_derivatives) +
          MemoryConsumption::memory_consumption (covariant) +
          MemoryConsumption::memory_consumption (contravariant) +
          MemoryConsumption::memory_consumption (unit_tangentials) +
          MemoryConsumption::memory_consumption (aux) +
          MemoryConsumption::memory_consumption (mapping_support_points) +
          MemoryConsumption::memory_consumption (cell_of_current_support_points) +
          MemoryConsumption::memory_consumption (volume_elements) +
          MemoryConsumption::memory_consumption (polynomial_degree) +
          MemoryConsumption::memory_consumption (n_shape_functions));
}


template <int dim, int spacedim>
void
MappingQGeneric<dim,spacedim>::InternalData::
initialize (const UpdateFlags      update_flags,
            const Quadrature<dim> &q,
            const unsigned int     n_original_q_points)
{
  // store the flags in the internal data object so we can access them
  // in fill_fe_*_values()
  this->update_each = update_flags;

  const unsigned int n_q_points = q.size();

  // see if we need the (transformation) shape function values
  // and/or gradients and resize the necessary arrays
  if (this->update_each & update_quadrature_points)
    shape_values.resize(n_shape_functions * n_q_points);

  if (this->update_each & (update_covariant_transformation
                           | update_contravariant_transformation
                           | update_JxW_values
                           | update_boundary_forms
                           | update_normal_vectors
                           | update_jacobians
                           | update_jacobian_grads
                           | update_inverse_jacobians
                           | update_jacobian_pushed_forward_grads
                           | update_jacobian_2nd_derivatives
                           | update_jacobian_pushed_forward_2nd_derivatives
                           | update_jacobian_3rd_derivatives
                           | update_jacobian_pushed_forward_3rd_derivatives))
    shape_derivatives.resize(n_shape_functions * n_q_points);

  if (this->update_each & update_covariant_transformation)
    covariant.resize(n_original_q_points);

  if (this->update_each & update_contravariant_transformation)
    contravariant.resize(n_original_q_points);

  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);

  const std::vector<Point<dim> > &ref_q_points = q.get_points();
  // now also fill the various fields with their correct values
  compute_shape_function_values (ref_q_points);

  tensor_product_quadrature = q.is_tensor_product();

  if (dim>1)
    {
      // find out if the one-dimensional formula is the same
      // in all directions
      if (tensor_product_quadrature)
        {
          const std::array<Quadrature<1>, dim> quad_array = q.get_tensor_basis();
          for (unsigned int i=1; i<dim && tensor_product_quadrature; ++i)
            {
              if (quad_array[i-1].size() != quad_array[i].size())
                {
                  tensor_product_quadrature = false;
                  break;
                }
              else
                {
                  const std::vector<Point<1>> &points_1 = quad_array[i-1].get_points();
                  const std::vector<Point<1>> &points_2 = quad_array[i].get_points();
                  const std::vector<double> &weights_1 = quad_array[i-1].get_weights();
                  const std::vector<double> &weights_2 = quad_array[i].get_weights();
                  for (unsigned int j=0; j<quad_array[i].size(); ++j)
                    {
                      if (std::abs(points_1[j][0]-points_2[j][0])>1.e-10
                          || std::abs(weights_1[j]-weights_2[j])>1.e-10)
                        {
                          tensor_product_quadrature = false;
                          break;
                        }

                    }
                }
            }

          if (tensor_product_quadrature)
            {
              const FE_Q<dim> fe(polynomial_degree);
              shape_info.reinit(q.get_tensor_basis()[0], fe);

              const unsigned int n_shape_values = fe.n_dofs_per_cell();
              const unsigned int max_size = std::max(n_q_points,n_shape_values);
              const unsigned int vec_length = ::VectorizedArray<double>::n_array_elements;
              const unsigned int n_comp = 1+ (spacedim-1)/vec_length;

              scratch.resize((dim-1)*max_size);
              values_dofs.resize(n_comp*n_shape_values);
            }
        }
    }
}



template <int dim, int spacedim>
void
MappingQGeneric<dim,spacedim>::InternalData::
initialize_face (const UpdateFlags      update_flags,
                 const Quadrature<dim> &q,
                 const unsigned int     n_original_q_points)
{
  initialize (update_flags, q, n_original_q_points);

  if (dim>1 && tensor_product_quadrature)
    {
      const unsigned int facedim = dim > 1 ? dim-1 : 1;
      const FE_Q<facedim> fe(polynomial_degree);
      shape_info.reinit(q.get_tensor_basis()[0], fe);

      const unsigned int n_shape_values = fe.n_dofs_per_cell();
      const unsigned int n_q_points = q.size();
      const unsigned int max_size = std::max(n_q_points,n_shape_values);
      const unsigned int vec_length = VectorizedArray<double>::n_array_elements;
      const unsigned int n_comp = 1+ (spacedim-1)/vec_length;

      scratch.resize((dim-1)*max_size);
      values_dofs.resize(n_comp*n_shape_values);
    }

  if (dim > 1)
    {
      if (this->update_each & (update_boundary_forms |
                               update_normal_vectors |
                               update_jacobians      |
                               update_JxW_values      |
                               update_inverse_jacobians))
        {
          aux.resize (dim-1, std::vector<Tensor<1,spacedim> > (n_original_q_points));

          // Compute tangentials to the unit cell.
          for (unsigned int i=0; i<unit_tangentials.size(); ++i)
            unit_tangentials[i].resize (n_original_q_points);
          switch (dim)
            {
            case 2:
            {
              // ensure a counterclockwise orientation of tangentials
              static const int tangential_orientation[4]= {-1,1,1,-1};
              for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
                {
                  Tensor<1,dim> tang;
                  tang[1-i/2] = tangential_orientation[i];
                  std::fill (unit_tangentials[i].begin(),
                             unit_tangentials[i].end(),
                             tang);
                }

              break;
            }

            case 3:
            {
              for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
                {
                  Tensor<1,dim> tang1, tang2;

                  const unsigned int nd=
                    GeometryInfo<dim>::unit_normal_direction[i];

                  // first tangential
                  // vector in direction
                  // of the (nd+1)%3 axis
                  // and inverted in case
                  // of unit inward normal
                  tang1[(nd+1)%dim]=GeometryInfo<dim>::unit_normal_orientation[i];
                  // second tangential
                  // vector in direction
                  // of the (nd+2)%3 axis
                  tang2[(nd+2)%dim]=1.;

                  // same unit tangents
                  // for all quadrature
                  // points on this face
                  std::fill (unit_tangentials[i].begin(),
                             unit_tangentials[i].end(),
                             tang1);
                  std::fill (unit_tangentials[GeometryInfo<dim>::faces_per_cell+i].begin(),
                             unit_tangentials[GeometryInfo<dim>::faces_per_cell+i].end(),
                             tang2);
                }

              break;
            }

            default:
              Assert (false, ExcNotImplemented());
            }
        }
    }
}





template <>
void
MappingQGeneric<1,1>::InternalData::
compute_shape_function_values (const std::vector<Point<1> > &unit_points)
{
  // if the polynomial degree is one, then we can simplify code a bit
  // by using hard-coded shape functions.
  if (polynomial_degree == 1)
    internal::MappingQ1::compute_shape_function_values_hardcode (n_shape_functions,
        unit_points, *this);
  else
    {
      // otherwise ask an object that describes the polynomial space
      internal::MappingQ1::compute_shape_function_values_general<1,1>(n_shape_functions,
          unit_points,*this);
    }
}

template <>
void
MappingQGeneric<2,2>::InternalData::
compute_shape_function_values (const std::vector<Point<2> > &unit_points)
{
  // if the polynomial degree is one, then we can simplify code a bit
  // by using hard-coded shape functions.
  if (polynomial_degree == 1)
    internal::MappingQ1::compute_shape_function_values_hardcode (n_shape_functions,
        unit_points, *this);
  else
    {
      // otherwise ask an object that describes the polynomial space
      internal::MappingQ1::compute_shape_function_values_general<2,2>(n_shape_functions,
          unit_points,*this);
    }
}

template <>
void
MappingQGeneric<3,3>::InternalData::
compute_shape_function_values (const std::vector<Point<3> > &unit_points)
{
  // if the polynomial degree is one, then we can simplify code a bit
  // by using hard-coded shape functions.
  if (polynomial_degree == 1)
    internal::MappingQ1::compute_shape_function_values_hardcode (n_shape_functions,
        unit_points, *this);
  else
    {
      // otherwise ask an object that describes the polynomial space
      internal::MappingQ1::compute_shape_function_values_general<3,3>(n_shape_functions,
          unit_points,*this);
    }
}

template <int dim, int spacedim>
void
MappingQGeneric<dim,spacedim>::InternalData::
compute_shape_function_values (const std::vector<Point<dim> > &unit_points)
{
  // for non-matching combinations of dim and spacedim, just run the general
  // case
  internal::MappingQ1::compute_shape_function_values_general<dim,spacedim>(n_shape_functions,
      unit_points,*this);
}


namespace internal
{
  namespace MappingQGenericImplementation
  {
    namespace
    {
      ::Table<2,double>
      compute_support_point_weights_on_quad(const unsigned int polynomial_degree)
      {
        ::Table<2,double> loqvs;

        // we are asked to compute weights for interior support points, but
        // there are no interior points if degree==1
        if (polynomial_degree == 1)
          return loqvs;

        const unsigned int M = polynomial_degree-1;
        const unsigned int n_inner_2d = M*M;
        const unsigned int n_outer_2d = 4 + 4*M;

        // set the weights of transfinite interpolation
        loqvs.reinit(n_inner_2d, n_outer_2d);
        QGaussLobatto<2> gl(polynomial_degree+1);
        for (unsigned int i=0; i<M; ++i)
          for (unsigned int j=0; j<M; ++j)
            {
              const Point<2> p = gl.point((i+1)*(polynomial_degree+1)+(j+1));
              const unsigned int index_table = i*M+j;
              for (unsigned int v=0; v<4; ++v)
                loqvs(index_table, v) =
                  -GeometryInfo<2>::d_linear_shape_function(p, v);
              loqvs(index_table, 4+i) = 1.-p[0];
              loqvs(index_table, 4+i+M) = p[0];
              loqvs(index_table, 4+j+2*M) = 1.-p[1];
              loqvs(index_table, 4+j+3*M) = p[1];
            }

        // the sum of weights of the points at the outer rim should be one. check
        // this
        for (unsigned int unit_point=0; unit_point<n_inner_2d; ++unit_point)
          Assert(std::fabs(std::accumulate(loqvs[unit_point].begin(),
                                           loqvs[unit_point].end(),0.) - 1)<1e-13*polynomial_degree,
                 ExcInternalError());

        return loqvs;
      }



      ::Table<2,double>
      compute_support_point_weights_on_hex(const unsigned int polynomial_degree)
      {
        ::Table<2,double> lohvs;

        // we are asked to compute weights for interior support points, but
        // there are no interior points if degree==1
        if (polynomial_degree == 1)
          return lohvs;

        const unsigned int M = polynomial_degree-1;

        const unsigned int n_inner = Utilities::fixed_power<3>(M);
        const unsigned int n_outer = 8+12*M+6*M*M;

        // set the weights of transfinite interpolation
        lohvs.reinit(n_inner, n_outer);
        QGaussLobatto<3> gl(polynomial_degree+1);
        for (unsigned int i=0; i<M; ++i)
          for (unsigned int j=0; j<M; ++j)
            for (unsigned int k=0; k<M; ++k)
              {
                const Point<3> p = gl.point((i+1)*(M+2)*(M+2)+(j+1)*(M+2)+(k+1));
                const unsigned int index_table = i*M*M+j*M+k;

                // vertices
                for (unsigned int v=0; v<8; ++v)
                  lohvs(index_table, v) =
                    GeometryInfo<3>::d_linear_shape_function(p, v);

                // lines
                {
                  constexpr std::array<unsigned int,4> line_coordinates_y
                  ({{0, 1, 4, 5}});
                  const Point<2> py(p[0], p[2]);
                  for (unsigned int l=0; l<4; ++l)
                    lohvs(index_table, 8+line_coordinates_y[l]*M+j) =
                      -GeometryInfo<2>::d_linear_shape_function(py, l);
                }

                {
                  constexpr std::array<unsigned int,4> line_coordinates_x
                  ({{2, 3, 6, 7}});
                  const Point<2> px(p[1], p[2]);
                  for (unsigned int l=0; l<4; ++l)
                    lohvs(index_table, 8+line_coordinates_x[l]*M+k) =
                      -GeometryInfo<2>::d_linear_shape_function(px, l);
                }

                {
                  constexpr std::array<unsigned int,4> line_coordinates_z
                  ({{8, 9, 10, 11}});
                  const Point<2> pz(p[0], p[1]);
                  for (unsigned int l=0; l<4; ++l)
                    lohvs(index_table, 8+line_coordinates_z[l]*M+i) =
                      -GeometryInfo<2>::d_linear_shape_function(pz, l);
                }

                // quads
                lohvs(index_table, 8+12*M+0*M*M+i*M+j) = 1.-p[0];
                lohvs(index_table, 8+12*M+1*M*M+i*M+j) = p[0];
                lohvs(index_table, 8+12*M+2*M*M+k*M+i) = 1.-p[1];
                lohvs(index_table, 8+12*M+3*M*M+k*M+i) = p[1];
                lohvs(index_table, 8+12*M+4*M*M+j*M+k) = 1.-p[2];
                lohvs(index_table, 8+12*M+5*M*M+j*M+k) = p[2];
              }

        // the sum of weights of the points at the outer rim should be one. check
        // this
        for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
          Assert(std::fabs(std::accumulate(lohvs[unit_point].begin(),
                                           lohvs[unit_point].end(),0.) - 1)<1e-13*polynomial_degree,
                 ExcInternalError());

        return lohvs;
      }



      std::vector<::Table<2,double> >
      compute_support_point_weights_perimeter_to_interior(const unsigned int polynomial_degree,
                                                          const unsigned int dim)
      {
        Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
        std::vector<::Table<2,double> > output(dim);
        if (polynomial_degree <= 1)
          return output;

        // fill the 1D interior weights
        QGaussLobatto<1> quadrature(polynomial_degree+1);
        output[0].reinit(polynomial_degree-1, GeometryInfo<1>::vertices_per_cell);
        for (unsigned int q=0; q<polynomial_degree-1; ++q)
          for (unsigned int i=0; i<GeometryInfo<1>::vertices_per_cell; ++i)
            output[0](q,i) = GeometryInfo<1>::d_linear_shape_function(quadrature.point(q+1),
                                                                      i);

        if (dim > 1)
          output[1] = compute_support_point_weights_on_quad(polynomial_degree);

        if (dim > 2)
          output[2] = compute_support_point_weights_on_hex(polynomial_degree);

        return output;
      }

      template <int dim>
      ::Table<2,double>
      compute_support_point_weights_cell(const unsigned int polynomial_degree)
      {
        Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
        if (polynomial_degree <= 1)
          return ::Table<2,double>();

        QGaussLobatto<dim> quadrature(polynomial_degree+1);
        std::vector<unsigned int> h2l(quadrature.size());
        FETools::hierarchic_to_lexicographic_numbering<dim>(polynomial_degree, h2l);

        ::Table<2,double> output(quadrature.size() - GeometryInfo<dim>::vertices_per_cell,
                                       GeometryInfo<dim>::vertices_per_cell);
        for (unsigned int q=0; q<output.size(0); ++q)
          for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
            output(q,i) = GeometryInfo<dim>::d_linear_shape_function(quadrature.point(h2l[q+GeometryInfo<dim>::vertices_per_cell]),
                                                                     i);

        return output;
      }



      template <int dim, int spacedim>
      Point<spacedim>
      compute_mapped_location_of_point
      (const typename ::MappingQGeneric<dim,spacedim>::InternalData &data)
      {
        AssertDimension (data.shape_values.size(),
                         data.mapping_support_points.size());

        // use now the InternalData to compute the point in real space.
        Point<spacedim> p_real;
        for (unsigned int i=0; i<data.mapping_support_points.size(); ++i)
          p_real += data.mapping_support_points[i] * data.shape(0,i);

        return p_real;
      }



      template <int dim>
      Point<dim>
      do_transform_real_to_unit_cell_internal
      (const typename ::Triangulation<dim,dim>::cell_iterator &cell,
       const Point<dim>                                             &p,
       const Point<dim>                                             &initial_p_unit,
       typename ::MappingQGeneric<dim,dim>::InternalData      &mdata)
      {
        const unsigned int spacedim = dim;

        const unsigned int n_shapes=mdata.shape_values.size();
        (void)n_shapes;
        Assert(n_shapes!=0, ExcInternalError());
        AssertDimension (mdata.shape_derivatives.size(), n_shapes);

        std::vector<Point<spacedim> > &points=mdata.mapping_support_points;
        AssertDimension (points.size(), n_shapes);


        // Newton iteration to solve
        //    f(x)=p(x)-p=0
        // where we are looking for 'x' and p(x) is the forward transformation
        // from unit to real cell. We solve this using a Newton iteration
        //    x_{n+1}=x_n-[f'(x)]^{-1}f(x)
        // The start value is set to be the linear approximation to the cell

        // The shape values and derivatives of the mapping at this point are
        // previously computed.

        Point<dim> p_unit = initial_p_unit;

        mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit));

        Point<spacedim> p_real = compute_mapped_location_of_point<dim,spacedim>(mdata);
        Tensor<1,spacedim> f = p_real-p;

        // early out if we already have our point
        if (f.norm_square() < 1e-24 * cell->diameter() * cell->diameter())
          return p_unit;

        // we need to compare the position of the computed p(x) against the given
        // point 'p'. We will terminate the iteration and return 'x' if they are
        // less than eps apart. The question is how to choose eps -- or, put maybe
        // more generally: in which norm we want these 'p' and 'p(x)' to be eps
        // apart.
        //
        // the question is difficult since we may have to deal with very elongated
        // cells where we may achieve 1e-12*h for the distance of these two points
        // in the 'long' direction, but achieving this tolerance in the 'short'
        // direction of the cell may not be possible
        //
        // what we do instead is then to terminate iterations if
        //    \| p(x) - p \|_A < eps
        // where the A-norm is somehow induced by the transformation of the cell.
        // in particular, we want to measure distances relative to the sizes of
        // the cell in its principal directions.
        //
        // to define what exactly A should be, note that to first order we have
        // the following (assuming that x* is the solution of the problem, i.e.,
        // p(x*)=p):
        //    p(x) - p = p(x) - p(x*)
        //             = -grad p(x) * (x*-x) + higher order terms
        // This suggest to measure with a norm that corresponds to
        //    A = {[grad p(x]^T [grad p(x)]}^{-1}
        // because then
        //    \| p(x) - p \|_A  \approx  \| x - x* \|
        // Consequently, we will try to enforce that
        //    \| p(x) - p \|_A  =  \| f \|  <=  eps
        //
        // Note that using this norm is a bit dangerous since the norm changes
        // in every iteration (A isn't fixed by depends on xk). However, if the
        // cell is not too deformed (it may be stretched, but not twisted) then
        // the mapping is almost linear and A is indeed constant or nearly so.
        const double eps = 1.e-11;
        const unsigned int newton_iteration_limit = 20;

        unsigned int newton_iteration = 0;
        double last_f_weighted_norm;
        do
          {
#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
            std::cout << "Newton iteration " << newton_iteration << std::endl;
#endif

            // f'(x)
            Tensor<2,spacedim> df;
            for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
              {
                const Tensor<1,dim> &grad_transform=mdata.derivative(0,k);
                const Point<spacedim> &point=points[k];

                for (unsigned int i=0; i<spacedim; ++i)
                  for (unsigned int j=0; j<dim; ++j)
                    df[i][j]+=point[i]*grad_transform[j];
              }

            // Solve  [f'(x)]d=f(x)
            AssertThrow(determinant(df) > 0,
                        (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
            Tensor<2,spacedim> df_inverse = invert(df);
            const Tensor<1,spacedim> delta = df_inverse * static_cast<const Tensor<1,spacedim>&>(f);

#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
            std::cout << "   delta=" << delta  << std::endl;
#endif

            // do a line search
            double step_length = 1;
            do
              {
                // update of p_unit. The spacedim-th component of transformed point
                // is simply ignored in codimension one case. When this component is
                // not zero, then we are projecting the point to the surface or
                // curve identified by the cell.
                Point<dim> p_unit_trial = p_unit;
                for (unsigned int i=0; i<dim; ++i)
                  p_unit_trial[i] -= step_length * delta[i];

                // shape values and derivatives
                // at new p_unit point
                mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit_trial));

                // f(x)
                Point<spacedim> p_real_trial = internal::MappingQGenericImplementation::compute_mapped_location_of_point<dim,spacedim>(mdata);
                const Tensor<1,spacedim> f_trial = p_real_trial-p;

#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
                std::cout << "     step_length=" << step_length << std::endl
                          << "       ||f ||   =" << f.norm() << std::endl
                          << "       ||f*||   =" << f_trial.norm() << std::endl
                          << "       ||f*||_A =" << (df_inverse * f_trial).norm() << std::endl;
#endif

                // see if we are making progress with the current step length
                // and if not, reduce it by a factor of two and try again
                //
                // strictly speaking, we should probably use the same norm as we use
                // for the outer algorithm. in practice, line search is just a
                // crutch to find a "reasonable" step length, and so using the l2
                // norm is probably just fine
                if (f_trial.norm() < f.norm())
                  {
                    p_real = p_real_trial;
                    p_unit = p_unit_trial;
                    f = f_trial;
                    break;
                  }
                else if (step_length > 0.05)
                  step_length /= 2;
                else
                  AssertThrow (false,
                               (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
              }
            while (true);

            ++newton_iteration;
            if (newton_iteration > newton_iteration_limit)
              AssertThrow (false,
                           (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
            last_f_weighted_norm = (df_inverse * f).norm();
          }
        while (last_f_weighted_norm > eps);

        return p_unit;
      }



      template <int dim>
      Point<dim>
      do_transform_real_to_unit_cell_internal_codim1
      (const typename ::Triangulation<dim,dim+1>::cell_iterator &cell,
       const Point<dim+1>                                             &p,
       const Point<dim>                                               &initial_p_unit,
       typename ::MappingQGeneric<dim,dim+1>::InternalData      &mdata)
      {
        const unsigned int spacedim = dim+1;

        const unsigned int n_shapes=mdata.shape_values.size();
        (void)n_shapes;
        Assert(n_shapes!=0, ExcInternalError());
        Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError());
        Assert(mdata.shape_second_derivatives.size()==n_shapes, ExcInternalError());

        std::vector<Point<spacedim> > &points=mdata.mapping_support_points;
        Assert(points.size()==n_shapes, ExcInternalError());

        Point<spacedim> p_minus_F;

        Tensor<1,spacedim>  DF[dim];
        Tensor<1,spacedim>  D2F[dim][dim];

        Point<dim> p_unit = initial_p_unit;
        Point<dim> f;
        Tensor<2,dim>  df;

        // Evaluate first and second derivatives
        mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit));

        for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
          {
            const Tensor<1,dim>   &grad_phi_k = mdata.derivative(0,k);
            const Tensor<2,dim>   &hessian_k  = mdata.second_derivative(0,k);
            const Point<spacedim> &point_k = points[k];

            for (unsigned int j=0; j<dim; ++j)
              {
                DF[j] += grad_phi_k[j] * point_k;
                for (unsigned int l=0; l<dim; ++l)
                  D2F[j][l] += hessian_k[j][l] * point_k;
              }
          }

        p_minus_F = p;
        p_minus_F -= compute_mapped_location_of_point<dim,spacedim>(mdata);


        for (unsigned int j=0; j<dim; ++j)
          f[j] = DF[j] * p_minus_F;

        for (unsigned int j=0; j<dim; ++j)
          {
            f[j] = DF[j] * p_minus_F;
            for (unsigned int l=0; l<dim; ++l)
              df[j][l] = -DF[j]*DF[l] + D2F[j][l] * p_minus_F;
          }


        const double eps = 1.e-12*cell->diameter();
        const unsigned int loop_limit = 10;

        unsigned int loop=0;

        while (f.norm()>eps && loop++<loop_limit)
          {
            // Solve  [df(x)]d=f(x)
            const Tensor<1,dim> d = invert(df) * static_cast<const Tensor<1,dim>&>(f);
            p_unit -= d;

            for (unsigned int j=0; j<dim; ++j)
              {
                DF[j].clear();
                for (unsigned int l=0; l<dim; ++l)
                  D2F[j][l].clear();
              }

            mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit));

            for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
              {
                const Tensor<1,dim>   &grad_phi_k = mdata.derivative(0,k);
                const Tensor<2,dim>   &hessian_k  = mdata.second_derivative(0,k);
                const Point<spacedim> &point_k = points[k];

                for (unsigned int j=0; j<dim; ++j)
                  {
                    DF[j] += grad_phi_k[j] * point_k;
                    for (unsigned int l=0; l<dim; ++l)
                      D2F[j][l] += hessian_k[j][l] * point_k;
                  }
              }

            //TODO: implement a line search here in much the same way as for
            // the corresponding function above that does so for dim==spacedim
            p_minus_F = p;
            p_minus_F -= compute_mapped_location_of_point<dim,spacedim>(mdata);

            for (unsigned int j=0; j<dim; ++j)
              {
                f[j] = DF[j] * p_minus_F;
                for (unsigned int l=0; l<dim; ++l)
                  df[j][l] = -DF[j]*DF[l] + D2F[j][l] * p_minus_F;
              }

          }


        // Here we check that in the last execution of while the first
        // condition was already wrong, meaning the residual was below
        // eps. Only if the first condition failed, loop will have been
        // increased and tested, and thus have reached the limit.
        AssertThrow (loop<loop_limit, (typename Mapping<dim,spacedim>::ExcTransformationFailed()));

        return p_unit;
      }

      template <int dim, int spacedim>
      void
      maybe_update_q_points_Jacobians_and_grads_tensor
      (const CellSimilarity::Similarity                                    cell_similarity,
       const typename ::MappingQGeneric<dim,spacedim>::InternalData &data,
       std::vector<Point<spacedim> >                                      &quadrature_points,
       std::vector<DerivativeForm<2,dim,spacedim> >                       &jacobian_grads)
      {
        const UpdateFlags update_flags = data.update_each;

        const unsigned int n_shape_values = data.n_shape_functions;
        const unsigned int n_q_points = data.shape_info.n_q_points;
        const unsigned int vec_length = VectorizedArray<double>::n_array_elements;
        const unsigned int n_comp = 1+ (spacedim-1)/vec_length;
        const unsigned int n_hessians = (dim*(dim+1))/2;

        const bool evaluate_values = update_flags & update_quadrature_points;
        const bool evaluate_gradients= (cell_similarity != CellSimilarity::translation)
                                       &&(update_flags & update_contravariant_transformation);
        const bool evaluate_hessians = (cell_similarity != CellSimilarity::translation)
                                       &&(update_flags & update_jacobian_grads);

        Assert (!evaluate_values || n_q_points > 0, ExcInternalError());
        Assert (!evaluate_values || n_q_points == quadrature_points.size(),
                ExcDimensionMismatch(n_q_points, quadrature_points.size()));
        Assert (!evaluate_gradients || data.n_shape_functions > 0, ExcInternalError());
        Assert (!evaluate_gradients || n_q_points == data.contravariant.size(),
                ExcDimensionMismatch(n_q_points, data.contravariant.size()));
        Assert (!evaluate_hessians || n_q_points == jacobian_grads.size(),
                ExcDimensionMismatch(n_q_points, jacobian_grads.size()));

        // prepare arrays
        if (evaluate_values || evaluate_gradients || evaluate_hessians)
          {
            data.values_dofs.resize(n_comp*n_shape_values);
            data.values_quad.resize(n_comp*n_q_points);
            data.gradients_quad.resize (n_comp*n_q_points*dim);

            if (evaluate_hessians)
              data.hessians_quad.resize(n_comp*n_q_points*n_hessians);

            const std::vector<unsigned int> &renumber_to_lexicographic
              = data.shape_info.lexicographic_numbering;
            for (unsigned int i=0; i<n_shape_values; ++i)
              for (unsigned int d=0; d<spacedim; ++d)
                {
                  const unsigned int in_comp = d%vec_length;
                  const unsigned int out_comp = d/vec_length;
                  data.values_dofs[out_comp*n_shape_values+i][in_comp]
                    = data.mapping_support_points[renumber_to_lexicographic[i]][d];
                }

            // do the actual tensorized evaluation
            SelectEvaluator<dim, -1, 0, n_comp, VectorizedArray<double> >::evaluate
            (data.shape_info, data.values_dofs.begin(), data.values_quad.begin(),
             data.gradients_quad.begin(), data.hessians_quad.begin(), data.scratch.begin(),
             evaluate_values, evaluate_gradients, evaluate_hessians);
          }

        // do the postprocessing
        if (evaluate_values)
          {
            for (unsigned int out_comp=0; out_comp<n_comp; ++out_comp)
              for (unsigned int i=0; i<n_q_points; ++i)
                for (unsigned int in_comp=0;
                     in_comp<vec_length && in_comp<spacedim-out_comp*vec_length; ++in_comp)
                  quadrature_points[i][out_comp*vec_length+in_comp]
                    = data.values_quad[out_comp*n_q_points+i][in_comp];
          }

        if (evaluate_gradients)
          {
            std::fill(data.contravariant.begin(), data.contravariant.end(),
                      DerivativeForm<1,dim,spacedim>());
            // We need to reinterpret the data after evaluate has been applied.
            for (unsigned int out_comp=0; out_comp<n_comp; ++out_comp)
              for (unsigned int point=0; point<n_q_points; ++point)
                for (unsigned int j=0; j<dim; ++j)
                  for (unsigned int in_comp=0;
                       in_comp<vec_length && in_comp<spacedim-out_comp*vec_length; ++in_comp)
                    {
                      const unsigned int total_number = point*dim+j;
                      const unsigned int new_comp = total_number/n_q_points;
                      const unsigned int new_point = total_number % n_q_points;
                      data.contravariant[new_point][out_comp*vec_length+in_comp][new_comp]
                        = data.gradients_quad[(out_comp*n_q_points+point)*dim+j][in_comp];
                    }
          }
        if (update_flags & update_covariant_transformation)
          if (cell_similarity != CellSimilarity::translation)
            for (unsigned int point=0; point<n_q_points; ++point)
              data.covariant[point] = (data.contravariant[point]).covariant_form();

        if (update_flags & update_volume_elements)
          if (cell_similarity != CellSimilarity::translation)
            for (unsigned int point=0; point<n_q_points; ++point)
              data.volume_elements[point] = data.contravariant[point].determinant();

        if (evaluate_hessians)
          {
            constexpr int desymmetrize_3d [6][2] = {{0,0},{1,1},{2,2},{0,1},{0,2},{1,2}};
            constexpr int desymmetrize_2d [3][2] = {{0,0},{1,1},{0,1}};

            // We need to reinterpret the data after evaluate has been applied.
            for (unsigned int out_comp=0; out_comp<n_comp; ++out_comp)
              for (unsigned int point=0; point<n_q_points; ++point)
                for (unsigned int j=0; j<n_hessians; ++j)
                  for (unsigned int in_comp=0;
                       in_comp<vec_length && in_comp<spacedim-out_comp*vec_length; ++in_comp)
                    {
                      const unsigned int total_number = point*n_hessians+j;
                      const unsigned int new_point = total_number % n_q_points;
                      const unsigned int new_hessian_comp = total_number/n_q_points;
                      const unsigned int new_hessian_comp_i = dim==2 ? desymmetrize_2d[new_hessian_comp][0]
                                                              : desymmetrize_3d[new_hessian_comp][0];
                      const unsigned int new_hessian_comp_j = dim==2 ? desymmetrize_2d[new_hessian_comp][1]
                                                              : desymmetrize_3d[new_hessian_comp][1];
                      const double value = data.hessians_quad[(out_comp*n_q_points+point)*n_hessians+j][in_comp];
                      jacobian_grads[new_point][out_comp*vec_length+in_comp][new_hessian_comp_i][new_hessian_comp_j] = value;
                      jacobian_grads[new_point][out_comp*vec_length+in_comp][new_hessian_comp_j][new_hessian_comp_i] = value;
                    }
          }
      }


      template <int dim, int spacedim>
      void
      maybe_compute_q_points
      (const typename QProjector<dim>::DataSetDescriptor                   data_set,
       const typename ::MappingQGeneric<dim,spacedim>::InternalData &data,
       std::vector<Point<spacedim> >                                      &quadrature_points)
      {
        const UpdateFlags update_flags = data.update_each;

        if (update_flags & update_quadrature_points)
          for (unsigned int point=0; point<quadrature_points.size(); ++point)
            {
              const double *shape = &data.shape(point+data_set,0);
              Point<spacedim> result = (shape[0] *
                                        data.mapping_support_points[0]);
              for (unsigned int k=1; k<data.n_shape_functions; ++k)
                for (unsigned int i=0; i<spacedim; ++i)
                  result[i] += shape[k] * data.mapping_support_points[k][i];
              quadrature_points[point] = result;
            }
      }



      template <int dim, int spacedim>
      void
      maybe_update_Jacobians
      (const CellSimilarity::Similarity                                    cell_similarity,
       const typename ::QProjector<dim>::DataSetDescriptor           data_set,
       const typename ::MappingQGeneric<dim,spacedim>::InternalData &data)
      {
        const UpdateFlags update_flags = data.update_each;

        if (update_flags & update_contravariant_transformation)
          // if the current cell is just a
          // translation of the previous one, no
          // need to recompute jacobians...
          if (cell_similarity != CellSimilarity::translation)
            {
              const unsigned int n_q_points = data.contravariant.size();

              std::fill(data.contravariant.begin(), data.contravariant.end(),
                        DerivativeForm<1,dim,spacedim>());

              Assert (data.n_shape_functions > 0, ExcInternalError());

              const Tensor<1,spacedim> *supp_pts =
                &data.mapping_support_points[0];

              for (unsigned int point=0; point<n_q_points; ++point)
                {
                  const Tensor<1,dim> *data_derv =
                    &data.derivative(point+data_set, 0);

                  double result [spacedim][dim];

                  // peel away part of sum to avoid zeroing the
                  // entries and adding for the first time
                  for (unsigned int i=0; i<spacedim; ++i)
                    for (unsigned int j=0; j<dim; ++j)
                      result[i][j] = data_derv[0][j] * supp_pts[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)
                        result[i][j] += data_derv[k][j] * supp_pts[k][i];

                  // write result into contravariant data. for
                  // j=dim in the case dim<spacedim, there will
                  // never be any nonzero data that arrives in
                  // here, so it is ok anyway because it was
                  // initialized to zero at the initialization
                  for (unsigned int i=0; i<spacedim; ++i)
                    for (unsigned int j=0; j<dim; ++j)
                      data.contravariant[point][i][j] = result[i][j];
                }
            }

        if (update_flags & update_covariant_transformation)
          if (cell_similarity != CellSimilarity::translation)
            {
              const unsigned int n_q_points = data.contravariant.size();
              for (unsigned int point=0; point<n_q_points; ++point)
                {
                  data.covariant[point] = (data.contravariant[point]).covariant_form();
                }
            }

        if (update_flags & update_volume_elements)
          if (cell_similarity != CellSimilarity::translation)
            {
              const unsigned int n_q_points = data.contravariant.size();
              for (unsigned int point=0; point<n_q_points; ++point)
                data.volume_elements[point] = data.contravariant[point].determinant();
            }

      }

      template <int dim, int spacedim>
      void
      maybe_update_jacobian_grads
      (const CellSimilarity::Similarity                                    cell_similarity,
       const typename QProjector<dim>::DataSetDescriptor                   data_set,
       const typename ::MappingQGeneric<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];
                }
          }
      }

      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 ::MappingQGeneric<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];
                              }

                          }
                  }
              }
          }
      }

      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 ::MappingQGeneric<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];
                  }
              }
          }
      }

      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 ::MappingQGeneric<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];
                            }
                  }
              }
          }
      }

      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 ::MappingQGeneric<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];
                  }
              }
          }
      }

      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 ::MappingQGeneric<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>
MappingQGeneric<dim,spacedim>::MappingQGeneric (const unsigned int p)
  :
  polynomial_degree(p),
  line_support_points(this->polynomial_degree+1),
  fe_q(dim == 3 ? new FE_Q<dim>(this->polynomial_degree) : nullptr),
  support_point_weights_perimeter_to_interior (internal::MappingQGenericImplementation::compute_support_point_weights_perimeter_to_interior(this->polynomial_degree, dim)),
  support_point_weights_cell (internal::MappingQGenericImplementation::compute_support_point_weights_cell<dim>(this->polynomial_degree))
{
  Assert (p >= 1, ExcMessage ("It only makes sense to create polynomial mappings "
                              "with a polynomial degree greater or equal to one."));
}



template <int dim, int spacedim>
MappingQGeneric<dim,spacedim>::MappingQGeneric (const MappingQGeneric<dim,spacedim> &mapping)
  :
  polynomial_degree(mapping.polynomial_degree),
  line_support_points(mapping.line_support_points),
  fe_q(dim == 3 ? new FE_Q<dim>(*mapping.fe_q) : nullptr),
  support_point_weights_perimeter_to_interior (mapping.support_point_weights_perimeter_to_interior),
  support_point_weights_cell (mapping.support_point_weights_cell)
{}




template <int dim, int spacedim>
std::unique_ptr<Mapping<dim,spacedim> >
MappingQGeneric<dim,spacedim>::clone () const
{
  return std_cxx14::make_unique<MappingQGeneric<dim,spacedim>>(*this);
}




template <int dim, int spacedim>
unsigned int
MappingQGeneric<dim,spacedim>::get_degree() const
{
  return polynomial_degree;
}



template <int dim, int spacedim>
Point<spacedim>
MappingQGeneric<dim,spacedim>::
transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                             const Point<dim> &p) const
{
  // set up the polynomial space
  const TensorProductPolynomials<dim>
  tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points()));
  Assert (tensor_pols.n() == Utilities::fixed_power<dim>(polynomial_degree+1),
          ExcInternalError());

  // then also construct the mapping from lexicographic to the Qp shape function numbering
  const std::vector<unsigned int>
  renumber (FETools::
            lexicographic_to_hierarchic_numbering
            (FiniteElementData<dim> (internal::MappingQGenericImplementation::get_dpo_vector<dim>
                                     (polynomial_degree), 1, polynomial_degree)));

  const std::vector<Point<spacedim> > support_points
    = this->compute_mapping_support_points(cell);

  Point<spacedim> mapped_point;
  for (unsigned int i=0; i<tensor_pols.n(); ++i)
    mapped_point += support_points[renumber[i]] * tensor_pols.compute_value (i, p);

  return mapped_point;
}


// In the code below, GCC tries to instantiate MappingQGeneric<3,4> when
// seeing which of the overloaded versions of
// do_transform_real_to_unit_cell_internal() to call. This leads to bad
// error messages and, generally, nothing very good. Avoid this by ensuring
// that this class exists, but does not have an inner InternalData
// type, thereby ruling out the codim-1 version of the function
// below when doing overload resolution.
template <>
class MappingQGeneric<3,4>
{};



// visual studio freaks out when trying to determine if
// do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
// candidate. So instead of letting the compiler pick the correct overload, we
// use template specialization to make sure we pick up the right function to
// call:

template <int dim, int spacedim>
Point<dim>
MappingQGeneric<dim,spacedim>::
transform_real_to_unit_cell_internal
(const typename Triangulation<dim,spacedim>::cell_iterator &,
 const Point<spacedim> &,
 const Point<dim> &) const
{
  // default implementation (should never be called)
  Assert(false, ExcInternalError());
  return Point<dim>();
}

template <>
Point<1>
MappingQGeneric<1,1>::
transform_real_to_unit_cell_internal
(const Triangulation<1,1>::cell_iterator &cell,
 const Point<1>                          &p,
 const Point<1>                          &initial_p_unit) const
{
  const int dim = 1;
  const int spacedim = 1;

  const Quadrature<dim> point_quadrature(initial_p_unit);

  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
  if (spacedim>dim)
    update_flags |= update_jacobian_grads;
  auto mdata = Utilities::dynamic_unique_cast<InternalData>(get_data(update_flags,
                                                            point_quadrature));

  mdata->mapping_support_points = this->compute_mapping_support_points (cell);

  // dispatch to the various specializations for spacedim=dim,
  // spacedim=dim+1, etc
  return internal::MappingQGenericImplementation::do_transform_real_to_unit_cell_internal<1>(cell, p, initial_p_unit, *mdata);
}

template <>
Point<2>
MappingQGeneric<2, 2>::
transform_real_to_unit_cell_internal
(const Triangulation<2, 2>::cell_iterator &cell,
 const Point<2>                           &p,
 const Point<2>                           &initial_p_unit) const
{
  const int dim = 2;
  const int spacedim = 2;

  const Quadrature<dim> point_quadrature(initial_p_unit);

  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
  if (spacedim>dim)
    update_flags |= update_jacobian_grads;
  auto mdata = Utilities::dynamic_unique_cast<InternalData>(get_data(update_flags,
                                                            point_quadrature));

  mdata->mapping_support_points = this->compute_mapping_support_points (cell);

  // dispatch to the various specializations for spacedim=dim,
  // spacedim=dim+1, etc
  return internal::MappingQGenericImplementation::do_transform_real_to_unit_cell_internal<2>(cell, p, initial_p_unit, *mdata);
}

template <>
Point<3>
MappingQGeneric<3, 3>::
transform_real_to_unit_cell_internal
(const Triangulation<3, 3>::cell_iterator &cell,
 const Point<3>                            &p,
 const Point<3>                                 &initial_p_unit) const
{
  const int dim = 3;
  const int spacedim = 3;

  const Quadrature<dim> point_quadrature(initial_p_unit);

  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
  if (spacedim>dim)
    update_flags |= update_jacobian_grads;
  auto mdata = Utilities::dynamic_unique_cast<InternalData>(get_data(update_flags,
                                                            point_quadrature));

  mdata->mapping_support_points = this->compute_mapping_support_points (cell);

  // dispatch to the various specializations for spacedim=dim,
  // spacedim=dim+1, etc
  return internal::MappingQGenericImplementation::do_transform_real_to_unit_cell_internal<3>(cell, p, initial_p_unit, *mdata);
}



template <>
Point<1>
MappingQGeneric<1, 2>::
transform_real_to_unit_cell_internal
(const Triangulation<1, 2>::cell_iterator &cell,
 const Point<2>                            &p,
 const Point<1>                                 &initial_p_unit) const
{
  const int dim = 1;
  const int spacedim = 2;

  const Quadrature<dim> point_quadrature(initial_p_unit);

  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
  if (spacedim>dim)
    update_flags |= update_jacobian_grads;
  auto mdata = Utilities::dynamic_unique_cast<InternalData>(get_data(update_flags,
                                                            point_quadrature));

  mdata->mapping_support_points = this->compute_mapping_support_points (cell);

  // dispatch to the various specializations for spacedim=dim,
  // spacedim=dim+1, etc
  return internal::MappingQGenericImplementation::do_transform_real_to_unit_cell_internal_codim1<1>(cell, p, initial_p_unit, *mdata);
}



template <>
Point<2>
MappingQGeneric<2, 3>::
transform_real_to_unit_cell_internal
(const Triangulation<2, 3>::cell_iterator &cell,
 const Point<3>                            &p,
 const Point<2>                                 &initial_p_unit) const
{
  const int dim = 2;
  const int spacedim = 3;

  const Quadrature<dim> point_quadrature(initial_p_unit);

  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
  if (spacedim>dim)
    update_flags |= update_jacobian_grads;
  auto mdata = Utilities::dynamic_unique_cast<InternalData>(get_data(update_flags,
                                                            point_quadrature));

  mdata->mapping_support_points = this->compute_mapping_support_points (cell);

  // dispatch to the various specializations for spacedim=dim,
  // spacedim=dim+1, etc
  return internal::MappingQGenericImplementation::do_transform_real_to_unit_cell_internal_codim1<2>(cell, p, initial_p_unit, *mdata);
}

template <>
Point<1>
MappingQGeneric<1, 3>::
transform_real_to_unit_cell_internal
(const Triangulation<1, 3>::cell_iterator &,
 const Point<3> &,
 const Point<1> &) const
{
  Assert (false, ExcNotImplemented());
  return Point<1>();
}



template <int dim, int spacedim>
Point<dim>
MappingQGeneric<dim,spacedim>::
transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                             const Point<spacedim>                            &p) const
{
  // Use an exact formula if one is available. this is only the case
  // for Q1 mappings in 1d, and in 2d if dim==spacedim
  if ((polynomial_degree == 1) &&
      ((dim == 1)
       ||
       ((dim == 2) && (dim == spacedim))))
    {
      // The dimension-dependent algorithms are much faster (about 25-45x in
      // 2D) but fail most of the time when the given point (p) is not in the
      // cell. The dimension-independent Newton algorithm given below is
      // slower, but more robust (though it still sometimes fails). Therefore
      // this function implements the following strategy based on the
      // p's dimension:
      //
      // * In 1D this mapping is linear, so the mapping is always invertible
      //   (and the exact formula is known) as long as the cell has non-zero
      //   length.
      // * In 2D the exact (quadratic) formula is called first. If either the
      //   exact formula does not succeed (negative discriminant in the
      //   quadratic formula) or succeeds but finds a solution outside of the
      //   unit cell, then the Newton solver is called. The rationale for the
      //   second choice is that the exact formula may provide two different
      //   answers when mapping a point outside of the real cell, but the
      //   Newton solver (if it converges) will only return one answer.
      //   Otherwise the exact formula successfully found a point in the unit
      //   cell and that value is returned.
      // * In 3D there is no (known to the authors) exact formula, so the Newton
      //   algorithm is used.
      const std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
      vertices = this->get_vertices(cell);
      try
        {
          switch (dim)
            {
            case 1:
            {
              // formula not subject to any issues in 1d
              if (spacedim == 1)
                return internal::MappingQ1::transform_real_to_unit_cell(vertices, p);
              else
                break;
            }

            case 2:
            {
              const Point<dim> point
                = internal::MappingQ1::transform_real_to_unit_cell(vertices, p);

              // formula not guaranteed to work for points outside of
              // the cell. only take the computed point if it lies
              // inside the reference cell
              const double eps = 1e-15;
              if (-eps <= point(1) && point(1) <= 1 + eps &&
                  -eps <= point(0) && point(0) <= 1 + eps)
                {
                  return point;
                }
              else
                break;
            }

            default:
            {
              // we should get here, based on the if-condition at the top
              Assert(false, ExcInternalError());
            }
            }
        }
      catch (const typename Mapping<spacedim,spacedim>::ExcTransformationFailed &)
        {
          // simply fall through and continue on to the standard Newton code
        }
    }
  else
    {
      // we can't use an explicit formula,
    }


  // Find the initial value for the Newton iteration by a normal
  // projection to the least square plane determined by the vertices
  // of the cell
  Point<dim> initial_p_unit;
  if (this->preserves_vertex_locations())
    initial_p_unit = cell->real_to_unit_cell_affine_approximation(p);
  else
    {
      // for the MappingQEulerian type classes, we want to still call the cell
      // iterator's affine approximation. do so by creating a dummy
      // triangulation with just the first vertices.
      //
      // we do this by first getting all support points, then
      // throwing away all but the vertices, and finally calling
      // the same function as above
      std::vector<Point<spacedim> > a
        = this->compute_mapping_support_points (cell);
      a.resize(GeometryInfo<dim>::vertices_per_cell);
      std::vector<CellData<dim> > cells(1);
      for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
        cells[0].vertices[i] = i;
      Triangulation<dim,spacedim> tria;
      tria.create_triangulation(a, cells, SubCellData());
      initial_p_unit = tria.begin_active()->real_to_unit_cell_affine_approximation(p);
    }
  // in 1d with spacedim > 1 the affine approximation is exact
  if (dim == 1 && polynomial_degree == 1)
    {
      return initial_p_unit;
    }
  else
    {
      // in case the function above should have given us something back that
      // lies outside the unit cell, then project it back into the reference
      // cell in hopes that this gives a better starting point to the
      // following iteration
      initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);

      // perform the Newton iteration and return the result. note that this
      // statement may throw an exception, which we simply pass up to the
      // caller
      return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
    }
}



template <int dim, int spacedim>
UpdateFlags
MappingQGeneric<dim,spacedim>::requires_update_flags (const UpdateFlags in) const
{
  // add flags if the respective quantities are necessary to compute
  // what we need. note that some flags appear in both the conditions
  // and in subsequent set operations. this leads to some circular
  // logic. the only way to treat this is to iterate. since there are
  // 5 if-clauses in the loop, it will take at most 5 iterations to
  // converge. do them:
  UpdateFlags out = in;
  for (unsigned int i=0; i<5; ++i)
    {
      // The following is a little incorrect:
      // If not applied on a face,
      // update_boundary_forms does not
      // make sense. On the other hand,
      // it is necessary on a
      // face. Currently,
      // update_boundary_forms is simply
      // ignored for the interior of a
      // cell.
      if (out & (update_JxW_values
                 | update_normal_vectors))
        out |= update_boundary_forms;

      if (out & (update_covariant_transformation
                 | update_JxW_values
                 | update_jacobians
                 | update_jacobian_grads
                 | update_boundary_forms
                 | update_normal_vectors))
        out |= update_contravariant_transformation;

      if (out & (update_inverse_jacobians
                 | update_jacobian_pushed_forward_grads
                 | update_jacobian_pushed_forward_2nd_derivatives
                 | update_jacobian_pushed_forward_3rd_derivatives) )
        out |= update_covariant_transformation;

      // The contravariant transformation is used in the Piola
      // transformation, which requires the determinant of the Jacobi
      // matrix of the transformation.  Because we have no way of
      // knowing here whether the finite element wants to use the
      // contravariant or the Piola transforms, we add the JxW values
      // to the list of flags to be updated for each cell.
      if (out & update_contravariant_transformation)
        out |= update_volume_elements;

      // the same is true when computing normal vectors: they require
      // the determinant of the Jacobian
      if (out & update_normal_vectors)
        out |= update_volume_elements;
    }

  return out;
}



template <int dim, int spacedim>
std::unique_ptr<typename Mapping<dim,spacedim>::InternalDataBase>
MappingQGeneric<dim,spacedim>::get_data (const UpdateFlags update_flags,
                                         const Quadrature<dim> &q) const
{
  auto data = std_cxx14::make_unique<InternalData>(polynomial_degree);
  data->initialize (this->requires_update_flags(update_flags), q, q.size());

  return std::move(data);
}



template <int dim, int spacedim>
std::unique_ptr<typename Mapping<dim,spacedim>::InternalDataBase>
MappingQGeneric<dim,spacedim>::get_face_data (const UpdateFlags        update_flags,
                                              const Quadrature<dim-1> &quadrature) const
{
  auto data = std_cxx14::make_unique<InternalData>(polynomial_degree);
  data->initialize_face (this->requires_update_flags(update_flags),
                         QProjector<dim>::project_to_all_faces(quadrature),
                         quadrature.size());

  return std::move(data);
}



template <int dim, int spacedim>
std::unique_ptr<typename Mapping<dim,spacedim>::InternalDataBase>
MappingQGeneric<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
                                                 const Quadrature<dim-1>& quadrature) const
{
  auto data = std_cxx14::make_unique<InternalData>(polynomial_degree);
  data->initialize_face (this->requires_update_flags(update_flags),
                         QProjector<dim>::project_to_all_subfaces(quadrature),
                         quadrature.size());

  return std::move(data);
}



template <int dim, int spacedim>
CellSimilarity::Similarity
MappingQGeneric<dim,spacedim>::
fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                const CellSimilarity::Similarity                           cell_similarity,
                const Quadrature<dim>                                     &quadrature,
                const typename Mapping<dim,spacedim>::InternalDataBase    &internal_data,
                internal::FEValuesImplementation::MappingRelatedData<dim,spacedim>      &output_data) const
{
  // ensure that the following static_cast is really correct:
  Assert (dynamic_cast<const InternalData *>(&internal_data) != nullptr,
          ExcInternalError());
  const InternalData &data = static_cast<const InternalData &>(internal_data);

  const unsigned int n_q_points=quadrature.size();

  // recompute the support points of the transformation of this
  // cell. we tried to be clever here in an earlier version of the
  // library by checking whether the cell is the same as the one we
  // had visited last, but it turns out to be difficult to determine
  // that because a cell for the purposes of a mapping is
  // characterized not just by its (triangulation, level, index)
  // triple, but also by the locations of its vertices, the manifold
  // object attached to the cell and all of its bounding faces/edges,
  // etc. to reliably test that the "cell" we are on is, therefore,
  // not easily done
  data.mapping_support_points = this->compute_mapping_support_points(cell);
  data.cell_of_current_support_points = cell;

  // if the order of the mapping is greater than 1, then do not reuse any cell
  // similarity information. This is necessary because the cell similarity
  // value is computed with just cell vertices and does not take into account
  // cell curvature.
  const CellSimilarity::Similarity computed_cell_similarity =
    (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);

  if (dim>1 && data.tensor_product_quadrature)
    {
      internal::MappingQGenericImplementation::maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>
      (computed_cell_similarity,
       data,
       output_data.quadrature_points,
       output_data.jacobian_grads);
    }
  else
    {
      internal::MappingQGenericImplementation::maybe_compute_q_points<dim,spacedim>
      (QProjector<dim>::DataSetDescriptor::cell (),
       data,
       output_data.quadrature_points);

      internal::MappingQGenericImplementation::maybe_update_Jacobians<dim,spacedim>
      (computed_cell_similarity,
       QProjector<dim>::DataSetDescriptor::cell (),
       data);

      internal::MappingQGenericImplementation::maybe_update_jacobian_grads<dim,spacedim>
      (computed_cell_similarity,
       QProjector<dim>::DataSetDescriptor::cell (),
       data,
       output_data.jacobian_grads);
    }

  internal::MappingQGenericImplementation::maybe_update_jacobian_pushed_forward_grads<dim,spacedim>
  (computed_cell_similarity,
   QProjector<dim>::DataSetDescriptor::cell (),
   data,
   output_data.jacobian_pushed_forward_grads);

  internal::MappingQGenericImplementation::maybe_update_jacobian_2nd_derivatives<dim,spacedim>
  (computed_cell_similarity,
   QProjector<dim>::DataSetDescriptor::cell (),
   data,
   output_data.jacobian_2nd_derivatives);

  internal::MappingQGenericImplementation::maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,spacedim>
  (computed_cell_similarity,
   QProjector<dim>::DataSetDescriptor::cell (),
   data,
   output_data.jacobian_pushed_forward_2nd_derivatives);

  internal::MappingQGenericImplementation::maybe_update_jacobian_3rd_derivatives<dim,spacedim>
  (computed_cell_similarity,
   QProjector<dim>::DataSetDescriptor::cell (),
   data,
   output_data.jacobian_3rd_derivatives);

  internal::MappingQGenericImplementation::maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,spacedim>
  (computed_cell_similarity,
   QProjector<dim>::DataSetDescriptor::cell (),
   data,
   output_data.jacobian_pushed_forward_3rd_derivatives);

  const UpdateFlags update_flags = data.update_each;
  const std::vector<double> &weights=quadrature.get_weights();

  // Multiply quadrature weights by absolute value of Jacobian determinants or
  // the area element g=sqrt(DX^t DX) in case of codim > 0

  if (update_flags & (update_normal_vectors
                      | update_JxW_values))
    {
      AssertDimension (output_data.JxW_values.size(), n_q_points);

      Assert( !(update_flags & update_normal_vectors ) ||
              (output_data.normal_vectors.size() == n_q_points),
              ExcDimensionMismatch(output_data.normal_vectors.size(), n_q_points));


      if (computed_cell_similarity != CellSimilarity::translation)
        for (unsigned int point=0; point<n_q_points; ++point)
          {

            if (dim == spacedim)
              {
                const double det = data.contravariant[point].determinant();

                // check for distorted cells.

                // TODO: this allows for anisotropies of up to 1e6 in 3D and
                // 1e12 in 2D. might want to find a finer
                // (dimension-independent) criterion
                Assert (det > 1e-12*Utilities::fixed_power<dim>(cell->diameter()/
                                                                std::sqrt(double(dim))),
                        (typename Mapping<dim,spacedim>::ExcDistortedMappedCell(cell->center(), det, point)));

                output_data.JxW_values[point] = weights[point] * det;
              }
            // if dim==spacedim, then there is no cell normal to
            // compute. since this is for FEValues (and not FEFaceValues),
            // there are also no face normals to compute
            else //codim>0 case
              {
                Tensor<1, spacedim> DX_t [dim];
                for (unsigned int i=0; i<spacedim; ++i)
                  for (unsigned int j=0; j<dim; ++j)
                    DX_t[j][i] = data.contravariant[point][i][j];

                Tensor<2, dim> G; //First fundamental form
                for (unsigned int i=0; i<dim; ++i)
                  for (unsigned int j=0; j<dim; ++j)
                    G[i][j] = DX_t[i] * DX_t[j];

                output_data.JxW_values[point]
                  = sqrt(determinant(G)) * weights[point];

                if (computed_cell_similarity == CellSimilarity::inverted_translation)
                  {
                    // we only need to flip the normal
                    if (update_flags & update_normal_vectors)
                      output_data.normal_vectors[point] *= -1.;
                  }
                else
                  {
                    if (update_flags & update_normal_vectors)
                      {
                        Assert(spacedim == dim+1,
                               ExcMessage("There is no (unique) cell normal for "
                                          + Utilities::int_to_string(dim) +
                                          "-dimensional cells in "
                                          + Utilities::int_to_string(spacedim) +
                                          "-dimensional space. This only works if the "
                                          "space dimension is one greater than the "
                                          "dimensionality of the mesh cells."));

                        if (dim==1)
                          output_data.normal_vectors[point] =
                            cross_product_2d(-DX_t[0]);
                        else //dim == 2
                          output_data.normal_vectors[point] =
                            cross_product_3d(DX_t[0], DX_t[1]);

                        output_data.normal_vectors[point] /= output_data.normal_vectors[point].norm();

                        if (cell->direction_flag() == false)
                          output_data.normal_vectors[point] *= -1.;
                      }

                  }
              } //codim>0 case

          }
    }



  // copy values from InternalData to vector given by reference
  if (update_flags & update_jacobians)
    {
      AssertDimension (output_data.jacobians.size(), n_q_points);
      if (computed_cell_similarity != CellSimilarity::translation)
        for (unsigned int point=0; point<n_q_points; ++point)
          output_data.jacobians[point] = data.contravariant[point];
    }

  // copy values from InternalData to vector given by reference
  if (update_flags & update_inverse_jacobians)
    {
      AssertDimension (output_data.inverse_jacobians.size(), n_q_points);
      if (computed_cell_similarity != CellSimilarity::translation)
        for (unsigned int point=0; point<n_q_points; ++point)
          output_data.inverse_jacobians[point] = data.covariant[point].transpose();
    }

  return computed_cell_similarity;
}






namespace internal
{
  namespace MappingQGenericImplementation
  {
    namespace
    {
      template <int dim, int spacedim>
      void
      maybe_compute_face_data (const ::MappingQGeneric<dim,spacedim> &mapping,
                               const typename ::Triangulation<dim,spacedim>::cell_iterator &cell,
                               const unsigned int               face_no,
                               const unsigned int               subface_no,
                               const unsigned int               n_q_points,
                               const std::vector<double>        &weights,
                               const typename ::MappingQGeneric<dim,spacedim>::InternalData &data,
                               internal::FEValuesImplementation::MappingRelatedData<dim,spacedim>         &output_data)
      {
        const UpdateFlags update_flags = data.update_each;

        if (update_flags & (update_boundary_forms |
                            update_normal_vectors |
                            update_jacobians      |
                            update_JxW_values      |
                            update_inverse_jacobians))
          {
            if (update_flags & update_boundary_forms)
              AssertDimension (output_data.boundary_forms.size(), n_q_points);
            if (update_flags & update_normal_vectors)
              AssertDimension (output_data.normal_vectors.size(), n_q_points);
            if (update_flags & update_JxW_values)
              AssertDimension (output_data.JxW_values.size(), n_q_points);

            Assert (data.aux.size()+1 >= dim, ExcInternalError());

            // first compute some common data that is used for evaluating
            // all of the flags below

            // map the unit tangentials to the real cell. checking for d!=dim-1
            // eliminates compiler warnings regarding unsigned int expressions <
            // 0.
            for (unsigned int d=0; d!=dim-1; ++d)
              {
                Assert (face_no+GeometryInfo<dim>::faces_per_cell*d <
                        data.unit_tangentials.size(),
                        ExcInternalError());
                Assert (data.aux[d].size() <=
                        data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d].size(),
                        ExcInternalError());

                mapping.transform (make_array_view(data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d]),
                                   mapping_contravariant,
                                   data,
                                   make_array_view(data.aux[d]));
              }

            if (update_flags & update_boundary_forms)
              {
                // if dim==spacedim, we can use the unit tangentials to compute the
                // boundary form by simply taking the cross product
                if (dim == spacedim)
                  {
                    for (unsigned int i=0; i<n_q_points; ++i)
                      switch (dim)
                        {
                        case 1:
                          // in 1d, we don't have access to any of the data.aux
                          // fields (because it has only dim-1 components), but we
                          // can still compute the boundary form by simply
                          // looking at the number of the face
                          output_data.boundary_forms[i][0] = (face_no == 0 ?
                                                              -1 : +1);
                          break;
                        case 2:
                          output_data.boundary_forms[i] =
                            cross_product_2d(data.aux[0][i]);
                          break;
                        case 3:
                          output_data.boundary_forms[i] =
                            cross_product_3d(data.aux[0][i], data.aux[1][i]);
                          break;
                        default:
                          Assert(false, ExcNotImplemented());
                        }
                  }
                else //(dim < spacedim)
                  {
                    // in the codim-one case, the boundary form results from the
                    // cross product of all the face tangential vectors and the cell
                    // normal vector
                    //
                    // to compute the cell normal, use the same method used in
                    // fill_fe_values for cells above
                    AssertDimension (data.contravariant.size(), n_q_points);

                    for (unsigned int point=0; point<n_q_points; ++point)
                      {
                        if (dim==1)
                          {
                            // J is a tangent vector
                            output_data.boundary_forms[point] = data.contravariant[point].transpose()[0];
                            output_data.boundary_forms[point] /=
                              (face_no == 0 ? -1. : +1.) * output_data.boundary_forms[point].norm();
                          }

                        if (dim==2)
                          {
                            const DerivativeForm<1,spacedim,dim> DX_t =
                              data.contravariant[point].transpose();

                            Tensor<1, spacedim> cell_normal =
                              cross_product_3d(DX_t[0], DX_t[1]);
                            cell_normal /= cell_normal.norm();

                            // then compute the face normal from the face tangent
                            // and the cell normal:
                            output_data.boundary_forms[point] =
                              cross_product_3d(data.aux[0][point], cell_normal);
                          }
                      }
                  }
              }

            if (update_flags & update_JxW_values)
              for (unsigned int i=0; i<output_data.boundary_forms.size(); ++i)
                {
                  output_data.JxW_values[i] = output_data.boundary_forms[i].norm() * weights[i];

                  if (subface_no != numbers::invalid_unsigned_int)
                    {
                      const double area_ratio = GeometryInfo<dim>::subface_ratio(cell->subface_case(face_no),
                                                                                 subface_no);
                      output_data.JxW_values[i] *= area_ratio;
                    }
                }

            if (update_flags & update_normal_vectors)
              for (unsigned int i=0; i<output_data.normal_vectors.size(); ++i)
                output_data.normal_vectors[i] = Point<spacedim>(output_data.boundary_forms[i] /
                                                                output_data.boundary_forms[i].norm());

            if (update_flags & update_jacobians)
              for (unsigned int point=0; point<n_q_points; ++point)
                output_data.jacobians[point] = data.contravariant[point];

            if (update_flags & update_inverse_jacobians)
              for (unsigned int point=0; point<n_q_points; ++point)
                output_data.inverse_jacobians[point] = data.covariant[point].transpose();
          }
      }


      template <int dim, int spacedim>
      void
      do_fill_fe_face_values (const ::MappingQGeneric<dim,spacedim>                             &mapping,
                              const typename ::Triangulation<dim,spacedim>::cell_iterator &cell,
                              const unsigned int                                                 face_no,
                              const unsigned int                                                 subface_no,
                              const typename QProjector<dim>::DataSetDescriptor                  data_set,
                              const Quadrature<dim-1>                                           &quadrature,
                              const typename ::MappingQGeneric<dim,spacedim>::InternalData      &data,
                              internal::FEValuesImplementation::MappingRelatedData<dim,spacedim>              &output_data)
      {
        if (dim>1 && data.tensor_product_quadrature)
          {
            maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>
            (CellSimilarity::none,
             data,
             output_data.quadrature_points,
             output_data.jacobian_grads);
          }
        else
          {
            maybe_compute_q_points<dim,spacedim> (data_set,
                                                  data,
                                                  output_data.quadrature_points);
            maybe_update_Jacobians<dim,spacedim> (CellSimilarity::none,
                                                  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,
                                 output_data);
      }
    }
  }
}



template <int dim, int spacedim>
void
MappingQGeneric<dim,spacedim>::
fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                     const unsigned int                                         face_no,
                     const Quadrature<dim-1>                                   &quadrature,
                     const typename Mapping<dim,spacedim>::InternalDataBase    &internal_data,
                     internal::FEValuesImplementation::MappingRelatedData<dim,spacedim>      &output_data) const
{
  // ensure that the following cast is really correct:
  Assert ((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
          ExcInternalError());
  const InternalData &data
    = static_cast<const InternalData &>(internal_data);

  // if necessary, recompute the support points of the transformation of this cell
  // (note that we need to first check the triangulation pointer, since otherwise
  // the second test might trigger an exception if the triangulations are not the
  // same)
  if ((data.mapping_support_points.size() == 0)
      ||
      (&cell->get_triangulation() !=
       &data.cell_of_current_support_points->get_triangulation())
      ||
      (cell != data.cell_of_current_support_points))
    {
      data.mapping_support_points = this->compute_mapping_support_points(cell);
      data.cell_of_current_support_points = cell;
    }

  internal::MappingQGenericImplementation::do_fill_fe_face_values
  (*this,
   cell, face_no, numbers::invalid_unsigned_int,
   QProjector<dim>::DataSetDescriptor::face (face_no,
                                             cell->face_orientation(face_no),
                                             cell->face_flip(face_no),
                                             cell->face_rotation(face_no),
                                             quadrature.size()),
   quadrature,
   data,
   output_data);
}



template <int dim, int spacedim>
void
MappingQGeneric<dim,spacedim>::
fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                        const unsigned int                                         face_no,
                        const unsigned int                                         subface_no,
                        const Quadrature<dim-1>                                   &quadrature,
                        const typename Mapping<dim,spacedim>::InternalDataBase    &internal_data,
                        internal::FEValuesImplementation::MappingRelatedData<dim,spacedim>      &output_data) const
{
  // ensure that the following cast is really correct:
  Assert ((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
          ExcInternalError());
  const InternalData &data
    = static_cast<const InternalData &>(internal_data);

  // if necessary, recompute the support points of the transformation of this cell
  // (note that we need to first check the triangulation pointer, since otherwise
  // the second test might trigger an exception if the triangulations are not the
  // same)
  if ((data.mapping_support_points.size() == 0)
      ||
      (&cell->get_triangulation() !=
       &data.cell_of_current_support_points->get_triangulation())
      ||
      (cell != data.cell_of_current_support_points))
    {
      data.mapping_support_points = this->compute_mapping_support_points(cell);
      data.cell_of_current_support_points = cell;
    }

  internal::MappingQGenericImplementation::do_fill_fe_face_values
  (*this,
   cell, face_no, subface_no,
   QProjector<dim>::DataSetDescriptor::subface (face_no, subface_no,
                                                cell->face_orientation(face_no),
                                                cell->face_flip(face_no),
                                                cell->face_rotation(face_no),
                                                quadrature.size(),
                                                cell->subface_case(face_no)),
   quadrature,
   data,
   output_data);
}



namespace internal
{
  namespace MappingQGenericImplementation
  {
    namespace
    {
      template <int dim, int spacedim, int rank>
      void
      transform_fields(const ArrayView<const Tensor<rank,dim> >               &input,
                       const MappingType                                       mapping_type,
                       const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
                       const ArrayView<Tensor<rank,spacedim> >                &output)
      {
        AssertDimension (input.size(), output.size());
        Assert ((dynamic_cast<const typename ::MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
                ExcInternalError());
        const typename ::MappingQGeneric<dim,spacedim>::InternalData
        &data = static_cast<const typename ::MappingQGeneric<dim,spacedim>::InternalData &>(mapping_data);

        switch (mapping_type)
          {
          case mapping_contravariant:
          {
            Assert (data.update_each & update_contravariant_transformation,
                    typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));

            for (unsigned int i=0; i<output.size(); ++i)
              output[i] = apply_transformation(data.contravariant[i], input[i]);

            return;
          }

          case mapping_piola:
          {
            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"));
            Assert (rank==1, ExcMessage("Only for rank 1"));
            if (rank!=1)
              return;

            for (unsigned int i=0; i<output.size(); ++i)
              {
                output[i] = apply_transformation(data.contravariant[i], input[i]);
                output[i] /= data.volume_elements[i];
              }
            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_covariant_transformation"));

            for (unsigned int i=0; i<output.size(); ++i)
              output[i] = apply_transformation(data.covariant[i], input[i]);

            return;
          }

          default:
            Assert(false, ExcNotImplemented());
          }
      }


      template <int dim, int spacedim, int rank>
      void
      transform_gradients(const ArrayView<const Tensor<rank,dim> >                &input,
                          const MappingType                                        mapping_type,
                          const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
                          const ArrayView<Tensor<rank,spacedim> >                 &output)
      {
        AssertDimension (input.size(), output.size());
        Assert ((dynamic_cast<const typename ::MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
                ExcInternalError());
        const typename ::MappingQGeneric<dim,spacedim>::InternalData
        &data = static_cast<const typename ::MappingQGeneric<dim,spacedim>::InternalData &>(mapping_data);

        switch (mapping_type)
          {
          case mapping_contravariant_gradient:
          {
            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)
              {
                DerivativeForm<1,spacedim,dim> A =
                  apply_transformation(data.contravariant[i], transpose(input[i]) );
                output[i] = apply_transformation(data.covariant[i], A.transpose() );
              }

            return;
          }

          case mapping_covariant_gradient:
          {
            Assert (data.update_each & update_covariant_transformation,
                    typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
            Assert (rank==2, ExcMessage("Only for rank 2"));

            for (unsigned int i=0; i<output.size(); ++i)
              {
                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,
                    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"));
            Assert (rank==2, ExcMessage("Only for rank 2"));

            for (unsigned int i=0; i<output.size(); ++i)
              {
                DerivativeForm<1,spacedim,dim> A =
                  apply_transformation(data.covariant[i], input[i] );
                Tensor<2,spacedim> T =
                  apply_transformation(data.contravariant[i], A.transpose() );

                output[i] = transpose(T);
                output[i] /= data.volume_elements[i];
              }

            return;
          }

          default:
            Assert(false, ExcNotImplemented());
          }
      }




      template <int dim, int spacedim>
      void
      transform_hessians(const ArrayView<const Tensor<3,dim> >                  &input,
                         const MappingType                                       mapping_type,
                         const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
                         const ArrayView<Tensor<3,spacedim> >                   &output)
      {
        AssertDimension (input.size(), output.size());
        Assert ((dynamic_cast<const typename ::MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
                ExcInternalError());
        const typename ::MappingQGeneric<dim,spacedim>::InternalData
        &data = static_cast<const typename ::MappingQGeneric<dim,spacedim>::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)
                {
                  double tmp1[dim][dim];
                  for (unsigned int J=0; J<dim; ++J)
                    for (unsigned int K=0; K<dim; ++K)
                      {
                        tmp1[J][K] = data.contravariant[q][i][0] * input[q][0][J][K];
                        for (unsigned int I=1; I<dim; ++I)
                          tmp1[J][K] += data.contravariant[q][i][I] * input[q][I][J][K];
                      }
                  for (unsigned int j=0; j<spacedim; ++j)
                    {
                      double tmp2[dim];
                      for (unsigned int K=0; K<dim; ++K)
                        {
                          tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
                          for (unsigned int J=1; J<dim; ++J)
                            tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
                        }
                      for (unsigned int k=0; k<spacedim; ++k)
                        {
                          output[q][i][j][k] = data.covariant[q][k][0] * tmp2[0];
                          for (unsigned int K=1; K<dim; ++K)
                            output[q][i][j][k] += data.covariant[q][k][K] * tmp2[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)
                {
                  double tmp1[dim][dim];
                  for (unsigned int J=0; J<dim; ++J)
                    for (unsigned int K=0; K<dim; ++K)
                      {
                        tmp1[J][K] = data.covariant[q][i][0] * input[q][0][J][K];
                        for (unsigned int I=1; I<dim; ++I)
                          tmp1[J][K] += data.covariant[q][i][I] * input[q][I][J][K];
                      }
                  for (unsigned int j=0; j<spacedim; ++j)
                    {
                      double tmp2[dim];
                      for (unsigned int K=0; K<dim; ++K)
                        {
                          tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
                          for (unsigned int J=1; J<dim; ++J)
                            tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
                        }
                      for (unsigned int k=0; k<spacedim; ++k)
                        {
                          output[q][i][j][k] = data.covariant[q][k][0] * tmp2[0];
                          for (unsigned int K=1; K<dim; ++K)
                            output[q][i][j][k] += data.covariant[q][k][K] * tmp2[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)
                {
                  double factor[dim];
                  for (unsigned int I=0; I<dim; ++I)
                    factor[I] = data.contravariant[q][i][I] / data.volume_elements[q];
                  double tmp1[dim][dim];
                  for (unsigned int J=0; J<dim; ++J)
                    for (unsigned int K=0; K<dim; ++K)
                      {
                        tmp1[J][K] = factor[0] * input[q][0][J][K];
                        for (unsigned int I=1; I<dim; ++I)
                          tmp1[J][K] += factor[I] * input[q][I][J][K];
                      }
                  for (unsigned int j=0; j<spacedim; ++j)
                    {
                      double tmp2[dim];
                      for (unsigned int K=0; K<dim; ++K)
                        {
                          tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
                          for (unsigned int J=1; J<dim; ++J)
                            tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
                        }
                      for (unsigned int k=0; k<spacedim; ++k)
                        {
                          output[q][i][j][k] = data.covariant[q][k][0] * tmp2[0];
                          for (unsigned int K=1; K<dim; ++K)
                            output[q][i][j][k] += data.covariant[q][k][K] * tmp2[K];
                        }
                    }
                }

            return;
          }

          default:
            Assert(false, ExcNotImplemented());
          }
      }




      template <int dim, int spacedim, int rank>
      void
      transform_differential_forms(const ArrayView<const DerivativeForm<rank, dim,spacedim> >   &input,
                                   const MappingType                                             mapping_type,
                                   const typename Mapping<dim,spacedim>::InternalDataBase       &mapping_data,
                                   const ArrayView<Tensor<rank+1, spacedim> >                   &output)
      {
        AssertDimension (input.size(), output.size());
        Assert ((dynamic_cast<const typename ::MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
                ExcInternalError());
        const typename ::MappingQGeneric<dim,spacedim>::InternalData
        &data = static_cast<const typename ::MappingQGeneric<dim,spacedim>::InternalData &>(mapping_data);

        switch (mapping_type)
          {
          case mapping_covariant:
          {
            Assert (data.update_each & 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]);

            return;
          }
          default:
            Assert(false, ExcNotImplemented());
          }
      }
    }
  }
}



template <int dim, int spacedim>
void
MappingQGeneric<dim,spacedim>::
transform (const ArrayView<const Tensor<1, dim> >                  &input,
           const MappingType                                        mapping_type,
           const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
           const ArrayView<Tensor<1, spacedim> >                   &output) const
{
  internal::MappingQGenericImplementation::transform_fields(input, mapping_type, mapping_data, output);
}



template <int dim, int spacedim>
void
MappingQGeneric<dim,spacedim>::
transform (const ArrayView<const DerivativeForm<1, dim,spacedim> >  &input,
           const MappingType                                         mapping_type,
           const typename Mapping<dim,spacedim>::InternalDataBase   &mapping_data,
           const ArrayView<Tensor<2, spacedim> >                    &output) const
{
  internal::MappingQGenericImplementation::transform_differential_forms(input, mapping_type, mapping_data, output);
}



template <int dim, int spacedim>
void
MappingQGeneric<dim,spacedim>::
transform (const ArrayView<const Tensor<2, dim> >                  &input,
           const MappingType                                        mapping_type,
           const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
           const ArrayView<Tensor<2, spacedim> >                   &output) const
{
  switch (mapping_type)
    {
    case mapping_contravariant:
      internal::MappingQGenericImplementation::transform_fields(input, mapping_type, mapping_data, output);
      return;

    case mapping_piola_gradient:
    case mapping_contravariant_gradient:
    case mapping_covariant_gradient:
      internal::MappingQGenericImplementation::transform_gradients(input, mapping_type, mapping_data, output);
      return;
    default:
      Assert(false, ExcNotImplemented());
    }
}



template <int dim, int spacedim>
void
MappingQGeneric<dim,spacedim>::
transform (const ArrayView<const  DerivativeForm<2, dim, spacedim> > &input,
           const MappingType                                          mapping_type,
           const typename Mapping<dim,spacedim>::InternalDataBase    &mapping_data,
           const ArrayView<Tensor<3,spacedim> >                      &output) const
{

  AssertDimension (input.size(), output.size());
  Assert (dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
          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)
            {
              double tmp[dim];
              for (unsigned int K=0; K<dim; ++K)
                {
                  tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
                  for (unsigned int J=1; J<dim; ++J)
                    tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
                }
              for (unsigned int k=0; k<spacedim; ++k)
                {
                  output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
                  for (unsigned int K=1; K<dim; ++K)
                    output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
                }
            }
      return;
    }

    default:
      Assert(false, ExcNotImplemented());
    }
}



template <int dim, int spacedim>
void
MappingQGeneric<dim,spacedim>::
transform (const ArrayView<const  Tensor<3,dim> >                  &input,
           const MappingType                                        mapping_type,
           const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
           const ArrayView<Tensor<3,spacedim> >                    &output) const
{
  switch (mapping_type)
    {
    case mapping_piola_hessian:
    case mapping_contravariant_hessian:
    case mapping_covariant_hessian:
      internal::MappingQGenericImplementation::transform_hessians(input, mapping_type, mapping_data, output);
      return;
    default:
      Assert(false, ExcNotImplemented());
    }
}



template <int dim, int spacedim>
void
MappingQGeneric<dim,spacedim>::
add_line_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                         std::vector<Point<spacedim> > &a) const
{
  // if we only need the midpoint, then ask for it.
  if (this->polynomial_degree==2)
    {
      for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
        {
          const typename Triangulation<dim,spacedim>::line_iterator line =
            (dim == 1  ?
             static_cast<typename Triangulation<dim,spacedim>::line_iterator>(cell) :
             cell->line(line_no));

          const Manifold<dim,spacedim> &manifold =
            ( ( line->manifold_id() == numbers::invalid_manifold_id ) &&
              ( dim < spacedim )
              ?
              cell->get_manifold()
              :
              line->get_manifold() );
          a.push_back(manifold.get_new_point_on_line(line));
        }
    }
  else
    // otherwise call the more complicated functions and ask for inner points
    // from the boundary description
    {
      std::vector<Point<spacedim> > tmp_points;
      // loop over each of the lines, and if it is at the boundary, then first
      // get the boundary description and second compute the points on it
      for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
        {
          const typename Triangulation<dim,spacedim>::line_iterator
          line = (dim == 1
                  ?
                  static_cast<typename Triangulation<dim,spacedim>::line_iterator>(cell)
                  :
                  cell->line(line_no));

          const Manifold<dim,spacedim> &manifold =
            ( ( line->manifold_id() == numbers::invalid_manifold_id ) &&
              ( dim < spacedim )
              ?
              cell->get_manifold() :
              line->get_manifold() );

          if (const Boundary<dim,spacedim> *boundary
              = dynamic_cast<const Boundary<dim,spacedim> *>(&manifold))
            {
              tmp_points.resize(this->polynomial_degree-1);
              boundary->get_intermediate_points_on_line(line, tmp_points);
              if (dim != 3 || cell->line_orientation(line_no))
                a.insert (a.end(), tmp_points.begin(), tmp_points.end());
              else
                a.insert (a.end(), tmp_points.rbegin(), tmp_points.rend());
            }
          else
            {
              const std::array<Point<spacedim>, 2> vertices
              {
                {
                  cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 0)),
                  cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 1))
                }
              };

              const std::size_t n_rows = support_point_weights_perimeter_to_interior[0].size(0);
              a.resize(a.size() + n_rows);
              auto a_view = make_array_view(a.end() - n_rows, a.end());
              manifold.get_new_points(make_array_view(vertices.begin(),
                                                      vertices.end()),
                                      support_point_weights_perimeter_to_interior[0],
                                      a_view);
            }
        }
    }
}



template <>
void
MappingQGeneric<3,3>::
add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell,
                        std::vector<Point<3> >                  &a) const
{
  const unsigned int faces_per_cell    = GeometryInfo<3>::faces_per_cell;

  // used if face quad at boundary or entirely in the interior of the domain
  std::vector<Point<3> > tmp_points;

  // loop over all faces and collect points on them
  for (unsigned int face_no=0; face_no<faces_per_cell; ++face_no)
    {
      const Triangulation<3>::face_iterator face = cell->face(face_no);

      // select the correct mappings for the present face
      const bool face_orientation = cell->face_orientation(face_no),
                 face_flip        = cell->face_flip       (face_no),
                 face_rotation    = cell->face_rotation   (face_no);

#ifdef DEBUG
      const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face,
                         lines_per_face    = GeometryInfo<3>::lines_per_face;

      // some sanity checks up front
      for (unsigned int i=0; i<vertices_per_face; ++i)
        Assert(face->vertex_index(i)==cell->vertex_index(
                 GeometryInfo<3>::face_to_cell_vertices(face_no, i,
                                                        face_orientation,
                                                        face_flip,
                                                        face_rotation)),
               ExcInternalError());

      // indices of the lines that bound a face are given by GeometryInfo<3>::
      // face_to_cell_lines
      for (unsigned int i=0; i<lines_per_face; ++i)
        Assert(face->line(i)==cell->line(GeometryInfo<3>::face_to_cell_lines(
                                           face_no, i, face_orientation, face_flip, face_rotation)),
               ExcInternalError());
#endif

      // On a quad, we have to check whether the manifold should determine the
      // point distribution from all surrounding points (new manifold code) or
      // the old-style Boundary code should simply return the intermediate
      // points. The second check is to find out whether the Boundary object
      // is actually a StraightBoundary (the default flat manifold assigned to
      // the triangulation if no manifold is assigned).
      const Boundary<3,3> *boundary =
        dynamic_cast<const Boundary<3,3> *>(&face->get_manifold());
      if (boundary != nullptr &&
          std::string(typeid(*boundary).name()).find("StraightBoundary") ==
          std::string::npos)
        {
          // ask the boundary/manifold object to return intermediate points on it
          tmp_points.resize((polynomial_degree-1)*(polynomial_degree-1));
          boundary->get_intermediate_points_on_quad(face, tmp_points);
          for (unsigned int i=0; i<tmp_points.size(); ++i)
            a.push_back(tmp_points[fe_q->adjust_quad_dof_index_for_face_orientation(i,
                                   face_orientation,
                                   face_flip,
                                   face_rotation)]);
        }
      else
        {
          // extract the points surrounding a quad from the points
          // already computed. First get the 4 vertices and then the points on
          // the four lines
          boost::container::small_vector<Point<3>, 200>
          tmp_points(GeometryInfo<2>::vertices_per_cell
                     + GeometryInfo<2>::lines_per_cell*(polynomial_degree-1));
          for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
            tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no,v)];
          if (polynomial_degree > 1)
            for (unsigned int line=0; line<GeometryInfo<2>::lines_per_cell; ++line)
              for (unsigned int i=0; i<polynomial_degree-1; ++i)
                tmp_points[4+line*(polynomial_degree-1)+i] =
                  a[GeometryInfo<3>::vertices_per_cell +
                    (polynomial_degree-1)*
                    GeometryInfo<3>::face_to_cell_lines(face_no,line) + i];

          const std::size_t n_rows = support_point_weights_perimeter_to_interior[1].size(0);
          a.resize(a.size() + n_rows);
          auto a_view = make_array_view(a.end() - n_rows, a.end());
          face->get_manifold().get_new_points (make_array_view(tmp_points.begin(),
                                                               tmp_points.end()),
                                               support_point_weights_perimeter_to_interior[1],
                                               a_view);
        }
    }
}



template <>
void
MappingQGeneric<2,3>::
add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell,
                        std::vector<Point<3> >                  &a) const
{
  if (const Boundary<2,3> *boundary =
        dynamic_cast<const Boundary<2,3> *>(&cell->get_manifold()))
    {
      std::vector<Point<3> > points((polynomial_degree-1)*(polynomial_degree-1));
      boundary->get_intermediate_points_on_quad(cell, points);
      a.insert(a.end(), points.begin(), points.end());
    }
  else
    {
      std::array<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
      for (unsigned int i=0; i<GeometryInfo<2>::vertices_per_cell; ++i)
        vertices[i] = cell->vertex(i);

      Table<2,double> weights(Utilities::fixed_power<2>(polynomial_degree-1),
                              GeometryInfo<2>::vertices_per_cell);
      for (unsigned int q=0, q2=0; q2<polynomial_degree-1; ++q2)
        for (unsigned int q1=0; q1<polynomial_degree-1; ++q1, ++q)
          {
            Point<2> point(line_support_points.point(q1+1)[0],
                           line_support_points.point(q2+1)[0]);
            for (unsigned int i=0; i<GeometryInfo<2>::vertices_per_cell; ++i)
              weights(q,i) = GeometryInfo<2>::d_linear_shape_function(point, i);
          }
      // TODO: use all surrounding points once Boundary path is removed
      const std::size_t n_rows = weights.size(0);
      a.resize(a.size() + n_rows);
      auto a_view = make_array_view(a.end() - n_rows, a.end());
      cell->get_manifold().get_new_points(make_array_view(vertices.begin(),
                                                          vertices.end()),
                                          weights,
                                          a_view);
    }
}



template <int dim, int spacedim>
void
MappingQGeneric<dim,spacedim>::
add_quad_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &,
                        std::vector<Point<spacedim> > &) const
{
  Assert (false, ExcInternalError());
}



template <int dim, int spacedim>
std::vector<Point<spacedim> >
MappingQGeneric<dim,spacedim>::
compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell) const
{
  // get the vertices first
  std::vector<Point<spacedim> > a;
  a.reserve(Utilities::fixed_power<dim>(polynomial_degree+1));
  for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
    a.push_back(cell->vertex(i));

  if (this->polynomial_degree > 1)
    {
      // check if all entities have the same manifold id which is when we can
      // simply ask the manifold for all points. the transfinite manifold can
      // do the interpolation better than this class, so if we detect that we
      // do not have to change anything here
      Assert(dim<=3, ExcImpossibleInDim(dim));
      bool all_manifold_ids_are_equal = (dim == spacedim);
      if (all_manifold_ids_are_equal &&
          dynamic_cast<const TransfiniteInterpolationManifold<dim,spacedim>*>(&cell->get_manifold()) == nullptr)
        {
          for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
            if (&cell->face(f)->get_manifold() != &cell->get_manifold())
              all_manifold_ids_are_equal = false;

          if (dim == 3)
            for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
              if (&cell->line(l)->get_manifold() != &cell->get_manifold())
                all_manifold_ids_are_equal = false;
        }

      if (all_manifold_ids_are_equal)
        {
          const std::size_t n_rows = support_point_weights_cell.size(0);
          a.resize(a.size() + n_rows);
          auto a_view = make_array_view(a.end() - n_rows, a.end());
          cell->get_manifold().get_new_points(make_array_view(a.begin(),
                                                              a.end() - n_rows),
                                              support_point_weights_cell,
                                              a_view);
        }
      else
        switch (dim)
          {
          case 1:
            add_line_support_points(cell, a);
            break;
          case 2:
            // in 2d, add the points on the four bounding lines to the exterior
            // (outer) points
            add_line_support_points(cell, a);

            // then get the interior support points
            if (dim != spacedim)
              add_quad_support_points(cell, a);
            else
              {
                const std::size_t n_rows = support_point_weights_perimeter_to_interior[1].size(0);
                a.resize(a.size() + n_rows);
                auto a_view = make_array_view(a.end() - n_rows, a.end());
                cell->get_manifold().get_new_points(make_array_view(a.begin(),
                                                                    a.end() - n_rows),
                                                    support_point_weights_perimeter_to_interior[1],
                                                    a_view);
              }
            break;

          case 3:
            // in 3d also add the points located on the boundary faces
            add_line_support_points (cell, a);
            add_quad_support_points (cell, a);

            // then compute the interior points
            {
              const std::size_t n_rows = support_point_weights_perimeter_to_interior[2].size(0);
              a.resize(a.size() + n_rows);
              auto a_view = make_array_view(a.end() - n_rows, a.end());
              cell->get_manifold().get_new_points(make_array_view(a.begin(),
                                                                  a.end() - n_rows),
                                                  support_point_weights_perimeter_to_interior[2],
                                                  a_view);
            }
            break;

          default:
            Assert(false, ExcNotImplemented());
            break;
          }
    }

  return a;
}



//--------------------------- Explicit instantiations -----------------------
#include "mapping_q_generic.inst"


DEAL_II_NAMESPACE_CLOSE