From 0668fd50b3c55b21510078a7eec494ac0c329845 Mon Sep 17 00:00:00 2001
From: Denis Davydov <davydden@gmail.com>
Date: Thu, 19 May 2016 12:36:08 +0200
Subject: [PATCH] move build cell tables to FETools

---
 include/deal.II/fe/fe_system.h |   5 -
 include/deal.II/fe/fe_tools.h  |  17 +++
 source/fe/fe_system.cc         | 198 ++-------------------------------
 source/fe/fe_tools.cc          | 182 ++++++++++++++++++++++++++++++
 source/fe/fe_tools.inst.in     |   8 ++
 5 files changed, 219 insertions(+), 191 deletions(-)

diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h
index 8c936a7e12..40323b7f48 100644
--- a/include/deal.II/fe/fe_system.h
+++ b/include/deal.II/fe/fe_system.h
@@ -982,11 +982,6 @@ private:
   void initialize (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
                    const std::vector<unsigned int> &multiplicities);
 
-  /**
-   * Used by @p initialize.
-   */
-  void build_cell_tables();
-
   /**
    * Used by @p initialize.
    */
diff --git a/include/deal.II/fe/fe_tools.h b/include/deal.II/fe/fe_tools.h
index d86bb1a2bc..f3e022d0d7 100644
--- a/include/deal.II/fe/fe_tools.h
+++ b/include/deal.II/fe/fe_tools.h
@@ -932,6 +932,23 @@ namespace FETools
                               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 <code>true</code>, the underlying finite element
