From 841783bd211631bc601ee716e61d907ea9e854e7 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Sun, 1 Nov 2015 09:56:00 -0600
Subject: [PATCH] Remove the first_cell/is_first_cell/current_update_flags()
 mechanism from FiniteElement.

What this patch does is address the issue that
`FiniteElement::get_data()` (as was previously the case with Mappings
as well) sets a flag in the returned `InternalData` object that
indicates that we are visiting a cell for the first time. `FEValues`
later sets this flag back once we're done with the first cell.

The way this is used is that the various finite element
implementations query `current_update_flags` in
`FE::fill_fe_*_values()`, which equals `update_each|update_once` if we
are on the first cell, and `update_each` on later cells. That means
that on the first cell we update fields in the output object that we
need to initialize only once, e.g., the values of shape functions on
quadrature points, but we won't do this again on later cells.

I don't particularly like this approach. It is opaque and took even me
a long time to reconstruct when I saw it (see #1305 for a
discussion). My preference is to do things we only need to do once in
`FE::get_data()`, and do things we need to do every time in
`FE::fill_fe_*_values()`. This kind of change was already implemented
for mappings in the patches references from #1305, and the current
pull request goes into the same direction.

The current PR only does one half of the necessary changes in order to
keep the patch reasonably self contained and readable. In particular,
what it really does is simply *always* do everything, by storing
`update_once|update_each` in `FE::InternalDataBase::update_each`. In
other words, it is the equivalent to telling finite element
implementations that we are always on the first cell. This requires
more work than in the previous state, although copying shape values
every time is arguably of small expense compared to transforming
gradients.

I will of course fix this up in upcoming patches. It requires
significant shuffling around. I'll open an issue in a minute
discussing the details. I'd simply like to propose this here already
as an incremental step forward, without making patch review too
onerous. The patch is mostly mechanical once you understand what it
does.
---
 include/deal.II/fe/fe.h                     | 30 -----------
 include/deal.II/fe/fe_face.h                | 10 +---
 include/deal.II/fe/fe_poly.h                | 30 ++++-------
 include/deal.II/fe/fe_poly.templates.h      |  6 +--
 include/deal.II/fe/fe_poly_face.h           | 14 ++---
 include/deal.II/fe/fe_poly_face.templates.h | 18 ++++---
 include/deal.II/fe/fe_poly_tensor.h         | 26 ++++-----
 include/deal.II/fe/fe_system.h              |  9 ----
 source/fe/fe.cc                             | 31 +----------
 source/fe/fe_dgp_nonparametric.cc           | 59 +++++++++------------
 source/fe/fe_face.cc                        |  4 +-
 source/fe/fe_poly.cc                        | 40 ++++++--------
 source/fe/fe_poly_tensor.cc                 | 33 +++++-------
 source/fe/fe_system.cc                      | 49 +++--------------
 source/fe/fe_values.cc                      |  6 ---
 15 files changed, 103 insertions(+), 262 deletions(-)

diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h
index c5a990ff67..add7dbc80c 100644
--- a/include/deal.II/fe/fe.h
+++ b/include/deal.II/fe/fe.h
@@ -519,16 +519,6 @@ public:
      */
     virtual ~InternalDataBase ();
 
-    /**
-     * Values updated by the constructor or by reinit.
-     */
-    UpdateFlags          update_flags;
-
-    /**
-     * Values computed by constructor.
-     */
-    UpdateFlags          update_once;
-
     /**
      * A set of update flags specifying the kind of information that
      * an implementation of the FiniteElement interface needs to compute on
@@ -547,30 +537,10 @@ public:
      */
     UpdateFlags          update_each;
 
-    /**
-     * If <tt>first_cell==true</tt> this function returns @p update_flags,
-     * i.e. <tt>update_once|update_each</tt>. If <tt>first_cell==false</tt> it
-     * returns @p update_each.
-     */
-    UpdateFlags  current_update_flags() const;
-
-    /**
-     * Set the @p first_cell flag to @p false. Used by the @p FEValues class
-     * to indicate that we have already done the work on the first cell.
-     */
-    virtual void clear_first_cell ();
-
     /**
      * Return an estimate (in bytes) or the memory consumption of this object.
      */
     virtual std::size_t memory_consumption () const;
-
-  private:
-    /**
-     * Initially set to true, but reset to false when clear_first_cell()
-     * is called.
-     */
-    bool first_cell;
   };
 
 public:
diff --git a/include/deal.II/fe/fe_face.h b/include/deal.II/fe/fe_face.h
index 352bc36399..ad9cfa624b 100644
--- a/include/deal.II/fe/fe_face.h
+++ b/include/deal.II/fe/fe_face.h
@@ -250,20 +250,14 @@ protected:
     // generate a new data object and initialize some fields
     typename FiniteElement<1,spacedim>::InternalDataBase *data =
       new typename FiniteElement<1,spacedim>::InternalDataBase;
+    data->update_each = update_once(update_flags) | update_each(update_flags);  // FIX: only update_each required
 
-    // check what needs to be initialized only once and what on every
-    // cell/face/subface we visit
-    data->update_once = update_once(update_flags);
-    data->update_each = update_each(update_flags);
-    data->update_flags = data->update_once | data->update_each;
-
-    const UpdateFlags flags(data->update_flags);
     const unsigned int n_q_points = quadrature.size();
     AssertDimension(n_q_points, 1);
     (void)n_q_points;
 
     // No derivatives of this element are implemented.
-    if (flags & update_gradients || flags & update_hessians)
+    if (data->update_each & update_gradients || data->update_each & update_hessians)
       {
         Assert(false, ExcNotImplemented());
       }
diff --git a/include/deal.II/fe/fe_poly.h b/include/deal.II/fe/fe_poly.h
index 79331d4e85..c73b58117c 100644
--- a/include/deal.II/fe/fe_poly.h
+++ b/include/deal.II/fe/fe_poly.h
@@ -222,16 +222,8 @@ protected:
     // generate a new data object and
     // initialize some fields
     InternalData *data = new InternalData;
