Implement proper indexing tables in FESystem::compute_fill().

This is about twice as fast as the old version.

diff --git a/doc/news/changes/minor/20220430DavidWells b/doc/news/changes/minor/20220430DavidWells
new file mode 100644 (file)
index 0000000..9f5fa04
--- /dev/null
+++ b/doc/news/changes/minor/20220430DavidWells
@@ -0,0 +1,4 @@
+Changed: Improved: The internal indexing in FESystem::fill_fe_values() has been
+improved and is now about twice as fast.
+<br>
+(David Wells, 2022/04/30)

diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h
index 276daaf03a6d42ea78076c45669f43b401876915..520ddf6369b7f22acd13a3825ec7cdb2fee63fd0 100644 (file)
--- a/include/deal.II/fe/fe_system.h
+++ b/include/deal.II/fe/fe_system.h
@@ -1049,6 +1049,39 @@ public:
   virtual std::size_t
   memory_consumption() const override;
 
+  /**
+   * Computing values on a cell requires copying values between the
+   * FEValuesImplementation::FiniteElementRelatedData objects owned by the
+   * present object and its subelements. The indexing for this is fairly
+   * complicated since each shape function in the nonprimitive base element may
+   * have a different number of nonzero components (e.g., a Taylor-Hood pair
+   * base element): we also need offsets into both the base element and present
+   * element's arrays. To this end, we will create a table with all the
+   * necessary information for each shape function in the nonprimitive base
+   * element.
+   */
+  struct BaseOffsets
+  {
+    /**
+     * Number of nonzero components for the current base element shape function.
+     */
+    unsigned int n_nonzero_components;
+
+    /**
+     * Index into the base element's
+     * FEValuesImplementation::FiniteElementRelatedData arrays for the current
+     * base element shape function.
+     */
+    unsigned int in_index;
+
+    /**
+     * Index into the present element's
+     * FEValuesImplementation::FiniteElementRelatedData arrays for the current
+     * base element shape function.
+     */
+    unsigned int out_index;
+  };
+
 protected:
   virtual std::unique_ptr<
     typename FiniteElement<dim, spacedim>::InternalDataBase>
@@ -1162,6 +1195,22 @@ private:
    */
   static const unsigned int invalid_face_number = numbers::invalid_unsigned_int;
 
+  /**
+   * Lookup tables for indexing primitive elements.
+   *
+   * These tables store, for each primitive base element, the array indices into
+   * which we will write values. One can think of these indices as a combination
+   * of the system index and the number of nonzero components per shape
+   * function. This data is needed to copy values from base element value
+   * reinitialization to the present element's value reinitialization.
+   */
+  std::vector<Table<2, unsigned int>> primitive_offset_tables;
+
+  /**
+   * Lookup tables for indexing nonprimitive elements.
+   */
+  std::vector<std::vector<BaseOffsets>> nonprimitive_offset_tables;
+
   /**
    * Pointers to underlying finite element objects.
    *

diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc
index c105825f55e3a1c990a1df1fd44d045408bcbccc..b54178bd9abf5dd1e05ed46e27ebf7168eff4f94 100644 (file)
--- a/source/fe/fe_system.cc
+++ b/source/fe/fe_system.cc
@@ -42,8 +42,187 @@ namespace
     });
   }
 } // namespace
-/* ----------------------- FESystem::InternalData ------------------- */
 
