]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Refactor MappingQ/Q1 into a common base class.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 4 Sep 2015 01:57:34 +0000 (20:57 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 6 Sep 2015 18:43:24 +0000 (13:43 -0500)
This patch separates out the common functionality for generic polynomial
mappings into a base class. It contains most of what used to be in MappingQ1
but really was not specific to p=1 mappings.

include/deal.II/fe/mapping_q.h
include/deal.II/fe/mapping_q1.h
include/deal.II/fe/mapping_q_generic.h [new file with mode: 0644]
source/fe/CMakeLists.txt
source/fe/mapping_q1.cc
source/fe/mapping_q_generic.cc [new file with mode: 0644]
source/fe/mapping_q_generic.inst.in [new file with mode: 0644]

index 470352e90d0514428a5ae95249d67b0232476c1b..7fb3cc036484fa6b8968a03736aade8700d43151 100644 (file)
@@ -20,6 +20,7 @@
 #include <deal.II/base/config.h>
 #include <deal.II/base/table.h>
 #include <deal.II/fe/mapping_q1.h>
+#include <deal.II/fe/mapping_q_generic.h>
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/grid/manifold.h>
 
@@ -56,8 +57,8 @@ class MappingQ : public MappingQ1<dim,spacedim>
 {
 public:
   /**
-   * Constructor.  @p p gives the degree of mapping polynomials on boundary
-   * cells.
+   * Constructor.  @p polynomial_degree denotes the polynomial degree
+   * of the polynomials that are used to map cells boundary.
    *
    * The second argument determines whether the higher order mapping should
    * also be used on interior cells. If its value is <code>false</code> (the
@@ -68,7 +69,7 @@ public:
    * would also like to use a higher order mapping in the interior. The
    * MappingQEulerian class is one such case.
    */
-  MappingQ (const unsigned int p,
+  MappingQ (const unsigned int polynomial_degree,
             const bool use_mapping_q_on_all_cells = false);
 
   /**
index bf07fa14b9ed5ad96dc3b269f8a7d3672378496b..bbdf4b4df8dcad725ae03200f4ca8c1fa87e293c 100644 (file)
@@ -23,7 +23,7 @@
 #include <deal.II/base/qprojector.h>
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/fe/mapping.h>
+#include <deal.II/fe/mapping_q_generic.h>
 
 #include <cmath>
 
@@ -37,7 +37,8 @@ DEAL_II_NAMESPACE_OPEN
  * Mapping of the reference to cell to a general
  * quadrilateral/hexahedra by $d$-linear shape functions.
  *
- * This function maps the reference (unit) cell to a general grid cell with
+ * This mapping implemented by this class maps the reference (unit) cell
+ * to a general grid cell with
  * straight lines in $d$ dimensions. (Note, however, that in 3D the
  * <i>faces</i> of a general, trilinearly mapped cell may be curved, even if the
  * edges are not). This is the standard mapping used for polyhedral domains. It
@@ -51,7 +52,7 @@ DEAL_II_NAMESPACE_OPEN
  * @author Guido Kanschat, 2000, 2001; Ralf Hartmann, 2000, 2001, 2005
  */
 template <int dim, int spacedim=dim>
-class MappingQ1 : public Mapping<dim,spacedim>
+class MappingQ1 : public MappingQGeneric<dim,spacedim>
 {
 public:
   /**
@@ -63,12 +64,6 @@ public:
   virtual
   Mapping<dim,spacedim> *clone () const;
 
-  /**
-   * Always returns @p true because MappingQ1 preserves vertex locations.
-   */
-  virtual
-  bool preserves_vertex_locations () const;
-
   /**
    * @name Mapping points between reference and real cells
    * @{
@@ -90,54 +85,6 @@ public:
    * @}
    */
 
-  /**
-   * @name Functions to transform tensors from reference to real coordinates
-   * @{
-   */
-
-  // for documentation, see the Mapping base class
-  virtual
-  void
-  transform (const VectorSlice<const std::vector<Tensor<1,dim> > >   input,
-             const MappingType                                       type,
-             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
-             VectorSlice<std::vector<Tensor<1,spacedim> > >          output) const;
-
-  // for documentation, see the Mapping base class
-  virtual
-  void
-  transform (const VectorSlice<const std::vector<DerivativeForm<1, dim,spacedim> > > input,
-             const MappingType                                                       type,
-             const typename Mapping<dim,spacedim>::InternalDataBase                 &internal,
-             VectorSlice<std::vector<Tensor<2,spacedim> > >                          output) const;
-
-  // for documentation, see the Mapping base class
-  virtual
-  void
-  transform (const VectorSlice<const std::vector<Tensor<2, dim> > >  input,
-             const MappingType                                       type,
-             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
-             VectorSlice<std::vector<Tensor<2,spacedim> > >          output) const;
-
-  // for documentation, see the Mapping base class
-  virtual
-  void
-  transform (const VectorSlice<const std::vector< DerivativeForm<2, dim, spacedim> > > input,
-             const MappingType                                                         type,
-             const typename Mapping<dim,spacedim>::InternalDataBase                   &internal,
-             VectorSlice<std::vector<Tensor<3,spacedim> > >                            output) const;
-
-  // for documentation, see the Mapping base class
-  virtual
-  void
-  transform (const VectorSlice<const std::vector<Tensor<3, dim> > >  input,
-             const MappingType                                       type,
-             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
-             VectorSlice<std::vector<Tensor<3,spacedim> > >          output) const;
-
-  /**
-   * @}
-   */
 
   /**
    * @name Interface with FEValues
@@ -146,298 +93,13 @@ public:
 
 public:
   /**
-   * Storage for internal data of d-linear mappings. See Mapping::InternalDataBase
-   * for an extensive description.
-   *
-   * For the current class MappingQ1, the InternalData class store
-   * data that is computed once when the object is created
-   * (in get_data()) as well as data the class wants to store from between
-   * the call to fill_fe_values(), fill_fe_face_values(), or
-   * fill_fe_subface_values() until possible later calls from the finite
-   * element to functions such as transform(). The latter class of
-   * member variables are marked as 'mutable'.
+   * Use the InternalData class of the base class without modification
+   * and additions.
    */
-  class InternalData : public Mapping<dim,spacedim>::InternalDataBase
-  {
-  public:
-    /**
-     * Constructor. The argument denotes the polynomial degree of
-     * the mapping to which this object will correspond.
-     */
-    InternalData(const unsigned int polynomial_degree);
-
-    /**
-     * Initialize the object's member variables related to cell data
-     * based on the given arguments.
-     *
-     * The function also calls compute_shape_function_values() to
-     * actually set the member variables related to the values and
-     * derivatives of the mapping shape functions.
-     */
-    void
-    initialize (const UpdateFlags      update_flags,
-                const Quadrature<dim> &quadrature,
-                const unsigned int     n_original_q_points);
-
-    /**
-     * Initialize the object's member variables related to cell and
-     * face data based on the given arguments. In order to initialize
-     * cell data, this function calls initialize().
-     */
-    void
-    initialize_face (const UpdateFlags      update_flags,
-                     const Quadrature<dim> &quadrature,
-                     const unsigned int     n_original_q_points);
-
-    /**
-     * Compute the values and/or derivatives of the shape functions
-     * used for the mapping.
-     *
-     * Which values, derivatives, or higher order derivatives are
-     * computed is determined by which of the member arrays have
-     * nonzero sizes. They are typically set to their appropriate
-     * sizes by the initialize() and initialize_face() functions,
-     * which indeed call this function internally. However, it is
-     * possible (and at times useful) to do the resizing by hand and
-     * then call this function directly. An example is in a Newton
-     * iteration where we update the location of a quadrature point
-     * (e.g., in MappingQ::transform_real_to_uni_cell()) and need to
-     * re-compute the mapping and its derivatives at this location,
-     * but have already sized all internal arrays correctly.
-     */
-    void compute_shape_function_values (const std::vector<Point<dim> > &unit_points);
-
-
-    /**
-     * Shape function at quadrature point. Shape functions are in tensor
-     * product order, so vertices must be reordered to obtain transformation.
-     */
-    const double &shape (const unsigned int qpoint,
-                         const unsigned int shape_nr) const;
-
-    /**
-     * Shape function at quadrature point. See above.
-     */
-    double &shape (const unsigned int qpoint,
-                   const unsigned int shape_nr);
-
-    /**
-     * Gradient of shape function in quadrature point. See above.
-     */
-    const Tensor<1,dim> &derivative (const unsigned int qpoint,
-                                     const unsigned int shape_nr) const;
-
-    /**
-     * Gradient of shape function in quadrature point. See above.
-     */
-    Tensor<1,dim> &derivative (const unsigned int qpoint,
-                               const unsigned int shape_nr);
-
-    /**
-     * Second derivative of shape function in quadrature point. See above.
-     */
-    const Tensor<2,dim> &second_derivative (const unsigned int qpoint,
-                                            const unsigned int shape_nr) const;
-
-    /**
-     * Second derivative of shape function in quadrature point. See above.
-     */
-    Tensor<2,dim> &second_derivative (const unsigned int qpoint,
-                                      const unsigned int shape_nr);
-
-    /**
-     * third derivative of shape function in quadrature point. See above.
-     */
-    const Tensor<3,dim> &third_derivative (const unsigned int qpoint,
-                                           const unsigned int shape_nr) const;
-
-    /**
-     * third derivative of shape function in quadrature point. See above.
-     */
-    Tensor<3,dim> &third_derivative (const unsigned int qpoint,
-                                     const unsigned int shape_nr);
-
-    /**
-     * fourth derivative of shape function in quadrature point. See above.
-     */
-    const Tensor<4,dim> &fourth_derivative (const unsigned int qpoint,
-                                            const unsigned int shape_nr) const;
-
-    /**
-     * fourth derivative of shape function in quadrature point. See above.
-     */
-    Tensor<4,dim> &fourth_derivative (const unsigned int qpoint,
-                                      const unsigned int shape_nr);
-
-    /**
-     * Return an estimate (in bytes) or the memory consumption of this object.
-     */
-    virtual std::size_t memory_consumption () const;
-
-    /**
-     * Values of shape functions. Access by function @p shape.
-     *
-     * Computed once.
-     */
-    std::vector<double> shape_values;
-
-    /**
-     * Values of shape function derivatives. Access by function @p derivative.
-     *
-     * Computed once.
-     */
-    std::vector<Tensor<1,dim> > shape_derivatives;
-
-    /**
-     * Values of shape function second derivatives. Access by function @p
-     * second_derivative.
-     *
-     * Computed once.
-     */
-    std::vector<Tensor<2,dim> > shape_second_derivatives;
-
-    /**
-     * Values of shape function third derivatives. Access by function @p
-     * second_derivative.
-     *
-     * Computed once.
-     */
-    std::vector<Tensor<3,dim> > shape_third_derivatives;
-
-    /**
-     * Values of shape function fourth derivatives. Access by function @p
-     * second_derivative.
-     *
-     * Computed once.
-     */
-    std::vector<Tensor<4,dim> > shape_fourth_derivatives;
-
-    /**
-     * Unit tangential vectors. Used for the computation of boundary forms and
-     * normal vectors.
-     *
-     * This vector has (dim-1)GeometryInfo::faces_per_cell entries. The first
-     * GeometryInfo::faces_per_cell contain the vectors in the first
-     * tangential direction for each face; the second set of
-     * GeometryInfo::faces_per_cell entries contain the vectors in the second
-     * tangential direction (only in 3d, since there we have 2 tangential
-     * directions per face), etc.
-     *
-     * Filled once.
-     */
-    std::vector<std::vector<Tensor<1,dim> > > unit_tangentials;
-
-    /**
-     * The polynomial degree of the mapping. Since the objects here
-     * are also used (with minor adjustments) by MappingQ, we need to
-     * store this.
-     */
-    unsigned int polynomial_degree;
-
-    /**
-     * Number of shape functions. If this is a Q1 mapping, then it is simply
-     * the number of vertices per cell. However, since also derived classes
-     * use this class (e.g. the Mapping_Q() class), the number of shape
-     * functions may also be different.
-     *
-     * In general, it is $(p+1)^\text{dim}$, where $p$ is the
-     * polynomial degree of the mapping.
-     */
-    const unsigned int n_shape_functions;
-
-    /**
-     * Tensors of covariant transformation at each of the quadrature points.
-     * The matrix stored is the Jacobian * G^{-1}, where G = Jacobian^{t} *
-     * Jacobian, is the first fundamental form of the map; if dim=spacedim
-     * then it reduces to the transpose of the inverse of the Jacobian matrix,
-     * which itself is stored in the @p contravariant field of this structure.
-     *
-     * Computed on each cell.
-     */
-    mutable std::vector<DerivativeForm<1,dim, spacedim > >  covariant;
-
-    /**
-     * Tensors of contravariant transformation at each of the quadrature
-     * points. The contravariant matrix is the Jacobian of the transformation,
-     * i.e. $J_{ij}=dx_i/d\hat x_j$.
-     *
-     * Computed on each cell.
-     */
-    mutable std::vector< DerivativeForm<1,dim,spacedim> > contravariant;
-
-    /**
-     * Auxiliary vectors for internal use.
-     */
-    mutable std::vector<std::vector<Tensor<1,spacedim> > > aux;
-
-    /**
-     * Stores the support points of the mapping shape functions on the @p
-     * cell_of_current_support_points.
-     */
-    mutable std::vector<Point<spacedim> > mapping_support_points;
-
-    /**
-     * Stores the cell of which the @p mapping_support_points are stored.
-     */
-    mutable typename Triangulation<dim,spacedim>::cell_iterator cell_of_current_support_points;
-
-    /**
-     * The determinant of the Jacobian in each quadrature point. Filled if
-     * #update_volume_elements.
-     */
-    mutable std::vector<double> volume_elements;
-  };
+  typedef typename MappingQGeneric<dim,spacedim>::InternalData InternalData;
 
 protected:
 
-  // documentation can be found in Mapping::requires_update_flags()
-  virtual
-  UpdateFlags
-  requires_update_flags (const UpdateFlags update_flags) const;
-
-  // documentation can be found in Mapping::get_data()
-  virtual
-  InternalData *
-  get_data (const UpdateFlags,
-            const Quadrature<dim> &quadrature) const;
-
-  // documentation can be found in Mapping::get_face_data()
-  virtual
-  typename Mapping<dim,spacedim>::InternalDataBase *
-  get_face_data (const UpdateFlags flags,
-                 const Quadrature<dim-1>& quadrature) const;
-
-  // documentation can be found in Mapping::get_subface_data()
-  virtual
-  typename Mapping<dim,spacedim>::InternalDataBase *
-  get_subface_data (const UpdateFlags flags,
-                    const Quadrature<dim-1>& quadrature) const;
-
-  // documentation can be found in Mapping::fill_fe_values()
-  virtual
-  CellSimilarity::Similarity
-  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,
-                  dealii::internal::FEValues::MappingRelatedData<dim, spacedim> &output_data) const;
-
-  // documentation can be found in Mapping::fill_fe_face_values()
-  virtual void
-  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,
-                       dealii::internal::FEValues::MappingRelatedData<dim, spacedim> &output_data) const;
-
-  // documentation can be found in Mapping::fill_fe_subface_values()
-  virtual void
-  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,
-                          dealii::internal::FEValues::MappingRelatedData<dim, spacedim> &output_data) const;
 
   /**
    * @}
@@ -480,15 +142,6 @@ protected:
                                              const Point<spacedim>                            &p) const;
 
 
-  /**
-   * Declare a convenience typedef for the class that describes offsets into
-   * quadrature formulas projected onto faces and subfaces.
-   */
-  typedef
-  typename QProjector<dim>::DataSetDescriptor
-  DataSetDescriptor;
-
-
   /**
    * Transforms a point @p p on the unit cell to the point @p p_real on the
    * real cell @p cell and returns @p p_real.
@@ -587,154 +240,6 @@ struct StaticMappingQ1
 
 /*@}*/
 
-/*----------------------------------------------------------------------*/
-
-#ifndef DOXYGEN
-
-template<int dim, int spacedim>
-inline
-const double &
-MappingQ1<dim,spacedim>::InternalData::shape (const unsigned int qpoint,
-                                              const unsigned int shape_nr) const
-{
-  Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(),
-         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
-                       shape_values.size()));
-  return shape_values [qpoint*n_shape_functions + shape_nr];
-}
-
-
-
-template<int dim, int spacedim>
-inline
-double &
-MappingQ1<dim,spacedim>::InternalData::shape (const unsigned int qpoint,
-                                              const unsigned int shape_nr)
-{
-  Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(),
-         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
-                       shape_values.size()));
-  return shape_values [qpoint*n_shape_functions + shape_nr];
-}
-
-
-template<int dim, int spacedim>
-inline
-const Tensor<1,dim> &
-MappingQ1<dim,spacedim>::InternalData::derivative (const unsigned int qpoint,
-                                                   const unsigned int shape_nr) const
-{
-  Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(),
-         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
-                       shape_derivatives.size()));
-  return shape_derivatives [qpoint*n_shape_functions + shape_nr];
-}
-
-
-
-template<int dim, int spacedim>
-inline
-Tensor<1,dim> &
-MappingQ1<dim,spacedim>::InternalData::derivative (const unsigned int qpoint,
-                                                   const unsigned int shape_nr)
-{
-  Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(),
-         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
-                       shape_derivatives.size()));
-  return shape_derivatives [qpoint*n_shape_functions + shape_nr];
-}
-
-
-template <int dim, int spacedim>
-inline
-const Tensor<2,dim> &
-MappingQ1<dim,spacedim>::InternalData::second_derivative (const unsigned int qpoint,
-                                                          const unsigned int shape_nr) const
-{
-  Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(),
-         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
-                       shape_second_derivatives.size()));
-  return shape_second_derivatives [qpoint*n_shape_functions + shape_nr];
-}
-
-
-template <int dim, int spacedim>
-inline
-Tensor<2,dim> &
-MappingQ1<dim,spacedim>::InternalData::second_derivative (const unsigned int qpoint,
-                                                          const unsigned int shape_nr)
-{
-  Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(),
-         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
-                       shape_second_derivatives.size()));
-  return shape_second_derivatives [qpoint*n_shape_functions + shape_nr];
-}
-
-template <int dim, int spacedim>
-inline
-const Tensor<3,dim> &
-MappingQ1<dim,spacedim>::InternalData::third_derivative (const unsigned int qpoint,
-                                                         const unsigned int shape_nr) const
-{
-  Assert(qpoint*n_shape_functions + shape_nr < shape_third_derivatives.size(),
-         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
-                       shape_third_derivatives.size()));
-  return shape_third_derivatives [qpoint*n_shape_functions + shape_nr];
-}
-
-
-template <int dim, int spacedim>
-inline
-Tensor<3,dim> &
-MappingQ1<dim,spacedim>::InternalData::third_derivative (const unsigned int qpoint,
-                                                         const unsigned int shape_nr)
-{
-  Assert(qpoint*n_shape_functions + shape_nr < shape_third_derivatives.size(),
-         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
-                       shape_third_derivatives.size()));
-  return shape_third_derivatives [qpoint*n_shape_functions + shape_nr];
-}
-
-
-template <int dim, int spacedim>
-inline
-const Tensor<4,dim> &
-MappingQ1<dim,spacedim>::InternalData::fourth_derivative (const unsigned int qpoint,
-                                                          const unsigned int shape_nr) const
-{
-  Assert(qpoint*n_shape_functions + shape_nr < shape_fourth_derivatives.size(),
-         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
-                       shape_fourth_derivatives.size()));
-  return shape_fourth_derivatives [qpoint*n_shape_functions + shape_nr];
-}
-
-
-template <int dim, int spacedim>
-inline
-Tensor<4,dim> &
-MappingQ1<dim,spacedim>::InternalData::fourth_derivative (const unsigned int qpoint,
-                                                          const unsigned int shape_nr)
-{
-  Assert(qpoint*n_shape_functions + shape_nr < shape_fourth_derivatives.size(),
-         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
-                       shape_fourth_derivatives.size()));
-  return shape_fourth_derivatives [qpoint*n_shape_functions + shape_nr];
-}
-
-
-
-template <int dim, int spacedim>
-inline
-bool
-MappingQ1<dim,spacedim>::preserves_vertex_locations () const
-{
-  return true;
-}
-
-#endif // DOXYGEN
-
-/* -------------- declaration of explicit specializations ------------- */
-
 
 DEAL_II_NAMESPACE_CLOSE
 
diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h
new file mode 100644 (file)
index 0000000..5a1e97b
--- /dev/null
@@ -0,0 +1,604 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2000 - 2015 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__mapping_q_generic_h
+#define dealii__mapping_q_generic_h
+
+
+#include <deal.II/base/derivative_form.h>
+#include <deal.II/base/config.h>
+#include <deal.II/base/table.h>
+#include <deal.II/base/qprojector.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/fe/mapping.h>
+
+#include <cmath>
+
+DEAL_II_NAMESPACE_OPEN
+
+/*!@addtogroup mapping */
+/*@{*/
+
+
+/**
+ * A base class for all polynomial mappings. In particular, this class
+ * provides the basis for the MappingQ1 and MappingQ classes that
+ * implement (bi-, tri-)linear mappings and higher order mappings,
+ * respectively.
+ *
+ *
+ * <h3>Implementation</h3>
+ *
+ * This class provides essentially the entire generic infrastructure
+ * for polynomial mappings. What it requires to work from derived
+ * classes is an implementation of the
+ * compute_mapping_support_points() class that provides a list of
+ * locations to which the support points of the mapping (e.g., the
+ * vertices of the cell in the lowest order case) should be mapped.
+ *
+ *
+ * @author Wolfgang Bangerth, 2015
+ */
+template <int dim, int spacedim=dim>
+class MappingQGeneric : public Mapping<dim,spacedim>
+{
+public:
+  /**
+   * Constructor.  @p polynomial_degree denotes the polynomial degree
+   * of the polynomials that are used to map cells from the reference
+   * to the real cell.
+   */
+  MappingQGeneric (const unsigned int polynomial_degree);
+
+  /**
+   * Always returns @p true because MappingQ1 preserves vertex locations.
+   */
+  virtual
+  bool preserves_vertex_locations () const;
+
+
+  /**
+   * @name Functions to transform tensors from reference to real coordinates
+   * @{
+   */
+
+  // for documentation, see the Mapping base class
+  virtual
+  void
+  transform (const VectorSlice<const std::vector<Tensor<1,dim> > >   input,
+             const MappingType                                       type,
+             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+             VectorSlice<std::vector<Tensor<1,spacedim> > >          output) const;
+
+  // for documentation, see the Mapping base class
+  virtual
+  void
+  transform (const VectorSlice<const std::vector<DerivativeForm<1, dim,spacedim> > > input,
+             const MappingType                                                       type,
+             const typename Mapping<dim,spacedim>::InternalDataBase                 &internal,
+             VectorSlice<std::vector<Tensor<2,spacedim> > >                          output) const;
+
+  // for documentation, see the Mapping base class
+  virtual
+  void
+  transform (const VectorSlice<const std::vector<Tensor<2, dim> > >  input,
+             const MappingType                                       type,
+             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+             VectorSlice<std::vector<Tensor<2,spacedim> > >          output) const;
+
+  // for documentation, see the Mapping base class
+  virtual
+  void
+  transform (const VectorSlice<const std::vector< DerivativeForm<2, dim, spacedim> > > input,
+             const MappingType                                                         type,
+             const typename Mapping<dim,spacedim>::InternalDataBase                   &internal,
+             VectorSlice<std::vector<Tensor<3,spacedim> > >                            output) const;
+
+  // for documentation, see the Mapping base class
+  virtual
+  void
+  transform (const VectorSlice<const std::vector<Tensor<3, dim> > >  input,
+             const MappingType                                       type,
+             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+             VectorSlice<std::vector<Tensor<3,spacedim> > >          output) const;
+
+  /**
+   * @}
+   */
+
+  /**
+   * @name Interface with FEValues
+   * @{
+   */
+
+public:
+  /**
+   * Storage for internal data of polynomial mappings. See
+   * Mapping::InternalDataBase for an extensive description.
+   *
+   * For the current class, the InternalData class stores
+   * data that is computed once when the object is created
+   * (in get_data()) as well as data the class wants to store from between
+   * the call to fill_fe_values(), fill_fe_face_values(), or
+   * fill_fe_subface_values() until possible later calls from the finite
+   * element to functions such as transform(). The latter class of
+   * member variables are marked as 'mutable'.
+   */
+  class InternalData : public Mapping<dim,spacedim>::InternalDataBase
+  {
+  public:
+    /**
+     * Constructor. The argument denotes the polynomial degree of
+     * the mapping to which this object will correspond.
+     */
+    InternalData(const unsigned int polynomial_degree);
+
+    /**
+     * Initialize the object's member variables related to cell data
+     * based on the given arguments.
+     *
+     * The function also calls compute_shape_function_values() to
+     * actually set the member variables related to the values and
+     * derivatives of the mapping shape functions.
+     */
+    void
+    initialize (const UpdateFlags      update_flags,
+                const Quadrature<dim> &quadrature,
+                const unsigned int     n_original_q_points);
+
+    /**
+     * Initialize the object's member variables related to cell and
+     * face data based on the given arguments. In order to initialize
+     * cell data, this function calls initialize().
+     */
+    void
+    initialize_face (const UpdateFlags      update_flags,
+                     const Quadrature<dim> &quadrature,
+                     const unsigned int     n_original_q_points);
+
+    /**
+     * Compute the values and/or derivatives of the shape functions
+     * used for the mapping.
+     *
+     * Which values, derivatives, or higher order derivatives are
+     * computed is determined by which of the member arrays have
+     * nonzero sizes. They are typically set to their appropriate
+     * sizes by the initialize() and initialize_face() functions,
+     * which indeed call this function internally. However, it is
+     * possible (and at times useful) to do the resizing by hand and
+     * then call this function directly. An example is in a Newton
+     * iteration where we update the location of a quadrature point
+     * (e.g., in MappingQ::transform_real_to_uni_cell()) and need to
+     * re-compute the mapping and its derivatives at this location,
+     * but have already sized all internal arrays correctly.
+     */
+    void compute_shape_function_values (const std::vector<Point<dim> > &unit_points);
+
+
+    /**
+     * Shape function at quadrature point. Shape functions are in tensor
+     * product order, so vertices must be reordered to obtain transformation.
+     */
+    const double &shape (const unsigned int qpoint,
+                         const unsigned int shape_nr) const;
+
+    /**
+     * Shape function at quadrature point. See above.
+     */
+    double &shape (const unsigned int qpoint,
+                   const unsigned int shape_nr);
+
+    /**
+     * Gradient of shape function in quadrature point. See above.
+     */
+    const Tensor<1,dim> &derivative (const unsigned int qpoint,
+                                     const unsigned int shape_nr) const;
+
+    /**
+     * Gradient of shape function in quadrature point. See above.
+     */
+    Tensor<1,dim> &derivative (const unsigned int qpoint,
+                               const unsigned int shape_nr);
+
+    /**
+     * Second derivative of shape function in quadrature point. See above.
+     */
+    const Tensor<2,dim> &second_derivative (const unsigned int qpoint,
+                                            const unsigned int shape_nr) const;
+
+    /**
+     * Second derivative of shape function in quadrature point. See above.
+     */
+    Tensor<2,dim> &second_derivative (const unsigned int qpoint,
+                                      const unsigned int shape_nr);
+
+    /**
+     * third derivative of shape function in quadrature point. See above.
+     */
+    const Tensor<3,dim> &third_derivative (const unsigned int qpoint,
+                                           const unsigned int shape_nr) const;
+
+    /**
+     * third derivative of shape function in quadrature point. See above.
+     */
+    Tensor<3,dim> &third_derivative (const unsigned int qpoint,
+                                     const unsigned int shape_nr);
+
+    /**
+     * fourth derivative of shape function in quadrature point. See above.
+     */
+    const Tensor<4,dim> &fourth_derivative (const unsigned int qpoint,
+                                            const unsigned int shape_nr) const;
+
+    /**
+     * fourth derivative of shape function in quadrature point. See above.
+     */
+    Tensor<4,dim> &fourth_derivative (const unsigned int qpoint,
+                                      const unsigned int shape_nr);
+
+    /**
+     * Return an estimate (in bytes) or the memory consumption of this object.
+     */
+    virtual std::size_t memory_consumption () const;
+
+    /**
+     * Values of shape functions. Access by function @p shape.
+     *
+     * Computed once.
+     */
+    std::vector<double> shape_values;
+
+    /**
+     * Values of shape function derivatives. Access by function @p derivative.
+     *
+     * Computed once.
+     */
+    std::vector<Tensor<1,dim> > shape_derivatives;
+
+    /**
+     * Values of shape function second derivatives. Access by function @p
+     * second_derivative.
+     *
+     * Computed once.
+     */
+    std::vector<Tensor<2,dim> > shape_second_derivatives;
+
+    /**
+     * Values of shape function third derivatives. Access by function @p
+     * second_derivative.
+     *
+     * Computed once.
+     */
+    std::vector<Tensor<3,dim> > shape_third_derivatives;
+
+    /**
+     * Values of shape function fourth derivatives. Access by function @p
+     * second_derivative.
+     *
+     * Computed once.
+     */
+    std::vector<Tensor<4,dim> > shape_fourth_derivatives;
+
+    /**
+     * Unit tangential vectors. Used for the computation of boundary forms and
+     * normal vectors.
+     *
+     * This vector has (dim-1)GeometryInfo::faces_per_cell entries. The first
+     * GeometryInfo::faces_per_cell contain the vectors in the first
+     * tangential direction for each face; the second set of
+     * GeometryInfo::faces_per_cell entries contain the vectors in the second
+     * tangential direction (only in 3d, since there we have 2 tangential
+     * directions per face), etc.
+     *
+     * Filled once.
+     */
+    std::vector<std::vector<Tensor<1,dim> > > unit_tangentials;
+
+    /**
+     * The polynomial degree of the mapping. Since the objects here
+     * are also used (with minor adjustments) by MappingQ, we need to
+     * store this.
+     */
+    unsigned int polynomial_degree;
+
+    /**
+     * Number of shape functions. If this is a Q1 mapping, then it is simply
+     * the number of vertices per cell. However, since also derived classes
+     * use this class (e.g. the Mapping_Q() class), the number of shape
+     * functions may also be different.
+     *
+     * In general, it is $(p+1)^\text{dim}$, where $p$ is the
+     * polynomial degree of the mapping.
+     */
+    const unsigned int n_shape_functions;
+
+    /**
+     * Tensors of covariant transformation at each of the quadrature points.
+     * The matrix stored is the Jacobian * G^{-1}, where G = Jacobian^{t} *
+     * Jacobian, is the first fundamental form of the map; if dim=spacedim
+     * then it reduces to the transpose of the inverse of the Jacobian matrix,
+     * which itself is stored in the @p contravariant field of this structure.
+     *
+     * Computed on each cell.
+     */
+    mutable std::vector<DerivativeForm<1,dim, spacedim > >  covariant;
+
+    /**
+     * Tensors of contravariant transformation at each of the quadrature
+     * points. The contravariant matrix is the Jacobian of the transformation,
+     * i.e. $J_{ij}=dx_i/d\hat x_j$.
+     *
+     * Computed on each cell.
+     */
+    mutable std::vector< DerivativeForm<1,dim,spacedim> > contravariant;
+
+    /**
+     * Auxiliary vectors for internal use.
+     */
+    mutable std::vector<std::vector<Tensor<1,spacedim> > > aux;
+
+    /**
+     * Stores the support points of the mapping shape functions on the @p
+     * cell_of_current_support_points.
+     */
+    mutable std::vector<Point<spacedim> > mapping_support_points;
+
+    /**
+     * Stores the cell of which the @p mapping_support_points are stored.
+     */
+    mutable typename Triangulation<dim,spacedim>::cell_iterator cell_of_current_support_points;
+
+    /**
+     * The determinant of the Jacobian in each quadrature point. Filled if
+     * #update_volume_elements.
+     */
+    mutable std::vector<double> volume_elements;
+  };
+
+protected:
+
+  // documentation can be found in Mapping::requires_update_flags()
+  virtual
+  UpdateFlags
+  requires_update_flags (const UpdateFlags update_flags) const;
+
+  // documentation can be found in Mapping::get_data()
+  virtual
+  InternalData *
+  get_data (const UpdateFlags,
+            const Quadrature<dim> &quadrature) const;
+
+  // documentation can be found in Mapping::get_face_data()
+  virtual
+  typename Mapping<dim,spacedim>::InternalDataBase *
+  get_face_data (const UpdateFlags flags,
+                 const Quadrature<dim-1>& quadrature) const;
+
+  // documentation can be found in Mapping::get_subface_data()
+  virtual
+  typename Mapping<dim,spacedim>::InternalDataBase *
+  get_subface_data (const UpdateFlags flags,
+                    const Quadrature<dim-1>& quadrature) const;
+
+  // documentation can be found in Mapping::fill_fe_values()
+  virtual
+  CellSimilarity::Similarity
+  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,
+                  dealii::internal::FEValues::MappingRelatedData<dim, spacedim> &output_data) const;
+
+  // documentation can be found in Mapping::fill_fe_face_values()
+  virtual void
+  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,
+                       dealii::internal::FEValues::MappingRelatedData<dim, spacedim> &output_data) const;
+
+  // documentation can be found in Mapping::fill_fe_subface_values()
+  virtual void
+  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,
+                          dealii::internal::FEValues::MappingRelatedData<dim, spacedim> &output_data) const;
+
+  /**
+   * @}
+   */
+
+
+  /**
+   * The degree of the polynomials used as shape functions for the mapping
+   * of cells.
+   */
+  const unsigned int polynomial_degree;
+
+  /**
+   * An interface that derived classes have to implement and that
+   * computes the locations of support points for the mapping. For
+   * example, for MappingQ1 these are the vertices. However, other
+   * classes may override this function differently. In particular,
+   * the MappingQ1Eulerian class does exactly this by not computing
+   * the support points from the geometry of the current cell but
+   * instead evaluating an externally given displacement field in
+   * addition to the geometry of the cell.
+   */
+  virtual
+  void
+  compute_mapping_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                                  std::vector<Point<spacedim> > &a) const = 0;
+};
+
+
+
+/*@}*/
+
+/*----------------------------------------------------------------------*/
+
+#ifndef DOXYGEN
+
+template<int dim, int spacedim>
+inline
+const double &
+MappingQGeneric<dim,spacedim>::InternalData::shape (const unsigned int qpoint,
+                                                    const unsigned int shape_nr) const
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_values.size()));
+  return shape_values [qpoint*n_shape_functions + shape_nr];
+}
+
+
+
+template<int dim, int spacedim>
+inline
+double &
+MappingQGeneric<dim,spacedim>::InternalData::shape (const unsigned int qpoint,
+                                                    const unsigned int shape_nr)
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_values.size()));
+  return shape_values [qpoint*n_shape_functions + shape_nr];
+}
+
+
+template<int dim, int spacedim>
+inline
+const Tensor<1,dim> &
+MappingQGeneric<dim,spacedim>::InternalData::derivative (const unsigned int qpoint,
+                                                         const unsigned int shape_nr) const
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_derivatives.size()));
+  return shape_derivatives [qpoint*n_shape_functions + shape_nr];
+}
+
+
+
+template<int dim, int spacedim>
+inline
+Tensor<1,dim> &
+MappingQGeneric<dim,spacedim>::InternalData::derivative (const unsigned int qpoint,
+                                                         const unsigned int shape_nr)
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_derivatives.size()));
+  return shape_derivatives [qpoint*n_shape_functions + shape_nr];
+}
+
+
+template <int dim, int spacedim>
+inline
+const Tensor<2,dim> &
+MappingQGeneric<dim,spacedim>::InternalData::second_derivative (const unsigned int qpoint,
+    const unsigned int shape_nr) const
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_second_derivatives.size()));
+  return shape_second_derivatives [qpoint*n_shape_functions + shape_nr];
+}
+
+
+template <int dim, int spacedim>
+inline
+Tensor<2,dim> &
+MappingQGeneric<dim,spacedim>::InternalData::second_derivative (const unsigned int qpoint,
+    const unsigned int shape_nr)
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_second_derivatives.size()));
+  return shape_second_derivatives [qpoint*n_shape_functions + shape_nr];
+}
+
+template <int dim, int spacedim>
+inline
+const Tensor<3,dim> &
+MappingQGeneric<dim,spacedim>::InternalData::third_derivative (const unsigned int qpoint,
+    const unsigned int shape_nr) const
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_third_derivatives.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_third_derivatives.size()));
+  return shape_third_derivatives [qpoint*n_shape_functions + shape_nr];
+}
+
+
+template <int dim, int spacedim>
+inline
+Tensor<3,dim> &
+MappingQGeneric<dim,spacedim>::InternalData::third_derivative (const unsigned int qpoint,
+    const unsigned int shape_nr)
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_third_derivatives.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_third_derivatives.size()));
+  return shape_third_derivatives [qpoint*n_shape_functions + shape_nr];
+}
+
+
+template <int dim, int spacedim>
+inline
+const Tensor<4,dim> &
+MappingQGeneric<dim,spacedim>::InternalData::fourth_derivative (const unsigned int qpoint,
+    const unsigned int shape_nr) const
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_fourth_derivatives.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_fourth_derivatives.size()));
+  return shape_fourth_derivatives [qpoint*n_shape_functions + shape_nr];
+}
+
+
+template <int dim, int spacedim>
+inline
+Tensor<4,dim> &
+MappingQGeneric<dim,spacedim>::InternalData::fourth_derivative (const unsigned int qpoint,
+    const unsigned int shape_nr)
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_fourth_derivatives.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_fourth_derivatives.size()));
+  return shape_fourth_derivatives [qpoint*n_shape_functions + shape_nr];
+}
+
+
+
+template <int dim, int spacedim>
+inline
+bool
+MappingQGeneric<dim,spacedim>::preserves_vertex_locations () const
+{
+  return true;
+}
+
+#endif // DOXYGEN
+
+/* -------------- declaration of explicit specializations ------------- */
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 228cd4e686c1170d200824500a49dc77dbdfe916..d00504a3c69a9d9194aba4e360235330e6d45af7 100644 (file)
@@ -52,6 +52,7 @@ SET(_src
   mapping_cartesian.cc
   mapping.cc
   mapping_fe_field.cc
+  mapping_q_generic.cc
   mapping_q1.cc
   mapping_q1_eulerian.cc
   mapping_q.cc
@@ -95,6 +96,7 @@ SET(_inst
   mapping_cartesian.inst.in
   mapping.inst.in
   mapping_fe_field.inst.in
+  mapping_q_generic.inst.in
   mapping_q1_eulerian.inst.in
   mapping_q1.inst.in
   mapping_q_eulerian.inst.in
index d4f83839d34ce0781ce72c18d772ae247cbf478d..2b522e6d0fa040c855e8b150b2e9c4338187e313 100644 (file)
 DEAL_II_NAMESPACE_OPEN
 
 
-template<int dim, int spacedim>
-MappingQ1<dim,spacedim>::InternalData::InternalData (const unsigned int polynomial_degree)
-  :
-  polynomial_degree (polynomial_degree),
-  n_shape_functions (Utilities::fixed_power<dim>(polynomial_degree+1))
-{}
-
-
-
-template<int dim, int spacedim>
-std::size_t
-MappingQ1<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
-MappingQ1<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);
-
-  // now also fill the various fields with their correct values
-  compute_shape_function_values (q.get_points());
-}
-
-
-
-template <int dim, int spacedim>
-void
-MappingQ1<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.
-          const unsigned int nfaces = GeometryInfo<dim>::faces_per_cell;
-          unit_tangentials.resize (nfaces*(dim-1),
-                                   std::vector<Tensor<1,dim> > (n_original_q_points));
-          if (dim==2)
-            {
-              // ensure a counterclockwise
-              // orientation of tangentials
-              static const int tangential_orientation[4]= {-1,1,1,-1};
-              for (unsigned int i=0; i<nfaces; ++i)
-                {
-                  Tensor<1,dim> tang;
-                  tang[1-i/2]=tangential_orientation[i];
-                  std::fill (unit_tangentials[i].begin(),
-                             unit_tangentials[i].end(), tang);
-                }
-            }
-          else if (dim==3)
-            {
-              for (unsigned int i=0; i<nfaces; ++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[nfaces+i].begin(),
-                             unit_tangentials[nfaces+i].end(), tang2);
-                }
-            }
-        }
-    }
-}
-
-
 
 
 template<int dim, int spacedim>
 MappingQ1<dim,spacedim>::MappingQ1 ()
