]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Convert data structures in MF::MappingInfo to a face-compatible variant.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 7 Apr 2018 21:15:28 +0000 (23:15 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 11 Apr 2018 11:43:20 +0000 (13:43 +0200)
include/deal.II/matrix_free/dof_info.templates.h
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/helper_functions.h
include/deal.II/matrix_free/mapping_data_on_the_fly.h
include/deal.II/matrix_free/mapping_info.h
include/deal.II/matrix_free/mapping_info.templates.h
include/deal.II/matrix_free/matrix_free.h
include/deal.II/matrix_free/matrix_free.templates.h
source/matrix_free/matrix_free.inst.in
tests/matrix_free/compress_mapping.cc

index 0f0062c3e5c14059bfe72259e60124ebe1050263..894167a1ca5a3937e8c49492b374d92d2f132d7c 100644 (file)
@@ -23,6 +23,7 @@
 #include <deal.II/lac/sparsity_pattern.h>
 #include <deal.II/matrix_free/dof_info.h>
 #include <deal.II/matrix_free/helper_functions.h>
+#include <deal.II/matrix_free/mapping_info.h>
 
 DEAL_II_NAMESPACE_OPEN
 
index d4e83961fb227957c5a1ff34d857d22d0ab9ccc5..8b6a6579fa54d301c73eeba7d9b5c738ae40e686 100644 (file)
@@ -135,17 +135,23 @@ public:
   void reinit (const typename Triangulation<dim>::cell_iterator &cell);
 
   /**
-   * For the transformation information stored in MappingInfo, this function
-   * returns the index which belongs to the current cell as specified in @p
-   * reinit. Note that MappingInfo has different fields for Cartesian cells,
-   * cells with affine mapping and with general mappings, so in order to
-   * access the correct data, this interface must be used together with
-   * get_cell_type.
+   * @deprecated Use get_mapping_data_index_offset() instead.
    */
+  DEAL_II_DEPRECATED
   unsigned int get_cell_data_number() const;
 
   /**
-   * Return the type of the cell the @p reinit function has been called for.
+   * Return the index offset within the geometry fields for the cell the @p
+   * reinit() function has been called for. This index can be used to access
+   * an index into a field that has the same compression behavior as the
+   * Jacobian of the geometry, e.g., to store an effective coefficient tensors
+   * that combines a coefficient with the geometry for lower memory transfer
+   * as the available data fields.
+   */
+  unsigned int get_mapping_data_index_offset() const;
+
+  /**
+   * Return the type of the cell the @p reinit() function has been called for.
    * Valid values are @p cartesian for Cartesian cells (which allows for
    * considerable data compression), @p affine for cells with affine mappings,
    * and @p general for general cells without any compressed storage applied.
@@ -481,6 +487,15 @@ public:
    */
   VectorizedArray<Number> JxW(const unsigned int q_point) const;
 
+  /**
+   * Gets the inverse and transposed version of Jacobian of the mapping
+   * between the unit to the real cell (representing the covariant
+   * transformation). This is exactly the matrix used internally to transform
+   * the unit cell gradients to gradients on the real cell.
+   */
+  Tensor<2,dim,VectorizedArray<Number> >
+  inverse_jacobian(const unsigned int q_index) const;
+
   //@}
 
   /**
@@ -800,12 +815,12 @@ protected:
   const internal::MatrixFreeFunctions::DoFInfo *dof_info;
 
   /**
-   * Stores a pointer to the underlying transformation data from unit to real
-   * cells for the given quadrature formula specified at construction. Also
-   * contained in matrix_info, but it simplifies code if we store a reference
-   * to it.
+   * Stores a pointer to the underlying transformation data from unit to
+   * real cells for the given quadrature formula specified at construction.
+   * Also contained in matrix_info, but it simplifies code if we store a
+   * reference to it.
    */
-  const internal::MatrixFreeFunctions::MappingInfo<dim,Number> *mapping_info;
+  const internal::MatrixFreeFunctions::MappingInfoStorage<dim,dim,Number> *mapping_data;
 
   /**
    * Stores a pointer to the unit cell shape data, i.e., values, gradients and
@@ -815,12 +830,6 @@ protected:
    */
   const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> *data;
 
-  /**
-   * A pointer to the Cartesian Jacobian information of the present cell. Only
-   * set to a useful value if on a Cartesian cell, otherwise zero.
-   */
-  const Tensor<1,dim,VectorizedArray<Number> > *cartesian_data;
-
   /**
    * A pointer to the Jacobian information of the present cell. Only set to a
    * useful value if on a non-Cartesian cell.
@@ -838,26 +847,7 @@ protected:
   /**
    * A pointer to the quadrature weights of the underlying quadrature formula.
    */
-  const VectorizedArray<Number> *quadrature_weights;
-
-  /**
-   * A pointer to the quadrature points on the present cell.
-   */
-  const Point<dim,VectorizedArray<Number> > *quadrature_points;
-
-  /**
-   * A pointer to the diagonal part of the Jacobian gradient on the present
-   * cell. Only set to a useful value if on a general cell with non-constant
-   * Jacobian.
-   */
-  const Tensor<2,dim,VectorizedArray<Number> > *jacobian_grad;
-
-  /**
-   * A pointer to the upper diagonal part of the Jacobian gradient on the
-   * present cell. Only set to a useful value if on a general cell with non-
-   * constant Jacobian.
-   */
-  const Tensor<1,(dim>1?dim*(dim-1)/2:1),Tensor<1,dim,VectorizedArray<Number> > > * jacobian_grad_upper;
+  const Number *quadrature_weights;
 
   /**
    * After a call to reinit(), stores the number of the cell we are currently
@@ -873,11 +863,6 @@ protected:
    */
   internal::MatrixFreeFunctions::CellType cell_type;
 
-  /**
-   * The stride to access the correct data in MappingInfo.
-   */
-  unsigned int cell_data_number;
-
   /**
    * Debug information to track whether dof values have been initialized
    * before accessed. Used to control exceptions when uninitialized data is
@@ -926,12 +911,6 @@ protected:
    */
   std::shared_ptr<internal::MatrixFreeFunctions::MappingDataOnTheFly<dim,Number> > mapped_geometry;
 
-  /**
-   * For use with on-the-fly evaluation, provide a data structure to store the
-   * global dof indices on the current cell from a reinit call.
-   */
-  std::vector<types::global_dof_index> old_style_dof_indices;
-
   /**
    * For a FiniteElement with more than one finite element, select at which
    * component this data structure should start.
@@ -2112,28 +2091,22 @@ FEEvaluationBase<dim,n_components_,Number>
                       :
                       0),
   active_quad_index  (fe_degree != numbers::invalid_unsigned_int ?
-                      data_in.get_mapping_info().
-                      mapping_data_gen[quad_no_in].
-                      quad_index_from_n_q_points(n_q_points)
+                      (data_in.get_mapping_info().cell_data[quad_no_in].
+                       quad_index_from_n_q_points(n_q_points))
                       :
                       0),
   matrix_info        (&data_in),
   dof_info           (&data_in.get_dof_info(fe_no_in)),
-  mapping_info       (&data_in.get_mapping_info()),
+  mapping_data       (internal::MatrixFreeFunctions::MappingInfoCellsOrFaces<dim,Number,false>::get
+                      (data_in.get_mapping_info(), quad_no)),
   data               (&data_in.get_shape_info
                       (fe_no_in, quad_no_in, active_fe_index,
                        active_quad_index)),
-  cartesian_data     (nullptr),
   jacobian           (nullptr),
   J_value            (nullptr),
-  quadrature_weights (mapping_info->mapping_data_gen[quad_no].
-                      quadrature_weights[active_quad_index].begin()),
-  quadrature_points  (nullptr),
-  jacobian_grad      (nullptr),
-  jacobian_grad_upper(nullptr),
+  quadrature_weights (mapping_data->descriptor[active_quad_index].quadrature_weights.begin()),
   cell               (numbers::invalid_unsigned_int),
-  cell_type          (internal::MatrixFreeFunctions::undefined),
-  cell_data_number   (numbers::invalid_unsigned_int),
+  cell_type          (internal::MatrixFreeFunctions::general),
   dof_values_initialized    (false),
   values_quad_initialized   (false),
   gradients_quad_initialized(false),
@@ -2150,7 +2123,7 @@ FEEvaluationBase<dim,n_components_,Number>
   AssertDimension (data->dofs_per_component_on_cell*n_fe_components,
                    dof_info->dofs_per_cell[active_fe_index]);
   AssertDimension (data->n_q_points,
-                   mapping_info->mapping_data_gen[quad_no].n_q_points[active_quad_index]);
+                   mapping_data->descriptor[active_quad_index].n_q_points);
   Assert (n_fe_components == 1 ||
           n_components == 1 ||
           n_components == n_fe_components,
@@ -2183,19 +2156,14 @@ FEEvaluationBase<dim,n_components_,Number>
   active_quad_index  (numbers::invalid_unsigned_int),
   matrix_info        (nullptr),
   dof_info           (nullptr),
-  mapping_info       (nullptr),
+  mapping_data       (nullptr),
   // select the correct base element from the given FE component
   data               (new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>(quadrature, fe, fe.component_to_base_index(first_selected_component).first)),
-  cartesian_data     (nullptr),
   jacobian           (nullptr),
   J_value            (nullptr),
   quadrature_weights (nullptr),
-  quadrature_points  (nullptr),
-  jacobian_grad      (nullptr),
-  jacobian_grad_upper(nullptr),
   cell               (0),
   cell_type          (internal::MatrixFreeFunctions::general),
-  cell_data_number   (numbers::invalid_unsigned_int),
   dof_values_initialized    (false),
   values_quad_initialized   (false),
   gradients_quad_initialized(false),
@@ -2218,9 +2186,10 @@ FEEvaluationBase<dim,n_components_,Number>
     mapped_geometry
       = std::make_shared<internal::MatrixFreeFunctions::MappingDataOnTheFly<dim,Number> >
         (mapping, quadrature, update_flags);
-  jacobian = mapped_geometry->get_inverse_jacobians().begin();
-  J_value = mapped_geometry->get_JxW_values().begin();
-  quadrature_points = mapped_geometry->get_quadrature_points().begin();
+
+  mapping_data = &mapped_geometry->get_data_storage();
+  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
+  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
 
   Assert(fe.element_multiplicity(base_element_number) == 1 ||
          fe.element_multiplicity(base_element_number)-first_selected_component >= n_components_,
@@ -2245,24 +2214,16 @@ FEEvaluationBase<dim,n_components_,Number>
   active_quad_index  (other.active_quad_index),
   matrix_info        (other.matrix_info),
   dof_info           (other.dof_info),
-  mapping_info       (other.mapping_info),
+  mapping_data       (other.mapping_data),
   data               (other.matrix_info == nullptr ?
                       new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>(*other.data) :
                       other.data),
-  cartesian_data     (nullptr),
   jacobian           (nullptr),
   J_value            (nullptr),
-  quadrature_weights (mapping_info != nullptr ?
-                      mapping_info->mapping_data_gen[quad_no].
-                      quadrature_weights[active_quad_index].begin()
-                      :
-                      nullptr),
-  quadrature_points  (nullptr),
-  jacobian_grad      (nullptr),
-  jacobian_grad_upper(nullptr),
+  quadrature_weights (other.matrix_info == nullptr ? nullptr :
+                      mapping_data->descriptor[active_quad_index].quadrature_weights.begin()),
   cell               (numbers::invalid_unsigned_int),
   cell_type          (internal::MatrixFreeFunctions::general),
-  cell_data_number   (numbers::invalid_unsigned_int),
   dof_values_initialized    (false),
   values_quad_initialized   (false),
   gradients_quad_initialized(false),
@@ -2281,10 +2242,11 @@ FEEvaluationBase<dim,n_components_,Number>
        MappingDataOnTheFly<dim,Number>(other.mapped_geometry->get_fe_values().get_mapping(),
                                        other.mapped_geometry->get_quadrature(),
                                        other.mapped_geometry->get_fe_values().get_update_flags()));
-      jacobian = mapped_geometry->get_inverse_jacobians().begin();
-      J_value = mapped_geometry->get_JxW_values().begin();
-      quadrature_points = mapped_geometry->get_quadrature_points().begin();
+      mapping_data = &mapped_geometry->get_data_storage();
       cell = 0;
+
+      jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
+      J_value = mapped_geometry->get_data_storage().JxW_values.begin();
     }
 }
 
@@ -2315,7 +2277,7 @@ FEEvaluationBase<dim,n_components_,Number>
 
   matrix_info = other.matrix_info;
   dof_info = other.dof_info;
-  mapping_info = other.mapping_info;
+  mapping_data = other.mapping_data;
   if (other.matrix_info == nullptr)
     {
       data = new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>(*other.data);
@@ -2328,20 +2290,14 @@ FEEvaluationBase<dim,n_components_,Number>
     }
   set_data_pointers();
 
-  cartesian_data = nullptr;
   jacobian = nullptr;
   J_value = nullptr;
-  quadrature_weights = mapping_info != nullptr ?
-                       mapping_info->mapping_data_gen[quad_no].
-                       quadrature_weights[active_quad_index].begin()
-                       :
-                       nullptr;
-  quadrature_points = nullptr;
-  jacobian_grad = nullptr;
-  jacobian_grad_upper = nullptr;
+  quadrature_weights = (mapping_data != nullptr ?
+                        mapping_data->descriptor[active_quad_index].quadrature_weights.begin()
+                        :
+                        nullptr);
   cell = numbers::invalid_unsigned_int;
   cell_type = internal::MatrixFreeFunctions::general;
-  cell_data_number = numbers::invalid_unsigned_int;
 
   // Create deep copy of mapped geometry for use in parallel...
   if (other.mapped_geometry.get() != nullptr)
@@ -2351,10 +2307,10 @@ FEEvaluationBase<dim,n_components_,Number>
        MappingDataOnTheFly<dim,Number>(other.mapped_geometry->get_fe_values().get_mapping(),
                                        other.mapped_geometry->get_quadrature(),
                                        other.mapped_geometry->get_fe_values().get_update_flags()));
-      jacobian = mapped_geometry->get_inverse_jacobians().begin();
-      J_value = mapped_geometry->get_JxW_values().begin();
-      quadrature_points = mapped_geometry->get_quadrature_points().begin();
       cell = 0;
+      mapping_data = &mapped_geometry->get_data_storage();
+      jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
+      J_value = mapped_geometry->get_data_storage().JxW_values.begin();
     }
 
   return *this;
@@ -2427,72 +2383,22 @@ FEEvaluationBase<dim,n_components_,Number>
 template <int dim, int n_components_, typename Number>
 inline
 void
-FEEvaluationBase<dim,n_components_,Number>::reinit (const unsigned int cell_in)
+FEEvaluationBase<dim,n_components_,Number>::reinit (const unsigned int cell_index)
 {
   Assert (mapped_geometry == nullptr,
           ExcMessage("FEEvaluation was initialized without a matrix-free object."
                      " Integer indexing is not possible"));
   if (mapped_geometry != nullptr)
     return;
-  Assert (dof_info != nullptr, ExcNotInitialized());
-  Assert (mapping_info != nullptr, ExcNotInitialized());
-  AssertIndexRange (cell_in, dof_info->row_starts.size()-1);
-  AssertDimension (((dof_info->cell_active_fe_index.size() > 0) ?
-                    dof_info->cell_active_fe_index[cell_in] : 0),
-                   active_fe_index);
-  cell = cell_in;
-  cell_type = mapping_info->get_cell_type(cell);
-  cell_data_number = mapping_info->get_cell_data_index(cell);
-
-  if (mapping_info->quadrature_points_initialized == true)
-    {
-      AssertIndexRange (cell_data_number, mapping_info->
-                        mapping_data_gen[quad_no].rowstart_q_points.size());
-      const unsigned int index = mapping_info->mapping_data_gen[quad_no].
-                                 rowstart_q_points[cell];
-      AssertIndexRange (index, mapping_info->mapping_data_gen[quad_no].
-                        quadrature_points.size());
-      quadrature_points =
-        &mapping_info->mapping_data_gen[quad_no].quadrature_points[index];
-    }
 
-  if (cell_type == internal::MatrixFreeFunctions::cartesian)
-    {
-      cartesian_data = &mapping_info->cartesian_data[cell_data_number].first;
-      J_value        = &mapping_info->cartesian_data[cell_data_number].second;
-    }
-  else if (cell_type == internal::MatrixFreeFunctions::affine)
-    {
-      jacobian  = &mapping_info->affine_data[cell_data_number].first;
-      J_value   = &mapping_info->affine_data[cell_data_number].second;
-    }
-  else
-    {
-      const unsigned int rowstart = mapping_info->
-                                    mapping_data_gen[quad_no].rowstart_jacobians[cell_data_number];
-      AssertIndexRange (rowstart, mapping_info->
-                        mapping_data_gen[quad_no].jacobians.size());
-      jacobian =
-        &mapping_info->mapping_data_gen[quad_no].jacobians[rowstart];
-      if (mapping_info->JxW_values_initialized == true)
-        {
-          AssertIndexRange (rowstart, mapping_info->
-                            mapping_data_gen[quad_no].JxW_values.size());
-          J_value = &(mapping_info->mapping_data_gen[quad_no].
-                      JxW_values[rowstart]);
-        }
-      if (mapping_info->second_derivatives_initialized == true)
-        {
-          AssertIndexRange(rowstart, mapping_info->
-                           mapping_data_gen[quad_no].jacobians_grad_diag.size());
-          jacobian_grad = &mapping_info->mapping_data_gen[quad_no].
-                          jacobians_grad_diag[rowstart];
-          AssertIndexRange(rowstart, mapping_info->
-                           mapping_data_gen[quad_no].jacobians_grad_upper.size());
-          jacobian_grad_upper = &mapping_info->mapping_data_gen[quad_no].
-                                jacobians_grad_upper[rowstart];
-        }
-    }
+  Assert (this->dof_info != nullptr, ExcNotInitialized());
+  Assert (this->mapping_data != nullptr, ExcNotInitialized());
+  this->cell = cell_index;
+  this->cell_type = this->matrix_info->get_mapping_info().get_cell_type(cell_index);
+
+  const unsigned int offsets = this->mapping_data->data_index_offsets[cell_index];
+  this->jacobian  = &this->mapping_data->jacobians[0][offsets];
+  this->J_value   = &this->mapping_data->JxW_values[offsets];
 
 #ifdef DEBUG
   dof_values_initialized      = false;
@@ -2550,8 +2456,24 @@ unsigned int
 FEEvaluationBase<dim,n_components_,Number>
 ::get_cell_data_number () const
 {
-  Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
-  return cell_data_number;
+  return get_mapping_data_index_offset();
+}
+
+
+
+template <int dim, int n_components_, typename Number>
+inline
+unsigned int
+FEEvaluationBase<dim,n_components_,Number>
+::get_mapping_data_index_offset () const
+{
+  if (matrix_info == 0)
+    return 0;
+  else
+    {
+      AssertIndexRange(cell, this->mapping_data->data_index_offsets.size());
+      return this->mapping_data->data_index_offsets[cell];
+    }
 }
 
 
@@ -2589,7 +2511,7 @@ FEEvaluationBase<dim,n_components_,Number>
   if (this->cell_type == internal::MatrixFreeFunctions::cartesian ||
       this->cell_type == internal::MatrixFreeFunctions::affine)
     {
-      Assert (this->mapping_info != nullptr, ExcNotImplemented());
+      Assert (this->mapping_data != nullptr, ExcNotImplemented());
       VectorizedArray<Number> J = this->J_value[0];
       for (unsigned int q=0; q<this->data->n_q_points; ++q)
         JxW_values[q] = J * this->quadrature_weights[q];
@@ -2606,11 +2528,12 @@ inline
 VectorizedArray<Number>
 FEEvaluationBase<dim,n_components_,Number>::JxW(const unsigned int q_point) const
 {
+  AssertIndexRange(q_point, data->n_q_points);
   Assert (this->J_value != nullptr, ExcNotInitialized());
   if (this->cell_type == internal::MatrixFreeFunctions::cartesian ||
       this->cell_type == internal::MatrixFreeFunctions::affine)
     {
-      Assert (this->mapping_info != nullptr, ExcInternalError());
+      Assert (this->quadrature_weights != nullptr, ExcInternalError());
       return this->J_value[0] * this->quadrature_weights[q_point];
     }
   else
@@ -2619,6 +2542,22 @@ FEEvaluationBase<dim,n_components_,Number>::JxW(const unsigned int q_point) cons
 
 
 
+template <int dim, int n_components_, typename Number>
+inline
+Tensor<2,dim,VectorizedArray<Number> >
+FEEvaluationBase<dim,n_components_,Number>
+::inverse_jacobian(const unsigned int q_index) const
+{
+  AssertIndexRange(q_index, data->n_q_points);
+  Assert (this->jacobian != nullptr, ExcNotImplemented());
+  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
+    return jacobian[0];
+  else
+    return jacobian[q_index];
+}
+
+
+
 namespace internal
 {
   // write access to generic vectors that have operator ().
@@ -3749,13 +3688,13 @@ FEEvaluationBase<dim,n_components_,Number>
       for (unsigned int comp=0; comp<n_components; comp++)
         for (unsigned int d=0; d<dim; ++d)
           grad_out[comp][d] = (this->gradients_quad[comp][d][q_point] *
-                               cartesian_data[0][d]);
+                               jacobian[0][d][d]);
     }
   // cell with general/affine Jacobian
   else
     {
       const Tensor<2,dim,VectorizedArray<Number> > &jac =
-        this->cell_type == internal::MatrixFreeFunctions::general ?
+        this->cell_type > internal::MatrixFreeFunctions::affine ?
         jacobian[q_point] : jacobian[0];
       for (unsigned int comp=0; comp<n_components; comp++)
         {
@@ -3841,32 +3780,35 @@ FEEvaluationBase<dim,n_components_,Number>
           internal::ExcAccessToUninitializedField());
   AssertIndexRange (q_point, this->data->n_q_points);
 
+  Assert(jacobian != nullptr, ExcNotImplemented());
+  const Tensor<2,dim,VectorizedArray<Number> > &jac =
+    jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point];
+
   Tensor<2,dim,VectorizedArray<Number> > hessian_out [n_components];
 
   // Cartesian cell
   if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
     {
-      const Tensor<1,dim,VectorizedArray<Number> > &jac = cartesian_data[0];
       for (unsigned int comp=0; comp<n_components; comp++)
         for (unsigned int d=0; d<dim; ++d)
           {
             hessian_out[comp][d][d] = (this->hessians_quad[comp][d][q_point] *
-                                       jac[d] * jac[d]);
+                                       jac[d][d] * jac[d][d]);
             switch (dim)
               {
               case 1:
                 break;
               case 2:
                 hessian_out[comp][0][1] = (this->hessians_quad[comp][2][q_point] *
-                                           jac[0] * jac[1]);
+                                           jac[0][0] * jac[1][1]);
                 break;
               case 3:
                 hessian_out[comp][0][1] = (this->hessians_quad[comp][3][q_point] *
-                                           jac[0] * jac[1]);
+                                           jac[0][0] * jac[1][1]);
                 hessian_out[comp][0][2] = (this->hessians_quad[comp][4][q_point] *
-                                           jac[0] * jac[2]);
+                                           jac[0][0] * jac[2][2]);
                 hessian_out[comp][1][2] = (this->hessians_quad[comp][5][q_point] *
-                                           jac[1] * jac[2]);
+                                           jac[1][1] * jac[2][2]);
                 break;
               default:
                 Assert (false, ExcNotImplemented());
@@ -3875,16 +3817,9 @@ FEEvaluationBase<dim,n_components_,Number>
               hessian_out[comp][e][d] = hessian_out[comp][d][e];
           }
     }
-  // cell with general Jacobian
-  else if (this->cell_type == internal::MatrixFreeFunctions::general)
+  // cell with general Jacobian, but constant within the cell
+  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
     {
-      Assert (this->mapping_info->second_derivatives_initialized == true,
-              ExcNotInitialized());
-      const Tensor<2,dim,VectorizedArray<Number> > &jac = jacobian[q_point];
-      const Tensor<2,dim,VectorizedArray<Number> > &jac_grad = jacobian_grad[q_point];
-      const Tensor<1,(dim>1?dim*(dim-1)/2:1),
-            Tensor<1,dim,VectorizedArray<Number> > >
-            & jac_grad_UT = jacobian_grad_upper[q_point];
       for (unsigned int comp=0; comp<n_components; comp++)
         {
           // compute laplacian before the gradient because it needs to access
@@ -3902,18 +3837,8 @@ FEEvaluationBase<dim,n_components_,Number>
                   hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
               }
 
-          // add diagonal part of J' * grad(u)
-          for (unsigned int d=0; d<dim; ++d)
-            for (unsigned int e=0; e<dim; ++e)
-              hessian_out[comp][d][d] += (jac_grad[d][e] *
-                                          this->gradients_quad[comp][e][q_point]);
-
-          // add off-diagonal part of J' * grad(u)
-          for (unsigned int d=0, count=0; d<dim; ++d)
-            for (unsigned int e=d+1; e<dim; ++e, ++count)
-              for (unsigned int f=0; f<dim; ++f)
-                hessian_out[comp][d][e] += (jac_grad_UT[count][f] *
-                                            this->gradients_quad[comp][f][q_point]);
+          // no J' * grad(u) part here because the Jacobian is constant
+          // throughout the cell and hence, its derivative is zero
 
           // take symmetric part
           for (unsigned int d=0; d<dim; ++d)
@@ -3921,10 +3846,11 @@ FEEvaluationBase<dim,n_components_,Number>
               hessian_out[comp][e][d] = hessian_out[comp][d][e];
         }
     }
-  // cell with general Jacobian, but constant within the cell
-  else // if (this->cell_type == internal::MatrixFreeFunctions::affine)
+  // cell with general Jacobian
+  else
     {
-      const Tensor<2,dim,VectorizedArray<Number> > &jac = jacobian[0];
+      const Tensor<1,dim*(dim+1)/2,Tensor<1,dim,VectorizedArray<Number> > > &jac_grad =
+        mapping_data->jacobian_gradients[0][this->get_mapping_data_index_offset()+q_point];
       for (unsigned int comp=0; comp<n_components; comp++)
         {
           // compute laplacian before the gradient because it needs to access
@@ -3942,8 +3868,18 @@ FEEvaluationBase<dim,n_components_,Number>
                   hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
               }
 
-          // no J' * grad(u) part here because the Jacobian is constant
-          // throughout the cell and hence, its derivative is zero
+          // add diagonal part of J' * grad(u)
+          for (unsigned int d=0; d<dim; ++d)
+            for (unsigned int e=0; e<dim; ++e)
+              hessian_out[comp][d][d] += (jac_grad[d][e] *
+                                          this->gradients_quad[comp][e][q_point]);
+
+          // add off-diagonal part of J' * grad(u)
+          for (unsigned int d=0, count=dim; d<dim; ++d)
+            for (unsigned int e=d+1; e<dim; ++e, ++count)
+              for (unsigned int f=0; f<dim; ++f)
+                hessian_out[comp][d][e] += (jac_grad[count][f] *
+                                            this->gradients_quad[comp][f][q_point]);
 
           // take symmetric part
           for (unsigned int d=0; d<dim; ++d)
@@ -3966,24 +3902,23 @@ FEEvaluationBase<dim,n_components_,Number>
           internal::ExcAccessToUninitializedField());
   AssertIndexRange (q_point, this->data->n_q_points);
 
+  Assert(jacobian != nullptr, ExcNotImplemented());
+  const Tensor<2,dim,VectorizedArray<Number> > &jac =
+    jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point];
+
   Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > > hessian_out;
 
   // Cartesian cell
   if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
     {
-      const Tensor<1,dim,VectorizedArray<Number> > &jac = cartesian_data[0];
       for (unsigned int comp=0; comp<n_components; comp++)
         for (unsigned int d=0; d<dim; ++d)
           hessian_out[comp][d] = (this->hessians_quad[comp][d][q_point] *
-                                  jac[d] * jac[d]);
+                                  jac[d][d] * jac[d][d]);
     }
-  // cell with general Jacobian
-  else if (this->cell_type == internal::MatrixFreeFunctions::general)
+  // cell with general Jacobian, but constant within the cell
+  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
     {
-      Assert (this->mapping_info->second_derivatives_initialized == true,
-              ExcNotInitialized());
-      const Tensor<2,dim,VectorizedArray<Number> > &jac = jacobian[q_point];
-      const Tensor<2,dim,VectorizedArray<Number> > &jac_grad = jacobian_grad[q_point];
       for (unsigned int comp=0; comp<n_components; comp++)
         {
           // compute laplacian before the gradient because it needs to access
@@ -4000,17 +3935,13 @@ FEEvaluationBase<dim,n_components_,Number>
               for (unsigned int f=1; f<dim; ++f)
                 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
             }
-
-          for (unsigned int d=0; d<dim; ++d)
-            for (unsigned int e=0; e<dim; ++e)
-              hessian_out[comp][d] += (jac_grad[d][e] *
-                                       this->gradients_quad[comp][e][q_point]);
         }
     }
-  // cell with general Jacobian, but constant within the cell
-  else // if (this->cell_type == internal::MatrixFreeFunctions::affine)
+  // cell with general Jacobian
+  else
     {
-      const Tensor<2,dim,VectorizedArray<Number> > &jac = jacobian[0];
+      const Tensor<1,dim*(dim+1)/2,Tensor<1,dim,VectorizedArray<Number> > > &jac_grad =
+        mapping_data->jacobian_gradients[0][this->get_mapping_data_index_offset()+q_point];
       for (unsigned int comp=0; comp<n_components; comp++)
         {
           // compute laplacian before the gradient because it needs to access
@@ -4027,6 +3958,11 @@ FEEvaluationBase<dim,n_components_,Number>
               for (unsigned int f=1; f<dim; ++f)
                 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
             }
+
+          for (unsigned int d=0; d<dim; ++d)
+            for (unsigned int e=0; e<dim; ++e)
+              hessian_out[comp][d] += (jac_grad[d][e] *
+                                       this->gradients_quad[comp][e][q_point]);
         }
     }
   return hessian_out;
@@ -4114,22 +4050,25 @@ FEEvaluationBase<dim,n_components_,Number>
   Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
   AssertIndexRange (q_point, this->data->n_q_points);
   this->gradients_quad_submitted = true;
+  Assert (this->J_value != nullptr, ExcNotInitialized());
+  Assert (this->jacobian != nullptr, ExcNotInitialized());
 #endif
+
   if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
     {
       const VectorizedArray<Number> JxW = J_value[0] * quadrature_weights[q_point];
       for (unsigned int comp=0; comp<n_components; comp++)
         for (unsigned int d=0; d<dim; ++d)
           this->gradients_quad[comp][d][q_point] = (grad_in[comp][d] *
-                                                    cartesian_data[0][d] * JxW);
+                                                    jacobian[0][d][d] * JxW);
     }
   else
     {
       const Tensor<2,dim,VectorizedArray<Number> > &jac =
-        this->cell_type == internal::MatrixFreeFunctions::general ?
+        this->cell_type > internal::MatrixFreeFunctions::affine ?
         jacobian[q_point] : jacobian[0];
       const VectorizedArray<Number> JxW =
-        this->cell_type == internal::MatrixFreeFunctions::general ?
+        this->cell_type > internal::MatrixFreeFunctions::affine ?
         J_value[q_point] : J_value[0] * quadrature_weights[q_point];
       for (unsigned int comp=0; comp<n_components; ++comp)
         for (unsigned int d=0; d<dim; ++d)
@@ -4320,26 +4259,27 @@ FEEvaluationAccess<dim,1,Number>
           internal::ExcAccessToUninitializedField());
   AssertIndexRange (q_point, this->data->n_q_points);
 
+  Assert (this->jacobian != nullptr, ExcNotInitialized());
+
   Tensor<1,dim,VectorizedArray<Number> > grad_out;
 
-  // Cartesian cell
   if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
     {
       for (unsigned int d=0; d<dim; ++d)
         grad_out[d] = (this->gradients_quad[0][d][q_point] *
-                       this->cartesian_data[0][d]);
+                       this->jacobian[0][d][d]);
     }
-  // cell with general/constant Jacobian
+  // cell with general/affine Jacobian
   else
     {
       const Tensor<2,dim,VectorizedArray<Number> > &jac =
-        this->cell_type == internal::MatrixFreeFunctions::general ?
-        this->jacobian[q_point] : this->jacobian[0];
+        this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
+                       q_point : 0];
       for (unsigned int d=0; d<dim; ++d)
         {
-          grad_out[d] = (jac[d][0] * this->gradients_quad[0][0][q_point]);
+          grad_out[d] = jac[d][0] * this->gradients_quad[0][0][q_point];
           for (unsigned int e=1; e<dim; ++e)
-            grad_out[d] += (jac[d][e] * this->gradients_quad[0][e][q_point]);
+            grad_out[d] += jac[d][e] * this->gradients_quad[0][e][q_point];
         }
     }
   return grad_out;
@@ -4427,36 +4367,40 @@ inline
 void
 FEEvaluationAccess<dim,1,Number>
 ::submit_gradient (const Tensor<1,dim,VectorizedArray<Number> > grad_in,
-                   const unsigned int q_point)
+                   const unsigned int q_index)
 {
 #ifdef DEBUG
   Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
-  AssertIndexRange (q_point, this->data->n_q_points);
+  AssertIndexRange (q_index, this->data->n_q_points);
   this->gradients_quad_submitted = true;
+  Assert (this->J_value != nullptr, ExcNotInitialized());
+  Assert (this->jacobian != nullptr, ExcNotInitialized());
 #endif
+
   if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
     {
-      const VectorizedArray<Number> JxW = this->J_value[0] * this->quadrature_weights[q_point];
+      const VectorizedArray<Number> JxW = this->J_value[0] *
+                                          this->quadrature_weights[q_index];
       for (unsigned int d=0; d<dim; ++d)
-        this->gradients_quad[0][d][q_point] = (grad_in[d] *
-                                               this->cartesian_data[0][d] *
+        this->gradients_quad[0][d][q_index] = (grad_in[d] *
+                                               this->jacobian[0][d][d] *
                                                JxW);
     }
   // general/affine cell type
   else
     {
       const Tensor<2,dim,VectorizedArray<Number> > &jac =
-        this->cell_type == internal::MatrixFreeFunctions::general ?
-        this->jacobian[q_point] : this->jacobian[0];
+        this->cell_type > internal::MatrixFreeFunctions::affine ?
+        this->jacobian[q_index] : this->jacobian[0];
       const VectorizedArray<Number> JxW =
-        this->cell_type == internal::MatrixFreeFunctions::general ?
-        this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point];
+        this->cell_type > internal::MatrixFreeFunctions::affine ?
+        this->J_value[q_index] : this->J_value[0] * this->quadrature_weights[q_index];
       for (unsigned int d=0; d<dim; ++d)
         {
           VectorizedArray<Number> new_val = jac[0][d] * grad_in[0];
           for (unsigned int e=1; e<dim; ++e)
             new_val += jac[e][d] * grad_in[e];
-          this->gradients_quad[0][d][q_point] = new_val * JxW;
+          this->gradients_quad[0][d][q_index] = new_val * JxW;
         }
     }
 }
@@ -4552,6 +4496,7 @@ FEEvaluationAccess<dim,dim,Number>
   Assert (this->gradients_quad_initialized==true,
           internal::ExcAccessToUninitializedField());
   AssertIndexRange (q_point, this->data->n_q_points);
+  Assert (this->jacobian != nullptr, ExcNotInitialized());
 
   VectorizedArray<Number> divergence;
 
@@ -4559,10 +4504,10 @@ FEEvaluationAccess<dim,dim,Number>
   if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
     {
       divergence = (this->gradients_quad[0][0][q_point] *
-                    this->cartesian_data[0][0]);
+                    this->jacobian[0][0][0]);
       for (unsigned int d=1; d<dim; ++d)
         divergence += (this->gradients_quad[d][d][q_point] *
-                       this->cartesian_data[0][d]);
+                       this->jacobian[0][d][d]);
     }
   // cell with general/constant Jacobian
   else
@@ -4714,7 +4659,10 @@ FEEvaluationAccess<dim,dim,Number>
   Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
   AssertIndexRange (q_point, this->data->n_q_points);
   this->gradients_quad_submitted = true;
+  Assert (this->J_value != nullptr, ExcNotInitialized());
+  Assert (this->jacobian != nullptr, ExcNotInitialized());
 #endif
+
   if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
     {
       const VectorizedArray<Number> fac = this->J_value[0] *
@@ -4722,7 +4670,7 @@ FEEvaluationAccess<dim,dim,Number>
       for (unsigned int d=0; d<dim; ++d)
         {
           this->gradients_quad[d][d][q_point] = (fac *
-                                                 this->cartesian_data[0][d]);
+                                                 this->jacobian[0][d][d]);
           for (unsigned int e=d+1; e<dim; ++e)
             {
               this->gradients_quad[d][e][q_point] = VectorizedArray<Number>();
@@ -4764,22 +4712,25 @@ FEEvaluationAccess<dim,dim,Number>
   Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
   AssertIndexRange (q_point, this->data->n_q_points);
   this->gradients_quad_submitted = true;
+  Assert (this->J_value != nullptr, ExcNotInitialized());
+  Assert (this->jacobian != nullptr, ExcNotInitialized());
 #endif
+
   if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
     {
       const VectorizedArray<Number> JxW = this->J_value[0] * this->quadrature_weights[q_point];
       for (unsigned int d=0; d<dim; ++d)
         this->gradients_quad[d][d][q_point] = (sym_grad.access_raw_entry(d) *
                                                JxW *
-                                               this->cartesian_data[0][d]);
+                                               this->jacobian[0][d][d]);
       for (unsigned int e=0, counter=dim; e<dim; ++e)
         for (unsigned int d=e+1; d<dim; ++d, ++counter)
           {
             const VectorizedArray<Number> value = sym_grad.access_raw_entry(counter) * JxW;
             this->gradients_quad[e][d][q_point] = (value *
-                                                   this->cartesian_data[0][d]);
+                                                   this->jacobian[0][d][d]);
             this->gradients_quad[d][e][q_point] = (value *
-                                                   this->cartesian_data[0][e]);
+                                                   this->jacobian[0][e][e]);
           }
     }
   // general/affine cell type
@@ -4943,23 +4894,13 @@ FEEvaluationAccess<1,1,Number>
           internal::ExcAccessToUninitializedField());
   AssertIndexRange (q_point, this->data->n_q_points);
 
-  Tensor<1,1,VectorizedArray<Number> > grad_out;
+  const Tensor<2,1,VectorizedArray<Number> > &jac =
+    this->cell_type == internal::MatrixFreeFunctions::general ?
+    this->jacobian[q_point] : this->jacobian[0];
 
-  // Cartesian cell
-  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
-    {
-      grad_out[0] = (this->gradients_quad[0][0][q_point] *
-                     this->cartesian_data[0][0]);
-    }
-  // cell with general/constant Jacobian
-  else
-    {
-      const Tensor<2,1,VectorizedArray<Number> > &jac =
-        this->cell_type == internal::MatrixFreeFunctions::general ?
-        this->jacobian[q_point] : this->jacobian[0];
+  Tensor<1,1,VectorizedArray<Number> > grad_out;
+  grad_out[0] = jac[0][0] * this->gradients_quad[0][0][q_point];
 
-      grad_out[0] = (jac[0][0] * this->gradients_quad[0][0][q_point]);
-    }
   return grad_out;
 }
 
@@ -5052,25 +4993,15 @@ FEEvaluationAccess<1,1,Number>
   AssertIndexRange (q_point, this->data->n_q_points);
   this->gradients_quad_submitted = true;
 #endif
-  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
-    {
-      const VectorizedArray<Number> JxW = this->J_value[0] * this->quadrature_weights[q_point];
-      this->gradients_quad[0][0][q_point] = (grad_in[0] *
-                                             this->cartesian_data[0][0] *
-                                             JxW);
-    }
-  // general/affine cell type
-  else
-    {
-      const Tensor<2,1,VectorizedArray<Number> > &jac =
-        this->cell_type == internal::MatrixFreeFunctions::general ?
-        this->jacobian[q_point] : this->jacobian[0];
-      const VectorizedArray<Number> JxW =
-        this->cell_type == internal::MatrixFreeFunctions::general ?
-        this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point];
 
-      this->gradients_quad[0][0][q_point] = jac[0][0] * grad_in[0] * JxW;
-    }
+  const Tensor<2,1,VectorizedArray<Number> > &jac =
+    this->cell_type > internal::MatrixFreeFunctions::affine ?
+    this->jacobian[q_point] : this->jacobian[0];
+  const VectorizedArray<Number> JxW =
+    this->cell_type > internal::MatrixFreeFunctions::affine ?
+    this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point];
+
+  this->gradients_quad[0][0][q_point] = jac[0][0] * grad_in[0] * JxW;
 }
 
 
@@ -5252,13 +5183,13 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
                   proposed_dof_comp = no;
                   break;
                 }
-          if (static_n_q_points ==
-              this->mapping_info->mapping_data_gen[this->quad_no].n_q_points[this->active_quad_index])
+          if (n_q_points ==
+              this->mapping_data->descriptor[this->active_quad_index].n_q_points)
             proposed_quad_comp = this->quad_no;
           else
-            for (unsigned int no=0; no<this->mapping_info->mapping_data_gen.size(); ++no)
-              if (this->mapping_info->mapping_data_gen[no].n_q_points[this->active_quad_index]
-                  == static_n_q_points)
+            for (unsigned int no=0; no<this->matrix_info->get_mapping_info().cell_data.size(); ++no)
+              if (this->matrix_info->get_mapping_info().cell_data[no].descriptor[this->active_quad_index].n_q_points
+                  == n_q_points)
                 {
                   proposed_quad_comp = no;
                   break;
@@ -5326,8 +5257,7 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
   if (fe_no != numbers::invalid_unsigned_int)
     {
       AssertDimension (n_q_points,
-                       this->mapping_info->mapping_data_gen[this->quad_no].
-                       n_q_points[this->active_quad_index]);
+                       this->mapping_data->descriptor[this->active_quad_index].n_q_points);
       AssertDimension (this->data->dofs_per_component_on_cell * this->n_fe_components,
                        this->dof_info->dofs_per_cell[this->active_fe_index]);
     }
@@ -5343,10 +5273,24 @@ Point<dim,VectorizedArray<Number> >
 FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 ::quadrature_point (const unsigned int q) const
 {
-  Assert (this->mapping_info->quadrature_points_initialized == true,
-          ExcNotInitialized());
-  Assert (this->quadrature_points != nullptr, ExcNotInitialized());
+  if (this->matrix_info == nullptr)
+    {
+      Assert ((this->mapped_geometry->get_fe_values().get_update_flags() |
+               update_quadrature_points),
+              ExcNotInitialized());
+    }
+  else
+    {
+      Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
+             ExcNotInitialized());
+    }
+
   AssertIndexRange (q, n_q_points);
+  const Point<dim,VectorizedArray<Number> > *quadrature_points = &this->mapping_data->
+      quadrature_points[this->mapping_data->quadrature_point_offsets[this->cell]];
+
+  const unsigned int n_q_points_1d_actual =
+    fe_degree == -1 ? this->data->n_q_points_1d : n_q_points_1d;
 
   // Cartesian mesh: not all quadrature points are stored, only the
   // diagonal. Hence, need to find the tensor product index and retrieve the
@@ -5357,15 +5301,15 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
       switch (dim)
         {
         case 1:
-          return this->quadrature_points[q];
+          return quadrature_points[q];
         case 2:
-          point[0] = this->quadrature_points[q%n_q_points_1d][0];
-          point[1] = this->quadrature_points[q/n_q_points_1d][1];
+          point[0] = quadrature_points[q%n_q_points_1d_actual][0];
+          point[1] = quadrature_points[q/n_q_points_1d_actual][1];
           return point;
         case 3:
-          point[0] = this->quadrature_points[q%n_q_points_1d][0];
-          point[1] = this->quadrature_points[(q/n_q_points_1d)%n_q_points_1d][1];
-          point[2] = this->quadrature_points[q/(n_q_points_1d*n_q_points_1d)][2];
+          point[0] = quadrature_points[q%n_q_points_1d_actual][0];
+          point[1] = quadrature_points[(q/n_q_points_1d_actual)%n_q_points_1d_actual][1];
+          point[2] = quadrature_points[q/(n_q_points_1d_actual*n_q_points_1d_actual)][2];
           return point;
         default:
           Assert (false, ExcNotImplemented());
@@ -5374,7 +5318,7 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
     }
   // all other cases: just return the respective data as it is fully stored
   else
-    return this->quadrature_points[q];
+    return quadrature_points[q];
 }
 
 
index 8cb01efebb52743201bc2db776ec86f658971ca8..1720820106bdadb4076666f1f7a59b2ee10a82bd 100644 (file)
@@ -132,64 +132,6 @@ namespace internal
       unsigned int n_procs;
     };
 
-    /**
-     * Data type to identify cell type.
-     */
-    enum CellType
-    {
-      /**
-       * The cell is Cartesian.
-       */
-      cartesian = 0,
-      /**
-       * The cell may be described with an affine mapping.
-       */
-      affine    = 1,
-      /**
-       * There is no special information available for compressing the
-       * representation of the cell.
-       */
-      general   = 2,
-      /**
-       * The cell type is undefined.
-       */
-      undefined = 3
-    };
-
-
-    /**
-     * A class that is used to compare floating point arrays (e.g.
-     * std::vectors, Tensor<1,dim>, etc.). The idea of this class is to
-     * consider two arrays as equal if they are the same within a given
-     * tolerance. We use this comparator class within an std::map<> of the
-     * given arrays. Note that this comparison operator does not satisfy all
-     * the mathematical properties one usually wants to have (consider e.g.
-     * the numbers a=0, b=0.1, c=0.2 with tolerance 0.15; the operator gives
-     * a<c, but neither of a<b? or b<c? is satisfied). This is not a problem
-     * in the use cases for this class, but be careful when using it in other
-     * contexts.
-     */
-    template <typename Number>
-    struct FPArrayComparator
-    {
-      FPArrayComparator (const Number scaling);
-
-      bool operator() (const std::vector<Number> &v1,
-                       const std::vector<Number> &v2) const;
-
-      template <int dim>
-      bool operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
-                       const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const;
-
-      template <int dim>
-      bool operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
-                       const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const;
-
-      Number tolerance;
-    };
-
-    // Note: Implementation in matrix_free.templates.h
-
   } // end of namespace MatrixFreeFunctions
 } // end of namespace internal
 