+    data->update_each = update_each(update_flags) | update_once(update_flags);  // FIX: only update_each required
 
-    // check what needs to be
-    // initialized only once and what
-    // on every cell/face/subface we
-    // visit
-    data->update_once = update_once(update_flags);
-    data->update_each = update_each(update_flags);
-    data->update_flags = data->update_once | data->update_each;
-
-    const UpdateFlags flags(data->update_flags);
     const unsigned int n_q_points = quadrature.size();
 
     // some scratch arrays
@@ -244,28 +236,28 @@ protected:
     // initialize fields only if really
     // necessary. otherwise, don't
     // allocate memory
-    if (flags & update_values)
+    if (update_flags & update_values)
       {
         values.resize (this->dofs_per_cell);
         data->shape_values.resize (this->dofs_per_cell,
                                    std::vector<double> (n_q_points));
       }
 
-    if (flags & update_gradients)
+    if (update_flags & update_gradients)
       {
         grads.resize (this->dofs_per_cell);
         data->shape_gradients.resize (this->dofs_per_cell,
                                       std::vector<Tensor<1,dim> > (n_q_points));
       }
 
-    if (flags & update_hessians)
+    if (update_flags & update_hessians)
       {
         grad_grads.resize (this->dofs_per_cell);
         data->shape_hessians.resize (this->dofs_per_cell,
                                      std::vector<Tensor<2,dim> > (n_q_points));
       }
 
