From 52bef21ac22c6af07dec9b6f6a0fcda0e53fb982 Mon Sep 17 00:00:00 2001
From: David Wells <wellsd2@rpi.edu>
Date: Sun, 9 Jul 2017 14:04:17 -0400
Subject: [PATCH] Put internal FE functions in specific namespaces.

This is the finite element equivalent to 6310b2dee3e and is also necessary for
the (possible) unity build patch.
---
 source/fe/fe_bdm.cc            |  79 ++++++-----
 source/fe/fe_dgp_monomial.cc   | 149 ++++++++++----------
 source/fe/fe_dgq.cc            |  28 ++--
 source/fe/fe_enriched.cc       | 162 +++++++++++-----------
 source/fe/fe_face.cc           |  30 ++--
 source/fe/fe_nedelec.cc        |  27 ++--
 source/fe/fe_nothing.cc        |  17 ++-
 source/fe/fe_poly_tensor.cc    | 181 ++++++++++++------------
 source/fe/fe_q.cc              |  28 ++--
 source/fe/fe_q_base.cc         | 121 ++++++++--------
 source/fe/fe_q_bubbles.cc      | 244 +++++++++++++++++----------------
 source/fe/fe_q_hierarchical.cc |  36 +++--
 source/fe/fe_system.cc         |  25 ++--
 13 files changed, 597 insertions(+), 530 deletions(-)

diff --git a/source/fe/fe_bdm.cc b/source/fe/fe_bdm.cc
index f0e0e22dc8..48255014e6 100644
--- a/source/fe/fe_bdm.cc
+++ b/source/fe/fe_bdm.cc
@@ -249,48 +249,54 @@ FE_BDM<dim>::get_ria_vector (const unsigned int deg)
 }
 
 
-namespace
+namespace internal
 {
-  // This function sets up the values of the polynomials we want to
-  // take moments with in the quadrature points. In fact, we multiply
-  // thos by the weights, such that the sum of function values and
-  // test_values over quadrature points yields the interpolated degree
-  // of freedom.
-  template <int dim>
-  void
-  initialize_test_values (std::vector<std::vector<double> > &test_values,
-                          const Quadrature<dim> &quadrature,
-                          const unsigned int deg)
+  namespace FE_BDM
   {
-    PolynomialsP<dim> poly(deg);
-    std::vector<Tensor<1,dim> > dummy1;
-    std::vector<Tensor<2,dim> > dummy2;
-    std::vector<Tensor<3,dim> > dummy3;
-    std::vector<Tensor<4,dim> > dummy4;
+    namespace
+    {
+      // This function sets up the values of the polynomials we want to
+      // take moments with in the quadrature points. In fact, we multiply
+      // thos by the weights, such that the sum of function values and
+      // test_values over quadrature points yields the interpolated degree
+      // of freedom.
+      template <int dim>
+      void
+      initialize_test_values (std::vector<std::vector<double> > &test_values,
+                              const Quadrature<dim> &quadrature,
+                              const unsigned int deg)
+      {
+        PolynomialsP<dim> poly(deg);
+        std::vector<Tensor<1,dim> > dummy1;
+        std::vector<Tensor<2,dim> > dummy2;
+        std::vector<Tensor<3,dim> > dummy3;
+        std::vector<Tensor<4,dim> > dummy4;
 
-    test_values.resize(quadrature.size());
+        test_values.resize(quadrature.size());
 
-    for (unsigned int k=0; k<quadrature.size(); ++k)
-      {
-        test_values[k].resize(poly.n());
-        poly.compute(quadrature.point(k), test_values[k], dummy1, dummy2,
-                     dummy3, dummy4);
-        for (unsigned int i=0; i < poly.n(); ++i)
+        for (unsigned int k=0; k<quadrature.size(); ++k)
           {
-            test_values[k][i] *= quadrature.weight(k);
+            test_values[k].resize(poly.n());
+            poly.compute(quadrature.point(k), test_values[k], dummy1, dummy2,
+                         dummy3, dummy4);
+            for (unsigned int i=0; i < poly.n(); ++i)
+              {
+                test_values[k][i] *= quadrature.weight(k);
+              }
           }
       }
-  }
 
-  // This specialization only serves to avoid error messages. Nothing
-  // useful can be computed in dimension zero and thus the vector
-  // length stays zero.
-  template <>
-  void
-  initialize_test_values (std::vector<std::vector<double> > &,
-                          const Quadrature<0> &,
-                          const unsigned int)
-  {}
+      // This specialization only serves to avoid error messages. Nothing
+      // useful can be computed in dimension zero and thus the vector
+      // length stays zero.
+      template <>
+      void
+      initialize_test_values (std::vector<std::vector<double> > &,
+                              const Quadrature<0> &,
+                              const unsigned int)
+      {}
+    }
+  }
 }
 
 
@@ -333,7 +339,7 @@ FE_BDM<dim>::initialize_support_points (const unsigned int deg)
   // point values on faces in 2D. In 3D, this is impossible, since the
   // moments are only taken with respect to PolynomialsP.
   if (dim>2)
-    initialize_test_values(test_values_face, face_points, deg);
+    internal::FE_BDM::initialize_test_values(test_values_face, face_points, deg);
 
   if (deg<=1) return;
 
@@ -347,7 +353,7 @@ FE_BDM<dim>::initialize_support_points (const unsigned int deg)
   // the test functions in the
   // interior quadrature points
 
-  initialize_test_values(test_values_cell, cell_points, deg-2);
+  internal::FE_BDM::initialize_test_values(test_values_cell, cell_points, deg-2);
 }
 
 
@@ -356,4 +362,3 @@ FE_BDM<dim>::initialize_support_points (const unsigned int deg)
 #include "fe_bdm.inst"
 
 DEAL_II_NAMESPACE_CLOSE