index a5c196a6057b7b2e665fd179c5dcf88bcfda1811..265731fae6e45e30f0e8c3c88ccba6e2da599b28 100644 (file)
@@ -27,6 +27,7 @@
 #include <deal.II/matrix_free/mapping_info.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/mapping_q1.h>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -101,36 +102,13 @@ namespace internal
       const dealii::FEValues<dim> &get_fe_values () const;
 
       /**
-       * Return a vector of inverse transpose Jacobians. For compatibility
-       * with FEEvaluation, it returns tensors of vectorized arrays, even
-       * though all components are equal.
+       * Return a reference to the underlying storage field of type
+       * MappingInfoStorage of the same format as the data fields in
+       * MappingInfo. This ensures compatibility with the precomputed data
+       * fields in the MappingInfo class.
        */
-      const AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > &
-      get_inverse_jacobians() const;
-
-      /**
-       * Return a vector of quadrature weights times the Jacobian determinant
-       * (JxW). For compatibility with FEEvaluation, it returns tensors of
-       * vectorized arrays, even though all components are equal.
-       */
-      const AlignedVector<VectorizedArray<Number> > &
-      get_JxW_values() const;
-
-      /**
-       * Return a vector of quadrature points in real space on the given cell.
-       * For compatibility with FEEvaluation, it returns tensors of vectorized
-       * arrays, even though all components are equal.
-       */
-      const AlignedVector<Point<dim,VectorizedArray<Number> > > &
-      get_quadrature_points() const;
-
-      /**
-       * Return a vector of normal vectors in real space on the given cell.
-       * For compatibility with FEEvaluation, it returns tensors of vectorized
-       * arrays, even though all components are equal.
-       */
-      const AlignedVector<Tensor<1,dim,VectorizedArray<Number> > > &
-      get_normal_vectors() const;
+      const MappingInfoStorage<dim,dim,Number> &
+      get_data_storage() const;
 
       /**
        * Return a reference to 1D quadrature underlying this object.
@@ -162,24 +140,10 @@ namespace internal
       const Quadrature<1> quadrature_1d;
 
       /**
-       * Inverse Jacobians, stored in vectorized array form.
-       */
-      AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > inverse_jacobians;
-
-      /**
-       * Stored Jacobian determinants and quadrature weights
-       */
-      AlignedVector<VectorizedArray<Number> > jxw_values;
-
-      /**
-       * Stored quadrature points
+       * The storage part created for a single cell and held in analogy to
+       * MappingInfo.
        */
