]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move functions into a separate namespace.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 25 May 2016 16:24:23 +0000 (11:24 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 25 May 2016 17:58:36 +0000 (12:58 -0500)
Also provide a higher level documentation.

doc/news/changes.h
include/deal.II/fe/fe_system.h
include/deal.II/fe/fe_tools.h
source/fe/fe_system.cc
source/fe/fe_tools.cc
source/fe/fe_tools.inst.in

index 87aa44cb90415d23e7d42a39affd7144d3f7d8d9..f4c40859da07fc9d1b54eea43a9adc985ead3efc 100644 (file)
@@ -269,11 +269,11 @@ inconvenience this causes.
  (Bruno Turcksin, 2016/05/22)
   </li>
 
- <li> New: There are now additional functions in the FETools namespace that build
+ <li> 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.
  <br>
- (Denis Davydov, 2016/05/20)
+ (Denis Davydov, Wolfgang Bangerth, 2016/05/20)
  </li>
 
  <li> New: Added PArpackSolver::reinit() when dealing with BlockVectors.
index 98064474d9cbae0d6f6e5ca396f98969d99b7196..0cd903b510ccf9d2002125cfcfa789e3d29b1c4d 100644 (file)
@@ -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
index 4c524b6cc40d6df2536f9a9afcffda623f0520e2..9dba8de61e652f6f54506bd8222552bae631508c 100644 (file)
@@ -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
+   * <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
index 9f1ff5a4dcfc192cc1a87e204a15b2be2891fb7b..f51088b4cf00a5fa6d218057b324cc862dc49321 100644 (file)
@@ -115,9 +115,9 @@ const unsigned int FESystem<dim,spacedim>::invalid_face_number;
 template <int dim, int spacedim>
 FESystem<dim,spacedim>::FESystem (const FiniteElement<dim,spacedim> &fe,
                                   const unsigned int n_elements) :
-  FiniteElement<dim,spacedim> (FETools::multiply_dof_numbers(&fe, n_elements),
-                               FETools::compute_restriction_is_additive_flags (&fe, n_elements),
-                               FETools::compute_nonzero_components(&fe, n_elements)),
+  FiniteElement<dim,spacedim> (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<const FiniteElement<dim,spacedim>*> fes;
@@ -134,10 +134,10 @@ FESystem<dim,spacedim>::FESystem (const FiniteElement<dim,spacedim> &fe1,
                                   const unsigned int        n1,
                                   const FiniteElement<dim,spacedim> &fe2,
                                   const unsigned int        n2) :
-  FiniteElement<dim,spacedim> (FETools::multiply_dof_numbers(&fe1, n1, &fe2, n2),
-                               FETools::compute_restriction_is_additive_flags (&fe1, n1,
+  FiniteElement<dim,spacedim> (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<dim,spacedim>::FESystem (const FiniteElement<dim,spacedim> &fe1,
                                   const unsigned int        n2,
                                   const FiniteElement<dim,spacedim> &fe3,
                                   const unsigned int        n3) :
-  FiniteElement<dim,spacedim> (FETools::multiply_dof_numbers(&fe1, n1,
-                                                             &fe2, n2,
-                                                             &fe3, n3),
-                               FETools::compute_restriction_is_additive_flags (&fe1, n1,
+  FiniteElement<dim,spacedim> (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<dim,spacedim>::FESystem (const FiniteElement<dim,spacedim> &fe1,
                                   const unsigned int        n3,
                                   const FiniteElement<dim,spacedim> &fe4,
                                   const unsigned int        n4) :
-  FiniteElement<dim,spacedim> (FETools::multiply_dof_numbers(&fe1, n1,
-                                                             &fe2, n2,
-                                                             &fe3, n3,
-                                                             &fe4, n4),
-                               FETools::compute_restriction_is_additive_flags (&fe1, n1,
+  FiniteElement<dim,spacedim> (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<dim,spacedim>::FESystem (const FiniteElement<dim,spacedim> &fe1,
                                   const unsigned int        n4,
                                   const FiniteElement<dim,spacedim> &fe5,
                                   const unsigned int        n5) :
-  FiniteElement<dim,spacedim> (FETools::multiply_dof_numbers(&fe1, n1,
-                                                             &fe2, n2,
-                                                             &fe3, n3,
-                                                             &fe4, n4,
-                                                             &fe5, n5),
-                               FETools::compute_restriction_is_additive_flags (&fe1, n1,
+  FiniteElement<dim,spacedim> (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<dim,spacedim>::FESystem (
   const std::vector<const FiniteElement<dim,spacedim>*>  &fes,
   const std::vector<unsigned int>                  &multiplicities)
   :
-  FiniteElement<dim,spacedim> (FETools::multiply_dof_numbers(fes, multiplicities),
-                               FETools::compute_restriction_is_additive_flags (fes, multiplicities),
-                               FETools::compute_nonzero_components(fes, multiplicities)),
+  FiniteElement<dim,spacedim> (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<dim,spacedim>::initialize (const std::vector<const FiniteElement<d
     this->system_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);
 
   }
 
index 7f87499c2b76f96382d261475e0c458a4784a27f..4c5922b10159c5c88f166abe72274d63d78a4dc3 100644 (file)
@@ -59,225 +59,159 @@ DEAL_II_NAMESPACE_OPEN
 
 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());
@@ -285,215 +219,188 @@ namespace FETools
                 }
         }
 
-    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(),
@@ -508,231 +415,213 @@ namespace FETools
                 }
         }
 
-    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))
                     {
@@ -749,139 +638,165 @@ namespace FETools
                     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]
@@ -902,11 +817,100 @@ namespace FETools
                     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());
+    }
   }
 
 
index 7a8aa495feefa3f9821640b2bddc9ef3e8b2e7dd..c45db192130b3ff87135ef2efdb7a4c27f0b9efd 100644 (file)
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 1998 - 2015 by the deal.II authors
+// Copyright (C) 1998 - 2016 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
 for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
   {
     namespace FETools
-      \{
+    \{
 #if deal_II_dimension <= deal_II_space_dimension
-          template
-          FiniteElementData<deal_II_dimension>
-       multiply_dof_numbers (const std::vector<const FiniteElement<deal_II_dimension,deal_II_space_dimension>*> &fes,
-                             const std::vector<unsigned int>                       &multiplicities,
-                             bool);
-                             
-       template
-       FiniteElementData<deal_II_dimension>
-       multiply_dof_numbers (const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe1,
-                             const unsigned int            N1,
-                             const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe2,
-                             const unsigned int            N2,
-                             const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe3,
-                             const unsigned int            N3,
-                             const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe4,
-                             const unsigned int            N4,
-                             const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe5,
-                             const unsigned int            N5);
-                             
-      template
-      std::vector<bool>
-      compute_restriction_is_additive_flags (const std::vector<const FiniteElement<deal_II_dimension,deal_II_space_dimension>*> &fes,
-                                             const std::vector<unsigned int>              &multiplicities);
-                                             
-      template
-      std::vector<bool>
-      compute_restriction_is_additive_flags (const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe1,
-                                             const unsigned int        N1,
-                                             const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe2,
-                                             const unsigned int        N2,
-                                             const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe3,
-                                             const unsigned int        N3,
-                                             const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe4,
-                                             const unsigned int        N4,
-                                             const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe5,
-                                             const unsigned int        N5);
-                                             
-      template
-      std::vector<ComponentMask>
-      compute_nonzero_components (const std::vector<const FiniteElement<deal_II_dimension,deal_II_space_dimension>*> &fes,
-                                  const std::vector<unsigned int>              &multiplicities,
+      namespace Compositing
+      \{
+         template
+         FiniteElementData<deal_II_dimension>
+         multiply_dof_numbers (const std::vector<const FiniteElement<deal_II_dimension,deal_II_space_dimension>*> &fes,
+                               const std::vector<unsigned int>                       &multiplicities,
+                               bool);
+
+         template
+         FiniteElementData<deal_II_dimension>
+         multiply_dof_numbers (const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe1,
+                               const unsigned int            N1,
+                               const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe2,
+                               const unsigned int            N2,
+                               const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe3,
+                               const unsigned int            N3,
+                               const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe4,
+                               const unsigned int            N4,
+                               const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe5,
+                               const unsigned int            N5);
+
+        template
+        std::vector<bool>
+        compute_restriction_is_additive_flags (const std::vector<const FiniteElement<deal_II_dimension,deal_II_space_dimension>*> &fes,
+                                               const std::vector<unsigned int>              &multiplicities);
+
+        template
+        std::vector<bool>
+        compute_restriction_is_additive_flags (const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe1,
+                                               const unsigned int        N1,
+                                               const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe2,
+                                               const unsigned int        N2,
+                                               const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe3,
+                                               const unsigned int        N3,
+                                               const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe4,
+                                               const unsigned int        N4,
+                                               const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe5,
+                                               const unsigned int        N5);
+
+        template
+        std::vector<ComponentMask>
+        compute_nonzero_components (const std::vector<const FiniteElement<deal_II_dimension,deal_II_space_dimension>*> &fes,
+                                    const std::vector<unsigned int>              &multiplicities,
                                   const bool do_tensor_product);
-                                  
-      template
-      std::vector<ComponentMask>
-      compute_nonzero_components (const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe1,
-                                  const unsigned int        N1,
-                                  const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe2,
-                                  const unsigned int        N2,
-                                  const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe3,
-                                  const unsigned int        N3,
-                                  const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe4,
-                                  const unsigned int        N4,
-                                  const FiniteElement<deal_II_dimension,deal_II_space_dimension> *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<deal_II_dimension,deal_II_space_dimension> &fe,
-                         const bool do_tensor_product);
-       
-       template                  
-       void
-       build_face_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &face_system_to_base_table,
-                         std::vector< std::pair< unsigned int, unsigned int > >                            &face_system_to_component_table,
-                         const FiniteElement<deal_II_dimension,deal_II_space_dimension> &fe,
-                         const bool do_tensor_product);
-                        
+
+        template
+        std::vector<ComponentMask>
+        compute_nonzero_components (const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe1,
+                                    const unsigned int        N1,
+                                    const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe2,
+                                    const unsigned int        N2,
+                                    const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe3,
+                                    const unsigned int        N3,
+                                    const FiniteElement<deal_II_dimension,deal_II_space_dimension> *fe4,
+                                    const unsigned int        N4,
+                                    const FiniteElement<deal_II_dimension,deal_II_space_dimension> *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<deal_II_dimension,deal_II_space_dimension> &fe,
+                           const bool do_tensor_product);
+
+         template                  
+         void
+         build_face_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &face_system_to_base_table,
+                           std::vector< std::pair< unsigned int, unsigned int > >                            &face_system_to_component_table,
+                           const FiniteElement<deal_II_dimension,deal_II_space_dimension> &fe,
+                           const bool do_tensor_product);
+
+      \}
+
       template
       void compute_block_renumbering (
         const FiniteElement<deal_II_dimension,deal_II_space_dimension> & ,
@@ -105,7 +109,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
       (const FiniteElement<deal_II_dimension,deal_II_space_dimension> &,
        std::vector<std::vector<FullMatrix<double> > > &, const bool, const double);
 #endif
-      \}
+   \}
   }
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.