+   * is assumed to be build using the tensor product rule. That is, the number of
+   * composite components is the sum of components in each finite element times
+   * multiplicity.
+   */
+  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);
+
   /**
    * Parse the name of a finite element and generate a finite element object
    * accordingly. The parser ignores space characters between words (things
diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc
index a3af7aa6b7..f7b6c34efb 100644
--- a/source/fe/fe_system.cc
+++ b/source/fe/fe_system.cc
@@ -1265,191 +1265,6 @@ compute_fill (const Mapping<dim,spacedim>                      &mapping,
 }
 
 
-
-template <int dim, int spacedim>
-void
-FESystem<dim,spacedim>::build_cell_tables()
-{
-  // If the system is not primitive, these have not been initialized by
-  // FiniteElement
-  this->system_to_component_table.resize(this->dofs_per_cell);
-  this->face_system_to_component_table.resize(this->dofs_per_face);
-
-  unsigned int total_index = 0;
-
-  for (unsigned int base=0; base < this->n_base_elements(); ++base)
-    for (unsigned int m = 0; m < this->element_multiplicity(base); ++m)
-      {
-        for (unsigned int k=0; k<base_element(base).n_components(); ++k)
-          this->component_to_base_table[total_index++]
-            = std::make_pair(std::make_pair(base,k), m);
-      }
-  Assert (total_index == this->component_to_base_table.size(),
-          ExcInternalError());
-
-  // 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<this->n_base_elements(); ++base)
-        for (unsigned int m=0; m<this->element_multiplicity(base);
-             ++m, comp_start+=base_element(base).n_components())
-          for (unsigned int local_index = 0;
-               local_index < base_element(base).dofs_per_vertex;
-               ++local_index, ++total_index)
-            {
-              const unsigned int index_in_base
-                = (base_element(base).dofs_per_vertex*vertex_number +
-                   local_index);
-
-              this->system_to_base_table[total_index]
-                = std::make_pair (std::make_pair(base, m), index_in_base);
-
-              if (base_element(base).is_primitive(index_in_base))
-                {
-                  const unsigned int comp_in_base
-                    = 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
-                    = base_element(base).system_to_component_index(index_in_base).second;
-                  this->system_to_component_table[total_index]
-                    = std::make_pair (comp, index_in_comp);
-                }
-              else
-                this->system_to_component_table[total_index] = non_primitive_index;
-            }
-    }
-
-  // 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<this->n_base_elements(); ++base)
-          for (unsigned int m=0; m<this->element_multiplicity(base);
-               ++m, comp_start+=base_element(base).n_components())
-            for (unsigned int local_index = 0;
-                 local_index < base_element(base).dofs_per_line;
-                 ++local_index, ++total_index)
-              {
-                const unsigned int index_in_base
-                  = (base_element(base).dofs_per_line*line_number +
-                     local_index +
-                     base_element(base).first_line_index);
-
-                this->system_to_base_table[total_index]
-                  = std::make_pair (std::make_pair(base,m), index_in_base);
-
-                if (base_element(base).is_primitive(index_in_base))
-                  {
-                    const unsigned int comp_in_base
-                      = 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
-                      = base_element(base).system_to_component_index(index_in_base).second;
-                    this->system_to_component_table[total_index]
-                      = std::make_pair (comp, index_in_comp);
-                  }
-                else
-                  this->system_to_component_table[total_index] = non_primitive_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<this->n_base_elements(); ++base)
-          for (unsigned int m=0; m<this->element_multiplicity(base);
-               ++m, comp_start += base_element(base).n_components())
-            for (unsigned int local_index = 0;
-                 local_index < base_element(base).dofs_per_quad;
-                 ++local_index, ++total_index)
-              {
-                const unsigned int index_in_base
-                  = (base_element(base).dofs_per_quad*quad_number +
-                     local_index +
-                     base_element(base).first_quad_index);
-
-                this->system_to_base_table[total_index]
-                  = std::make_pair (std::make_pair(base,m), index_in_base);
-
-                if (base_element(base).is_primitive(index_in_base))
-                  {
-                    const unsigned int comp_in_base
-                      = 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
-                      = base_element(base).system_to_component_index(index_in_base).second;
-                    this->system_to_component_table[total_index]
-                      = std::make_pair (comp, index_in_comp);
-                  }
-                else
-                  this->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<this->n_base_elements(); ++base)
-          for (unsigned int m=0; m<this->element_multiplicity(base);
-               ++m, comp_start+=base_element(base).n_components())
-            for (unsigned int local_index = 0;
-                 local_index < base_element(base).dofs_per_hex;
-                 ++local_index, ++total_index)
-              {
-                const unsigned int index_in_base
-                  = (base_element(base).dofs_per_hex*hex_number +
-                     local_index +
-                     base_element(base).first_hex_index);
-
-                this->system_to_base_table[total_index]
-                  = std::make_pair (std::make_pair(base,m), index_in_base);
-
-                if (base_element(base).is_primitive(index_in_base))
-                  {
-                    const unsigned int comp_in_base
-                      = 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
-                      = base_element(base).system_to_component_index(index_in_base).second;
-                    this->system_to_component_table[total_index]
-                      = std::make_pair (comp, index_in_comp);
-                  }
-                else
-                  this->system_to_component_table[total_index] = non_primitive_index;
-              }
-      }
-}
-
-
-
 template <int dim, int spacedim>
 void
 FESystem<dim,spacedim>::build_face_tables()
@@ -1831,7 +1646,18 @@ void FESystem<dim,spacedim>::initialize (const std::vector<const FiniteElement<d
 
   Assert(ind>0, ExcInternalError());
 
-  build_cell_tables();
+  {
+    // If the system is not primitive, these have not been initialized by
+    // FiniteElement
+    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);
+
+  }
   build_face_tables();
 
   // restriction and prolongation matrices are build on demand
diff --git a/source/fe/fe_tools.cc b/source/fe/fe_tools.cc
index 06edebbab9..c2ebb2a222 100644
--- a/source/fe/fe_tools.cc
+++ b/source/fe/fe_tools.cc
@@ -558,6 +558,188 @@ namespace FETools
     return compute_nonzero_components (fe_list, 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;
+
+    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());
+
+
+    // 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())
+            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);
+
+                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;
+              }
+      }
+
+    // 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())
+              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;
+                }
+        }
+
+    // 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())
+              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);
+
+                  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)
+        {
+          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())
+              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);
+
+                  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;
+                }
+        }
+  }
+
   // Not implemented in the general case.
   template <class FE>
   FiniteElement<FE::dimension, FE::space_dimension> *
diff --git a/source/fe/fe_tools.inst.in b/source/fe/fe_tools.inst.in
index 90535f80ef..6c7a389d51 100644
--- a/source/fe/fe_tools.inst.in
+++ b/source/fe/fe_tools.inst.in
@@ -74,6 +74,14 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
                                   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 compute_block_renumbering (
-- 
2.39.5