From 9b4a09104caea3841abe91392e1d6a0adfd444a6 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Thu, 19 May 2016 13:39:08 +0200 Subject: [PATCH] same for fe_face_tables --- include/deal.II/fe/fe_system.h | 5 - include/deal.II/fe/fe_tools.h | 14 +++ source/fe/fe_system.cc | 162 +-------------------------------- source/fe/fe_tools.cc | 159 ++++++++++++++++++++++++++++++++ source/fe/fe_tools.inst.in | 7 ++ 5 files changed, 184 insertions(+), 163 deletions(-) diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h index 40323b7f48..98064474d9 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_face_tables(); - /** * Used by @p initialize. */ diff --git a/include/deal.II/fe/fe_tools.h b/include/deal.II/fe/fe_tools.h index f3e022d0d7..14fb9995e7 100644 --- a/include/deal.II/fe/fe_tools.h +++ b/include/deal.II/fe/fe_tools.h @@ -949,6 +949,20 @@ namespace FETools const FiniteElement &finite_element, const bool do_tensor_product = true); + /** + * build... + * @param face_system_to_base_table + * @param face_system_to_component_table + * @param finite_element + * @param do_tensor_product + */ + template + void + build_face_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &face_system_to_base_table, + std::vector< std::pair< unsigned int, unsigned int > > &face_system_to_component_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 f7b6c34efb..9f1ff5a4dc 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -1265,163 +1265,6 @@ compute_fill (const Mapping &mapping, } -template -void -FESystem::build_face_tables() -{ - // Initialize index tables. do this in the same way as done for the cell - // tables, except that we now loop over the objects of faces - - // For non-primitive shape functions, have a special invalid index - const std::pair - non_primitive_index (numbers::invalid_unsigned_int, - numbers::invalid_unsigned_int); - - // 1. Vertices - unsigned int total_index = 0; - for (unsigned int vertex_number=0; - vertex_number::vertices_per_face; - ++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) - { - // get (cell) index of this shape function inside the base - // element to see whether the shape function is primitive - // (assume that all shape functions on vertices share the same - // primitivity property; assume likewise for all shape functions - // located on lines, quads, etc. this way, we can ask for - // primitivity of only _one_ shape function, which is taken as - // representative for all others located on the same type of - // object): - const unsigned int index_in_base - = (base_element(base).dofs_per_vertex*vertex_number + - local_index); - - const unsigned int face_index_in_base - = (base_element(base).dofs_per_vertex*vertex_number + - local_index); - - this->face_system_to_base_table[total_index] - = std::make_pair (std::make_pair(base,m), face_index_in_base); - - if (base_element(base).is_primitive(index_in_base)) - { - const unsigned int comp_in_base - = base_element(base).face_system_to_component_index(face_index_in_base).first; - const unsigned int comp - = comp_start + comp_in_base; - const unsigned int face_index_in_comp - = base_element(base).face_system_to_component_index(face_index_in_base).second; - this->face_system_to_component_table[total_index] - = std::make_pair (comp, face_index_in_comp); - } - else - this->face_system_to_component_table[total_index] = non_primitive_index; - } - } - - // 2. Lines - if (GeometryInfo::lines_per_face > 0) - for (unsigned int line_number= 0; - line_number != GeometryInfo::lines_per_face; - ++line_number) - { - unsigned int comp_start = 0; - for (unsigned int base = 0; base < this->n_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) - { - // do everything alike for this type of object - const unsigned int index_in_base - = (base_element(base).dofs_per_line*line_number + - local_index + - base_element(base).first_line_index); - - const unsigned int face_index_in_base - = (base_element(base).first_face_line_index + - base_element(base).dofs_per_line * line_number + - local_index); - - this->face_system_to_base_table[total_index] - = std::make_pair (std::make_pair(base,m), face_index_in_base); - - if (base_element(base).is_primitive(index_in_base)) - { - const unsigned int comp_in_base - = base_element(base).face_system_to_component_index(face_index_in_base).first; - const unsigned int comp - = comp_start + comp_in_base; - const unsigned int face_index_in_comp - = base_element(base).face_system_to_component_index(face_index_in_base).second; - this->face_system_to_component_table[total_index] - = std::make_pair (comp, face_index_in_comp); - } - else - this->face_system_to_component_table[total_index] = non_primitive_index; - } - } - - // 3. Quads - if (GeometryInfo::quads_per_face > 0) - for (unsigned int quad_number= 0; - quad_number != GeometryInfo::quads_per_face; - ++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) - { - // do everything alike for this type of object - const unsigned int index_in_base - = (base_element(base).dofs_per_quad*quad_number + - local_index + - base_element(base).first_quad_index); - - const unsigned int face_index_in_base - = (base_element(base).first_face_quad_index + - base_element(base).dofs_per_quad * quad_number + - local_index); - - this->face_system_to_base_table[total_index] - = std::make_pair (std::make_pair(base,m), face_index_in_base); - - if (base_element(base).is_primitive(index_in_base)) - { - const unsigned int comp_in_base - = base_element(base).face_system_to_component_index(face_index_in_base).first; - const unsigned int comp - = comp_start + comp_in_base; - const unsigned int face_index_in_comp - = base_element(base).face_system_to_component_index(face_index_in_base).second; - this->face_system_to_component_table[total_index] - = std::make_pair (comp, face_index_in_comp); - } - else - this->face_system_to_component_table[total_index] = non_primitive_index; - } - } - Assert (total_index == this->dofs_per_face, ExcInternalError()); - Assert (total_index == this->face_system_to_component_table.size(), - ExcInternalError()); - Assert (total_index == this->face_system_to_base_table.size(), - ExcInternalError()); -} - - - template void FESystem::build_interface_constraints () { @@ -1657,8 +1500,11 @@ void FESystem::initialize (const std::vectorcomponent_to_base_table, *this); + FETools::build_face_tables(this->face_system_to_base_table, + this->face_system_to_component_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 c2ebb2a222..8f29669e06 100644 --- a/source/fe/fe_tools.cc +++ b/source/fe/fe_tools.cc @@ -740,6 +740,165 @@ namespace FETools } } + template + void + build_face_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &face_system_to_base_table, + std::vector< std::pair< unsigned int, unsigned int > > &face_system_to_component_table, + const FiniteElement &fe, + const bool do_tensor_product) + { + // Initialize index tables. do this in the same way as done for the cell + // tables, except that we now loop over the objects of faces + + // For non-primitive shape functions, have a special invalid index + const std::pair + non_primitive_index (numbers::invalid_unsigned_int, + numbers::invalid_unsigned_int); + + // 1. Vertices + unsigned int total_index = 0; + for (unsigned int vertex_number=0; + vertex_number::vertices_per_face; + ++vertex_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; base::lines_per_face > 0) + for (unsigned int line_number= 0; + line_number != GeometryInfo::lines_per_face; + ++line_number) + { + unsigned int comp_start = 0; + for (unsigned int base = 0; base < fe.n_base_elements(); ++base) + for (unsigned int m=0; m::quads_per_face > 0) + for (unsigned int quad_number= 0; + quad_number != GeometryInfo::quads_per_face; + ++quad_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 6c7a389d51..7a8aa495fe 100644 --- a/source/fe/fe_tools.inst.in +++ b/source/fe/fe_tools.inst.in @@ -82,6 +82,13 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS 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 + build_face_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &face_system_to_base_table, + std::vector< std::pair< unsigned int, unsigned int > > &face_system_to_component_table, + const FiniteElement &fe, + const bool do_tensor_product); template void compute_block_renumbering ( -- 2.39.5