mapping_manifold.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 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/std_cxx14/memory.h>
#include <deal.II/base/tensor_product_polynomials.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/lac/full_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/dofs/dof_accessor.h>
#include <deal.II/grid/manifold.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_manifold.h>
#include <deal.II/fe/mapping_q1.h>

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


DEAL_II_NAMESPACE_OPEN

template <int dim, int spacedim>
std::size_t
MappingManifold<dim,spacedim>::InternalData::memory_consumption () const
{
  return (Mapping<dim,spacedim>::InternalDataBase::memory_consumption() +
          MemoryConsumption::memory_consumption (vertices) +
          MemoryConsumption::memory_consumption (cell) +
          MemoryConsumption::memory_consumption (quad) +
          MemoryConsumption::memory_consumption (cell_manifold_quadrature_weights) +
          MemoryConsumption::memory_consumption (vertex_weights) +
          MemoryConsumption::memory_consumption (unit_tangentials) +
          MemoryConsumption::memory_consumption (covariant) +
          MemoryConsumption::memory_consumption (contravariant) +
          MemoryConsumption::memory_consumption (aux) +
          MemoryConsumption::memory_consumption (volume_elements) +
          MemoryConsumption::memory_consumption (manifold) );
}


template <int dim, int spacedim>
void
MappingManifold<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;

  // Store the quadrature
  this->quad = q;

  // Resize the weights
  this->vertex_weights.resize(GeometryInfo<dim>::vertices_per_cell);

  // see if we need the (transformation) shape function values
  // and/or gradients and resize the necessary arrays
  if (this->update_each & (update_quadrature_points | update_contravariant_transformation) )
    compute_manifold_quadrature_weights(q);

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


template <int dim, int spacedim>
void
MappingManifold<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)
    {
      if (this->update_each & update_boundary_forms)
        {
          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 <int dim, int spacedim>
MappingManifold<dim,spacedim>::MappingManifold (const MappingManifold<dim,spacedim> &)
{}



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



template <int dim, int spacedim>
Point<dim>
MappingManifold<dim,spacedim>::
transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_iterator &,
                             const Point<spacedim> &) const
{
  Assert(false, ExcNotImplemented());
  return Point<dim>();
}



template <int dim, int spacedim>
Point<spacedim>
MappingManifold<dim,spacedim>::
transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                             const Point<dim> &p) const
{
  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
  std::array<double, GeometryInfo<dim>::vertices_per_cell> weights;

  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
    {
      vertices[v] = cell->vertex(v);
      weights[v] = GeometryInfo<dim>::d_linear_shape_function(p, v);
    }
  return cell->get_manifold().get_new_point(make_array_view(vertices.begin(),
                                                            vertices.end()),
                                            make_array_view(weights.begin(),
                                                            weights.end()));
}



// In the code below, GCC tries to instantiate MappingManifold<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 MappingManifold<3,4>
{};



template <int dim, int spacedim>
UpdateFlags
MappingManifold<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 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 elements wants to use the
      // contravariant of 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_JxW_values;

      if (out & update_normal_vectors)
        out |= update_JxW_values;
    }

  // Now throw an exception if we stumble upon something that was not
  // implemented yet
  Assert(!(out & (
             update_jacobian_grads |
             update_jacobian_pushed_forward_grads |
             update_jacobian_pushed_forward_2nd_derivatives |
             update_jacobian_pushed_forward_3rd_derivatives
           )), ExcNotImplemented());

  return out;
}



