From b2aa0a4827d8d8227334b35b0743914711ef95b9 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 30 Apr 2022 19:41:19 -0400 Subject: [PATCH] Implement proper indexing tables in FESystem::compute_fill(). This is about twice as fast as the old version. --- doc/news/changes/minor/20220430DavidWells | 4 + include/deal.II/fe/fe_system.h | 49 ++++ source/fe/fe_system.cc | 302 ++++++++++++++++------ 3 files changed, 276 insertions(+), 79 deletions(-) create mode 100644 doc/news/changes/minor/20220430DavidWells diff --git a/doc/news/changes/minor/20220430DavidWells b/doc/news/changes/minor/20220430DavidWells new file mode 100644 index 0000000000..9f5fa04ace --- /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. +
+(David Wells, 2022/04/30) diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h index 276daaf03a..520ddf6369 100644 --- 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::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> primitive_offset_tables; + + /** + * Lookup tables for indexing nonprimitive elements. + */ + std::vector> nonprimitive_offset_tables; + /** * Pointers to underlying finite element objects. * diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc index c105825f55..b54178bd9a 100644 --- 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 + Table<2, unsigned int> + setup_primitive_offset_table(const FESystem &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 + std::vector::BaseOffsets> + setup_nonprimitive_offset_table(const FESystem &fe, + const unsigned int base_no) + { + std::vector::BaseOffsets> table; + const FiniteElement &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 + void + copy_primitive_base_element_values( + const FESystem &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 + &base_data, + FEValuesImplementation::FiniteElementRelatedData + &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 + void + copy_nonprimitive_base_element_values( + const FESystem &fe, + const unsigned int base_no, + const unsigned int n_q_points, + const UpdateFlags base_flags, + const std::vector::BaseOffsets> &offsets, + const FEValuesImplementation::FiniteElementRelatedData + &base_data, + FEValuesImplementation::FiniteElementRelatedData + &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 FESystem::InternalData::InternalData( @@ -1320,86 +1499,33 @@ FESystem::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::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([&]() { -- 2.39.5