-
diff --git a/source/fe/fe_dgp_monomial.cc b/source/fe/fe_dgp_monomial.cc
index 27316601e1..9a5f0686aa 100644
--- a/source/fe/fe_dgp_monomial.cc
+++ b/source/fe/fe_dgp_monomial.cc
@@ -23,88 +23,93 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-// namespace for some functions that are used in this file.
-namespace
+namespace internal
 {
-  // storage of hand-chosen support
-  // points
-  //
-  // For dim=2, dofs_per_cell of
-  // FE_DGPMonomial(k) is given by
-  // 0.5(k+1)(k+2), i.e.
-  //
-  // k    0  1  2  3  4  5  6  7
-  // dofs 1  3  6 10 15 21 28 36
-  //
-  // indirect access of unit points:
-  // the points for degree k are
-  // located at
-  //
-  // points[start_index[k]..start_index[k+1]-1]
-  const unsigned int start_index2d[6]= {0,1,4,10,20,35};
-  const double points2d[35][2]=
+  namespace FE_DGPMonomial
   {
-    {0,0},
-    {0,0},{1,0},{0,1},
-    {0,0},{1,0},{0,1},{1,1},{0.5,0},{0,0.5},
-    {0,0},{1,0},{0,1},{1,1},{1./3.,0},{2./3.,0},{0,1./3.},{0,2./3.},{0.5,1},{1,0.5},
-    {0,0},{1,0},{0,1},{1,1},{0.25,0},{0.5,0},{0.75,0},{0,0.25},{0,0.5},{0,0.75},{1./3.,1},{2./3.,1},{1,1./3.},{1,2./3.},{0.5,0.5}
-  };
-
-  // For dim=3, dofs_per_cell of
-  // FE_DGPMonomial(k) is given by
-  // 1./6.(k+1)(k+2)(k+3), i.e.
-  //
-  // k    0  1  2  3  4  5  6   7
-  // dofs 1  4 10 20 35 56 84 120
-  const unsigned int start_index3d[6]= {0,1,5,15/*,35*/};
-  const double points3d[35][3]=
-  {
-    {0,0,0},
-    {0,0,0},{1,0,0},{0,1,0},{0,0,1},
-    {0,0,0},{1,0,0},{0,1,0},{0,0,1},{0.5,0,0},{0,0.5,0},{0,0,0.5},{1,1,0},{1,0,1},{0,1,1}
-  };
+    namespace
+    {
+      // storage of hand-chosen support
+      // points
+      //
+      // For dim=2, dofs_per_cell of
+      // FE_DGPMonomial(k) is given by
+      // 0.5(k+1)(k+2), i.e.
+      //
+      // k    0  1  2  3  4  5  6  7
+      // dofs 1  3  6 10 15 21 28 36
+      //
+      // indirect access of unit points:
+      // the points for degree k are
+      // located at
+      //
+      // points[start_index[k]..start_index[k+1]-1]
+      const unsigned int start_index2d[6]= {0,1,4,10,20,35};
+      const double points2d[35][2]=
+      {
+        {0,0},
+        {0,0},{1,0},{0,1},
+        {0,0},{1,0},{0,1},{1,1},{0.5,0},{0,0.5},
+        {0,0},{1,0},{0,1},{1,1},{1./3.,0},{2./3.,0},{0,1./3.},{0,2./3.},{0.5,1},{1,0.5},
+        {0,0},{1,0},{0,1},{1,1},{0.25,0},{0.5,0},{0.75,0},{0,0.25},{0,0.5},{0,0.75},{1./3.,1},{2./3.,1},{1,1./3.},{1,2./3.},{0.5,0.5}
+      };
+
+      // For dim=3, dofs_per_cell of
+      // FE_DGPMonomial(k) is given by
+      // 1./6.(k+1)(k+2)(k+3), i.e.
+      //
+      // k    0  1  2  3  4  5  6   7
+      // dofs 1  4 10 20 35 56 84 120
+      const unsigned int start_index3d[6]= {0,1,5,15/*,35*/};
+      const double points3d[35][3]=
+      {
+        {0,0,0},
+        {0,0,0},{1,0,0},{0,1,0},{0,0,1},
+        {0,0,0},{1,0,0},{0,1,0},{0,0,1},{0.5,0,0},{0,0.5,0},{0,0,0.5},{1,1,0},{1,0,1},{0,1,1}
+      };
 
 
-  template <int dim>
-  void generate_unit_points (const unsigned int,
-                             std::vector<Point<dim> > &);
+      template <int dim>
+      void generate_unit_points (const unsigned int,
+                                 std::vector<Point<dim> > &);
 
-  template <>
-  void generate_unit_points (const unsigned int k,
-                             std::vector<Point<1> > &p)
-  {
-    Assert(p.size()==k+1, ExcDimensionMismatch(p.size(), k+1));
-    const double h = 1./k;
-    for (unsigned int i=0; i<p.size(); ++i)
-      p[i](0)=i*h;
-  }
+      template <>
+      void generate_unit_points (const unsigned int k,
+                                 std::vector<Point<1> > &p)
+      {
+        Assert(p.size()==k+1, ExcDimensionMismatch(p.size(), k+1));
+        const double h = 1./k;
+        for (unsigned int i=0; i<p.size(); ++i)
+          p[i](0)=i*h;
+      }
 
-  template <>
-  void generate_unit_points (const unsigned int k,
-                             std::vector<Point<2> > &p)
-  {
-    Assert(k<=4, ExcNotImplemented());
-    Assert(p.size()==start_index2d[k+1]-start_index2d[k], ExcInternalError());
-    for (unsigned int i=0; i<p.size(); ++i)
+      template <>
+      void generate_unit_points (const unsigned int k,
+                                 std::vector<Point<2> > &p)
       {
-        p[i](0)=points2d[start_index2d[k]+i][0];
-        p[i](1)=points2d[start_index2d[k]+i][1];
+        Assert(k<=4, ExcNotImplemented());
+        Assert(p.size()==start_index2d[k+1]-start_index2d[k], ExcInternalError());
+        for (unsigned int i=0; i<p.size(); ++i)
+          {
+            p[i](0)=points2d[start_index2d[k]+i][0];
+            p[i](1)=points2d[start_index2d[k]+i][1];
+          }
       }
-  }
 
-  template <>
-  void generate_unit_points (const unsigned int k,
-                             std::vector<Point<3> > &p)
-  {
-    Assert(k<=2, ExcNotImplemented());
-    Assert(p.size()==start_index3d[k+1]-start_index3d[k], ExcInternalError());
-    for (unsigned int i=0; i<p.size(); ++i)
+      template <>
+      void generate_unit_points (const unsigned int k,
+                                 std::vector<Point<3> > &p)
       {
-        p[i](0)=points3d[start_index3d[k]+i][0];
-        p[i](1)=points3d[start_index3d[k]+i][1];
-        p[i](2)=points3d[start_index3d[k]+i][2];
+        Assert(k<=2, ExcNotImplemented());
+        Assert(p.size()==start_index3d[k+1]-start_index3d[k], ExcInternalError());
+        for (unsigned int i=0; i<p.size(); ++i)
+          {
+            p[i](0)=points3d[start_index3d[k]+i][0];
+            p[i](1)=points3d[start_index3d[k]+i][1];
+            p[i](2)=points3d[start_index3d[k]+i][2];
+          }
       }
+    }
   }
 }
 
@@ -199,7 +204,7 @@ get_interpolation_matrix (const FiniteElement<dim> &source_fe,
   else
     {
       std::vector<Point<dim> > unit_points(this->dofs_per_cell);
-      generate_unit_points(this->degree, unit_points);
+      internal::FE_DGPMonomial::generate_unit_points(this->degree, unit_points);
 
       FullMatrix<double> source_fe_matrix(unit_points.size(), source_fe.dofs_per_cell);
       for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
diff --git a/source/fe/fe_dgq.cc b/source/fe/fe_dgq.cc
index 489ec3cab7..c61689604c 100644
--- a/source/fe/fe_dgq.cc
+++ b/source/fe/fe_dgq.cc
@@ -30,15 +30,21 @@
 DEAL_II_NAMESPACE_OPEN
 
 
-namespace
+namespace internal
 {
-  std::vector<Point<1> >
-  get_QGaussLobatto_points (const unsigned int degree)
+  namespace FE_DGQ
   {
-    if (degree > 0)
-      return QGaussLobatto<1>(degree+1).get_points();
-    else
-      return std::vector<Point<1> >(1, Point<1>(0.5));
+    namespace
+    {
+      std::vector<Point<1> >
+      get_QGaussLobatto_points (const unsigned int degree)
+      {
+        if (degree > 0)
+          return QGaussLobatto<1>(degree+1).get_points();
+        else
+          return std::vector<Point<1> >(1, Point<1>(0.5));
+      }
+    }
   }
 }
 
@@ -48,14 +54,15 @@ template <int dim, int spacedim>
 FE_DGQ<dim, spacedim>::FE_DGQ (const unsigned int degree)
   :
   FE_Poly<TensorProductPolynomials<dim>, dim, spacedim>
-  (TensorProductPolynomials<dim>(Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))),
+  (TensorProductPolynomials<dim>(Polynomials::generate_complete_Lagrange_basis
+                                 (internal::FE_DGQ::get_QGaussLobatto_points(degree))),
    FiniteElementData<dim>(get_dpo_vector(degree), 1, degree, FiniteElementData<dim>::L2),
    std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(degree),1, degree).dofs_per_cell, true),
    std::vector<ComponentMask>(FiniteElementData<dim>(get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector<bool>(1,true)))
 {
   // Compute support points, which are the tensor product of the Lagrange
   // interpolation points in the constructor.
-  Quadrature<dim> support_quadrature(get_QGaussLobatto_points(degree));
+  Quadrature<dim> support_quadrature(internal::FE_DGQ::get_QGaussLobatto_points(degree));
   Assert (support_quadrature.get_points().size() > 0,
           (typename FiniteElement<dim, spacedim>::ExcFEHasNoSupportPoints ()));
   this->unit_support_points = support_quadrature.get_points();
@@ -872,7 +879,8 @@ FE_DGQLegendre<dim,spacedim>::clone() const
 template <int dim, int spacedim>
 FE_DGQHermite<dim,spacedim>::FE_DGQHermite (const unsigned int degree)
   : FE_DGQ<dim,spacedim>(degree < 3 ?
-                         Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))
+                         Polynomials::generate_complete_Lagrange_basis
+                         (internal::FE_DGQ::get_QGaussLobatto_points(degree))
                          :
                          Polynomials::HermiteInterpolation::generate_complete_basis(degree))
 {}
