/**
- * 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
+ * <code>do_tensor_product</code>:
+ *
+ * <ol>
+ * <li> Tensor product construction (<code>do_tensor_product=true</code>):
+ * 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.
+ *
+ * <li> Combining shape functions
+ * (<code>do_tensor_product=false</code>): In contrast to the
+ * previous strategy, combining shape functions simply takes
+ * <i>all</i> 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.
+ *
+ * </ol>
*/
- template <int dim, int spacedim>
- FiniteElementData<dim>
- multiply_dof_numbers (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
- const std::vector<unsigned int> &multiplicities,
- const bool do_tensor_product = true);
+ namespace Compositing
+ {
- /**
- * Same as above but for a specific number of sub-elements.
- */
- template <int dim, int spacedim>
- FiniteElementData<dim>
- multiply_dof_numbers (const FiniteElement<dim,spacedim> *fe1,
- const unsigned int N1,
- const FiniteElement<dim,spacedim> *fe2=NULL,
- const unsigned int N2=0,
- const FiniteElement<dim,spacedim> *fe3=NULL,
- const unsigned int N3=0,
- const FiniteElement<dim,spacedim> *fe4=NULL,
- const unsigned int N4=0,
- const FiniteElement<dim,spacedim> *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 <int dim, int spacedim>
+ FiniteElementData<dim>
+ multiply_dof_numbers (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
+ const std::vector<unsigned int> &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 <int dim, int spacedim>
- std::vector<bool>
- compute_restriction_is_additive_flags (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
- const std::vector<unsigned int> &multiplicities);
+ /**
+ * Same as above but for a specific number of sub-elements.
+ */
+ template <int dim, int spacedim>
+ FiniteElementData<dim>
+ multiply_dof_numbers (const FiniteElement<dim,spacedim> *fe1,
+ const unsigned int N1,
+ const FiniteElement<dim,spacedim> *fe2=NULL,
+ const unsigned int N2=0,
+ const FiniteElement<dim,spacedim> *fe3=NULL,
+ const unsigned int N3=0,
+ const FiniteElement<dim,spacedim> *fe4=NULL,
+ const unsigned int N4=0,
+ const FiniteElement<dim,spacedim> *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 <int dim, int spacedim>
- std::vector<bool>
- compute_restriction_is_additive_flags (const FiniteElement<dim,spacedim> *fe1,
- const unsigned int N1,
- const FiniteElement<dim,spacedim> *fe2=NULL,
- const unsigned int N2=0,
- const FiniteElement<dim,spacedim> *fe3=NULL,
- const unsigned int N3=0,
- const FiniteElement<dim,spacedim> *fe4=NULL,
- const unsigned int N4=0,
- const FiniteElement<dim,spacedim> *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 <int dim, int spacedim>
+ std::vector<bool>
+ compute_restriction_is_additive_flags (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
+ const std::vector<unsigned int> &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 <int dim, int spacedim>
- std::vector<ComponentMask>
- compute_nonzero_components (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
- const std::vector<unsigned int> &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 <int dim, int spacedim>
+ std::vector<bool>
+ compute_restriction_is_additive_flags (const FiniteElement<dim,spacedim> *fe1,
+ const unsigned int N1,
+ const FiniteElement<dim,spacedim> *fe2=NULL,
+ const unsigned int N2=0,
+ const FiniteElement<dim,spacedim> *fe3=NULL,
+ const unsigned int N3=0,
+ const FiniteElement<dim,spacedim> *fe4=NULL,
+ const unsigned int N4=0,
+ const FiniteElement<dim,spacedim> *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 <int dim, int spacedim>
- std::vector<ComponentMask>
- compute_nonzero_components (const FiniteElement<dim,spacedim> *fe1,
- const unsigned int N1,
- const FiniteElement<dim,spacedim> *fe2=NULL,
- const unsigned int N2=0,
- const FiniteElement<dim,spacedim> *fe3=NULL,
- const unsigned int N3=0,
- const FiniteElement<dim,spacedim> *fe4=NULL,
- const unsigned int N4=0,
- const FiniteElement<dim,spacedim> *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 <int dim, int spacedim>
+ std::vector<ComponentMask>
+ compute_nonzero_components (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
+ const std::vector<unsigned int> &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 <int dim, int spacedim>
- 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<dim,spacedim> &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 <int dim, int spacedim>
+ std::vector<ComponentMask>
+ compute_nonzero_components (const FiniteElement<dim,spacedim> *fe1,
+ const unsigned int N1,
+ const FiniteElement<dim,spacedim> *fe2=NULL,
+ const unsigned int N2=0,
+ const FiniteElement<dim,spacedim> *fe3=NULL,
+ const unsigned int N3=0,
+ const FiniteElement<dim,spacedim> *fe4=NULL,
+ const unsigned int N4=0,
+ const FiniteElement<dim,spacedim> *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 <int dim, int spacedim>
+ 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<dim,spacedim> &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 <int dim, int spacedim>
+ void
+ build_face_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &face_system_to_base_table,
+ std::vector< std::pair< unsigned int, unsigned int > > &face_system_to_component_table,
+ const FiniteElement<dim,spacedim> &finite_element,
+ const bool do_tensor_product = true);
+
+ }
- /**
- * 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 <int dim, int spacedim>
- void
- build_face_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &face_system_to_base_table,
- std::vector< std::pair< unsigned int, unsigned int > > &face_system_to_component_table,
- const FiniteElement<dim,spacedim> &finite_element,
- const bool do_tensor_product = true);
/**
* Parse the name of a finite element and generate a finite element object
namespace FETools
{
- template <int dim, int spacedim>
- FiniteElementData<dim>
- multiply_dof_numbers (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
- const std::vector<unsigned int> &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; i<fes.size(); i++)
- if (multiplicities[i]>0)
- {
- n_components = fes[i]->n_components();
- break;
- }
+ template <int dim, int spacedim>
+ FiniteElementData<dim>
+ multiply_dof_numbers (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
+ const std::vector<unsigned int> &multiplicities,
+ const bool do_tensor_product)
+ {
+ AssertDimension(fes.size(), multiplicities.size());
- for (unsigned int i=0; i<fes.size(); i++)
- if (multiplicities[i]>0)
- {
- 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<dim>::Conformity total_conformity
- = typename FiniteElementData<dim>::Conformity();
- {
- unsigned int index = 0;
- for (index=0; index<fes.size(); ++index)
- if (multiplicities[index]>0)
+ unsigned int n_components = 0;
+ // Get the number of components from the first given finite element.
+ for (unsigned int i=0; i<fes.size(); i++)
+ if (multiplicities[i]>0)
{
- total_conformity = fes[index]->conforming_space;
+ n_components = fes[i]->n_components();
break;
}
- for (; index<fes.size(); ++index)
- if (multiplicities[index]>0)
- total_conformity =
- typename FiniteElementData<dim>::Conformity(total_conformity
- &
- fes[index]->conforming_space);
- }
-
- std::vector<unsigned int> 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; i<fes.size(); i++)
+ if (multiplicities[i]>0)
+ {
+ 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<dim> (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 <int dim, int spacedim>
- FiniteElementData<dim>
- multiply_dof_numbers (const FiniteElement<dim,spacedim> *fe1,
- const unsigned int N1,
- const FiniteElement<dim,spacedim> *fe2,
- const unsigned int N2,
- const FiniteElement<dim,spacedim> *fe3,
- const unsigned int N3,
- const FiniteElement<dim,spacedim> *fe4,
- const unsigned int N4,
- const FiniteElement<dim,spacedim> *fe5,
- const unsigned int N5)
- {
- std::vector<const FiniteElement<dim,spacedim>*> fes;
- fes.push_back(fe1);
- fes.push_back(fe2);
- fes.push_back(fe3);
- fes.push_back(fe4);
- fes.push_back(fe5);
-
- std::vector<unsigned int> 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 <int dim, int spacedim>
- std::vector<bool>
- compute_restriction_is_additive_flags (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
- const std::vector<unsigned int> &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; i<fes.size(); ++i)
- if (multiplicities[i]>0) // 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<bool> 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<GeometryInfo<dim>::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<dim>::Conformity total_conformity
+ = typename FiniteElementData<dim>::Conformity();
{
- for (unsigned int base=0; base<fes.size(); ++base)
- for (unsigned int m=0; m<multiplicities[base]; ++m)
- 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);
+ unsigned int index = 0;
+ for (index=0; index<fes.size(); ++index)
+ if (multiplicities[index]>0)
+ {
+ 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 (; index<fes.size(); ++index)
+ if (multiplicities[index]>0)
+ total_conformity =
+ typename FiniteElementData<dim>::Conformity(total_conformity
+ &
+ fes[index]->conforming_space);
}
- // 2. Lines
- if (GeometryInfo<dim>::lines_per_cell > 0)
- for (unsigned int line_number= 0;
- line_number != GeometryInfo<dim>::lines_per_cell;
- ++line_number)
- {
- for (unsigned int base=0; base<fes.size(); ++base)
- for (unsigned int m=0; m<multiplicities[base]; ++m)
- 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);
+ std::vector<unsigned int> 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<dim>::quads_per_cell > 0)
- for (unsigned int quad_number= 0;
- quad_number != GeometryInfo<dim>::quads_per_cell;
- ++quad_number)
- {
- for (unsigned int base=0; base<fes.size(); ++base)
- for (unsigned int m=0; m<multiplicities[base]; ++m)
- 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);
+ 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<dim> (dpo,
+ (do_tensor_product ? multiplied_n_components : n_components),
+ degree,
+ total_conformity,
+ block_indices);
+ }
- // 4. Hexes
- if (GeometryInfo<dim>::hexes_per_cell > 0)
- for (unsigned int hex_number= 0;
- hex_number != GeometryInfo<dim>::hexes_per_cell;
- ++hex_number)
+ template <int dim, int spacedim>
+ FiniteElementData<dim>
+ multiply_dof_numbers (const FiniteElement<dim,spacedim> *fe1,
+ const unsigned int N1,
+ const FiniteElement<dim,spacedim> *fe2,
+ const unsigned int N2,
+ const FiniteElement<dim,spacedim> *fe3,
+ const unsigned int N3,
+ const FiniteElement<dim,spacedim> *fe4,
+ const unsigned int N4,
+ const FiniteElement<dim,spacedim> *fe5,
+ const unsigned int N5)
+ {
+ std::vector<const FiniteElement<dim,spacedim>*> fes;
+ fes.push_back(fe1);
+ fes.push_back(fe2);
+ fes.push_back(fe3);
+ fes.push_back(fe4);
+ fes.push_back(fe5);
+
+ std::vector<unsigned int> 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 <int dim, int spacedim>
+ std::vector<bool>
+ compute_restriction_is_additive_flags (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
+ const std::vector<unsigned int> &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; i<fes.size(); ++i)
+ if (multiplicities[i]>0) // 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<bool> 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<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex_number)
{
for (unsigned int base=0; base<fes.size(); ++base)
for (unsigned int m=0; m<multiplicities[base]; ++m)
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 (index_in_base < fes[base]->dofs_per_cell,
ExcInternalError());
}
}
- Assert (total_index == n_shape_functions, ExcInternalError());
+ // 2. Lines
+ if (GeometryInfo<dim>::lines_per_cell > 0)
+ for (unsigned int line_number= 0;
+ line_number != GeometryInfo<dim>::lines_per_cell;
+ ++line_number)
+ {
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base]; ++m)
+ 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);
- 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<dim>::quads_per_cell > 0)
+ for (unsigned int quad_number= 0;
+ quad_number != GeometryInfo<dim>::quads_per_cell;
+ ++quad_number)
+ {
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base]; ++m)
+ 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 (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 <int dim, int spacedim>
- std::vector<bool>
- compute_restriction_is_additive_flags (const FiniteElement<dim,spacedim> *fe1,
- const unsigned int N1,
- const FiniteElement<dim,spacedim> *fe2,
- const unsigned int N2,
- const FiniteElement<dim,spacedim> *fe3,
- const unsigned int N3,
- const FiniteElement<dim,spacedim> *fe4,
- const unsigned int N4,
- const FiniteElement<dim,spacedim> *fe5,
- const unsigned int N5)
- {
- std::vector<const FiniteElement<dim,spacedim>*> fe_list;
- std::vector<unsigned int> multiplicities;
+ // 4. Hexes
+ if (GeometryInfo<dim>::hexes_per_cell > 0)
+ for (unsigned int hex_number= 0;
+ hex_number != GeometryInfo<dim>::hexes_per_cell;
+ ++hex_number)
+ {
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base]; ++m)
+ 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 (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 <int dim, int spacedim>
+ std::vector<bool>
+ compute_restriction_is_additive_flags (const FiniteElement<dim,spacedim> *fe1,
+ const unsigned int N1,
+ const FiniteElement<dim,spacedim> *fe2,
+ const unsigned int N2,
+ const FiniteElement<dim,spacedim> *fe3,
+ const unsigned int N3,
+ const FiniteElement<dim,spacedim> *fe4,
+ const unsigned int N4,
+ const FiniteElement<dim,spacedim> *fe5,
+ const unsigned int N5)
+ {
+ std::vector<const FiniteElement<dim,spacedim>*> fe_list;
+ std::vector<unsigned int> multiplicities;
+ fe_list.push_back (fe1);
+ multiplicities.push_back (N1);
- template <int dim, int spacedim>
- std::vector<ComponentMask>
- compute_nonzero_components (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
- const std::vector<unsigned int> &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; i<fes.size(); ++i)
- if (multiplicities[i]>0) //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; i<fes.size(); ++i)
- if (multiplicities[i]>0) //needed because fe might be NULL
- n_components += fes[i]->n_components() * multiplicities[i];
- }
- else
- {
- for (unsigned int i=0; i<fes.size(); ++i)
- if (multiplicities[i]>0) //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; i<fes.size(); ++i)
- if (multiplicities[i]>0) //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<std::vector<bool> >
- retval (n_shape_functions, std::vector<bool> (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<GeometryInfo<dim>::vertices_per_cell;
- ++vertex_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fes.size(); ++base)
- for (unsigned int m=0; m<multiplicities[base];
- ++m, comp_start+=fes[base]->n_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; c<fes[base]->n_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<dim>::lines_per_cell > 0)
- for (unsigned int line_number= 0;
- line_number != GeometryInfo<dim>::lines_per_cell;
- ++line_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fes.size(); ++base)
- for (unsigned int m=0; m<multiplicities[base];
- ++m, comp_start+=fes[base]->n_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; c<fes[base]->n_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 <int dim, int spacedim>
+ std::vector<ComponentMask>
+ compute_nonzero_components (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
+ const std::vector<unsigned int> &multiplicities,
+ const bool do_tensor_product)
+ {
+ AssertDimension(fes.size(), multiplicities.size());
- // 3. Quads
- if (GeometryInfo<dim>::quads_per_cell > 0)
- for (unsigned int quad_number= 0;
- quad_number != GeometryInfo<dim>::quads_per_cell;
- ++quad_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fes.size(); ++base)
- for (unsigned int m=0; m<multiplicities[base];
- ++m, comp_start+=fes[base]->n_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; i<fes.size(); ++i)
+ if (multiplicities[i]>0) //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; c<fes[base]->n_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; i<fes.size(); ++i)
+ if (multiplicities[i]>0) //needed because fe might be NULL
+ n_components += fes[i]->n_components() * multiplicities[i];
+ }
+ else
+ {
+ for (unsigned int i=0; i<fes.size(); ++i)
+ if (multiplicities[i]>0) //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; i<fes.size(); ++i)
+ if (multiplicities[i]>0) //needed because fe might be NULL
+ Assert (n_components == fes[i]->n_components(),
+ ExcDimensionMismatch(n_components,fes[i]->n_components()));
}
- // 4. Hexes
- if (GeometryInfo<dim>::hexes_per_cell > 0)
- for (unsigned int hex_number= 0;
- hex_number != GeometryInfo<dim>::hexes_per_cell;
- ++hex_number)
+ // generate the array that will hold the output
+ std::vector<std::vector<bool> >
+ retval (n_shape_functions, std::vector<bool> (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<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex_number)
{
unsigned int comp_start = 0;
for (unsigned int base=0; base<fes.size(); ++base)
for (unsigned int m=0; m<multiplicities[base];
++m, comp_start+=fes[base]->n_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(),
}
}
- Assert (total_index == n_shape_functions, ExcInternalError());
+ // 2. Lines
+ if (GeometryInfo<dim>::lines_per_cell > 0)
+ for (unsigned int line_number= 0;
+ line_number != GeometryInfo<dim>::lines_per_cell;
+ ++line_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base];
+ ++m, comp_start+=fes[base]->n_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<vector<bool> > into a vector<ComponentMask>.
- // 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<bool> anyway.
- std::vector<ComponentMask> xretval (retval.size());
- for (unsigned int i=0; i<retval.size(); ++i)
- xretval[i] = ComponentMask(retval[i]);
- return xretval;
- }
+ Assert (comp_start+fes[base]->n_components() <=
+ retval[total_index].size(),
+ ExcInternalError());
+ for (unsigned int c=0; c<fes[base]->n_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<dim>::quads_per_cell > 0)
+ for (unsigned int quad_number= 0;
+ quad_number != GeometryInfo<dim>::quads_per_cell;
+ ++quad_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base];
+ ++m, comp_start+=fes[base]->n_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 <int dim, int spacedim>
- std::vector<ComponentMask>
- compute_nonzero_components (const FiniteElement<dim,spacedim> *fe1,
- const unsigned int N1,
- const FiniteElement<dim,spacedim> *fe2,
- const unsigned int N2,
- const FiniteElement<dim,spacedim> *fe3,
- const unsigned int N3,
- const FiniteElement<dim,spacedim> *fe4,
- const unsigned int N4,
- const FiniteElement<dim,spacedim> *fe5,
- const unsigned int N5)
- {
- std::vector<const FiniteElement<dim,spacedim>*> fe_list;
- std::vector<unsigned int> multiplicities;
+ Assert (comp_start+fes[base]->n_components() <=
+ retval[total_index].size(),
+ ExcInternalError());
+ for (unsigned int c=0; c<fes[base]->n_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<dim>::hexes_per_cell > 0)
+ for (unsigned int hex_number= 0;
+ hex_number != GeometryInfo<dim>::hexes_per_cell;
+ ++hex_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base];
+ ++m, comp_start+=fes[base]->n_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; c<fes[base]->n_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<vector<bool> > into a vector<ComponentMask>.
+ // 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<bool> anyway.
+ std::vector<ComponentMask> xretval (retval.size());
+ for (unsigned int i=0; i<retval.size(); ++i)
+ xretval[i] = ComponentMask(retval[i]);
+ return xretval;
+ }
- fe_list.push_back (fe5);
- multiplicities.push_back (N5);
- return compute_nonzero_components (fe_list, multiplicities);
- }
+ /**
+ * Compute the non-zero vector components of a composed finite element.
+ */
+ template <int dim, int spacedim>
+ std::vector<ComponentMask>
+ compute_nonzero_components (const FiniteElement<dim,spacedim> *fe1,
+ const unsigned int N1,
+ const FiniteElement<dim,spacedim> *fe2,
+ const unsigned int N2,
+ const FiniteElement<dim,spacedim> *fe3,
+ const unsigned int N3,
+ const FiniteElement<dim,spacedim> *fe4,
+ const unsigned int N4,
+ const FiniteElement<dim,spacedim> *fe5,
+ const unsigned int N5)
+ {
+ std::vector<const FiniteElement<dim,spacedim>*> fe_list;
+ std::vector<unsigned int> multiplicities;
- template <int dim, int spacedim>
- 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<dim,spacedim> &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<fe.base_element(base).n_components(); ++k)
- component_to_base_table[total_index++]
- = std::make_pair(std::make_pair(base,k), m);
- }
- Assert (total_index == component_to_base_table.size(),
- ExcInternalError());
- }
- else
- {
- // The base element establishing a component does not make sense in this case.
- // Set up to something meaningless:
- for (unsigned int i = 0; i < component_to_base_table.size(); i++)
- component_to_base_table[i] = std::make_pair(std::make_pair(numbers::invalid_unsigned_int,numbers::invalid_unsigned_int), numbers::invalid_unsigned_int);
+ 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);
- // 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<unsigned int, unsigned int>
- 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<GeometryInfo<dim>::vertices_per_cell;
- ++vertex_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fe.n_base_elements(); ++base)
- for (unsigned int m=0; m<fe.element_multiplicity(base);
- ++m, comp_start+=fe.base_element(base).n_components() * do_tensor_product)
- for (unsigned int local_index = 0;
- local_index < fe.base_element(base).dofs_per_vertex;
- ++local_index, ++total_index)
- {
- const unsigned int index_in_base
- = (fe.base_element(base).dofs_per_vertex*vertex_number +
- local_index);
+ fe_list.push_back (fe5);
+ multiplicities.push_back (N5);
- system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base, m), index_in_base);
+ return compute_nonzero_components (fe_list, multiplicities);
+ }
- 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 <int dim, int spacedim>
+ 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<dim,spacedim> &fe,
+ const bool do_tensor_product)
+ {
+ unsigned int total_index = 0;
- // 2. Lines
- if (GeometryInfo<dim>::lines_per_cell > 0)
- for (unsigned int line_number= 0;
- line_number != GeometryInfo<dim>::lines_per_cell;
- ++line_number)
+ if (do_tensor_product)
{
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fe.n_base_elements(); ++base)
- for (unsigned int m=0; m<fe.element_multiplicity(base);
- ++m, comp_start+=fe.base_element(base).n_components() * do_tensor_product)
- for (unsigned int local_index = 0;
- local_index < fe.base_element(base).dofs_per_line;
- ++local_index, ++total_index)
- {
- const unsigned int index_in_base
- = (fe.base_element(base).dofs_per_line*line_number +
- local_index +
- fe.base_element(base).first_line_index);
-
- system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base,m), index_in_base);
-
- 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;
- }
+ 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<fe.base_element(base).n_components(); ++k)
+ component_to_base_table[total_index++]
+ = std::make_pair(std::make_pair(base,k), m);
+ }
+ Assert (total_index == component_to_base_table.size(),
+ ExcInternalError());
}
-
- // 3. Quads
- if (GeometryInfo<dim>::quads_per_cell > 0)
- for (unsigned int quad_number= 0;
- quad_number != GeometryInfo<dim>::quads_per_cell;
- ++quad_number)
+ else
{
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fe.n_base_elements(); ++base)
- for (unsigned int m=0; m<fe.element_multiplicity(base);
- ++m, comp_start += fe.base_element(base).n_components() * do_tensor_product)
- for (unsigned int local_index = 0;
- local_index < fe.base_element(base).dofs_per_quad;
- ++local_index, ++total_index)
- {
- const unsigned int index_in_base
- = (fe.base_element(base).dofs_per_quad*quad_number +
- local_index +
- fe.base_element(base).first_quad_index);
+ // The base element establishing a component does not make sense in this case.
+ // Set up to something meaningless:
+ for (unsigned int i = 0; i < component_to_base_table.size(); i++)
+ component_to_base_table[i] = std::make_pair(std::make_pair(numbers::invalid_unsigned_int,numbers::invalid_unsigned_int), numbers::invalid_unsigned_int);
- system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base,m), index_in_base);
-
- 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;
- }
}
- // 4. Hexes
- if (GeometryInfo<dim>::hexes_per_cell > 0)
- for (unsigned int hex_number= 0;
- hex_number != GeometryInfo<dim>::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<unsigned int, unsigned int>
+ 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<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex_number)
{
unsigned int comp_start = 0;
for (unsigned int base=0; base<fe.n_base_elements(); ++base)
for (unsigned int m=0; m<fe.element_multiplicity(base);
++m, comp_start+=fe.base_element(base).n_components() * do_tensor_product)
for (unsigned int local_index = 0;
- local_index < fe.base_element(base).dofs_per_hex;
+ local_index < fe.base_element(base).dofs_per_vertex;
++local_index, ++total_index)
{
const unsigned int index_in_base
- = (fe.base_element(base).dofs_per_hex*hex_number +
- local_index +
- fe.base_element(base).first_hex_index);
+ = (fe.base_element(base).dofs_per_vertex*vertex_number +
+ local_index);
system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base,m), index_in_base);
+ = std::make_pair (std::make_pair(base, m), index_in_base);
if (fe.base_element(base).is_primitive(index_in_base))
{
system_to_component_table[total_index] = non_primitive_index;
}
}
- }
- template <int dim, int spacedim>
- void
- build_face_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &face_system_to_base_table,
- std::vector< std::pair< unsigned int, unsigned int > > &face_system_to_component_table,
- const FiniteElement<dim,spacedim> &fe,
- const bool do_tensor_product)
- {
- // Initialize index tables. do this in the same way as done for the cell
- // tables, except that we now loop over the objects of faces
-
- // For non-primitive shape functions, have a special invalid index
- const std::pair<unsigned int, unsigned int>
- non_primitive_index (numbers::invalid_unsigned_int,
- numbers::invalid_unsigned_int);
-
- // 1. Vertices
- unsigned int total_index = 0;
- for (unsigned int vertex_number=0;
- vertex_number<GeometryInfo<dim>::vertices_per_face;
- ++vertex_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fe.n_base_elements(); ++base)
- for (unsigned int m=0; m<fe.element_multiplicity(base);
- ++m, comp_start += fe.base_element(base).n_components() * do_tensor_product)
- for (unsigned int local_index = 0;
- local_index < fe.base_element(base).dofs_per_vertex;
- ++local_index, ++total_index)
- {
- // get (cell) index of this shape function inside the base
- // element to see whether the shape function is primitive
- // (assume that all shape functions on vertices share the same
- // primitivity property; assume likewise for all shape functions
- // located on lines, quads, etc. this way, we can ask for
- // primitivity of only _one_ shape function, which is taken as
- // representative for all others located on the same type of
- // object):
- const unsigned int index_in_base
- = (fe.base_element(base).dofs_per_vertex*vertex_number +
- local_index);
-
- const unsigned int face_index_in_base
- = (fe.base_element(base).dofs_per_vertex*vertex_number +
- local_index);
-
- face_system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base,m), face_index_in_base);
-
- if (fe.base_element(base).is_primitive(index_in_base))
+ // 2. Lines
+ if (GeometryInfo<dim>::lines_per_cell > 0)
+ for (unsigned int line_number= 0;
+ line_number != GeometryInfo<dim>::lines_per_cell;
+ ++line_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fe.n_base_elements(); ++base)
+ for (unsigned int m=0; m<fe.element_multiplicity(base);
+ ++m, comp_start+=fe.base_element(base).n_components() * do_tensor_product)
+ for (unsigned int local_index = 0;
+ local_index < fe.base_element(base).dofs_per_line;
+ ++local_index, ++total_index)
{
- const unsigned int comp_in_base
- = fe.base_element(base).face_system_to_component_index(face_index_in_base).first;
- const unsigned int comp
- = comp_start + comp_in_base;
- const unsigned int face_index_in_comp
- = fe.base_element(base).face_system_to_component_index(face_index_in_base).second;
- face_system_to_component_table[total_index]
- = std::make_pair (comp, face_index_in_comp);
+ const unsigned int index_in_base
+ = (fe.base_element(base).dofs_per_line*line_number +
+ local_index +
+ fe.base_element(base).first_line_index);
+
+ system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), index_in_base);
+
+ 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;
}
- else
- face_system_to_component_table[total_index] = non_primitive_index;
- }
- }
+ }
- // 2. Lines
- if (GeometryInfo<dim>::lines_per_face > 0)
- for (unsigned int line_number= 0;
- line_number != GeometryInfo<dim>::lines_per_face;
- ++line_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
- for (unsigned int m=0; m<fe.element_multiplicity(base);
- ++m, comp_start += fe.base_element(base).n_components() * do_tensor_product)
- for (unsigned int local_index = 0;
- local_index < fe.base_element(base).dofs_per_line;
- ++local_index, ++total_index)
- {
- // do everything alike for this type of object
- const unsigned int index_in_base
- = (fe.base_element(base).dofs_per_line*line_number +
- local_index +
- fe.base_element(base).first_line_index);
+ // 3. Quads
+ if (GeometryInfo<dim>::quads_per_cell > 0)
+ for (unsigned int quad_number= 0;
+ quad_number != GeometryInfo<dim>::quads_per_cell;
+ ++quad_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fe.n_base_elements(); ++base)
+ for (unsigned int m=0; m<fe.element_multiplicity(base);
+ ++m, comp_start += fe.base_element(base).n_components() * do_tensor_product)
+ for (unsigned int local_index = 0;
+ local_index < fe.base_element(base).dofs_per_quad;
+ ++local_index, ++total_index)
+ {
+ const unsigned int index_in_base
+ = (fe.base_element(base).dofs_per_quad*quad_number +
+ local_index +
+ fe.base_element(base).first_quad_index);
- const unsigned int face_index_in_base
- = (fe.base_element(base).first_face_line_index +
- fe.base_element(base).dofs_per_line * line_number +
- local_index);
+ system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), index_in_base);
- face_system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base,m), face_index_in_base);
+ if (fe.base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base
+ = fe.base_element(base).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;
+ }
+ }
- if (fe.base_element(base).is_primitive(index_in_base))
- {
- const unsigned int comp_in_base
- = fe.base_element(base).face_system_to_component_index(face_index_in_base).first;
- const unsigned int comp
- = comp_start + comp_in_base;
- const unsigned int face_index_in_comp
- = fe.base_element(base).face_system_to_component_index(face_index_in_base).second;
- face_system_to_component_table[total_index]
- = std::make_pair (comp, face_index_in_comp);
- }
- else
- face_system_to_component_table[total_index] = non_primitive_index;
- }
- }
+ // 4. Hexes
+ if (GeometryInfo<dim>::hexes_per_cell > 0)
+ for (unsigned int hex_number= 0;
+ hex_number != GeometryInfo<dim>::hexes_per_cell;
+ ++hex_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fe.n_base_elements(); ++base)
+ for (unsigned int m=0; m<fe.element_multiplicity(base);
+ ++m, comp_start+=fe.base_element(base).n_components() * do_tensor_product)
+ for (unsigned int local_index = 0;
+ local_index < fe.base_element(base).dofs_per_hex;
+ ++local_index, ++total_index)
+ {
+ const unsigned int index_in_base
+ = (fe.base_element(base).dofs_per_hex*hex_number +
+ local_index +
+ fe.base_element(base).first_hex_index);
+
+ system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), index_in_base);
- // 3. Quads
- if (GeometryInfo<dim>::quads_per_face > 0)
- for (unsigned int quad_number= 0;
- quad_number != GeometryInfo<dim>::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 <int dim, int spacedim>
+ void
+ build_face_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &face_system_to_base_table,
+ std::vector< std::pair< unsigned int, unsigned int > > &face_system_to_component_table,
+ const FiniteElement<dim,spacedim> &fe,
+ const bool do_tensor_product)
+ {
+ // Initialize index tables. do this in the same way as done for the cell
+ // tables, except that we now loop over the objects of faces
+
+ // For non-primitive shape functions, have a special invalid index
+ const std::pair<unsigned int, unsigned int>
+ non_primitive_index (numbers::invalid_unsigned_int,
+ numbers::invalid_unsigned_int);
+
+ // 1. Vertices
+ unsigned int total_index = 0;
+ for (unsigned int vertex_number=0;
+ vertex_number<GeometryInfo<dim>::vertices_per_face;
+ ++vertex_number)
{
unsigned int comp_start = 0;
for (unsigned int base=0; base<fe.n_base_elements(); ++base)
for (unsigned int m=0; m<fe.element_multiplicity(base);
++m, comp_start += fe.base_element(base).n_components() * do_tensor_product)
for (unsigned int local_index = 0;
- local_index < fe.base_element(base).dofs_per_quad;
+ local_index < fe.base_element(base).dofs_per_vertex;
++local_index, ++total_index)
{
- // do everything alike for this type of object
+ // get (cell) index of this shape function inside the base
+ // element to see whether the shape function is primitive
+ // (assume that all shape functions on vertices share the same
+ // primitivity property; assume likewise for all shape functions
+ // located on lines, quads, etc. this way, we can ask for
+ // primitivity of only _one_ shape function, which is taken as
+ // representative for all others located on the same type of
+ // object):
const unsigned int index_in_base
- = (fe.base_element(base).dofs_per_quad*quad_number +
- local_index +
- fe.base_element(base).first_quad_index);
+ = (fe.base_element(base).dofs_per_vertex*vertex_number +
+ local_index);
const unsigned int face_index_in_base
- = (fe.base_element(base).first_face_quad_index +
- fe.base_element(base).dofs_per_quad * quad_number +
+ = (fe.base_element(base).dofs_per_vertex*vertex_number +
local_index);
face_system_to_base_table[total_index]
face_system_to_component_table[total_index] = non_primitive_index;
}
}
- Assert (total_index == fe.dofs_per_face, ExcInternalError());
- Assert (total_index == face_system_to_component_table.size(),
- ExcInternalError());
- Assert (total_index == face_system_to_base_table.size(),
- ExcInternalError());
+
+ // 2. Lines
+ if (GeometryInfo<dim>::lines_per_face > 0)
+ for (unsigned int line_number= 0;
+ line_number != GeometryInfo<dim>::lines_per_face;
+ ++line_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
+ for (unsigned int m=0; m<fe.element_multiplicity(base);
+ ++m, comp_start += fe.base_element(base).n_components() * do_tensor_product)
+ for (unsigned int local_index = 0;
+ local_index < fe.base_element(base).dofs_per_line;
+ ++local_index, ++total_index)
+ {
+ // do everything alike for this type of object
+ const unsigned int index_in_base
+ = (fe.base_element(base).dofs_per_line*line_number +
+ local_index +
+ fe.base_element(base).first_line_index);
+
+ const unsigned int face_index_in_base
+ = (fe.base_element(base).first_face_line_index +
+ fe.base_element(base).dofs_per_line * line_number +
+ local_index);
+
+ face_system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), face_index_in_base);
+
+ if (fe.base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base
+ = fe.base_element(base).face_system_to_component_index(face_index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int face_index_in_comp
+ = fe.base_element(base).face_system_to_component_index(face_index_in_base).second;
+ face_system_to_component_table[total_index]
+ = std::make_pair (comp, face_index_in_comp);
+ }
+ else
+ face_system_to_component_table[total_index] = non_primitive_index;
+ }
+ }
+
+ // 3. Quads
+ if (GeometryInfo<dim>::quads_per_face > 0)
+ for (unsigned int quad_number= 0;
+ quad_number != GeometryInfo<dim>::quads_per_face;
+ ++quad_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fe.n_base_elements(); ++base)
+ for (unsigned int m=0; m<fe.element_multiplicity(base);
+ ++m, comp_start += fe.base_element(base).n_components() * do_tensor_product)
+ for (unsigned int local_index = 0;
+ local_index < fe.base_element(base).dofs_per_quad;
+ ++local_index, ++total_index)
+ {
+ // do everything alike for this type of object
+ const unsigned int index_in_base
+ = (fe.base_element(base).dofs_per_quad*quad_number +
+ local_index +
+ fe.base_element(base).first_quad_index);
+
+ const unsigned int face_index_in_base
+ = (fe.base_element(base).first_face_quad_index +
+ fe.base_element(base).dofs_per_quad * quad_number +
+ local_index);
+
+ face_system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), face_index_in_base);
+
+ if (fe.base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base
+ = fe.base_element(base).face_system_to_component_index(face_index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int face_index_in_comp
+ = fe.base_element(base).face_system_to_component_index(face_index_in_base).second;
+ face_system_to_component_table[total_index]
+ = std::make_pair (comp, face_index_in_comp);
+ }
+ else
+ face_system_to_component_table[total_index] = non_primitive_index;
+ }
+ }
+ Assert (total_index == fe.dofs_per_face, ExcInternalError());
+ Assert (total_index == face_system_to_component_table.size(),
+ ExcInternalError());
+ Assert (total_index == face_system_to_base_table.size(),
+ ExcInternalError());
+ }
}
projects / dealii.git / commitdiff