-      AlignedVector<Point<dim,VectorizedArray<Number> > > quadrature_points;
-
-      /**
-       * Stored normal vectors (for face integration)
-       */
-      AlignedVector<Tensor<1,dim,VectorizedArray<Number> > > normal_vectors;
+      MappingInfoStorage<dim,dim,Number> mapping_info_storage;
     };
 
 
@@ -192,13 +156,24 @@ namespace internal
                                                           const UpdateFlags update_flags)
       :
       fe_values(mapping, fe_dummy, Quadrature<dim>(quadrature),
-                internal::MatrixFreeFunctions::MappingInfo<dim,Number>::compute_update_flags(update_flags)),
-      quadrature_1d(quadrature),
-      inverse_jacobians(fe_values.get_quadrature().size()),
-      jxw_values(fe_values.get_quadrature().size()),
-      quadrature_points(fe_values.get_quadrature().size()),
-      normal_vectors(fe_values.get_quadrature().size())
+                MappingInfo<dim,Number>::compute_update_flags(update_flags)),
+      quadrature_1d(quadrature)
     {
+      mapping_info_storage.descriptor.resize(1);
+      mapping_info_storage.descriptor[0].initialize(quadrature);
+      mapping_info_storage.data_index_offsets.resize(1);
+      mapping_info_storage.JxW_values.resize(fe_values.n_quadrature_points);
+      mapping_info_storage.jacobians[0].resize(fe_values.n_quadrature_points);
+      if (update_flags & update_quadrature_points)
+        {
+          mapping_info_storage.quadrature_point_offsets.resize(1, 0);
+          mapping_info_storage.quadrature_points.resize(fe_values.n_quadrature_points);
+        }
+      if (fe_values.get_update_flags() & update_normal_vectors)
+        {
+          mapping_info_storage.normal_vectors.resize(fe_values.n_quadrature_points);
+          mapping_info_storage.normals_times_jacobians[0].resize(fe_values.n_quadrature_points);
+        }
       Assert(!(fe_values.get_update_flags() & update_jacobian_grads),
              ExcNotImplemented());
     }
@@ -210,17 +185,9 @@ namespace internal
     MappingDataOnTheFly<dim,Number>::MappingDataOnTheFly (const Quadrature<1> &quadrature,
                                                           const UpdateFlags update_flags)
       :
-      fe_values(fe_dummy, Quadrature<dim>(quadrature),
-                internal::MatrixFreeFunctions::MappingInfo<dim,Number>::compute_update_flags(update_flags)),
-      quadrature_1d(quadrature),
-      inverse_jacobians(fe_values.get_quadrature().size()),
-      jxw_values(fe_values.get_quadrature().size()),
-      quadrature_points(fe_values.get_quadrature().size()),
-      normal_vectors(fe_values.get_quadrature().size())
-    {
-      Assert(!(fe_values.get_update_flags() & update_jacobian_grads),
-             ExcNotImplemented());
-    }
+      MappingDataOnTheFly(::dealii::StaticMappingQ1<dim,dim>::mapping,
+                          quadrature, update_flags)
+    {}
 
 
 
@@ -235,18 +202,27 @@ namespace internal
       fe_values.reinit(present_cell);
       for (unsigned int q=0; q<fe_values.get_quadrature().size(); ++q)
         {
-          if (fe_values.get_update_flags() & update_inverse_jacobians)
-            for (unsigned int d=0; d<dim; ++d)
-              for (unsigned int e=0; e<dim; ++e)
-                inverse_jacobians[q][d][e] = fe_values.inverse_jacobian(q)[e][d];
+          if (fe_values.get_update_flags() & update_JxW_values)
+            mapping_info_storage.JxW_values[q] = fe_values.JxW(q);
+          if (fe_values.get_update_flags() & update_jacobians)
+            {
+              Tensor<2,dim> jac = fe_values.jacobian(q);
+              jac = invert(transpose(jac));
+              for (unsigned int d=0; d<dim; ++d)
+                for (unsigned int e=0; e<dim; ++e)
+                  mapping_info_storage.jacobians[0][q][d][e] = jac[d][e];
+            }
           if (fe_values.get_update_flags() & update_quadrature_points)
             for (unsigned int d=0; d<dim; ++d)
-              quadrature_points[q][d] = fe_values.quadrature_point(q)[d];
+              mapping_info_storage.quadrature_points[q][d] = fe_values.quadrature_point(q)[d];
           if (fe_values.get_update_flags() & update_normal_vectors)
-            for (unsigned int d=0; d<dim; ++d)
-              normal_vectors[q][d] = fe_values.normal_vector(q)[d];
-          if (fe_values.get_update_flags() & update_JxW_values)
-            jxw_values[q] = fe_values.JxW(q);
+            {
+              for (unsigned int d=0; d<dim; ++d)
+                mapping_info_storage.normal_vectors[q][d] = fe_values.normal_vector(q)[d];
+              mapping_info_storage.normals_times_jacobians[0][q] =
+                mapping_info_storage.normal_vectors[q] *
+                mapping_info_storage.jacobians[0][q];
+            }
         }
     }
 
@@ -284,40 +260,10 @@ namespace internal
 
     template <int dim, typename Number>
     inline
-    const AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > &
-    MappingDataOnTheFly<dim,Number>::get_inverse_jacobians() const
-    {
-      return inverse_jacobians;
-    }
-
-
-
-    template <int dim, typename Number>
-    inline
-    const AlignedVector<Tensor<1,dim,VectorizedArray<Number> > > &
-    MappingDataOnTheFly<dim,Number>::get_normal_vectors() const
-    {
-      return normal_vectors;
-    }
-
-
-
-    template <int dim, typename Number>
-    inline
-    const AlignedVector<Point<dim,VectorizedArray<Number> > > &
-    MappingDataOnTheFly<dim,Number>::get_quadrature_points() const
-    {
-      return quadrature_points;
-    }
-
-
-
-    template <int dim, typename Number>
-    inline
-    const AlignedVector<VectorizedArray<Number> > &
-    MappingDataOnTheFly<dim,Number>::get_JxW_values() const
+    const MappingInfoStorage<dim,dim,Number> &
+    MappingDataOnTheFly<dim,Number>::get_data_storage() const
     {
-      return jxw_values;
+      return mapping_info_storage;
     }
 
 
