From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Fri, 6 Nov 2015 16:02:01 +0000 (-0600)
Subject: Re-inline a function previously broken out into its own context.
X-Git-Tag: v8.4.0-rc2~243^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F1841%2Fhead;p=dealii.git

Re-inline a function previously broken out into its own context.
---

diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h
index 5547211bd2..ace5ac4e14 100644
--- a/include/deal.II/fe/fe_system.h
+++ b/include/deal.II/fe/fe_system.h
@@ -1089,24 +1089,6 @@ private:
   std::vector<std::pair<unsigned int, unsigned int> >
   hp_object_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;
 
-  /**
-   * Compute the equivalent of compute_fill(), but only for the base element
-   * specified by the second-to-last argument. Some elements are grouped
-   * together to stay within the limit of 8 function arguments of boost.
-   */
-  template <int dim_1>
-  void
-  compute_fill_one_base (const Mapping<dim,spacedim>                      &mapping,
-                         const std::pair<typename Triangulation<dim,spacedim>::cell_iterator,
-                         CellSimilarity::Similarity>       cell_and_similarity,
-                         const std::pair<unsigned int,unsigned int>        face_sub_no,
-                         const Quadrature<dim_1>                          &quadrature,
-                         const std::pair<const typename Mapping<dim,spacedim>::InternalDataBase *,
-                         const typename FiniteElement<dim,spacedim>::InternalDataBase *> mapping_and_fe_internal,
-                         const unsigned int                                base_element,
-                         const internal::FEValues::MappingRelatedData<dim,spacedim> &mapping_data,
-                         internal::FEValues::FiniteElementRelatedData<dim,spacedim> &output_data) const;
-
   /**
    * Usually: Fields of cell-independent data.
    *
diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc
index 4570406cfa..0dd710f662 100644
--- a/source/fe/fe_system.cc
+++ b/source/fe/fe_system.cc
@@ -1057,151 +1057,6 @@ fill_fe_subface_values (const Mapping<dim,spacedim>                      &mappin
 
 
 
-template <int dim, int spacedim>
-template <int dim_1>
-void
-FESystem<dim,spacedim>::
-compute_fill_one_base (const Mapping<dim,spacedim>                      &mapping,
-                       const std::pair<typename Triangulation<dim,spacedim>::cell_iterator,
-                       CellSimilarity::Similarity>       cell_and_similarity,
-                       const std::pair<unsigned int, unsigned int>       face_sub_no,
-                       const Quadrature<dim_1>                          &quadrature,
-                       const std::pair<const typename Mapping<dim,spacedim>::InternalDataBase *,
-                       const typename FiniteElement<dim,spacedim>::InternalDataBase *> mapping_and_fe_internal,
-                       const unsigned int                                base_no,
-                       const internal::FEValues::MappingRelatedData<dim,spacedim> &mapping_data,
-                       internal::FEValues::FiniteElementRelatedData<dim,spacedim> &output_data) const
-{
-  const typename Triangulation<dim,spacedim>::cell_iterator cell            = cell_and_similarity.first;
-  const CellSimilarity::Similarity                          cell_similarity = cell_and_similarity.second;
-
-  const InternalData &fe_data = static_cast<const InternalData &> (*mapping_and_fe_internal.second);
-  const unsigned int n_q_points = quadrature.size();
-
-  const FiniteElement<dim,spacedim> &
-  base_fe      = base_element(base_no);
-  typename FiniteElement<dim,spacedim>::InternalDataBase &
-  base_fe_data = fe_data.get_fe_data(base_no);
-  internal::FEValues::FiniteElementRelatedData<dim,spacedim> &
-  base_data    = fe_data.get_fe_output_object(base_no);
-
-  // fill_fe_face_values needs argument Quadrature<dim-1> for both cases
-  // dim_1==dim-1 and dim_1=dim. Hence the following workaround
-  const Quadrature<dim>   *cell_quadrature = 0;
-  const Quadrature<dim-1> *face_quadrature = 0;
-
-  // static cast to the common base class of quadrature being either
-  // Quadrature<dim> or Quadrature<dim-1>:
-  const Subscriptor *quadrature_base_pointer = &quadrature;
-
-  if (face_sub_no.first==invalid_face_number)
-    {
-      Assert(dim_1==dim, ExcDimensionMismatch(dim_1,dim));
-      Assert (dynamic_cast<const Quadrature<dim> *>(quadrature_base_pointer) != 0,
-              ExcInternalError());
-
-      cell_quadrature
-        = static_cast<const Quadrature<dim> *>(quadrature_base_pointer);
-    }
-  else
-    {
-      Assert(dim_1==dim-1, ExcDimensionMismatch(dim_1,dim-1));
-      Assert (dynamic_cast<const Quadrature<dim-1> *>(quadrature_base_pointer) != 0,
-              ExcInternalError());
-
-      face_quadrature
-        = static_cast<const Quadrature<dim-1> *>(quadrature_base_pointer);
-    }
-
-
-  // Make sure that in the case of fill_fe_values the data is only
-  // copied from base_data to data if base_data is changed. therefore
-  // use fe_fe_data.current_update_flags()
-  //
-  // for the case of fill_fe_(sub)face_values the data needs to be
-  // copied from base_data to data on each face, therefore use
-  // base_fe_data.update_flags.
-  if (face_sub_no.first==invalid_face_number)
-    base_fe.fill_fe_values(mapping, cell, *cell_quadrature, *mapping_and_fe_internal.first,
-                           base_fe_data, mapping_data, base_data, cell_similarity);
-  else if (face_sub_no.second==invalid_face_number)
-    base_fe.fill_fe_face_values(mapping, cell, face_sub_no.first,
-                                *face_quadrature, *mapping_and_fe_internal.first, base_fe_data, mapping_data, base_data);
-  else
-    base_fe.fill_fe_subface_values(mapping, cell, face_sub_no.first, face_sub_no.second,
-                                   *face_quadrature, *mapping_and_fe_internal.first, base_fe_data, mapping_data, base_data);
-
-  // now data has been generated, so copy it. we used to work by
-  // looping over all base elements (i.e. this outer loop), then over
-  // multiplicity, then over the shape functions from that base
-  // element, but that requires that we can infer the global number of
-  // a shape function from its number in the base element. for that we
-  // had the component_to_system_table.
-  //
-  // however, this does of course no longer work since we have
-  // non-primitive elements. so we go the other way round: loop over
-  // 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 = 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
-  // enter this section
-  if (cell_similarity != CellSimilarity::translation)
-    for (unsigned int system_index=0; system_index<this->dofs_per_cell;
-         ++system_index)
-      if (this->system_to_base_table[system_index].first.first == base_no)
-        {
-          const unsigned int
-          base_index = this->system_to_base_table[system_index].second;
-          Assert (base_index<base_fe.dofs_per_cell, ExcInternalError());
-
-          // now copy. if the shape function is primitive, then there
-          // is only one value to be copied, but for non-primitive
-          // elements, there might be more values to be copied
-          //
-          // so, find out from which index to take this one value, and
-          // to which index to put
-          unsigned int out_index = 0;
-          for (unsigned int i=0; i<system_index; ++i)
-            out_index += this->n_nonzero_components(i);
-          unsigned int in_index = 0;
-          for (unsigned int i=0; i<base_index; ++i)
-            in_index += base_fe.n_nonzero_components(i);
-
-          // then loop over the number of components to be copied
-          Assert (this->n_nonzero_components(system_index) ==
-                  base_fe.n_nonzero_components(base_index),
-                  ExcInternalError());
-          for (unsigned int s=0; s<this->n_nonzero_components(system_index); ++s)
-            {
-              if (base_flags & update_values)
-                for (unsigned int q=0; q<n_q_points; ++q)
-                  output_data.shape_values[out_index+s][q] =
-                    base_data.shape_values(in_index+s,q);
-
-              if (base_flags & update_gradients)
-                for (unsigned int q=0; q<n_q_points; ++q)
-                  output_data.shape_gradients[out_index+s][q] =
-                    base_data.shape_gradients[in_index+s][q];
-
-              if (base_flags & update_hessians)
-                for (unsigned int q=0; q<n_q_points; ++q)
-                  output_data.shape_hessians[out_index+s][q] =
-                    base_data.shape_hessians[in_index+s][q];
-
-              if (base_flags & update_3rd_derivatives)
-                for (unsigned int q=0; q<n_q_points; ++q)
-                  output_data.shape_3rd_derivatives[out_index+s][q] =
-                    base_data.shape_3rd_derivatives[in_index+s][q];
-
-            }
-        }
-}
-
-
-
 template <int dim, int spacedim>
 template <int dim_1>
 void
@@ -1211,7 +1066,7 @@ compute_fill (const Mapping<dim,spacedim>                      &mapping,
               const unsigned int                                face_no,
               const unsigned int                                sub_no,
               const Quadrature<dim_1>                          &quadrature,
-              const CellSimilarity::Similarity                 ,
+              const CellSimilarity::Similarity                  cell_similarity,
               const typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
               const typename FiniteElement<dim,spacedim>::InternalDataBase &fedata,
               const internal::FEValues::MappingRelatedData<dim,spacedim> &mapping_data,
@@ -1231,18 +1086,141 @@ compute_fill (const Mapping<dim,spacedim>                      &mapping,
   const UpdateFlags flags = fe_data.update_each;
 
 
+  // loop over the base elements, let them compute what they need to compute,
+  // and then copy what is necessary.
+  //
+  // one may think that it would be a good idea to parallelize this over
+  // base elements, but it turns out to be not worthwhile: doing so lets
+  // multiple threads access data objects that were created by the current
+  // thread, leading to many NUMA memory access inefficiencies. we specifically
+  // want to avoid this if this class is called in a WorkStream context where
+  // we very carefully allocate objects only on the thread where they
+  // will actually be used; spawning new tasks here would be counterproductive
   if (flags & (update_values | update_gradients
                | update_hessians | update_3rd_derivatives ))
     for (unsigned int base_no=0; base_no<this->n_base_elements(); ++base_no)
       {
-        compute_fill_one_base<dim_1> (mapping,
-                                      std::make_pair(cell, CellSimilarity::none),
-                                      std::make_pair(face_no, sub_no),
-                                      quadrature,
-                                      std::make_pair(&mapping_internal, &fedata),
-                                      base_no,
-                                      mapping_data,
-                                      output_data);
+        const FiniteElement<dim,spacedim> &
+        base_fe      = base_element(base_no);
+        typename FiniteElement<dim,spacedim>::InternalDataBase &
+        base_fe_data = fe_data.get_fe_data(base_no);
+        internal::FEValues::FiniteElementRelatedData<dim,spacedim> &
+        base_data    = fe_data.get_fe_output_object(base_no);
+
+        // fill_fe_face_values needs argument Quadrature<dim-1> for both cases
+        // dim_1==dim-1 and dim_1=dim. Hence the following workaround
+        const Quadrature<dim>   *cell_quadrature = 0;
+        const Quadrature<dim-1> *face_quadrature = 0;
+        const unsigned int n_q_points = quadrature.size();
+
+        // static cast to the common base class of quadrature being either
+        // Quadrature<dim> or Quadrature<dim-1>:
+        const Subscriptor *quadrature_base_pointer = &quadrature;
+
+        if (face_no==invalid_face_number)
+          {
+            Assert(dim_1==dim, ExcDimensionMismatch(dim_1,dim));
+            Assert (dynamic_cast<const Quadrature<dim> *>(quadrature_base_pointer) != 0,
+                    ExcInternalError());
+
+            cell_quadrature
+              = static_cast<const Quadrature<dim> *>(quadrature_base_pointer);
+          }
+        else
+          {
+            Assert(dim_1==dim-1, ExcDimensionMismatch(dim_1,dim-1));
+            Assert (dynamic_cast<const Quadrature<dim-1> *>(quadrature_base_pointer) != 0,
+                    ExcInternalError());
+
+            face_quadrature
+              = static_cast<const Quadrature<dim-1> *>(quadrature_base_pointer);
+          }
+
+
+        // Make sure that in the case of fill_fe_values the data is only
+        // copied from base_data to data if base_data is changed. therefore
+        // use fe_fe_data.current_update_flags()
+        //
+        // for the case of fill_fe_(sub)face_values the data needs to be
+        // copied from base_data to data on each face, therefore use
+        // base_fe_data.update_flags.
+        if (face_no==invalid_face_number)
+          base_fe.fill_fe_values(mapping, cell, *cell_quadrature, mapping_internal,
+                                 base_fe_data, mapping_data, base_data, cell_similarity);
+        else if (sub_no==invalid_face_number)
+          base_fe.fill_fe_face_values(mapping, cell, face_no,
+                                      *face_quadrature, mapping_internal, base_fe_data, mapping_data, base_data);
+        else
+          base_fe.fill_fe_subface_values(mapping, cell, face_no, sub_no,
+                                         *face_quadrature, mapping_internal, base_fe_data, mapping_data, base_data);
+
+        // now data has been generated, so copy it. we used to work by
+        // looping over all base elements (i.e. this outer loop), then over
+        // multiplicity, then over the shape functions from that base
+        // element, but that requires that we can infer the global number of
+        // a shape function from its number in the base element. for that we
+        // had the component_to_system_table.
+        //
+        // however, this does of course no longer work since we have
+        // non-primitive elements. so we go the other way round: loop over
+        // 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 = 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
+        // enter this section
+        if (cell_similarity != CellSimilarity::translation)
+          for (unsigned int system_index=0; system_index<this->dofs_per_cell;
+               ++system_index)
+            if (this->system_to_base_table[system_index].first.first == base_no)
+              {
+                const unsigned int
+                base_index = this->system_to_base_table[system_index].second;
+                Assert (base_index<base_fe.dofs_per_cell, ExcInternalError());
+
+                // now copy. if the shape function is primitive, then there
+                // is only one value to be copied, but for non-primitive
+                // elements, there might be more values to be copied
+                //
+                // so, find out from which index to take this one value, and
+                // to which index to put
+                unsigned int out_index = 0;
+                for (unsigned int i=0; i<system_index; ++i)
+                  out_index += this->n_nonzero_components(i);
+                unsigned int in_index = 0;
+                for (unsigned int i=0; i<base_index; ++i)
+                  in_index += base_fe.n_nonzero_components(i);
+
+                // then loop over the number of components to be copied
+                Assert (this->n_nonzero_components(system_index) ==
+                        base_fe.n_nonzero_components(base_index),
+                        ExcInternalError());
+                for (unsigned int s=0; s<this->n_nonzero_components(system_index); ++s)
+                  {
+                    if (base_flags & update_values)
+                      for (unsigned int q=0; q<n_q_points; ++q)
+                        output_data.shape_values[out_index+s][q] =
+                          base_data.shape_values(in_index+s,q);
+
+                    if (base_flags & update_gradients)
+                      for (unsigned int q=0; q<n_q_points; ++q)
+                        output_data.shape_gradients[out_index+s][q] =
+                          base_data.shape_gradients[in_index+s][q];
+
+                    if (base_flags & update_hessians)
+                      for (unsigned int q=0; q<n_q_points; ++q)
+                        output_data.shape_hessians[out_index+s][q] =
+                          base_data.shape_hessians[in_index+s][q];
+
+                    if (base_flags & update_3rd_derivatives)
+                      for (unsigned int q=0; q<n_q_points; ++q)
+                        output_data.shape_3rd_derivatives[out_index+s][q] =
+                          base_data.shape_3rd_derivatives[in_index+s][q];
+
+                  }
+              }
       }
 }