From 9b4a09104caea3841abe91392e1d6a0adfd444a6 Mon Sep 17 00:00:00 2001 From: Denis Davydov <davydden@gmail.com> 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<const FiniteElement<dim,spacedim>*> &fes, const std::vector<unsigned int> &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<dim,spacedim> &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 <int dim, int spacedim> + 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<dim,spacedim> &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<dim,spacedim> &mapping, } -template <int dim, int spacedim> -void -FESystem<dim,spacedim>::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<unsigned int, unsigned int> - 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<GeometryInfo<dim>::vertices_per_face; - ++vertex_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; base<this->n_base_elements(); ++base) - for (unsigned int m=0; m<this->element_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<dim>::lines_per_face > 0) - for (unsigned int line_number= 0; - line_number != GeometryInfo<dim>::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; m<this->element_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<dim>::quads_per_face > 0) - for (unsigned int quad_number= 0; - quad_number != GeometryInfo<dim>::quads_per_face; - ++quad_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; base<this->n_base_elements(); ++base) - for (unsigned int m=0; m<this->element_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 <int dim, int spacedim> void FESystem<dim,spacedim>::build_interface_constraints () { @@ -1657,8 +1500,11 @@ void FESystem<dim,spacedim>::initialize (const std::vector<const FiniteElement<d this->component_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 <int dim, int spacedim> + 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<dim,spacedim> &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<unsigned int, unsigned int> + 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<GeometryInfo<dim>::vertices_per_face; + ++vertex_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; base<fe.n_base_elements(); ++base) + for (unsigned int m=0; m<fe.element_multiplicity(base); + ++m, comp_start += fe.base_element(base).n_components()) + for (unsigned int local_index = 0; + local_index < fe.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 + = (fe.base_element(base).dofs_per_vertex*vertex_number + + local_index); + + const unsigned int face_index_in_base + = (fe.base_element(base).dofs_per_vertex*vertex_number + + local_index); + + face_system_to_base_table[total_index] + = std::make_pair (std::make_pair(base,m), face_index_in_base); + + if (fe.base_element(base).is_primitive(index_in_base)) + { + const unsigned int comp_in_base + = fe.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 + = fe.base_element(base).face_system_to_component_index(face_index_in_base).second; + face_system_to_component_table[total_index] + = std::make_pair (comp, face_index_in_comp); + } + else + face_system_to_component_table[total_index] = non_primitive_index; + } + } + + // 2. Lines + if (GeometryInfo<dim>::lines_per_face > 0) + for (unsigned int line_number= 0; + line_number != GeometryInfo<dim>::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<fe.element_multiplicity(base); + ++m, comp_start += fe.base_element(base).n_components()) + for (unsigned int local_index = 0; + local_index < fe.base_element(base).dofs_per_line; + ++local_index, ++total_index) + { + // do everything alike for this type of object + const unsigned int index_in_base + = (fe.base_element(base).dofs_per_line*line_number + + local_index + + fe.base_element(base).first_line_index); + + const unsigned int face_index_in_base + = (fe.base_element(base).first_face_line_index + + fe.base_element(base).dofs_per_line * line_number + + local_index); + + face_system_to_base_table[total_index] + = std::make_pair (std::make_pair(base,m), face_index_in_base); + + if (fe.base_element(base).is_primitive(index_in_base)) + { + const unsigned int comp_in_base + = fe.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 + = fe.base_element(base).face_system_to_component_index(face_index_in_base).second; + face_system_to_component_table[total_index] + = std::make_pair (comp, face_index_in_comp); + } + else + face_system_to_component_table[total_index] = non_primitive_index; + } + } + + // 3. Quads + if (GeometryInfo<dim>::quads_per_face > 0) + for (unsigned int quad_number= 0; + quad_number != GeometryInfo<dim>::quads_per_face; + ++quad_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; base<fe.n_base_elements(); ++base) + for (unsigned int m=0; m<fe.element_multiplicity(base); + ++m, comp_start += fe.base_element(base).n_components()) + for (unsigned int local_index = 0; + local_index < fe.base_element(base).dofs_per_quad; + ++local_index, ++total_index) + { + // do everything alike for this type of object + const unsigned int index_in_base + = (fe.base_element(base).dofs_per_quad*quad_number + + local_index + + fe.base_element(base).first_quad_index); + + const unsigned int face_index_in_base + = (fe.base_element(base).first_face_quad_index + + fe.base_element(base).dofs_per_quad * quad_number + + local_index); + + face_system_to_base_table[total_index] + = std::make_pair (std::make_pair(base,m), face_index_in_base); + + if (fe.base_element(base).is_primitive(index_in_base)) + { + const unsigned int comp_in_base + = fe.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 + = fe.base_element(base).face_system_to_component_index(face_index_in_base).second; + face_system_to_component_table[total_index] + = std::make_pair (comp, face_index_in_comp); + } + else + face_system_to_component_table[total_index] = non_primitive_index; + } + } + Assert (total_index == fe.dofs_per_face, ExcInternalError()); + Assert (total_index == face_system_to_component_table.size(), + ExcInternalError()); + Assert (total_index == face_system_to_base_table.size(), + ExcInternalError()); + } + + // Not implemented in the general case. template <class FE> FiniteElement<FE::dimension, FE::space_dimension> * 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<deal_II_dimension,deal_II_space_dimension> &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<deal_II_dimension,deal_II_space_dimension> &fe, + const bool do_tensor_product); template void compute_block_renumbering ( -- 2.39.5