Do some work already in get_data().

Do not delay until fill_fe_values(), if something can already be done early on.

diff --git a/include/deal.II/fe/fe_poly.h b/include/deal.II/fe/fe_poly.h
index 134b5ef3a06c5c9d900b30381a4b7414f655c5f3..01db90689f36c714c13451b58ad77f297182547d 100644 (file)
--- a/include/deal.II/fe/fe_poly.h
+++ b/include/deal.II/fe/fe_poly.h
@@ -223,7 +223,7 @@ protected:
   get_data(const UpdateFlags                                                    update_flags,
            const Mapping<dim,spacedim>                                         &/*mapping*/,
            const Quadrature<dim>                                               &quadrature,
-           dealii::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &/*output_data*/) const
+           dealii::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &output_data) const
   {
     // generate a new data object and
     // initialize some fields
@@ -232,51 +232,58 @@ protected:
 
     const unsigned int n_q_points = quadrature.size();
 
-    // some scratch arrays
-    std::vector<double> values(0);
-    std::vector<Tensor<1,dim> > grads(0);
-    std::vector<Tensor<2,dim> > grad_grads(0);
-    std::vector<Tensor<3,dim> > third_derivatives(0);
-    std::vector<Tensor<4,dim> > fourth_derivatives(0);
-
-    // initialize fields only if really
-    // necessary. otherwise, don't
-    // allocate memory
-    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));
-      }
+    // initialize some scratch arrays. we need them for the underlying
+    // polynomial to put the values and derivatives of shape functions
+    // to put there, depending on what the user requested
+    std::vector<double> values(update_flags & update_values ?
+                               this->dofs_per_cell : 0);
+    std::vector<Tensor<1,dim> > grads(update_flags & update_gradients ?
+                                      this->dofs_per_cell : 0);
+    std::vector<Tensor<2,dim> > grad_grads(update_flags & update_hessians ?
+                                           this->dofs_per_cell : 0);
+    std::vector<Tensor<3,dim> > third_derivatives(update_flags & update_3rd_derivatives ?
+                                                  this->dofs_per_cell : 0);
+    std::vector<Tensor<4,dim> > fourth_derivatives;   // won't be needed, so leave empty
+
+    // now also initialize fields the fields of this class's own
+    // temporary storage, depending on what we need for the given
+    // update flags.
+    //
+    // there is one exception from the rule: if we are dealing with
+    // cells (i.e., if this function is not called via
+    // get_(sub)face_data()), then we can already store things in the
+    // final location where FEValues::reinit() later wants to see
+    // things. we then don't need the intermediate space. we determine
+    // whether we are on a cell by asking whether the number of
+    // elements in the output array equals the number of quadrature
+    // points (yes, it's a cell) or not (because in that case the
+    // number of quadrature points we use here equals the number of
+    // quadrature points summed over *all* faces or subfaces, whereas
+    // the number of output slots equals the number of quadrature
+    // points on only *one* face)
+    if ((update_flags & update_values)
+        &&
+        !((output_data.shape_values.n_rows() > 0)
+          &&
+          (output_data.shape_values.n_cols() == n_q_points)))
+      data->shape_values.resize (this->dofs_per_cell,
+                                 std::vector<double> (n_q_points));
 
     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));
-      }
+      data->shape_gradients.resize (this->dofs_per_cell,
+                                    std::vector<Tensor<1,dim> > (n_q_points));
 
     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));
-      }
+      data->shape_hessians.resize (this->dofs_per_cell,
+                                   std::vector<Tensor<2,dim> > (n_q_points));
 
     if (update_flags & update_3rd_derivatives)
-      {
-        third_derivatives.resize (this->dofs_per_cell);
-        data->shape_3rd_derivatives.resize (this->dofs_per_cell,
-                                            std::vector<Tensor<3,dim> > (n_q_points));
-      }
-
-    // next already fill those fields
-    // of which we have information by
-    // now. note that the shape
-    // gradients are only those on the
-    // unit cell, and need to be
-    // transformed when visiting an
-    // actual cell
+      data->shape_3rd_derivatives.resize (this->dofs_per_cell,
+                                          std::vector<Tensor<3,dim> > (n_q_points));
+
+    // next already fill those fields of which we have information by
+    // now. note that the shape gradients are only those on the unit
+    // cell, and need to be transformed when visiting an actual cell
     if (update_flags & (update_values | update_gradients
                         | update_hessians | update_3rd_derivatives) )
       for (unsigned int i=0; i<n_q_points; ++i)
@@ -286,10 +293,27 @@ protected:
                              third_derivatives,
                              fourth_derivatives);
 
+          // the values of shape functions at quadrature points don't change.
+          // consequently, write these values right into the output array if
+          // we can, i.e., if the output array has the correct size. this is
+          // the case on cells. on faces, we already precompute data on *all*
+          // faces and subfaces, but we later on copy only a portion of it
+          // into the output object; in that case, copy the data from all
+          // faces into the scratch object
           if (update_flags & update_values)