index d4afeaea8a31a1839fab03e225a12e42c77db222..74da08ff9904b4d9092392c77b11b15766a02a49 100644 (file)
@@ -37,71 +37,268 @@ namespace internal
   namespace MatrixFreeFunctions
   {
     /**
-     * The class that stores all geometry-dependent data related with cell
-     * interiors for use in the matrix-free class.
+     * An enum to identify various types of cells. The most general type is
+     * what we typically compute in the FEValues context but for many
+     * geometries we can save significant storage.
+     */
+    enum CellType : unsigned char
+    {
+      /**
+       * The cell is Cartesian.
+       */
+      cartesian  = 0,
+      /**
+       * The cell may be described with an affine mapping.
+       */
+      affine     = 1,
+      /**
+       * The current face is flat, i.e., the normal factor on a face is the
+       * same on all quadrature points.
+       */
+      flat_faces = 2,
+      /**
+       * There is no special information available for compressing the
+       * representation of the cell.
+       */
+      general    = 3
+    };
+
+
+
+    /**
+     * Definition of a structure that stores all cached data related to the
+     * evaluated geometry from the mapping. In order to support hp-adaptivity
+     * and compressed storage (in particular for Jacobians, JxW values, and
+     * normals), storage length can be different for different rows. Thus, it
+     * allows to jump at the data of individual rows similar to compressed row
+     * storage in sparse matrices. We have two different start indices for
+     * fields with different sizes. The first category of offsets are the
+     * indices for Jacobians of the transformation from unit to real cell (we
+     * store the inverse Jacobian), second derivatives, JxW values, and normal
+     * vectors. We keep separate arrays for all these data structures because
+     * a user code might access only some of them. In such a case, one array
+     * will be gone through in a contiguous order with access to all entries,
+     * which makes it easy for the processor to prefetch data. Having all data
+     * in a single array would require some strides in the access pattern,
+     * which is much more complicated for the processor to predict (and indeed
+     * leads to prefetching of data that does not get used on Intel processors
+     * such as BroadwellEP).
      *
-     * @author Katharina Kormann and Martin Kronbichler, 2010, 2011
+     * The second category of indices are the offsets for the quadrature
+     * points. Quadrature points can be compressed less than the other fields
+     * and thus need longer fields. Quadrature point indices are often used in
+     * other contexts such as evaluation of right hand sides.
+     *
+     * The third component is a descriptor of data from the unit cells, called
+     * QuadratureDescriptor, which contains the quadrature weights and
+     * permutations of how to go through quadrature points in case of face
+     * data. The latter comes in a vector for the support of hp adaptivity,
+     * with several data fields for the individual quadrature formulas.
+     *
+     * @author Katharina Kormann, Martin Kronbichler, 2018
      */
-    template <int dim, typename Number>
-    struct MappingInfo
+    template <int structdim, int spacedim, typename Number>
+    struct MappingInfoStorage
     {
+      struct QuadratureDescriptor
+      {
+        /**
+         * Constructor. Does nothing.
+         */
+        QuadratureDescriptor();
+
+        /**
+         * Set up the lengths in the various members of this struct.
+         */
+        void initialize(const Quadrature<1> &quadrature_1d,
+                        const UpdateFlags    update_flags_inner_faces = update_default);
+
+        /**
+         * Returns the memory consumption in bytes.
+         */
+        std::size_t memory_consumption() const;
+
+        /**
+         * Number of quadrature points applied on the given cell or face.
+         */
+        unsigned int n_q_points;
+
+        /**
+         * Quadrature formula applied on the given cell or face.
+         */
+        Quadrature<structdim> quadrature;
+
+        /**
+         * Quadrature weights separated by dimension for use in specific
+         * situations.
+         */
+        std::array<AlignedVector<Number>, structdim> tensor_quadrature_weights;
+
+        /**
+         * A cached vector of quadrature weights in the given number format
+         * (non-vectorized, as it is cheap to broadcast the value to all lanes
+         * when it is used in a vectorized context).
+         */
+        AlignedVector<Number> quadrature_weights;
+
+        /**
+         * For quadrature on faces, the evaluation of basis functions is not
+         * in the correct order if a face is not in the standard orientation
+         * to a given element. This data structure is used to re-order the
+         * data evaluated on quadrature points to represent the correct order.
+         */
+        dealii::Table<2,unsigned int> face_orientations;
+      };
+
       /**
-       * Determines how many bits of an unsigned int are used to distinguish
-       * the cell types (Cartesian, with constant Jacobian, or general)
+       * A class describing the layout of the sections in the @p data_storage
+       * field and also includes some data that depends on the number of
+       * quadrature points in the hp context such as the inner quadrature
+       * formula and re-indexing for faces that are not in the standard
+       * orientation.
        */
-      static const std::size_t  n_cell_type_bits = 2;
+      std::vector<QuadratureDescriptor> descriptor;
 
       /**
-       * Determines how many types of different cells can be detected at most.
-       * Corresponds to the number of bits we reserved for it.
+       * Stores the index offset into the arrays @p jxw_values, @p jacobians,
+       * @p normal_vectors and the second derivatives. Note that affine cells
+       * have shorter fields of length 1, where the others have lengths equal
+       * to the numer of quadrature points of the given cell.
        */
-      static const unsigned int n_cell_types = 1U<<n_cell_type_bits;
+      AlignedVector<unsigned int> data_index_offsets;
 
       /**
-       * An abbreviation for the length of vector lines of the current data
-       * type.
+       * The storage of the Jacobian determinant (times the quadrature weight
+       * in case the transformation is non-affine) on quadrature
+       * points.
+       *
+       * Indexed by @p data_index_offsets.
        */
-      static const unsigned int n_vector_elements = VectorizedArray<Number>::n_array_elements;
+      AlignedVector<VectorizedArray<Number> > JxW_values;
 
+      /**
+       * Stores the normal vectors.
+       *
+       * Indexed by @p data_index_offsets.
+       */
+      AlignedVector<Tensor<1,spacedim,VectorizedArray<Number> > > normal_vectors;
+
+      /**
+       * The storage of covariant transformation on quadrature points, i.e.,
+       * the inverse and transposed Jacobians of the transformation from the
+       * unit to the real cell.
+       *
+       * Indexed by @p data_index_offsets.
+       *
+       * Contains two fields for access from both sides for interior faces,
+       * but the default case (cell integrals or boundary integrals) only
+       * fills the zeroth component and ignores the first one.
+       */
+      AlignedVector<Tensor<2,spacedim,VectorizedArray<Number> > > jacobians[2];
+
+      /**
+       * The storage of the gradients of the inverse Jacobian
+       * transformation. Because of symmetry, only the upper diagonal and
+       * diagonal part are needed. The first index runs through the
+       * derivatives, starting with the diagonal and then continuing row-wise,
+       * i.e., $\partial^2/\partial x_1 \partial x_2$ first, then
+       * $\partial^2/\partial x_1 \partial x_3$, and so on. The second index
+       * is the spatial coordinate.
+       *
+       * Indexed by @p data_index_offsets.
+       *
+       * Contains two fields for access from both sides for interior faces,
+       * but the default case (cell integrals or boundary integrals) only
+       * fills the zeroth component and ignores the first one.
+       */
+      AlignedVector<Tensor<1,spacedim *(spacedim+1)/2,
+                    Tensor<1,spacedim,VectorizedArray<Number> > > > jacobian_gradients[2];
+
+      /**
+       * Stores the Jacobian transformations times the normal vector (this
+       * represents a shortcut that is accessed often and can thus get higher
+       * performance).
+       *
+       * Indexed by @p data_index_offsets.
+       */
+      AlignedVector<Tensor<1,spacedim,VectorizedArray<Number> > > normals_times_jacobians [2];
+
+      /**
+       * Stores the index offset of a particular cell into the quadrature
+       * points array in real coordinates. Note that Cartesian cells have
+       * shorter fields (length is @p n_q_points_1d) than non-Cartesian cells
+       * (length is @p n_q_points) or faces.
+       */
+      AlignedVector<unsigned int> quadrature_point_offsets;
+
+      /**
+       * Stores the quadrature points in real coordinates, including a
+       * compression scheme for Cartesian cells where we do not need to store
+       * the full data on all points.
+       *
+       * Indexed by @p quadrature_point_offsets.
+       */
+      AlignedVector<Point<spacedim,VectorizedArray<Number> > > quadrature_points;
+
+      /**
+       * Returns the quadrature index for a given number of quadrature
+       * points. If not in hp mode or if the index is not found, this
+       * function always returns index 0. Hence, this function does not
+       * check whether the given degree is actually present.
+       */
+      unsigned int
+      quad_index_from_n_q_points (const unsigned int n_q_points) const;
+
+      /**
+       * Prints a detailed summary of memory consumption in the different
+       * structures of this class to the given output stream.
+       */
+      template <typename StreamType>
+      void print_memory_consumption(StreamType     &out,
+                                    const SizeInfo &task_info) const;
+
+      /**
+       * Returns the memory consumption in bytes.
+       */
+      std::size_t memory_consumption () const;
+    };
+
+
+
+    /**
+     * The class that stores all geometry-dependent data related with cell
+     * interiors for use in the matrix-free class.
+     *
+     * @author Katharina Kormann and Martin Kronbichler, 2010, 2011, 2017
+     */
+    template <int dim, typename Number>
+    struct MappingInfo
+    {
       /**
        * Empty constructor.
        */
       MappingInfo();
 
       /**
-       * Compute the information in the given cells. The cells are specified
-       * by the level and the index within the level (as given by
+       * Compute the information in the given cells and faces. The cells are
+       * specified by the level and the index within the level (as given by
        * CellIterator::level() and CellIterator::index(), in order to allow
        * for different kinds of iterators, e.g. standard DoFHandler,
        * multigrid, etc.)  on a fixed Triangulation. In addition, a mapping
        * and several quadrature formulas are given.
        */
       void initialize (const dealii::Triangulation<dim>                &tria,
-                       const std::vector<std::pair<unsigned int,unsigned int> > &cells,
+                       const std::vector<std::pair<unsigned int,unsigned int> >  &cells,
                        const std::vector<unsigned int>         &active_fe_index,
                        const Mapping<dim>                      &mapping,
                        const std::vector<dealii::hp::QCollection<1> >  &quad,
-                       const UpdateFlags                        update_flags);
-
-      /**
-       * Helper function to determine which update flags must be set in the
-       * internal functions to initialize all data as requested by the user.
-       */
-      static UpdateFlags
-      compute_update_flags (const UpdateFlags                        update_flags,
-                            const std::vector<dealii::hp::QCollection<1> >  &quad =
-                              std::vector<dealii::hp::QCollection<1> >());
+                       const UpdateFlags                        update_flags_cells);
 
       /**
        * Return the type of a given cell as detected during initialization.
        */
       CellType get_cell_type (const unsigned int cell_chunk_no) const;
 
-      /**
-       * Return the type of a given cell as detected during initialization.
-       */
-      unsigned int get_cell_data_index (const unsigned int cell_chunk_no) const;
-
       /**
        * Clear all data fields in this class.
        */
@@ -118,232 +315,143 @@ namespace internal
        */
       template <typename StreamType>
       void print_memory_consumption(StreamType     &out,
-                                    const SizeInfo &size_info) const;
+                                    const SizeInfo &task_info) const;
 
       /**
-       * Stores whether a cell is Cartesian, has constant transform data
-       * (Jacobians) or is general. cell_type % 4 gives this information (0:
-       * Cartesian, 1: constant Jacobian throughout cell, 2: general cell),
-       * and cell_type / 4 gives the index in the data field of where to find
-       * the information in the fields Jacobian and JxW values (except for
-       * quadrature points, for which the index runs as usual).
+       * Stores whether a cell is Cartesian (cell type 0), has constant
+       * transform data (Jacobians) (cell type 1), or is general (cell type
+       * 3). Type 2 is only used for faces and no cells are assigned this
+       * value.
        */
-      std::vector<unsigned int> cell_type;
+      std::vector<CellType> cell_type;
 
       /**
-       * The first field stores the inverse Jacobian for Cartesian cells:
-       * There, it is a diagonal rank-2 tensor, so we actually just store a
-       * rank-1 tensor. It is the same on all cells, therefore we only store
-       * it once per cell, and use similarities from one cell to another, too
-       * (on structured meshes, there are usually many cells with the same
-       * Jacobian).
-       *
-       * The second field stores the Jacobian determinant for Cartesian cells
-       * (without the quadrature weight, which depends on the quadrature
-       * point, whereas the determinant is the same on each quadrature point).
+       * Stores whether a face (and both cells adjacent to the face) is
+       * Cartesian (face type 0), whether it represents an affine situation
+       * (face type 1), whether it is a flat face where the normal vector is
+       * the same throughout the face (face type 2), or is general (face type
+       * 3).
        */
-      AlignedVector<std::pair<Tensor<1,dim,VectorizedArray<Number> >,
-                    VectorizedArray<Number> > > cartesian_data;
+      std::vector<CellType> face_type;
 
       /**
-       * The first field stores the Jacobian for non-Cartesian cells where all
-       * the Jacobians on the cell are the same (i.e., constant, which comes
-       * from a linear transformation from unit to real cell). Also use
-       * similarities from one cell to another (on structured meshes, there
-       * are usually many cells with the same Jacobian).
-       *
-       * The second field stores the Jacobian determinant for non-Cartesian
-       * cells with constant Jacobian throughout the cell (without the
-       * quadrature weight, which depends on the quadrature point, whereas the
-       * determinant is the same on each quadrature point).
+       * The data cache for the cells.
        */
-      AlignedVector<std::pair<Tensor<2,dim,VectorizedArray<Number> >,
-                    VectorizedArray<Number> > > affine_data;
+      std::vector<MappingInfoStorage<dim,dim,Number> > cell_data;
 
       /**
-       * Definition of a structure that stores data that depends on the
-       * quadrature formula (if we have more than one quadrature formula on a
-       * given problem, these fields will be different)
+       * The data cache for the faces.
        */
-      struct MappingInfoDependent
-      {
-        /**
-         * This field stores the row starts for the inverse Jacobian
-         * transformations, quadrature weights and second derivatives.
-         */
-        std::vector<unsigned int> rowstart_jacobians;
-
-        /**
-         * This field stores the inverse Jacobian transformation from unit to
-         * real cell, which is needed for most gradient transformations
-         * (corresponds to FEValues::inverse_jacobian) for general cells.
-         */
-        AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > jacobians;
-
-        /**
-         * This field stores the Jacobian determinant times the quadrature
-         * weights (JxW in deal.II speak) for general cells.
-         */
-        AlignedVector<VectorizedArray<Number> > JxW_values;
-
-        /**
-         * Stores the diagonal part of the gradient of the inverse Jacobian
-         * transformation. The first index runs over the derivatives
-         * $\partial^2/\partial x_i^2$, the second over the space coordinate.
-         * Needed for computing the Laplacian of FE functions on the real
-         * cell. Uses a separate storage from the off-diagonal part
-         * $\partial^2/\partial x_i \partial x_j, i\neq j$ because that is
-         * only needed for computing a full Hessian.
-         */
-        AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > jacobians_grad_diag;
-
-        /**
-         * Stores the off-diagonal part of the gradient of the inverse
-         * Jacobian transformation. Because of symmetry, only the upper
-         * diagonal part is needed. The first index runs through the
-         * derivatives row-wise, i.e., $\partial^2/\partial x_1 \partial x_2$
-         * first, then $\partial^2/\partial x_1 \partial x_3$, and so on. The
-         * second index is the spatial coordinate. Not filled currently.
-         */
-        AlignedVector<Tensor<1,(dim>1?dim*(dim-1)/2:1),
-                      Tensor<1,dim,VectorizedArray<Number> > > > jacobians_grad_upper;
+      std::vector<MappingInfoStorage<dim-1,dim,Number> > face_data;
 
-        /**
-         * Stores the row start for quadrature points in real coordinates for
-         * both types of cells. Note that Cartesian cells will have shorter
-         * fields (length is @p n_q_points_1d) than non-Cartesian cells
-         * (length is @p n_q_points).
-         */
-        std::vector<unsigned int> rowstart_q_points;
-
-        /**
-         * Stores the quadrature points in real coordinates for Cartesian
-         * cells (does not need to store the full data on all points)
-         */
-        AlignedVector<Point<dim,VectorizedArray<Number> > > quadrature_points;
-
-        /**
-         * The dim-dimensional quadrature formula underlying the problem
-         * (constructed from a 1D tensor product quadrature formula).
-         */
-        dealii::hp::QCollection<dim>    quadrature;
+      /**
+       * The data cache for the face-associated-with-cell topology, following
+       * the @p cell_type variable for the cell types.
+       */
+      std::vector<MappingInfoStorage<dim-1,dim,Number> > face_data_by_cells;
 
-        /**
-         * The (dim-1)-dimensional quadrature formula corresponding to face
-         * evaluation (constructed from a 1D tensor product quadrature
-         * formula).
-         */
-        dealii::hp::QCollection<dim-1>  face_quadrature;
+      /**
+       * Computes the information in the given cells, called within
+       * initialize.
+       */
+      void initialize_cells (const dealii::Triangulation<dim>                &tria,
+                             const std::vector<std::pair<unsigned int,unsigned int> > &cells,
+                             const std::vector<unsigned int>         &active_fe_index,
+                             const Mapping<dim>                      &mapping,
+                             const std::vector<dealii::hp::QCollection<1> >  &quad,
+                             const UpdateFlags                        update_flags_cells);
 
-        /**
-         * The number of quadrature points for the current quadrature formula.
-         */
-        std::vector<unsigned int> n_q_points;
+      /**
+       * Helper function to determine which update flags must be set in the
+       * internal functions to initialize all data as requested by the user.
+       */
+      static UpdateFlags
+      compute_update_flags (const UpdateFlags                        update_flags,
+                            const std::vector<dealii::hp::QCollection<1> > &quad =
+                              std::vector<dealii::hp::QCollection<1> >());
+    };
 
-        /**
-         * The number of quadrature points for the current quadrature formula
-         * when applied to a face. Only set if the quadrature formula is
-         * derived from a tensor product, since it is not defined from the
-         * full quadrature formula otherwise.
-         */
-        std::vector<unsigned int> n_q_points_face;
 
-        /**
-         * The quadrature weights (vectorized data format) on the unit cell.
-         */
-        std::vector<AlignedVector<VectorizedArray<Number> > > quadrature_weights;
 
-        /**
-         * This variable stores the number of quadrature points for all
-         * quadrature indices in the underlying element for easier access to
-         * data in the hp case.
-         */
-        std::vector<unsigned int> quad_index_conversion;
+    /**
+     * A helper class to extract either cell or face data from mapping info
+     * for use in FEEvaluationBase.
+     *
+     * @author Katharina Kormann, Martin Kronbichler, 2018
+     */
+    template <int, typename, bool> struct MappingInfoCellsOrFaces;
 
-        /**
-         * Return the quadrature index for a given number of quadrature
-         * points. If not in hp mode or if the index is not found, this
-         * function always returns index 0. Hence, this function does not
-         * check whether the given degree is actually present.
-         */
-        unsigned int
-        quad_index_from_n_q_points (const unsigned int n_q_points) const;
+    template <int dim, typename Number>
+    struct MappingInfoCellsOrFaces<dim,Number,false>
+    {
+      static const MappingInfoStorage<dim,dim,Number> *
+      get(const MappingInfo<dim,Number> &mapping_info,
+          const unsigned int quad_no)
+      {
+        AssertIndexRange(quad_no, mapping_info.cell_data.size());
+        return &mapping_info.cell_data[quad_no];
+      }
+    };
 
+    template <int dim, typename Number>
+    struct MappingInfoCellsOrFaces<dim,Number,true>
+    {
+      static const MappingInfoStorage<dim-1,dim,Number> *
+      get(const MappingInfo<dim,Number> &mapping_info,
+          const unsigned int quad_no)
+      {
+        AssertIndexRange(quad_no, mapping_info.face_data.size());
+        return &mapping_info.face_data[quad_no];
+      }
+    };
 
-        /**
-         * Prints a detailed summary of memory consumption in the different
-         * structures of this class to the given output stream.
-         */
-        template <typename StreamType>
-        void print_memory_consumption(StreamType     &out,
-                                      const SizeInfo &size_info) const;
 
-        /**
-         * Return the memory consumption in bytes.
-         */
-        std::size_t memory_consumption () const;
-      };
 
-      /**
-       * Contains all the stuff that depends on the quadrature formula
-       */
-      std::vector<MappingInfoDependent> mapping_data_gen;
+    /**
+     * A class that is used to compare floating point arrays (e.g. std::vectors,
+     * Tensor<1,dim>, etc.). The idea of this class is to consider two arrays as
+     * equal if they are the same within a given tolerance. We use this
+     * comparator class within an std::map<> of the given arrays. Note that this
+     * comparison operator does not satisfy all the mathematical properties one
+     * usually wants to have (consider e.g. the numbers a=0, b=0.1, c=0.2 with
+     * tolerance 0.15; the operator gives a<c, but neither of a<b? or b<c? is
+     * satisfied). This is not a problem in the use cases for this class, but be
+     * careful when using it in other contexts.
+     */
+    template <typename Number>
+    struct FPArrayComparator
+    {
+      FPArrayComparator (const Number scaling);
 
-      /**
-       * Stores whether JxW values have been initialized
-       */
-      bool JxW_values_initialized;
+      bool operator() (const std::vector<Number> &v1,
+                       const std::vector<Number> &v2) const;
 
-      /**
-       * Stores whether we computed second derivatives.
-       */
-      bool second_derivatives_initialized;
+      bool operator ()(const Tensor<1,VectorizedArray<Number>::n_array_elements,Number> &t1,
+                       const Tensor<1,VectorizedArray<Number>::n_array_elements,Number> &t2) const;
 
-      /**
-       * Stores whether we computed quadrature points.
-       */
-      bool quadrature_points_initialized;
+      template <int dim>
+      bool operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
+                       const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const;
 
-      /**
-       * Internal temporary data used for the initialization.
-       */
-      struct CellData
-      {
-        CellData (const double jac_size);
-        void resize (const unsigned int size);
-
-        AlignedVector<Tensor<1,dim,VectorizedArray<Number> > >  quadrature_points;
-        AlignedVector<Tensor<2,dim,VectorizedArray<Number> > >  general_jac;
-        AlignedVector<Tensor<3,dim,VectorizedArray<Number> > >  general_jac_grad;
-        Tensor<2,dim,VectorizedArray<Number> > const_jac;
-        const double                           jac_size;
-      };
+      template <int dim>
+      bool operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
+                       const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const;
 
-      /**
-       * Helper function called internally during the initialize function.
-       */
-      void evaluate_on_cell (const dealii::Triangulation<dim> &tria,
-                             const std::pair<unsigned int,unsigned int> *cells,
-                             const unsigned int  cell,
-                             const unsigned int  my_q,
-                             CellType (&cell_t_prev)[n_vector_elements],
-                             CellType (&cell_t)[n_vector_elements],
-                             dealii::FEValues<dim,dim> &fe_values,
-                             CellData          &cell_data) const;
+      Number tolerance;
     };
 
 
 
     /* ------------------- inline functions ----------------------------- */
 
-    template <int dim, typename Number>
+    template <int structdim, int spacedim, typename Number>
     inline
     unsigned int
-    MappingInfo<dim,Number>::MappingInfoDependent::
-    quad_index_from_n_q_points (const unsigned int n_q_points) const
+    MappingInfoStorage<structdim,spacedim,Number>
+    ::quad_index_from_n_q_points (const unsigned int n_q_points) const
     {
-      for (unsigned int i=0; i<quad_index_conversion.size(); ++i)
-        if (n_q_points == quad_index_conversion[i])
+      for (unsigned int i=0; i<descriptor.size(); ++i)
+        if (n_q_points == descriptor[i].n_q_points)
           return i;
       return 0;
     }
@@ -356,20 +464,7 @@ namespace internal
     MappingInfo<dim,Number>::get_cell_type (const unsigned int cell_no) const
     {
       AssertIndexRange (cell_no, cell_type.size());
-      CellType enum_cell_type = (CellType)(cell_type[cell_no] % n_cell_types);
-      Assert(enum_cell_type != undefined, ExcInternalError());
-      return enum_cell_type;
-    }
-
-
-
-    template <int dim, typename Number>
-    inline
-    unsigned int
-    MappingInfo<dim,Number>::get_cell_data_index (const unsigned int cell_no) const
-    {
-      AssertIndexRange (cell_no, cell_type.size());
-      return cell_type[cell_no] >> n_cell_type_bits;
+      return cell_type[cell_no];
     }
 
   } // end of namespace MatrixFreeFunctions
index 8c79097a3f7a2ca32001ce1806b0cf2347c732b2..f08e728bafccd45156ed6b37d32bf55c0c47dfc8 100644 (file)
@@ -18,6 +18,8 @@
 
 #include <deal.II/base/utilities.h>
 #include <deal.II/base/memory_consumption.h>
+#include <deal.II/base/multithread_info.h>
+#include <deal.II/base/thread_management.h>
 #include <deal.II/fe/fe_nothing.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/mapping_q1.h>