diff --git a/source/fe/fe_enriched.cc b/source/fe/fe_enriched.cc
index b98f2e6e86..4463d3c4b2 100644
--- a/source/fe/fe_enriched.cc
+++ b/source/fe/fe_enriched.cc
@@ -22,93 +22,99 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-namespace
+namespace internal
 {
-  /**
-   * Auxiliary function to create multiplicity vector from input enrichment functions.
-   */
-  template <typename T>
-  std::vector<unsigned int>
-  build_multiplicities(const std::vector<std::vector<T > > &functions )
+  namespace FE_Enriched
   {
-    std::vector<unsigned int> multiplicities;
-    multiplicities.push_back(1); // the first one is non-enriched FE
-    for (unsigned int i = 0; i < functions.size(); i++)
-      multiplicities.push_back(functions[i].size());
-
-    return multiplicities;
-  }
-
-
-  /**
-   * Auxiliary function to build FiniteElement's vector
-   */
-  template <int dim, int spacedim>
-  std::vector< const FiniteElement< dim, spacedim > * >
-  build_fes(const FiniteElement<dim,spacedim> *fe_base,
-            const std::vector<const FiniteElement<dim,spacedim> * > &fe_enriched)
-  {
-    std::vector< const FiniteElement< dim, spacedim > * > fes;
-    fes.push_back(fe_base);
-    for (unsigned int i = 0; i < fe_enriched.size(); i++)
-      fes.push_back(fe_enriched[i]);
-
-    return fes;
-  }
+    namespace
+    {
+      /**
+       * Auxiliary function to create multiplicity vector from input enrichment functions.
+       */
+      template <typename T>
+      std::vector<unsigned int>
+      build_multiplicities(const std::vector<std::vector<T > > &functions )
+      {
+        std::vector<unsigned int> multiplicities;
+        multiplicities.push_back(1); // the first one is non-enriched FE
+        for (unsigned int i = 0; i < functions.size(); i++)
+          multiplicities.push_back(functions[i].size());
 
+        return multiplicities;
+      }
 
-  /**
-   * Auxiliary function which check consistency of the input parameters.
-   * Returns true if everything is ok.
-   */
-  template <int dim, int spacedim>
-  bool
-  consistency_check (const std::vector< const FiniteElement< dim, spacedim > * > &fes,
-                     const std::vector< unsigned int > &multiplicities,
-                     const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > &functions)
-  {
-    AssertThrow(fes.size() > 0,
-                ExcMessage("FEs size should be >=1"));
-    AssertThrow(fes.size() == multiplicities.size(),
-                ExcMessage("FEs and multiplicities should have the same size"));
 
-    AssertThrow (functions.size() == fes.size() - 1,
-                 ExcDimensionMismatch(functions.size(), fes.size()-1));
+      /**
+       * Auxiliary function to build FiniteElement's vector
+       */
+      template <int dim, int spacedim>
+      std::vector< const FiniteElement< dim, spacedim > * >
+      build_fes(const FiniteElement<dim,spacedim> *fe_base,
+                const std::vector<const FiniteElement<dim,spacedim> * > &fe_enriched)
+      {
+        std::vector< const FiniteElement< dim, spacedim > * > fes;
+        fes.push_back(fe_base);
+        for (unsigned int i = 0; i < fe_enriched.size(); i++)
+          fes.push_back(fe_enriched[i]);
 
-    AssertThrow(multiplicities[0] == 1,
-                ExcMessage("First multiplicity should be 1"));
+        return fes;
+      }
 
-    const unsigned int n_comp_base = fes[0]->n_components();
 
-    // start from fe=1 as 0th is always non-enriched FE.
-    for (unsigned int fe=1; fe < fes.size(); fe++)
+      /**
+       * Auxiliary function which check consistency of the input parameters.
+       * Returns true if everything is ok.
+       */
+      template <int dim, int spacedim>
+      bool
+      consistency_check (const std::vector< const FiniteElement< dim, spacedim > * > &fes,
+                         const std::vector< unsigned int > &multiplicities,
+                         const std::vector<std::vector<std::function<const Function<spacedim> *(const typename dealii::Triangulation<dim, spacedim>::cell_iterator &) > > > &functions)
       {
-        const FE_Nothing<dim> *fe_nothing = dynamic_cast<const FE_Nothing<dim>*>(fes[fe]);
-        if (fe_nothing)
-          AssertThrow (fe_nothing->is_dominating(),
-                       ExcMessage("Only dominating FE_Nothing can be used in FE_Enriched"));
-
-        AssertThrow (fes[fe]->n_components() == n_comp_base,
-                     ExcMessage("All elements must have the same number of components"));
+        AssertThrow(fes.size() > 0,
+                    ExcMessage("FEs size should be >=1"));
+        AssertThrow(fes.size() == multiplicities.size(),
+                    ExcMessage("FEs and multiplicities should have the same size"));
+
+        AssertThrow (functions.size() == fes.size() - 1,
+                     ExcDimensionMismatch(functions.size(), fes.size()-1));
+
+        AssertThrow(multiplicities[0] == 1,
+                    ExcMessage("First multiplicity should be 1"));
+
+        const unsigned int n_comp_base = fes[0]->n_components();
+
+        // start from fe=1 as 0th is always non-enriched FE.
+        for (unsigned int fe=1; fe < fes.size(); fe++)
+          {
+            const FE_Nothing<dim> *fe_nothing = dynamic_cast<const FE_Nothing<dim>*>(fes[fe]);
+            if (fe_nothing)
+              AssertThrow (fe_nothing->is_dominating(),
+                           ExcMessage("Only dominating FE_Nothing can be used in FE_Enriched"));
+
+            AssertThrow (fes[fe]->n_components() == n_comp_base,
+                         ExcMessage("All elements must have the same number of components"));
+          }
+        return true;
       }
-    return true;
-  }
 
 
-  /**
-   * Auxiliary function which determines whether the FiniteElement will be enriched.
-   */
-  template <int dim, int spacedim>
-  bool
-  check_if_enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes)
-  {
-    // start from fe=1 as 0th is always non-enriched FE.
-    for (unsigned int fe=1; fe < fes.size(); fe++)
-      if (dynamic_cast<const FE_Nothing<dim>*>(fes[fe]) == nullptr)
-        // this is not FE_Nothing => there will be enrichment
-        return true;
+      /**
+       * Auxiliary function which determines whether the FiniteElement will be enriched.
+       */
+      template <int dim, int spacedim>
+      bool
+      check_if_enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes)
+      {
+        // start from fe=1 as 0th is always non-enriched FE.
+        for (unsigned int fe=1; fe < fes.size(); fe++)
+          if (dynamic_cast<const FE_Nothing<dim>*>(fes[fe]) == nullptr)
+            // this is not FE_Nothing => there will be enrichment
+            return true;
 
-    return false;
+        return false;
+      }
+    }
   }
 }
 
@@ -149,8 +155,8 @@ FE_Enriched<dim,spacedim>::FE_Enriched (const FiniteElement<dim,spacedim> *fe_ba
                                         const std::vector<const FiniteElement<dim,spacedim> * > &fe_enriched,
                                         const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > &functions)
   :