-            for (unsigned int k=0; k<this->dofs_per_cell; ++k)
-              data->shape_values[k][i] = values[k];
-
+            if (output_data.shape_values.n_rows() > 0)
+              {
+                if (output_data.shape_values.n_cols() == n_q_points)
+                  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+                    output_data.shape_values[k][i] = values[k];
+                else
+                  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+                    data->shape_values[k][i] = values[k];
+              }
+
+          // for everything else, derivatives need to be transformed,
+          // so we write them into our scratch space and only later
+          // copy stuff into where FEValues wants it
           if (update_flags & update_gradients)
             for (unsigned int k=0; k<this->dofs_per_cell; ++k)
               data->shape_gradients[k][i] = grads[k];

diff --git a/include/deal.II/fe/fe_poly.templates.h b/include/deal.II/fe/fe_poly.templates.h
index c55bf59053b5417b78b6c139c39b03cb22bdfb41..bc2bf1d9d6e075c0edb5c4817f18480e91e60e8f 100644 (file)
--- a/include/deal.II/fe/fe_poly.templates.h
+++ b/include/deal.II/fe/fe_poly.templates.h
@@ -254,11 +254,9 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &,
 
   const UpdateFlags flags(fe_data.update_each);
 
-  if (flags & update_values)
-    for (unsigned int k=0; k<this->dofs_per_cell; ++k)
-      for (unsigned int i=0; i<quadrature.size(); ++i)
-        output_data.shape_values(k,i) = fe_data.shape_values[k][i];
-
+  // transform gradients and higher derivatives. there is nothing to do
+  // for values since we already emplaced them into output_data when
+  // we were in get_data()
   if (flags & update_gradients && cell_similarity != CellSimilarity::translation)
     for (unsigned int k=0; k<this->dofs_per_cell; ++k)
       mapping.transform (fe_data.shape_gradients[k],
@@ -329,6 +327,9 @@ fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator
 
   const UpdateFlags flags(fe_data.update_each);
 
+  // transform gradients and higher derivatives. we also have to copy
+  // the values (unlike in the case of fill_fe_values()) since
+  // we need to take into account the offsets
   if (flags & update_values)
     for (unsigned int k=0; k<this->dofs_per_cell; ++k)
       for (unsigned int i=0; i<quadrature.size(); ++i)
@@ -408,6 +409,9 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato
 
   const UpdateFlags flags(fe_data.update_each);
 
+  // transform gradients and higher derivatives. we also have to copy
+  // the values (unlike in the case of fill_fe_values()) since
+  // we need to take into account the offsets
   if (flags & update_values)
     for (unsigned int k=0; k<this->dofs_per_cell; ++k)
       for (unsigned int i=0; i<quadrature.size(); ++i)

diff --git a/source/fe/fe_poly.cc b/source/fe/fe_poly.cc
index bd52bab07826f5ffe6550bfa35bedf4a00cddf41..583c23a5fc671e7eb6a48d83a9577a5aa0d24ede 100644 (file)
--- a/source/fe/fe_poly.cc
+++ b/source/fe/fe_poly.cc
@@ -49,10 +49,9 @@ fill_fe_values (const Triangulation<1,2>::cell_iterator &,
 
   for (unsigned int k=0; k<this->dofs_per_cell; ++k)
     {
-      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];
-
+      // transform gradients and higher derivatives. there is nothing to do
+      // for values since we already emplaced them into output_data when
+      // we were in get_data()
       if (fe_data.update_each & update_gradients && cell_similarity != CellSimilarity::translation)
         mapping.transform (fe_data.shape_gradients[k],
                            mapping_covariant,
@@ -107,10 +106,9 @@ fill_fe_values (const Triangulation<2,3>::cell_iterator &,
 
   for (unsigned int k=0; k<this->dofs_per_cell; ++k)
     {
-      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];
-
+      // transform gradients and higher derivatives. there is nothing to do
+      // for values since we already emplaced them into output_data when
+      // we were in get_data()
       if (fe_data.update_each & update_gradients && cell_similarity != CellSimilarity::translation)
         mapping.transform (fe_data.shape_gradients[k],
                            mapping_covariant,
@@ -166,10 +164,9 @@ fill_fe_values (const Triangulation<1,2>::cell_iterator &,
 
   for (unsigned int k=0; k<this->dofs_per_cell; ++k)
     {
-      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];
-
+      // transform gradients and higher derivatives. there is nothing to do
+      // for values since we already emplaced them into output_data when
+      // we were in get_data()
       if (fe_data.update_each & update_gradients && cell_similarity != CellSimilarity::translation)
         mapping.transform (fe_data.shape_gradients[k],
                            mapping_covariant,
@@ -220,11 +217,9 @@ fill_fe_values (const Triangulation<2,3>::cell_iterator &,
 
   for (unsigned int k=0; k<this->dofs_per_cell; ++k)
     {
-      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];
-
-
+      // transform gradients and higher derivatives. there is nothing to do
+      // for values since we already emplaced them into output_data when
+      // we were in get_data()
       if (fe_data.update_each & update_gradients && cell_similarity != CellSimilarity::translation)
         mapping.transform (fe_data.shape_gradients[k],
                            mapping_covariant,