@@ -32,14 +34,165 @@ namespace internal
 {
   namespace MatrixFreeFunctions
   {
-    // ----------------- actual MappingInfo functions -------------------------
+
+    /* ------------------------ MappingInfoStorage implementation ---------- */
+
+    template <int structdim, int spacedim, typename Number>
+    MappingInfoStorage<structdim, spacedim, Number>::QuadratureDescriptor
+    ::QuadratureDescriptor()
+      :
+      n_q_points (numbers::invalid_unsigned_int)
+    {
+    }
+
+
+
+    template <int structdim, int spacedim, typename Number>
+    void
+    MappingInfoStorage<structdim, spacedim, Number>::QuadratureDescriptor
+    ::initialize(const Quadrature<1> &quadrature_1d,
+                 const UpdateFlags    update_flags_inner_faces)
+    {
+      Assert(structdim+1 <= spacedim ||
+             update_flags_inner_faces == update_default,
+             ExcMessage("Volume cells do not allow for setting inner faces"));
+      quadrature = Quadrature<structdim>(quadrature_1d);
+      n_q_points = quadrature.size();
+      quadrature_weights.resize(n_q_points);
+      for (unsigned int i=0; i<n_q_points; ++i)
+        quadrature_weights[i] = quadrature.weight(i);
+
+      for (unsigned int d=0; d<structdim; ++d)
+        {
+          tensor_quadrature_weights[d].resize(quadrature_1d.size());
+          for (unsigned int i=0; i<quadrature_1d.size(); ++i)
+            tensor_quadrature_weights[d][i] = quadrature_1d.weight(i);
+        }
+
+      // face orientation for faces in 3D
+      if (structdim == spacedim-1 && spacedim == 3 &&
+          update_flags_inner_faces != update_default)
+        {
+          const unsigned int n=quadrature_1d.size();
+          face_orientations.reinit(8, n*n);
+          for (unsigned int j=0, i=0; j<n; ++j)
+            for (unsigned int k=0; k<n; ++k, ++i)
+              {
+                // face_orientation=true,  face_flip=false, face_rotation=false
+                face_orientations[0][i] = i;
+                // face_orientation=false, face_flip=false, face_rotation=false
+                face_orientations[1][i] = j       + k      *n;
+                // face_orientation=true,  face_flip=true,  face_rotation=false
+                face_orientations[2][i] = (n-1-k) + (n-1-j)*n;
+                // face_orientation=false, face_flip=true,  face_rotation=false
+                face_orientations[3][i] = (n-1-j) + (n-1-k)*n;
+                // face_orientation=true,  face_flip=false, face_rotation=true
+                face_orientations[4][i] = j       + (n-1-k)*n;
+                // face_orientation=false, face_flip=false, face_rotation=true
+                face_orientations[5][i] = k       + (n-1-j)*n;
+                // face_orientation=true,  face_flip=true,  face_rotation=true
+                face_orientations[6][i] = (n-1-j) + k      *n;
+                // face_orientation=false, face_flip=true,  face_rotation=true
+                face_orientations[7][i] = (n-1-k) + j      *n;
+              }
+        }
+    }
+
+
+
+    template <int structdim, int spacedim, typename Number>
+    std::size_t
+    MappingInfoStorage<structdim, spacedim, Number>::QuadratureDescriptor
+    ::memory_consumption() const
+    {
+      std::size_t memory =
+        sizeof (this) +
+        quadrature.memory_consumption() +
+        quadrature_weights.memory_consumption() +
+        face_orientations.memory_consumption();
+      for (unsigned int d=0; d<structdim; ++d)
+        memory += tensor_quadrature_weights[d].memory_consumption();
+      return memory;
+    }
+
+
+
+    template <int structdim, int spacedim, typename Number>
+    std::size_t
+    MappingInfoStorage<structdim,spacedim,Number>::memory_consumption() const
+    {
+      return
+        MemoryConsumption::memory_consumption (descriptor) +
+        MemoryConsumption::memory_consumption (data_index_offsets) +
+        MemoryConsumption::memory_consumption (JxW_values) +
+        MemoryConsumption::memory_consumption (normal_vectors) +
+        MemoryConsumption::memory_consumption (jacobians[0]) +
+        MemoryConsumption::memory_consumption (jacobians[1]) +
+        MemoryConsumption::memory_consumption (jacobian_gradients[0]) +
+        MemoryConsumption::memory_consumption (jacobian_gradients[1]) +
+        MemoryConsumption::memory_consumption (normals_times_jacobians[0]) +
+        MemoryConsumption::memory_consumption (normals_times_jacobians[1]) +
+        MemoryConsumption::memory_consumption (quadrature_point_offsets) +
+        MemoryConsumption::memory_consumption (quadrature_points);
+    }
+
+
+
+    template <int structdim, int spacedim, typename Number>
+    template <typename StreamType>
+    void
+    MappingInfoStorage<structdim,spacedim,Number>::print_memory_consumption
+    (StreamType     &out,
+     const SizeInfo &task_info) const
+    {
+      // print_memory_statistics involves global communication, so we can
+      // disable the check here only if no processor has any such data
+      const std::size_t size = Utilities::MPI::sum(jacobians[0].size(),
+                                                   task_info.communicator);
+      if (size > 0)
+        {
+          out << "      Memory JxW data:               ";
+          task_info.print_memory_statistics
+          (out, MemoryConsumption::memory_consumption (data_index_offsets) +
+           MemoryConsumption::memory_consumption (JxW_values));
+          out << "      Memory Jacobian data:          ";
+          task_info.print_memory_statistics
+          (out, MemoryConsumption::memory_consumption (jacobians[0]) +
+           MemoryConsumption::memory_consumption (jacobians[1]));
+          out << "      Memory second derivative data: ";
+          task_info.print_memory_statistics
+          (out, MemoryConsumption::memory_consumption (jacobian_gradients[0]) +
+           MemoryConsumption::memory_consumption (jacobian_gradients[1]));
+        }
+      const std::size_t normal_size = Utilities::MPI::sum(normal_vectors.size(),
+                                                          task_info.communicator);
+      if (normal_size > 0)
+        {
+          out << "      Memory normal vectors data:    ";
+          task_info.print_memory_statistics
+          (out, MemoryConsumption::memory_consumption (normal_vectors) +
+           MemoryConsumption::memory_consumption (normals_times_jacobians[0]) +
+           MemoryConsumption::memory_consumption (normals_times_jacobians[1]));
+        }
+
+      const std::size_t quad_size =
+        Utilities::MPI::sum(quadrature_points.size(),
+                            task_info.communicator);
+      if (quad_size > 0)
+        {
+          out << "      Memory quadrature points:      ";
+          task_info.print_memory_statistics
+          (out, MemoryConsumption::memory_consumption (quadrature_point_offsets) +
+           MemoryConsumption::memory_consumption (quadrature_points));
+        }
+    }
+
+
+
+    /* ------------------------ MappingInfo implementation ----------------- */
 
     template <int dim, typename Number>
     MappingInfo<dim,Number>::MappingInfo()
-      :
-      JxW_values_initialized (false),
-      second_derivatives_initialized (false),
-      quadrature_points_initialized (false)
     {}
 
 
@@ -48,13 +201,11 @@ namespace internal
     void
     MappingInfo<dim,Number>::clear ()
     {
-      JxW_values_initialized = false;
-      quadrature_points_initialized = false;
-      second_derivatives_initialized = false;
-      mapping_data_gen.clear();
+      cell_data.clear();
+      face_data.clear();
+      face_data_by_cells.clear();
       cell_type.clear();
-      cartesian_data.clear();
-      affine_data.clear();
+      face_type.clear();
     }
 
 
@@ -65,17 +216,10 @@ namespace internal
     compute_update_flags (const UpdateFlags update_flags,
                           const std::vector<dealii::hp::QCollection<1> > &quad)
     {
-      // this class is build around the evaluation this class is build around
-      // the evaluation of inverse gradients, so compute them in any case
-      UpdateFlags new_flags = update_inverse_jacobians;
-
-      // if the user requested gradients, need inverse Jacobians
-      if (update_flags & update_gradients || update_flags & update_inverse_jacobians)
-        new_flags |= update_inverse_jacobians;
-
-      // for JxW, would only need JxW values.
-      if (update_flags & update_JxW_values)
-        new_flags |= update_JxW_values;
+      // this class is build around the evaluation of jacobians, so compute
+      // them in any case. The Jacobians will be inverted manually. Since we
+      // always do support integration, we also include the JxW values
+      UpdateFlags new_flags = update_jacobians | update_JxW_values;
 
       // for Hessian information, need inverse Jacobians and the derivative of
       // Jacobians (these two together will give use the gradients of the
@@ -112,7 +256,30 @@ namespace internal
 
 
 
-    namespace internal
+    template <int dim, typename Number>
+    void
+    MappingInfo<dim,Number>::initialize
+    (const dealii::Triangulation<dim>                          &tria,
+     const std::vector<std::pair<unsigned int,unsigned int> >  &cells,
+     const std::vector<unsigned int>                           &active_fe_index,
+     const Mapping<dim>                                        &mapping,
+     const std::vector<dealii::hp::QCollection<1> >            &quad,
+     const UpdateFlags                                          update_flags_cells)
+    {
+      clear();
+
+      // Could call these functions in parallel, but not useful because the
+      // work inside is nicely split up already
+      initialize_cells(tria, cells, active_fe_index, mapping, quad, update_flags_cells);
+    }
+
+
+
+    /* ------------------------- initialization of cells ------------------- */
+
+    // Anonymous namespace with implementation of extraction of values on cell
+    // range
+    namespace
     {
       template <int dim>
       double get_jacobian_size (const dealii::Triangulation<dim> &tria)
@@ -121,331 +288,438 @@ namespace internal
           return 1;
         else return tria.begin()->diameter();
       }
-    }
 
 
 
-    template <int dim, typename Number>
-    void
-    MappingInfo<dim,Number>::initialize
-    (const dealii::Triangulation<dim>                         &tria,
-     const std::vector<std::pair<unsigned int,unsigned int> > &cells,
-     const std::vector<unsigned int>                          &active_fe_index,
-     const Mapping<dim>                                       &mapping,
-     const std::vector<dealii::hp::QCollection<1> >           &quad,
-     const UpdateFlags                                         update_flags_input)
-    {
-      clear();
-      const unsigned int n_quads = quad.size();
-      const unsigned int n_cells = cells.size();
-      const unsigned int vectorization_length =
-        VectorizedArray<Number>::n_array_elements;
-      Assert (n_cells%vectorization_length == 0, ExcInternalError());
-      const unsigned int n_macro_cells = n_cells/vectorization_length;
-      mapping_data_gen.resize (n_quads);
-      cell_type.resize (n_macro_cells);
+      template <int dim, typename Number>
+      struct CompressedCellData
+      {
+        CompressedCellData(const double expected_size)
+          :
+          data(FPArrayComparator<Number>(expected_size))
+        {}
+
+        std::map<Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> >
+        , unsigned int, FPArrayComparator<Number> > data;
+      };
+
+      /**
+       * Internal temporary data used for the initialization.
+       */
+      template <int dim, typename Number>
+      struct LocalData
+      {
+        LocalData (const double jac_size);
+        void resize (const unsigned int size);
 
-      // dummy FE that is used to set up an FEValues object. Do not need the
-      // actual finite element because we will only evaluate quantities for
-      // the mapping that are independent of the FE
-      FE_Nothing<dim> dummy_fe;
-      UpdateFlags update_flags = compute_update_flags (update_flags_input, quad);
+        AlignedVector<Point<dim,VectorizedArray<Number> > > quadrature_points;
+        AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > general_jac;
+        AlignedVector<VectorizedArray<Number> > JxW_values;
+        AlignedVector<Tensor<3,dim,VectorizedArray<Number> > > general_jac_grad;
+        AlignedVector<Tensor<1,dim,VectorizedArray<Number> > > normal_vectors;
+        Tensor<2,dim,VectorizedArray<Number> > const_jac;
+        const double                           jac_size;
+      };
 
-      if (update_flags & update_JxW_values)
-        JxW_values_initialized = true;
-      if (update_flags & update_jacobian_grads)
-        second_derivatives_initialized = true;
-      if (update_flags & update_quadrature_points)
-        quadrature_points_initialized = true;
-
-      // when we make comparisons about the size of Jacobians we need to know
-      // the approximate size of typical entries in Jacobians. We need to fix
-      // the Jacobian size once and for all. We choose the diameter of the
-      // first cell (on level zero, which is the best accuracy we can hope
-      // for, since diameters on finer levels are computed by differences of
-      // nearby cells). If the mesh extends over a certain domain, the
-      // precision of double values is essentially limited by this precision.
-      const double jacobian_size = internal::get_jacobian_size(tria);
-
-      // objects that hold the data for up to vectorization_length cells while
-      // we fill them up. Only after all vectorization_length cells have been
-      // processed, we can insert the data into the data structures of this
-      // class
-      CellData data (jacobian_size);
 
-      for (unsigned int my_q=0; my_q<n_quads; ++my_q)
-        {
-          MappingInfoDependent &current_data = mapping_data_gen[my_q];
-          const unsigned int n_hp_quads = quad[my_q].size();
-          AssertIndexRange (0, n_hp_quads);
-          current_data.n_q_points.reserve (n_hp_quads);
-          current_data.n_q_points_face.reserve (n_hp_quads);
-          current_data.quadrature_weights.resize (n_hp_quads);
-          std::vector<unsigned int> n_q_points_1d (n_hp_quads),
-              step_size_cartesian (n_hp_quads);
-          if (n_hp_quads > 1)
-            current_data.quad_index_conversion.resize(n_hp_quads);
-          for (unsigned int q=0; q<n_hp_quads; ++q)
-            {
-              n_q_points_1d[q] = quad[my_q][q].size();
-              const unsigned int n_q_points =
-                Utilities::fixed_power<dim>(n_q_points_1d[q]);
-              current_data.n_q_points.push_back (n_q_points);
-
-              current_data.n_q_points_face.push_back
-              (dim>1 ? Utilities::fixed_power<dim-1>(n_q_points_1d[q]) : 1);
-              current_data.quadrature.push_back
-              (Quadrature<dim>(quad[my_q][q]));
-              current_data.face_quadrature.push_back
-              (Quadrature<dim-1>(quad[my_q][q]));
-
-              // set quadrature weights in vectorized form
-              current_data.quadrature_weights[q].resize(n_q_points);
-              for (unsigned int i=0; i<n_q_points; ++i)
-                current_data.quadrature_weights[q][i] =
-                  current_data.quadrature[q].get_weights()[i];
-
-              if (n_hp_quads > 1)
-                current_data.quad_index_conversion[q] = n_q_points;
-
-              // To walk on the diagonal for lexicographic ordering, we have
-              // to jump one index ahead in each direction. For direction 0,
-              // this is just the next point, for direction 1, it means adding
-              // n_q_points_1d, and so on.
-              step_size_cartesian[q] = 0;
-              unsigned int factor = 1;
-              for (unsigned int d=0; d<dim; ++d)
+
+      template <int dim, typename Number>
+      LocalData<dim,Number>::LocalData (const double jac_size_in)
+        :
+        jac_size (jac_size_in)
+      {}
+
+
+
+      template <int dim, typename Number>
+      void
+      LocalData<dim,Number>::resize (const unsigned int size)
+      {
+        if (JxW_values.size() != size)
+          {
+            quadrature_points.resize_fast(size);
+            general_jac.resize_fast(size*2);
+            JxW_values.resize_fast(size);
+            general_jac_grad.resize_fast(size*2);
+            normal_vectors.resize_fast(size);
+          }
+      }
+
+      /**
+       * Helper function called internally during the initialize function.
+       */
+      template <int dim, typename Number>
+      void evaluate_on_cell (const dealii::Triangulation<dim> &tria,
+                             const std::pair<unsigned int,unsigned int> *cells,
+                             const unsigned int         my_q,
+                             CellType                  &cell_t_prev,
+                             CellType (&cell_t)[VectorizedArray<Number>::n_array_elements],
+                             dealii::FEValues<dim,dim> &fe_val,
+                             LocalData<dim,Number>     &cell_data)
+      {
+        const unsigned int n_q_points = fe_val.n_quadrature_points;
+        const UpdateFlags update_flags = fe_val.get_update_flags();
+
+        cell_data.const_jac = Tensor<2,dim,VectorizedArray<Number> >();
+
+        // this should be the same value as used in HashValue::scaling (but we
+        // not have that field here)
+        const double zero_tolerance_double = cell_data.jac_size *
+                                             std::numeric_limits<double>::epsilon() * 1024.;
+        for (unsigned int j=0; j<VectorizedArray<Number>::n_array_elements; ++j)
+          {
+            typename dealii::Triangulation<dim>::cell_iterator
+            cell_it (&tria, cells[j].first, cells[j].second);
+            fe_val.reinit(cell_it);
+            cell_t[j] = general;
+
+            // extract quadrature points and store them temporarily. if we have
+            // Cartesian cells, we can compress the indices
+            if (update_flags & update_quadrature_points)
+              for (unsigned int q=0; q<n_q_points; ++q)
                 {
-                  step_size_cartesian[q] += factor;
-                  factor *= n_q_points_1d[q];
+                  const Point<dim> &point = fe_val.quadrature_point(q);
+                  for (unsigned int d=0; d<dim; ++d)
+                    cell_data.quadrature_points[q][d][j] = point[d];
                 }
-            }
 
-          // if there are no cells, there is nothing to do
-          if (cells.size() == 0)
-            continue;
-
-          Tensor<3,dim,VectorizedArray<Number> > jac_grad, grad_jac_inv;
-          Tensor<1,dim,VectorizedArray<Number> > tmp;
-
-          // encodes the cell types of the current cell. Since several cells
-          // must be considered together, this variable holds the individual
-          // info of the last chunk of cells
-          CellType cell_t [vectorization_length],
-                   cell_t_prev [vectorization_length];
-          for (unsigned int j=0; j<vectorization_length; ++j)
-            cell_t_prev[j] = undefined;
-
-          // fe_values object that is used to compute the mapping data. for
-          // the hp case there might be more than one finite element. since we
-          // manually select the active FE index and not via a
-          // hp::DoFHandler<dim>::active_cell_iterator, we need to manually
-          // select the correct finite element, so just hold a vector of
-          // FEValues
-          std::vector<std::shared_ptr<dealii::FEValues<dim> > >
-          fe_values (current_data.quadrature.size());
-          UpdateFlags update_flags_feval =
-            ((update_flags & update_inverse_jacobians) ? update_jacobians : update_default) |
-            ((update_flags & update_jacobian_grads) ? update_jacobian_grads : update_default) |
-            ((update_flags & update_quadrature_points) ? update_quadrature_points : update_default);
-
-          // resize the fields that have fixed size or for which we know
-          // something from an earlier loop
-          current_data.rowstart_q_points.resize (n_macro_cells+1);
-          if (my_q > 0)
-            {
-              const unsigned int n_cells_var =
-                mapping_data_gen[0].rowstart_jacobians.size()-1;
-              current_data.rowstart_jacobians.reserve (n_cells_var+1);
-              const unsigned int reserve_size = n_cells_var *
-                                                current_data.n_q_points[0];
-              if (mapping_data_gen[0].jacobians.size() > 0)
-                current_data.jacobians.reserve (reserve_size);
-              if (mapping_data_gen[0].JxW_values.size() > 0)
-                current_data.jacobians.reserve (reserve_size);
-              if (mapping_data_gen[0].jacobians_grad_diag.size() > 0)
-                current_data.jacobians_grad_diag.reserve (reserve_size);
-              if (mapping_data_gen[0].jacobians_grad_upper.size() > 0)
-                current_data.jacobians_grad_upper.reserve (reserve_size);
-            }
+            // if this is not the first quadrature formula and we already have
+            // determined that this cell is either Cartesian or with constant
+            // Jacobian, we have nothing more to do.
+            if (my_q > 0 && cell_t_prev <= affine)
+              continue;
+
+            // first round: if the transformation is detected to be the same as
+            // on the old cell, we only need to copy over the data.
+            if (fe_val.get_cell_similarity() == CellSimilarity::translation
+                &&
+                my_q == 0)
+              {
+                if (j==0)
+                  cell_t[j] = cell_t_prev;
+                else
+                  cell_t[j] = cell_t[j-1];
+              }
+
+            const DerivativeForm<1,dim,dim> &jac_0 = fe_val.jacobian(0);
+
+            if (my_q == 0)
+              {
+                // check whether the Jacobian is constant on this cell the first
+                // time we come around here
+                if (cell_t[j] == general)
+                  {
+                    bool jacobian_constant = true;
+                    for (unsigned int q=1; q<n_q_points; ++q)
+                      {
+                        const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(q);
+                        for (unsigned int d=0; d<dim; ++d)
+                          for (unsigned int e=0; e<dim; ++e)
+                            if (std::fabs(jac_0[d][e]-jac[d][e]) >
+                                zero_tolerance_double)
+                              jacobian_constant = false;
+                        if (jacobian_constant == false)
+                          break;
+                      }
+
+                    // check whether the Jacobian is diagonal to machine
+                    // accuracy
+                    bool cell_cartesian = jacobian_constant;
+                    for (unsigned int d=0; d<dim; ++d)
+                      for (unsigned int e=0; e<dim; ++e)
+                        if (d!=e)
+                          if (std::fabs(jac_0[d][e]) >
+                              zero_tolerance_double)
+                            {
+                              cell_cartesian=false;
+                              break;
+                            }
+
+                    // in case we have only one quadrature point, we can have
+                    // non-constant Jacobians, but we cannot detect it by
+                    // comparison from one quadrature point to the next: in that
+                    // case, need to look at second derivatives and see whether
+                    // there are some non-zero entries (this is necessary since
+                    // we determine the constness of the Jacobian for the first
+                    // quadrature formula and might not look at them any more
+                    // for the second, third quadrature formula). in any case,
+                    // the flag update_jacobian_grads will be set in that case
+                    if (cell_cartesian == false && n_q_points == 1 &&
+                        update_flags & update_jacobian_grads)
+                      {
+                        const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(0);
+                        const DerivativeForm<2,dim,dim> &jacobian_grad =
+                          fe_val.jacobian_grad(0);
+                        for (unsigned int d=0; d<dim; ++d)
+                          for (unsigned int e=0; e<dim; ++e)
+                            for (unsigned int f=0; f<dim; ++f)
+                              {
+                                double jac_grad_comp = (jac[f][0] *
+                                                        jacobian_grad[d][e][0]);
+                                for (unsigned int g=1; g<dim; ++g)
+                                  jac_grad_comp += (jac[f][g] *
+                                                    jacobian_grad[d][e][g]);
+                                if (std::fabs(jac_grad_comp) >
+                                    zero_tolerance_double)
+                                  jacobian_constant = false;
+                              }
+                      }
+                    // set cell type
+                    if (cell_cartesian == true)
+                      cell_t[j] = cartesian;
+                    else if (jacobian_constant == true)
+                      cell_t[j] = affine;
+                    else
+                      cell_t[j] = general;
+                  }
+
+                // Cartesian cell
+                if (cell_t[j] == cartesian)
+                  {
+                    // set Jacobian into diagonal (off-diagonal part is already
+                    // zeroed out)
+                    for (unsigned int d=0; d<dim; ++d)
+                      cell_data.const_jac[d][d][j] = jac_0[d][d];
+                    continue;
+                  }
 
-          // we would like to put a Tensor<1,dim,VectorizedArray<Number> > as
-          // key into the std::map, but std::map allocation does not align the
-          // allocated memory correctly, so put it into a tensor of the
-          // correct length instead
-          FPArrayComparator<Number> comparator(jacobian_size);
-          typedef Tensor<1,VectorizedArray<Number>::n_array_elements,Number> VEC_ARRAY;
-          std::map<Tensor<1,dim,VEC_ARRAY>, unsigned int,
-              FPArrayComparator<Number> > cartesians(comparator);
-          std::map<Tensor<2,dim,VEC_ARRAY>, unsigned int,
-              FPArrayComparator<Number> > affines(comparator);
-
-          // loop over all cells
-          for (unsigned int cell=0; cell<n_macro_cells; ++cell)
+                // cell with affine mapping
+                else if (cell_t[j] == affine)
+                  {
+                    // compress out very small values
+                    for (unsigned int d=0; d<dim; ++d)
+                      for (unsigned int e=0; e<dim; ++e)
+                        if (std::fabs(jac_0[d][e]))
+                          cell_data.const_jac[d][e][j] = jac_0[d][e];
+                    continue;
+                  }
+              }
+
+            // general cell case
+
+            // go through all quadrature points and fill in the data into the
+            // temporary data structures with slots for the vectorized data
+            // types
+            for (unsigned int q=0; q<n_q_points; ++q)
+              {
+                // compress out very small numbers which are only noise. Then it
+                // is cleaner to use zero straight away (though it does not save
+                // any memory)
+                const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(q);
+                for (unsigned int d=0; d<dim; ++d)
+                  for (unsigned int e=0; e<dim; ++e)
+                    cell_data.general_jac[q][d][e][j] =
+                      std::fabs(jac[d][e]) < zero_tolerance_double ? 0. : jac[d][e];
+
+                // need to do some calculus based on the gradient of the
+                // Jacobian, in order to find the gradient of the inverse
+                // Jacobian which is needed in user code. however, we would like
+                // to perform that on vectorized data types instead of doubles
+                // or floats. to this end, copy the gradients first
+                if (update_flags & update_jacobian_grads)
+                  {
+                    const DerivativeForm<2,dim,dim> &jacobian_grad = fe_val.jacobian_grad(q);
+                    for (unsigned int d=0; d<dim; ++d)
+                      for (unsigned int e=0; e<dim; ++e)
+                        for (unsigned int f=0; f<dim; ++f)
+                          cell_data.general_jac_grad[q][d][e][f][j] = jacobian_grad[d][e][f];
+                  }
+              }
+          } // end loop over entries of vectorization (n_array_elements cells)
+
+        // set information for next cell
+        cell_t_prev = cell_t[VectorizedArray<Number>::n_array_elements-1];
+      }
+
+
+
+      template <int dim, typename Number>
+      void
+      initialize_cell_range
+      (const std::pair<unsigned int,unsigned int>                cell_range,
+       const dealii::Triangulation<dim>                         &tria,
+       const std::vector<std::pair<unsigned int,unsigned int> > &cells,
+       const std::vector<unsigned int>                          &active_fe_index,
+       const Mapping<dim>                                       &mapping,
+       const std::vector<dealii::hp::QCollection<1> >           &quad,
+       const UpdateFlags                                         update_flags,
+       MappingInfo<dim,Number>                                  &mapping_info,
+       std::pair<std::vector<MappingInfoStorage<dim,dim,Number> >,
+       CompressedCellData<dim,Number> >                         &data)
+      {
+        FE_Nothing<dim> dummy_fe;
+
+        Tensor<3,dim,VectorizedArray<Number> > jac_grad, grad_jac_inv;
+        Tensor<1,dim,VectorizedArray<Number> > tmp;
+
+        // when we make comparisons about the size of Jacobians we need to
+        // know the approximate size of typical entries in Jacobians. We need
+        // to fix the Jacobian size once and for all. We choose the diameter
+        // of the first cell (on level zero, which is the best accuracy we can
+        // hope for, since diameters on finer levels are computed by
+        // differences of nearby cells) as the order of magnitude by which we
+        // make comparisons "relative."
+        const double jacobian_size = get_jacobian_size(tria);
+
+        // objects that hold the data for up to vectorization_width cells while
+        // we fill them up. Only after all vectorization_width cells have been
+        // processed, we can insert the data into the data structures of this
+        // class
+        LocalData<dim,Number> cell_data (jacobian_size);
+
+        // encodes the cell types of the current cell. Since several cells
+        // must be considered together, this variable holds the individual
+        // info of the last chunk of cells
+        CellType cell_t [VectorizedArray<Number>::n_array_elements];
+        CellType cell_t_prev = general;
+
+        // fe_values object that is used to compute the mapping data. for
+        // the hp case there might be more than one finite element. since we
+        // manually select the active FE index and not via a
+        // hp::DoFHandler<dim>::active_cell_iterator, we need to manually
+        // select the correct finite element, so just hold a vector of
+        // FEValues
+        std::vector<std::vector<std::shared_ptr<dealii::FEValues<dim> > > >
+        fe_values (mapping_info.cell_data.size());
+        for (unsigned int i=0; i<fe_values.size(); ++i)
+          fe_values[i].resize(mapping_info.cell_data[i].descriptor.size());
+        UpdateFlags update_flags_feval =
+          (update_flags & update_jacobians ? update_jacobians : update_default) |
+          (update_flags & update_jacobian_grads ? update_jacobian_grads : update_default) |
+          (update_flags & update_quadrature_points ? update_quadrature_points : update_default);
+
+        std::vector<std::vector<unsigned int> > n_q_points_1d (quad.size()),
+            step_size_cartesian (quad.size());
+        for (unsigned int my_q=0; my_q<quad.size(); ++my_q)
+          {
+            n_q_points_1d[my_q].resize(quad[my_q].size());
+            step_size_cartesian[my_q].resize(quad[my_q].size());
+            for (unsigned int hpq=0; hpq<quad[my_q].size(); ++hpq)
+              {
+                n_q_points_1d[my_q][hpq] = quad[my_q][hpq].size();
+
+                // To walk on the diagonal for lexicographic ordering, we have
+                // to jump one index ahead in each direction. For direction 0,
+                // this is just the next point, for direction 1, it means adding
+                // n_q_points_1d, and so on.
+                step_size_cartesian[my_q][hpq] = 0;
+                unsigned int factor = 1;
+                for (unsigned int d=0; d<dim; ++d)
+                  {
+                    step_size_cartesian[my_q][hpq] += factor;
+                    factor *= n_q_points_1d[my_q][hpq];
+                  }
+              }
+          }
+
+        const unsigned int end_cell = std::min(mapping_info.cell_type.size(),
+                                               std::size_t(cell_range.second));
+        // loop over given cells
+        for (unsigned int cell=cell_range.first; cell<end_cell; ++cell)
+          for (unsigned int my_q=0; my_q<mapping_info.cell_data.size(); ++my_q)
             {
               // GENERAL OUTLINE: First generate the data in format "number"
-              // for vectorization_length cells, and then find the most
+              // for vectorization_width cells, and then find the most
               // general type of cell for appropriate vectorized formats. then
               // fill this data in
               const unsigned int fe_index = active_fe_index.size() > 0 ?
                                             active_fe_index[cell] : 0;
-              const unsigned int n_q_points = current_data.n_q_points[fe_index];
-              if (fe_values[fe_index].get() == nullptr)
-                fe_values[fe_index].reset
-                (new dealii::FEValues<dim> (mapping, dummy_fe,
-                                            current_data.quadrature[fe_index],
+              const unsigned int n_q_points =
+                mapping_info.cell_data[my_q].descriptor[fe_index].n_q_points;
+              if (fe_values[my_q][fe_index].get() == nullptr)
+                fe_values[my_q][fe_index].reset
+                (new dealii::FEValues<dim> (mapping, dummy_fe, mapping_info.cell_data[my_q].
+                                            descriptor[fe_index].quadrature,
                                             update_flags_feval));
-              dealii::FEValues<dim> &fe_val = *fe_values[fe_index];
-              data.resize (n_q_points);
+              dealii::FEValues<dim> &fe_val = *fe_values[my_q][fe_index];
+              cell_data.resize (n_q_points);
 
               // if the fe index has changed from the previous cell, set the
               // old cell type to invalid (otherwise, we might detect
               // similarity due to some cells further ahead)
-              if (cell > 0 && active_fe_index.size() > 0 &&
-                  active_fe_index[cell] != active_fe_index[cell-1])
-                cell_t_prev[vectorization_length-1] = undefined;
-              evaluate_on_cell (tria, &cells[cell*vectorization_length],
-                                cell, my_q, cell_t_prev, cell_t, fe_val, data);
+              if (my_q > 0)
+                cell_t_prev = CellType(mapping_info.cell_type[cell]);
+              else if (cell > cell_range.first && active_fe_index.size() > 0 &&
+                       active_fe_index[cell] != active_fe_index[cell-1])
+                cell_t_prev = general;
+
+              evaluate_on_cell (tria, &cells[cell*VectorizedArray<Number>::n_array_elements],
+                                my_q, cell_t_prev, cell_t, fe_val,
+                                cell_data);
 
               // now reorder the data into vectorized types. if we are here
               // for the first time, we need to find out whether the Jacobian
               // allows for some simplification (Cartesian, affine) taking
-              // vectorization_length cell together and we have to insert that
-              // data into the respective fields. Also, we have to compress
-              // different cell indicators into one structure.
+              // vectorization_width cell together
 
               if (my_q == 0)
                 {
-                  // find the most general cell type (most general type is 2
+                  // find the most general cell type (most general type is 3
                   // (general cell))
                   CellType most_general_type = cartesian;
-                  for (unsigned int j=0; j<vectorization_length; ++j)
+                  for (unsigned int j=0; j<VectorizedArray<Number>::n_array_elements; ++j)
                     if (cell_t[j] > most_general_type)
                       most_general_type = cell_t[j];
-                  AssertIndexRange (most_general_type, 3);
-                  unsigned int insert_position = numbers::invalid_unsigned_int;
-
-                  // Cartesian cell with diagonal Jacobian: only insert the
-                  // diagonal of the inverse and the Jacobian determinant. We
-                  // do this by using an std::map that collects pointers to
-                  // all Cartesian Jacobians. We need a pointer in the
-                  // std::map because it cannot store data based on
-                  // VectorizedArray (alignment issue). We circumvent the
-                  // problem by temporarily filling the next value into the
-                  // cartesian_data field and, in case we did an insertion,
-                  // the data is already in the correct place.
-                  if (most_general_type == cartesian)
-                    {
-                      std::pair<Tensor<1,dim,VEC_ARRAY>,unsigned int> new_entry;
-                      new_entry.second = cartesians.size();
-                      for (unsigned int d=0; d<dim; ++d)
-                        for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
-                          new_entry.first[d][v] = data.const_jac[d][d][v];
-
-                      std::pair<typename std::map<Tensor<1,dim,VEC_ARRAY>,
-                          unsigned int, FPArrayComparator<Number> >::iterator,
-                          bool> it = cartesians.insert(new_entry);
-                      if (it.second == false)
-                        insert_position = it.first->second;
-                      else
-                        insert_position = new_entry.second;
-                    }
+                  AssertIndexRange ((unsigned int)most_general_type, 4U);
+                  mapping_info.cell_type[cell] = most_general_type;
+                }
+
+              AssertThrow(data.first[my_q].JxW_values.size() <
+                          static_cast<std::size_t>(std::numeric_limits<unsigned int>::max()),
+                          ExcMessage("Index overflow. Cannot fit data in 32 bit integers"));
 
-                  // Constant Jacobian case. same strategy as before, but with
-                  // other data fields
-                  else if (most_general_type == affine)
+              unsigned int insert_position = data.first[my_q].JxW_values.size();
+              // Cartesian/affine cell with constant Jacobians throughout the
+              // cell. We need to store the data in another data field because
+              // std::map cannot store data based on VectorizedArray directly
+              // (alignment issue).
+              if (mapping_info.cell_type[cell] <= affine)
+                {
+                  if (my_q == 0)
                     {
-                      std::pair<Tensor<2,dim,VEC_ARRAY>,unsigned int> new_entry;
-                      new_entry.second = affines.size();
+                      std::pair<Tensor<2,dim,Tensor<1,VectorizedArray<Number>::
+                      n_array_elements,Number> >,unsigned int> new_entry;
+                      // This number overlaps with the general data but we
+                      // take care of that when we merge data from different
+                      // threads
+                      new_entry.second = data.second.data.size();
                       for (unsigned int d=0; d<dim; ++d)
                         for (unsigned int e=0; e<dim; ++e)
                           for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
-                            new_entry.first[d][e][v] = data.const_jac[d][e][v];
-
-                      std::pair<typename std::map<Tensor<2,dim,VEC_ARRAY>,
-                          unsigned int, FPArrayComparator<Number> >::iterator,
-                          bool> it = affines.insert(new_entry);
-                      if (it.second == false)
-                        insert_position = it.first->second;
-                      else
-                        insert_position = new_entry.second;
-                    }
+                            new_entry.first[d][e][v] = cell_data.const_jac[d][e][v];
 
-                  // general cell case: first resize the data field to fit the
-                  // new data. if we are here the first time, assume that
-                  // there are many general cells to come, so reserve some
-                  // memory in order to not have too many reallocations and
-                  // memcpy's. The scheme used here involves at most one
-                  // reallocation.
-                  else
-                    {
-                      Assert (most_general_type == general, ExcInternalError());
-                      insert_position = current_data.rowstart_jacobians.size();
-                      if (current_data.rowstart_jacobians.size() == 0)
-                        {
-                          unsigned int reserve_size = (n_macro_cells-cell+1)/2;
-                          current_data.rowstart_jacobians.reserve
-                          (reserve_size);
-                          reserve_size *= n_q_points;
-                          current_data.jacobians.reserve (reserve_size);
-                          if (update_flags & update_JxW_values)
-                            current_data.JxW_values.reserve (reserve_size);
-                          if (update_flags & update_jacobian_grads)
-                            {
-                              current_data.jacobians_grad_diag.reserve (reserve_size);
-                              current_data.jacobians_grad_upper.reserve (reserve_size);
-                            }
-                        }
+                      insert_position = data.second.data.insert(new_entry).first->second;
                     }
-
-                  cell_type[cell] = ((insert_position << n_cell_type_bits) +
-                                     (unsigned int)most_general_type);
-
-                } // end if (my_q == 0)
+                  else
+                    insert_position = data.first[0].data_index_offsets[cell-cell_range.first];
+                }
 
               // general cell case: now go through all quadrature points and
               // collect the data. done for all different quadrature formulas,
               // so do it outside the above loop.
-              if (get_cell_type(cell) == general)
+              data.first[my_q].data_index_offsets.push_back(insert_position);
+              if (mapping_info.get_cell_type(cell) == general)
                 {
-                  const unsigned int previous_size =
-                    current_data.jacobians.size();
-                  current_data.rowstart_jacobians.push_back (previous_size);
-                  if (update_flags & update_JxW_values)
-                    {
-                      AssertDimension (previous_size,
-                                       current_data.JxW_values.size());
-                    }
-                  if (update_flags & update_jacobian_grads)
-                    {
-                      AssertDimension (previous_size,
-                                       current_data.jacobians_grad_diag.size());
-                      AssertDimension (previous_size,
-                                       current_data.jacobians_grad_upper.size());
-                    }
                   for (unsigned int q=0; q<n_q_points; ++q)
                     {
-                      Tensor<2,dim,VectorizedArray<Number> > &jac = data.general_jac[q];
-                      Tensor<3,dim,VectorizedArray<Number> > &jacobian_grad = data.general_jac_grad[q];
-                      for (unsigned int j=0; j<vectorization_length; ++j)
-                        if (cell_t[j] == cartesian || cell_t[j] == affine)
+                      Tensor<2,dim,VectorizedArray<Number> > &jac = cell_data.general_jac[q];
+                      Tensor<3,dim,VectorizedArray<Number> > &jacobian_grad = cell_data.general_jac_grad[q];
+                      for (unsigned int j=0; j<VectorizedArray<Number>::n_array_elements; ++j)
+                        if (cell_t[j] < general)
                           {
                             for (unsigned int d=0; d<dim; ++d)
                               for (unsigned int e=0; e<dim; ++e)
                                 {
-                                  jac[d][e][j] = data.const_jac[d][e][j];
+                                  jac[d][e][j] = cell_data.const_jac[d][e][j];
                                   for (unsigned int f=0; f<dim; ++f)
                                     jacobian_grad[d][e][f][j] = 0.;
                                 }
                           }
 
-                      const VectorizedArray<Number> det = determinant (jac);
-                      current_data.jacobians.push_back (transpose(invert(jac)));
-                      const Tensor<2,dim,VectorizedArray<Number> > &inv_jac = current_data.jacobians.back();
-
-                      if (update_flags & update_JxW_values)
-                        current_data.JxW_values.push_back
-                        (det * current_data.quadrature_weights[fe_index][q]);
+                      data.first[my_q].JxW_values.push_back(determinant(jac)*
+                                                            fe_val.get_quadrature().weight(q));
+                      Tensor<2,dim,VectorizedArray<Number> > inv_jac = transpose(invert(jac));
+                      data.first[my_q].jacobians[0].push_back(inv_jac);
 
                       if (update_flags & update_jacobian_grads)
                         {
@@ -498,22 +772,18 @@ namespace internal
                                   }
                               }
 
-                          {
-                            VectorizedArray<Number> grad_diag[dim][dim];
-                            for (unsigned int d=0; d<dim; ++d)
-                              for (unsigned int e=0; e<dim; ++e)
-                                grad_diag[d][e] = grad_jac_inv[d][d][e];
-                            current_data.jacobians_grad_diag.push_back
-                            (Tensor<2,dim,VectorizedArray<Number> >(grad_diag));
-                          }
+                          // the diagonal part of Jacobian gradient comes first
+                          Tensor<1,dim*(dim+1)/2,Tensor<1,dim,VectorizedArray<Number> > > final_grad;
+                          for (unsigned int d=0; d<dim; ++d)
+                            for (unsigned int e=0; e<dim; ++e)
+                              final_grad[d][e] = grad_jac_inv[d][d][e];
 
-                          // sets upper-diagonal part of Jacobian
-                          Tensor<1,(dim>1?dim*(dim-1)/2:1),Tensor<1,dim,VectorizedArray<Number> > > grad_upper;
+                          // then the upper-diagonal part
                           for (unsigned int d=0, count=0; d<dim; ++d)
                             for (unsigned int e=d+1; e<dim; ++e, ++count)
                               for (unsigned int f=0; f<dim; ++f)
-                                grad_upper[count][f] = grad_jac_inv[d][e][f];
-                          current_data.jacobians_grad_upper.push_back(grad_upper);
+                                final_grad[dim+count][f] = grad_jac_inv[d][e][f];
+                          data.first[my_q].jacobian_gradients[0].push_back(final_grad);
                         }
                     }
                 }
@@ -531,397 +801,389 @@ namespace internal
                   // (with some little indexing) the location of all
                   // quadrature points.
                   const unsigned int old_size =
-                    current_data.quadrature_points.size();
-                  current_data.rowstart_q_points[cell] = old_size;
+                    data.first[my_q].quadrature_points.size();
+                  data.first[my_q].quadrature_point_offsets.push_back(old_size);
 
-                  Tensor<1,dim,VectorizedArray<Number> > quad_point;
-
-                  if (get_cell_type(cell) == cartesian)
+                  if (mapping_info.get_cell_type(cell) == cartesian)
                     {
-                      current_data.quadrature_points.resize (old_size+
-                                                             n_q_points_1d[fe_index]);
-                      for (unsigned int q=0; q<n_q_points_1d[fe_index]; ++q)
-                        for (unsigned int d=0; d<dim; ++d)
-                          current_data.quadrature_points[old_size+q][d] =
-                            data.quadrature_points[q*step_size_cartesian[fe_index]][d];
+                      for (unsigned int q=0; q<n_q_points_1d[my_q][fe_index]; ++q)
+                        {
+                          Point<dim,VectorizedArray<Number> > quad_point;
+                          for (unsigned int d=0; d<dim; ++d)
+                            quad_point[d] =
+                              cell_data.quadrature_points[q*step_size_cartesian[my_q][fe_index]][d];
+                          data.first[my_q].quadrature_points.push_back(quad_point);
+                        }
+
                     }
                   else
                     {
-                      current_data.quadrature_points.resize (old_size + n_q_points);
                       for (unsigned int q=0; q<n_q_points; ++q)
-                        for (unsigned int d=0; d<dim; ++d)
-                          current_data.quadrature_points[old_size+q][d] =
-                            data.quadrature_points[q][d];
+                        data.first[my_q].quadrature_points.push_back
+                        (cell_data.quadrature_points[q]);
                     }
                 }
-            } // end for ( cell < n_macro_cells )
-          current_data.rowstart_jacobians.push_back
-          (current_data.jacobians.size());
-          current_data.rowstart_q_points[n_macro_cells] =
-            current_data.quadrature_points.size();
-
-          // finally, fill the accumulated data for Cartesian and affine cells
-          //  into cartesian_data and affine_data, invert and transpose the
-          //  Jacobians, and compute the JxW value.
-          if (my_q == 0)
-            {
-              cartesian_data.resize(cartesians.size());
-              for (typename std::map<Tensor<1,dim,VEC_ARRAY>,
-                   unsigned int, FPArrayComparator<Number> >::iterator
-                   it = cartesians.begin(); it != cartesians.end(); ++it)
-                {
-                  VectorizedArray<Number> det = make_vectorized_array<Number>(1.);
-                  for (unsigned int d=0; d<dim; ++d)
-                    {
-                      VectorizedArray<Number> jac_d;
-                      for (unsigned int v=0;
-                           v<VectorizedArray<Number>::n_array_elements; ++v)
-                        jac_d[v] = it->first[d][v];
-                      cartesian_data[it->second].first[d] = 1./jac_d;
-                      det *= jac_d;
-                    }
-                  cartesian_data[it->second].second = det;
-                }
-              affine_data.resize(affines.size());
-              for (typename std::map<Tensor<2,dim,VEC_ARRAY>,
-                   unsigned int, FPArrayComparator<Number> >::iterator
-                   it = affines.begin(); it != affines.end(); ++it)
-                {
-                  Tensor<2,dim,VectorizedArray<Number> > jac;
-                  for (unsigned int d=0; d<dim; ++d)
-                    for (unsigned int e=0; e<dim; ++e)
-                      for (unsigned int v=0;
-                           v<VectorizedArray<Number>::n_array_elements; ++v)
-                        jac[d][e][v] = it->first[d][e][v];
+            } // end for ( cell < end_cells )
+      }
 
-                  affine_data[it->second].second = determinant(jac);
-                  affine_data[it->second].first = transpose(invert(jac));
-                }
-            }
-        }
-    }
+
+
+      template <typename CONTAINER>
+      void
+      merge_compressed_data(const CONTAINER &source,
+                            CONTAINER &destination,
+                            std::vector<unsigned int> &indices)
+      {
+        indices.resize(source.size());
+        typename CONTAINER::iterator lookup = destination.begin();
+        for (typename CONTAINER::const_iterator it = source.begin();
+             it != source.end(); ++it)
+          {
+            typename CONTAINER::value_type entry = *it;
+            entry.second = destination.size();
+            lookup = destination.insert(lookup, entry);
+            AssertIndexRange(it->second, indices.size());
+            indices[it->second] = lookup->second;
+            // best guess for insert position of next item
+            ++lookup;
+          }
+      }
+
+
+
+      template <int structdim, int dim, typename Number>
+      void
+      copy_data (const unsigned int                first_cell,
+                 const std::array<std::size_t,2>  &data_shift,
+                 const std::vector<unsigned int>  &indices_compressed,
+                 const std::vector<CellType>      &cell_type,
+                 MappingInfoStorage<structdim,dim,Number> &data_cells_local,
+                 MappingInfoStorage<structdim,dim,Number> &data_cells)
+      {
+        // Copy the index offsets and shift by the appropriate value
+        for (unsigned int lcell=0;
+             lcell<data_cells_local.data_index_offsets.size(); ++lcell)
+          {
+            const unsigned int cell = lcell + first_cell;
+            data_cells.data_index_offsets[cell]
+              = cell_type[cell] <= static_cast<unsigned int>(affine) ?
+                indices_compressed[data_cells_local.data_index_offsets[lcell]]
+                :
+                data_cells_local.data_index_offsets[lcell] + data_shift[0];
+            if (data_cells_local.quadrature_point_offsets.size()>lcell)
+              data_cells.quadrature_point_offsets[cell] =
+                data_cells_local.quadrature_point_offsets[lcell] +
+                data_shift[1];
+          }
+
+        // Copy quadrature points
+        if (data_cells.quadrature_point_offsets.empty() == false)
+          {
+            Point<dim,VectorizedArray<Number> > *out_point =
+              &data_cells.quadrature_points[data_shift[1]];
+            for (const Point<dim,VectorizedArray<Number> > *point =
+                   data_cells_local.quadrature_points.begin(); point !=
+                 data_cells_local.quadrature_points.end(); ++point, ++out_point)
+              *out_point = *point;
+            data_cells_local.quadrature_points.clear();
+          }
+
+        // If we have collected Jacobian data, copy Jacobians, JxW values,
+        // Jacobian gradients
+        if (data_cells_local.JxW_values.empty())
+          return;
+
+        std::copy(data_cells_local.JxW_values.begin(),
+                  data_cells_local.JxW_values.end(),
+                  data_cells.JxW_values.begin()+data_shift[0]);
+        data_cells_local.JxW_values.clear();
+        std::copy(data_cells_local.normal_vectors.begin(),
+                  data_cells_local.normal_vectors.end(),
+                  data_cells.normal_vectors.begin()+data_shift[0]);
+        data_cells_local.normal_vectors.clear();
+        for (unsigned int i=0; i<2; ++i)
+          {
+            std::copy(data_cells_local.jacobians[i].begin(),
+                      data_cells_local.jacobians[i].end(),
+                      data_cells.jacobians[i].begin()+data_shift[0]);
+            data_cells_local.jacobians[i].clear();
+            std::copy(data_cells_local.jacobian_gradients[i].begin(),
+                      data_cells_local.jacobian_gradients[i].end(),
+                      data_cells.jacobian_gradients[i].begin()+data_shift[0]);
+            data_cells_local.jacobian_gradients[i].clear();
+            std::copy(data_cells_local.normals_times_jacobians[i].begin(),
+                      data_cells_local.normals_times_jacobians[i].end(),
+                      data_cells.normals_times_jacobians[i].begin()+data_shift[0]);
+            data_cells_local.normals_times_jacobians[i].clear();
+          }
+      }
+
+    } // end of anonymous namespace
 
 
 
     template <int dim, typename Number>
     void