-    if (flags & update_3rd_derivatives)
+    if (update_flags & update_3rd_derivatives)
       {
         third_derivatives.resize (this->dofs_per_cell);
         data->shape_3rd_derivatives.resize (this->dofs_per_cell,
@@ -279,8 +271,8 @@ protected:
     // unit cell, and need to be
     // transformed when visiting an
     // actual cell
-    if (flags & (update_values | update_gradients
-                 | update_hessians | update_3rd_derivatives) )
+    if (update_flags & (update_values | update_gradients
+                        | update_hessians | update_3rd_derivatives) )
       for (unsigned int i=0; i<n_q_points; ++i)
         {
           poly_space.compute(quadrature.point(i),
@@ -288,19 +280,19 @@ protected:
                              third_derivatives,
                              fourth_derivatives);
 
-          if (flags & update_values)
+          if (update_flags & update_values)
             for (unsigned int k=0; k<this->dofs_per_cell; ++k)
               data->shape_values[k][i] = values[k];
 
-          if (flags & update_gradients)
+          if (update_flags & update_gradients)
             for (unsigned int k=0; k<this->dofs_per_cell; ++k)
               data->shape_gradients[k][i] = grads[k];
 
-          if (flags & update_hessians)
+          if (update_flags & update_hessians)
             for (unsigned int k=0; k<this->dofs_per_cell; ++k)
               data->shape_hessians[k][i] = grad_grads[k];
 
-          if (flags & update_3rd_derivatives)
+          if (update_flags & update_3rd_derivatives)
             for (unsigned int k=0; k<this->dofs_per_cell; ++k)
               data->shape_3rd_derivatives[k][i] = third_derivatives[k];
         }
diff --git a/include/deal.II/fe/fe_poly.templates.h b/include/deal.II/fe/fe_poly.templates.h
index bdaa97d2e1..74519fda03 100644
--- a/include/deal.II/fe/fe_poly.templates.h
+++ b/include/deal.II/fe/fe_poly.templates.h
@@ -240,7 +240,7 @@ fill_fe_values (const Mapping<dim,spacedim>                                  &ma
   Assert (dynamic_cast<const InternalData *> (&fedata) != 0, ExcInternalError());
   const InternalData &fe_data = static_cast<const InternalData &> (fedata);
 
-  const UpdateFlags flags(fe_data.current_update_flags());
+  const UpdateFlags flags(fe_data.update_each);
 
   for (unsigned int k=0; k<this->dofs_per_cell; ++k)
     {
@@ -312,7 +312,7 @@ fill_fe_face_values (const Mapping<dim,spacedim>
                                                 cell->face_rotation(face),
                                                 quadrature.size());
 
-  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
+  const UpdateFlags flags(fe_data.update_each);
 
   for (unsigned int k=0; k<this->dofs_per_cell; ++k)
     {
@@ -388,7 +388,7 @@ fill_fe_subface_values (const Mapping<dim,spacedim>
                                                    quadrature.size(),
                                                    cell->subface_case(face));
 
-  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
+  const UpdateFlags flags(fe_data.update_each);
 
   for (unsigned int k=0; k<this->dofs_per_cell; ++k)
     {
diff --git a/include/deal.II/fe/fe_poly_face.h b/include/deal.II/fe/fe_poly_face.h
index 697d5e801c..31d85ea509 100644
--- a/include/deal.II/fe/fe_poly_face.h
+++ b/include/deal.II/fe/fe_poly_face.h
@@ -99,16 +99,8 @@ protected:
     // generate a new data object and
     // initialize some fields
     InternalData *data = new InternalData;
+    data->update_each = update_once(update_flags) | update_each(update_flags);  // FIX: only update_each required
 
-    // check what needs to be
-    // initialized only once and what
-    // on every cell/face/subface we
-    // visit
-    data->update_once = update_once(update_flags);
-    data->update_each = update_each(update_flags);
-    data->update_flags = data->update_once | data->update_each;
-
-    const UpdateFlags flags(data->update_flags);
     const unsigned int n_q_points = quadrature.size();
 
     // some scratch arrays
@@ -121,7 +113,7 @@ protected:
     // initialize fields only if really
     // necessary. otherwise, don't
     // allocate memory
-    if (flags & update_values)
+    if (data->update_each & update_values)
       {
         values.resize (poly_space.n());
         data->shape_values.resize (poly_space.n(),
@@ -139,7 +131,7 @@ protected:
       }
     // No derivatives of this element
     // are implemented.
-    if (flags & update_gradients || flags & update_hessians)
+    if (data->update_each & update_gradients || data->update_each & update_hessians)
       {
         Assert(false, ExcNotImplemented());
       }
diff --git a/include/deal.II/fe/fe_poly_face.templates.h b/include/deal.II/fe/fe_poly_face.templates.h
index 786f0bd3e9..5497997a91 100644
--- a/include/deal.II/fe/fe_poly_face.templates.h
+++ b/include/deal.II/fe/fe_poly_face.templates.h
@@ -118,9 +118,7 @@ fill_fe_face_values (const Mapping<dim,spacedim> &,
   Assert (dynamic_cast<const InternalData *> (&fedata) != 0, ExcInternalError());
   const InternalData &fe_data = static_cast<const InternalData &> (fedata);
 
-  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
-
-  if (flags & update_values)
+  if (fe_data.update_each & update_values)
     for (unsigned int i=0; i<quadrature.size(); ++i)
       {
         for (unsigned int k=0; k<this->dofs_per_cell; ++k)
@@ -136,7 +134,10 @@ fill_fe_face_values (const Mapping<dim,spacedim> &,
                 for (unsigned int k=0; k<this->dofs_per_quad; ++k)
                   output_data.shape_values(foffset+k,i) = fe_data.shape_values[k+this->first_face_quad_index][i];
               }
+
+            // fall through...
           }
+
           case 2:
           {
             // Fill data for line shape functions
@@ -150,7 +151,10 @@ fill_fe_face_values (const Mapping<dim,spacedim> &,
                         = fe_data.shape_values[k+(line*this->dofs_per_line)+this->first_face_line_index][i];
                   }
               }
+
+            // fall through...
           }
+
           case 1:
           {
             // Fill data for vertex shape functions
@@ -185,12 +189,10 @@ fill_fe_subface_values (const Mapping<dim,spacedim> &,
   Assert (dynamic_cast<const InternalData *> (&fedata) != 0, ExcInternalError());
   const InternalData &fe_data = static_cast<const InternalData &> (fedata);
 
-  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
-
   const unsigned int foffset = fe_data.shape_values.size() * face;
   const unsigned int offset = subface*quadrature.size();
 
-  if (flags & update_values)
+  if (fe_data.update_each & update_values)
     for (unsigned int i=0; i<quadrature.size(); ++i)
       {
         for (unsigned int k=0; k<this->dofs_per_cell; ++k)
@@ -199,8 +201,8 @@ fill_fe_subface_values (const Mapping<dim,spacedim> &,
           output_data.shape_values(foffset+k,i) = fe_data.shape_values[k][i+offset];
       }
 
-  Assert (!(flags & update_gradients), ExcNotImplemented());
-  Assert (!(flags & update_hessians), ExcNotImplemented());
+  Assert (!(fe_data.update_each & update_gradients), ExcNotImplemented());
+  Assert (!(fe_data.update_each & update_hessians), ExcNotImplemented());
 }
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/include/deal.II/fe/fe_poly_tensor.h b/include/deal.II/fe/fe_poly_tensor.h
index 78bb7dbfc1..fefaba3f58 100644
--- a/include/deal.II/fe/fe_poly_tensor.h
+++ b/include/deal.II/fe/fe_poly_tensor.h
@@ -188,16 +188,8 @@ protected:
     // generate a new data object and
     // initialize some fields
     InternalData *data = new InternalData;
+    data->update_each = update_each(update_flags) | update_once(update_flags);  // FIX: only update_each required
 
-    // check what needs to be
-    // initialized only once and what
-    // on every cell/face/subface we
-    // visit
-    data->update_once = update_once(update_flags);
-    data->update_each = update_each(update_flags);
-    data->update_flags = data->update_once | data->update_each;
-
-    const UpdateFlags flags(data->update_flags);
     const unsigned int n_q_points = quadrature.size();
 
     // some scratch arrays
@@ -207,13 +199,13 @@ protected:
     std::vector<Tensor<4,dim> > third_derivatives(0);
     std::vector<Tensor<5,dim> > fourth_derivatives(0);
 
-    if (flags & (update_values | update_gradients | update_hessians) )
+    if (update_flags & (update_values | update_gradients | update_hessians) )
       data->sign_change.resize (this->dofs_per_cell);
 
     // initialize fields only if really
     // necessary. otherwise, don't
     // allocate memory
-    if (flags & update_values)
+    if (update_flags & update_values)
       {
         values.resize (this->dofs_per_cell);
         data->shape_values.resize (this->dofs_per_cell);
@@ -223,7 +215,7 @@ protected:
           data->transformed_shape_values.resize (n_q_points);
       }
 
-    if (flags & update_gradients)
+    if (update_flags & update_gradients)
       {
         grads.resize (this->dofs_per_cell);
         data->shape_grads.resize (this->dofs_per_cell);
@@ -241,7 +233,7 @@ protected:
           data->untransformed_shape_grads.resize(n_q_points);
       }
 
-    if (flags & update_hessians)
+    if (update_flags & update_hessians)
       {
         grad_grads.resize (this->dofs_per_cell);
         data->shape_grad_grads.resize (this->dofs_per_cell);
@@ -259,7 +251,7 @@ protected:
     // node values N_i holds
     // N_i(v_j)=\delta_ij for all basis
     // functions v_j
-    if (flags & (update_values | update_gradients))
+    if (update_flags & (update_values | update_gradients))
       for (unsigned int k=0; k<n_q_points; ++k)
         {
           poly_space.compute(quadrature.point(k),
@@ -267,7 +259,7 @@ protected:
                              third_derivatives,
                              fourth_derivatives);
 
-          if (flags & update_values)
+          if (update_flags & update_values)
             {
               if (inverse_node_matrix.n_cols() == 0)
                 for (unsigned int i=0; i<this->dofs_per_cell; ++i)
@@ -282,7 +274,7 @@ protected:
                   }
             }
 
-          if (flags & update_gradients)
+          if (update_flags & update_gradients)
             {
               if (inverse_node_matrix.n_cols() == 0)
                 for (unsigned int i=0; i<this->dofs_per_cell; ++i)
@@ -297,7 +289,7 @@ protected:
                   }
             }
 
-          if (flags & update_hessians)
+          if (update_flags & update_hessians)
             {
               if (inverse_node_matrix.n_cols() == 0)
                 for (unsigned int i=0; i<this->dofs_per_cell; ++i)
diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h
index 999de4106a..ff256ef3a6 100644
--- a/include/deal.II/fe/fe_system.h
+++ b/include/deal.II/fe/fe_system.h
@@ -1145,15 +1145,6 @@ private:
     internal::FEValues::FiniteElementRelatedData<dim,spacedim> &
     get_fe_output_object (const unsigned int base_no) const;
 
-    /**
-     * Set the @p first_cell flag to @p false. Used by the @p FEValues class
-     * to indicate that we have already done the work on the first cell.
-     *
-     * In addition to calling the respective function of the base class, this
-     * function also calls the functions of the sub-data objects.
-     */
-    virtual void clear_first_cell ();
-
   private:
 
     /**
diff --git a/source/fe/fe.cc b/source/fe/fe.cc
index 4bf0d9cc12..8107ac7e5f 100644
--- a/source/fe/fe.cc
+++ b/source/fe/fe.cc
@@ -36,10 +36,7 @@ DEAL_II_NAMESPACE_OPEN
 
 template <int dim, int spacedim>
 FiniteElement<dim, spacedim>::InternalDataBase::InternalDataBase ():
-  update_flags(update_default),
-  update_once(update_default),
-  update_each(update_default),
-  first_cell(true)
+  update_each(update_default)
 {}
 
 
@@ -59,32 +56,6 @@ FiniteElement<dim, spacedim>::InternalDataBase::memory_consumption () const
 
 
 
-template <int dim, int spacedim>
-UpdateFlags
-FiniteElement<dim,spacedim>::InternalDataBase::current_update_flags () const
-{
-  if (first_cell)
-    {
-      Assert (update_flags==(update_once|update_each),
-              ExcInternalError());
-      return update_flags;
-    }
-  else
-    return update_each;
-}
-
-
-
-template <int dim, int spacedim>
-void
-FiniteElement<dim,spacedim>::InternalDataBase::clear_first_cell ()
-{
-  first_cell = false;
-}
-
-
-
-
 template <int dim, int spacedim>
 FiniteElement<dim,spacedim>::
 FiniteElement (const FiniteElementData<dim> &fe_data,
diff --git a/source/fe/fe_dgp_nonparametric.cc b/source/fe/fe_dgp_nonparametric.cc
index 4588b0c4b8..00045beebd 100644
--- a/source/fe/fe_dgp_nonparametric.cc
+++ b/source/fe/fe_dgp_nonparametric.cc
@@ -255,13 +255,7 @@ FE_DGPNonparametric<dim,spacedim>::get_data (
 {
   // generate a new data object
   typename FiniteElement<dim,spacedim>::InternalDataBase *data = new typename FiniteElement<dim,spacedim>::InternalDataBase;
-  // check what needs to be
-  // initialized only once and what
-  // on every cell/face/subface we
-  // visit
-  data->update_once = update_once(update_flags);
-  data->update_each = update_each(update_flags);
-  data->update_flags = data->update_once | data->update_each;
+  data->update_each = update_once(update_flags) | update_each(update_flags);   // FIX: only update_each required
 
   // other than that, there is nothing we can add here as discussed
   // in the general documentation of this class
@@ -287,18 +281,17 @@ fill_fe_values (const Mapping<dim,spacedim> &,
                 internal::FEValues::FiniteElementRelatedData<dim,spacedim> &output_data,
                 const CellSimilarity::Similarity /*cell_similarity*/) const
 {
-  const UpdateFlags flags(fe_data.current_update_flags());
-  Assert (flags & update_quadrature_points, ExcInternalError());
+  Assert (fe_data.update_each & update_quadrature_points, ExcInternalError());
 
   const unsigned int n_q_points = mapping_data.quadrature_points.size();
 
-  std::vector<double> values(flags & update_values ? this->dofs_per_cell : 0);
-  std::vector<Tensor<1,dim> > grads(flags & update_gradients ? this->dofs_per_cell : 0);
-  std::vector<Tensor<2,dim> > grad_grads(flags & update_hessians ? this->dofs_per_cell : 0);
+  std::vector<double> values(fe_data.update_each & update_values ? this->dofs_per_cell : 0);
+  std::vector<Tensor<1,dim> > grads(fe_data.update_each & update_gradients ? this->dofs_per_cell : 0);
+  std::vector<Tensor<2,dim> > grad_grads(fe_data.update_each & update_hessians ? this->dofs_per_cell : 0);
   std::vector<Tensor<3,dim> > empty_vector_of_3rd_order_tensors;
   std::vector<Tensor<4,dim> > empty_vector_of_4th_order_tensors;
 
-  if (flags & (update_values | update_gradients))
+  if (fe_data.update_each & (update_values | update_gradients))
     for (unsigned int i=0; i<n_q_points; ++i)
       {
         polynomial_space.compute(mapping_data.quadrature_points[i],
@@ -307,11 +300,11 @@ fill_fe_values (const Mapping<dim,spacedim> &,
                                  empty_vector_of_4th_order_tensors);
         for (unsigned int k=0; k<this->dofs_per_cell; ++k)
           {
-            if (flags & update_values)
+            if (fe_data.update_each & update_values)
               output_data.shape_values[k][i] = values[k];
-            if (flags & update_gradients)
+            if (fe_data.update_each & update_gradients)
               output_data.shape_gradients[k][i] = grads[k];
-            if (flags & update_hessians)
+            if (fe_data.update_each & update_hessians)
               output_data.shape_hessians[k][i] = grad_grads[k];
           }
       }
@@ -331,18 +324,17 @@ fill_fe_face_values (const Mapping<dim,spacedim> &,
                      const internal::FEValues::MappingRelatedData<dim,spacedim> &mapping_data,
                      internal::FEValues::FiniteElementRelatedData<dim,spacedim> &output_data) const
 {
-  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
-  Assert (flags & update_quadrature_points, ExcInternalError());
+  Assert (fe_data.update_each & update_quadrature_points, ExcInternalError());
 
   const unsigned int n_q_points = mapping_data.quadrature_points.size();
 
-  std::vector<double> values(flags & update_values ? this->dofs_per_cell : 0);
-  std::vector<Tensor<1,dim> > grads(flags & update_gradients ? this->dofs_per_cell : 0);
-  std::vector<Tensor<2,dim> > grad_grads(flags & update_hessians ? this->dofs_per_cell : 0);
+  std::vector<double> values(fe_data.update_each & update_values ? this->dofs_per_cell : 0);
+  std::vector<Tensor<1,dim> > grads(fe_data.update_each & update_gradients ? this->dofs_per_cell : 0);
+  std::vector<Tensor<2,dim> > grad_grads(fe_data.update_each & update_hessians ? this->dofs_per_cell : 0);
   std::vector<Tensor<3,dim> > empty_vector_of_3rd_order_tensors;
   std::vector<Tensor<4,dim> > empty_vector_of_4th_order_tensors;
 
-  if (flags & (update_values | update_gradients))
+  if (fe_data.update_each & (update_values | update_gradients))
     for (unsigned int i=0; i<n_q_points; ++i)
       {
         polynomial_space.compute(mapping_data.quadrature_points[i],
@@ -351,11 +343,11 @@ fill_fe_face_values (const Mapping<dim,spacedim> &,
                                  empty_vector_of_4th_order_tensors);
         for (unsigned int k=0; k<this->dofs_per_cell; ++k)
           {
-            if (flags & update_values)
+            if (fe_data.update_each & update_values)
               output_data.shape_values[k][i] = values[k];
-            if (flags & update_gradients)
+            if (fe_data.update_each & update_gradients)
               output_data.shape_gradients[k][i] = grads[k];
-            if (flags & update_hessians)
+            if (fe_data.update_each & update_hessians)
               output_data.shape_hessians[k][i] = grad_grads[k];
           }
       }
@@ -376,18 +368,17 @@ fill_fe_subface_values (const Mapping<dim,spacedim> &,
                         const internal::FEValues::MappingRelatedData<dim,spacedim> &mapping_data,
                         internal::FEValues::FiniteElementRelatedData<dim,spacedim> &output_data) const
 {
-  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
-  Assert (flags & update_quadrature_points, ExcInternalError());
+  Assert (fe_data.update_each & update_quadrature_points, ExcInternalError());
 
   const unsigned int n_q_points = mapping_data.quadrature_points.size();
 
-  std::vector<double> values(flags & update_values ? this->dofs_per_cell : 0);
-  std::vector<Tensor<1,dim> > grads(flags & update_gradients ? this->dofs_per_cell : 0);
-  std::vector<Tensor<2,dim> > grad_grads(flags & update_hessians ? this->dofs_per_cell : 0);
+  std::vector<double> values(fe_data.update_each & update_values ? this->dofs_per_cell : 0);
+  std::vector<Tensor<1,dim> > grads(fe_data.update_each & update_gradients ? this->dofs_per_cell : 0);
+  std::vector<Tensor<2,dim> > grad_grads(fe_data.update_each & update_hessians ? this->dofs_per_cell : 0);
   std::vector<Tensor<3,dim> > empty_vector_of_3rd_order_tensors;
   std::vector<Tensor<4,dim> > empty_vector_of_4th_order_tensors;
 
-  if (flags & (update_values | update_gradients))
+  if (fe_data.update_each & (update_values | update_gradients))
     for (unsigned int i=0; i<n_q_points; ++i)
       {
         polynomial_space.compute(mapping_data.quadrature_points[i],
@@ -396,11 +387,11 @@ fill_fe_subface_values (const Mapping<dim,spacedim> &,
                                  empty_vector_of_4th_order_tensors);
         for (unsigned int k=0; k<this->dofs_per_cell; ++k)
           {
-            if (flags & update_values)
+            if (fe_data.update_each & update_values)
               output_data.shape_values[k][i] = values[k];
-            if (flags & update_gradients)
+            if (fe_data.update_each & update_gradients)
               output_data.shape_gradients[k][i] = grads[k];
-            if (flags & update_hessians)
+            if (fe_data.update_each & update_hessians)
               output_data.shape_hessians[k][i] = grad_grads[k];
           }
       }
diff --git a/source/fe/fe_face.cc b/source/fe/fe_face.cc
index 9732fab950..7f3c8b9458 100644
--- a/source/fe/fe_face.cc
+++ b/source/fe/fe_face.cc
@@ -494,10 +494,8 @@ fill_fe_face_values (const Mapping<1,spacedim> &,
                      const internal::FEValues::MappingRelatedData<1,spacedim> &,
                      internal::FEValues::FiniteElementRelatedData<1,spacedim> &output_data) const
 {
-  const UpdateFlags flags(fedata.update_once | fedata.update_each);
-
   const unsigned int foffset = face;
-  if (flags & update_values)
+  if (fedata.update_each & update_values)
     {
       for (unsigned int k=0; k<this->dofs_per_cell; ++k)
         output_data.shape_values(k,0) = 0.;
diff --git a/source/fe/fe_poly.cc b/source/fe/fe_poly.cc
index 7981da4946..c4c49b6ddd 100644
--- a/source/fe/fe_poly.cc
+++ b/source/fe/fe_poly.cc
@@ -47,21 +47,19 @@ fill_fe_values (const Mapping<1,2>                                &mapping,
   Assert (dynamic_cast<const InternalData *> (&fedata) != 0, ExcInternalError());
   const InternalData &fe_data = static_cast<const InternalData &> (fedata);
 
-  const UpdateFlags flags(fe_data.current_update_flags());
-
   for (unsigned int k=0; k<this->dofs_per_cell; ++k)
     {
-      if (flags & update_values)
+      if (fe_data.update_each & update_values)
         for (unsigned int i=0; i<quadrature.size(); ++i)
           output_data.shape_values(k,i) = fe_data.shape_values[k][i];
 
-      if (flags & update_gradients && cell_similarity != CellSimilarity::translation)
+      if (fe_data.update_each & update_gradients && cell_similarity != CellSimilarity::translation)
         mapping.transform (fe_data.shape_gradients[k],
                            mapping_covariant,
                            mapping_internal,
                            output_data.shape_gradients[k]);
 
-      if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+      if (fe_data.update_each & update_hessians && cell_similarity != CellSimilarity::translation)
         {
           mapping.transform (fe_data.shape_hessians[k],
                              mapping_covariant_gradient,
@@ -75,7 +73,7 @@ fill_fe_values (const Mapping<1,2>                                &mapping,
                 * output_data.shape_gradients[k][i][j];
         }
 
-      if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+      if (fe_data.update_each & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
         {
           mapping.transform (fe_data.shape_3rd_derivatives[k],
                              mapping_covariant_hessian,
@@ -107,21 +105,19 @@ fill_fe_values (const Mapping<2,3>                                &mapping,
   Assert (dynamic_cast<const InternalData *> (&fedata) != 0, ExcInternalError());
   const InternalData &fe_data = static_cast<const InternalData &> (fedata);
 
-  const UpdateFlags flags(fe_data.current_update_flags());
-
   for (unsigned int k=0; k<this->dofs_per_cell; ++k)
     {
-      if (flags & update_values)
+      if (fe_data.update_each & update_values)
         for (unsigned int i=0; i<quadrature.size(); ++i)
           output_data.shape_values(k,i) = fe_data.shape_values[k][i];
 
-      if (flags & update_gradients && cell_similarity != CellSimilarity::translation)
+      if (fe_data.update_each & update_gradients && cell_similarity != CellSimilarity::translation)
         mapping.transform (fe_data.shape_gradients[k],
                            mapping_covariant,
                            mapping_internal,
                            output_data.shape_gradients[k]);
 
-      if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+      if (fe_data.update_each & update_hessians && cell_similarity != CellSimilarity::translation)
         {
           mapping.transform (fe_data.shape_hessians[k],
                              mapping_covariant_gradient,
@@ -135,7 +131,7 @@ fill_fe_values (const Mapping<2,3>                                &mapping,
                 * output_data.shape_gradients[k][i][j];
         }
 
-      if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+      if (fe_data.update_each & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
         {
           mapping.transform (fe_data.shape_3rd_derivatives[k],
                              mapping_covariant_hessian,
@@ -168,21 +164,19 @@ fill_fe_values (const Mapping<1,2>                                &mapping,
   Assert (dynamic_cast<const InternalData *> (&fedata) != 0, ExcInternalError());
   const InternalData &fe_data = static_cast<const InternalData &> (fedata);
 
-  const UpdateFlags flags(fe_data.current_update_flags());
-
   for (unsigned int k=0; k<this->dofs_per_cell; ++k)
     {
-      if (flags & update_values)
+      if (fe_data.update_each & update_values)
         for (unsigned int i=0; i<quadrature.size(); ++i)
           output_data.shape_values(k,i) = fe_data.shape_values[k][i];
 
-      if (flags & update_gradients && cell_similarity != CellSimilarity::translation)
+      if (fe_data.update_each & update_gradients && cell_similarity != CellSimilarity::translation)
         mapping.transform (fe_data.shape_gradients[k],
                            mapping_covariant,
                            mapping_internal,
                            output_data.shape_gradients[k]);
 
-      if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+      if (fe_data.update_each & update_hessians && cell_similarity != CellSimilarity::translation)
         {
           mapping.transform (fe_data.shape_hessians[k],
                              mapping_covariant_gradient,
@@ -196,7 +190,7 @@ fill_fe_values (const Mapping<1,2>                                &mapping,
                 * output_data.shape_gradients[k][i][j];
         }
 
-      if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+      if (fe_data.update_each & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
         {
           mapping.transform (fe_data.shape_3rd_derivatives[k],
                              mapping_covariant_hessian,
@@ -224,22 +218,20 @@ fill_fe_values (const Mapping<2,3>                                &mapping,
   Assert (dynamic_cast<const InternalData *> (&fedata) != 0, ExcInternalError());
   const InternalData &fe_data = static_cast<const InternalData &> (fedata);
 
-  const UpdateFlags flags(fe_data.current_update_flags());
-
   for (unsigned int k=0; k<this->dofs_per_cell; ++k)
     {
-      if (flags & update_values)
+      if (fe_data.update_each & update_values)
         for (unsigned int i=0; i<quadrature.size(); ++i)
           output_data.shape_values(k,i) = fe_data.shape_values[k][i];
 
 
-      if (flags & update_gradients && cell_similarity != CellSimilarity::translation)
+      if (fe_data.update_each & update_gradients && cell_similarity != CellSimilarity::translation)
         mapping.transform (fe_data.shape_gradients[k],
                            mapping_covariant,
                            mapping_internal,
                            output_data.shape_gradients[k]);
 
-      if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+      if (fe_data.update_each & update_hessians && cell_similarity != CellSimilarity::translation)
         {
           mapping.transform (fe_data.shape_hessians[k],
                              mapping_covariant_gradient,
@@ -253,7 +245,7 @@ fill_fe_values (const Mapping<2,3>                                &mapping,
                 * output_data.shape_gradients[k][i][j];
         }
 
-      if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+      if (fe_data.update_each & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
         {
           mapping.transform (fe_data.shape_3rd_derivatives[k],
                              mapping_covariant_hessian,
diff --git a/source/fe/fe_poly_tensor.cc b/source/fe/fe_poly_tensor.cc
index 8939e9d05d..cd116087e8 100644
--- a/source/fe/fe_poly_tensor.cc
+++ b/source/fe/fe_poly_tensor.cc
@@ -307,11 +307,10 @@ fill_fe_values (const Mapping<dim,spacedim>                                  &ma
   const InternalData &fe_data = static_cast<const InternalData &> (fedata);
 
   const unsigned int n_q_points = quadrature.size();
-  const UpdateFlags flags(fe_data.current_update_flags());
 
-  Assert(!(flags & update_values) || fe_data.shape_values.size() == this->dofs_per_cell,
+  Assert(!(fe_data.update_each & update_values) || fe_data.shape_values.size() == this->dofs_per_cell,
          ExcDimensionMismatch(fe_data.shape_values.size(), this->dofs_per_cell));
-  Assert(!(flags & update_values) || fe_data.shape_values[0].size() == n_q_points,
+  Assert(!(fe_data.update_each & update_values) || fe_data.shape_values[0].size() == n_q_points,
          ExcDimensionMismatch(fe_data.shape_values[0].size(), n_q_points));
 
   // Create table with sign changes, due to the special structure of the RT elements.
@@ -338,7 +337,7 @@ fill_fe_values (const Mapping<dim,spacedim>                                  &ma
       // the previous one; or, even if it is a translation, if we use mappings
       // other than the standard mappings that require us to recompute values
       // and derivatives because of possible sign changes
-      if (flags & update_values &&
+      if (fe_data.update_each & update_values &&
           ((cell_similarity != CellSimilarity::translation)
            ||
            ((mapping_type == mapping_piola) || (mapping_type == mapping_raviart_thomas)
@@ -404,7 +403,7 @@ fill_fe_values (const Mapping<dim,spacedim>                                  &ma
         }
 
       // update gradients. apply the same logic as above
-      if (flags & update_gradients
+      if (fe_data.update_each & update_gradients
           &&
           ((cell_similarity != CellSimilarity::translation)
            ||
@@ -532,7 +531,7 @@ fill_fe_values (const Mapping<dim,spacedim>                                  &ma
         }
 
       // update hessians. apply the same logic as above
-      if (flags & update_hessians
+      if (fe_data.update_each & update_hessians
           &&
           ((cell_similarity != CellSimilarity::translation)
            ||
@@ -742,7 +741,7 @@ fill_fe_values (const Mapping<dim,spacedim>                                  &ma
         }
 
       // third derivatives are not implemented
-      if (flags & update_3rd_derivatives
+      if (fe_data.update_each & update_3rd_derivatives
           &&
           ((cell_similarity != CellSimilarity::translation)
            ||
@@ -788,8 +787,6 @@ fill_fe_face_values (const Mapping<dim,spacedim>
                                                 cell->face_rotation(face),
                                                 n_q_points);
 
-  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
-
 //TODO: Size assertions
 
 // Create table with sign changes, due to the special structure of the RT elements.
@@ -810,7 +807,7 @@ fill_fe_face_values (const Mapping<dim,spacedim>
       const unsigned int first = output_data.shape_function_to_row_table[i * this->n_components() +
                                  this->get_nonzero_components(i).first_selected_component()];
 
-      if (flags & update_values)
+      if (fe_data.update_each & update_values)
         {
           switch (mapping_type)
             {
@@ -880,7 +877,7 @@ fill_fe_face_values (const Mapping<dim,spacedim>
             }
         }
 
-      if (flags & update_gradients)
+      if (fe_data.update_each & update_gradients)
         {
           VectorSlice< std::vector<Tensor<2,spacedim> > >
           transformed_shape_grads (fe_data.transformed_shape_grads, offset, n_q_points);
@@ -1004,7 +1001,7 @@ fill_fe_face_values (const Mapping<dim,spacedim>
             }
         }
 
-      if (flags & update_hessians)
+      if (fe_data.update_each & update_hessians)
 
         {
 
@@ -1208,7 +1205,7 @@ fill_fe_face_values (const Mapping<dim,spacedim>
         }
 
       // third derivatives are not implemented
-      if (flags & update_3rd_derivatives)
+      if (fe_data.update_each & update_3rd_derivatives)
         {
           Assert(false, ExcNotImplemented())
         }
@@ -1251,8 +1248,6 @@ fill_fe_subface_values (const Mapping<dim,spacedim>
                                                    n_q_points,
                                                    cell->subface_case(face));
 
-  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
-
 //   Assert(mapping_type == independent
 //       || ( mapping_type == independent_on_cartesian
 //            && dynamic_cast<const MappingCartesian<dim>*>(&mapping) != 0),
@@ -1276,7 +1271,7 @@ fill_fe_subface_values (const Mapping<dim,spacedim>
       const unsigned int first = output_data.shape_function_to_row_table[i * this->n_components() +
                                  this->get_nonzero_components(i).first_selected_component()];
 
-      if (flags & update_values)
+      if (fe_data.update_each & update_values)
         {
           switch (mapping_type)
             {
@@ -1343,7 +1338,7 @@ fill_fe_subface_values (const Mapping<dim,spacedim>
             }
         }
 
-      if (flags & update_gradients)
+      if (fe_data.update_each & update_gradients)
         {
           VectorSlice< std::vector<Tensor<2, spacedim > > >
           transformed_shape_grads (fe_data.transformed_shape_grads, offset, n_q_points);
@@ -1469,7 +1464,7 @@ fill_fe_subface_values (const Mapping<dim,spacedim>
             }
         }
 
-      if (flags & update_hessians)
+      if (fe_data.update_each & update_hessians)
 
         {
 
@@ -1672,7 +1667,7 @@ fill_fe_subface_values (const Mapping<dim,spacedim>
         }
 
       // third derivatives are not implemented
-      if (flags & update_3rd_derivatives)
+      if (fe_data.update_each & update_3rd_derivatives)
         {
           Assert(false, ExcNotImplemented())
         }
diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc
index afd837eca0..95f097d156 100644
--- a/source/fe/fe_system.cc
+++ b/source/fe/fe_system.cc
@@ -104,20 +104,6 @@ InternalData::get_fe_output_object (const unsigned int base_no) const
 
 
 
-template <int dim, int spacedim>
-void
-FESystem<dim,spacedim>::InternalData::clear_first_cell ()
-{
-  // call respective function of base
-  // class
-  FiniteElement<dim,spacedim>::InternalDataBase::clear_first_cell ();
-  // then the functions of all the
-  // sub-objects
-  for (unsigned int i=0; i<base_fe_datas.size(); ++i)
-    base_fe_datas[i]->clear_first_cell ();
-}
-
-
 /* ---------------------------------- FESystem ------------------- */
 
 
@@ -893,10 +879,7 @@ FESystem<dim,spacedim>::get_data (const UpdateFlags      flags_,
 {
   UpdateFlags flags = flags_;
   InternalData *data = new InternalData(this->n_base_elements());
-
-  data->update_once = update_once (flags);
-  data->update_each = update_each (flags);
-  flags = data->update_once | data->update_each;
+  data->update_each = update_each(flags) | update_once(flags);
 
   // get data objects from each of the base elements and store
   // them. do the creation of these objects in parallel as their
@@ -922,13 +905,11 @@ FESystem<dim,spacedim>::get_data (const UpdateFlags      flags_,
       internal::FEValues::FiniteElementRelatedData<dim,spacedim> &base_fe_output_object
         = data->get_fe_output_object(base_no);
       base_fe_output_object.initialize (quadrature.size(), base_element(base_no),
-                                        flags | base_fe_data->update_flags);
+                                        flags | base_fe_data->update_each);
 
       // then store the pointer to the base internal object
       data->set_fe_data(base_no, base_fe_data);
     }
-  data->update_flags = data->update_once |
-                       data->update_each;
   return data;
 }
 
@@ -944,10 +925,7 @@ FESystem<dim,spacedim>::get_face_data (
 {
   UpdateFlags flags = flags_;
   InternalData *data = new InternalData(this->n_base_elements());
-
-  data->update_once = update_once (flags);
-  data->update_each = update_each (flags);
-  flags = data->update_once | data->update_each;
+  data->update_each = update_each(flags) | update_once(flags);
 
   // get data objects from each of the base elements and store
   // them. do the creation of these objects in parallel as their
@@ -973,13 +951,11 @@ FESystem<dim,spacedim>::get_face_data (
       internal::FEValues::FiniteElementRelatedData<dim,spacedim> &base_fe_output_object
         = data->get_fe_output_object(base_no);
       base_fe_output_object.initialize (quadrature.size(), base_element(base_no),
-                                        flags | base_fe_data->update_flags);
+                                        flags | base_fe_data->update_each);
 
       // then store the pointer to the base internal object
       data->set_fe_data(base_no, base_fe_data);
     }
-  data->update_flags = data->update_once |
-                       data->update_each;
   return data;
 }
 
@@ -997,10 +973,7 @@ FESystem<dim,spacedim>::get_subface_data (
 {
   UpdateFlags flags = flags_;
   InternalData *data = new InternalData(this->n_base_elements());
-
-  data->update_once = update_once (flags);
-  data->update_each = update_each (flags);
-  flags = data->update_once | data->update_each;
+  data->update_each = update_each(flags) | update_once(flags);
 
   // get data objects from each of the base elements and store
   // them. do the creation of these objects in parallel as their
@@ -1026,13 +999,11 @@ FESystem<dim,spacedim>::get_subface_data (
       internal::FEValues::FiniteElementRelatedData<dim,spacedim> &base_fe_output_object
         = data->get_fe_output_object(base_no);
       base_fe_output_object.initialize (quadrature.size(), base_element(base_no),
-                                        flags | base_fe_data->update_flags);
+                                        flags | base_fe_data->update_each);
 
       // then store the pointer to the base internal object
       data->set_fe_data(base_no, base_fe_data);
     }
-  data->update_flags = data->update_once |
-                       data->update_each;
   return data;
 }
 
@@ -1183,9 +1154,7 @@ compute_fill_one_base (const Mapping<dim,spacedim>                      &mapping
   // all shape functions of the composed element, and here only treat
   // those shape functions that belong to a given base element
   //TODO: Introduce the needed table and loop only over base element shape functions. This here is not efficient at all AND very bad style
-  const UpdateFlags base_flags(dim_1==dim ?
-                               base_fe_data.current_update_flags() :
-                               base_fe_data.update_flags);
+  const UpdateFlags base_flags = base_fe_data.update_each;
 
   // if the current cell is just a translation of the previous one,
   // the underlying data has not changed, and we don't even need to
@@ -1270,9 +1239,7 @@ compute_fill (const Mapping<dim,spacedim>                      &mapping,
   // (fill_fe_values) or dim_1==dim-1
   // (fill_fe_(sub)face_values)
   Assert(dim_1==dim || dim_1==dim-1, ExcInternalError());
-  const UpdateFlags flags(dim_1==dim ?
-                          fe_data.current_update_flags() :
-                          fe_data.update_flags);
+  const UpdateFlags flags = fe_data.update_each;
 
 
   if (flags & (update_values | update_gradients
diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc
index f5341ebee6..2457767b9b 100644
--- a/source/fe/fe_values.cc
+++ b/source/fe/fe_values.cc
@@ -3768,8 +3768,6 @@ void FEValues<dim,spacedim>::do_reinit ()
                                 this->mapping_output,
                                 this->finite_element_output,
                                 this->cell_similarity);
-
-  this->fe_data->clear_first_cell ();
 }
 
 
@@ -3974,8 +3972,6 @@ void FEFaceValues<dim,spacedim>::do_reinit (const unsigned int face_no)
                                      *this->fe_data,
                                      this->mapping_output,
                                      this->finite_element_output);
-
-  this->fe_data->clear_first_cell ();
 }
 
 
@@ -4217,8 +4213,6 @@ void FESubfaceValues<dim,spacedim>::do_reinit (const unsigned int face_no,
                                         *this->fe_data,
                                         this->mapping_output,
                                         this->finite_element_output);
-
-  this->fe_data->clear_first_cell ();
 }
 
 
-- 
2.39.5