From 0668fd50b3c55b21510078a7eec494ac0c329845 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Thu, 19 May 2016 12:36:08 +0200 Subject: [PATCH] move build cell tables to FETools --- include/deal.II/fe/fe_system.h | 5 - include/deal.II/fe/fe_tools.h | 17 +++ source/fe/fe_system.cc | 198 ++------------------------------- source/fe/fe_tools.cc | 182 ++++++++++++++++++++++++++++++ source/fe/fe_tools.inst.in | 8 ++ 5 files changed, 219 insertions(+), 191 deletions(-) diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h index 8c936a7e12..40323b7f48 100644 --- a/include/deal.II/fe/fe_system.h +++ b/include/deal.II/fe/fe_system.h @@ -982,11 +982,6 @@ private: void initialize (const std::vector*> &fes, const std::vector &multiplicities); - /** - * Used by @p initialize. - */ - void build_cell_tables(); - /** * Used by @p initialize. */ diff --git a/include/deal.II/fe/fe_tools.h b/include/deal.II/fe/fe_tools.h index d86bb1a2bc..f3e022d0d7 100644 --- a/include/deal.II/fe/fe_tools.h +++ b/include/deal.II/fe/fe_tools.h @@ -932,6 +932,23 @@ namespace FETools const FiniteElement *fe5=NULL, const unsigned int N5=0); + /** + * For a given (composite) @p finite_element build @p system_to_component_table, + * @p system_to_base_table and @p component_to_base_table. + * + * If @p do_tensor_product is true, the underlying finite element + * is assumed to be build using the tensor product rule. That is, the number of + * composite components is the sum of components in each finite element times + * multiplicity. + */ + template + void + build_cell_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &system_to_base_table, + std::vector< std::pair< unsigned int, unsigned int > > &system_to_component_table, + std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &component_to_base_table, + const FiniteElement &finite_element, + const bool do_tensor_product = true); + /** * Parse the name of a finite element and generate a finite element object * accordingly. The parser ignores space characters between words (things diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc index a3af7aa6b7..f7b6c34efb 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -1265,191 +1265,6 @@ compute_fill (const Mapping &mapping, } - -template -void -FESystem::build_cell_tables() -{ - // If the system is not primitive, these have not been initialized by - // FiniteElement - this->system_to_component_table.resize(this->dofs_per_cell); - this->face_system_to_component_table.resize(this->dofs_per_face); - - unsigned int total_index = 0; - - for (unsigned int base=0; base < this->n_base_elements(); ++base) - for (unsigned int m = 0; m < this->element_multiplicity(base); ++m) - { - for (unsigned int k=0; kcomponent_to_base_table[total_index++] - = std::make_pair(std::make_pair(base,k), m); - } - Assert (total_index == this->component_to_base_table.size(), - ExcInternalError()); - - // Initialize index tables. Multi-component base elements have to be - // thought of. For non-primitive shape functions, have a special invalid - // index. - const std::pair - non_primitive_index (numbers::invalid_unsigned_int, - numbers::invalid_unsigned_int); - - // First enumerate vertex indices, where we first enumerate all indices on - // the first vertex in the order of the base elements, then of the second - // vertex, etc - total_index = 0; - for (unsigned int vertex_number=0; - vertex_number::vertices_per_cell; - ++vertex_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; basen_base_elements(); ++base) - for (unsigned int m=0; melement_multiplicity(base); - ++m, comp_start+=base_element(base).n_components()) - for (unsigned int local_index = 0; - local_index < base_element(base).dofs_per_vertex; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (base_element(base).dofs_per_vertex*vertex_number + - local_index); - - this->system_to_base_table[total_index] - = std::make_pair (std::make_pair(base, m), index_in_base); - - if (base_element(base).is_primitive(index_in_base)) - { - const unsigned int comp_in_base - = base_element(base).system_to_component_index(index_in_base).first; - const unsigned int comp - = comp_start + comp_in_base; - const unsigned int index_in_comp - = base_element(base).system_to_component_index(index_in_base).second; - this->system_to_component_table[total_index] - = std::make_pair (comp, index_in_comp); - } - else - this->system_to_component_table[total_index] = non_primitive_index; - } - } - - // 2. Lines - if (GeometryInfo::lines_per_cell > 0) - for (unsigned int line_number= 0; - line_number != GeometryInfo::lines_per_cell; - ++line_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; basen_base_elements(); ++base) - for (unsigned int m=0; melement_multiplicity(base); - ++m, comp_start+=base_element(base).n_components()) - for (unsigned int local_index = 0; - local_index < base_element(base).dofs_per_line; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (base_element(base).dofs_per_line*line_number + - local_index + - base_element(base).first_line_index); - - this->system_to_base_table[total_index] - = std::make_pair (std::make_pair(base,m), index_in_base); - - if (base_element(base).is_primitive(index_in_base)) - { - const unsigned int comp_in_base - = base_element(base).system_to_component_index(index_in_base).first; - const unsigned int comp - = comp_start + comp_in_base; - const unsigned int index_in_comp - = base_element(base).system_to_component_index(index_in_base).second; - this->system_to_component_table[total_index] - = std::make_pair (comp, index_in_comp); - } - else - this->system_to_component_table[total_index] = non_primitive_index; - } - } - - // 3. Quads - if (GeometryInfo::quads_per_cell > 0) - for (unsigned int quad_number= 0; - quad_number != GeometryInfo::quads_per_cell; - ++quad_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; basen_base_elements(); ++base) - for (unsigned int m=0; melement_multiplicity(base); - ++m, comp_start += base_element(base).n_components()) - for (unsigned int local_index = 0; - local_index < base_element(base).dofs_per_quad; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (base_element(base).dofs_per_quad*quad_number + - local_index + - base_element(base).first_quad_index); - - this->system_to_base_table[total_index] - = std::make_pair (std::make_pair(base,m), index_in_base); - - if (base_element(base).is_primitive(index_in_base)) - { - const unsigned int comp_in_base - = base_element(base).system_to_component_index(index_in_base).first; - const unsigned int comp - = comp_start + comp_in_base; - const unsigned int index_in_comp - = base_element(base).system_to_component_index(index_in_base).second; - this->system_to_component_table[total_index] - = std::make_pair (comp, index_in_comp); - } - else - this->system_to_component_table[total_index] = non_primitive_index; - } - } - - // 4. Hexes - if (GeometryInfo::hexes_per_cell > 0) - for (unsigned int hex_number= 0; - hex_number != GeometryInfo::hexes_per_cell; - ++hex_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; basen_base_elements(); ++base) - for (unsigned int m=0; melement_multiplicity(base); - ++m, comp_start+=base_element(base).n_components()) - for (unsigned int local_index = 0; - local_index < base_element(base).dofs_per_hex; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (base_element(base).dofs_per_hex*hex_number + - local_index + - base_element(base).first_hex_index); - - this->system_to_base_table[total_index] - = std::make_pair (std::make_pair(base,m), index_in_base); - - if (base_element(base).is_primitive(index_in_base)) - { - const unsigned int comp_in_base - = base_element(base).system_to_component_index(index_in_base).first; - const unsigned int comp - = comp_start + comp_in_base; - const unsigned int index_in_comp - = base_element(base).system_to_component_index(index_in_base).second; - this->system_to_component_table[total_index] - = std::make_pair (comp, index_in_comp); - } - else - this->system_to_component_table[total_index] = non_primitive_index; - } - } -} - - - template void FESystem::build_face_tables() @@ -1831,7 +1646,18 @@ void FESystem::initialize (const std::vector0, ExcInternalError()); - build_cell_tables(); + { + // If the system is not primitive, these have not been initialized by + // FiniteElement + this->system_to_component_table.resize(this->dofs_per_cell); + this->face_system_to_component_table.resize(this->dofs_per_face); + + FETools::build_cell_tables(this->system_to_base_table, + this->system_to_component_table, + this->component_to_base_table, + *this); + + } build_face_tables(); // restriction and prolongation matrices are build on demand diff --git a/source/fe/fe_tools.cc b/source/fe/fe_tools.cc index 06edebbab9..c2ebb2a222 100644 --- a/source/fe/fe_tools.cc +++ b/source/fe/fe_tools.cc @@ -558,6 +558,188 @@ namespace FETools return compute_nonzero_components (fe_list, multiplicities); } + template + void + build_cell_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &system_to_base_table, + std::vector< std::pair< unsigned int, unsigned int > > &system_to_component_table, + std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &component_to_base_table, + const FiniteElement &fe, + const bool do_tensor_product) + { + unsigned int total_index = 0; + + for (unsigned int base=0; base < fe.n_base_elements(); ++base) + for (unsigned int m = 0; m < fe.element_multiplicity(base); ++m) + { + for (unsigned int k=0; k + non_primitive_index (numbers::invalid_unsigned_int, + numbers::invalid_unsigned_int); + + // First enumerate vertex indices, where we first enumerate all indices on + // the first vertex in the order of the base elements, then of the second + // vertex, etc + total_index = 0; + for (unsigned int vertex_number=0; + vertex_number::vertices_per_cell; + ++vertex_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; base::lines_per_cell > 0) + for (unsigned int line_number= 0; + line_number != GeometryInfo::lines_per_cell; + ++line_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; base::quads_per_cell > 0) + for (unsigned int quad_number= 0; + quad_number != GeometryInfo::quads_per_cell; + ++quad_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; base::hexes_per_cell > 0) + for (unsigned int hex_number= 0; + hex_number != GeometryInfo::hexes_per_cell; + ++hex_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; base FiniteElement * diff --git a/source/fe/fe_tools.inst.in b/source/fe/fe_tools.inst.in index 90535f80ef..6c7a389d51 100644 --- a/source/fe/fe_tools.inst.in +++ b/source/fe/fe_tools.inst.in @@ -74,6 +74,14 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS const unsigned int N4, const FiniteElement *fe5, const unsigned int N5); + + template + void + build_cell_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &system_to_base_table, + std::vector< std::pair< unsigned int, unsigned int > > &system_to_component_table, + std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &component_to_base_table, + const FiniteElement &fe, + const bool do_tensor_product); template void compute_block_renumbering ( -- 2.39.5