-    MappingInfo<dim,Number>::evaluate_on_cell (const dealii::Triangulation<dim> &tria,
-                                               const std::pair<unsigned int,unsigned int> *cells,
-                                               const unsigned int  cell,
-                                               const unsigned int  my_q,
-                                               CellType (&cell_t_prev)[n_vector_elements],
-                                               CellType (&cell_t)[n_vector_elements],
-                                               dealii::FEValues<dim,dim> &fe_val,
-                                               CellData          &data) const
+    MappingInfo<dim,Number>::initialize_cells
+    (const dealii::Triangulation<dim>                         &tria,
+     const std::vector<std::pair<unsigned int,unsigned int> > &cells,
+     const std::vector<unsigned int>                          &active_fe_index,
+     const Mapping<dim>                                       &mapping,
+     const std::vector<dealii::hp::QCollection<1> >           &quad,
+     const UpdateFlags                                         update_flags_input)
     {
-      const unsigned int n_q_points = fe_val.n_quadrature_points;
-      const UpdateFlags update_flags = fe_val.get_update_flags();
-
-      // this should be the same value as used in HashValue::scaling (but we
-      // not have that field here)
-      const double zero_tolerance_double = data.jac_size *
-                                           std::numeric_limits<double>::epsilon() * 1024.;
-      for (unsigned int j=0; j<n_vector_elements; ++j)
-        {
-          typename dealii::Triangulation<dim>::cell_iterator
-          cell_it (&tria, cells[j].first, cells[j].second);
-          fe_val.reinit(cell_it);
-          cell_t[j] = undefined;
+      const unsigned int n_quads = quad.size();
+      const unsigned int n_cells = cells.size();
+      const unsigned int vectorization_width =
+        VectorizedArray<Number>::n_array_elements;
+      Assert (n_cells%vectorization_width == 0, ExcInternalError());
+      const unsigned int n_macro_cells = n_cells/vectorization_width;
+      cell_data.resize (n_quads);
+      cell_type.resize (n_macro_cells);
 
-          // extract quadrature points and store them temporarily. if we have
-          // Cartesian cells, we can compress the indices
-          if (update_flags & update_quadrature_points)
-            for (unsigned int q=0; q<n_q_points; ++q)
-              {
-                const Point<dim> &point = fe_val.quadrature_point(q);
-                for (unsigned int d=0; d<dim; ++d)
-                  data.quadrature_points[q][d][j] = point[d];
-              }
+      // dummy FE that is used to set up an FEValues object. Do not need the
+      // actual finite element because we will only evaluate quantities for
+      // the mapping that are independent of the FE
+      UpdateFlags update_flags = compute_update_flags (update_flags_input, quad);
 
-          // if this is not the first quadrature formula and we already have
-          // determined that this cell is either Cartesian or with constant
-          // Jacobian, we have nothing more to do.
-          if (my_q > 0 && (get_cell_type(cell) == cartesian
-                           || get_cell_type(cell) == affine) )
-            continue;
-
-          // first round: if the transformation is detected to be the same as
-          // on the old cell, we only need to copy over the data.
-          if (fe_val.get_cell_similarity() == CellSimilarity::translation
-              &&
-              my_q == 0)
-            {
-              if (j==0)
-                {
-                  Assert (cell>0, ExcInternalError());
-                  cell_t[j] = cell_t_prev[n_vector_elements-1];
-                }
-              else
-                cell_t[j] = cell_t[j-1];
-            }
+      for (unsigned int my_q=0; my_q<n_quads; ++my_q)
+        {
+          const unsigned int n_hp_quads = quad[my_q].size();
+          AssertIndexRange (0, n_hp_quads);
+          cell_data[my_q].descriptor.resize(n_hp_quads);
+          for (unsigned int q=0; q<n_hp_quads; ++q)
+            cell_data[my_q].descriptor[q].initialize(quad[my_q][q],
+                                                     update_default);
+        }
 
-          const DerivativeForm<1,dim,dim> &jac_0 = fe_val.jacobian(0);
+      if (n_macro_cells == 0)
+        return;
 
-          if (my_q == 0)
-            {
-              // check whether the Jacobian is constant on this cell the first
-              // time we come around here
-              if (cell_t[j] == undefined)
-                {
-                  bool jacobian_constant = true;
-                  for (unsigned int q=1; q<n_q_points; ++q)
-                    {
-                      const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(q);
-                      for (unsigned int d=0; d<dim; ++d)
-                        for (unsigned int e=0; e<dim; ++e)
-                          if (std::fabs(jac_0[d][e]-jac[d][e]) >
-                              zero_tolerance_double)
-                            jacobian_constant = false;
-                      if (jacobian_constant == false)
-                        break;
-                    }
+      // Create as many chunks of cells as we have threads and spawn the work
+      unsigned int work_per_chunk =
+        std::max(8U, (n_macro_cells + MultithreadInfo::n_threads() - 1) /
+                 MultithreadInfo::n_threads());
 
-                  // check whether the Jacobian is diagonal to machine
-                  // accuracy
-                  bool cell_cartesian = jacobian_constant;
-                  for (unsigned int d=0; d<dim; ++d)
-                    for (unsigned int e=0; e<dim; ++e)
-                      if (d!=e)
-                        if (std::fabs(jac_0[d][e]) >
-                            zero_tolerance_double)
-                          {
-                            cell_cartesian=false;
-                            break;
-                          }
+      std::vector<std::pair<std::vector<MappingInfoStorage<dim,dim,Number> >,
+          CompressedCellData<dim,Number> > > data_cells_local;
+      // Reserve enough space to avoid re-allocation (which would break the
+      // references to the data fields passed to the tasks!)
+      data_cells_local.reserve(MultithreadInfo::n_threads());
 
-                  // in case we have only one quadrature point, we can have
-                  // non-constant Jacobians, but we cannot detect it by
-                  // comparison from one quadrature point to the next: in that
-                  // case, need to look at second derivatives and see whether
-                  // there are some non-zero entries (this is necessary since
-                  // we determine the constness of the Jacobian for the first
-                  // quadrature formula and might not look at them any more
-                  // for the second, third quadrature formula). in any case,
-                  // the flag update_jacobian_grads will be set in that case
-                  if (cell_cartesian == false && n_q_points == 1 &&
-                      update_flags & update_jacobian_grads)
-                    {
-                      const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(0);
-                      const DerivativeForm<2,dim,dim> &jacobian_grad =
-                        fe_val.jacobian_grad(0);
-                      for (unsigned int d=0; d<dim; ++d)
-                        for (unsigned int e=0; e<dim; ++e)
-                          for (unsigned int f=0; f<dim; ++f)
-                            {
-                              double jac_grad_comp = (jac[f][0] *
-                                                      jacobian_grad[d][e][0]);
-                              for (unsigned int g=1; g<dim; ++g)
-                                jac_grad_comp += (jac[f][g] *
-                                                  jacobian_grad[d][e][g]);
-                              if (std::fabs(jac_grad_comp) >
-                                  zero_tolerance_double)
-                                jacobian_constant = false;
-                            }
-                    }
-                  // set cell type
-                  if (cell_cartesian == true)
-                    cell_t[j] = cartesian;
-                  else if (jacobian_constant == true)
-                    cell_t[j] = affine;
-                  else
-                    cell_t[j] = general;
-                }
+      {
+        Threads::TaskGroup<> tasks;
+        std::pair<unsigned int,unsigned int> cell_range(0U, work_per_chunk);
+        while (cell_range.first < n_macro_cells)
+          {
+            data_cells_local.push_back
+            (std::make_pair (std::vector<MappingInfoStorage<dim,dim,Number> >(n_quads),
+                             CompressedCellData<dim,Number>(get_jacobian_size(tria))));
+            tasks += Threads::new_task(&initialize_cell_range<dim,Number>,
+                                       cell_range, tria,
+                                       cells, active_fe_index, mapping, quad,
+                                       update_flags, *this,
+                                       data_cells_local.back());
+            cell_range.first = cell_range.second;
+            cell_range.second += work_per_chunk;
+          }
+        tasks.join_all();
+      }
 
-              // Cartesian cell
-              if (cell_t[j] == cartesian)
-                {
-                  // set Jacobian into diagonal and clear off-diagonal part
-                  for (unsigned int d=0; d<dim; ++d)
-                    {
-                      data.const_jac[d][d][j] = jac_0[d][d];
-                      for (unsigned int e=d+1; e<dim; ++e)
-                        {
-                          data.const_jac[d][e][j] = 0.;
-                          data.const_jac[e][d][j] = 0.;
-                        }
-                    }
-                  continue;
-                }
+      // Fill in each thread's constant Jacobians into the data of the zeroth
+      // chunk in serial
+      std::vector<std::vector<unsigned int> > indices_compressed(data_cells_local.size());
+      for (unsigned int i=0; i<data_cells_local.size(); ++i)
+        merge_compressed_data(data_cells_local[i].second.data,
+                              data_cells_local[0].second.data,
+                              indices_compressed[i]);
+
+      // Collect all data in the final data fields.
+      // First allocate the memory
+      const unsigned int n_constant_jacobians = data_cells_local[0].second.data.size();
+      for (unsigned int my_q=0; my_q<cell_data.size(); ++my_q)
+        {
+          cell_data[my_q].data_index_offsets.resize(cell_type.size());
+          std::vector<std::array<std::size_t,2> > shift(data_cells_local.size());
+          shift[0][0] = n_constant_jacobians;
+          shift[0][1] = 0;
+          for (unsigned int i=1; i<data_cells_local.size(); ++i)
+            {
+              shift[i][0] = shift[i-1][0] + data_cells_local[i-1].first[my_q].JxW_values.size();
+              shift[i][1] = shift[i-1][1] + data_cells_local[i-1].first[my_q].quadrature_points.size();
+            }
+          cell_data[my_q].JxW_values.
+          resize_fast(shift.back()[0] + data_cells_local.back().first[my_q].
+                      JxW_values.size());
+          cell_data[my_q].jacobians[0].resize_fast(cell_data[my_q].JxW_values.size());
+          if (update_flags & update_jacobian_grads)
+            cell_data[my_q].jacobian_gradients[0].resize_fast(cell_data[my_q].JxW_values.size());
+          if (update_flags & update_quadrature_points)
+            {
+              cell_data[my_q].quadrature_point_offsets.resize(cell_type.size());
+              cell_data[my_q].quadrature_points.
+              resize_fast(shift.back()[1] + data_cells_local.back().first[my_q].
+                          quadrature_points.size());
+            }
 
-              // cell with affine mapping
-              else if (cell_t[j] == affine)
+          // Start tasks that copy the local data
+          Threads::TaskGroup<> tasks;
+          for (unsigned int i=0; i<data_cells_local.size(); ++i)
+            tasks += Threads::new_task(&copy_data<dim,dim,Number>,
+                                       work_per_chunk * i, shift[i],
+                                       indices_compressed[i], cell_type,
+                                       data_cells_local[i].first[my_q],
+                                       cell_data[my_q]);
+
+          // finally, insert the constant cell data at the beginning (the
+          // other tasks can already start copying the non-constant data)
+          if (my_q == 0)
+            {
+              for (auto &it : data_cells_local[0].second.data)
                 {
-                  // compress out very small values
+                  Tensor<2,dim,VectorizedArray<Number> > jac;
                   for (unsigned int d=0; d<dim; ++d)
                     for (unsigned int e=0; e<dim; ++e)
-                      data.const_jac[d][e][j] =
-                        std::fabs(jac_0[d][e]) < zero_tolerance_double ?
-                        0 : jac_0[d][e];
-                  continue;
+                      for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+                        jac[d][e][v] = it.first[d][e][v];
+                  AssertIndexRange(it.second, n_constant_jacobians);
+                  const std::size_t index = it.second;
+                  cell_data[my_q].JxW_values[index] = determinant(jac);
+                  // invert and transpose jac
+                  cell_data[my_q].jacobians[0][index] = transpose(invert(jac));
+                  // second derivative of transformation is zero on affine cells
                 }
             }
-
-          // general cell case
-
-          // go through all quadrature points and fill in the data into the
-          // temporary data structures with slots for the vectorized data
-          // types
-          for (unsigned int q=0; q<n_q_points; ++q)
+          else
             {
-
-              // compress out very small numbers which are only noise. Then it
-              // is cleaner to use zero straight away (though it does not save
-              // any memory)
-              const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(q);
-              for (unsigned int d=0; d<dim; ++d)
-                for (unsigned int e=0; e<dim; ++e)
-                  data.general_jac[q][d][e][j] =
-                    std::fabs(jac[d][e]) < zero_tolerance_double ? 0. : jac[d][e];
-
-              // need to do some calculus based on the gradient of the
-              // Jacobian, in order to find the gradient of the inverse
-              // Jacobian which is needed in user code. however, we would like
-              // to perform that on vectorized data types instead of doubles
-              // or floats. to this end, copy the gradients first
-              if (update_flags & update_jacobian_grads)
+              for (unsigned int i=0; i<n_constant_jacobians; ++i)
                 {
-                  const DerivativeForm<2,dim,dim> &jacobian_grad = fe_val.jacobian_grad(q);
-                  for (unsigned int d=0; d<dim; ++d)
-                    for (unsigned int e=0; e<dim; ++e)
-                      for (unsigned int f=0; f<dim; ++f)
-                        data.general_jac_grad[q][d][e][f][j] = jacobian_grad[d][e][f];
+                  cell_data[my_q].JxW_values[i]   = cell_data[0].JxW_values[i];
+                  cell_data[my_q].jacobians[0][i] = cell_data[0].jacobians[0][i];
                 }
             }
-        } // end loop over all entries in vectorization (n_vector_elements
-      // cells)
 
-      // set information for next cell
-      for (unsigned int j=0; j<n_vector_elements; ++j)
-        cell_t_prev[j] = cell_t[j];
+          // ... wait for the parallel work to finish
+          tasks.join_all();
+        }
     }
 
 
+
     template <int dim, typename Number>
-    MappingInfo<dim,Number>::CellData::CellData (const double jac_size_in)
-      :
-      jac_size (jac_size_in)
-    {}
+    std::size_t MappingInfo<dim,Number>::memory_consumption() const
+    {
+      std::size_t
+      memory  = MemoryConsumption::memory_consumption (cell_data);
+      memory += MemoryConsumption::memory_consumption (face_data);
+      memory += cell_type.capacity()*sizeof(CellType);
+      memory += face_type.capacity()*sizeof(CellType);
+      memory += sizeof (*this);
+      return memory;
+    }
 
 
 
     template <int dim, typename Number>
-    void
-    MappingInfo<dim,Number>::CellData::resize (const unsigned int size)
+    template <typename StreamType>
+    void MappingInfo<dim,Number>::print_memory_consumption(StreamType     &out,
+                                                           const SizeInfo &task_info) const
     {
-      if (general_jac.size() != size)
+      out << "    Cell types:                      ";
+      task_info.print_memory_statistics
+      (out, cell_type.capacity()*sizeof(CellType));
+      out << "    Face types:                      ";
+      task_info.print_memory_statistics
+      (out, face_type.capacity()*sizeof(CellType));
+      for (unsigned int j=0; j<cell_data.size(); ++j)
         {
-          quadrature_points.resize(size);
-          general_jac.resize(size);
-          general_jac_grad.resize(size);
+          out << "    Data component " << j << std::endl;
+          cell_data[j].print_memory_consumption(out, task_info);
+          face_data[j].print_memory_consumption(out, task_info);
         }
     }
 
 
 
-    template <int dim, typename Number>
-    std::size_t MappingInfo<dim,Number>::MappingInfoDependent::memory_consumption() const
+    /* ------------------------------------------------------------------ */
+
+    template <typename Number>
+    FPArrayComparator<Number>::FPArrayComparator (const Number scaling)
+      :
+      tolerance (scaling * std::numeric_limits<double>::epsilon() * 1024.)
+    {}
+
+
+
+    template <typename Number>
+    bool
+    FPArrayComparator<Number>::operator() (const std::vector<Number> &v1,
+                                           const std::vector<Number> &v2) const
     {
-      std::size_t
-      memory = MemoryConsumption::memory_consumption (jacobians);
-      memory += MemoryConsumption::memory_consumption (JxW_values);
-      memory += MemoryConsumption::memory_consumption (jacobians_grad_diag);
-      memory += MemoryConsumption::memory_consumption (jacobians_grad_upper);
-      memory += MemoryConsumption::memory_consumption (rowstart_q_points);
-      memory += MemoryConsumption::memory_consumption (quadrature_points);
-      memory += MemoryConsumption::memory_consumption (quadrature);
-      memory += MemoryConsumption::memory_consumption (face_quadrature);
-      memory += MemoryConsumption::memory_consumption (quadrature_weights);
-      memory += MemoryConsumption::memory_consumption (n_q_points);
-      memory += MemoryConsumption::memory_consumption (n_q_points_face);
-      memory += MemoryConsumption::memory_consumption (quad_index_conversion);
-      return memory;
+      const unsigned int s1 = v1.size(), s2 = v2.size();
+      if (s1 < s2)
+        return true;
+      else if (s1 > s2)
+        return false;
+      else
+        for (unsigned int i=0; i<s1; ++i)
+          if (v1[i] < v2[i] - tolerance)
+            return true;
+          else if (v1[i] > v2[i] + tolerance)
+            return false;
+      return false;
     }
 
 
 
-    template <int dim, typename Number>
-    std::size_t MappingInfo<dim,Number>::memory_consumption() const
+    template <typename Number>
+    bool
+    FPArrayComparator<Number>::
+    operator ()(const Tensor<1,VectorizedArray<Number>::n_array_elements,Number> &t1,
+                const Tensor<1,VectorizedArray<Number>::n_array_elements,Number> &t2) const
     {
-      std::size_t
-      memory= MemoryConsumption::memory_consumption (mapping_data_gen);
-      memory += MemoryConsumption::memory_consumption (affine_data);
-      memory += MemoryConsumption::memory_consumption (cartesian_data);
-      memory += MemoryConsumption::memory_consumption (cell_type);
-      memory += sizeof (*this);
-      return memory;
+      for (unsigned int k=0; k<VectorizedArray<Number>::n_array_elements; ++k)
+        if (t1[k] < t2[k] - tolerance)
+          return true;
+        else if (t1[k] > t2[k] + tolerance)
+          return false;
+      return false;
     }
 
 
 
-    template <int dim, typename Number>
-    template <typename StreamType>
-    void MappingInfo<dim,Number>::MappingInfoDependent::print_memory_consumption
-    (StreamType     &out,
-     const SizeInfo &size_info) const
+    template <typename Number>
+    template <int dim>
+    bool
+    FPArrayComparator<Number>::
+    operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
+                const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const
     {
-      // print_memory_statistics involves global communication, so we can
-      // disable the check here only if no processor has any such data
-#ifdef DEAL_II_WITH_MPI
-      unsigned int general_size_glob = 0, general_size_loc = jacobians.size();
-      int ierr = MPI_Allreduce (&general_size_loc, &general_size_glob, 1,
-                                MPI_UNSIGNED, MPI_MAX, size_info.communicator);
-      AssertThrowMPI (ierr);
-#else
-      unsigned int general_size_glob = jacobians.size();
-#endif
-      if (general_size_glob > 0)
-        {
-          out << "      Memory Jacobian data:          ";
-          size_info.print_memory_statistics
-          (out, MemoryConsumption::memory_consumption (jacobians) +
-           MemoryConsumption::memory_consumption (JxW_values));
-          out << "      Memory second derivative data: ";
-          size_info.print_memory_statistics
-          (out,MemoryConsumption::memory_consumption (jacobians_grad_diag) +
-           MemoryConsumption::memory_consumption (jacobians_grad_upper));
-        }
-
-#ifdef DEAL_II_WITH_MPI
-      unsigned int quad_size_glob = 0, quad_size_loc = quadrature_points.size();
-      ierr = MPI_Allreduce (&quad_size_loc, &quad_size_glob, 1, MPI_UNSIGNED,
-                            MPI_MAX, size_info.communicator);
-      AssertThrowMPI (ierr);
-#else
-      unsigned int quad_size_glob = quadrature_points.size();
-#endif
-      if (quad_size_glob > 0)
-        {
-          out << "      Memory quadrature points:      ";
-          size_info.print_memory_statistics
-          (out, MemoryConsumption::memory_consumption (rowstart_q_points) +
-           MemoryConsumption::memory_consumption (quadrature_points));
-        }
+      for (unsigned int d=0; d<dim; ++d)
+        for (unsigned int k=0; k<VectorizedArray<Number>::n_array_elements; ++k)
+          if (t1[d][k] < t2[d][k] - tolerance)
+            return true;
+          else if (t1[d][k] > t2[d][k] + tolerance)
+            return false;
+      return false;
     }
 
 
 
-    template <int dim, typename Number>
-    template <typename StreamType>
-    void MappingInfo<dim,Number>::print_memory_consumption(StreamType     &out,
-                                                           const SizeInfo &size_info) const
+    template <typename Number>
+    template <int dim>
+    bool
+    FPArrayComparator<Number>::
+    operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
+                const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const
     {
-      out << "    Cell types:                      ";
-      size_info.print_memory_statistics
-      (out, MemoryConsumption::memory_consumption (cell_type));
-      out << "    Memory transformations compr:    ";
-      size_info.print_memory_statistics
-      (out, MemoryConsumption::memory_consumption (affine_data) +
-       MemoryConsumption::memory_consumption (cartesian_data));
-      for (unsigned int j=0; j<mapping_data_gen.size(); ++j)
-        {
-          out << "    Data component " << j << std::endl;
-          mapping_data_gen[j].print_memory_consumption(out, size_info);
-        }
+      for (unsigned int d=0; d<dim; ++d)
+        for (unsigned int e=0; e<dim; ++e)
+          for (unsigned int k=0; k<VectorizedArray<Number>::n_array_elements; ++k)
+            if (t1[d][e][k] < t2[d][e][k] - tolerance)
+              return true;
+            else if (t1[d][e][k] > t2[d][e][k] + tolerance)
+              return false;
+      return false;
     }
 
   } // end of namespace MatrixFreeFunctions