template <int dim, int spacedim>
std::unique_ptr<typename Mapping<dim,spacedim>::InternalDataBase>
MappingManifold<dim,spacedim>::get_data (const UpdateFlags update_flags,
                                         const Quadrature<dim> &q) const
{
  auto data = std_cxx14::make_unique<InternalData>();
  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>
MappingManifold<dim,spacedim>::get_face_data (const UpdateFlags        update_flags,
                                              const Quadrature<dim-1> &quadrature) const
{
  auto data = std_cxx14::make_unique<InternalData>();
  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>
MappingManifold<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
                                                 const Quadrature<dim-1>& quadrature) const
{
  auto data = std_cxx14::make_unique<InternalData>();
  data->initialize_face (this->requires_update_flags(update_flags),
                         QProjector<dim>::project_to_all_subfaces(quadrature),
                         quadrature.size());

  return std::move(data);
}



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

        AssertDimension(data.vertices.size(), GeometryInfo<dim>::vertices_per_cell);

        if (update_flags & update_quadrature_points)
          {
            for (unsigned int point=0; point<quadrature_points.size(); ++point)
              {
                quadrature_points[point] = data.manifold->get_new_point
                                           (make_array_view(data.vertices),
                                            make_array_view(data.cell_manifold_quadrature_weights[point+data_set]));
              }
          }
      }



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

        if (update_flags & update_contravariant_transformation)
          {
            const unsigned int n_q_points = data.contravariant.size();

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

            AssertDimension(GeometryInfo<dim>::vertices_per_cell,
                            data.vertices.size());
            for (unsigned int point=0; point<n_q_points; ++point)
              {
                // Start by figuring out how to compute the direction in
                // the reference space:
                const Point<dim> &p = data.quad.point(point+data_set);

                // And get its image on the manifold:
                const Point<spacedim> P = data.manifold->get_new_point
                                          (make_array_view(data.vertices),
                                           make_array_view(data.cell_manifold_quadrature_weights[point+data_set]));

                // To compute the Jacobian, we choose dim points aligned
                // with the dim reference axes, which are still in the
                // given cell, and ask for the tangent vector in these
                // directions. Choosing the points is somewhat arbitrary,
                // so we try to be smart and we pick points which are
                // on the opposite quadrant w.r.t. the evaluation
                // point.
                for (unsigned int i=0; i<dim; ++i)
                  {
                    const Point<dim> ei = Point<dim>::unit_vector(i);
                    const double pi = p[i];
                    Assert(pi >=0 && pi <= 1.0,
                           ExcInternalError("Was expecting a quadrature point "
                                            "inside the unit reference element."));

                    // In the length L, we store also the direction sign,
                    // which is positive, if the coordinate is < .5,
                    const double L = pi > .5 ? -pi: 1-pi;

                    const Point<dim> np(p + L*ei);

                    // Get the weights to compute the np point in real space
                    for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
                      data.vertex_weights[j] = GeometryInfo<dim>::d_linear_shape_function(np, j);

                    const Point<spacedim> NP=
                      data.manifold->get_new_point(make_array_view(data.vertices),
                                                   make_array_view(data.vertex_weights));

                    const Tensor<1,spacedim> T = data.manifold->get_tangent_vector(P, NP);

                    for (unsigned int d=0; d<spacedim; ++d)
                      data.contravariant[point][d][i] = T[d]/L;
                  }
              }

            if (update_flags & update_covariant_transformation)
              {
                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)
              {
                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>
CellSimilarity::Similarity
MappingManifold<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();

  data.store_vertices(cell);
  data.manifold = &(cell->get_manifold());

  internal::MappingManifoldImplementation::maybe_compute_q_points<dim,spacedim>
  (QProjector<dim>::DataSetDescriptor::cell (),
   data,
   output_data.quadrature_points);

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

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


      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]
                = std::sqrt(determinant(G)) * weights[point];

              if (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 (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 (cell_similarity != CellSimilarity::translation)
        for (unsigned int point=0; point<n_q_points; ++point)
          output_data.inverse_jacobians[point] = data.covariant[point].transpose();
    }

  return cell_similarity;
}






namespace internal
{
  namespace MappingManifoldImplementation
  {
    namespace
    {
      template <int dim, int spacedim>
      void
      maybe_compute_face_data
      (const ::MappingManifold<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 ::MappingManifold<dim,spacedim>::InternalData &data,
       internal::FEValuesImplementation::MappingRelatedData<dim,spacedim>               &output_data)
      {
        const UpdateFlags update_flags = data.update_each;

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

            // 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 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)
                  {
                    switch (dim)
                      {
                      case 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();

                        break;
                      }

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

                        break;
                      }

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

            if (update_flags & (update_normal_vectors
                                | update_JxW_values))
              for (unsigned int i=0; i<output_data.boundary_forms.size(); ++i)
                {
                  if (update_flags & update_JxW_values)
                    {
                      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)
                    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 ::MappingManifold<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 ::MappingManifold<dim,spacedim>::InternalData &data,
       internal::FEValuesImplementation::MappingRelatedData<dim,spacedim>               &output_data)
      {
        data.store_vertices(cell);

        data.manifold = &cell->face(face_no)->get_manifold();

        maybe_compute_q_points<dim,spacedim> (data_set,
                                              data,
                                              output_data.quadrature_points);
        maybe_update_Jacobians<dim,spacedim> (data_set,
                                              data);

        maybe_compute_face_data (mapping,
                                 cell, face_no, subface_no, quadrature.size(),
                                 quadrature.get_weights(), data,
                                 output_data);
      }

      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 ::MappingManifold<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
                ExcInternalError());
        const typename ::MappingManifold<dim,spacedim>::InternalData
        &data = static_cast<const typename ::MappingManifold<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 ::MappingManifold<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
                ExcInternalError());
        const typename ::MappingManifold<dim,spacedim>::InternalData
        &data = static_cast<const typename ::MappingManifold<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 ::MappingManifold<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
                ExcInternalError());
        const typename ::MappingManifold<dim,spacedim>::InternalData
        &data = static_cast<const typename ::MappingManifold<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 ::MappingManifold<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
                ExcInternalError());
        const typename ::MappingManifold<dim,spacedim>::InternalData
        &data = static_cast<const typename ::MappingManifold<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
MappingManifold<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);

  internal::MappingManifoldImplementation::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
MappingManifold<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);

  internal::MappingManifoldImplementation::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);
}



template <int dim, int spacedim>
void
MappingManifold<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::MappingManifoldImplementation::transform_fields(input, mapping_type, mapping_data, output);
}



template <int dim, int spacedim>
void
MappingManifold<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::MappingManifoldImplementation::transform_differential_forms(input, mapping_type, mapping_data, output);
}



template <int dim, int spacedim>
void
MappingManifold<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::MappingManifoldImplementation::transform_fields(input, mapping_type, mapping_data, output);
      return;

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



template <int dim, int spacedim>
void
MappingManifold<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
MappingManifold<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::MappingManifoldImplementation::transform_hessians(input, mapping_type, mapping_data, output);
      return;
    default:
      Assert(false, ExcNotImplemented());
    }
}

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


DEAL_II_NAMESPACE_CLOSE