From 9b4a09104caea3841abe91392e1d6a0adfd444a6 Mon Sep 17 00:00:00 2001
From: Denis Davydov <davydden@gmail.com>
Date: Thu, 19 May 2016 13:39:08 +0200
Subject: [PATCH] same for fe_face_tables

---
 include/deal.II/fe/fe_system.h |   5 -
 include/deal.II/fe/fe_tools.h  |  14 +++
 source/fe/fe_system.cc         | 162 +--------------------------------
 source/fe/fe_tools.cc          | 159 ++++++++++++++++++++++++++++++++
 source/fe/fe_tools.inst.in     |   7 ++
 5 files changed, 184 insertions(+), 163 deletions(-)

diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h
index 40323b7f48..98064474d9 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_face_tables();
-
   /**
    * Used by @p initialize.
    */
diff --git a/include/deal.II/fe/fe_tools.h b/include/deal.II/fe/fe_tools.h
index f3e022d0d7..14fb9995e7 100644
--- a/include/deal.II/fe/fe_tools.h
+++ b/include/deal.II/fe/fe_tools.h
@@ -949,6 +949,20 @@ namespace FETools
                     const FiniteElement<dim,spacedim> &finite_element,
                     const bool do_tensor_product = true);
 
+  /**
+   * build...
+   * @param face_system_to_base_table
+   * @param face_system_to_component_table
+   * @param finite_element
+   * @param do_tensor_product
+   */
+  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
    * accordingly. The parser ignores space characters between words (things
diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc
index f7b6c34efb..9f1ff5a4dc 100644
--- a/source/fe/fe_system.cc
+++ b/source/fe/fe_system.cc
@@ -1265,163 +1265,6 @@ compute_fill (const Mapping<dim,spacedim>                      &mapping,
 }
 
 
-template <int dim, int spacedim>
-void
-FESystem<dim,spacedim>::build_face_tables()
-{
-  // 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<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)
-            {
-              // 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
-                = (base_element(base).dofs_per_vertex*vertex_number +
-                   local_index);
-
-              const unsigned int face_index_in_base
-                = (base_element(base).dofs_per_vertex*vertex_number +
-                   local_index);
-
-              this->face_system_to_base_table[total_index]
-                = std::make_pair (std::make_pair(base,m), face_index_in_base);
-
-              if (base_element(base).is_primitive(index_in_base))
-                {
-                  const unsigned int comp_in_base
-                    = 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
-                    = base_element(base).face_system_to_component_index(face_index_in_base).second;
-                  this->face_system_to_component_table[total_index]
-                    = std::make_pair (comp, face_index_in_comp);
-                }
-              else
-                this->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 < 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)
-              {
-                // do everything alike for this type of object
-                const unsigned int index_in_base
-                  = (base_element(base).dofs_per_line*line_number +
-                     local_index +
-                     base_element(base).first_line_index);
-
-                const unsigned int face_index_in_base
-                  = (base_element(base).first_face_line_index +
-                     base_element(base).dofs_per_line * line_number +
-                     local_index);
-
-                this->face_system_to_base_table[total_index]
-                  = std::make_pair (std::make_pair(base,m), face_index_in_base);
-
-                if (base_element(base).is_primitive(index_in_base))
-                  {
-                    const unsigned int comp_in_base
-                      = 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
-                      = base_element(base).face_system_to_component_index(face_index_in_base).second;
-                    this->face_system_to_component_table[total_index]
-                      = std::make_pair (comp, face_index_in_comp);
-                  }
-                else
-                  this->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<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)
-              {
-                // do everything alike for this type of object
-                const unsigned int index_in_base
-                  = (base_element(base).dofs_per_quad*quad_number +
-                     local_index +
-                     base_element(base).first_quad_index);
-
-                const unsigned int face_index_in_base
-                  = (base_element(base).first_face_quad_index +
-                     base_element(base).dofs_per_quad * quad_number +
-                     local_index);
-
-                this->face_system_to_base_table[total_index]
-                  = std::make_pair (std::make_pair(base,m), face_index_in_base);
-
-                if (base_element(base).is_primitive(index_in_base))
-                  {
-                    const unsigned int comp_in_base
-                      = 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
-                      = base_element(base).face_system_to_component_index(face_index_in_base).second;
-                    this->face_system_to_component_table[total_index]
-                      = std::make_pair (comp, face_index_in_comp);
-                  }
-                else
-                  this->face_system_to_component_table[total_index] = non_primitive_index;
-              }
-      }
-  Assert (total_index == this->dofs_per_face, ExcInternalError());
-  Assert (total_index == this->face_system_to_component_table.size(),
-          ExcInternalError());
-  Assert (total_index == this->face_system_to_base_table.size(),
-          ExcInternalError());
-}
-
-
-
 template <int dim, int spacedim>
 void FESystem<dim,spacedim>::build_interface_constraints ()
 {
@@ -1657,8 +1500,11 @@ void FESystem<dim,spacedim>::initialize (const std::vector<const FiniteElement<d
                                this->component_to_base_table,
                                *this);
 
+    FETools::build_face_tables(this->face_system_to_base_table,
+                               this->face_system_to_component_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 c2ebb2a222..8f29669e06 100644
--- a/source/fe/fe_tools.cc
+++ b/source/fe/fe_tools.cc
@@ -740,6 +740,165 @@ namespace FETools
         }
   }
 
+  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())
+            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))
+                  {
+                    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;
+              }
+      }
+
+    // 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())
+              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())
+              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());
+  }
+
+
   // 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 6c7a389d51..7a8aa495fe 100644
--- a/source/fe/fe_tools.inst.in
+++ b/source/fe/fe_tools.inst.in
@@ -82,6 +82,13 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
       	                 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 (
-- 
2.39.5