index f94a1bf8655be64f06a7e0fc0b34127e7ebbc360..4d3daa83edd36f9051d6a4f839d1a74f8adce30d 100644 (file)
@@ -1397,9 +1397,8 @@ unsigned int
 MatrixFree<dim,Number>::get_n_q_points(const unsigned int quad_index,
                                        const unsigned int active_fe_index) const
 {
-  AssertIndexRange (quad_index,
-                    mapping_info.mapping_data_gen.size());
-  return mapping_info.mapping_data_gen[quad_index].n_q_points[active_fe_index];
+  AssertIndexRange (quad_index, mapping_info.cell_data.size());
+  return mapping_info.cell_data[quad_index].descriptor[active_fe_index].n_q_points;
 }
 
 
@@ -1422,9 +1421,8 @@ unsigned int
 MatrixFree<dim,Number>::get_n_q_points_face(const unsigned int quad_index,
                                             const unsigned int active_fe_index) const
 {
-  AssertIndexRange (quad_index,
-                    mapping_info.mapping_data_gen.size());
-  return mapping_info.mapping_data_gen[quad_index].n_q_points_face[active_fe_index];
+  AssertIndexRange (quad_index, mapping_info.face_data.size());
+  return mapping_info.face_data[quad_index].descriptor[active_fe_index].n_q_points;
 }
 
 