-  FE_Enriched<dim,spacedim> (build_fes(fe_base,fe_enriched),
-                             build_multiplicities(functions),
+  FE_Enriched<dim,spacedim> (internal::FE_Enriched::build_fes(fe_base,fe_enriched),
+                             internal::FE_Enriched::build_multiplicities(functions),
                              functions)
 {}
 
@@ -164,11 +170,11 @@ FE_Enriched<dim,spacedim>::FE_Enriched (const std::vector< const FiniteElement<
                                FETools::Compositing::compute_restriction_is_additive_flags(fes,multiplicities),
                                FETools::Compositing::compute_nonzero_components(fes,multiplicities,false)),
   enrichments(functions),
-  is_enriched(check_if_enriched(fes)),
+  is_enriched(internal::FE_Enriched::check_if_enriched(fes)),
   fe_system(fes,multiplicities)
 {
   // descriptive error are thrown within the function.
-  Assert(consistency_check(fes,multiplicities,functions),
+  Assert(internal::FE_Enriched::consistency_check(fes,multiplicities,functions),
          ExcInternalError());
 
   initialize(fes, multiplicities);
diff --git a/source/fe/fe_face.cc b/source/fe/fe_face.cc
index dad2e319f2..448d2b3b6c 100644
--- a/source/fe/fe_face.cc
+++ b/source/fe/fe_face.cc
@@ -27,25 +27,31 @@
 DEAL_II_NAMESPACE_OPEN
 
 
-namespace
+namespace internal
 {
-  std::vector<Point<1> >
-  get_QGaussLobatto_points (const unsigned int degree)
+  namespace FE_FaceQ
   {
-    if (degree > 0)
-      return QGaussLobatto<1>(degree+1).get_points();
-    else
-      return std::vector<Point<1> >(1, Point<1>(0.5));
+    namespace
+    {
+      std::vector<Point<1> >
+      get_QGaussLobatto_points (const unsigned int degree)
+      {
+        if (degree > 0)
+          return QGaussLobatto<1>(degree+1).get_points();
+        else
+          return std::vector<Point<1> >(1, Point<1>(0.5));
+      }
+    }
   }
 }
 
 template <int dim, int spacedim>
 FE_FaceQ<dim,spacedim>::FE_FaceQ (const unsigned int degree)
   :
-  FE_PolyFace<TensorProductPolynomials<dim-1>, dim, spacedim> (
-    TensorProductPolynomials<dim-1>(Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))),
-    FiniteElementData<dim>(get_dpo_vector(degree), 1, degree, FiniteElementData<dim>::L2),
-    std::vector<bool>(1,true))
+  FE_PolyFace<TensorProductPolynomials<dim-1>, dim, spacedim>
+  (TensorProductPolynomials<dim-1>(Polynomials::generate_complete_Lagrange_basis(internal::FE_FaceQ::get_QGaussLobatto_points(degree))),
+   FiniteElementData<dim>(get_dpo_vector(degree), 1, degree, FiniteElementData<dim>::L2),
+   std::vector<bool>(1,true))
 {
   // initialize unit face support points
   const unsigned int codim = dim-1;
@@ -56,7 +62,7 @@ FE_FaceQ<dim,spacedim>::FE_FaceQ (const unsigned int degree)
       this->unit_face_support_points[0][d] = 0.5;
   else
     {
-      std::vector<Point<1> > points = get_QGaussLobatto_points(degree);
+      std::vector<Point<1> > points = internal::FE_FaceQ::get_QGaussLobatto_points(degree);
 
       unsigned int k=0;
       for (unsigned int iz=0; iz <= ((codim>2) ? this->degree : 0) ; ++iz)
diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc
index d9c557968a..d3c5d7eea6 100644
--- a/source/fe/fe_nedelec.cc
+++ b/source/fe/fe_nedelec.cc
@@ -44,17 +44,20 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace internal
 {
-  namespace
+  namespace FE_Nedelec
   {
-    double
-    get_embedding_computation_tolerance(const unsigned int p)
+    namespace
     {
-      // This heuristic was computed by monitoring the worst residual
-      // resulting from the least squares computation when computing
-      // the face embedding matrices in the FE_Nedelec constructor.
-      // The residual growth is exponential, but is bounded by this
-      // function up to degree 12.
-      return 1.e-15*std::exp(std::pow(p,1.075));
+      double
+      get_embedding_computation_tolerance(const unsigned int p)
+      {
+        // This heuristic was computed by monitoring the worst residual
+        // resulting from the least squares computation when computing
+        // the face embedding matrices in the FE_Nedelec constructor.
+        // The residual growth is exponential, but is bounded by this
+        // function up to degree 12.
+        return 1.e-15*std::exp(std::pow(p,1.075));
+      }
     }
   }
 }
@@ -107,7 +110,7 @@ FE_Nedelec<dim>::FE_Nedelec (const unsigned int order)
 
   FETools::compute_face_embedding_matrices<dim,double>
   (*this, face_embeddings, 0, 0,
-   internal::get_embedding_computation_tolerance(order));
+   internal::FE_Nedelec::get_embedding_computation_tolerance(order));
 
   switch (dim)
     {
@@ -3022,7 +3025,7 @@ FE_Nedelec<dim>
       this_nonconst.reinit_restriction_and_prolongation_matrices ();
       // Fill prolongation matrices with embedding operators
       FETools::compute_embedding_matrices (this_nonconst, this_nonconst.prolongation, true,
-                                           internal::get_embedding_computation_tolerance(this->degree));
+                                           internal::FE_Nedelec::get_embedding_computation_tolerance(this->degree));
 #ifdef DEBUG_NEDELEC
       deallog << "Restriction" << std::endl;
 #endif
@@ -3073,7 +3076,7 @@ FE_Nedelec<dim>
       this_nonconst.reinit_restriction_and_prolongation_matrices ();
       // Fill prolongation matrices with embedding operators
       FETools::compute_embedding_matrices (this_nonconst, this_nonconst.prolongation, true,
-                                           internal::get_embedding_computation_tolerance(this->degree));
+                                           internal::FE_Nedelec::get_embedding_computation_tolerance(this->degree));
 #ifdef DEBUG_NEDELEC
       deallog << "Restriction" << std::endl;
 #endif
diff --git a/source/fe/fe_nothing.cc b/source/fe/fe_nothing.cc
index 71216f1a05..ba059c1a62 100644
--- a/source/fe/fe_nothing.cc
+++ b/source/fe/fe_nothing.cc
@@ -19,10 +19,16 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-namespace
+namespace internal
 {
-  const char *
-  zero_dof_message = "This element has no shape functions.";
+  namespace FE_Nothing
+  {
+    namespace
+    {
+      const char *
+      zero_dof_message = "This element has no shape functions.";
+    }
+  }
 }
 
 
@@ -85,8 +91,8 @@ double
 FE_Nothing<dim,spacedim>::shape_value (const unsigned int /*i*/,
                                        const Point<dim> & /*p*/) const
 {
-  (void)zero_dof_message;
-  Assert(false,ExcMessage(zero_dof_message));
+  (void)internal::FE_Nothing::zero_dof_message;
+  Assert(false,ExcMessage(internal::FE_Nothing::zero_dof_message));
   return 0;
 }
 
@@ -306,4 +312,3 @@ get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> & /*source_f
 
 
 DEAL_II_NAMESPACE_CLOSE
-
diff --git a/source/fe/fe_poly_tensor.cc b/source/fe/fe_poly_tensor.cc
index a88e6155d7..d3d9626b44 100644
--- a/source/fe/fe_poly_tensor.cc
+++ b/source/fe/fe_poly_tensor.cc
@@ -23,109 +23,116 @@
 #include <deal.II/fe/fe_poly_tensor.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/mapping_cartesian.h>
+#include <deal.II/grid/tria.h>
 
 DEAL_II_NAMESPACE_OPEN
 
-namespace
+namespace internal
 {
-//---------------------------------------------------------------------------
-// Utility method, which is used to determine the change of sign for
-// the DoFs on the faces of the given cell.
-//---------------------------------------------------------------------------
-
-  /**
-   * On noncartesian grids, the sign of the DoFs associated with the faces of
-   * the elements has to be changed in some cases.  This procedure implements an
-   * algorithm, which determines the DoFs, which need this sign change for a
-   * given cell.
-   */
-  void
-  get_face_sign_change_rt (const Triangulation<1>::cell_iterator &,
-                           const unsigned int,
-                           std::vector<double>                   &face_sign)
+  namespace FE_PolyTensor
   {
-    // nothing to do in 1d
-    std::fill (face_sign.begin (), face_sign.end (), 1.0);
-  }
-
-
+    namespace
+    {
+      //---------------------------------------------------------------------------
+      // Utility method, which is used to determine the change of sign for
+      // the DoFs on the faces of the given cell.
+      //---------------------------------------------------------------------------
+
+      /**
+       * On noncartesian grids, the sign of the DoFs associated with the faces of
+       * the elements has to be changed in some cases.  This procedure implements an
+       * algorithm, which determines the DoFs, which need this sign change for a
+       * given cell.
+       */
+      void
+      get_face_sign_change_rt (const dealii::Triangulation<1>::cell_iterator &,
+                               const unsigned int,
+                               std::vector<double>                           &face_sign)
+      {
+        // nothing to do in 1d
+        std::fill (face_sign.begin (), face_sign.end (), 1.0);
+      }
 
-  void
-  get_face_sign_change_rt (const Triangulation<2>::cell_iterator &cell,
-                           const unsigned int                     dofs_per_face,
-                           std::vector<double>                   &face_sign)
-  {
-    const unsigned int dim = 2;
-    const unsigned int spacedim = 2;
 
-    // Default is no sign
-    // change. I.e. multiply by one.
-    std::fill (face_sign.begin (), face_sign.end (), 1.0);
 
-    for (unsigned int f = GeometryInfo<dim>::faces_per_cell / 2;
-         f < GeometryInfo<dim>::faces_per_cell; ++f)
+      void
+      get_face_sign_change_rt (const dealii::Triangulation<2>::cell_iterator &cell,
+                               const unsigned int                             dofs_per_face,
+                               std::vector<double>                           &face_sign)
       {
-        Triangulation<dim,spacedim>::face_iterator face = cell->face (f);
-        if (!face->at_boundary ())
-          {
-            const unsigned int nn = cell->neighbor_face_no(f);
+        const unsigned int dim = 2;
+        const unsigned int spacedim = 2;
 
-            if (nn < GeometryInfo<dim>::faces_per_cell / 2)
-              for (unsigned int j = 0; j < dofs_per_face; ++j)
-                {
-                  Assert (f * dofs_per_face + j < face_sign.size(),
-                          ExcInternalError());
+        // Default is no sign
+        // change. I.e. multiply by one.
+        std::fill (face_sign.begin (), face_sign.end (), 1.0);
 
-//TODO: This is probably only going to work for those elements for which all dofs are face dofs
-                  face_sign[f * dofs_per_face + j] = -1.0;
-                }
+        for (unsigned int f = GeometryInfo<dim>::faces_per_cell / 2;
+             f < GeometryInfo<dim>::faces_per_cell; ++f)
+          {
+            dealii::Triangulation<dim,spacedim>::face_iterator face = cell->face (f);
+            if (!face->at_boundary ())
+              {
+                const unsigned int nn = cell->neighbor_face_no(f);
+
+                if (nn < GeometryInfo<dim>::faces_per_cell / 2)
+                  for (unsigned int j = 0; j < dofs_per_face; ++j)
+                    {
+                      Assert (f * dofs_per_face + j < face_sign.size(),
+                              ExcInternalError());
+
+                      //TODO: This is probably only going to work for those elements for which all dofs are face dofs
+                      face_sign[f * dofs_per_face + j] = -1.0;
+                    }
+              }
           }
       }
-  }
 
 
 
-  void
-  get_face_sign_change_rt (const Triangulation<3>::cell_iterator &/*cell*/,
-                           const unsigned int                     /*dofs_per_face*/,
-                           std::vector<double>                   &face_sign)
-  {
-    std::fill (face_sign.begin (), face_sign.end (), 1.0);
-//TODO: think about what it would take here
-  }
+      void
+      get_face_sign_change_rt (const dealii::Triangulation<3>::cell_iterator &/*cell*/,
+                               const unsigned int                             /*dofs_per_face*/,
+                               std::vector<double>                           &face_sign)
+      {
+        std::fill (face_sign.begin (), face_sign.end (), 1.0);
+        //TODO: think about what it would take here
+      }
 
-  void
-  get_face_sign_change_nedelec (const Triangulation<1>::cell_iterator &/*cell*/,
-                                const unsigned int                     /*dofs_per_face*/,
-                                std::vector<double>                   &face_sign)
-  {
-    // nothing to do in 1d
-    std::fill (face_sign.begin (), face_sign.end (), 1.0);
-  }
+      void
+      get_face_sign_change_nedelec (const dealii::Triangulation<1>::cell_iterator &/*cell*/,
+                                    const unsigned int                             /*dofs_per_face*/,
+                                    std::vector<double>                           &face_sign)
+      {
+        // nothing to do in 1d
+        std::fill (face_sign.begin (), face_sign.end (), 1.0);
+      }
 
 
 
-  void
-  get_face_sign_change_nedelec (const Triangulation<2>::cell_iterator &/*cell*/,
-                                const unsigned int                     /*dofs_per_face*/,
-                                std::vector<double>                   &face_sign)
-  {
-    std::fill (face_sign.begin (), face_sign.end (), 1.0);
-//TODO: think about what it would take here
-  }
+      void
+      get_face_sign_change_nedelec (const dealii::Triangulation<2>::cell_iterator &/*cell*/,
+                                    const unsigned int                             /*dofs_per_face*/,
+                                    std::vector<double>                           &face_sign)
+      {
+        std::fill (face_sign.begin (), face_sign.end (), 1.0);
+        //TODO: think about what it would take here
+      }
 
 
-  void
-  get_face_sign_change_nedelec (const Triangulation<3>::cell_iterator &cell,
-                                const unsigned int                     /*dofs_per_face*/,
-                                std::vector<double>                   &face_sign)
-  {
-    const unsigned int dim = 3;
-    std::fill (face_sign.begin (), face_sign.end (), 1.0);
-//TODO: This is probably only going to work for those elements for which all dofs are face dofs
-    for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
-      if (!(cell->line_orientation (l)))
-        face_sign[l] = -1.0;
+      void
+      get_face_sign_change_nedelec (const dealii::Triangulation<3>::cell_iterator &cell,
+                                    const unsigned int                             /*dofs_per_face*/,
+                                    std::vector<double>                           &face_sign)
+      {
+        const unsigned int dim = 3;
+        std::fill (face_sign.begin (), face_sign.end (), 1.0);
+        //TODO: This is probably only going to work for those elements for which all dofs are face dofs
+        for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+          if (!(cell->line_orientation (l)))
+            face_sign[l] = -1.0;
+      }
+    }
   }
 }
 
@@ -335,9 +342,9 @@ fill_fe_values
   std::fill( fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0 );
 
   if (mapping_type == mapping_raviart_thomas)
-    get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change);
+    internal::FE_PolyTensor::get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change);
   else if (mapping_type == mapping_nedelec)
