From: Wolfgang Bangerth Date: Wed, 25 May 2016 16:24:23 +0000 (-0500) Subject: Move functions into a separate namespace. X-Git-Tag: v8.5.0-rc1~1005^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d995b84852ac4f84c12ecb9bbe9536733b8d868a;p=dealii.git Move functions into a separate namespace. Also provide a higher level documentation. --- diff --git a/doc/news/changes.h b/doc/news/changes.h index 87aa44cb90..f4c40859da 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -269,11 +269,11 @@ inconvenience this causes. (Bruno Turcksin, 2016/05/22) -
  • New: There are now additional functions in the FETools namespace that build +
  • New: There are now additional functions in the FETools::Compositing namespace that build finite elements out of simpler finite elements, either by forming tensor products or by combining the set of shape functions.
    - (Denis Davydov, 2016/05/20) + (Denis Davydov, Wolfgang Bangerth, 2016/05/20)
  • New: Added PArpackSolver::reinit() when dealing with BlockVectors. diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h index 98064474d9..0cd903b510 100644 --- a/include/deal.II/fe/fe_system.h +++ b/include/deal.II/fe/fe_system.h @@ -34,7 +34,10 @@ DEAL_II_NAMESPACE_OPEN * one. To the outside world, the resulting object looks just like a usual * finite element object, which is composed of several other finite elements * that are possibly of different type. The result is then a vector-valued - * finite element. %Vector valued elements are discussed in a number of + * finite element. An example is given in the documentation of namespace + * FETools::Compositing, when using the "tensor product" strategy. + * + * %Vector valued elements are discussed in a number of * tutorial programs, for example step-8, step-20, step-21, and in particular * in the * @ref vector_valued diff --git a/include/deal.II/fe/fe_tools.h b/include/deal.II/fe/fe_tools.h index 4c524b6cc4..9dba8de61e 100644 --- a/include/deal.II/fe/fe_tools.h +++ b/include/deal.II/fe/fe_tools.h @@ -841,148 +841,248 @@ namespace FETools /** - * 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 number of components - * returned in the FiniteElementData object is the sum over the - * product of the number of components in each of the finite - * elements times the corresponding multiplicity. Otherwise the - * number of components is taken from the first finite element with - * non-zero multiplicity, and all other elements with non-zero - * multiplicities need to have the same number of vector components. + * A namespace that contains functions that help build more + * complex finite elements from simpler ("base") elements. + * + * There are generally two ways in which one can build more complex + * elements, and this is reflected by several of the functions in + * this namespace having arguments called + * do_tensor_product: + * + *
      + *
    1. Tensor product construction (do_tensor_product=true): + * The tensor product construction, in the simplest case, builds a + * vector-valued element from scalar elements (see + * @ref vector_valued "this documentation module" and + * @ref GlossComponent "this glossary entry" for more information). + * To give an example, consider creating a vector-valued element with + * two vector components, where the first should have linear shape + * functions and the second quadratic shape functions. In 1d, the + * shape functions (on the reference cell) of the base elements are then + * @f{align*} + * Q_1 &= \{ 1-x, x \}, + * \\ Q_2 &= \{ 2(\frac 12 - x)(1-x), 2(x - \frac 12)x, 4x(1-x) \}, + * @f} + * where shape functions are ordered in the usual way (first on the + * first vertex, then on the second vertex, then in the interior of + * the cell). The tensor product construction will create an element with + * the following shape functions: + * @f{align*} + * Q_1 \times Q_2 &= + * \left\{ + * \begin{pmatrix} 1-x \\ 0 \end{pmatrix}, + * \begin{pmatrix} 0 \\ 2(\frac 12 - x)(1-x) \end{pmatrix}, + * \begin{pmatrix} x \\ 0 \end{pmatrix}, + * \begin{pmatrix} 0 \\ 2(x - \frac 12)x \end{pmatrix}, + * \begin{pmatrix} 0 \\ 4x(1-x) \end{pmatrix} + * \right\}. + * @f} + * The list here is again in standard order. + * + * Of course, the procedure also works if the base elements are + * already vector valued themselves: in that case, the composed + * element simply has as many vector components as the base elements + * taken together. + * + *
    2. Combining shape functions + * (do_tensor_product=false): In contrast to the + * previous strategy, combining shape functions simply takes + * all of the shape functions together. In the case above, + * this would yield the following element: + * @f{align*} + * Q_1 + Q_2 &= \{ 1-x, 2(\frac 12 - x)(1-x), + * x, 2(x - \frac 12)x, 4x(1-x) \}. + * @f} + * In other words, if the base elements are scalar, the resulting + * element will also be. In general, the base elements all will + * have to have the same number of vector components. + * + * The element constructed above of course no longer has a linearly + * independent set of shape functions. As a consequence, any matrix + * one creates by treating all shape functions the same will be + * singular. In practice, this strategy is therefore typically used + * in situations where one explicitly makes sure that certain shape + * functions are treated differently (e.g., by multiplying them with + * weight factors), or in cases where the shape functions one + * combines are not linearly dependent. + * + *
    */ - template - FiniteElementData - multiply_dof_numbers (const std::vector*> &fes, - const std::vector &multiplicities, - const bool do_tensor_product = true); + namespace Compositing + { - /** - * 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); + /** + * 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 number of components + * returned in the FiniteElementData object is the sum over the + * product of the number of components in each of the finite + * elements times the corresponding multiplicity. Otherwise the + * number of components is taken from the first finite element with + * non-zero multiplicity, and all other elements with non-zero + * multiplicities need to have the same number of vector components. + * + * See the documentation of namespace FETools::Compositing for more + * information about the @p do_tensor_product argument. + */ + template + FiniteElementData + multiply_dof_numbers (const std::vector*> &fes, + const std::vector &multiplicities, + const bool do_tensor_product = true); - /** - * Compute the "restriction is additive" flags (see the - * documentation of the FiniteElement class) for a list of finite - * elements with multiplicities given in the second argument. - */ - template - std::vector - compute_restriction_is_additive_flags (const std::vector*> &fes, - const std::vector &multiplicities); + /** + * 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); - /** - * Take a @p FiniteElement object and return a boolean vector - * describing the @p restriction_is_additive_flags (see the - * documentation of the FiniteElement class) for each shape function - * of the mixed element consisting of @p N1, @p N2, ... copies of - * the sub-elements @p fe1, @p fe2, ... - */ - 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 "restriction is additive" flags (see the + * documentation of the FiniteElement class) for a list of finite + * elements with multiplicities given in the second argument. + * + * The "restriction is additive" flags are properties of + * individual shape functions that do not depend on whether the + * composed element uses the tensor product or combination + * strategy outlined in the documentation of the + * FETools::Composition namespace. Consequently, this function + * does not have a @p do_tensor_product argument. + */ + template + std::vector + compute_restriction_is_additive_flags (const std::vector*> &fes, + const std::vector &multiplicities); - /** - * Compute the nonzero components for each shape function of a - * composed finite element described by a list of finite elements - * with multiplicities given in the second argument. - * - * If @p do_tensor_product is true, the number of components (and - * thus the size of the ComponentMask objects) is the sum over the - * product of the number of components in each of the finite - * elements times the corresponding multiplicity. Otherwise the - * number of components is taken from the first finite element with - * non-zero multiplicity, and all other elements with non-zero - * multiplicities need to have the same number of vector components. - */ - template - std::vector - compute_nonzero_components (const std::vector*> &fes, - const std::vector &multiplicities, - const bool do_tensor_product = true); + /** + * Take a @p FiniteElement object and return a boolean vector + * describing the @p restriction_is_additive_flags (see the + * documentation of the FiniteElement class) for each shape function + * of the mixed element consisting of @p N1, @p N2, ... copies of + * the sub-elements @p fe1, @p fe2, ... + * + * The "restriction is additive" flags are properties of + * individual shape functions that do not depend on whether the + * composed element uses the tensor product or combination + * strategy outlined in the documentation of the + * FETools::Composition namespace. Consequently, this function + * does not have a @p do_tensor_product argument. + */ + 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 non-zero vector components of a composed finite - * element. This function is similar to the previous one, except - * that the pointers indicate the elements to be composed, and the - * arguments @p N1, @p N2, ... the multiplicities. Null pointers - * indicate that an argument is to be skipped. - */ - 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); + /** + * Compute the nonzero components for each shape function of a + * composed finite element described by a list of finite elements + * with multiplicities given in the second argument. + * + * If @p do_tensor_product is true, the number of components (and + * thus the size of the ComponentMask objects) is the sum over the + * product of the number of components in each of the finite + * elements times the corresponding multiplicity. Otherwise the + * number of components is taken from the first finite element with + * non-zero multiplicity, and all other elements with non-zero + * multiplicities need to have the same number of vector components. + * + * See the documentation of namespace FETools::Compositing for more + * information about the @p do_tensor_product argument. + */ + template + std::vector + compute_nonzero_components (const std::vector*> &fes, + const std::vector &multiplicities, + const bool do_tensor_product = true); - /** - * 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 number of components - * used for the composite element is the sum over the - * product of the number of components in each of the finite - * elements times the corresponding multiplicity. Otherwise the - * number of components is taken from the first finite element with - * non-zero multiplicity, and all other elements with non-zero - * multiplicities need to have the same number of vector components. - */ - 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); + /** + * Compute the non-zero vector components of a composed finite + * element. This function is similar to the previous one, except + * that the pointers indicate the elements to be composed, and the + * arguments @p N1, @p N2, ... the multiplicities. Null pointers + * indicate that an argument is to be skipped. + */ + 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); + + /** + * 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 number of components + * used for the composite element is the sum over the + * product of the number of components in each of the finite + * elements times the corresponding multiplicity. Otherwise the + * number of components is taken from the first finite element with + * non-zero multiplicity, and all other elements with non-zero + * multiplicities need to have the same number of vector components. + * + * See the documentation of namespace FETools::Compositing for more + * information about the @p do_tensor_product argument. + */ + 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); + + /** + * For a given (composite) @p finite_element build @p face_system_to_base_table, + * and @p face_system_to_component_table. + * + * If @p do_tensor_product is true, the number of components + * used for the composite element is the sum over the + * product of the number of components in each of the finite + * elements times the corresponding multiplicity. Otherwise the + * number of components is taken from the first finite element with + * non-zero multiplicity, and all other elements with non-zero + * multiplicities need to have the same number of vector components. + * + * See the documentation of namespace FETools::Compositing for more + * information about the @p do_tensor_product argument. + */ + 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); + + } - /** - * For a given (composite) @p finite_element build @p face_system_to_base_table, - * and @p face_system_to_component_table. - * - * If @p do_tensor_product is true, the number of components - * used for the composite element is the sum over the - * product of the number of components in each of the finite - * elements times the corresponding multiplicity. Otherwise the - * number of components is taken from the first finite element with - * non-zero multiplicity, and all other elements with non-zero - * multiplicities need to have the same number of vector components. - */ - 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 diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc index 9f1ff5a4dc..f51088b4cf 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -115,9 +115,9 @@ const unsigned int FESystem::invalid_face_number; template FESystem::FESystem (const FiniteElement &fe, const unsigned int 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)), + FiniteElement (FETools::Compositing::multiply_dof_numbers(&fe, n_elements), + FETools::Compositing::compute_restriction_is_additive_flags (&fe, n_elements), + FETools::Compositing::compute_nonzero_components(&fe, n_elements)), base_elements((n_elements>0)) { std::vector*> fes; @@ -134,10 +134,10 @@ FESystem::FESystem (const FiniteElement &fe1, const unsigned int n1, const FiniteElement &fe2, const unsigned int n2) : - FiniteElement (FETools::multiply_dof_numbers(&fe1, n1, &fe2, n2), - FETools::compute_restriction_is_additive_flags (&fe1, n1, + FiniteElement (FETools::Compositing::multiply_dof_numbers(&fe1, n1, &fe2, n2), + FETools::Compositing::compute_restriction_is_additive_flags (&fe1, n1, &fe2, n2), - FETools::compute_nonzero_components(&fe1, n1, + FETools::Compositing::compute_nonzero_components(&fe1, n1, &fe2, n2)), base_elements((n1>0)+(n2>0)) { @@ -159,13 +159,13 @@ FESystem::FESystem (const FiniteElement &fe1, const unsigned int n2, const FiniteElement &fe3, const unsigned int n3) : - FiniteElement (FETools::multiply_dof_numbers(&fe1, n1, - &fe2, n2, - &fe3, n3), - FETools::compute_restriction_is_additive_flags (&fe1, n1, + FiniteElement (FETools::Compositing::multiply_dof_numbers(&fe1, n1, + &fe2, n2, + &fe3, n3), + FETools::Compositing::compute_restriction_is_additive_flags (&fe1, n1, &fe2, n2, &fe3, n3), - FETools::compute_nonzero_components(&fe1, n1, + FETools::Compositing::compute_nonzero_components(&fe1, n1, &fe2, n2, &fe3, n3)), base_elements((n1>0)+(n2>0)+(n3>0)) @@ -192,15 +192,15 @@ FESystem::FESystem (const FiniteElement &fe1, const unsigned int n3, const FiniteElement &fe4, const unsigned int n4) : - FiniteElement (FETools::multiply_dof_numbers(&fe1, n1, - &fe2, n2, - &fe3, n3, - &fe4, n4), - FETools::compute_restriction_is_additive_flags (&fe1, n1, + FiniteElement (FETools::Compositing::multiply_dof_numbers(&fe1, n1, + &fe2, n2, + &fe3, n3, + &fe4, n4), + FETools::Compositing::compute_restriction_is_additive_flags (&fe1, n1, &fe2, n2, &fe3, n3, &fe4, n4), - FETools::compute_nonzero_components(&fe1, n1, + FETools::Compositing::compute_nonzero_components(&fe1, n1, &fe2, n2, &fe3, n3, &fe4 ,n4)), @@ -232,17 +232,17 @@ FESystem::FESystem (const FiniteElement &fe1, const unsigned int n4, const FiniteElement &fe5, const unsigned int n5) : - FiniteElement (FETools::multiply_dof_numbers(&fe1, n1, - &fe2, n2, - &fe3, n3, - &fe4, n4, - &fe5, n5), - FETools::compute_restriction_is_additive_flags (&fe1, n1, + FiniteElement (FETools::Compositing::multiply_dof_numbers(&fe1, n1, + &fe2, n2, + &fe3, n3, + &fe4, n4, + &fe5, n5), + FETools::Compositing::compute_restriction_is_additive_flags (&fe1, n1, &fe2, n2, &fe3, n3, &fe4, n4, &fe5, n5), - FETools::compute_nonzero_components(&fe1, n1, + FETools::Compositing::compute_nonzero_components(&fe1, n1, &fe2, n2, &fe3, n3, &fe4 ,n4, @@ -271,9 +271,9 @@ FESystem::FESystem ( const std::vector*> &fes, const std::vector &multiplicities) : - FiniteElement (FETools::multiply_dof_numbers(fes, multiplicities), - FETools::compute_restriction_is_additive_flags (fes, multiplicities), - FETools::compute_nonzero_components(fes, multiplicities)), + FiniteElement (FETools::Compositing::multiply_dof_numbers(fes, multiplicities), + FETools::Compositing::compute_restriction_is_additive_flags (fes, multiplicities), + FETools::Compositing::compute_nonzero_components(fes, multiplicities)), base_elements(count_nonzeros(multiplicities)) { initialize(fes, multiplicities); @@ -1495,14 +1495,14 @@ void FESystem::initialize (const std::vectorsystem_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); + FETools::Compositing::build_cell_tables(this->system_to_base_table, + this->system_to_component_table, + this->component_to_base_table, + *this); - FETools::build_face_tables(this->face_system_to_base_table, - this->face_system_to_component_table, - *this); + FETools::Compositing::build_face_tables(this->face_system_to_base_table, + this->face_system_to_component_table, + *this); } diff --git a/source/fe/fe_tools.cc b/source/fe/fe_tools.cc index 7f87499c2b..4c5922b101 100644 --- a/source/fe/fe_tools.cc +++ b/source/fe/fe_tools.cc @@ -59,225 +59,159 @@ DEAL_II_NAMESPACE_OPEN namespace FETools { - template - FiniteElementData - multiply_dof_numbers (const std::vector*> &fes, - const std::vector &multiplicities, - const bool do_tensor_product) + namespace Compositing { - 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; - } + template + FiniteElementData + multiply_dof_numbers (const std::vector*> &fes, + const std::vector &multiplicities, + const bool do_tensor_product) + { + AssertDimension(fes.size(), multiplicities.size()); - 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]; + 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; - multiplied_n_components+=fes[i]->n_components() * multiplicities[i]; + unsigned int multiplied_n_components = 0; - Assert (do_tensor_product || (n_components == fes[i]->n_components()), - ExcDimensionMismatch(n_components, fes[i]->n_components())); + unsigned int degree = 0; // degree is the maximal degree of the 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) + unsigned int n_components = 0; + // Get the number of components from the first given finite element. + for (unsigned int i=0; i0) { - total_conformity = fes[index]->conforming_space; + n_components = fes[i]->n_components(); 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 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]; - 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); + multiplied_n_components+=fes[i]->n_components() * multiplicities[i]; - return FiniteElementData (dpo, - (do_tensor_product ? multiplied_n_components : n_components), - degree, - total_conformity, - block_indices); - } + Assert (do_tensor_product || (n_components == fes[i]->n_components()), + ExcDimensionMismatch(n_components, fes[i]->n_components())); - 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); - } + degree = std::max(degree, fes[i]->tensor_degree() ); + } - 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) + // 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(); { - 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); + unsigned int index = 0; + for (index=0; index0) + { + total_conformity = fes[index]->conforming_space; + break; + } - Assert (index_in_base < fes[base]->dofs_per_cell, - ExcInternalError()); - retval[total_index] = fes[base]->restriction_is_additive(index_in_base); - } + for (; index0) + total_conformity = + typename FiniteElementData::Conformity(total_conformity + & + fes[index]->conforming_space); } - // 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); + 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); - Assert (index_in_base < fes[base]->dofs_per_cell, - ExcInternalError()); - retval[total_index] = fes[base]->restriction_is_additive(index_in_base); - } - } + BlockIndices block_indices (0,0); - // 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); + 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); - Assert (index_in_base < fes[base]->dofs_per_cell, - ExcInternalError()); - retval[total_index] = fes[base]->restriction_is_additive(index_in_base); - } - } + return FiniteElementData (dpo, + (do_tensor_product ? multiplied_n_components : n_components), + degree, + total_conformity, + block_indices); + } - // 4. Hexes - if (GeometryInfo::hexes_per_cell > 0) - for (unsigned int hex_number= 0; - hex_number != GeometryInfo::hexes_per_cell; - ++hex_number) + 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_hex; + local_index < fes[base]->dofs_per_vertex; ++local_index, ++total_index) { const unsigned int index_in_base - = (fes[base]->dofs_per_hex*hex_number + - local_index + - fes[base]->first_hex_index); + = (fes[base]->dofs_per_vertex*vertex_number + + local_index); Assert (index_in_base < fes[base]->dofs_per_cell, ExcInternalError()); @@ -285,215 +219,188 @@ namespace FETools } } - Assert (total_index == n_shape_functions, ExcInternalError()); + // 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); - return retval; - } + 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); + } + } - /** - * 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; + // 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); - fe_list.push_back (fe1); - multiplicities.push_back (N1); + Assert (index_in_base < fes[base]->dofs_per_cell, + ExcInternalError()); + retval[total_index] = fes[base]->restriction_is_additive(index_in_base); + } + } - fe_list.push_back (fe2); - multiplicities.push_back (N2); + Assert (total_index == n_shape_functions, ExcInternalError()); - fe_list.push_back (fe3); - multiplicities.push_back (N3); + return retval; + } - 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); - } + /** + * 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); - template - std::vector - compute_nonzero_components (const std::vector*> &fes, - const std::vector &multiplicities, - const bool do_tensor_product) - { - AssertDimension(fes.size(), multiplicities.size()); + fe_list.push_back (fe2); + multiplicities.push_back (N2); - // 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]; + fe_list.push_back (fe3); + multiplicities.push_back (N3); - 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())); - } + fe_list.push_back (fe4); + multiplicities.push_back (N4); - // 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); + fe_list.push_back (fe5); + multiplicities.push_back (N5); + return compute_restriction_is_additive_flags (fe_list, multiplicities); + } - 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]; - } - } - } + template + std::vector + compute_nonzero_components (const std::vector*> &fes, + const std::vector &multiplicities, + const bool do_tensor_product) + { + AssertDimension(fes.size(), multiplicities.size()); - // 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); + // 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]; - 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]; - } - } + 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())); } - // 4. Hexes - if (GeometryInfo::hexes_per_cell > 0) - for (unsigned int hex_number= 0; - hex_number != GeometryInfo::hexes_per_cell; - ++hex_number) + // 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_hex; + local_index < fes[base]->dofs_per_vertex; ++local_index, ++total_index) { const unsigned int index_in_base - = (fes[base]->dofs_per_hex*hex_number + - local_index + - fes[base]->first_hex_index); + = (fes[base]->dofs_per_vertex*vertex_number + + local_index); Assert (comp_start+fes[base]->n_components() <= retval[total_index].size(), @@ -508,231 +415,213 @@ namespace FETools } } - Assert (total_index == n_shape_functions, ExcInternalError()); + // 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); - // 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; in_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); - /** - * 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, - 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; + 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]; + } + } + } - fe_list.push_back (fe1); - multiplicities.push_back (N1); + // 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); - fe_list.push_back (fe2); - multiplicities.push_back (N2); + 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]; + } + } + } - fe_list.push_back (fe3); - multiplicities.push_back (N3); + Assert (total_index == n_shape_functions, ExcInternalError()); - fe_list.push_back (fe4); - multiplicities.push_back (N4); + // 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; - 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; + fe_list.push_back (fe1); + multiplicities.push_back (N1); - if (do_tensor_product) - { - 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 + 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; - // 2. Lines - if (GeometryInfo::lines_per_cell > 0) - for (unsigned int line_number= 0; - line_number != GeometryInfo::lines_per_cell; - ++line_number) + if (do_tensor_product) { - 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) + else { - 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) + + // 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; base - 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_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::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_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::quads_per_face > 0) - for (unsigned int quad_number= 0; - quad_number != GeometryInfo::quads_per_face; - ++quad_number) + if (fe.base_element(base).is_primitive(index_in_base)) + { + const unsigned int comp_in_base + = fe.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 + = fe.base_element(base).system_to_component_index(index_in_base).second; + system_to_component_table[total_index] + = std::make_pair (comp, index_in_comp); + } + else + system_to_component_table[total_index] = non_primitive_index; + } + } + } + + 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 - 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, + namespace Compositing + \{ + 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 - 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 - 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 + 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 + 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 + 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 ( const FiniteElement & , @@ -105,7 +109,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS (const FiniteElement &, std::vector > > &, const bool, const double); #endif - \} + \} }