@@ -1475,9 +1473,8 @@ const Quadrature<dim> &
 MatrixFree<dim,Number>::get_quadrature (const unsigned int quad_index,
                                         const unsigned int active_fe_index) const
 {
-  AssertIndexRange (quad_index, mapping_info.mapping_data_gen.size());
-  return mapping_info.mapping_data_gen[quad_index].
-         quadrature[active_fe_index];
+  AssertIndexRange (quad_index, mapping_info.cell_data.size());
+  return mapping_info.cell_data[quad_index].descriptor[active_fe_index].quadrature;
 }
 
 
@@ -1488,9 +1485,8 @@ const Quadrature<dim-1> &
 MatrixFree<dim,Number>::get_face_quadrature (const unsigned int quad_index,
                                              const unsigned int active_fe_index) const
 {
-  AssertIndexRange (quad_index, mapping_info.mapping_data_gen.size());
-  return mapping_info.mapping_data_gen[quad_index].
-         face_quadrature[active_fe_index];
+  AssertIndexRange (quad_index, mapping_info.face_data.size());
+  return mapping_info.face_data[quad_index].descriptor[active_fe_index].quadrature;
 }
 
 
index 22b8eade937b24cf87b3316b7290a99fe0475594..0864dc8f028cc9edb972dc65e23624cc327faa9c 100644 (file)
@@ -1034,73 +1034,6 @@ namespace internal
         boundary_cells_start = boundary_cells_end = n_macro_cells;
     }
 
-
-
-    /* ------------------------------------------------------------------ */
-
-    template <typename Number>
-    FPArrayComparator<Number>::FPArrayComparator (const Number scaling)
-      :
-      tolerance (scaling * std::numeric_limits<double>::epsilon() * 1024.)
-    {}
-
-
-
-    template <typename Number>
-    bool
-    FPArrayComparator<Number>::operator() (const std::vector<Number> &v1,
-                                           const std::vector<Number> &v2) const
-    {
-      const unsigned int s1 = v1.size(), s2 = v2.size();
-      if (s1 < s2)
-        return true;
-      else if (s1 > s2)
-        return false;
-      else
-        for (unsigned int i=0; i<s1; ++i)
-          if (v1[i] < v2[i] - tolerance)
-            return true;
-          else if (v1[i] > v2[i] + tolerance)
-            return false;
-      return false;
-    }
-
-
-
-    template <typename Number>
-    template <int dim>
-    bool
-    FPArrayComparator<Number>::
-    operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
-                const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const
-    {
-      for (unsigned int d=0; d<dim; ++d)
-        for (unsigned int k=0; k<VectorizedArray<Number>::n_array_elements; ++k)
-          if ((t1)[d][k] < (t2)[d][k] - tolerance)
-            return true;
-          else if ((t1)[d][k] > (t2)[d][k] + tolerance)
-            return false;
-      return false;
-    }
-
-
-
-    template <typename Number>
-    template <int dim>
-    bool
-    FPArrayComparator<Number>::
-    operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
-                const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const
-    {
-      for (unsigned int d=0; d<dim; ++d)
-        for (unsigned int e=0; e<dim; ++e)
-          for (unsigned int k=0; k<VectorizedArray<Number>::n_array_elements; ++k)
-            if ((t1)[d][e][k] < (t2)[d][e][k] - tolerance)
-              return true;
-            else if ((t1)[d][e][k] > (t2)[d][e][k] + tolerance)
-              return false;
-      return false;
-    }
   }
 }
 
index a79ed61004bdb6d95bd83972f6d4936337c45c34..b15a166af640f898dd0f8c33cce0bacd4b8847f3 100644 (file)
@@ -32,6 +32,11 @@ for (deal_II_dimension : DIMENSIONS)
     template struct internal::MatrixFreeFunctions::MappingInfo<deal_II_dimension,double>;
     template struct internal::MatrixFreeFunctions::MappingInfo<deal_II_dimension,float>;
 
+    template struct internal::MatrixFreeFunctions::MappingInfoStorage<deal_II_dimension,deal_II_dimension,double>;
+    template struct internal::MatrixFreeFunctions::MappingInfoStorage<deal_II_dimension,deal_II_dimension,float>;
+    template struct internal::MatrixFreeFunctions::MappingInfoStorage<deal_II_dimension-1,deal_II_dimension,double>;
+    template struct internal::MatrixFreeFunctions::MappingInfoStorage<deal_II_dimension-1,deal_II_dimension,float>;
+
 #ifndef DEAL_II_MSVC
     template
     void
index 777fa0f01e60636635c5da7968a2b4aea1e7bc8a..15070956665c226257ee347794eab643b91fdd61 100644 (file)
@@ -79,14 +79,15 @@ void test ()
   mf.reinit (dof, constraints, quad, data);
 
   const unsigned int n_macro_cells = mf.n_macro_cells();
-  const unsigned int n_cartesian = mf.get_mapping_info().cartesian_data.size();
-  const unsigned int n_affine = mf.get_mapping_info().affine_data.size();
-  const unsigned int n_general = mf.get_mapping_info().mapping_data_gen[0].rowstart_jacobians.size()-1;
+  std::vector<unsigned int> n_cell_types(4, 0);
+  for (unsigned int i=0; i<n_macro_cells; ++i)
+    n_cell_types[mf.get_mapping_info().get_cell_type(i)]++;
 
   // should do at least some compression
-  Assert(n_cartesian+n_affine+n_general < n_macro_cells, ExcInternalError());
-  Assert(n_cartesian * 5 < n_macro_cells, ExcInternalError());
-  Assert(n_affine * 10 < n_macro_cells, ExcInternalError());
+  Assert(n_cell_types[0]+n_cell_types[1] > 0, ExcInternalError());
+  Assert(mf.get_mapping_info().cell_data[0].jacobians[0].size() <
+         (n_cell_types[3]*quad.size()+n_macro_cells-n_cell_types[3]),
+         ExcInternalError());
   deallog << "OK" << std::endl;
 }
 
@@ -115,15 +116,14 @@ void test_cube ()
   mf.reinit (dof, constraints, quad, data);
 
   const unsigned int n_macro_cells = mf.n_macro_cells();
-  const unsigned int n_cartesian = mf.get_mapping_info().cartesian_data.size();
-  const unsigned int n_affine = mf.get_mapping_info().affine_data.size();
-  const unsigned int n_general = mf.get_mapping_info().mapping_data_gen[0].rowstart_jacobians.size()-1;
+  std::vector<unsigned int> n_cell_types(4, 0);
+  for (unsigned int i=0; i<n_macro_cells; ++i)
+    n_cell_types[mf.get_mapping_info().get_cell_type(i)]++;
 
   // should have one Cartesian cell and no other
   // cell type
-  AssertDimension(n_cartesian, 1);
-  AssertDimension(n_affine, 0);
-  AssertDimension(n_general, 0);
+  AssertDimension(n_cell_types[0], n_macro_cells);
+  AssertDimension(mf.get_mapping_info().cell_data[0].jacobians[0].size(), 1);
   Assert(n_macro_cells > 1, ExcInternalError());
   deallog << "OK" << std::endl;
 }
@@ -160,15 +160,14 @@ void test_parallelogram ()
   mf.reinit (dof, constraints, quad, data);
 
   const unsigned int n_macro_cells = mf.n_macro_cells();
-  const unsigned int n_cartesian = mf.get_mapping_info().cartesian_data.size();
-  const unsigned int n_affine = mf.get_mapping_info().affine_data.size();
-  const unsigned int n_general = mf.get_mapping_info().mapping_data_gen[0].rowstart_jacobians.size()-1;
+  std::vector<unsigned int> n_cell_types(4, 0);
+  for (unsigned int i=0; i<n_macro_cells; ++i)
+    n_cell_types[mf.get_mapping_info().get_cell_type(i)]++;
 
   // should have one affine cell and no other
   // cell type
-  AssertDimension(n_cartesian, 0);
-  AssertDimension(n_affine, 1);
-  AssertDimension(n_general, 0);
+  AssertDimension(n_cell_types[1], n_macro_cells);
+  AssertDimension(mf.get_mapping_info().cell_data[0].jacobians[0].size(), 1);
   Assert(n_macro_cells > 1, ExcInternalError());
   deallog << "OK" << std::endl;
 }

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.