-    get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change);
+    internal::FE_PolyTensor::get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change);
 
 
   for (unsigned int i=0; i<this->dofs_per_cell; ++i)
@@ -816,10 +823,10 @@ fill_fe_face_values
   std::fill( fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0 );
 
   if (mapping_type == mapping_raviart_thomas)
-    get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change);
+    internal::FE_PolyTensor::get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change);
 
   else if (mapping_type == mapping_nedelec)
-    get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change);
+    internal::FE_PolyTensor::get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change);
 
   for (unsigned int i=0; i<this->dofs_per_cell; ++i)
     {
@@ -1295,10 +1302,10 @@ fill_fe_subface_values
   std::fill( fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0 );
 
   if (mapping_type == mapping_raviart_thomas)
-    get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change);
+    internal::FE_PolyTensor::get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change);
 
   else if (mapping_type == mapping_nedelec)
-    get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change);
+    internal::FE_PolyTensor::get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change);
 
   for (unsigned int i=0; i<this->dofs_per_cell; ++i)
     {
diff --git a/source/fe/fe_q.cc b/source/fe/fe_q.cc
index 339205589a..a6941c5f76 100644
--- a/source/fe/fe_q.cc
+++ b/source/fe/fe_q.cc
@@ -25,19 +25,25 @@
 DEAL_II_NAMESPACE_OPEN
 
 
-namespace
+namespace internal
 {
-  std::vector<Point<1> >
-  get_QGaussLobatto_points (const unsigned int degree)
+  namespace FE_Q
   {
-    if (degree > 0)
-      return QGaussLobatto<1>(degree+1).get_points();
-    else
+    namespace
+    {
+      std::vector<Point<1> >
+      get_QGaussLobatto_points (const unsigned int degree)
       {
-        typedef FE_Q_Base<TensorProductPolynomials<1>, 1, 1> FEQ;
-        AssertThrow(false, FEQ::ExcFEQCannotHaveDegree0());
+        if (degree > 0)
+          return QGaussLobatto<1>(degree+1).get_points();
+        else
+          {
+            typedef FE_Q_Base<TensorProductPolynomials<1>, 1, 1> FEQ;
+            AssertThrow(false, FEQ::ExcFEQCannotHaveDegree0());
+          }
+        return std::vector<Point<1> >();
       }
-    return std::vector<Point<1> >();
+    }
   }
 }
 
@@ -47,13 +53,13 @@ template <int dim, int spacedim>
 FE_Q<dim,spacedim>::FE_Q (const unsigned int degree)
   :
   FE_Q_Base<TensorProductPolynomials<dim>, dim, spacedim>
-  (TensorProductPolynomials<dim>(Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))),
+  (TensorProductPolynomials<dim>(Polynomials::generate_complete_Lagrange_basis(internal::FE_Q::get_QGaussLobatto_points(degree))),
    FiniteElementData<dim>(this->get_dpo_vector(degree),
                           1, degree,
                           FiniteElementData<dim>::H1),
    std::vector<bool> (1, false))
 {
-  this->initialize(get_QGaussLobatto_points(degree));
+  this->initialize(internal::FE_Q::get_QGaussLobatto_points(degree));
 }
 
 
diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc
index b958a14177..ad3cc46af1 100644
--- a/source/fe/fe_q_base.cc
+++ b/source/fe/fe_q_base.cc
@@ -35,66 +35,69 @@
 DEAL_II_NAMESPACE_OPEN
 
 
-namespace FE_Q_Helper
+namespace internal
 {
-  namespace
+  namespace FE_Q_Base
   {
-    // get the renumbering for faces
-    template <int dim>
-    inline
-    std::vector<unsigned int>
-    face_lexicographic_to_hierarchic_numbering (const unsigned int degree)
+    namespace
     {
-      std::vector<unsigned int> dpo(dim, 1U);
-      for (unsigned int i=1; i<dpo.size(); ++i)
-        dpo[i]=dpo[i-1]*(degree-1);
-      const dealii::FiniteElementData<dim-1> face_data(dpo,1,degree);
-      std::vector<unsigned int> face_renumber (face_data.dofs_per_cell);
-      FETools::lexicographic_to_hierarchic_numbering (face_data, face_renumber);
-      return face_renumber;
-    }
+      // get the renumbering for faces
+      template <int dim>
+      inline
+      std::vector<unsigned int>
+      face_lexicographic_to_hierarchic_numbering (const unsigned int degree)
+      {
+        std::vector<unsigned int> dpo(dim, 1U);
+        for (unsigned int i=1; i<dpo.size(); ++i)
+          dpo[i]=dpo[i-1]*(degree-1);
+        const dealii::FiniteElementData<dim-1> face_data(dpo,1,degree);
+        std::vector<unsigned int> face_renumber (face_data.dofs_per_cell);
+        FETools::lexicographic_to_hierarchic_numbering (face_data, face_renumber);
+        return face_renumber;
+      }
 
-    // dummy specialization for dim == 1 to avoid linker errors
-    template <>
-    inline
-    std::vector<unsigned int>
-    face_lexicographic_to_hierarchic_numbering<1> (const unsigned int)
-    {
-      return std::vector<unsigned int>();
-    }
+      // dummy specialization for dim == 1 to avoid linker errors
+      template <>
+      inline
+      std::vector<unsigned int>
+      face_lexicographic_to_hierarchic_numbering<1> (const unsigned int)
+      {
+        return std::vector<unsigned int>();
+      }
 
 
 
-    // in get_restriction_matrix() and get_prolongation_matrix(), want to undo
-    // tensorization on inner loops for performance reasons. this clears a
-    // dim-array
-    template <int dim>
-    inline
-    void
-    zero_indices (unsigned int (&indices)[dim])
-    {
-      for (unsigned int d=0; d<dim; ++d)
-        indices[d] = 0;
-    }
+      // in get_restriction_matrix() and get_prolongation_matrix(), want to undo
+      // tensorization on inner loops for performance reasons. this clears a
+      // dim-array
+      template <int dim>
+      inline
+      void
+      zero_indices (unsigned int (&indices)[dim])
+      {
+        for (unsigned int d=0; d<dim; ++d)
+          indices[d] = 0;
+      }
 
 
 
-    // in get_restriction_matrix() and get_prolongation_matrix(), want to undo
-    // tensorization on inner loops for performance reasons. this increments
-    // tensor product indices
-    template <int dim>
-    inline
-    void
-    increment_indices (unsigned int       (&indices)[dim],
-                       const unsigned int   dofs1d)
-    {
-      ++indices[0];
-      for (int d=0; d<dim-1; ++d)
-        if (indices[d]==dofs1d)
-          {
-            indices[d] = 0;
-            indices[d+1]++;
-          }
+      // in get_restriction_matrix() and get_prolongation_matrix(), want to undo
+      // tensorization on inner loops for performance reasons. this increments
+      // tensor product indices
+      template <int dim>
+      inline
+      void
+      increment_indices (unsigned int       (&indices)[dim],
+                         const unsigned int   dofs1d)
+      {
+        ++indices[0];
+        for (int d=0; d<dim-1; ++d)
+          if (indices[d]==dofs1d)
+            {
+              indices[d] = 0;
+              indices[d+1]++;
+            }
+      }
     }
   }
 }
@@ -205,7 +208,7 @@ struct FE_Q_Base<PolynomialType,xdim,xspacedim>::Implementation
     const std::vector<unsigned int> &index_map_inverse =
       fe.poly_space.get_numbering_inverse();
     const std::vector<unsigned int> face_index_map =
-      FE_Q_Helper::face_lexicographic_to_hierarchic_numbering<dim>(q_deg);
+      internal::FE_Q_Base::face_lexicographic_to_hierarchic_numbering<dim>(q_deg);
     Assert(std::abs(fe.poly_space.compute_value(index_map_inverse[0],Point<dim>())
                     - 1.) < 1e-14,
            ExcInternalError());
@@ -329,7 +332,7 @@ struct FE_Q_Base<PolynomialType,xdim,xspacedim>::Implementation
     const std::vector<unsigned int> &index_map_inverse =
       fe.poly_space.get_numbering_inverse();
     const std::vector<unsigned int> face_index_map =
-      FE_Q_Helper::face_lexicographic_to_hierarchic_numbering<dim>(q_deg);
+      internal::FE_Q_Base::face_lexicographic_to_hierarchic_numbering<dim>(q_deg);
     Assert(std::abs(fe.poly_space.compute_value(index_map_inverse[0],Point<dim>())
                     - 1.) < 1e-14,
            ExcInternalError());
@@ -928,7 +931,7 @@ void FE_Q_Base<PolynomialType,dim,spacedim>::initialize_unit_face_support_points
 
   // find renumbering of faces and assign from values of quadrature
   std::vector<unsigned int> face_index_map =
-    FE_Q_Helper::face_lexicographic_to_hierarchic_numbering<dim>(q_degree);
+    internal::FE_Q_Base::face_lexicographic_to_hierarchic_numbering<dim>(q_degree);
   Quadrature<1> support_1d(points);
   Quadrature<codim> support_quadrature(support_1d);
   this->unit_face_support_points.resize(support_quadrature.size());
@@ -1284,11 +1287,11 @@ FE_Q_Base<PolynomialType,dim,spacedim>
       // now expand from 1D info. block innermost dimension (x_0) in order to
       // avoid difficult checks at innermost loop
       unsigned int j_indices[dim];
-      FE_Q_Helper::zero_indices<dim> (j_indices);
+      internal::FE_Q_Base::zero_indices<dim> (j_indices);
       for (unsigned int j=0; j<q_dofs_per_cell; j+=dofs1d)
         {
           unsigned int i_indices[dim];
-          FE_Q_Helper::zero_indices<dim> (i_indices);
+          internal::FE_Q_Base::zero_indices<dim> (i_indices);
           for (unsigned int i=0; i<q_dofs_per_cell; i+=dofs1d)
             {
               double val_extra_dim = 1.;
@@ -1308,10 +1311,10 @@ FE_Q_Base<PolynomialType,dim,spacedim>
 
               // update indices that denote the tensor product position. a bit
               // fuzzy and therefore not done for innermost x_0 direction
-              FE_Q_Helper::increment_indices<dim> (i_indices, dofs1d);
+              internal::FE_Q_Base::increment_indices<dim> (i_indices, dofs1d);
             }
           Assert (i_indices[dim-1] == 1, ExcInternalError());
-          FE_Q_Helper::increment_indices<dim> (j_indices, dofs1d);
+          internal::FE_Q_Base::increment_indices<dim> (j_indices, dofs1d);
         }
 
       // the discontinuous node is simply mapped on the discontinuous node on
@@ -1428,7 +1431,7 @@ FE_Q_Base<PolynomialType,dim,spacedim>
                       this->poly_space.compute_value(index_map_inverse[j], point);
                   }
               unsigned int j_indices[dim];