+  :
+  MappingQGeneric<dim,spacedim> (1)
 {}
 
 
@@ -322,1030 +169,6 @@ namespace internal
       Assert(false, ExcInternalError());
       return Point<3>();
     }
-
-
-
-    template <int spacedim>
-    void
-    compute_shape_function_values (const unsigned int            n_shape_functions,
-                                   const std::vector<Point<1> > &unit_points,
-                                   typename dealii::MappingQ1<1,spacedim>::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)
-            {
-              // the following may or may not
-              // work if dim != spacedim
-              Assert (spacedim == 1, ExcNotImplemented());
-
-              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)
-            {
-              // if lower order derivative don't work, neither should this
-              Assert (spacedim == 1, ExcNotImplemented());
-
-              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)
-            {
-              // if lower order derivative don't work, neither should this
-              Assert (spacedim == 1, ExcNotImplemented());
-
-              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;
-            }
-        }
-    }
-
-
-    template <int spacedim>
-    void
-    compute_shape_function_values (const unsigned int            n_shape_functions,
-                                   const std::vector<Point<2> > &unit_points,
-                                   typename dealii::MappingQ1<2,spacedim>::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;
-            }
-        }
-    }
-
-
-
-    template <int spacedim>
-    void
-    compute_shape_function_values (const unsigned int            n_shape_functions,
-                                   const std::vector<Point<3> > &unit_points,
-                                   typename dealii::MappingQ1<3,spacedim>::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)
-            {
-              // the following may or may not
-              // work if dim != spacedim
-              Assert (spacedim == 3, ExcNotImplemented());
-
-              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)
-            {
-              // if lower order derivative don't work, neither should this
-              Assert (spacedim == 3, ExcNotImplemented());
-
-              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)
-            {
-              // if lower order derivative don't work, neither should this
-              Assert (spacedim == 3, ExcNotImplemented());
-
-              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;
-            }
-        }
-    }
-  }
-}
-
-
-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;
-  }
-}
-
-
-
-template<int dim, int spacedim>
-void
-MappingQ1<dim,spacedim>::InternalData::
-compute_shape_function_values (const std::vector<Point<dim> > &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)
-      &&
-      (dim == spacedim))
-    internal::MappingQ1::compute_shape_function_values<spacedim> (n_shape_functions,
-        unit_points, *this);
-  else
-    // otherwise ask an object that describes the polynomial space
-    {
-      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 QGaussLobatto<1> line_support_points (polynomial_degree + 1);
-      const TensorProductPolynomials<dim>
-      tensor_pols (Polynomials::generate_complete_Lagrange_basis(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> (get_dpo_vector<dim>(polynomial_degree), 1,
-                                          polynomial_degree)));
-
-      std::vector<double> values;
-      std::vector<Tensor<1,dim> > grads;
-      if (shape_values.size()!=0)
-        {
-          Assert(shape_values.size()==n_shape_functions*n_points,
-                 ExcInternalError());
-          values.resize(n_shape_functions);
-        }
-      if (shape_derivatives.size()!=0)
-        {
-          Assert(shape_derivatives.size()==n_shape_functions*n_points,
-                 ExcInternalError());
-          grads.resize(n_shape_functions);
-        }
-
-      std::vector<Tensor<2,dim> > grad2;
-      if (shape_second_derivatives.size()!=0)
-        {
-          Assert(shape_second_derivatives.size()==n_shape_functions*n_points,
-                 ExcInternalError());
-          grad2.resize(n_shape_functions);
-        }
-
-      std::vector<Tensor<3,dim> > grad3;
-      if (shape_third_derivatives.size()!=0)
-        {
-          Assert(shape_third_derivatives.size()==n_shape_functions*n_points,
-                 ExcInternalError());
-          grad3.resize(n_shape_functions);
-        }
-
-      std::vector<Tensor<4,dim> > grad4;
-      if (shape_fourth_derivatives.size()!=0)
-        {
-          Assert(shape_fourth_derivatives.size()==n_shape_functions*n_points,
-                 ExcInternalError());
-          grad4.resize(n_shape_functions);
-        }
-
-
-      if (shape_values.size()!=0 ||
-          shape_derivatives.size()!=0 ||
-          shape_second_derivatives.size()!=0 ||
-          shape_third_derivatives.size()!=0 ||
-          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 (shape_values.size()!=0)
-              for (unsigned int i=0; i<n_shape_functions; ++i)
-                shape(point,renumber[i]) = values[i];
-
-            if (shape_derivatives.size()!=0)
-              for (unsigned int i=0; i<n_shape_functions; ++i)
-                derivative(point,renumber[i]) = grads[i];
-
-            if (shape_second_derivatives.size()!=0)
-              for (unsigned int i=0; i<n_shape_functions; ++i)
-                second_derivative(point,renumber[i]) = grad2[i];
-
-            if (shape_third_derivatives.size()!=0)
-              for (unsigned int i=0; i<n_shape_functions; ++i)
-                third_derivative(point,renumber[i]) = grad3[i];
-
-            if (shape_fourth_derivatives.size()!=0)
-              for (unsigned int i=0; i<n_shape_functions; ++i)
-                fourth_derivative(point,renumber[i]) = grad4[i];
-          }
-    }
-}
-
-
-
-template<int dim, int spacedim>
-UpdateFlags
-MappingQ1<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;
-    }
-
-  return out;
-}
-
-
-
-template<int dim, int spacedim>
-typename MappingQ1<dim,spacedim>::InternalData *
-MappingQ1<dim,spacedim>::get_data (const UpdateFlags update_flags,
-                                   const Quadrature<dim> &q) const
-{
-  InternalData *data = new InternalData(1);
-  data->initialize (requires_update_flags(update_flags), q, q.size());
-
-  return data;
-}
-
-
-
-template<int dim, int spacedim>
-typename Mapping<dim,spacedim>::InternalDataBase *
-MappingQ1<dim,spacedim>::get_face_data (const UpdateFlags        update_flags,
-                                        const Quadrature<dim-1> &quadrature) const
-{
-  InternalData *data = new InternalData(1);
-  data->initialize_face (requires_update_flags(update_flags),
-                         QProjector<dim>::project_to_all_faces(quadrature),
-                         quadrature.size());
-
-  return data;
-}
-
-
-
-template<int dim, int spacedim>
-typename Mapping<dim,spacedim>::InternalDataBase *
-MappingQ1<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
-                                           const Quadrature<dim-1>& quadrature) const
-{
-  InternalData *data = new InternalData(1);
-  data->initialize_face (requires_update_flags(update_flags),
-                         QProjector<dim>::project_to_all_subfaces(quadrature),
-                         quadrature.size());
-
-  return data;
-}
-
-
-
-namespace internal
-{
-  namespace
-  {
-    /**
-     * Compute the locations of quadrature points on the object described by
-     * the first argument (and the cell for which the mapping support points
-     * have already been set), but only if the update_flags of the @p data
-     * argument indicate so.
-     */
-    template <int dim, int spacedim>
-    void
-    maybe_compute_q_points (const typename QProjector<dim>::DataSetDescriptor                 data_set,
-                            const typename dealii::MappingQ1<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;
-            }
-        }
-    }
-
-
-    /**
-     * Update the co- and contravariant matrices as well as their determinant, for the cell
-     * described stored in the data object, but only if the update_flags of the @p data
-     * argument indicate so.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim>
-    void
-    maybe_update_Jacobians (const CellSimilarity::Similarity                                   cell_similarity,
-                            const typename dealii::QProjector<dim>::DataSetDescriptor          data_set,
-                            const typename dealii::MappingQ1<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();
-          }
-
-    }
-
-    /**
-     * Update the Hessian of the transformation from unit to real cell, the
-     * Jacobian gradients.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim>
-    void
-    maybe_update_jacobian_grads (const CellSimilarity::Similarity                                   cell_similarity,
-                                 const typename QProjector<dim>::DataSetDescriptor                  data_set,
-                                 const typename dealii::MappingQ1<dim,spacedim>::InternalData      &data,
-                                 std::vector<DerivativeForm<2,dim,spacedim> >                      &jacobian_grads)
-    {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_grads)
-        {
-          const unsigned int n_q_points = jacobian_grads.size();
-
-          if (cell_similarity != CellSimilarity::translation)
-            {
-              for (unsigned int point=0; point<n_q_points; ++point)
-                {
-                  const Tensor<2,dim> *second =
-                    &data.second_derivative(point+data_set, 0);
-                  double result [spacedim][dim][dim];
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        result[i][j][l] = (second[0][j][l] *
-                                           data.mapping_support_points[0][i]);
-                  for (unsigned int k=1; k<data.n_shape_functions; ++k)
-                    for (unsigned int i=0; i<spacedim; ++i)
-                      for (unsigned int j=0; j<dim; ++j)
-                        for (unsigned int l=0; l<dim; ++l)
-                          result[i][j][l]
-                          += (second[k][j][l]
-                              *
-                              data.mapping_support_points[k][i]);
-
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        jacobian_grads[point][i][j][l] = result[i][j][l];
-                }
-            }
-        }
-    }
-
-    /**
-     * Update the Hessian of the transformation from unit to real cell, the
-     * Jacobian gradients, pushed forward to the real cell coordinates.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim>
-    void
-    maybe_update_jacobian_pushed_forward_grads (const CellSimilarity::Similarity                                   cell_similarity,
-                                                const typename QProjector<dim>::DataSetDescriptor                  data_set,
-                                                const typename dealii::MappingQ1<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)
-            {
-              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]);
-
-                  // pushing forward the derivative coordinates
-                  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] = result[i][0][0] *
-                                                                          data.covariant[point][j][0] *
-                                                                          data.covariant[point][l][0];
-                          for (unsigned int jr=0; jr<dim; ++jr)
-                            {
-                              const unsigned int lr_start = jr==0? 1:0;
-                              for (unsigned int lr=lr_start; lr<dim; ++lr)
-                                jacobian_pushed_forward_grads[point][i][j][l] += result[i][jr][lr] *
-                                                                                 data.covariant[point][j][jr] *
-                                                                                 data.covariant[point][l][lr];
-                            }
-                        }
-                }
-            }
-        }
-    }
-
-    /**
-         * Update the third derivatives of the transformation from unit to real cell, the
-         * Jacobian hessians.
-         *
-         * Skip the computation if possible as indicated by the first argument.
-         */
-    template <int dim, int spacedim>
-    void
-    maybe_update_jacobian_2nd_derivatives (const CellSimilarity::Similarity                              cell_similarity,
-                                           const typename QProjector<dim>::DataSetDescriptor             data_set,
-                                           const typename dealii::MappingQ1<dim,spacedim>::InternalData &data,
-                                           std::vector<DerivativeForm<3,dim,spacedim> >                 &jacobian_2nd_derivatives)
-    {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_2nd_derivatives)
-        {
-          const unsigned int n_q_points = jacobian_2nd_derivatives.size();
-
-          if (cell_similarity != CellSimilarity::translation)
-            {
-              for (unsigned int point=0; point<n_q_points; ++point)
-                {
-                  const Tensor<3,dim> *third =
-                    &data.third_derivative(point+data_set, 0);
-                  double result [spacedim][dim][dim][dim];
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          result[i][j][l][m] = (third[0][j][l][m] *
-                                                data.mapping_support_points[0][i]);
-                  for (unsigned int k=1; k<data.n_shape_functions; ++k)
-                    for (unsigned int i=0; i<spacedim; ++i)
-                      for (unsigned int j=0; j<dim; ++j)
-                        for (unsigned int l=0; l<dim; ++l)
-                          for (unsigned int m=0; m<dim; ++m)
-                            result[i][j][l][m]
-                            += (third[k][j][l][m]
-                                *
-                                data.mapping_support_points[k][i]);
-
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          jacobian_2nd_derivatives[point][i][j][l][m] = result[i][j][l][m];
-                }
-            }
-        }
-    }
-
-    /**
-     * Update the Hessian of the Hessian of the transformation from unit
-     * to real cell, the Jacobian Hessian gradients, pushed forward to the
-     * real cell coordinates.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim>
-    void
-    maybe_update_jacobian_pushed_forward_2nd_derivatives (const CellSimilarity::Similarity                                   cell_similarity,
-                                                          const typename QProjector<dim>::DataSetDescriptor                  data_set,
-                                                          const typename dealii::MappingQ1<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)
-            {
-              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<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]
-                              = result[i][0][0][0]*
-                                data.covariant[point][j][0] *
-                                data.covariant[point][l][0] *
-                                data.covariant[point][m][0];
-                            for (unsigned int jr=0; jr<dim; ++jr)
-                              for (unsigned int lr=0; lr<dim; ++lr)
-                                {
-                                  const unsigned int mr0 = (jr+lr == 0)? 1:0;
-                                  for (unsigned int mr=mr0; mr<dim; ++mr)
-                                    jacobian_pushed_forward_2nd_derivatives[point][i][j][l][m]
-                                    += result[i][jr][lr][mr] *
-                                       data.covariant[point][j][jr] *
-                                       data.covariant[point][l][lr] *
-                                       data.covariant[point][m][mr];
-                                }
-                          }
-                }
-            }
-        }
-    }
-
-    /**
-         * Update the fourth derivatives of the transformation from unit to real cell, the
-         * Jacobian hessian gradients.
-         *
-         * Skip the computation if possible as indicated by the first argument.
-         */
-    template <int dim, int spacedim>
-    void
-    maybe_update_jacobian_3rd_derivatives (const CellSimilarity::Similarity                              cell_similarity,
-                                           const typename QProjector<dim>::DataSetDescriptor             data_set,
-                                           const typename dealii::MappingQ1<dim,spacedim>::InternalData &data,
-                                           std::vector<DerivativeForm<4,dim,spacedim> >                 &jacobian_3rd_derivatives)
-    {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_3rd_derivatives)
-        {
-          const unsigned int n_q_points = jacobian_3rd_derivatives.size();
-
-          if (cell_similarity != CellSimilarity::translation)
-            {
-              for (unsigned int point=0; point<n_q_points; ++point)
-                {
-                  const Tensor<4,dim> *fourth =
-                    &data.fourth_derivative(point+data_set, 0);
-                  double result [spacedim][dim][dim][dim][dim];
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          for (unsigned int n=0; n<dim; ++n)
-                            result[i][j][l][m][n] = (fourth[0][j][l][m][n] *
-                                                     data.mapping_support_points[0][i]);
-                  for (unsigned int k=1; k<data.n_shape_functions; ++k)
-                    for (unsigned int i=0; i<spacedim; ++i)
-                      for (unsigned int j=0; j<dim; ++j)
-                        for (unsigned int l=0; l<dim; ++l)
-                          for (unsigned int m=0; m<dim; ++m)
-                            for (unsigned int n=0; n<dim; ++n)
-                              result[i][j][l][m][n]
-                              += (fourth[k][j][l][m][n]
-                                  *
-                                  data.mapping_support_points[k][i]);
-
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          for (unsigned int n=0; n<dim; ++n)
-                            jacobian_3rd_derivatives[point][i][j][l][m][n] = result[i][j][l][m][n];
-                }
-            }
-        }
-    }
-
-    /**
-     * Update the Hessian gradient of the transformation from unit to real cell, the
-     * Jacobian Hessians, pushed forward to the real cell coordinates.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim>
-    void
-    maybe_update_jacobian_pushed_forward_3rd_derivatives (const CellSimilarity::Similarity                                   cell_similarity,
-                                                          const typename QProjector<dim>::DataSetDescriptor                  data_set,
-                                                          const typename dealii::MappingQ1<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)
-            {
-              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<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]
-                                = result[i][0][0][0][0] *
-                                  data.covariant[point][j][0] *
-                                  data.covariant[point][l][0] *
-                                  data.covariant[point][m][0] *
-                                  data.covariant[point][n][0];
-                              for (unsigned int jr=0; jr<dim; ++jr)
-                                for (unsigned int lr=0; lr<dim; ++lr)
-                                  for (unsigned int mr=0; mr<dim; ++mr)
-                                    {
-                                      const unsigned int nr0 = (jr+lr+mr==0)? 1:0;
-                                      for (unsigned int nr=nr0; nr<dim; ++nr)
-                                        jacobian_pushed_forward_3rd_derivatives[point][i][j][l][m][n]
-                                        += result[i][jr][lr][mr][nr] *
-                                           data.covariant[point][j][jr] *
-                                           data.covariant[point][l][lr] *
-                                           data.covariant[point][m][mr] *
-                                           data.covariant[point][n][nr];
-                                    }
-                            }
-                }
-            }
-        }
-    }
   }
 }
 
@@ -1367,904 +190,6 @@ MappingQ1<dim,spacedim>::compute_mapping_support_points(
 
 
 
-template<int dim, int spacedim>
-CellSimilarity::Similarity
-MappingQ1<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::FEValues::MappingRelatedData<dim,spacedim>      &output_data) const
-{
-  // ensure that the following static_cast is really correct:
-  Assert (dynamic_cast<const InternalData *>(&internal_data) != 0,
-          ExcInternalError());
-  const InternalData &data = static_cast<const InternalData &>(internal_data);
-
-  const unsigned int n_q_points=quadrature.size();
-
-  // 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))
-    {
-      compute_mapping_support_points(cell, data.mapping_support_points);
-      data.cell_of_current_support_points = cell;
-    }
-
-  internal::maybe_compute_q_points<dim,spacedim> (DataSetDescriptor::cell (),
-                                                  data,
-                                                  output_data.quadrature_points);
-  internal::maybe_update_Jacobians<dim,spacedim> (cell_similarity,
-                                                  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));
-
-
-      if (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 (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
-                  {
-                    const unsigned int codim = spacedim-dim;
-                    (void)codim;
-
-                    if (update_flags & update_normal_vectors)
-                      {
-                        Assert( codim==1 , ExcMessage("There is no cell normal in codim 2."));
-
-                        if (dim==1)
-                          cross_product(output_data.normal_vectors[point],
-                                        -DX_t[0]);
-                        else //dim == 2
-                          cross_product(output_data.normal_vectors[point],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();
-    }
-
-  internal::maybe_update_jacobian_grads<dim,spacedim> (cell_similarity,
-                                                       DataSetDescriptor::cell (),
-                                                       data,
-                                                       output_data.jacobian_grads);
-
-  internal::maybe_update_jacobian_pushed_forward_grads<dim,spacedim> (cell_similarity,
-      DataSetDescriptor::cell (),
-      data,
-      output_data.jacobian_pushed_forward_grads);
-
-  internal::maybe_update_jacobian_2nd_derivatives<dim,spacedim> (cell_similarity,
-      DataSetDescriptor::cell (),
-      data,
-      output_data.jacobian_2nd_derivatives);
-
-  internal::maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,spacedim> (cell_similarity,
-      DataSetDescriptor::cell (),
-      data,
-      output_data.jacobian_pushed_forward_2nd_derivatives);
-
-  internal::maybe_update_jacobian_3rd_derivatives<dim,spacedim> (cell_similarity,
-      DataSetDescriptor::cell (),
-      data,
-      output_data.jacobian_3rd_derivatives);
-
-  internal::maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,spacedim> (cell_similarity,
-      DataSetDescriptor::cell (),
-      data,
-      output_data.jacobian_pushed_forward_3rd_derivatives);
-
-  return cell_similarity;
-}
-
-
-
-
-
-
-namespace internal
-{
-  namespace
-  {
-    /**
-     * Depending on what information is called for in the update flags of the
-     * @p data object, compute the various pieces of information that is required
-     * by the fill_fe_face_values() and fill_fe_subface_values() functions.
-     * This function simply unifies the work that would be done by
-     * those two functions.
-     *
-     * The resulting data is put into the @p output_data argument.
-     */
-    template <int dim, int spacedim>
-    void
-    maybe_compute_face_data (const dealii::MappingQ1<dim,spacedim> &mapping,
-                             const typename dealii::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 dealii::MappingQ1<dim,spacedim>::InternalData &data,
-                             internal::FEValues::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 (data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d],
-                                 mapping_contravariant,
-                                 data,
-                                 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:
-                    cross_product (output_data.boundary_forms[i], data.aux[0][i]);
-                    break;
-                  case 3:
-                    cross_product (output_data.boundary_forms[i], 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)
-                    {
-                      Tensor<1,spacedim> cell_normal;
-                      const DerivativeForm<1,spacedim,dim> DX_t =
-                        data.contravariant[point].transpose();
-                      cross_product(cell_normal,DX_t[0],DX_t[1]);
-                      cell_normal /= cell_normal.norm();
-
-                      // then compute the face normal from the face tangent
-                      // and the cell normal:
-                      cross_product (output_data.boundary_forms[point],
-                                     data.aux[0][point], cell_normal);
-                    }
-                }
-            }
-
-          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();
-        }
-    }
-
-
-    /**
-     * Do the work of MappingQ1::fill_fe_face_values() and
-     * MappingQ1::fill_fe_subface_values() in a generic way,
-     * using the 'data_set' to differentiate whether we will
-     * work on a face (and if so, which one) or subface.
-     */
-    template<int dim, int spacedim>
-    void
-    do_fill_fe_face_values (const dealii::MappingQ1<dim,spacedim>                             &mapping,
-                            const typename dealii::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 dealii::MappingQ1<dim,spacedim>::InternalData      &data,
-                            internal::FEValues::MappingRelatedData<dim,spacedim>              &output_data)
-    {
-      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
-MappingQ1<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::FEValues::MappingRelatedData<dim,spacedim>      &output_data) const
-{
-  // ensure that the following cast is really correct:
-  Assert ((dynamic_cast<const InternalData *>(&internal_data) != 0),
-          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))
-    {
-      compute_mapping_support_points(cell, data.mapping_support_points);
-      data.cell_of_current_support_points = cell;
-    }
-
-  internal::do_fill_fe_face_values (*this,
-                                    cell, face_no, numbers::invalid_unsigned_int,
-                                    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
-MappingQ1<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::FEValues::MappingRelatedData<dim,spacedim>      &output_data) const
-{
-  // ensure that the following cast is really correct:
-  Assert ((dynamic_cast<const InternalData *>(&internal_data) != 0),
-          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))
-    {
-      compute_mapping_support_points(cell, data.mapping_support_points);
-      data.cell_of_current_support_points = cell;
-    }
-
-  internal::do_fill_fe_face_values (*this,
-                                    cell, face_no, subface_no,
-                                    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
-{
-  template <int dim, int spacedim, int rank>
-  void
-  transform_fields(const VectorSlice<const std::vector<Tensor<rank,dim> > > input,
-                   const MappingType                                        mapping_type,
-                   const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
-                   VectorSlice<std::vector<Tensor<rank,spacedim> > >        output)
-  {
-    AssertDimension (input.size(), output.size());
-    Assert ((dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *>(&mapping_data) != 0),
-            ExcInternalError());
-    const typename MappingQ1<dim,spacedim>::InternalData
-    &data = static_cast<const typename MappingQ1<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 VectorSlice<const std::vector<Tensor<rank,dim> > > input,
-                      const MappingType                                        mapping_type,
-                      const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
-                      VectorSlice<std::vector<Tensor<rank,spacedim> > >        output)
-  {
-    AssertDimension (input.size(), output.size());
-    Assert ((dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *>(&mapping_data) != 0),
-            ExcInternalError());
-    const typename MappingQ1<dim,spacedim>::InternalData
-    &data = static_cast<const typename MappingQ1<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 VectorSlice<const std::vector<Tensor<3,dim> > >   input,
-                     const MappingType                                       mapping_type,
-                     const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
-                     VectorSlice<std::vector<Tensor<3,spacedim> > >          output)
-  {
-    AssertDimension (input.size(), output.size());
-    Assert ((dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *>(&mapping_data) != 0),
-            ExcInternalError());
-    const typename MappingQ1<dim,spacedim>::InternalData
-    &data = static_cast<const typename MappingQ1<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)
-            for (unsigned int j=0; j<spacedim; ++j)
-              for (unsigned int k=0; k<spacedim; ++k)
-                {
-                  output[q][i][j][k] =    data.contravariant[q][i][0]
-                                          * data.covariant[q][j][0]
-                                          * data.covariant[q][k][0]
-                                          * input[q][0][0][0];
-                  for (unsigned int I=0; I<dim; ++I)
-                    for (unsigned int J=0; J<dim; ++J)
-                      {
-                        const unsigned int K0 = (0==(I+J))? 1 : 0;
-                        for (unsigned int K=K0; K<dim; ++K)
-                          output[q][i][j][k] +=    data.contravariant[q][i][I]
-                                                   * data.covariant[q][j][J]
-                                                   * data.covariant[q][k][K]
-                                                   * input[q][I][J][K];
-                      }
-
-                }
-        return;
-      }
-
-      case mapping_covariant_hessian:
-      {
-        Assert (data.update_each & update_covariant_transformation,
-                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
-
-        for (unsigned int q=0; q<output.size(); ++q)
-          for (unsigned int i=0; i<spacedim; ++i)
-            for (unsigned int j=0; j<spacedim; ++j)
-              for (unsigned int k=0; k<spacedim; ++k)
-                {
-                  output[q][i][j][k] =    data.covariant[q][i][0]
-                                          * data.covariant[q][j][0]
-                                          * data.covariant[q][k][0]
-                                          * input[q][0][0][0];
-                  for (unsigned int I=0; I<dim; ++I)
-                    for (unsigned int J=0; J<dim; ++J)
-                      {
-                        const unsigned int K0 = (0==(I+J))? 1 : 0;
-                        for (unsigned int K=K0; K<dim; ++K)
-                          output[q][i][j][k] +=   data.covariant[q][i][I]
-                                                  * data.covariant[q][j][J]
-                                                  * data.covariant[q][k][K]
-                                                  * input[q][I][J][K];
-                      }
-
-                }
-
-        return;
-      }
-
-      case mapping_piola_hessian:
-      {
-        Assert (data.update_each & update_covariant_transformation,
-                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
-        Assert (data.update_each & update_contravariant_transformation,
-                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
-        Assert (data.update_each & update_volume_elements,
-                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_volume_elements"));
-
-        for (unsigned int q=0; q<output.size(); ++q)
-          for (unsigned int i=0; i<spacedim; ++i)
-            for (unsigned int j=0; j<spacedim; ++j)
-              for (unsigned int k=0; k<spacedim; ++k)
-                {
-                  output[q][i][j][k] =    data.contravariant[q][i][0]
-                                          / data.volume_elements[q]
-                                          * data.covariant[q][j][0]
-                                          * data.covariant[q][k][0]
-                                          * input[q][0][0][0];
-                  for (unsigned int I=0; I<dim; ++I)
-                    for (unsigned int J=0; J<dim; ++J)
-                      {
-                        const unsigned int K0 = (0==(I+J))? 1 : 0;
-                        for (unsigned int K=K0; K<dim; ++K)
-                          output[q][i][j][k] +=    data.contravariant[q][i][I]
-                                                   / data.volume_elements[q]
-                                                   * data.covariant[q][j][J]
-                                                   * data.covariant[q][k][K]
-                                                   * input[q][I][J][K];
-                      }
-
-                }
-
-        return;
-      }
-
-      default:
-        Assert(false, ExcNotImplemented());
-      }
-  }
-
-
-
-
-  template<int dim, int spacedim, int rank>
-  void
-  transform_differential_forms(const VectorSlice<const std::vector<DerivativeForm<rank, dim,spacedim> > > input,
-                               const MappingType                                                          mapping_type,
-                               const typename Mapping<dim,spacedim>::InternalDataBase                    &mapping_data,
-                               VectorSlice<std::vector<Tensor<rank+1, spacedim> > >                       output)
-  {
-    AssertDimension (input.size(), output.size());
-    Assert ((dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *>(&mapping_data) != 0),
-            ExcInternalError());
-    const typename MappingQ1<dim,spacedim>::InternalData
-    &data = static_cast<const typename MappingQ1<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
-MappingQ1<dim,spacedim>::
-transform (const VectorSlice<const std::vector<Tensor<1, dim> > >   input,
-           const MappingType                                        mapping_type,
-           const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
-           VectorSlice<std::vector<Tensor<1, spacedim> > >          output) const
-{
-  transform_fields(input, mapping_type, mapping_data, output);
-}
-
-
-
-template<int dim, int spacedim>
-void
-MappingQ1<dim,spacedim>::
-transform (const VectorSlice<const std::vector<DerivativeForm<1, dim,spacedim> > > input,
-           const MappingType                                                       mapping_type,
-           const typename Mapping<dim,spacedim>::InternalDataBase                 &mapping_data,
-           VectorSlice<std::vector<Tensor<2, spacedim> > >                         output) const
-{
-  transform_differential_forms(input, mapping_type, mapping_data, output);
-}
-
-
-
-template<int dim, int spacedim>
-void
-MappingQ1<dim,spacedim>::
-transform (const VectorSlice<const std::vector<Tensor<2, dim> > >   input,
-           const MappingType                                        mapping_type,
-           const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
-           VectorSlice<std::vector<Tensor<2, spacedim> > >          output) const
-{
-  switch (mapping_type)
-    {
-    case mapping_contravariant:
-      transform_fields(input, mapping_type, mapping_data, output);
-      return;
-
-    case mapping_piola_gradient:
-    case mapping_contravariant_gradient:
-    case mapping_covariant_gradient:
-      transform_gradients(input, mapping_type, mapping_data, output);
-      return;
-    default:
-      Assert(false, ExcNotImplemented());
-    }
-}
-
-
-
-template<int dim, int spacedim>
-void
-MappingQ1<dim,spacedim>::
-transform (const VectorSlice<const std::vector< DerivativeForm<2, dim, spacedim> > > input,
-           const MappingType                                                         mapping_type,
-           const typename Mapping<dim,spacedim>::InternalDataBase                   &mapping_data,
-           VectorSlice<std::vector<Tensor<3,spacedim> > >                            output) const
-{
-
-  AssertDimension (input.size(), output.size());
-  Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
-          ExcInternalError());
-  const InternalData &data = static_cast<const InternalData &>(mapping_data);
-
-  switch (mapping_type)
-    {
-    case mapping_covariant_gradient:
-    {
-      Assert (data.update_each & update_contravariant_transformation,
-              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
-
-      for (unsigned int q=0; q<output.size(); ++q)
-        for (unsigned int i=0; i<spacedim; ++i)
-          for (unsigned int j=0; j<spacedim; ++j)
-            for (unsigned int k=0; k<spacedim; ++k)
-              {
-                output[q][i][j][k] = data.covariant[q][j][0]
-                                     * data.covariant[q][k][0]
-                                     * input[q][i][0][0];
-                for (unsigned int J=0; J<dim; ++J)
-                  {
-                    const unsigned int K0 = (0==J)? 1 : 0;
-                    for (unsigned int K=K0; K<dim; ++K)
-                      output[q][i][j][k] += data.covariant[q][j][J]
-                                            * data.covariant[q][k][K]
-                                            * input[q][i][J][K];
-                  }
-
-              }
-      return;
-    }
-
-    default:
-      Assert(false, ExcNotImplemented());
-    }
-}
-
-
-
-template<int dim, int spacedim>
-void
-MappingQ1<dim,spacedim>::
-transform (const VectorSlice<const std::vector< Tensor<3,dim> > >   input,
-           const MappingType                                        mapping_type,
-           const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
-           VectorSlice<std::vector<Tensor<3,spacedim> > >           output) const
-{
-
-  switch (mapping_type)
-    {
-    case mapping_piola_hessian:
-    case mapping_contravariant_hessian:
-    case mapping_covariant_hessian:
-      transform_hessians(input, mapping_type, mapping_data, output);
-      return;
-    default:
-      Assert(false, ExcNotImplemented());
-    }
-
-}
-
-
 
 template<int dim, int spacedim>
 Point<spacedim>
diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc
new file mode 100644 (file)
index 0000000..3d8df89
--- /dev/null
@@ -0,0 +1,2259 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2000 - 2015 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/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/tensor_product_polynomials.h>
+#include <deal.II/base/memory_consumption.h>
+#include <deal.II/base/std_cxx11/array.h>
+#include <deal.II/base/std_cxx11/unique_ptr.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.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_q1.h>
+
+#include <cmath>
+#include <algorithm>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+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))
+{}
+
+
+
+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);
+
+  // now also fill the various fields with their correct values
+  compute_shape_function_values (q.get_points());
+}
+
+
+
+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)
+    {
+      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.
+          const unsigned int nfaces = GeometryInfo<dim>::faces_per_cell;
+          unit_tangentials.resize (nfaces*(dim-1),
+                                   std::vector<Tensor<1,dim> > (n_original_q_points));
+          if (dim==2)
+            {
+              // ensure a counterclockwise
+              // orientation of tangentials
+              static const int tangential_orientation[4]= {-1,1,1,-1};
+              for (unsigned int i=0; i<nfaces; ++i)
+                {
+                  Tensor<1,dim> tang;
+                  tang[1-i/2]=tangential_orientation[i];
+                  std::fill (unit_tangentials[i].begin(),
+                             unit_tangentials[i].end(), tang);
+                }
+            }
+          else if (dim==3)
+            {
+              for (unsigned int i=0; i<nfaces; ++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[nfaces+i].begin(),
+                             unit_tangentials[nfaces+i].end(), tang2);
+                }
+            }
+        }
+    }
+}
+
+
+
+namespace internal
+{
+  namespace MappingQGeneric
+  {
+    // 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_cxx11::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_cxx11::array<Point<spacedim>, GeometryInfo<2>::vertices_per_cell> &vertices,
+     const Point<spacedim> &p)
+    {
+      Assert(spacedim == 2, ExcInternalError());
+      const double x = p(0);
+      const double y = p(1);
+
+      const double x0 = vertices[0](0);
+      const double x1 = vertices[1](0);
+      const double x2 = vertices[2](0);
+      const double x3 = vertices[3](0);
+
+      const double y0 = vertices[0](1);
+      const double y1 = vertices[1](1);
+      const double y2 = vertices[2](1);
+      const double y3 = vertices[3](1);
+
+      const double a = (x1 - x3)*(y0 - y2) - (x0 - x2)*(y1 - y3);
+      const double b = -(x0 - x1 - x2 + x3)*y + (x - 2*x1 + x3)*y0 - (x - 2*x0 + x2)*y1
+                       - (x - x1)*y2 + (x - x0)*y3;
+      const double c = (x0 - x1)*y - (x - x1)*y0 + (x - x0)*y1;
+
+      const double discriminant = b*b - 4*a*c;
+      // fast exit if the point is not in the cell (this is the only case
+      // where the discriminant is negative)
+      if (discriminant < 0.0)
+        {
+          return Point<2>(2, 2);
+        }
+
+      double eta1;
+      double eta2;
+      // special case #1: if a is zero, then use the linear formula
+      if (a == 0.0 && b != 0.0)
+        {
+          eta1 = -c/b;
+          eta2 = -c/b;
+        }
+      // special case #2: if c is very small:
+      else if (std::abs(c/b) < 1e-12)
+        {
+          eta1 = (-b - std::sqrt(discriminant)) / (2*a);
+          eta2 = (-b + std::sqrt(discriminant)) / (2*a);
+        }
+      // finally, use the numerically stable version of the quadratic formula:
+      else
+        {
+          eta1 = 2*c / (-b - std::sqrt(discriminant));
+          eta2 = 2*c / (-b + std::sqrt(discriminant));
+        }
+      // pick the one closer to the center of the cell.
+      const 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 double subexpr0 = -eta*x2 + x0*(eta - 1);
+      const double xi_denominator0 = eta*x3 - x1*(eta - 1) + subexpr0;
+      const 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 double max_y = std::max(std::max(std::abs(y0), std::abs(y1)),
+                                        std::max(std::abs(y2), std::abs(y3)));
+          const double subexpr1 = -eta*y2 + y0*(eta - 1);
+          const 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
+            {
+              return Point<2>(2, 2);
+            }
+        }
+    }
+
+
+
+    template<int spacedim>
+    Point<3>
+    transform_real_to_unit_cell
+    (const std_cxx11::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 spacedim>
+    void
+    compute_shape_function_values (const unsigned int            n_shape_functions,
+                                   const std::vector<Point<1> > &unit_points,
+                                   typename dealii::MappingQGeneric<1,spacedim>::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)
+            {
+              // the following may or may not
+              // work if dim != spacedim
+              Assert (spacedim == 1, ExcNotImplemented());
+
+              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)
+            {
+              // if lower order derivative don't work, neither should this
+              Assert (spacedim == 1, ExcNotImplemented());
+
+              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)
+            {
+              // if lower order derivative don't work, neither should this
+              Assert (spacedim == 1, ExcNotImplemented());
+
+              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;
+            }
+        }
+    }
+
+
+    template <int spacedim>
+    void
+    compute_shape_function_values (const unsigned int            n_shape_functions,
+                                   const std::vector<Point<2> > &unit_points,
+                                   typename dealii::MappingQGeneric<2,spacedim>::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;
+            }
+        }
+    }
+
+
+
+    template <int spacedim>
+    void
+    compute_shape_function_values (const unsigned int            n_shape_functions,
+                                   const std::vector<Point<3> > &unit_points,
+                                   typename dealii::MappingQGeneric<3,spacedim>::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)
+            {
+              // the following may or may not
+              // work if dim != spacedim
+              Assert (spacedim == 3, ExcNotImplemented());
+
+              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)
+            {
+              // if lower order derivative don't work, neither should this
+              Assert (spacedim == 3, ExcNotImplemented());
+
+              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)
+            {
+              // if lower order derivative don't work, neither should this
+              Assert (spacedim == 3, ExcNotImplemented());
+
+              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;
+            }
+        }
+    }
+  }
+}
+
+
+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;
+  }
+}
+
+
+
+template<int dim, int spacedim>
+void
+MappingQGeneric<dim,spacedim>::InternalData::
+compute_shape_function_values (const std::vector<Point<dim> > &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)
+      &&
+      (dim == spacedim))
+    internal::MappingQGeneric::compute_shape_function_values<spacedim> (n_shape_functions,
+        unit_points, *this);
+  else
+    // otherwise ask an object that describes the polynomial space
+    {
+      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 QGaussLobatto<1> line_support_points (polynomial_degree + 1);
+      const TensorProductPolynomials<dim>
+      tensor_pols (Polynomials::generate_complete_Lagrange_basis(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> (get_dpo_vector<dim>(polynomial_degree), 1,
+                                          polynomial_degree)));
+
+      std::vector<double> values;
+      std::vector<Tensor<1,dim> > grads;
+      if (shape_values.size()!=0)
+        {
+          Assert(shape_values.size()==n_shape_functions*n_points,
+                 ExcInternalError());
+          values.resize(n_shape_functions);
+        }
+      if (shape_derivatives.size()!=0)
+        {
+          Assert(shape_derivatives.size()==n_shape_functions*n_points,
+                 ExcInternalError());
+          grads.resize(n_shape_functions);
+        }
+
+      std::vector<Tensor<2,dim> > grad2;
+      if (shape_second_derivatives.size()!=0)
+        {
+          Assert(shape_second_derivatives.size()==n_shape_functions*n_points,
+                 ExcInternalError());
+          grad2.resize(n_shape_functions);
+        }
+
+      std::vector<Tensor<3,dim> > grad3;
+      if (shape_third_derivatives.size()!=0)
+        {
+          Assert(shape_third_derivatives.size()==n_shape_functions*n_points,
+                 ExcInternalError());
+          grad3.resize(n_shape_functions);
+        }
+
+      std::vector<Tensor<4,dim> > grad4;
+      if (shape_fourth_derivatives.size()!=0)
+        {
+          Assert(shape_fourth_derivatives.size()==n_shape_functions*n_points,
+                 ExcInternalError());
+          grad4.resize(n_shape_functions);
+        }
+
+
+      if (shape_values.size()!=0 ||
+          shape_derivatives.size()!=0 ||
+          shape_second_derivatives.size()!=0 ||
+          shape_third_derivatives.size()!=0 ||
+          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 (shape_values.size()!=0)
+              for (unsigned int i=0; i<n_shape_functions; ++i)
+                shape(point,renumber[i]) = values[i];
+
+            if (shape_derivatives.size()!=0)
+              for (unsigned int i=0; i<n_shape_functions; ++i)
+                derivative(point,renumber[i]) = grads[i];
+
+            if (shape_second_derivatives.size()!=0)
+              for (unsigned int i=0; i<n_shape_functions; ++i)
+                second_derivative(point,renumber[i]) = grad2[i];
+
+            if (shape_third_derivatives.size()!=0)
+              for (unsigned int i=0; i<n_shape_functions; ++i)
+                third_derivative(point,renumber[i]) = grad3[i];
+
+            if (shape_fourth_derivatives.size()!=0)
+              for (unsigned int i=0; i<n_shape_functions; ++i)
+                fourth_derivative(point,renumber[i]) = grad4[i];
+          }
+    }
+}
+
+
+
+template<int dim, int spacedim>
+MappingQGeneric<dim,spacedim>::MappingQGeneric (const unsigned int p)
+  :
+  polynomial_degree(p)
+{}
+
+
+
+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
+      // 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;
+    }
+
+  return out;
+}
+
+
+
+template<int dim, int spacedim>
+typename MappingQGeneric<dim,spacedim>::InternalData *
+MappingQGeneric<dim,spacedim>::get_data (const UpdateFlags update_flags,
+                                         const Quadrature<dim> &q) const
+{
+  InternalData *data = new InternalData(polynomial_degree);
+  data->initialize (this->requires_update_flags(update_flags), q, q.size());
+
+  return data;
+}
+
+
+
+template<int dim, int spacedim>
+typename Mapping<dim,spacedim>::InternalDataBase *
+MappingQGeneric<dim,spacedim>::get_face_data (const UpdateFlags        update_flags,
+                                              const Quadrature<dim-1> &quadrature) const
+{
+  InternalData *data = new InternalData(polynomial_degree);
+  data->initialize_face (this->requires_update_flags(update_flags),
+                         QProjector<dim>::project_to_all_faces(quadrature),
+                         quadrature.size());
+
+  return data;
+}
+
+
+
+template<int dim, int spacedim>
+typename Mapping<dim,spacedim>::InternalDataBase *
+MappingQGeneric<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
+                                                 const Quadrature<dim-1>& quadrature) const
+{
+  InternalData *data = new InternalData(polynomial_degree);
+  data->initialize_face (this->requires_update_flags(update_flags),
+                         QProjector<dim>::project_to_all_subfaces(quadrature),
+                         quadrature.size());
+
+  return data;
+}
+
+
+
+namespace internal
+{
+  namespace
+  {
+    /**
+     * Compute the locations of quadrature points on the object described by
+     * the first argument (and the cell for which the mapping support points
+     * have already been set), but only if the update_flags of the @p data
+     * argument indicate so.
+     */
+    template <int dim, int spacedim>
+    void
+    maybe_compute_q_points (const typename QProjector<dim>::DataSetDescriptor                 data_set,
+                            const typename dealii::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;
+            }
+        }
+    }
+
+
+    /**
+     * Update the co- and contravariant matrices as well as their determinant, for the cell
+     * described stored in the data object, but only if the update_flags of the @p data
+     * argument indicate so.
+     *
+     * Skip the computation if possible as indicated by the first argument.
+     */
+    template <int dim, int spacedim>
+    void
+    maybe_update_Jacobians (const CellSimilarity::Similarity                                   cell_similarity,
+                            const typename dealii::QProjector<dim>::DataSetDescriptor          data_set,
+                            const typename dealii::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();
+          }
+
+    }
+
+    /**
+     * Update the Hessian of the transformation from unit to real cell, the
+     * Jacobian gradients.
+     *
+     * Skip the computation if possible as indicated by the first argument.
+     */
+    template <int dim, int spacedim>
+    void
+    maybe_update_jacobian_grads (const CellSimilarity::Similarity                                   cell_similarity,
+                                 const typename QProjector<dim>::DataSetDescriptor                  data_set,
+                                 const typename dealii::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];
+                }
+            }
+        }
+    }
+
+    /**
+     * Update the Hessian of the transformation from unit to real cell, the
+     * Jacobian gradients, pushed forward to the real cell coordinates.
+     *
+     * Skip the computation if possible as indicated by the first argument.
+     */
+    template <int dim, int spacedim>
+    void
+    maybe_update_jacobian_pushed_forward_grads (const CellSimilarity::Similarity                                   cell_similarity,
+                                                const typename QProjector<dim>::DataSetDescriptor                  data_set,
+                                                const typename dealii::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)
+            {
+              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]);
+
+                  // pushing forward the derivative coordinates
+                  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] = result[i][0][0] *
+                                                                          data.covariant[point][j][0] *
+                                                                          data.covariant[point][l][0];
+                          for (unsigned int jr=0; jr<dim; ++jr)
+                            {
+                              const unsigned int lr_start = jr==0? 1:0;
+                              for (unsigned int lr=lr_start; lr<dim; ++lr)
+                                jacobian_pushed_forward_grads[point][i][j][l] += result[i][jr][lr] *
+                                                                                 data.covariant[point][j][jr] *
+                                                                                 data.covariant[point][l][lr];
+                            }
+                        }
+                }
+            }
+        }
+    }
+
+    /**
+     * Update the third derivatives of the transformation from unit to real cell, the
+     * Jacobian hessians.
+     *
+     * Skip the computation if possible as indicated by the first argument.
+     */
+    template <int dim, int spacedim>
+    void
+    maybe_update_jacobian_2nd_derivatives (const CellSimilarity::Similarity                              cell_similarity,
+                                           const typename QProjector<dim>::DataSetDescriptor             data_set,
+                                           const typename dealii::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];
+                }
+            }
+        }
+    }
+
+    /**
+     * Update the Hessian of the Hessian of the transformation from unit
+     * to real cell, the Jacobian Hessian gradients, pushed forward to the
+     * real cell coordinates.
+     *
+     * Skip the computation if possible as indicated by the first argument.
+     */
+    template <int dim, int spacedim>
+    void
+    maybe_update_jacobian_pushed_forward_2nd_derivatives (const CellSimilarity::Similarity                                   cell_similarity,
+                                                          const typename QProjector<dim>::DataSetDescriptor                  data_set,
+                                                          const typename dealii::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)
+            {
+              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<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]
+                              = result[i][0][0][0]*
+                                data.covariant[point][j][0] *
+                                data.covariant[point][l][0] *
+                                data.covariant[point][m][0];
+                            for (unsigned int jr=0; jr<dim; ++jr)
+                              for (unsigned int lr=0; lr<dim; ++lr)
+                                {
+                                  const unsigned int mr_start = (jr+lr == 0)? 1:0;
+                                  for (unsigned int mr=0; mr<mr_start; ++mr)
+                                    jacobian_pushed_forward_2nd_derivatives[point][i][j][l][m]
+                                    += result[i][jr][lr][mr] *
+                                       data.covariant[point][j][jr] *
+                                       data.covariant[point][l][lr] *
+                                       data.covariant[point][m][mr];
+                                }
+                          }
+                }
+            }
+        }
+    }
+
+    /**
+         * Update the fourth derivatives of the transformation from unit to real cell, the
+         * Jacobian hessian gradients.
+         *
+         * Skip the computation if possible as indicated by the first argument.
+         */
+    template <int dim, int spacedim>
+    void
+    maybe_update_jacobian_3rd_derivatives (const CellSimilarity::Similarity                              cell_similarity,
+                                           const typename QProjector<dim>::DataSetDescriptor             data_set,
+                                           const typename dealii::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];
+                }
+            }
+        }
+    }
+
+    /**
+     * Update the Hessian gradient of the transformation from unit to real cell, the
+     * Jacobian Hessians, pushed forward to the real cell coordinates.
+     *
+     * Skip the computation if possible as indicated by the first argument.
+     */
+    template <int dim, int spacedim>
+    void
+    maybe_update_jacobian_pushed_forward_3rd_derivatives (const CellSimilarity::Similarity                                   cell_similarity,
+                                                          const typename QProjector<dim>::DataSetDescriptor                  data_set,
+                                                          const typename dealii::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)
+            {
+              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<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]
+                                = result[i][0][0][0][0] *
+                                  data.covariant[point][j][0] *
+                                  data.covariant[point][l][0] *
+                                  data.covariant[point][m][0] *
+                                  data.covariant[point][n][0];
+                              for (unsigned int jr=0; jr<dim; ++jr)
+                                for (unsigned int lr=0; lr<dim; ++lr)
+                                  for (unsigned int mr=0; mr<dim; ++mr)
+                                    {
+                                      const unsigned int nr_start = (jr+lr+mr==0)? 1:0;
+                                      for (unsigned int nr=0; nr<nr_start; ++nr)
+                                        jacobian_pushed_forward_3rd_derivatives[point][i][j][l][m][n]
+                                        += result[i][jr][lr][mr][nr] *
+                                           data.covariant[point][j][jr] *
+                                           data.covariant[point][l][lr] *
+                                           data.covariant[point][m][mr] *
+                                           data.covariant[point][n][nr];
+                                    }
+                            }
+                }
+            }
+        }
+    }
+  }
+}
+
+
+
+
+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::FEValues::MappingRelatedData<dim,spacedim>      &output_data) const
+{
+  // ensure that the following static_cast is really correct:
+  Assert (dynamic_cast<const InternalData *>(&internal_data) != 0,
+          ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(internal_data);
+
+  const unsigned int n_q_points=quadrature.size();
+
+  // 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))
+    {
+      compute_mapping_support_points(cell, data.mapping_support_points);
+      data.cell_of_current_support_points = cell;
+    }
+
+  internal::maybe_compute_q_points<dim,spacedim> (QProjector<dim>::DataSetDescriptor::cell (),
+                                                  data,
+                                                  output_data.quadrature_points);
+  internal::maybe_update_Jacobians<dim,spacedim> (cell_similarity,
+                                                  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));
+
+
+      if (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 (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
+                  {
+                    const unsigned int codim = spacedim-dim;
+                    (void)codim;
+
+                    if (update_flags & update_normal_vectors)
+                      {
+                        Assert( codim==1 , ExcMessage("There is no cell normal in codim 2."));
+
+                        if (dim==1)
+                          cross_product(output_data.normal_vectors[point],
+                                        -DX_t[0]);
+                        else //dim == 2
+                          cross_product(output_data.normal_vectors[point],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();
+    }
+
+  internal::maybe_update_jacobian_grads<dim,spacedim> (cell_similarity,
+                                                       QProjector<dim>::DataSetDescriptor::cell (),
+                                                       data,
+                                                       output_data.jacobian_grads);
+
+  internal::maybe_update_jacobian_pushed_forward_grads<dim,spacedim> (cell_similarity,
+      QProjector<dim>::DataSetDescriptor::cell (),
+      data,
+      output_data.jacobian_pushed_forward_grads);
+
+  internal::maybe_update_jacobian_2nd_derivatives<dim,spacedim> (cell_similarity,
+      QProjector<dim>::DataSetDescriptor::cell (),
+      data,
+      output_data.jacobian_2nd_derivatives);
+
+  internal::maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,spacedim> (cell_similarity,
+      QProjector<dim>::DataSetDescriptor::cell (),
+      data,
+      output_data.jacobian_pushed_forward_2nd_derivatives);
+
+  internal::maybe_update_jacobian_3rd_derivatives<dim,spacedim> (cell_similarity,
+      QProjector<dim>::DataSetDescriptor::cell (),
+      data,
+      output_data.jacobian_3rd_derivatives);
+
+  internal::maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,spacedim> (cell_similarity,
+      QProjector<dim>::DataSetDescriptor::cell (),
+      data,
+      output_data.jacobian_pushed_forward_3rd_derivatives);
+
+  return cell_similarity;
+}
+
+
+
+
+
+
+namespace internal
+{
+  namespace
+  {
+    /**
+     * Depending on what information is called for in the update flags of the
+     * @p data object, compute the various pieces of information that is required
+     * by the fill_fe_face_values() and fill_fe_subface_values() functions.
+     * This function simply unifies the work that would be done by
+     * those two functions.
+     *
+     * The resulting data is put into the @p output_data argument.
+     */
+    template <int dim, int spacedim>
+    void
+    maybe_compute_face_data (const dealii::MappingQGeneric<dim,spacedim> &mapping,
+                             const typename dealii::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 dealii::MappingQGeneric<dim,spacedim>::InternalData &data,
+                             internal::FEValues::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 (data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d],
+                                 mapping_contravariant,
+                                 data,
+                                 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:
+                    cross_product (output_data.boundary_forms[i], data.aux[0][i]);
+                    break;
+                  case 3:
+                    cross_product (output_data.boundary_forms[i], 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)
+                    {
+                      Tensor<1,spacedim> cell_normal;
+                      const DerivativeForm<1,spacedim,dim> DX_t =
+                        data.contravariant[point].transpose();
+                      cross_product(cell_normal,DX_t[0],DX_t[1]);
+                      cell_normal /= cell_normal.norm();
+
+                      // then compute the face normal from the face tangent
+                      // and the cell normal:
+                      cross_product (output_data.boundary_forms[point],
+                                     data.aux[0][point], cell_normal);
+                    }
+                }
+            }
+
+          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();
+        }
+    }
+
+
+    /**
+     * Do the work of MappingQGeneric::fill_fe_face_values() and
+     * MappingQGeneric::fill_fe_subface_values() in a generic way,
+     * using the 'data_set' to differentiate whether we will
+     * work on a face (and if so, which one) or subface.
+     */
+    template<int dim, int spacedim>
+    void
+    do_fill_fe_face_values (const dealii::MappingQGeneric<dim,spacedim>                             &mapping,
+                            const typename dealii::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 dealii::MappingQGeneric<dim,spacedim>::InternalData      &data,
+                            internal::FEValues::MappingRelatedData<dim,spacedim>              &output_data)
+    {
+      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::FEValues::MappingRelatedData<dim,spacedim>      &output_data) const
+{
+  // ensure that the following cast is really correct:
+  Assert ((dynamic_cast<const InternalData *>(&internal_data) != 0),
+          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))
+    {
+      compute_mapping_support_points(cell, data.mapping_support_points);
+      data.cell_of_current_support_points = cell;
+    }
+
+  internal::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::FEValues::MappingRelatedData<dim,spacedim>      &output_data) const
+{
+  // ensure that the following cast is really correct:
+  Assert ((dynamic_cast<const InternalData *>(&internal_data) != 0),
+          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))
+    {
+      compute_mapping_support_points(cell, data.mapping_support_points);
+      data.cell_of_current_support_points = cell;
+    }
+
+  internal::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
+{
+  template <int dim, int spacedim, int rank>
+  void
+  transform_fields(const VectorSlice<const std::vector<Tensor<rank,dim> > > input,
+                   const MappingType                                        mapping_type,
+                   const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
+                   VectorSlice<std::vector<Tensor<rank,spacedim> > >        output)
+  {
+    AssertDimension (input.size(), output.size());
+    Assert ((dynamic_cast<const typename MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != 0),
+            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 VectorSlice<const std::vector<Tensor<rank,dim> > > input,
+                      const MappingType                                        mapping_type,
+                      const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
+                      VectorSlice<std::vector<Tensor<rank,spacedim> > >        output)
+  {
+    AssertDimension (input.size(), output.size());
+    Assert ((dynamic_cast<const typename MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != 0),
+            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 VectorSlice<const std::vector<Tensor<3,dim> > >   input,
+                     const MappingType                                       mapping_type,
+                     const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+                     VectorSlice<std::vector<Tensor<3,spacedim> > >          output)
+  {
+    AssertDimension (input.size(), output.size());
+    Assert ((dynamic_cast<const typename MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != 0),
+            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)
+            for (unsigned int j=0; j<spacedim; ++j)
+              for (unsigned int k=0; k<spacedim; ++k)
+                {
+                  output[q][i][j][k] =    data.contravariant[q][i][0]
+                                          * data.covariant[q][j][0]
+                                          * data.covariant[q][k][0]
+                                          * input[q][0][0][0];
+                  for (unsigned int I=0; I<dim; ++I)
+                    for (unsigned int J=0; J<dim; ++J)
+                      {
+                        const unsigned int K0 = (0==(I+J))? 1 : 0;
+                        for (unsigned int K=K0; K<dim; ++K)
+                          output[q][i][j][k] +=    data.contravariant[q][i][I]
+                                                   * data.covariant[q][j][J]
+                                                   * data.covariant[q][k][K]
+                                                   * input[q][I][J][K];
+                      }
+
+                }
+        return;
+      }
+
+      case mapping_covariant_hessian:
+      {
+        Assert (data.update_each & update_covariant_transformation,
+                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+
+        for (unsigned int q=0; q<output.size(); ++q)
+          for (unsigned int i=0; i<spacedim; ++i)
+            for (unsigned int j=0; j<spacedim; ++j)
+              for (unsigned int k=0; k<spacedim; ++k)
+                {
+                  output[q][i][j][k] =    data.covariant[q][i][0]
+                                          * data.covariant[q][j][0]
+                                          * data.covariant[q][k][0]
+                                          * input[q][0][0][0];
+                  for (unsigned int I=0; I<dim; ++I)
+                    for (unsigned int J=0; J<dim; ++J)
+                      {
+                        const unsigned int K0 = (0==(I+J))? 1 : 0;
+                        for (unsigned int K=K0; K<dim; ++K)
+                          output[q][i][j][k] +=   data.covariant[q][i][I]
+                                                  * data.covariant[q][j][J]
+                                                  * data.covariant[q][k][K]
+                                                  * input[q][I][J][K];
+                      }
+
+                }
+
+        return;
+      }
+
+      case mapping_piola_hessian:
+      {
+        Assert (data.update_each & update_covariant_transformation,
+                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+        Assert (data.update_each & update_contravariant_transformation,
+                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+        Assert (data.update_each & update_volume_elements,
+                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_volume_elements"));
+
+        for (unsigned int q=0; q<output.size(); ++q)
+          for (unsigned int i=0; i<spacedim; ++i)
+            for (unsigned int j=0; j<spacedim; ++j)
+              for (unsigned int k=0; k<spacedim; ++k)
+                {
+                  output[q][i][j][k] =    data.contravariant[q][i][0]
+                                          / data.volume_elements[q]
+                                          * data.covariant[q][j][0]
+                                          * data.covariant[q][k][0]
+                                          * input[q][0][0][0];
+                  for (unsigned int I=0; I<dim; ++I)
+                    for (unsigned int J=0; J<dim; ++J)
+                      {
+                        const unsigned int K0 = (0==(I+J))? 1 : 0;
+                        for (unsigned int K=K0; K<dim; ++K)
+                          output[q][i][j][k] +=    data.contravariant[q][i][I]
+                                                   / data.volume_elements[q]
+                                                   * data.covariant[q][j][J]
+                                                   * data.covariant[q][k][K]
+                                                   * input[q][I][J][K];
+                      }
+
+                }
+
+        return;
+      }
+
+      default:
+        Assert(false, ExcNotImplemented());
+      }
+  }
+
+
+
+
+  template<int dim, int spacedim, int rank>
+  void
+  transform_differential_forms(const VectorSlice<const std::vector<DerivativeForm<rank, dim,spacedim> > > input,
+                               const MappingType                                                          mapping_type,
+                               const typename Mapping<dim,spacedim>::InternalDataBase                    &mapping_data,
+                               VectorSlice<std::vector<Tensor<rank+1, spacedim> > >                       output)
+  {
+    AssertDimension (input.size(), output.size());
+    Assert ((dynamic_cast<const typename MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != 0),
+            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 VectorSlice<const std::vector<Tensor<1, dim> > >   input,
+           const MappingType                                        mapping_type,
+           const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
+           VectorSlice<std::vector<Tensor<1, spacedim> > >          output) const
+{
+  transform_fields(input, mapping_type, mapping_data, output);
+}
+
+
+
+template<int dim, int spacedim>
+void
+MappingQGeneric<dim,spacedim>::
+transform (const VectorSlice<const std::vector<DerivativeForm<1, dim,spacedim> > > input,
+           const MappingType                                                       mapping_type,
+           const typename Mapping<dim,spacedim>::InternalDataBase                 &mapping_data,
+           VectorSlice<std::vector<Tensor<2, spacedim> > >                         output) const
+{
+  transform_differential_forms(input, mapping_type, mapping_data, output);
+}
+
+
+
+template<int dim, int spacedim>
+void
+MappingQGeneric<dim,spacedim>::
+transform (const VectorSlice<const std::vector<Tensor<2, dim> > >   input,
+           const MappingType                                        mapping_type,
+           const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
+           VectorSlice<std::vector<Tensor<2, spacedim> > >          output) const
+{
+  switch (mapping_type)
+    {
+    case mapping_contravariant:
+      transform_fields(input, mapping_type, mapping_data, output);
+      return;
+
+    case mapping_piola_gradient:
+    case mapping_contravariant_gradient:
+    case mapping_covariant_gradient:
+      transform_gradients(input, mapping_type, mapping_data, output);
+      return;
+    default:
+      Assert(false, ExcNotImplemented());
+    }
+}
+
+
+
+template<int dim, int spacedim>
+void
+MappingQGeneric<dim,spacedim>::
+transform (const VectorSlice<const std::vector< DerivativeForm<2, dim, spacedim> > > input,
+           const MappingType                                                         mapping_type,
+           const typename Mapping<dim,spacedim>::InternalDataBase                   &mapping_data,
+           VectorSlice<std::vector<Tensor<3,spacedim> > >                            output) const
+{
+
+  AssertDimension (input.size(), output.size());
+  Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
+          ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(mapping_data);
+
+  switch (mapping_type)
+    {
+    case mapping_covariant_gradient:
+    {
+      Assert (data.update_each & update_contravariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+
+      for (unsigned int q=0; q<output.size(); ++q)
+        for (unsigned int i=0; i<spacedim; ++i)
+          for (unsigned int j=0; j<spacedim; ++j)
+            for (unsigned int k=0; k<spacedim; ++k)
+              {
+                output[q][i][j][k] = data.covariant[q][j][0]
+                                     * data.covariant[q][k][0]
+                                     * input[q][i][0][0];
+                for (unsigned int J=0; J<dim; ++J)
+                  {
+                    const unsigned int K0 = (0==J)? 1 : 0;
+                    for (unsigned int K=K0; K<dim; ++K)
+                      output[q][i][j][k] += data.covariant[q][j][J]
+                                            * data.covariant[q][k][K]
+                                            * input[q][i][J][K];
+                  }
+
+              }
+      return;
+    }
+
+    default:
+      Assert(false, ExcNotImplemented());
+    }
+}
+
+
+
+template<int dim, int spacedim>
+void
+MappingQGeneric<dim,spacedim>::
+transform (const VectorSlice<const std::vector< Tensor<3,dim> > >   input,
+           const MappingType                                        mapping_type,
+           const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
+           VectorSlice<std::vector<Tensor<3,spacedim> > >           output) const
+{
+
+  switch (mapping_type)
+    {
+    case mapping_piola_hessian:
+    case mapping_contravariant_hessian:
+    case mapping_covariant_hessian:
+      transform_hessians(input, mapping_type, mapping_data, output);
+      return;
+    default:
+      Assert(false, ExcNotImplemented());
+    }
+
+}
+
+
+
+//--------------------------- Explicit instantiations -----------------------
+#include "mapping_q_generic.inst"
+
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/fe/mapping_q_generic.inst.in b/source/fe/mapping_q_generic.inst.in
new file mode 100644 (file)
index 0000000..06156a2
--- /dev/null
@@ -0,0 +1,23 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 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.
+//
+// ---------------------------------------------------------------------
+
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+  {
+#if deal_II_dimension <= deal_II_space_dimension       
+    template class MappingQGeneric<deal_II_dimension, deal_II_space_dimension>;
+#endif
+  }
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.