From db5c8069debb8195839c4a108ac6b6f553688ede Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sat, 23 Jan 2016 11:43:40 -0600 Subject: [PATCH] Move some internal functions. FESystem has a number of functions that are called in the member initializer list of the constructors. These are all static functions, of course, and don't access any member variables. This patch moves them ouf of the class altogether. This has the advantage that it reduces the size of the .h file, and that it reduces the number of functions exported by the shared libraries. Rather, these functions may as well live inside an anonymous namespace in the .cc file. There are no functional changes, nor anything that would be visible to the user. --- include/deal.II/fe/fe_system.h | 82 --- source/fe/fe_system.cc | 955 +++++++++++++++++---------------- 2 files changed, 485 insertions(+), 552 deletions(-) diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h index bcb53a1e10..cb00921b29 100644 --- a/include/deal.II/fe/fe_system.h +++ b/include/deal.II/fe/fe_system.h @@ -983,88 +983,6 @@ private: * of the FiniteElement class. Called from the constructor. */ void initialize_quad_dof_index_permutation (); - - /** - * Helper function used in the constructor: take a @p FiniteElementData - * object and return an object of the same type with the number of degrees - * of freedom per vertex, line, etc. multiplied by @p n. Don't touch the - * number of functions for the transformation from unit to real cell. - */ - static FiniteElementData - multiply_dof_numbers (const FiniteElement *fe1, - const unsigned int N1, - const FiniteElement *fe2=NULL, - const unsigned int N2=0, - const FiniteElement *fe3=NULL, - const unsigned int N3=0, - const FiniteElement *fe4=NULL, - const unsigned int N4=0, - const FiniteElement *fe5=NULL, - const unsigned int N5=0); - - /** - * Same as above but for any number of sub-elements. - */ - static FiniteElementData - multiply_dof_numbers (const std::vector*> &fes, - const std::vector &multiplicities); - - - - /** - * Helper function used in the constructor: takes a @p FiniteElement object - * and returns an boolean vector including the @p - * restriction_is_additive_flags of the mixed element consisting of @p N - * elements of the sub-element @p fe. - */ - static std::vector - compute_restriction_is_additive_flags ( - const FiniteElement *fe1, - const unsigned int N1, - const FiniteElement *fe2=NULL, - const unsigned int N2=0, - const FiniteElement *fe3=NULL, - const unsigned int N3=0, - const FiniteElement *fe4=NULL, - const unsigned int N4=0, - const FiniteElement *fe5=NULL, - const unsigned int N5=0); - - /** - * Compute the named flags for a list of finite elements with multiplicities - * given in the second argument. This function is called from all the above - * functions. - */ - static std::vector - compute_restriction_is_additive_flags ( - const std::vector*> &fes, - const std::vector &multiplicities); - - - /** - * Compute the non-zero vector components of a composed finite element. - */ - static std::vector - compute_nonzero_components (const FiniteElement *fe1, - const unsigned int N1, - const FiniteElement *fe2=NULL, - const unsigned int N2=0, - const FiniteElement *fe3=NULL, - const unsigned int N3=0, - const FiniteElement *fe4=NULL, - const unsigned int N4=0, - const FiniteElement *fe5=NULL, - const unsigned int N5=0); - - /** - * Compute the nonzero components of a list of finite elements with - * multiplicities given in the second argument. This function is called from - * all the above functions. - */ - static std::vector - compute_nonzero_components (const std::vector*> &fes, - const std::vector &multiplicities); - /** * This function is simply singled out of the constructors since there are * several of them. It sets up the index table for the system as well as @p diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc index 2d1265f22f..aefbfa6849 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -110,6 +110,491 @@ InternalData::get_fe_output_object (const unsigned int base_no) const template const unsigned int FESystem::invalid_face_number; +namespace +{ + + /** + * Take vectors of finite elements and multiplicities and multiply out + * how many degrees of freedom the composed element has per vertex, + * line, etc. + */ + template + FiniteElementData + multiply_dof_numbers (const std::vector*> &fes, + const std::vector &multiplicities) + { + Assert (fes.size() == multiplicities.size(), ExcDimensionMismatch (fes.size(), multiplicities.size())); + + unsigned int multiplied_dofs_per_vertex = 0; + unsigned int multiplied_dofs_per_line = 0; + unsigned int multiplied_dofs_per_quad = 0; + unsigned int multiplied_dofs_per_hex = 0; + + unsigned int multiplied_n_components = 0; + + unsigned int degree = 0; // degree is the maximal degree of the components + + for (unsigned int i=0; i0) + { + multiplied_dofs_per_vertex += fes[i]->dofs_per_vertex * multiplicities[i]; + multiplied_dofs_per_line += fes[i]->dofs_per_line * multiplicities[i]; + multiplied_dofs_per_quad += fes[i]->dofs_per_quad * multiplicities[i]; + multiplied_dofs_per_hex += fes[i]->dofs_per_hex * multiplicities[i]; + + multiplied_n_components+=fes[i]->n_components() * multiplicities[i]; + + degree = std::max(degree, fes[i]->tensor_degree() ); + } + + // assume conformity of the first finite element and then take away + // bits as indicated by the base elements. if all multiplicities + // happen to be zero, then it doesn't matter what we set it to. + typename FiniteElementData::Conformity total_conformity + = typename FiniteElementData::Conformity(); + { + unsigned int index = 0; + for (index=0; index0) + { + total_conformity = fes[index]->conforming_space; + break; + } + + for (; index0) + total_conformity = + typename FiniteElementData::Conformity(total_conformity + & + fes[index]->conforming_space); + } + + std::vector dpo; + dpo.push_back(multiplied_dofs_per_vertex); + dpo.push_back(multiplied_dofs_per_line); + if (dim>1) dpo.push_back(multiplied_dofs_per_quad); + if (dim>2) dpo.push_back(multiplied_dofs_per_hex); + + return FiniteElementData (dpo, + multiplied_n_components, + degree, + total_conformity); + } + + /** + * Same as above but for a specific number of sub-elements. + */ + template + FiniteElementData + multiply_dof_numbers (const FiniteElement *fe1, + const unsigned int N1, + const FiniteElement *fe2=NULL, + const unsigned int N2=0, + const FiniteElement *fe3=NULL, + const unsigned int N3=0, + const FiniteElement *fe4=NULL, + const unsigned int N4=0, + const FiniteElement *fe5=NULL, + const unsigned int N5=0) + { + std::vector*> fes; + fes.push_back(fe1); + fes.push_back(fe2); + fes.push_back(fe3); + fes.push_back(fe4); + fes.push_back(fe5); + + std::vector mult; + mult.push_back(N1); + mult.push_back(N2); + mult.push_back(N3); + mult.push_back(N4); + mult.push_back(N5); + return multiply_dof_numbers(fes, mult); + } + + + + /** + * Compute the named flags for a list of finite elements with multiplicities + * given in the second argument. This function is called from all the above + * functions. + */ + template + std::vector + compute_restriction_is_additive_flags (const std::vector*> &fes, + const std::vector &multiplicities) + { + Assert (fes.size() == multiplicities.size(), ExcInternalError()); + + // first count the number of dofs and components that will emerge from the + // given FEs + unsigned int n_shape_functions = 0; + for (unsigned int i=0; i0) // check needed as fe might be NULL + n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i]; + + // generate the array that will hold the output + std::vector retval (n_shape_functions, false); + + // finally go through all the shape functions of the base elements, and copy + // their flags. this somehow copies the code in build_cell_table, which is + // not nice as it uses too much implicit knowledge about the layout of the + // individual bases in the composed FE, but there seems no way around... + // + // for each shape function, copy the flags from the base element to this + // one, taking into account multiplicities, and other complications + unsigned int total_index = 0; + for (unsigned int vertex_number=0; + vertex_number::vertices_per_cell; + ++vertex_number) + { + for (unsigned int base=0; basedofs_per_vertex; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_vertex*vertex_number + + local_index); + + Assert (index_in_base < fes[base]->dofs_per_cell, + ExcInternalError()); + retval[total_index] = fes[base]->restriction_is_additive(index_in_base); + } + } + + // 2. Lines + if (GeometryInfo::lines_per_cell > 0) + for (unsigned int line_number= 0; + line_number != GeometryInfo::lines_per_cell; + ++line_number) + { + for (unsigned int base=0; basedofs_per_line; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_line*line_number + + local_index + + fes[base]->first_line_index); + + Assert (index_in_base < fes[base]->dofs_per_cell, + ExcInternalError()); + retval[total_index] = fes[base]->restriction_is_additive(index_in_base); + } + } + + // 3. Quads + if (GeometryInfo::quads_per_cell > 0) + for (unsigned int quad_number= 0; + quad_number != GeometryInfo::quads_per_cell; + ++quad_number) + { + for (unsigned int base=0; basedofs_per_quad; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_quad*quad_number + + local_index + + fes[base]->first_quad_index); + + Assert (index_in_base < fes[base]->dofs_per_cell, + ExcInternalError()); + retval[total_index] = fes[base]->restriction_is_additive(index_in_base); + } + } + + // 4. Hexes + if (GeometryInfo::hexes_per_cell > 0) + for (unsigned int hex_number= 0; + hex_number != GeometryInfo::hexes_per_cell; + ++hex_number) + { + for (unsigned int base=0; basedofs_per_hex; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_hex*hex_number + + local_index + + fes[base]->first_hex_index); + + Assert (index_in_base < fes[base]->dofs_per_cell, + ExcInternalError()); + retval[total_index] = fes[base]->restriction_is_additive(index_in_base); + } + } + + Assert (total_index == n_shape_functions, ExcInternalError()); + + return retval; + } + + + + /** + * Take a @p FiniteElement object + * and return an boolean vector including the @p + * restriction_is_additive_flags of the mixed element consisting of @p N + * elements of the sub-element @p fe. + */ + template + std::vector + compute_restriction_is_additive_flags (const FiniteElement *fe1, + const unsigned int N1, + const FiniteElement *fe2=NULL, + const unsigned int N2=0, + const FiniteElement *fe3=NULL, + const unsigned int N3=0, + const FiniteElement *fe4=NULL, + const unsigned int N4=0, + const FiniteElement *fe5=NULL, + const unsigned int N5=0) + { + std::vector*> fe_list; + std::vector multiplicities; + + fe_list.push_back (fe1); + multiplicities.push_back (N1); + + fe_list.push_back (fe2); + multiplicities.push_back (N2); + + fe_list.push_back (fe3); + multiplicities.push_back (N3); + + fe_list.push_back (fe4); + multiplicities.push_back (N4); + + fe_list.push_back (fe5); + multiplicities.push_back (N5); + return compute_restriction_is_additive_flags (fe_list, multiplicities); + } + + + + /** + * Compute the nonzero components of a list of finite elements with + * multiplicities given in the second argument. + */ + template + std::vector + compute_nonzero_components (const std::vector*> &fes, + const std::vector &multiplicities) + { + Assert (fes.size() == multiplicities.size(), ExcInternalError()); + + // first count the number of dofs and components that will emerge from the + // given FEs + unsigned int n_shape_functions = 0; + for (unsigned int i=0; i0) //needed because fe might be NULL + n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i]; + + unsigned int n_components = 0; + for (unsigned int i=0; i0) //needed because fe might be NULL + n_components += fes[i]->n_components() * multiplicities[i]; + + // generate the array that will hold the output + std::vector > + retval (n_shape_functions, std::vector (n_components, false)); + + // finally go through all the shape functions of the base elements, and copy + // their flags. this somehow copies the code in build_cell_table, which is + // not nice as it uses too much implicit knowledge about the layout of the + // individual bases in the composed FE, but there seems no way around... + // + // for each shape function, copy the non-zero flags from the base element to + // this one, taking into account multiplicities, multiple components in base + // elements, and other complications + unsigned int 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_components()) + for (unsigned int local_index = 0; + local_index < fes[base]->dofs_per_vertex; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_vertex*vertex_number + + local_index); + + Assert (comp_start+fes[base]->n_components() <= + retval[total_index].size(), + ExcInternalError()); + for (unsigned int c=0; cn_components(); ++c) + { + Assert (c < fes[base]->get_nonzero_components(index_in_base).size(), + ExcInternalError()); + retval[total_index][comp_start+c] + = fes[base]->get_nonzero_components(index_in_base)[c]; + } + } + } + + // 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_components()) + for (unsigned int local_index = 0; + local_index < fes[base]->dofs_per_line; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_line*line_number + + local_index + + fes[base]->first_line_index); + + Assert (comp_start+fes[base]->n_components() <= + retval[total_index].size(), + ExcInternalError()); + for (unsigned int c=0; cn_components(); ++c) + { + Assert (c < fes[base]->get_nonzero_components(index_in_base).size(), + ExcInternalError()); + retval[total_index][comp_start+c] + = fes[base]->get_nonzero_components(index_in_base)[c]; + } + } + } + + // 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_components()) + for (unsigned int local_index = 0; + local_index < fes[base]->dofs_per_quad; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_quad*quad_number + + local_index + + fes[base]->first_quad_index); + + Assert (comp_start+fes[base]->n_components() <= + retval[total_index].size(), + ExcInternalError()); + for (unsigned int c=0; cn_components(); ++c) + { + Assert (c < fes[base]->get_nonzero_components(index_in_base).size(), + ExcInternalError()); + retval[total_index][comp_start+c] + = fes[base]->get_nonzero_components(index_in_base)[c]; + } + } + } + + // 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_components()) + for (unsigned int local_index = 0; + local_index < fes[base]->dofs_per_hex; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_hex*hex_number + + local_index + + fes[base]->first_hex_index); + + Assert (comp_start+fes[base]->n_components() <= + retval[total_index].size(), + ExcInternalError()); + for (unsigned int c=0; cn_components(); ++c) + { + Assert (c < fes[base]->get_nonzero_components(index_in_base).size(), + ExcInternalError()); + retval[total_index][comp_start+c] + = fes[base]->get_nonzero_components(index_in_base)[c]; + } + } + } + + Assert (total_index == n_shape_functions, ExcInternalError()); + + // now copy the vector > into a vector. + // this appears complicated but we do it this way since it's just + // awkward to generate ComponentMasks directly and so we need the + // recourse of the inner vector anyway. + std::vector xretval (retval.size()); + for (unsigned int i=0; i + std::vector + compute_nonzero_components (const FiniteElement *fe1, + const unsigned int N1, + const FiniteElement *fe2=NULL, + const unsigned int N2=0, + const FiniteElement *fe3=NULL, + const unsigned int N3=0, + const FiniteElement *fe4=NULL, + const unsigned int N4=0, + const FiniteElement *fe5=NULL, + const unsigned int N5=0) + { + std::vector*> fe_list; + std::vector multiplicities; + + fe_list.push_back (fe1); + multiplicities.push_back (N1); + + fe_list.push_back (fe2); + multiplicities.push_back (N2); + + fe_list.push_back (fe3); + multiplicities.push_back (N3); + + fe_list.push_back (fe4); + multiplicities.push_back (N4); + + fe_list.push_back (fe5); + multiplicities.push_back (N5); + + return compute_nonzero_components (fe_list, multiplicities); + } +} + + template FESystem::FESystem (const FiniteElement &fe, @@ -2401,476 +2886,6 @@ compare_for_face_domination (const FiniteElement &fe_other) const -template -FiniteElementData -FESystem::multiply_dof_numbers (const FiniteElement *fe1, - const unsigned int N1, - const FiniteElement *fe2, - const unsigned int N2, - const FiniteElement *fe3, - const unsigned int N3, - const FiniteElement *fe4, - const unsigned int N4, - const FiniteElement *fe5, - const unsigned int N5) -{ - std::vector*> fes; - fes.push_back(fe1); - fes.push_back(fe2); - fes.push_back(fe3); - fes.push_back(fe4); - fes.push_back(fe5); - - std::vector mult; - mult.push_back(N1); - mult.push_back(N2); - mult.push_back(N3); - mult.push_back(N4); - mult.push_back(N5); - return multiply_dof_numbers(fes, mult); -} - - - -template -FiniteElementData -FESystem::multiply_dof_numbers ( - const std::vector*> &fes, - const std::vector &multiplicities) -{ - Assert (fes.size() == multiplicities.size(), ExcDimensionMismatch (fes.size(), multiplicities.size())); - - unsigned int multiplied_dofs_per_vertex = 0; - unsigned int multiplied_dofs_per_line = 0; - unsigned int multiplied_dofs_per_quad = 0; - unsigned int multiplied_dofs_per_hex = 0; - - unsigned int multiplied_n_components = 0; - - unsigned int degree = 0; // degree is the maximal degree of the components - - for (unsigned int i=0; i0) - { - multiplied_dofs_per_vertex+=fes[i]->dofs_per_vertex * multiplicities[i]; - multiplied_dofs_per_line+=fes[i]->dofs_per_line * multiplicities[i]; - multiplied_dofs_per_quad+=fes[i]->dofs_per_quad * multiplicities[i]; - multiplied_dofs_per_hex+=fes[i]->dofs_per_hex * multiplicities[i]; - - multiplied_n_components+=fes[i]->n_components() * multiplicities[i]; - - degree = std::max(degree, fes[i]->tensor_degree() ); - } - - // assume conformity of the first finite element and then take away - // bits as indicated by the base elements. if all multiplicities - // happen to be zero, then it doesn't matter what we set it to. - typename FiniteElementData::Conformity total_conformity - = typename FiniteElementData::Conformity(); - { - unsigned int index = 0; - for (index=0; index0) - { - total_conformity = fes[index]->conforming_space; - break; - } - - for (; index0) - total_conformity = - typename FiniteElementData::Conformity(total_conformity - & - fes[index]->conforming_space); - } - - std::vector dpo; - dpo.push_back(multiplied_dofs_per_vertex); - dpo.push_back(multiplied_dofs_per_line); - if (dim>1) dpo.push_back(multiplied_dofs_per_quad); - if (dim>2) dpo.push_back(multiplied_dofs_per_hex); - - return FiniteElementData (dpo, - multiplied_n_components, - degree, - total_conformity); -} - - - -template -std::vector -FESystem::compute_restriction_is_additive_flags (const FiniteElement *fe1, - const unsigned int N1, - const FiniteElement *fe2, - const unsigned int N2, - const FiniteElement *fe3, - const unsigned int N3, - const FiniteElement *fe4, - const unsigned int N4, - const FiniteElement *fe5, - const unsigned int N5) -{ - std::vector*> fe_list; - std::vector multiplicities; - - fe_list.push_back (fe1); - multiplicities.push_back (N1); - - fe_list.push_back (fe2); - multiplicities.push_back (N2); - - fe_list.push_back (fe3); - multiplicities.push_back (N3); - - fe_list.push_back (fe4); - multiplicities.push_back (N4); - - fe_list.push_back (fe5); - multiplicities.push_back (N5); - return compute_restriction_is_additive_flags (fe_list, multiplicities); -} - - -template -std::vector -FESystem:: -compute_restriction_is_additive_flags (const std::vector*> &fes, - const std::vector &multiplicities) -{ - Assert (fes.size() == multiplicities.size(), ExcInternalError()); - - // first count the number of dofs and components that will emerge from the - // given FEs - unsigned int n_shape_functions = 0; - for (unsigned int i=0; i0) // check needed as fe might be NULL - n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i]; - - // generate the array that will hold the output - std::vector retval (n_shape_functions, false); - - // finally go through all the shape functions of the base elements, and copy - // their flags. this somehow copies the code in build_cell_table, which is - // not nice as it uses too much implicit knowledge about the layout of the - // individual bases in the composed FE, but there seems no way around... - // - // for each shape function, copy the flags from the base element to this - // one, taking into account multiplicities, and other complications - unsigned int total_index = 0; - for (unsigned int vertex_number=0; - vertex_number::vertices_per_cell; - ++vertex_number) - { - for (unsigned int base=0; basedofs_per_vertex; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_vertex*vertex_number + - local_index); - - Assert (index_in_base < fes[base]->dofs_per_cell, - ExcInternalError()); - retval[total_index] = fes[base]->restriction_is_additive(index_in_base); - } - } - - // 2. Lines - if (GeometryInfo::lines_per_cell > 0) - for (unsigned int line_number= 0; - line_number != GeometryInfo::lines_per_cell; - ++line_number) - { - for (unsigned int base=0; basedofs_per_line; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_line*line_number + - local_index + - fes[base]->first_line_index); - - Assert (index_in_base < fes[base]->dofs_per_cell, - ExcInternalError()); - retval[total_index] = fes[base]->restriction_is_additive(index_in_base); - } - } - - // 3. Quads - if (GeometryInfo::quads_per_cell > 0) - for (unsigned int quad_number= 0; - quad_number != GeometryInfo::quads_per_cell; - ++quad_number) - { - for (unsigned int base=0; basedofs_per_quad; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_quad*quad_number + - local_index + - fes[base]->first_quad_index); - - Assert (index_in_base < fes[base]->dofs_per_cell, - ExcInternalError()); - retval[total_index] = fes[base]->restriction_is_additive(index_in_base); - } - } - - // 4. Hexes - if (GeometryInfo::hexes_per_cell > 0) - for (unsigned int hex_number= 0; - hex_number != GeometryInfo::hexes_per_cell; - ++hex_number) - { - for (unsigned int base=0; basedofs_per_hex; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_hex*hex_number + - local_index + - fes[base]->first_hex_index); - - Assert (index_in_base < fes[base]->dofs_per_cell, - ExcInternalError()); - retval[total_index] = fes[base]->restriction_is_additive(index_in_base); - } - } - - Assert (total_index == n_shape_functions, ExcInternalError()); - - return retval; -} - - - -template -std::vector -FESystem::compute_nonzero_components (const FiniteElement *fe1, - const unsigned int N1, - const FiniteElement *fe2, - const unsigned int N2, - const FiniteElement *fe3, - const unsigned int N3, - const FiniteElement *fe4, - const unsigned int N4, - const FiniteElement *fe5, - const unsigned int N5) -{ - std::vector*> fe_list; - std::vector multiplicities; - - fe_list.push_back (fe1); - multiplicities.push_back (N1); - - fe_list.push_back (fe2); - multiplicities.push_back (N2); - - fe_list.push_back (fe3); - multiplicities.push_back (N3); - - fe_list.push_back (fe4); - multiplicities.push_back (N4); - - fe_list.push_back (fe5); - multiplicities.push_back (N5); - - return compute_nonzero_components (fe_list, multiplicities); -} - - - -template -std::vector -FESystem:: -compute_nonzero_components (const std::vector*> &fes, - const std::vector &multiplicities) -{ - Assert (fes.size() == multiplicities.size(), ExcInternalError()); - - // first count the number of dofs and components that will emerge from the - // given FEs - unsigned int n_shape_functions = 0; - for (unsigned int i=0; i0) //needed because fe might be NULL - n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i]; - - unsigned int n_components = 0; - for (unsigned int i=0; i0) //needed because fe might be NULL - n_components += fes[i]->n_components() * multiplicities[i]; - - // generate the array that will hold the output - std::vector > - retval (n_shape_functions, std::vector (n_components, false)); - - // finally go through all the shape functions of the base elements, and copy - // their flags. this somehow copies the code in build_cell_table, which is - // not nice as it uses too much implicit knowledge about the layout of the - // individual bases in the composed FE, but there seems no way around... - // - // for each shape function, copy the non-zero flags from the base element to - // this one, taking into account multiplicities, multiple components in base - // elements, and other complications - unsigned int 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_components()) - for (unsigned int local_index = 0; - local_index < fes[base]->dofs_per_vertex; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_vertex*vertex_number + - local_index); - - Assert (comp_start+fes[base]->n_components() <= - retval[total_index].size(), - ExcInternalError()); - for (unsigned int c=0; cn_components(); ++c) - { - Assert (index_in_base < fes[base]->nonzero_components.size(), - ExcInternalError()); - Assert (c < fes[base]->nonzero_components[index_in_base].size(), - ExcInternalError()); - retval[total_index][comp_start+c] - = fes[base]->nonzero_components[index_in_base][c]; - }; - } - } - - // 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_components()) - for (unsigned int local_index = 0; - local_index < fes[base]->dofs_per_line; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_line*line_number + - local_index + - fes[base]->first_line_index); - - Assert (comp_start+fes[base]->n_components() <= - retval[total_index].size(), - ExcInternalError()); - for (unsigned int c=0; cn_components(); ++c) - { - Assert (index_in_base < fes[base]->nonzero_components.size(), - ExcInternalError()); - Assert (c < fes[base]->nonzero_components[index_in_base].size(), - ExcInternalError()); - retval[total_index][comp_start+c] - = fes[base]->nonzero_components[index_in_base][c]; - }; - } - } - - // 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_components()) - for (unsigned int local_index = 0; - local_index < fes[base]->dofs_per_quad; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_quad*quad_number + - local_index + - fes[base]->first_quad_index); - - Assert (comp_start+fes[base]->n_components() <= - retval[total_index].size(), - ExcInternalError()); - for (unsigned int c=0; cn_components(); ++c) - { - Assert (index_in_base < fes[base]->nonzero_components.size(), - ExcInternalError()); - Assert (c < fes[base]->nonzero_components[index_in_base].size(), - ExcInternalError()); - retval[total_index][comp_start+c] - = fes[base]->nonzero_components[index_in_base][c]; - }; - } - } - - // 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_components()) - for (unsigned int local_index = 0; - local_index < fes[base]->dofs_per_hex; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_hex*hex_number + - local_index + - fes[base]->first_hex_index); - - Assert (comp_start+fes[base]->n_components() <= - retval[total_index].size(), - ExcInternalError()); - for (unsigned int c=0; cn_components(); ++c) - { - Assert (index_in_base < fes[base]->nonzero_components.size(), - ExcInternalError()); - Assert (c < fes[base]->nonzero_components[index_in_base].size(), - ExcInternalError()); - retval[total_index][comp_start+c] - = fes[base]->nonzero_components[index_in_base][c]; - }; - } - } - - Assert (total_index == n_shape_functions, ExcInternalError()); - - // now copy the vector > into a vector. - // this appears complicated but we do it this way since it's just - // awkward to generate ComponentMasks directly and so we need the - // recourse of the inner vector anyway. - std::vector xretval (retval.size()); - for (unsigned int i=0; i const FiniteElement & FESystem::base_element (const unsigned int index) const -- 2.39.5