+namespace internal
+{
+  /**
+   * Setup a table of offsets for a primitive FE. Unlike the nonprimitive
+   * case, here the number of nonzero components per shape function is always
+   * 1 and the number of components in the FE is always the multiplicity.
+   */
+  template <int dim, int spacedim = dim>
+  Table<2, unsigned int>
+  setup_primitive_offset_table(const FESystem<dim, spacedim> &fe,
+                               const unsigned int             base_no)
+  {
+    Assert(fe.base_element(base_no).is_primitive(), ExcInternalError());
+    Table<2, unsigned int> table(fe.element_multiplicity(base_no),
+                                 fe.base_element(base_no).n_dofs_per_cell());
+    // 0 is a bad default value since it is a valid index
+    table.fill(numbers::invalid_unsigned_int);
+
+    unsigned int out_index = 0;
+    for (unsigned int system_index = 0; system_index < fe.n_dofs_per_cell();
+         ++system_index)
+      {
+        if (fe.system_to_base_index(system_index).first.first == base_no)
+          {
+            Assert(fe.n_nonzero_components(system_index) == 1,
+                   ExcInternalError());
+            const unsigned int base_component =
+              fe.system_to_base_index(system_index).first.second;
+            const unsigned int base_index =
+              fe.system_to_base_index(system_index).second;
+            Assert(base_index < fe.base_element(base_no).n_dofs_per_cell(),
+                   ExcInternalError());
+
+            table[base_component][base_index] = out_index;
+          }
+        out_index += fe.n_nonzero_components(system_index);
+      }
+
+    return table;
+  }
+
+  /**
+   * Setup a table of offsets for a nonprimitive FE.
+   */
+  template <int dim, int spacedim = dim>
+  std::vector<typename FESystem<dim, spacedim>::BaseOffsets>
+  setup_nonprimitive_offset_table(const FESystem<dim, spacedim> &fe,
+                                  const unsigned int             base_no)
+  {
+    std::vector<typename FESystem<dim, spacedim>::BaseOffsets> table;
+    const FiniteElement<dim, spacedim> &base_fe = fe.base_element(base_no);
+
+    unsigned int out_index = 0;
+    for (unsigned int system_index = 0; system_index < fe.n_dofs_per_cell();
+         ++system_index)
+      {
+        if (fe.system_to_base_index(system_index).first.first == base_no)
+          {
+            const unsigned int base_index =
+              fe.system_to_base_index(system_index).second;
+            Assert(base_index < base_fe.n_dofs_per_cell(), ExcInternalError());
+            table.emplace_back();
+
+            table.back().n_nonzero_components =
+              fe.n_nonzero_components(system_index);
+            unsigned int in_index = 0;
+            for (unsigned int i = 0; i < base_index; ++i)
+              in_index += base_fe.n_nonzero_components(i);
+
+            table.back().in_index  = in_index;
+            table.back().out_index = out_index;
+          }
+        out_index += fe.n_nonzero_components(system_index);
+      }
+
+    Assert(table.size() ==
+             base_fe.n_dofs_per_cell() * fe.element_multiplicity(base_no),
+           ExcInternalError());
+    return table;
+  }
+
+  /**
+   * Copy data between internal FEValues objects from a primitive FE to the
+   * current FE.
+   */
+  template <int dim, int spacedim = dim>
+  void
+  copy_primitive_base_element_values(
+    const FESystem<dim, spacedim> &fe,
+    const unsigned int             base_no,
+    const unsigned int             n_q_points,
+    const UpdateFlags              base_flags,
+    const Table<2, unsigned int> & base_to_system_table,
+    const FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>
+      &base_data,
+    FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>
+      &output_data)
+  {
+    Assert(fe.base_element(base_no).is_primitive(), ExcInternalError());
+    const unsigned int n_components = fe.element_multiplicity(base_no);
+    const unsigned int n_dofs_per_cell =
+      fe.base_element(base_no).n_dofs_per_cell();
+    for (unsigned int component = 0; component < n_components; ++component)
+      for (unsigned int b = 0; b < n_dofs_per_cell; ++b)
+        {
+          const unsigned int out_index = base_to_system_table[component][b];
+
+          if (base_flags & update_values)
+            for (unsigned int q = 0; q < n_q_points; ++q)
+              output_data.shape_values[out_index][q] =
+                base_data.shape_values[b][q];
+
+          if (base_flags & update_gradients)
+            for (unsigned int q = 0; q < n_q_points; ++q)
+              output_data.shape_gradients[out_index][q] =
+                base_data.shape_gradients[b][q];
+
+          if (base_flags & update_hessians)
+            for (unsigned int q = 0; q < n_q_points; ++q)
+              output_data.shape_hessians[out_index][q] =
+                base_data.shape_hessians[b][q];
+
+          if (base_flags & update_3rd_derivatives)
+            for (unsigned int q = 0; q < n_q_points; ++q)
+              output_data.shape_3rd_derivatives[out_index][q] =
+                base_data.shape_3rd_derivatives[b][q];
+        }
+  }
+
+  /**
+   * Copy data between internal FEValues objects from a nonprimitive FE to the
+   * current FE.
+   */
+  template <int dim, int spacedim = dim>
+  void
+  copy_nonprimitive_base_element_values(
+    const FESystem<dim, spacedim> &fe,
+    const unsigned int             base_no,
+    const unsigned int             n_q_points,
+    const UpdateFlags              base_flags,
+    const std::vector<typename FESystem<dim, spacedim>::BaseOffsets> &offsets,
+    const FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>
+      &base_data,
+    FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>
+      &output_data)
+  {
+    (void)fe;
+    (void)base_no;
+    Assert(!fe.base_element(base_no).is_primitive(), ExcInternalError());
+
+    for (const auto &offset : offsets)
+      {
+        if (base_flags & update_values)
+          for (unsigned int s = 0; s < offset.n_nonzero_components; ++s)
+            for (unsigned int q = 0; q < n_q_points; ++q)
+              output_data.shape_values[offset.out_index + s][q] =
+                base_data.shape_values[offset.in_index + s][q];
+
+        if (base_flags & update_gradients)
+          for (unsigned int s = 0; s < offset.n_nonzero_components; ++s)
+            for (unsigned int q = 0; q < n_q_points; ++q)
+              output_data.shape_gradients[offset.out_index + s][q] =
+                base_data.shape_gradients[offset.in_index + s][q];
+
+        if (base_flags & update_hessians)
+          for (unsigned int s = 0; s < offset.n_nonzero_components; ++s)
+            for (unsigned int q = 0; q < n_q_points; ++q)
+              output_data.shape_hessians[offset.out_index + s][q] =
+                base_data.shape_hessians[offset.in_index + s][q];
+
+        if (base_flags & update_3rd_derivatives)
+          for (unsigned int s = 0; s < offset.n_nonzero_components; ++s)
+            for (unsigned int q = 0; q < n_q_points; ++q)
+              output_data.shape_3rd_derivatives[offset.out_index + s][q] =
+                base_data.shape_3rd_derivatives[offset.in_index + s][q];
+      }
+  }
+} // namespace internal
+
+/* ----------------------- FESystem::InternalData ------------------- */
 
 template <int dim, int spacedim>
 FESystem<dim, spacedim>::InternalData::InternalData(
@@ -1320,86 +1499,33 @@ FESystem<dim, spacedim>::compute_fill(
                                          base_fe_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
+        // now data has been generated, so copy it. This procedure is different
+        // for primitive and non-primitive base elements, so at this point we
+        // dispatch to helper functions.
         const UpdateFlags base_flags = base_fe_data.update_each;
 
-        // some base element might involve values that depend on the shape
-        // of the geometry, so we always need to copy the shape values around
-        // also in case we detected a cell similarity (but no heavy work will
-        // be done inside the individual elements in case we have a
-        // translation and simple elements).
-        for (unsigned int system_index = 0;
-             system_index < this->n_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.n_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());
-
-              if (base_flags & update_values)
-                for (unsigned int s = 0;
-                     s < this->n_nonzero_components(system_index);
-                     ++s)
-                  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 s = 0;
-                     s < this->n_nonzero_components(system_index);
-                     ++s)
-                  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 s = 0;
-                     s < this->n_nonzero_components(system_index);
-                     ++s)
-                  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 s = 0;
-                     s < this->n_nonzero_components(system_index);
-                     ++s)
-                  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];
-            }
+        if (base_fe.is_primitive())
+          {
+            internal::copy_primitive_base_element_values(
+              *this,
+              base_no,
+              n_q_points,
+              base_flags,
+              primitive_offset_tables[base_no],
+              base_data,
+              output_data);
+          }
+        else
+          {
+            internal::copy_nonprimitive_base_element_values(
+              *this,
+              base_no,
+              n_q_points,
+              base_flags,
+              nonprimitive_offset_tables[base_no],
+              base_data,
+              output_data);
+          }
       }
 }
 
@@ -1732,6 +1858,24 @@ FESystem<dim, spacedim>::initialize(
       }
   });
 
+  init_tasks += Threads::new_task([&]() {
+    primitive_offset_tables.resize(this->n_base_elements());
+
+    for (unsigned int base_no = 0; base_no < this->n_base_elements(); ++base_no)
+      if (base_element(base_no).is_primitive())
+        primitive_offset_tables[base_no] =
+          internal::setup_primitive_offset_table(*this, base_no);
+  });
+
+  init_tasks += Threads::new_task([&]() {
+    nonprimitive_offset_tables.resize(this->n_base_elements());
+
+    for (unsigned int base_no = 0; base_no < this->n_base_elements(); ++base_no)
+      if (!base_element(base_no).is_primitive())
+        nonprimitive_offset_tables[base_no] =
+          internal::setup_nonprimitive_offset_table(*this, base_no);
+  });
+
   // initialize face support points (for dim==2,3). same procedure as above
   if (dim > 1)
     init_tasks += Threads::new_task([&]() {