]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move some internal functions. 2106/head
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 23 Jan 2016 17:43:40 +0000 (11:43 -0600)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 23 Jan 2016 18:25:11 +0000 (12:25 -0600)
FESystem has a number of functions that are called in the member initializer list
of the constructors. These are all static functions, of course, and don't access
any member variables.

This patch moves them ouf of the class altogether. This has the advantage that
it reduces the size of the .h file, and that it reduces the number of functions
exported by the shared libraries. Rather, these functions may as well live
inside an anonymous namespace in the .cc file.

There are no functional changes, nor anything that would be visible to the user.

include/deal.II/fe/fe_system.h
source/fe/fe_system.cc

index bcb53a1e1088bab154014561948e3d3cb5fae9a2..cb00921b29db4124842c9f2eb33745dd31748662 100644 (file)
@@ -983,88 +983,6 @@ private:
    * of the FiniteElement class. Called from the constructor.
    */
   void initialize_quad_dof_index_permutation ();
-
-  /**
-   * Helper function used in the constructor: take a @p FiniteElementData
-   * object and return an object of the same type with the number of degrees
-   * of freedom per vertex, line, etc.  multiplied by @p n. Don't touch the
-   * number of functions for the transformation from unit to real cell.
-   */
-  static FiniteElementData<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);
-
-  /**
-   * Same as above but for any number of sub-elements.
-   */
-  static FiniteElementData<dim>
-  multiply_dof_numbers (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
-                        const std::vector<unsigned int>                       &multiplicities);
-
-
-
-  /**
-   * Helper function used in the constructor: takes a @p FiniteElement object
-   * and returns an boolean vector including the @p
-   * restriction_is_additive_flags of the mixed element consisting of @p N
-   * elements of the sub-element @p fe.
-   */
-  static std::vector<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 named flags for a list of finite elements with multiplicities
-   * given in the second argument. This function is called from all the above
-   * functions.
-   */
-  static std::vector<bool>
-  compute_restriction_is_additive_flags (
-    const std::vector<const FiniteElement<dim,spacedim>*> &fes,
-    const std::vector<unsigned int>              &multiplicities);
-
-
-  /**
-   * Compute the non-zero vector components of a composed finite element.
-   */
-  static 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 of a list of finite elements with
-   * multiplicities given in the second argument. This function is called from
-   * all the above functions.
-   */
-  static std::vector<ComponentMask>
-  compute_nonzero_components (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
-                              const std::vector<unsigned int>              &multiplicities);
-
   /**
    * This function is simply singled out of the constructors since there are
    * several of them. It sets up the index table for the system as well as @p
index 2d1265f22f5bb115e0fe5ce014aacface0faf793..aefbfa6849bbd548e511c4a4f150d6ebfbd94ead 100644 (file)
@@ -110,6 +110,491 @@ InternalData::get_fe_output_object (const unsigned int base_no) const
 template <int dim, int spacedim>
 const unsigned int FESystem<dim,spacedim>::invalid_face_number;
 
+namespace
+{
+
+  /**
+   * Take vectors of finite elements and multiplicities and multiply out
+   * how many degrees of freedom the composed element has per vertex,
+   * line, etc.
+   */
+  template <int dim, int spacedim>
+  FiniteElementData<dim>
+  multiply_dof_numbers (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
+                        const std::vector<unsigned int>                       &multiplicities)
+  {
+    Assert (fes.size() == multiplicities.size(), ExcDimensionMismatch (fes.size(), multiplicities.size()));
+
+    unsigned int multiplied_dofs_per_vertex = 0;
+    unsigned int multiplied_dofs_per_line = 0;
+    unsigned int multiplied_dofs_per_quad = 0;
+    unsigned int multiplied_dofs_per_hex = 0;
+
+    unsigned int multiplied_n_components = 0;
+
+    unsigned int degree = 0; // degree is the maximal degree of the components
+
+    for (unsigned int i=0; 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];
+
+          multiplied_n_components+=fes[i]->n_components() * multiplicities[i];
+
+          degree = std::max(degree, fes[i]->tensor_degree() );
+        }
+
+    // assume conformity of the first finite element and then take away
+    // bits as indicated by the base elements. if all multiplicities
+    // happen to be zero, then it doesn't matter what we set it to.
+    typename FiniteElementData<dim>::Conformity total_conformity
+      = typename FiniteElementData<dim>::Conformity();
+    {
+      unsigned int index = 0;
+      for (index=0; index<fes.size(); ++index)
+        if (multiplicities[index]>0)
+          {
+            total_conformity = fes[index]->conforming_space;
+            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);
+
+    return FiniteElementData<dim> (dpo,
+                                   multiplied_n_components,
+                                   degree,
+                                   total_conformity);
+  }
+
+  /**
+   * 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)
+  {
+    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);
+  }
+
+
+
+  /**
+   * Compute the named flags for a list of finite elements with multiplicities
+   * given in the second argument. This function is called from all the above
+   * functions.
+   */
+  template <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)
+  {
+    Assert (fes.size() == multiplicities.size(), ExcInternalError());
+
+    // first count the number of dofs and components that will emerge from the
+    // given FEs
+    unsigned int n_shape_functions = 0;
+    for (unsigned int i=0; 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_vertex;
+                 ++local_index, ++total_index)
+              {
+                const unsigned int index_in_base
+                  = (fes[base]->dofs_per_vertex*vertex_number +
+                     local_index);
+
+                Assert (index_in_base < fes[base]->dofs_per_cell,
+                        ExcInternalError());
+                retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
+              }
+      }
+
+    // 2. Lines
+    if (GeometryInfo<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);
+
+                  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);
+                }
+        }
+
+    // 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);
+
+                  Assert (index_in_base < fes[base]->dofs_per_cell,
+                          ExcInternalError());
+                  retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
+                }
+        }
+
+    Assert (total_index == n_shape_functions, ExcInternalError());
+
+    return retval;
+  }
+
+
+
+  /**
+   * Take a @p FiniteElement object
+   * and return an boolean vector including the @p
+   * restriction_is_additive_flags of the mixed element consisting of @p N
+   * elements of the sub-element @p fe.
+   */
+  template <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)
+  {
+    std::vector<const FiniteElement<dim,spacedim>*> fe_list;
+    std::vector<unsigned int>              multiplicities;
+
+    fe_list.push_back (fe1);
+    multiplicities.push_back (N1);
+
+    fe_list.push_back (fe2);
+    multiplicities.push_back (N2);
+
+    fe_list.push_back (fe3);
+    multiplicities.push_back (N3);
+
+    fe_list.push_back (fe4);
+    multiplicities.push_back (N4);
+
+    fe_list.push_back (fe5);
+    multiplicities.push_back (N5);
+    return compute_restriction_is_additive_flags (fe_list, multiplicities);
+  }
+
+
+
+  /**
+   * Compute the nonzero components of a list of finite elements with
+   * multiplicities given in the second argument.
+   */
+  template <int dim, int spacedim>
+  std::vector<ComponentMask>
+  compute_nonzero_components (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
+                              const std::vector<unsigned int>              &multiplicities)
+  {
+    Assert (fes.size() == multiplicities.size(), ExcInternalError());
+
+    // first count the number of dofs and components that will emerge from the
+    // given FEs
+    unsigned int n_shape_functions = 0;
+    for (unsigned int i=0; i<fes.size(); ++i)
+      if (multiplicities[i]>0) //needed because fe might be NULL
+        n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i];
+
+    unsigned int n_components = 0;
+    for (unsigned int i=0; i<fes.size(); ++i)
+      if (multiplicities[i]>0) //needed because fe might be NULL
+        n_components += fes[i]->n_components() * multiplicities[i];
+
+    // 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())
+            for (unsigned int local_index = 0;
+                 local_index < fes[base]->dofs_per_vertex;
+                 ++local_index, ++total_index)
+              {
+                const unsigned int index_in_base
+                  = (fes[base]->dofs_per_vertex*vertex_number +
+                     local_index);
+
+                Assert (comp_start+fes[base]->n_components() <=
+                        retval[total_index].size(),
+                        ExcInternalError());
+                for (unsigned int c=0; 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())
+              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];
+                    }
+                }
+        }
+
+    // 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())
+              for (unsigned int local_index = 0;
+                   local_index < fes[base]->dofs_per_quad;
+                   ++local_index, ++total_index)
+                {
+                  const unsigned int index_in_base
+                    = (fes[base]->dofs_per_quad*quad_number +
+                       local_index +
+                       fes[base]->first_quad_index);
+
+                  Assert (comp_start+fes[base]->n_components() <=
+                          retval[total_index].size(),
+                          ExcInternalError());
+                  for (unsigned int c=0; 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];
+                    }
+                }
+        }
+
+    // 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())
+              for (unsigned int local_index = 0;
+                   local_index < fes[base]->dofs_per_hex;
+                   ++local_index, ++total_index)
+                {
+                  const unsigned int index_in_base
+                    = (fes[base]->dofs_per_hex*hex_number +
+                       local_index +
+                       fes[base]->first_hex_index);
+
+                  Assert (comp_start+fes[base]->n_components() <=
+                          retval[total_index].size(),
+                          ExcInternalError());
+                  for (unsigned int c=0; 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];
+                    }
+                }
+        }
+
+    Assert (total_index == n_shape_functions, ExcInternalError());
+
+    // 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;
+  }
+
+
+  /**
+   * 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=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)
+  {
+    std::vector<const FiniteElement<dim,spacedim>*> fe_list;
+    std::vector<unsigned int>              multiplicities;
+
+    fe_list.push_back (fe1);
+    multiplicities.push_back (N1);
+
+    fe_list.push_back (fe2);
+    multiplicities.push_back (N2);
+
+    fe_list.push_back (fe3);
+    multiplicities.push_back (N3);
+
+    fe_list.push_back (fe4);
+    multiplicities.push_back (N4);
+
+    fe_list.push_back (fe5);
+    multiplicities.push_back (N5);
+
+    return compute_nonzero_components (fe_list, multiplicities);
+  }
+}
+
+
 
 template <int dim, int spacedim>
 FESystem<dim,spacedim>::FESystem (const FiniteElement<dim,spacedim> &fe,
@@ -2401,476 +2886,6 @@ compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const
 
 
 
-template <int dim, int spacedim>
-FiniteElementData<dim>
-FESystem<dim,spacedim>::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>
-FiniteElementData<dim>
-FESystem<dim,spacedim>::multiply_dof_numbers (
-  const std::vector<const FiniteElement<dim,spacedim>*>   &fes,
-  const std::vector<unsigned int>                   &multiplicities)
-{
-  Assert (fes.size() == multiplicities.size(), ExcDimensionMismatch (fes.size(), multiplicities.size()));
-
-  unsigned int multiplied_dofs_per_vertex = 0;
-  unsigned int multiplied_dofs_per_line = 0;
-  unsigned int multiplied_dofs_per_quad = 0;
-  unsigned int multiplied_dofs_per_hex = 0;
-
-  unsigned int multiplied_n_components = 0;
-
-  unsigned int degree = 0; // degree is the maximal degree of the components
-
-  for (unsigned int i=0; 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];
-
-        multiplied_n_components+=fes[i]->n_components() * multiplicities[i];
-
-        degree = std::max(degree, fes[i]->tensor_degree() );
-      }
-
-  // assume conformity of the first finite element and then take away
-  // bits as indicated by the base elements. if all multiplicities
-  // happen to be zero, then it doesn't matter what we set it to.
-  typename FiniteElementData<dim>::Conformity total_conformity
-    = typename FiniteElementData<dim>::Conformity();
-  {
-    unsigned int index = 0;
-    for (index=0; index<fes.size(); ++index)
-      if (multiplicities[index]>0)
-        {
-          total_conformity = fes[index]->conforming_space;
-          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);
-
-  return FiniteElementData<dim> (dpo,
-                                 multiplied_n_components,
-                                 degree,
-                                 total_conformity);
-}
-
-
-
-template <int dim, int spacedim>
-std::vector<bool>
-FESystem<dim,spacedim>::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);
-
-  fe_list.push_back (fe2);
-  multiplicities.push_back (N2);
-
-  fe_list.push_back (fe3);
-  multiplicities.push_back (N3);
-
-  fe_list.push_back (fe4);
-  multiplicities.push_back (N4);
-
-  fe_list.push_back (fe5);
-  multiplicities.push_back (N5);
-  return compute_restriction_is_additive_flags (fe_list, multiplicities);
-}
-
-
-template <int dim, int spacedim>
-std::vector<bool>
-FESystem<dim,spacedim>::
-compute_restriction_is_additive_flags (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
-                                       const std::vector<unsigned int>              &multiplicities)
-{
-  Assert (fes.size() == multiplicities.size(), ExcInternalError());
-
-  // first count the number of dofs and components that will emerge from the
-  // given FEs
-  unsigned int n_shape_functions = 0;
-  for (unsigned int i=0; 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_vertex;
-               ++local_index, ++total_index)
-            {
-              const unsigned int index_in_base
-                = (fes[base]->dofs_per_vertex*vertex_number +
-                   local_index);
-
-              Assert (index_in_base < fes[base]->dofs_per_cell,
-                      ExcInternalError());
-              retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
-            }
-    }
-
-  // 2. Lines
-  if (GeometryInfo<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);
-
-                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);
-              }
-      }
-
-  // 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);
-
-                Assert (index_in_base < fes[base]->dofs_per_cell,
-                        ExcInternalError());
-                retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
-              }
-      }
-
-  Assert (total_index == n_shape_functions, ExcInternalError());
-
-  return retval;
-}
-
-
-
-template <int dim, int spacedim>
-std::vector<ComponentMask>
-FESystem<dim,spacedim>::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;
-
-  fe_list.push_back (fe1);
-  multiplicities.push_back (N1);
-
-  fe_list.push_back (fe2);
-  multiplicities.push_back (N2);
-
-  fe_list.push_back (fe3);
-  multiplicities.push_back (N3);
-
-  fe_list.push_back (fe4);
-  multiplicities.push_back (N4);
-
-  fe_list.push_back (fe5);
-  multiplicities.push_back (N5);
-
-  return compute_nonzero_components (fe_list, multiplicities);
-}
-
-
-
-template <int dim, int spacedim>
-std::vector<ComponentMask>
-FESystem<dim,spacedim>::
-compute_nonzero_components (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
-                            const std::vector<unsigned int>              &multiplicities)
-{
-  Assert (fes.size() == multiplicities.size(), ExcInternalError());
-
-  // first count the number of dofs and components that will emerge from the
-  // given FEs
-  unsigned int n_shape_functions = 0;
-  for (unsigned int i=0; i<fes.size(); ++i)
-    if (multiplicities[i]>0) //needed because fe might be NULL
-      n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i];
-
-  unsigned int n_components = 0;
-  for (unsigned int i=0; i<fes.size(); ++i)
-    if (multiplicities[i]>0) //needed because fe might be NULL
-      n_components += fes[i]->n_components() * multiplicities[i];
-
-  // 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())
-          for (unsigned int local_index = 0;
-               local_index < fes[base]->dofs_per_vertex;
-               ++local_index, ++total_index)
-            {
-              const unsigned int index_in_base
-                = (fes[base]->dofs_per_vertex*vertex_number +
-                   local_index);
-
-              Assert (comp_start+fes[base]->n_components() <=
-                      retval[total_index].size(),
-                      ExcInternalError());
-              for (unsigned int c=0; c<fes[base]->n_components(); ++c)
-                {
-                  Assert (index_in_base < fes[base]->nonzero_components.size(),
-                          ExcInternalError());
-                  Assert (c < fes[base]->nonzero_components[index_in_base].size(),
-                          ExcInternalError());
-                  retval[total_index][comp_start+c]
-                    = fes[base]->nonzero_components[index_in_base][c];
-                };
-            }
-    }
-
-  // 2. Lines
-  if (GeometryInfo<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())
-            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 (index_in_base < fes[base]->nonzero_components.size(),
-                            ExcInternalError());
-                    Assert (c < fes[base]->nonzero_components[index_in_base].size(),
-                            ExcInternalError());
-                    retval[total_index][comp_start+c]
-                      = fes[base]->nonzero_components[index_in_base][c];
-                  };
-              }
-      }
-
-  // 3. Quads
-  if (GeometryInfo<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())
-            for (unsigned int local_index = 0;
-                 local_index < fes[base]->dofs_per_quad;
-                 ++local_index, ++total_index)
-              {
-                const unsigned int index_in_base
-                  = (fes[base]->dofs_per_quad*quad_number +
-                     local_index +
-                     fes[base]->first_quad_index);
-
-                Assert (comp_start+fes[base]->n_components() <=
-                        retval[total_index].size(),
-                        ExcInternalError());
-                for (unsigned int c=0; c<fes[base]->n_components(); ++c)
-                  {
-                    Assert (index_in_base < fes[base]->nonzero_components.size(),
-                            ExcInternalError());
-                    Assert (c < fes[base]->nonzero_components[index_in_base].size(),
-                            ExcInternalError());
-                    retval[total_index][comp_start+c]
-                      = fes[base]->nonzero_components[index_in_base][c];
-                  };
-              }
-      }
-
-  // 4. Hexes
-  if (GeometryInfo<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())
-            for (unsigned int local_index = 0;
-                 local_index < fes[base]->dofs_per_hex;
-                 ++local_index, ++total_index)
-              {
-                const unsigned int index_in_base
-                  = (fes[base]->dofs_per_hex*hex_number +
-                     local_index +
-                     fes[base]->first_hex_index);
-
-                Assert (comp_start+fes[base]->n_components() <=
-                        retval[total_index].size(),
-                        ExcInternalError());
-                for (unsigned int c=0; c<fes[base]->n_components(); ++c)
-                  {
-                    Assert (index_in_base < fes[base]->nonzero_components.size(),
-                            ExcInternalError());
-                    Assert (c < fes[base]->nonzero_components[index_in_base].size(),
-                            ExcInternalError());
-                    retval[total_index][comp_start+c]
-                      = fes[base]->nonzero_components[index_in_base][c];
-                  };
-              }
-      }
-
-  Assert (total_index == n_shape_functions, ExcInternalError());
-
-  // now copy the vector<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;
-}
-
-
-
-
 template <int dim, int spacedim>
 const FiniteElement<dim,spacedim> &
 FESystem<dim,spacedim>::base_element (const unsigned int index) const

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.