From 43da333198f4d5d0665d30d564c1e58954cd90ed Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Wed, 18 May 2016 16:56:09 +0200 Subject: [PATCH] move FESystem's auxiliary functions to FETools and extend them to cover the non-tensor product case of the enriched finite element --- include/deal.II/fe/fe_tools.h | 94 ++++++ source/fe/fe_system.cc | 567 +++------------------------------- source/fe/fe_tools.cc | 499 ++++++++++++++++++++++++++++++ source/fe/fe_tools.inst.in | 56 ++++ 4 files changed, 687 insertions(+), 529 deletions(-) diff --git a/include/deal.II/fe/fe_tools.h b/include/deal.II/fe/fe_tools.h index dab7a60736..d86bb1a2bc 100644 --- a/include/deal.II/fe/fe_tools.h +++ b/include/deal.II/fe/fe_tools.h @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -838,6 +839,99 @@ namespace FETools std::vector lexicographic_to_hierarchic_numbering (const FiniteElementData &fe_data); + + /** + * Take vectors of finite elements and multiplicities and multiply out + * how many degrees of freedom the composed element has per vertex, + * line, etc. + * + * If @p do_tensor_product is true, the returned number of components is sum of + * products of number of components in each finite elements times @p multiplicities. + * Otherwise the number of components is take from the first fine element + * with non-zero multiplicity. + */ + template + FiniteElementData + multiply_dof_numbers (const std::vector*> &fes, + const std::vector &multiplicities, + const bool do_tensor_product = true); + + /** + * 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); + + /** + * 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); + + /** + * 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); + + /** + * Compute the nonzero components of a list of finite elements with + * multiplicities given in the second argument. + * + * If @p do_tensor_product is false, the number of components of the resulting + * component mask is the same as the first finite element with non-zero multiplicity. + * Otherwise the number of components equals the sum of number of components + * of a each finite element times the multiplicity. + */ + template + std::vector + compute_nonzero_components (const std::vector*> &fes, + const std::vector &multiplicities, + const bool do_tensor_product = true); + + /** + * Compute the non-zero vector components of a composed finite element. + */ + template + 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); + /** * 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 2245f46400..a3af7aa6b7 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -22,6 +22,7 @@ #include #include #include +#include #include @@ -110,505 +111,13 @@ 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) - { - AssertDimension(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); - - BlockIndices block_indices (0,0); - - for (unsigned int base=0; base < fes.size(); ++base) - for (unsigned int m = 0; m < multiplicities[base]; ++m) - block_indices.push_back(fes[base]->dofs_per_cell); - - return FiniteElementData (dpo, - multiplied_n_components, - degree, - total_conformity, - block_indices); - } - - /** - * 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) - { - AssertDimension(fes.size(), multiplicities.size()); - - // 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) - { - AssertDimension(fes.size(), multiplicities.size()); - - // 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, const unsigned int n_elements) : - FiniteElement (multiply_dof_numbers(&fe, n_elements), - compute_restriction_is_additive_flags (&fe, n_elements), - compute_nonzero_components(&fe, n_elements)), + FiniteElement (FETools::multiply_dof_numbers(&fe, n_elements), + FETools::compute_restriction_is_additive_flags (&fe, n_elements), + FETools::compute_nonzero_components(&fe, n_elements)), base_elements((n_elements>0)) { std::vector*> fes; @@ -625,11 +134,11 @@ FESystem::FESystem (const FiniteElement &fe1, const unsigned int n1, const FiniteElement &fe2, const unsigned int n2) : - FiniteElement (multiply_dof_numbers(&fe1, n1, &fe2, n2), - compute_restriction_is_additive_flags (&fe1, n1, + FiniteElement (FETools::multiply_dof_numbers(&fe1, n1, &fe2, n2), + FETools::compute_restriction_is_additive_flags (&fe1, n1, &fe2, n2), - compute_nonzero_components(&fe1, n1, - &fe2, n2)), + FETools::compute_nonzero_components(&fe1, n1, + &fe2, n2)), base_elements((n1>0)+(n2>0)) { std::vector*> fes; @@ -650,15 +159,15 @@ FESystem::FESystem (const FiniteElement &fe1, const unsigned int n2, const FiniteElement &fe3, const unsigned int n3) : - FiniteElement (multiply_dof_numbers(&fe1, n1, - &fe2, n2, - &fe3, n3), - compute_restriction_is_additive_flags (&fe1, n1, + FiniteElement (FETools::multiply_dof_numbers(&fe1, n1, + &fe2, n2, + &fe3, n3), + FETools::compute_restriction_is_additive_flags (&fe1, n1, &fe2, n2, &fe3, n3), - compute_nonzero_components(&fe1, n1, - &fe2, n2, - &fe3, n3)), + FETools::compute_nonzero_components(&fe1, n1, + &fe2, n2, + &fe3, n3)), base_elements((n1>0)+(n2>0)+(n3>0)) { std::vector*> fes; @@ -683,18 +192,18 @@ FESystem::FESystem (const FiniteElement &fe1, const unsigned int n3, const FiniteElement &fe4, const unsigned int n4) : - FiniteElement (multiply_dof_numbers(&fe1, n1, - &fe2, n2, - &fe3, n3, - &fe4, n4), - compute_restriction_is_additive_flags (&fe1, n1, + FiniteElement (FETools::multiply_dof_numbers(&fe1, n1, + &fe2, n2, + &fe3, n3, + &fe4, n4), + FETools::compute_restriction_is_additive_flags (&fe1, n1, &fe2, n2, &fe3, n3, &fe4, n4), - compute_nonzero_components(&fe1, n1, - &fe2, n2, - &fe3, n3, - &fe4 ,n4)), + FETools::compute_nonzero_components(&fe1, n1, + &fe2, n2, + &fe3, n3, + &fe4 ,n4)), base_elements((n1>0)+(n2>0)+(n3>0)+(n4>0)) { std::vector*> fes; @@ -723,21 +232,21 @@ FESystem::FESystem (const FiniteElement &fe1, const unsigned int n4, const FiniteElement &fe5, const unsigned int n5) : - FiniteElement (multiply_dof_numbers(&fe1, n1, - &fe2, n2, - &fe3, n3, - &fe4, n4, - &fe5, n5), - compute_restriction_is_additive_flags (&fe1, n1, + FiniteElement (FETools::multiply_dof_numbers(&fe1, n1, + &fe2, n2, + &fe3, n3, + &fe4, n4, + &fe5, n5), + FETools::compute_restriction_is_additive_flags (&fe1, n1, &fe2, n2, &fe3, n3, &fe4, n4, &fe5, n5), - compute_nonzero_components(&fe1, n1, - &fe2, n2, - &fe3, n3, - &fe4 ,n4, - &fe5, n5)), + FETools::compute_nonzero_components(&fe1, n1, + &fe2, n2, + &fe3, n3, + &fe4 ,n4, + &fe5, n5)), base_elements((n1>0)+(n2>0)+(n3>0)+(n4>0)+(n5>0)) { std::vector*> fes; @@ -762,9 +271,9 @@ FESystem::FESystem ( const std::vector*> &fes, const std::vector &multiplicities) : - FiniteElement (multiply_dof_numbers(fes, multiplicities), - compute_restriction_is_additive_flags (fes, multiplicities), - compute_nonzero_components(fes, multiplicities)), + FiniteElement (FETools::multiply_dof_numbers(fes, multiplicities), + FETools::compute_restriction_is_additive_flags (fes, multiplicities), + FETools::compute_nonzero_components(fes, multiplicities)), base_elements(count_nonzeros(multiplicities)) { initialize(fes, multiplicities); diff --git a/source/fe/fe_tools.cc b/source/fe/fe_tools.cc index dc82e00ba6..06edebbab9 100644 --- a/source/fe/fe_tools.cc +++ b/source/fe/fe_tools.cc @@ -59,6 +59,505 @@ DEAL_II_NAMESPACE_OPEN namespace FETools { + template + FiniteElementData + multiply_dof_numbers (const std::vector*> &fes, + const std::vector &multiplicities, + const bool do_tensor_product) + { + AssertDimension(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 + + unsigned int n_components = 0; + // Get the number of components from the first given finite element. + for (unsigned int i=0; i0) + { + n_components = fes[i]->n_components(); + break; + } + + 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]; + + Assert (do_tensor_product || (n_components == fes[i]->n_components()), + ExcDimensionMismatch(n_components, fes[i]->n_components())); + + 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); + + BlockIndices block_indices (0,0); + + for (unsigned int base=0; base < fes.size(); ++base) + for (unsigned int m = 0; m < multiplicities[base]; ++m) + block_indices.push_back(fes[base]->dofs_per_cell); + + return FiniteElementData (dpo, + (do_tensor_product ? multiplied_n_components : n_components), + degree, + total_conformity, + block_indices); + } + + template + FiniteElementData + 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 + std::vector + compute_restriction_is_additive_flags (const std::vector*> &fes, + const std::vector &multiplicities) + { + AssertDimension(fes.size(), multiplicities.size()); + + // 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, + 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 + compute_nonzero_components (const std::vector*> &fes, + const std::vector &multiplicities, + const bool do_tensor_product) + { + AssertDimension(fes.size(), multiplicities.size()); + + // 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; + if (do_tensor_product) + { + for (unsigned int i=0; i0) //needed because fe might be NULL + n_components += fes[i]->n_components() * multiplicities[i]; + } + else + { + for (unsigned int i=0; i0) //needed because fe might be NULL + { + n_components = fes[i]->n_components(); + break; + } + // Now check that all FEs have the same number of components: + for (unsigned int i=0; i0) //needed because fe might be NULL + Assert (n_components == fes[i]->n_components(), + ExcDimensionMismatch(n_components,fes[i]->n_components())); + } + + // 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() * do_tensor_product) + 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() * do_tensor_product) + 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() * do_tensor_product) + 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() * do_tensor_product) + 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, + 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); + } + // Not implemented in the general case. template FiniteElement * diff --git a/source/fe/fe_tools.inst.in b/source/fe/fe_tools.inst.in index 8c6e7b4d14..90535f80ef 100644 --- a/source/fe/fe_tools.inst.in +++ b/source/fe/fe_tools.inst.in @@ -19,6 +19,62 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS namespace FETools \{ #if deal_II_dimension <= deal_II_space_dimension + template + FiniteElementData + multiply_dof_numbers (const std::vector*> &fes, + const std::vector &multiplicities, + bool); + + template + FiniteElementData + 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); + + template + std::vector + compute_restriction_is_additive_flags (const std::vector*> &fes, + const std::vector &multiplicities); + + template + std::vector + 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); + + template + std::vector + compute_nonzero_components (const std::vector*> &fes, + const std::vector &multiplicities, + const bool do_tensor_product); + + template + std::vector + 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); + template void compute_block_renumbering ( const FiniteElement & , -- 2.39.5