-              FE_Q_Helper::zero_indices<dim> (j_indices);
+              internal::FE_Q_Base::zero_indices<dim> (j_indices);
               double sum_check = 0;
               for (unsigned int j = 0; j<q_dofs_per_cell; j += dofs1d)
                 {
@@ -1454,7 +1457,7 @@ FE_Q_Base<PolynomialType,dim,spacedim>
                         my_restriction(mother_dof,child_dof)=val;
                       sum_check += val;
                     }
-                  FE_Q_Helper::increment_indices<dim> (j_indices, dofs1d);
+                  internal::FE_Q_Base::increment_indices<dim> (j_indices, dofs1d);
                 }
               Assert (std::fabs(sum_check-1) <
                       std::max(eps, 5e-16*std::sqrt(this->dofs_per_cell)),
diff --git a/source/fe/fe_q_bubbles.cc b/source/fe/fe_q_bubbles.cc
index c82b9cd83d..4312fb0542 100644
--- a/source/fe/fe_q_bubbles.cc
+++ b/source/fe/fe_q_bubbles.cc
@@ -38,131 +38,135 @@
 DEAL_II_NAMESPACE_OPEN
 
 
-namespace FE_Q_Bubbles_Helper
+namespace internal
 {
-  namespace
+  namespace FE_Q_Bubbles
   {
-    template <int dim, int spacedim>
-    inline
-    void
-    compute_embedding_matrices(const FE_Q_Bubbles<dim, spacedim> &fe,
-                               std::vector<std::vector<FullMatrix<double> > > &matrices,
-                               const bool isotropic_only)
+    namespace
     {
-      const unsigned int dpc    = fe.dofs_per_cell;
-      const unsigned int degree = fe.degree;
+      template <int dim, int spacedim>
+      inline
+      void
+      compute_embedding_matrices(const dealii::FE_Q_Bubbles<dim, spacedim> &fe,
+                                 std::vector<std::vector<FullMatrix<double> > > &matrices,
+                                 const bool                                 isotropic_only)
+      {
+        const unsigned int dpc    = fe.dofs_per_cell;
+        const unsigned int degree = fe.degree;
 
-      // Initialize quadrature formula on fine cells
-      std::unique_ptr<Quadrature<dim> > q_fine;
-      Quadrature<1> q_dummy(std::vector<Point<1> >(1), std::vector<double> (1,1.));
-      switch (dim)
-        {
-        case 1:
-          if (spacedim==1)
-            q_fine.reset(new QGauss<dim> (degree+1));
-          else if (spacedim==2)
-            q_fine.reset(new QAnisotropic<dim>(QGauss<1>(degree+1), q_dummy));
-          else
-            q_fine.reset(new QAnisotropic<dim>(QGauss<1>(degree+1), q_dummy, q_dummy));
-          break;
-        case 2:
-          if (spacedim==2)
+        // Initialize quadrature formula on fine cells
+        std::unique_ptr<Quadrature<dim> > q_fine;
+        Quadrature<1> q_dummy(std::vector<Point<1> >(1), std::vector<double> (1,1.));
+        switch (dim)
+          {
+          case 1:
+            if (spacedim==1)
+              q_fine.reset(new QGauss<dim> (degree+1));
+            else if (spacedim==2)
+              q_fine.reset(new QAnisotropic<dim>(QGauss<1>(degree+1), q_dummy));
+            else
+              q_fine.reset(new QAnisotropic<dim>(QGauss<1>(degree+1), q_dummy, q_dummy));
+            break;
+          case 2:
+            if (spacedim==2)
+              q_fine.reset(new QGauss<dim> (degree+1));
+            else
+              q_fine.reset(new QAnisotropic<dim>(QGauss<1>(degree+1), QGauss<1>(degree+1), q_dummy));
+            break;
+          case 3:
             q_fine.reset(new QGauss<dim> (degree+1));
-          else
-            q_fine.reset(new QAnisotropic<dim>(QGauss<1>(degree+1), QGauss<1>(degree+1), q_dummy));
-          break;
-        case 3:
-          q_fine.reset(new QGauss<dim> (degree+1));
-          break;
-        default:
-          Assert(false, ExcInternalError());
-        }
+            break;
+          default:
+            Assert(false, ExcInternalError());
+          }
 
-      Assert(q_fine.get() != nullptr, ExcInternalError());
-      const unsigned int nq = q_fine->size();
+        Assert(q_fine.get() != nullptr, ExcInternalError());
+        const unsigned int nq = q_fine->size();
 
-      // loop over all possible refinement cases
-      unsigned int ref_case = (isotropic_only)
-                              ? RefinementCase<dim>::isotropic_refinement
-                              : RefinementCase<dim>::cut_x;
-      for (; ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
-        {
-          const unsigned int nc
-            = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
-
-          for (unsigned int i=0; i<nc; ++i)
-            {
-              Assert(matrices[ref_case-1][i].n() == dpc,
-                     ExcDimensionMismatch(matrices[ref_case-1][i].n(),dpc));
-              Assert(matrices[ref_case-1][i].m() == dpc,
-                     ExcDimensionMismatch(matrices[ref_case-1][i].m(),dpc));
-            }
-
-          // create a respective refinement on the triangulation
-          Triangulation<dim, spacedim> tr;
-          GridGenerator::hyper_cube (tr, 0, 1);
-          tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
-          tr.execute_coarsening_and_refinement();
-
-          DoFHandler<dim, spacedim> dh(tr);
-          dh.distribute_dofs(fe);
-
-          FEValues<dim, spacedim> fine (StaticMappingQ1<dim,spacedim>::mapping, fe, *q_fine,
-                                        update_quadrature_points
-                                        | update_JxW_values | update_values);
-
-          const unsigned int n_dofs = dh.n_dofs();
-
-          FullMatrix<double> fine_mass(n_dofs);
-          FullMatrix<double> coarse_rhs_matrix(n_dofs, dpc);
-
-          std::vector<std::vector<types::global_dof_index> > child_ldi
-          (nc, std::vector<types::global_dof_index>(fe.dofs_per_cell));
-
-          //now create the mass matrix and all the right_hand sides
-          unsigned int child_no = 0;
-          typename dealii::DoFHandler<dim>::active_cell_iterator cell
-            = dh.begin_active();
-          for (; cell!=dh.end(); ++cell, ++child_no)
-            {
-              fine.reinit(cell);
-              cell->get_dof_indices(child_ldi[child_no]);
-
-              for (unsigned int q=0; q<nq; ++q)
-                for (unsigned int i=0; i<dpc; ++i)
-                  for (unsigned int j=0; j<dpc; ++j)
-                    {
-                      const unsigned int gdi=child_ldi[child_no][i];
-                      const unsigned int gdj=child_ldi[child_no][j];
-                      fine_mass(gdi, gdj)+=fine.shape_value(i,q)
-                                           *fine.shape_value(j,q)
-                                           *fine.JxW(q);
-                      Point<dim> quad_tmp;
-                      for (unsigned int k=0; k<dim; ++k)
-                        quad_tmp(k) = fine.quadrature_point(q)(k);
-                      coarse_rhs_matrix(gdi, j)
-                      +=fine.shape_value(i,q)
-                        *fe.shape_value(j, quad_tmp)
-                        *fine.JxW(q);
-                    }
-            }
-
-          //now solve for all right-hand sides simultaneously
-          dealii::FullMatrix<double> solution (n_dofs, dpc);
-          fine_mass.gauss_jordan();
-          fine_mass.mmult(solution, coarse_rhs_matrix);
-
-          //and distribute to the fine cell matrices
-          for (unsigned int child_no=0; child_no<nc; ++child_no)
-            for (unsigned int i=0; i<dpc; ++i)
-              for (unsigned int j=0; j<dpc; ++j)
-                {
-                  const unsigned int gdi=child_ldi[child_no][i];
-                  //remove small entries
-                  if (std::fabs(solution(gdi, j)) > 1.e-12)
-                    matrices[ref_case-1][child_no](i,j)=solution(gdi, j);
-                }
-        }
+        // loop over all possible refinement cases
+        unsigned int ref_case = (isotropic_only)
+                                ? RefinementCase<dim>::isotropic_refinement
+                                : RefinementCase<dim>::cut_x;
+        for (; ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
+          {
+            const unsigned int nc
+              = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
+
+            for (unsigned int i=0; i<nc; ++i)
+              {
+                Assert(matrices[ref_case-1][i].n() == dpc,
+                       ExcDimensionMismatch(matrices[ref_case-1][i].n(),dpc));
+                Assert(matrices[ref_case-1][i].m() == dpc,
+                       ExcDimensionMismatch(matrices[ref_case-1][i].m(),dpc));
+              }
+
+            // create a respective refinement on the triangulation
+            dealii::Triangulation<dim, spacedim> tr;
+            GridGenerator::hyper_cube (tr, 0, 1);
+            tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
+            tr.execute_coarsening_and_refinement();
+
+            dealii::DoFHandler<dim, spacedim> dh(tr);
+            dh.distribute_dofs(fe);
+
+            dealii::FEValues<dim, spacedim> fine
+            (StaticMappingQ1<dim,spacedim>::mapping, fe, *q_fine,
+             update_quadrature_points
+             | update_JxW_values | update_values);
+
+            const unsigned int n_dofs = dh.n_dofs();
+
+            FullMatrix<double> fine_mass(n_dofs);
+            FullMatrix<double> coarse_rhs_matrix(n_dofs, dpc);
+
+            std::vector<std::vector<types::global_dof_index> > child_ldi
+            (nc, std::vector<types::global_dof_index>(fe.dofs_per_cell));
+
+            //now create the mass matrix and all the right_hand sides
+            unsigned int child_no = 0;
+            typename dealii::DoFHandler<dim>::active_cell_iterator cell
+              = dh.begin_active();
+            for (; cell!=dh.end(); ++cell, ++child_no)
+              {
+                fine.reinit(cell);
+                cell->get_dof_indices(child_ldi[child_no]);
+
+                for (unsigned int q=0; q<nq; ++q)
+                  for (unsigned int i=0; i<dpc; ++i)
+                    for (unsigned int j=0; j<dpc; ++j)
+                      {
+                        const unsigned int gdi=child_ldi[child_no][i];
+                        const unsigned int gdj=child_ldi[child_no][j];
+                        fine_mass(gdi, gdj)+=fine.shape_value(i,q)
+                                             *fine.shape_value(j,q)
+                                             *fine.JxW(q);
+                        Point<dim> quad_tmp;
+                        for (unsigned int k=0; k<dim; ++k)
+                          quad_tmp(k) = fine.quadrature_point(q)(k);
+                        coarse_rhs_matrix(gdi, j)
+                        +=fine.shape_value(i,q)
+                          *fe.shape_value(j, quad_tmp)
+                          *fine.JxW(q);
+                      }
+              }
+
+            //now solve for all right-hand sides simultaneously
+            dealii::FullMatrix<double> solution (n_dofs, dpc);
+            fine_mass.gauss_jordan();
+            fine_mass.mmult(solution, coarse_rhs_matrix);
+
+            //and distribute to the fine cell matrices
+            for (unsigned int child_no=0; child_no<nc; ++child_no)
+              for (unsigned int i=0; i<dpc; ++i)
+                for (unsigned int j=0; j<dpc; ++j)
+                  {
+                    const unsigned int gdi=child_ldi[child_no][i];
+                    //remove small entries
+                    if (std::fabs(solution(gdi, j)) > 1.e-12)
+                      matrices[ref_case-1][child_no](i,j)=solution(gdi, j);
+                  }
+          }
+      }
     }
   }
 }
@@ -196,7 +200,7 @@ FE_Q_Bubbles<dim,spacedim>::FE_Q_Bubbles (const unsigned int q_degree)
   this->reinit_restriction_and_prolongation_matrices();
   if (dim == spacedim)
     {
-      FE_Q_Bubbles_Helper::compute_embedding_matrices
+      internal::FE_Q_Bubbles::compute_embedding_matrices
       (*this, this->prolongation, false);
       // Fill restriction matrices with L2-projection
       FETools::compute_projection_matrices (*this, this->restriction);
@@ -234,7 +238,7 @@ FE_Q_Bubbles<dim,spacedim>::FE_Q_Bubbles (const Quadrature<1> &points)
   this->reinit_restriction_and_prolongation_matrices();
   if (dim == spacedim)
     {
-      FE_Q_Bubbles_Helper::compute_embedding_matrices
+      internal::FE_Q_Bubbles::compute_embedding_matrices
       (*this, this->prolongation, false);
       // Fill restriction matrices with L2-projection
       FETools::compute_projection_matrices (*this, this->restriction);
diff --git a/source/fe/fe_q_hierarchical.cc b/source/fe/fe_q_hierarchical.cc
index 2e55a25841..184e0a5c71 100644
--- a/source/fe/fe_q_hierarchical.cc
+++ b/source/fe/fe_q_hierarchical.cc
@@ -28,23 +28,29 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-namespace
+namespace internal
 {
-  /**
-   * A function which maps  in[i] to i,i.e. output[in[i]] = i;
-   */
-  inline
-  std::vector<unsigned int>
-  invert_numbering (const std::vector<unsigned int> &in)
+  namespace FE_Q_Hierarchical
   {
-    std::vector<unsigned int> out (in.size());
-    for (unsigned int i=0; i<in.size(); ++i)
+    namespace
+    {
+      /**
+       * A function which maps  in[i] to i,i.e. output[in[i]] = i;
+       */
+      inline
+      std::vector<unsigned int>
+      invert_numbering (const std::vector<unsigned int> &in)
       {
-        Assert (in[i] < out.size(),
-                dealii::ExcIndexRange(in[i],0,out.size()));
-        out[in[i]]=i;
+        std::vector<unsigned int> out (in.size());
+        for (unsigned int i=0; i<in.size(); ++i)
+          {
+            Assert (in[i] < out.size(),
+                    dealii::ExcIndexRange(in[i],0,out.size()));
+            out[in[i]]=i;
+          }
+        return out;
       }
-    return out;
+    }
   }
 }
 
@@ -1812,8 +1818,8 @@ FE_Q_Hierarchical<dim>::
 face_fe_q_hierarchical_to_hierarchic_numbering (const unsigned int degree)
 {
   FiniteElementData<dim-1> fe_data(FE_Q_Hierarchical<dim-1>::get_dpo_vector(degree),1,degree);
-  return invert_numbering(FE_Q_Hierarchical<dim-1>::
-                          hierarchic_to_fe_q_hierarchical_numbering (fe_data));
+  return internal::FE_Q_Hierarchical::invert_numbering
+         (FE_Q_Hierarchical<dim-1>::hierarchic_to_fe_q_hierarchical_numbering (fe_data));
 }
 
 
diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc
index 0f93679872..40ee23cd78 100644
--- a/source/fe/fe_system.cc
+++ b/source/fe/fe_system.cc
@@ -30,18 +30,21 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-namespace
+namespace internal
 {
-  bool IsNonZero (unsigned int i)
+  namespace FESystem
   {
-    return i>0;
-  }
-
-  unsigned int count_nonzeros(const std::vector<unsigned int> &vec)
-  {
-    return std::count_if(vec.begin(), vec.end(), IsNonZero);
+    namespace
+    {
+      unsigned int count_nonzeros(const std::vector<unsigned int> &vec)
+      {
+        return std::count_if(vec.begin(), vec.end(), [](const unsigned int i)
+        {
+          return i > 0;
+        });
+      }
+    }
   }
-
 }
 /* ----------------------- FESystem::InternalData ------------------- */
 
@@ -275,7 +278,7 @@ FESystem<dim,spacedim>::FESystem (
   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))
+  base_elements(internal::FESystem::count_nonzeros(multiplicities))
 {
   initialize(fes, multiplicities);
 }
@@ -1458,7 +1461,7 @@ void FESystem<dim,spacedim>::initialize (const std::vector<const FiniteElement<d
           ExcDimensionMismatch (fes.size(), multiplicities.size()) );
   Assert (fes.size() > 0,
           ExcMessage ("Need to pass at least one finite element."));
-  Assert (count_nonzeros(multiplicities) > 0,
+  Assert (internal::FESystem::count_nonzeros(multiplicities) > 0,
           ExcMessage("You only passed FiniteElements with multiplicity 0."));
 
   // Note that we need to skip every fe with multiplicity 0 in the following block of code
-- 
2.39.5