]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Replaced projection-based interpolation by explicit interpolation formulas wherever...
authorbuerg <buerg@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 22 Jul 2012 15:21:35 +0000 (15:21 +0000)
committerbuerg <buerg@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 22 Jul 2012 15:21:35 +0000 (15:21 +0000)
git-svn-id: https://svn.dealii.org/trunk@25717 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/source/fe/fe_nedelec.cc

index e1c05bb37bc219c9f79e2929b04476c87fb98933..b1da9ccd6e9e7665ef91b45b64a7d4a913ca5248 100644 (file)
 #include <sstream>
 #include <iostream>
 
-//TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
-//adjust_line_dof_index_for_line_orientation_table fields, and write tests
-//similar to bits/face_orientation_and_fe_q_*
-
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -161,7 +157,6 @@ FE_PolyTensor<PolynomialsNedelec<dim>, dim>
       default:
         Assert (false, ExcNotImplemented ());
     }
-
 }
 
 
@@ -497,6 +492,8 @@ FE_Nedelec<1>::initialize_restriction ()
 
 
 
+
+
                    // Restriction operator
 template <int dim>
 void
@@ -515,6 +512,9 @@ FE_Nedelec<dim>::initialize_restriction ()
     n_edge_quadrature_points = edge_quadrature.size ();
   const unsigned int
     index = RefinementCase<dim>::isotropic_refinement - 1;
+  const unsigned int deg = this->degree-1;
+  const std::vector<Polynomials::Polynomial<double> >&
+    legendre_polynomials = Polynomials::Legendre::generate_complete_basis (deg);
 
   switch (dim)
     {
@@ -522,232 +522,83 @@ FE_Nedelec<dim>::initialize_restriction ()
         {
                         // First interpolate the shape
                         // functions of the child cells
-                        // to the lowest order shape
-                        // functions of the parent cell.
+                        // to the face shape functions
+                        // of the parent cell.
           for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
             for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
                  ++q_point)
               {
-                const double weight = 2.0 * edge_quadrature.weight (q_point);
-
                 if (edge_quadrature_points[q_point] (0) < 0.5)
-                  {
-                    Point<dim> quadrature_point (0.0,
-                                                 2.0 * edge_quadrature_points[q_point] (0));
-
-                    this->restriction[index][0] (0, dof) += weight
-                                                         * this->shape_value_component
-                                                           (dof,
-                                                            quadrature_point,
-                                                            1);
-                    quadrature_point (0) = 1.0;
-                    this->restriction[index][1] (this->degree, dof)
-                      += weight * this->shape_value_component (dof,
-                                                               quadrature_point,
-                                                               1);
-                    quadrature_point (0) = quadrature_point (1);
-                    quadrature_point (1) = 0.0;
-                    this->restriction[index][0] (2 * this->degree, dof)
-                      += weight * this->shape_value_component (dof,
-                                                               quadrature_point,
-                                                               0);
-                    quadrature_point (1) = 1.0;
-                    this->restriction[index][2] (3 * this->degree, dof)
-                      += weight * this->shape_value_component (dof,
-                                                               quadrature_point,
-                                                               0);
-                  }
+                  for (unsigned int i = 0; i <= deg; ++i)
+                    {
+                      const double weight
+                        = 2.0 * edge_quadrature.weight (q_point)
+                              * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0));
+                      Point<dim> quadrature_point (0.0,
+                                                   2.0 * edge_quadrature_points[q_point] (0));
+
+                      this->restriction[index][0] (i, dof) += weight
+                                                           * this->shape_value_component
+                                                             (dof,
+                                                              quadrature_point,
+                                                              1);
+                      quadrature_point (0) = 1.0;
+                      this->restriction[index][1] (i + this->degree, dof)
+                        += weight * this->shape_value_component (dof,
+                                                                 quadrature_point,
+                                                                 1);
+                      quadrature_point (0) = quadrature_point (1);
+                      quadrature_point (1) = 0.0;
+                      this->restriction[index][0] (i + 2 * this->degree, dof)
+                        += weight * this->shape_value_component (dof,
+                                                                 quadrature_point,
+                                                                 0);
+                      quadrature_point (1) = 1.0;
+                      this->restriction[index][2] (i + 3 * this->degree, dof)
+                        += weight * this->shape_value_component (dof,
+                                                                 quadrature_point,
+                                                                 0);
+                    }
 
                 else
-                  {
-                    Point<dim> quadrature_point (0.0,
-                                                 2.0 * edge_quadrature_points[q_point] (0)
-                                                     - 1.0);
-
-                    this->restriction[index][2] (0, dof) += weight
-                                                         * this->shape_value_component
-                                                           (dof,
-                                                            quadrature_point,
-                                                            1);
-                    quadrature_point (0) = 1.0;
-                    this->restriction[index][3] (this->degree, dof)
-                      += weight * this->shape_value_component (dof,
-                                                               quadrature_point,
-                                                               1);
-                    quadrature_point (0) = quadrature_point (1);
-                    quadrature_point (1) = 0.0;
-                    this->restriction[index][1] (2 * this->degree, dof)
-                      += weight * this->shape_value_component (dof,
-                                                               quadrature_point,
-                                                               0);
-                    quadrature_point (1) = 1.0;
-                    this->restriction[index][3] (3 * this->degree, dof)
-                      += weight * this->shape_value_component (dof,
-                                                               quadrature_point,
-                                                               0);
-                  }
+                  for (unsigned int i = 0; i <= deg; ++i)
+                    {
+                      const double weight
+                        = 2.0 * edge_quadrature.weight (q_point)
+                              * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0));
+                      Point<dim> quadrature_point (0.0,
+                                                   2.0 * edge_quadrature_points[q_point] (0)
+                                                       - 1.0);
+
+                      this->restriction[index][2] (i, dof) += weight
+                                                           * this->shape_value_component
+                                                             (dof,
+                                                              quadrature_point,
+                                                              1);
+                      quadrature_point (0) = 1.0;
+                      this->restriction[index][3] (i + this->degree, dof)
+                        += weight * this->shape_value_component (dof,
+                                                                 quadrature_point,
+                                                                 1);
+                      quadrature_point (0) = quadrature_point (1);
+                      quadrature_point (1) = 0.0;
+                      this->restriction[index][1] (i + 2 * this->degree, dof)
+                        += weight * this->shape_value_component (dof,
+                                                                 quadrature_point,
+                                                                 0);
+                      quadrature_point (1) = 1.0;
+                      this->restriction[index][3] (i + 3 * this->degree, dof)
+                        += weight * this->shape_value_component (dof,
+                                                                 quadrature_point,
+                                                                 0);
+                    }
               }
 
-                        // Then project the shape functions
-                        // of the child cells to the higher
-                        // order shape functions of the
-                        // parent cell.
+                        // Then interpolate the shape functions
+                        // of the child cells to the interior
+                        // shape functions of the parent cell.
           if (this->degree > 1)
             {
-              const unsigned int deg = this->degree-1;
-              const std::vector<Polynomials::Polynomial<double> >&
-                legendre_polynomials
-                  = Polynomials::Legendre::generate_complete_basis (deg);
-              FullMatrix<double> system_matrix_inv (deg, deg);
-
-              {
-                FullMatrix<double> assembling_matrix (deg,
-                                                      n_edge_quadrature_points);
-
-                for (unsigned int q_point = 0;
-                     q_point < n_edge_quadrature_points; ++q_point)
-                  {
-                    const double weight
-                      = std::sqrt (edge_quadrature.weight (q_point));
-
-                    for (unsigned int i = 0; i < deg; ++i)
-                      assembling_matrix (i, q_point) = weight
-                                                       * legendre_polynomials[i + 1].value
-                                                         (edge_quadrature_points[q_point] (0));
-                  }
-
-                FullMatrix<double> system_matrix (deg, deg);
-
-                assembling_matrix.mTmult (system_matrix, assembling_matrix);
-                system_matrix_inv.invert (system_matrix);
-              }
-
-              FullMatrix<double> solution (this->degree-1, 4);
-              FullMatrix<double> system_rhs (this->degree-1, 4);
-              Vector<double> tmp (4);
-
-              for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
-                for (unsigned int i = 0; i < 2; ++i)
-                  {
-                    system_rhs = 0.0;
-
-                    for (unsigned int q_point = 0;
-                         q_point < n_edge_quadrature_points; ++q_point)
-                      {
-                        const double weight
-                          = edge_quadrature.weight (q_point);
-                        const Point<dim> quadrature_point_0 (i,
-                                                                 edge_quadrature_points[q_point] (0));
-                        const Point<dim> quadrature_point_1
-                                             (edge_quadrature_points[q_point] (0),
-                                              i);
-
-                        if (edge_quadrature_points[q_point] (0) < 0.5)
-                          {
-                            Point<dim> quadrature_point_2 (i,
-                                                                   2.0 * edge_quadrature_points[q_point] (0));
-
-                            tmp (0) = weight
-                                      * (2.0 * this->shape_value_component
-                                                   (dof, quadrature_point_2, 1)
-                                                 - this->restriction[index][i]
-                                                   (i * this->degree, dof)
-                                                 * this->shape_value_component
-                                                   (i * this->degree,
-                                                    quadrature_point_0, 1));
-                            tmp (1) = -1.0 * weight
-                                               * this->restriction[index][i + 2]
-                                                 (i * this->degree, dof)
-                                               * this->shape_value_component
-                                                 (i * this->degree,
-                                                  quadrature_point_0, 1);
-                            quadrature_point_2
-                              = Point<dim> (2.0 * edge_quadrature_points[q_point] (0),
-                                                i);
-                            tmp (2) = weight
-                                      * (2.0 * this->shape_value_component
-                                                   (dof, quadrature_point_2, 0)
-                                                 - this->restriction[index][2 * i]
-                                                   ((i + 2) * this->degree, dof)
-                                                 * this->shape_value_component
-                                                   ((i + 2) * this->degree,
-                                                    quadrature_point_1, 0));
-                            tmp (3) = -1.0 * weight
-                                               * this->restriction[index][2 * i + 1]
-                                                 ((i + 2) * this->degree, dof)
-                                               * this->shape_value_component
-                                                 ((i + 2) * this->degree,
-                                                  quadrature_point_1, 0);
-                          }
-
-                        else
-                          {
-                            tmp (0) = -1.0 * weight
-                                           * this->restriction[index][i]
-                                                 (i * this->degree, dof)
-                                               * this->shape_value_component
-                                                 (i * this->degree,
-                                                  quadrature_point_0, 1);
-
-                            Point<dim> quadrature_point_2 (i,
-                                                                   2.0 * edge_quadrature_points[q_point] (0)
-                                                                       - 1.0);
-
-                            tmp (1) = weight
-                                          * (2.0 * this->shape_value_component
-                                                   (dof, quadrature_point_2, 1)
-                                                 - this->restriction[index][i + 2]
-                                                   (i * this->degree, dof)
-                                                 * this->shape_value_component
-                                                   (i * this->degree,
-                                                    quadrature_point_0, 1));
-                            tmp (2) = -1.0 * weight
-                                               * this->restriction[index][2 * i]
-                                                 ((i + 2) * this->degree, dof)
-                                               * this->shape_value_component
-                                                 ((i + 2) * this->degree,
-                                                  quadrature_point_1, 0);
-                            quadrature_point_2
-                              = Point<dim> (2.0 * edge_quadrature_points[q_point] (0)
-                                                    - 1.0, i);
-                            tmp (3) = weight
-                                          * (2.0 * this->shape_value_component
-                                                   (dof, quadrature_point_2, 0)
-                                                 - this->restriction[index][2 * i + 1]
-                                                   ((i + 2) * this->degree, dof)
-                                                 * this->shape_value_component
-                                                   ((i + 2) * this->degree,
-                                                    quadrature_point_1, 0));
-                          }
-
-                        for (unsigned int j = 0; j < this->degree-1; ++j)
-                          {
-                            const double L_j
-                              = legendre_polynomials[j + 1].value
-                                (edge_quadrature_points[q_point] (0));
-
-                            for (unsigned int k = 0; k < tmp.size (); ++k)
-                              system_rhs (j, k) += tmp (k) * L_j;
-                          }
-                      }
-
-                    system_matrix_inv.mmult (solution, system_rhs);
-
-                    for (unsigned int j = 0; j < this->degree-1; ++j)
-                      for (unsigned int k = 0; k < 2; ++k)
-                        {
-                          if (std::abs (solution (j, k)) > 1e-14)
-                            this->restriction[index][i + 2 * k]
-                            (i * this->degree + j + 1, dof)
-                              = solution (j, k);
-
-                          if (std::abs (solution (j, k + 2)) > 1e-14)
-                            this->restriction[index][2 * i + k]
-                            ((i + 2) * this->degree + j + 1, dof)
-                              = solution (j, k + 2);
-                        }
-                  }
-
               const QGauss<dim> quadrature (2 * this->degree);
               const std::vector<Point<dim> >&
                 quadrature_points = quadrature.get_points ();
@@ -757,185 +608,71 @@ FE_Nedelec<dim>::initialize_restriction ()
                     (this->degree);
               const unsigned int n_boundary_dofs
                 = GeometryInfo<dim>::faces_per_cell * this->degree;
-              const unsigned int& n_quadrature_points = quadrature.size ();
 
-              {
-                FullMatrix<double> assembling_matrix ((this->degree-1) * this->degree,
-                                                      n_quadrature_points);
-
-                for (unsigned int q_point = 0; q_point < n_quadrature_points;
+              for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
+                for (unsigned int q_point = 0; q_point < quadrature.size ();
                      ++q_point)
                   {
-                    const double weight
-                      = std::sqrt (quadrature.weight (q_point));
+                    const double weight = 2.0 * quadrature.weight (q_point);
+                    Point<dim> quadrature_point (2.0 * quadrature_points[q_point] (0),
+                                                 2.0 * quadrature_points[q_point] (1));
+                    unsigned int k;
 
-                    for (unsigned int i = 0; i < this->degree; ++i)
+                    if (quadrature_points[q_point] (0) < 0.5)
                       {
-                        const double L_i = weight
-                                           * legendre_polynomials[i].value
-                                             (quadrature_points[q_point] (0));
-
-                        for (unsigned int j = 0; j < this->degree-1; ++j)
-                          assembling_matrix (i * (this->degree-1) + j, q_point)
-                            = L_i * lobatto_polynomials[j + 2].value
-                                    (quadrature_points[q_point] (1));
-                      }
-                  }
-
-                FullMatrix<double> system_matrix (assembling_matrix.m (),
-                                                  assembling_matrix.m ());
-
-                assembling_matrix.mTmult (system_matrix, assembling_matrix);
-                system_matrix_inv.reinit (system_matrix.m (), system_matrix.m ());
-                system_matrix_inv.invert (system_matrix);
-              }
-
-              solution.reinit (system_matrix_inv.m (), 8);
-              system_rhs.reinit (system_matrix_inv.m (), 8);
-              tmp.reinit (8);
-
-              for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
-                {
-                  system_rhs = 0.0;
-
-                  for (unsigned int q_point = 0;
-                       q_point < n_quadrature_points; ++q_point)
-                    {
-                      tmp = 0.0;
-
-                      if (quadrature_points[q_point] (0) < 0.5)
-                        {
-                          if (quadrature_points[q_point] (1) < 0.5)
-                            {
-                              const Point<dim> quadrature_point
-                                               (2.0 * quadrature_points[q_point] (0),
-                                                2.0 * quadrature_points[q_point] (1));
-
-                              tmp (0) += 2.0 * this->shape_value_component
-                                               (dof, quadrature_point, 0);
-                              tmp (1) += 2.0 * this->shape_value_component
-                                               (dof, quadrature_point, 1);
-                            }
-
-                          else
-                            {
-                              const Point<dim> quadrature_point
-                                               (2.0 * quadrature_points[q_point] (0),
-                                                2.0 * quadrature_points[q_point] (1)
-                                                    - 1.0);
-
-                              tmp (4) += 2.0 * this->shape_value_component
-                                               (dof, quadrature_point, 0);
-                              tmp (5) += 2.0 * this->shape_value_component
-                                               (dof, quadrature_point, 1);
-                            }
-                        }
-
-                      else
                         if (quadrature_points[q_point] (1) < 0.5)
-                          {
-                            const Point<dim> quadrature_point
-                                             (2.0 * quadrature_points[q_point] (0)
-                                                  - 1.0,
-                                              2.0 * quadrature_points[q_point] (1));
-
-                            tmp (2) += 2.0 * this->shape_value_component
-                                             (dof, quadrature_point, 0);
-                            tmp (3) += 2.0 * this->shape_value_component
-                                             (dof, quadrature_point, 1);
-                          }
-
+                          k = 0;
+                        
                         else
                           {
-                            const Point<dim> quadrature_point
-                                             (2.0 * quadrature_points[q_point] (0)
-                                                  - 1.0,
-                                              2.0 * quadrature_points[q_point] (1)
-                                                  - 1.0);
-
-                            tmp (6) += 2.0 * this->shape_value_component
-                                             (dof, quadrature_point, 0);
-                            tmp (7) += 2.0 * this->shape_value_component
-                                             (dof, quadrature_point, 1);
-                          }
-
-                      for (unsigned int i = 0; i < 2; ++i)
-                        for (unsigned int j = 0; j < this->degree; ++j)
-                          {
-                            tmp (2 * i) -= this->restriction[index][i]
-                                           (j + 2 * this->degree, dof)
-                                        * this->shape_value_component
-                                          (j + 2 * this->degree,
-                                           quadrature_points[q_point], 0);
-                            tmp (2 * i + 1) -= this->restriction[index][i]
-                                               (i * this->degree + j, dof)
-                                            * this->shape_value_component
-                                              (i * this->degree + j,
-                                               quadrature_points[q_point], 1);
-                            tmp (2 * (i + 2)) -= this->restriction[index][i + 2]
-                                                 (j + 3 * this->degree, dof)
-                                              * this->shape_value_component
-                                                (j + 3 * this->degree,
-                                                 quadrature_points[q_point],
-                                                 0);
-                            tmp (2 * i + 5) -= this->restriction[index][i + 2]
-                                               (i * this->degree + j, dof)
-                                            * this->shape_value_component
-                                              (i * this->degree + j,
-                                               quadrature_points[q_point], 1);
+                            quadrature_point (1) -= 1.0;
+                            k = 1;
                           }
-
-                      tmp *= quadrature.weight (q_point);
-
-                      for (unsigned int i = 0; i < this->degree; ++i)
+                      }
+                    
+                    else
+                      if (quadrature_points[q_point] (1) < 0.5)
                         {
-                          const double L_i_0
-                            = legendre_polynomials[i].value
-                              (quadrature_points[q_point] (0));
-                          const double L_i_1
-                            = legendre_polynomials[i].value
-                              (quadrature_points[q_point] (1));
-
-                          for (unsigned int j = 0; j < this->degree-1; ++j)
-                            {
-                              const double l_j_0
-                                = L_i_0 * lobatto_polynomials[j + 2].value
-                                          (quadrature_points[q_point] (1));
-                              const double l_j_1
-                                = L_i_1 * lobatto_polynomials[j + 2].value
-                                          (quadrature_points[q_point] (0));
-
-                              for (unsigned int k = 0; k < 4; ++k)
-                                {
-                                  system_rhs (i * (this->degree-1) + j, 2 * k)
-                                    += tmp (2 * k) * l_j_0;
-                                  system_rhs (i * (this->degree-1) + j, 2 * k + 1)
-                                    += tmp (2 * k + 1) * l_j_1;
-                                }
-                            }
+                          quadrature_point (0) -= 1.0;
+                          k = 2;
                         }
-                    }
-
-                  system_matrix_inv.mmult (solution, system_rhs);
-
-                  for (unsigned int i = 0; i < this->degree; ++i)
-                    for (unsigned int j = 0; j < this->degree-1; ++j)
-                      for (unsigned int k = 0; k < 4; ++k)
+                      
+                      else
                         {
-                          if (std::abs (solution (i * (this->degree-1) + j, 2 * k))
-                                > 1e-14)
-                            this->restriction[index][k]
-                              (i * (this->degree-1) + j + n_boundary_dofs, dof)
-                              = solution (i * (this->degree-1) + j, 2 * k);
-
-                          if (std::abs (solution (i * (this->degree-1) + j, 2 * k + 1))
-                                > 1e-14)
-                            this->restriction[index][k]
-                            (i + (this->degree-1 + j) * this->degree + n_boundary_dofs,
-                             dof)
-                              = solution (i * (this->degree-1) + j, 2 * k + 1);
+                          quadrature_point (0) -= 1.0;
+                          quadrature_point (1) -= 1.0;
+                          k = 3;
                         }
-                }
+                                  
+                    for (unsigned int i = 0; i < this->degree; ++i)
+                      {
+                        const double L_i_0 = weight
+                                             * legendre_polynomials[i].value
+                                               (quadrature_points[q_point] (0));
+                        const double L_i_1 = weight
+                                             * legendre_polynomials[i].value
+                                               (quadrature_points[q_point] (1));
+
+                        for (unsigned int j = 0; j < deg; ++j)
+                          {
+                            this->restriction[index][k] (i * deg + j + n_boundary_dofs, dof)
+                              += L_i_0 * (this->shape_grad_component (dof, quadrature_point, 0)[1]
+                                          * legendre_polynomials[j + 1].value
+                                            (quadrature_points[q_point] (1))
+                                          + this->shape_value_component (dof, quadrature_point, 0)
+                                          * lobatto_polynomials[j + 2].value
+                                            (quadrature_points[q_point] (1)));
+                            this->restriction[index][k] (i + (deg + j) * this->degree
+                                                           + n_boundary_dofs, dof)
+                              += L_i_1 * (this->shape_grad_component (dof, quadrature_point, 1)[0]
+                                          * legendre_polynomials[j + 1].value
+                                            (quadrature_points[q_point] (0))
+                                          + this->shape_value_component (dof, quadrature_point, 1)
+                                          * lobatto_polynomials[j + 2].value
+                                            (quadrature_points[q_point] (0)));
+                          }
+                      }
+                  }
             }
 
           break;
@@ -945,292 +682,96 @@ FE_Nedelec<dim>::initialize_restriction ()
         {
                         // First interpolate the shape
                         // functions of the child cells
-                        // to the lowest order shape
-                        // functions of the parent cell.
+                        // to the edge shape functions
+                        // of the parent cell.
           for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
             for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
                  ++q_point)
               {
                 const double weight = 2.0 * edge_quadrature.weight (q_point);
-
+                
                 if (edge_quadrature_points[q_point] (0) < 0.5)
-                  for (unsigned int i = 0; i < 2; ++i)
-                    for (unsigned int j = 0; j < 2; ++j)
-                      {
-                        Point<dim> quadrature_point (i,
-                                                     2.0 * edge_quadrature_points[q_point] (0),
-                                                     j);
-
-                        this->restriction[index][i + 4 * j]
-                        ((i + 4 * j) * this->degree, dof)
-                          += weight * this->shape_value_component (dof,
-                                                                   quadrature_point,
-                                                                   1);
-                        quadrature_point
-                          = Point<dim> (2.0 * edge_quadrature_points[q_point] (0),
-                                        i, j);
-                        this->restriction[index][2 * (i + 2 * j)]
-                        ((i + 4 * j + 2) * this->degree, dof)
-                          += weight * this->shape_value_component (dof,
-                                                                   quadrature_point,
-                                                                   0);
-                        quadrature_point = Point<dim> (i, j,
-                                                       2.0 * edge_quadrature_points[q_point] (0));
-                        this->restriction[index][i + 2 * j]
-                        ((i + 2 * (j + 4)) * this->degree, dof)
-                          += weight * this->shape_value_component (dof,
-                                                                   quadrature_point,
-                                                                   2);
-                      }
+                  for (unsigned int i = 0; i < this->degree; ++i)
+                    {
+                      const double L_i
+                        = weight * legendre_polynomials[i].value
+                                   (edge_quadrature_points[q_point] (0));
 
-                    else
-                      for (unsigned int i = 0; i < 2; ++i)
-                        for (unsigned int j = 0; j < 2; ++j)
+                      for (unsigned int j = 0; j < 2; ++j)
+                        for (unsigned int k = 0; k < 2; ++k)
+                          {
+                            Point<dim> quadrature_point (j,
+                                                         2.0 * edge_quadrature_points[q_point] (0),
+                                                         k);
+
+                            this->restriction[index][j + 4 * k]
+                            (i + (j + 4 * k) * this->degree, dof)
+                              += L_i * this->shape_value_component (dof,
+                                                                    quadrature_point,
+                                                                    1);
+                            quadrature_point
+                              = Point<dim> (2.0 * edge_quadrature_points[q_point] (0),
+                                            j, k);
+                            this->restriction[index][2 * (j + 2 * k)]
+                            (i + (j + 4 * k + 2) * this->degree, dof)
+                              += L_i * this->shape_value_component (dof,
+                                                                    quadrature_point,
+                                                                    0);
+                            quadrature_point = Point<dim> (j, k,
+                                                           2.0 * edge_quadrature_points[q_point] (0));
+                            this->restriction[index][j + 2 * k]
+                            (i + (j + 2 * (k + 4)) * this->degree, dof)
+                              += L_i * this->shape_value_component (dof,
+                                                                    quadrature_point,
+                                                                    2);
+                          }
+                    }
+                  
+                else
+                  for (unsigned int i = 0; i < this->degree; ++i)
+                    {
+                      const double L_i
+                        = weight * legendre_polynomials[i].value
+                                   (edge_quadrature_points[q_point] (0));
+
+                      for (unsigned int j = 0; j < 2; ++j)
+                        for (unsigned int k = 0; k < 2; ++k)
                           {
-                            Point<dim> quadrature_point (i,
+                            Point<dim> quadrature_point (j,
                                                          2.0 * edge_quadrature_points[q_point] (0)
-                                                             - 1.0, j);
+                                                             - 1.0, k);
 
-                            this->restriction[index][i + 4 * j + 2]
-                            ((i + 4 * j) * this->degree, dof)
-                              += weight * this->shape_value_component (dof,
-                                                                       quadrature_point,
-                                                                       1);
+                            this->restriction[index][j + 4 * k + 2]
+                            (i + (j + 4 * k) * this->degree, dof)
+                              += L_i * this->shape_value_component (dof,
+                                                                    quadrature_point,
+                                                                    1);
                             quadrature_point
-                              = Point<dim> (2.0 * edge_quadrature_points[q_point] (0)
-                                                - 1.0, i, j);
-                            this->restriction[index][2 * (i + 2 * j) + 1]
-                            ((i + 4 * j + 2) * this->degree, dof)
-                              += weight * this->shape_value_component (dof,
-                                                                       quadrature_point,
-                                                                       0);
-                            quadrature_point = Point<dim> (i, j,
+                              = Point<dim> (2.0 * edge_quadrature_points[q_point] (0) - 1.0,
+                                            j, k);
+                            this->restriction[index][2 * (j + 2 * k) + 1]
+                            (i + (j + 4 * k + 2) * this->degree, dof)
+                              += L_i * this->shape_value_component (dof,
+                                                                    quadrature_point,
+                                                                    0);
+                            quadrature_point = Point<dim> (j, k,
                                                            2.0 * edge_quadrature_points[q_point] (0)
                                                                - 1.0);
-                            this->restriction[index][i + 2 * (j + 2)]
-                            ((i + 2 * (j + 4)) * this->degree, dof)
-                              += weight * this->shape_value_component (dof,
-                                                                       quadrature_point,
-                                                                       2);
+                            this->restriction[index][j + 2 * (k + 2)]
+                            (i + (j + 2 * (k + 4)) * this->degree, dof)
+                              += L_i * this->shape_value_component (dof,
+                                                                    quadrature_point,
+                                                                    2);
                           }
+                    }
               }
 
-                        // Then project the shape functions
-                        // of the child cells to the higher
-                        // order shape functions of the
-                        // parent cell.
+                        // Then interpolate the shape functions
+                        // of the child cells to the face
+                        // and interior shape functions of
+                        // the parent cell.
           if (this->degree > 1)
             {
-              const unsigned int deg = this->degree-1;
-              const std::vector<Polynomials::Polynomial<double> >&
-                legendre_polynomials
-                  = Polynomials::Legendre::generate_complete_basis (deg);
-              FullMatrix<double> system_matrix_inv (deg, deg);
-
-              {
-                FullMatrix<double> assembling_matrix (deg,
-                                                      n_edge_quadrature_points);
-
-                for (unsigned int q_point = 0;
-                     q_point < n_edge_quadrature_points; ++q_point)
-                  {
-                    const double weight = std::sqrt (edge_quadrature.weight
-                                                     (q_point));
-
-                    for (unsigned int i = 0; i < deg; ++i)
-                      assembling_matrix (i, q_point) = weight
-                                                       * legendre_polynomials[i + 1].value
-                                                         (edge_quadrature_points[q_point] (0));
-                  }
-
-                FullMatrix<double> system_matrix (deg, deg);
-
-                assembling_matrix.mTmult (system_matrix, assembling_matrix);
-                system_matrix_inv.invert (system_matrix);
-              }
-
-              FullMatrix<double> solution (deg, 6);
-              FullMatrix<double> system_rhs (deg, 6);
-              Vector<double> tmp (6);
-
-              for (unsigned int i = 0; i < 2; ++i)
-                for (unsigned int j = 0; j < 2; ++j)
-                  for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
-                    {
-                      system_rhs = 0.0;
-
-                      for (unsigned int q_point = 0;
-                           q_point < n_edge_quadrature_points; ++q_point)
-                        {
-                          const double weight = edge_quadrature.weight
-                                                (q_point);
-                          const Point<dim> quadrature_point_0 (i,
-                                                               edge_quadrature_points[q_point] (0),
-                                                               j);
-                          const Point<dim>
-                            quadrature_point_1
-                            (edge_quadrature_points[q_point] (0), i, j);
-                          const Point<dim> quadrature_point_2 (i, j,
-                                                               edge_quadrature_points[q_point] (0));
-
-                          if (edge_quadrature_points[q_point] (0) < 0.5)
-                            {
-                              Point<dim> quadrature_point_3 (i,
-                                                             2.0 * edge_quadrature_points[q_point] (0),
-                                                             j);
-
-                              tmp (0) = weight
-                                        * (2.0 * this->shape_value_component
-                                                 (dof, quadrature_point_3, 1)
-                                               - this->restriction[index][i + 4 * j]
-                                                 ((i + 4 * j) * this->degree,
-                                                  dof)
-                                               * this->shape_value_component
-                                                 ((i + 4 * j) * this->degree,
-                                                  quadrature_point_0, 1));
-                              tmp (1) = -1.0 * weight
-                                             * this->restriction[index][i + 4 * j + 2]
-                                               ((i + 4 * j) * this->degree,
-                                                dof)
-                                             * this->shape_value_component
-                                               ((i + 4 * j) * this->degree,
-                                                quadrature_point_0, 1);
-                              quadrature_point_3
-                                = Point<dim> (2.0 * edge_quadrature_points[q_point] (0),
-                                              i, j);
-                              tmp (2) = weight
-                                        * (2.0 * this->shape_value_component
-                                                 (dof, quadrature_point_3, 0)
-                                               - this->restriction[index][2 * (i + 2 * j)]
-                                                 ((i + 4 * j + 2) * this->degree,
-                                                  dof)
-                                               * this->shape_value_component
-                                                 ((i + 4 * j + 2) * this->degree,
-                                                  quadrature_point_1, 0));
-                              tmp (3) = -1.0 * weight
-                                             * this->restriction[index][2 * (i + 2 * j) + 1]
-                                               ((i + 4 * j + 2) * this->degree,
-                                                dof)
-                                             * this->shape_value_component
-                                               ((i + 4 * j + 2) * this->degree,
-                                                quadrature_point_1, 0);
-                              quadrature_point_3 = Point<dim> (i, j,
-                                                               2.0 * edge_quadrature_points[q_point] (0));
-                              tmp (4) = weight
-                                        * (2.0 * this->shape_value_component
-                                                 (dof, quadrature_point_3, 2)
-                                               - this->restriction[index][i + 2 * j]
-                                                 ((i + 2 * (j + 4)) * this->degree,
-                                                  dof)
-                                               * this->shape_value_component
-                                                 ((i + 2 * (j + 4)) * this->degree,
-                                                  quadrature_point_2, 2));
-                              tmp (5) = -1.0 * weight
-                                             * this->restriction[index][i + 2 * (j + 2)]
-                                               ((i + 2 * (j + 4)) * this->degree,
-                                                dof)
-                                             * this->shape_value_component
-                                               ((i + 2 * (j + 4)) * this->degree,
-                                                quadrature_point_2, 2);
-                            }
-
-                          else
-                            {
-                              tmp (0) = -1.0 * weight
-                                             * this->restriction[index][i + 4 * j]
-                                               ((i + 4 * j) * this->degree,
-                                                dof)
-                                             * this->shape_value_component
-                                               ((i + 4 * j) * this->degree,
-                                                quadrature_point_0, 1);
-
-                              Point<dim> quadrature_point_3 (i,
-                                                             2.0 * edge_quadrature_points[q_point] (0)
-                                                                 - 1.0, j);
-
-                              tmp (1) = weight
-                                        * (2.0 * this->shape_value_component
-                                                 (dof, quadrature_point_3, 1)
-                                               - this->restriction[index][i + 4 * j + 2]
-                                                 ((i + 4 * j) * this->degree,
-                                                  dof)
-                                               * this->shape_value_component
-                                                 ((i + 4 * j) * this->degree,
-                                                  quadrature_point_0, 1));
-                              tmp (2) = -1.0 * weight
-                                             * this->restriction[index][2 * (i + 2 * j)]
-                                               ((i + 4 * j + 2) * this->degree,
-                                                dof)
-                                             * this->shape_value_component
-                                               ((i + 4 * j + 2) * this->degree,
-                                                quadrature_point_1, 0);
-                              quadrature_point_3
-                                = Point<dim> (2.0 * edge_quadrature_points[q_point] (0)
-                                                  - 1.0, i, j);
-                              tmp (3) = weight
-                                        * (2.0 * this->shape_value_component
-                                                 (dof, quadrature_point_3, 0)
-                                               - this->restriction[index][2 * (i + 2 * j) + 1]
-                                                 ((i + 4 * j + 2) * this->degree,
-                                                  dof)
-                                               * this->shape_value_component
-                                                 ((i + 4 * j + 2) * this->degree,
-                                                  quadrature_point_1, 0));
-                              tmp (4) = -1.0 * weight
-                                             * this->restriction[index][i + 2 * j]
-                                               ((i + 2 * (j + 4)) * this->degree,
-                                                dof)
-                                             * this->shape_value_component
-                                               ((i + 2 * (j + 4)) * this->degree,
-                                                quadrature_point_2, 2);
-                              quadrature_point_3 = Point<dim> (i, j,
-                                                               2.0 * edge_quadrature_points[q_point] (0)
-                                                                   - 1.0);
-                              tmp (5) = weight
-                                        * (2.0 * this->shape_value_component
-                                                 (dof, quadrature_point_3, 2)
-                                               - this->restriction[index][i + 2 * (j + 2)]
-                                                 ((i + 2 * (j + 4)) * this->degree,
-                                                  dof)
-                                               * this->shape_value_component
-                                                 ((i + 2 * (j + 4)) * this->degree,
-                                                  quadrature_point_2, 2));
-                            }
-
-                          for (unsigned int k = 0; k < deg; ++k)
-                            {
-                              const double L_k
-                                = legendre_polynomials[k + 1].value
-                                  (edge_quadrature_points[q_point] (0));
-
-                              for (unsigned int l = 0; l < tmp.size (); ++l)
-                                system_rhs (k, l) += tmp (l) * L_k;
-                            }
-                        }
-
-                      system_matrix_inv.mmult (solution, system_rhs);
-
-                      for (unsigned int k = 0; k < 2; ++k)
-                        for (unsigned int l = 0; l < deg; ++l)
-                          {
-                            if (std::abs (solution (l, k)) > 1e-14)
-                              this->restriction[index][i + 2 * (2 * j + k)]
-                              ((i + 4 * j) * this->degree + l + 1, dof)
-                                = solution (l, k);
-
-                            if (std::abs (solution (l, k + 2)) > 1e-14)
-                              this->restriction[index][2 * (i + 2 * j) + k]
-                              ((i + 4 * j + 2) * this->degree + l + 1, dof)
-                                = solution (l, k + 2);
-
-                            if (std::abs (solution (l, k + 4)) > 1e-14)
-                              this->restriction[index][i + 2 * (j + 2 * k)]
-                              ((i + 2 * (j + 4)) * this->degree + l + 1,
-                               dof)
-                                = solution (l, k + 4);
-                          }
-                    }
-
               const QGauss<2> face_quadrature (2 * this->degree);
               const std::vector<Point<2> >& face_quadrature_points
                 = face_quadrature.get_points ();
@@ -1238,724 +779,268 @@ FE_Nedelec<dim>::initialize_restriction ()
                 lobatto_polynomials
                   = Polynomials::Lobatto::generate_complete_basis
                     (this->degree);
+              const QGauss<dim> quadrature (2 * this->degree);
+              const std::vector<Point<dim> >&
+                quadrature_points = quadrature.get_points ();
               const unsigned int n_edge_dofs
                 = GeometryInfo<dim>::lines_per_cell * this->degree;
+              const unsigned int n_boundary_dofs
+                = 2 * GeometryInfo<dim>::faces_per_cell * deg * this->degree
+                    + n_edge_dofs;
               const unsigned int& n_face_quadrature_points
                 = face_quadrature.size ();
-
-              {
-                FullMatrix<double> assembling_matrix
-                                   (deg * this->degree,
-                                    n_face_quadrature_points);
-
-                for (unsigned int q_point = 0;
-                     q_point < n_face_quadrature_points; ++q_point)
-                  {
-                    const double weight
-                      = std::sqrt (face_quadrature.weight (q_point));
-
-                    for (unsigned int i = 0; i <= deg; ++i)
-                      {
-                        const double L_i = weight
-                                           * legendre_polynomials[i].value
-                                             (face_quadrature_points[q_point] (0));
-
-                        for (unsigned int j = 0; j < deg; ++j)
-                          assembling_matrix (i * deg + j, q_point)
-                            = L_i * lobatto_polynomials[j + 2].value
-                                    (face_quadrature_points[q_point] (1));
-                      }
-                  }
-
-                FullMatrix<double> system_matrix (assembling_matrix.m (),
-                                                  assembling_matrix.m ());
-
-                assembling_matrix.mTmult (system_matrix,
-                                          assembling_matrix);
-                system_matrix_inv.reinit (system_matrix.m (),
-                                          system_matrix.m ());
-                system_matrix_inv.invert (system_matrix);
-              }
-
-              solution.reinit (system_matrix_inv.m (), 24);
-              system_rhs.reinit (system_matrix_inv.m (), 24);
-              tmp.reinit (24);
-
-              for (unsigned int i = 0; i < 2; ++i)
-                for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
-                  {
-                    system_rhs = 0.0;
-
-                    for (unsigned int q_point = 0;
-                         q_point < n_face_quadrature_points; ++q_point)
-                      {
-                        tmp = 0.0;
-
-                        if (face_quadrature_points[q_point] (0) < 0.5)
-                          {
-                            if (face_quadrature_points[q_point] (1) < 0.5)
-                              {
-                                Point<dim> quadrature_point_0 (i,
-                                                                   2.0 * face_quadrature_points[q_point] (0),
-                                                                   2.0 * face_quadrature_points[q_point] (1));
-
-                                tmp (0) += 2.0 * this->shape_value_component
-                                                 (dof, quadrature_point_0, 1);
-                                tmp (1) += 2.0 * this->shape_value_component
-                                                 (dof, quadrature_point_0, 2);
-                                quadrature_point_0
-                                  = Point<dim> (2.0 * face_quadrature_points[q_point] (0),
-                                                i,
-                                                2.0 * face_quadrature_points[q_point] (1));
-                                tmp (8) += 2.0 * this->shape_value_component
-                                                 (dof, quadrature_point_0, 2);
-                                tmp (9) += 2.0 * this->shape_value_component
-                                                 (dof, quadrature_point_0, 0);
-                                quadrature_point_0
-                                  = Point<dim> (2.0 * face_quadrature_points[q_point] (0),
-                                                2.0 * face_quadrature_points[q_point] (1),
-                                                i);
-                                tmp (16) += 2.0 * this->shape_value_component
-                                                  (dof, quadrature_point_0, 0);
-                                tmp (17) += 2.0 * this->shape_value_component
-                                                  (dof, quadrature_point_0, 1);
-                              }
-
-                            else
-                              {
-                                Point<dim> quadrature_point_0 (i,
-                                                                   2.0 * face_quadrature_points[q_point] (0),
-                                                                   2.0 * face_quadrature_points[q_point] (1)
-                                                                       - 1.0);
-
-                                tmp (2) += 2.0 * this->shape_value_component
-                                                 (dof, quadrature_point_0, 1);
-                                tmp (3) += 2.0 * this->shape_value_component
-                                                 (dof, quadrature_point_0, 2);
-                                quadrature_point_0
-                                  = Point<dim> (2.0 * face_quadrature_points[q_point] (0),
-                                                i,
-                                                2.0 * face_quadrature_points[q_point] (1)
-                                                    - 1.0);
-                                tmp (10) += 2.0 * this->shape_value_component
-                                                  (dof, quadrature_point_0, 2);
-                                tmp (11) += 2.0 * this->shape_value_component
-                                                  (dof, quadrature_point_0, 0);
-                                quadrature_point_0
-                                  = Point<dim> (2.0 * face_quadrature_points[q_point] (0),
-                                                2.0 * face_quadrature_points[q_point] (1)
-                                                    - 1.0, i);
-                                tmp (18) += 2.0 * this->shape_value_component
-                                                  (dof, quadrature_point_0, 0);
-                                tmp (19) += 2.0 * this->shape_value_component
-                                                  (dof, quadrature_point_0, 1);
-                              }
-                          }
-
-                        else
-                          if (face_quadrature_points[q_point] (1) < 0.5)
+              const unsigned int& n_quadrature_points = quadrature.size ();
+              
+                            // First, the interpolation to
+                            // the face shape functions.
+              for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
+                {
+                  for (unsigned int q_point = 0; q_point < n_face_quadrature_points;
+                       ++q_point)
+                    {
+                      const double weight = 2.0 * face_quadrature.weight (q_point);
+                      
+                      for (unsigned int i = 0; i < 2; ++i)
+                        {
+                          Point<dim> quadrature_point (i,
+                                                       2.0 * face_quadrature_points[q_point] (0),
+                                                       2.0 * face_quadrature_points[q_point] (1));
+                          unsigned int l;
+                          unsigned int m;
+                          
+                          if (face_quadrature_points[q_point] (0) < 0.5)
                             {
-                              Point<dim> quadrature_point_0 (i,
-                                                             2.0 * face_quadrature_points[q_point] (0)
-                                                                 - 1.0,
-                                                             2.0 * face_quadrature_points[q_point] (1));
-
-                              tmp (4) += 2.0 * this->shape_value_component
-                                               (dof, quadrature_point_0, 1);
-                              tmp (5) += 2.0 * this->shape_value_component
-                                               (dof, quadrature_point_0, 2);
-                              quadrature_point_0
-                                = Point<dim> (2.0 * face_quadrature_points[q_point] (0)
-                                                  - 1.0, i,
-                                              2.0 * face_quadrature_points[q_point] (1));
-                              tmp (12) += 2.0 * this->shape_value_component
-                                                (dof, quadrature_point_0, 2);
-                              tmp (13) += 2.0 * this->shape_value_component
-                                                (dof, quadrature_point_0, 0);
-                              quadrature_point_0
-                                = Point<dim> (2.0 * face_quadrature_points[q_point] (0)
-                                                  - 1.0,
-                                              2.0 * face_quadrature_points[q_point] (1),
-                                              i);
-                              tmp (20) += 2.0 * this->shape_value_component
-                                                (dof, quadrature_point_0, 0);
-                              tmp (21) += 2.0 * this->shape_value_component
-                                                (dof, quadrature_point_0, 1);
+                              m = 0;
+                              
+                              if (face_quadrature_points[q_point] (1) < 0.5)
+                                l = 0;
+                              
+                              else
+                                {
+                                  quadrature_point (2) -= 1.0;
+                                  l = 1;
+                                }
                             }
-
+                        
                           else
                             {
-                              Point<dim> quadrature_point_0 (i,
-                                                             2.0 * face_quadrature_points[q_point] (0)
-                                                                 - 1.0,
-                                                             2.0 * face_quadrature_points[q_point] (1)
-                                                                 - 1.0);
-
-                              tmp (6) += 2.0 * this->shape_value_component
-                                               (dof, quadrature_point_0, 1);
-                              tmp (7) += 2.0 * this->shape_value_component
-                                               (dof, quadrature_point_0, 2);
-                              quadrature_point_0
-                                = Point<dim> (2.0 * face_quadrature_points[q_point] (0)
-                                                  - 1.0, i,
-                                              2.0 * face_quadrature_points[q_point] (1)
-                                                  - 1.0);
-                              tmp (14) += 2.0 * this->shape_value_component
-                                                (dof, quadrature_point_0, 2);
-                              tmp (15) += 2.0 * this->shape_value_component
-                                                (dof, quadrature_point_0, 0);
-                              quadrature_point_0
-                                = Point<dim> (2.0 * face_quadrature_points[q_point] (0)
-                                                  - 1.0,
-                                              2.0 * face_quadrature_points[q_point] (1)
-                                                  - 1.0, i);
-                              tmp (22) += 2.0 * this->shape_value_component
-                                                (dof, quadrature_point_0, 0);
-                              tmp (23) += 2.0 * this->shape_value_component
-                                                (dof, quadrature_point_0, 1);
+                              quadrature_point (1) -= 1.0;
+                              m = 1;
+                              
+                              if (face_quadrature_points[q_point] (1) < 0.5)
+                                l = 0;
+                              
+                              else
+                                {
+                                  quadrature_point (2) -= 1.0;
+                                  l = 1;
+                                }
                             }
-
-                        const Point<dim> quadrature_point_0 (i,
-                                                             face_quadrature_points[q_point] (0),
-                                                             face_quadrature_points[q_point] (1));
-                        const Point<dim> quadrature_point_1
-                                         (face_quadrature_points[q_point] (0),
-                                          i,
-                                          face_quadrature_points[q_point] (1));
-                        const Point<dim> quadrature_point_2
-                                         (face_quadrature_points[q_point] (0),
-                                          face_quadrature_points[q_point] (1),
-                                          i);
-
-                        for (unsigned int j = 0; j < 2; ++j)
-                          for (unsigned int k = 0; k < 2; ++k)
-                            for (unsigned int l = 0; l <= deg; ++l)
-                              {
-                                tmp (2 * (j + 2 * k))
-                                  -= this->restriction[index][i + 2 * (2 * j + k)]
-                                     ((i + 4 * j) * this->degree + l, dof)
-                                  * this->shape_value_component
-                                    ((i + 4 * j) * this->degree + l,
-                                     quadrature_point_0, 1);
-                                tmp (2 * (j + 2 * k) + 1)
-                                  -= this->restriction[index][i + 2 * (2 * j + k)]
-                                     ((i + 2 * (k + 4)) * this->degree + l,
-                                      dof)
-                                  * this->shape_value_component
-                                    ((i + 2 * (k + 4)) * this->degree + l,
-                                     quadrature_point_0, 2);
-                                tmp (2 * (j + 2 * (k + 2)))
-                                  -= this->restriction[index][2 * (i + 2 * j) + k]
-                                     ((2 * (i + 4) + k) * this->degree + l,
-                                      dof)
-                                  * this->shape_value_component
-                                    ((2 * (i + 4) + k) * this->degree + l,
-                                     quadrature_point_1, 2);
-                                tmp (2 * (j + 2 * k) + 9)
-                                  -= this->restriction[index][2 * (i + 2 * j) + k]
-                                     ((i + 4 * j + 2) * this->degree + l,
-                                      dof)
-                                  * this->shape_value_component
-                                    ((i + 4 * j + 2) * this->degree + l,
-                                     quadrature_point_1, 0);
-                                tmp (2 * (j + 2 * (k + 4)))
-                                  -= this->restriction[index][2 * (2 * i + j) + k]
-                                     ((4 * i + j + 2) * this->degree + l,
-                                      dof)
-                                  * this->shape_value_component
-                                    ((4 * i + j + 2) * this->degree + l,
-                                     quadrature_point_2, 0);
-                                tmp (2 * (j + 2 * k) + 17)
-                                  -= this->restriction[index][2 * (2 * i + j) + k]
-                                     ((4 * i + k) * this->degree + l, dof)
-                                  * this->shape_value_component
-                                    ((4 * i + k) * this->degree + l,
-                                     quadrature_point_2, 1);
-                              }
-
-                        tmp *= face_quadrature.weight (q_point);
-
-                        for (unsigned int j = 0; j <= deg; ++j)
-                          {
-                            const double L_j_0
-                              = legendre_polynomials[j].value
-                                (face_quadrature_points[q_point] (0));
-                            const double L_j_1
-                              = legendre_polynomials[j].value
-                                (face_quadrature_points[q_point] (1));
-
-                            for (unsigned int k = 0; k < deg; ++k)
-                              {
-                                const double l_k_0
-                                  = L_j_0 * lobatto_polynomials[k + 2].value
-                                            (face_quadrature_points[q_point] (1));
-                                const double l_k_1
-                                  = L_j_1 * lobatto_polynomials[k + 2].value
-                                            (face_quadrature_points[q_point] (0));
-
-                                for (unsigned int l = 0; l < 4; ++l)
-                                  {
-                                    system_rhs (j * deg + k, 2 * l)
-                                      += tmp (2 * l) * l_k_0;
-                                    system_rhs (j * deg + k, 2 * l + 1)
-                                      += tmp (2 * l + 1) * l_k_1;
-                                    system_rhs (j * deg + k, 2 * (l + 4))
-                                      += tmp (2 * (l + 4)) * l_k_1;
-                                    system_rhs (j * deg + k, 2 * l + 9)
-                                      += tmp (2 * l + 9) * l_k_0;
-                                    system_rhs (j * deg + k, 2 * (l + 8))
-                                      += tmp (2 * (l + 8)) * l_k_0;
-                                    system_rhs (j * deg + k, 2 * l + 17)
-                                      += tmp (2 * l + 17) * l_k_1;
-                                  }
-                              }
-                          }
-                      }
-
-                    system_matrix_inv.mmult (solution, system_rhs);
-
-                    for (unsigned int j = 0; j < 2; ++j)
-                      for (unsigned int k = 0; k < 2; ++k)
-                        for (unsigned int l = 0; l <= deg; ++l)
-                          for (unsigned int m = 0; m < deg; ++m)
+                          
+                          for (unsigned int j = 0; j < this->degree; ++j)
                             {
-                              if (std::abs (solution (l * deg + m,
-                                                      2 * (j + 2 * k)))
-                                    > 1e-14)
-                                this->restriction[index][i + 2 * (2 * j + k)]
-                                ((2 * i * this->degree + l) * deg + m
-                                                            + n_edge_dofs,
-                                 dof) = solution (l * deg + m,
-                                                  2 * (j + 2 * k));
-
-                               if (std::abs (solution (l * deg + m,
-                                                       2 * (j + 2 * k) + 1))
-                                     > 1e-14)
-                                 this->restriction[index][i + 2 * (2 * j + k)]
-                                 (((2 * i + 1) * deg + m) * this->degree + l
-                                                          + n_edge_dofs, dof)
-                                   = solution (l * deg + m,
-                                               2 * (j + 2 * k) + 1);
-
-                               if (std::abs (solution (l * deg + m,
-                                                       2 * (j + 2 * (k + 2))))
-                                     > 1e-14)
-                                 this->restriction[index][2 * (i + 2 * j) + k]
-                                 ((2 * (i + 2) * this->degree + l) * deg + m
-                                                                   + n_edge_dofs,
-                                  dof) = solution (l * deg + m,
-                                                   2 * (j + 2 * (k + 2)));
-
-                               if (std::abs (solution (l * deg + m,
-                                                       2 * (j + 2 * k) + 9))
-                                     > 1e-14)
-                                 this->restriction[index][2 * (i + 2 * j) + k]
-                                 (((2 * i + 5) * deg + m) * this->degree + l
-                                                          + n_edge_dofs, dof)
-                                   = solution (l * deg + m,
-                                               2 * (j + 2 * k) + 9);
-
-                               if (std::abs (solution (l * deg + m,
-                                                       2 * (j + 2 * (k + 4))))
-                                     > 1e-14)
-                                 this->restriction[index][2 * (2 * i + j) + k]
-                                 ((2 * (i + 4) * this->degree + l) * deg + m
-                                                                   + n_edge_dofs,
-                                  dof) = solution (l * deg + m,
-                                                   2 * (j + 2 * (k + 4)));
-
-                               if (std::abs (solution (l * deg + m,
-                                                       2 * (j + 2 * k) + 17))
-                                     > 1e-14)
-                                 this->restriction[index][2 * (2 * i + j) + k]
-                                 (((2 * i + 9) * deg + m) * this->degree + l
-                                                          + n_edge_dofs, dof)
-                                   = solution (l * deg + m,
-                                               2 * (j + 2 * k) + 17);
-                             }
-                  }
-
-              const QGauss<dim> quadrature (2 * this->degree);
-              const std::vector<Point<dim> >&
-                quadrature_points = quadrature.get_points ();
-              const unsigned int n_boundary_dofs
-                = 2 * GeometryInfo<dim>::faces_per_cell * deg * this->degree
-                    + n_edge_dofs;
-              const unsigned int& n_quadrature_points = quadrature.size ();
-
-              {
-                FullMatrix<double>
-                  assembling_matrix (deg * deg * this->degree,
-                                     n_quadrature_points);
-
-                for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                     ++q_point)
-                  {
-                    const double weight = std::sqrt (quadrature.weight
-                                                     (q_point));
-
-                    for (unsigned int i = 0; i <= deg; ++i)
-                      {
-                        const double L_i = weight
-                                           * legendre_polynomials[i].value
-                                             (quadrature_points[q_point] (0));
-
-                        for (unsigned int j = 0; j < deg; ++j)
-                          {
-                            const double l_j
-                              = L_i * lobatto_polynomials[j + 2].value
-                                      (quadrature_points[q_point] (1));
-
-                            for (unsigned int k = 0; k < deg; ++k)
-                              assembling_matrix ((i * deg + j) * deg + k,
-                                                 q_point)
-                                = l_j * lobatto_polynomials[k + 2].value
-                                        (quadrature_points[q_point] (2));
-                          }
-                      }
-                  }
-
-                FullMatrix<double> system_matrix (assembling_matrix.m (),
-                                                  assembling_matrix.m ());
-
-                assembling_matrix.mTmult (system_matrix, assembling_matrix);
-                system_matrix_inv.reinit (system_matrix.m (),
-                                          system_matrix.m ());
-                system_matrix_inv.invert (system_matrix);
-              }
-
-              solution.reinit (system_matrix_inv.m (), 24);
-              system_rhs.reinit (system_matrix_inv.m (), 24);
-              tmp.reinit (24);
-
-              for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
-                {
-                  system_rhs = 0.0;
-
-                  for (unsigned int q_point = 0;
-                       q_point < n_quadrature_points; ++q_point)
+                              const double L_j_0
+                                = weight * legendre_polynomials[j].value
+                                           (face_quadrature_points[q_point] (0));
+                              const double L_j_1
+                                = weight * legendre_polynomials[j].value
+                                           (face_quadrature_points[q_point] (1));
+                            
+                              for (unsigned int k = 0; k < deg; ++k)
+                                {
+                                  const double Le_k_0
+                                    = legendre_polynomials[k + 1].value
+                                      (face_quadrature_points[q_point] (0));
+                                  const double Le_k_1
+                                    = legendre_polynomials[k + 1].value
+                                      (face_quadrature_points[q_point] (1));
+                                  const double lo_k_0
+                                    = lobatto_polynomials[k + 2].value
+                                      (face_quadrature_points[q_point] (0));
+                                  const double lo_k_1
+                                    = lobatto_polynomials[k + 2].value
+                                      (face_quadrature_points[q_point] (1));
+                                
+                                  this->restriction[index][i + 2 * (2 * l + m)]
+                                  ((2 * i * this->degree + j) * deg + k + n_edge_dofs, dof)
+                                    += L_j_0 * (this->shape_grad_component (dof, quadrature_point, 1)[2]
+                                                * Le_k_1
+                                                + this->shape_value_component (dof, quadrature_point, 1)
+                                                * lo_k_1);
+                                  this->restriction[index][i + 2 * (2 * l + m)]
+                                  (((2 * i + 1) * deg + k) * this->degree + j + n_edge_dofs, dof)
+                                    += L_j_1 * (this->shape_grad_component (dof, quadrature_point, 2)[1]
+                                                * Le_k_0
+                                                + this->shape_value_component (dof, quadrature_point, 2)
+                                                * lo_k_0);
+                                  this->restriction[index][2 * (i + 2 * l) + m]
+                                  ((2 * (i + 2) * this->degree + j) * deg + k + n_edge_dofs, dof)
+                                    += L_j_0 * (this->shape_grad_component (dof, quadrature_point, 2)[0]
+                                                * Le_k_1
+                                                + this->shape_value_component (dof, quadrature_point, 2)
+                                                * lo_k_1);
+                                  this->restriction[index][2 * (i + 2 * l) + m]
+                                  (((2 * i + 5) * deg + k) * this->degree + j + n_edge_dofs, dof)
+                                    += L_j_1 * (this->shape_grad_component (dof, quadrature_point, 0)[2]
+                                                * Le_k_0
+                                                + this->shape_value_component (dof, quadrature_point, 0)
+                                                * lo_k_0);
+                                  this->restriction[index][2 * (2 * i + l) + m]
+                                  ((2 * (i + 4) * this->degree + j) * deg + k + n_edge_dofs, dof)
+                                    += L_j_0 * (this->shape_grad_component (dof, quadrature_point, 0)[1]
+                                                * Le_k_1
+                                                + this->shape_value_component (dof, quadrature_point, 0)
+                                                * lo_k_1);
+                                  this->restriction[index][2 * (2 * i + l) + m]
+                                  (((2 * i + 9) * deg + k) * this->degree + j + n_edge_dofs, dof)
+                                    += L_j_1 * (this->shape_grad_component (dof, quadrature_point, 1)[0]
+                                                * Le_k_0
+                                                + this->shape_value_component (dof, quadrature_point, 1)
+                                                * lo_k_0);
+                                }
+                            }
+                        }
+                    }
+                  
+                            // Then, the interpolation to
+                            // the interior shape functions.
+                  for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                       ++q_point)
                     {
-                      tmp = 0.0;
-
+                      Point<dim> quadrature_point = 2.0 * quadrature_points[q_point];
+                      unsigned int l;
+                      unsigned int m;
+                      unsigned int n;
+                      
                       if (quadrature_points[q_point] (0) < 0.5)
                         {
+                          n = 0;
+                          
                           if (quadrature_points[q_point] (1) < 0.5)
                             {
+                              m = 0;
+                              
                               if (quadrature_points[q_point] (2) < 0.5)
+                                l = 0;
+                              
+                              else
                                 {
-                                  const Point<dim> quadrature_point
-                                                   (2.0 * quadrature_points[q_point] (0),
-                                                    2.0 * quadrature_points[q_point] (1),
-                                                    2.0 * quadrature_points[q_point] (2));
-
-                                  tmp (0) += 2.0 * this->shape_value_component
-                                                   (dof, quadrature_point, 0);
-                                  tmp (1) += 2.0 * this->shape_value_component
-                                                   (dof, quadrature_point, 1);
-                                  tmp (2) += 2.0 * this->shape_value_component
-                                                   (dof, quadrature_point, 2);
+                                  quadrature_point (2) -= 1.0;
+                                  l = 1;
                                 }
-
+                            }
+                          
+                          else
+                            {
+                              quadrature_point (1) -= 1.0;
+                              m = 1;
+                              
+                              if (quadrature_points[q_point] (2) < 0.5)
+                                l = 0;
+                              
                               else
                                 {
-                                  const Point<dim> quadrature_point
-                                                   (2.0 * quadrature_points[q_point] (0),
-                                                    2.0 * quadrature_points[q_point] (1),
-                                                    2.0 * quadrature_points[q_point] (2)
-                                                        - 1.0);
-
-                                  tmp (3) += 2.0 * this->shape_value_component
-                                                   (dof, quadrature_point, 0);
-                                  tmp (4) += 2.0 * this->shape_value_component
-                                                   (dof, quadrature_point, 1);
-                                  tmp (5) += 2.0 * this->shape_value_component
-                                                   (dof, quadrature_point, 2);
+                                  quadrature_point (2) -= 1.0;
+                                  l = 1;
                                 }
                             }
-
-                          else
-                            if (quadrature_points[q_point] (2) < 0.5)
-                              {
-                                const Point<dim> quadrature_point
-                                                 (2.0 * quadrature_points[q_point] (0),
-                                                  2.0 * quadrature_points[q_point] (1)
-                                                      - 1.0,
-                                                  2.0 * quadrature_points[q_point] (2));
-
-                                tmp (6) += 2.0 * this->shape_value_component
-                                                 (dof, quadrature_point, 0);
-                                tmp (7) += 2.0 * this->shape_value_component
-                                                 (dof, quadrature_point, 1);
-                                tmp (8) += 2.0 * this->shape_value_component
-                                                 (dof, quadrature_point, 2);
-                              }
-
-                            else
-                              {
-                                const Point<dim> quadrature_point
-                                                 (2.0 * quadrature_points[q_point] (0),
-                                                  2.0 * quadrature_points[q_point] (1)
-                                                      - 1.0,
-                                                  2.0 * quadrature_points[q_point] (2)
-                                                      - 1.0);
-
-                                tmp (9) += 2.0 * this->shape_value_component
-                                                 (dof, quadrature_point, 0);
-                                tmp (10) += 2.0 * this->shape_value_component
-                                                  (dof, quadrature_point, 1);
-                                tmp (11) += 2.0 * this->shape_value_component
-                                                  (dof, quadrature_point, 2);
-                              }
                         }
-
+                      
                       else
-                        if (quadrature_points[q_point] (1) < 0.5)
-                          {
-                            if (quadrature_points[q_point] (2) < 0.5)
-                              {
-                                const Point<dim> quadrature_point
-                                                 (2.0 * quadrature_points[q_point] (0)
-                                                      - 1.0,
-                                                  2.0 * quadrature_points[q_point] (1),
-                                                  2.0 * quadrature_points[q_point] (2));
-
-                                tmp (12) += 2.0 * this->shape_value_component
-                                                  (dof, quadrature_point, 0);
-                                tmp (13) += 2.0 * this->shape_value_component
-                                                  (dof, quadrature_point, 1);
-                                tmp (14) += 2.0 * this->shape_value_component
-                                                  (dof, quadrature_point, 2);
-                              }
-
-                            else
-                              {
-                                const Point<dim> quadrature_point
-                                                 (2.0 * quadrature_points[q_point] (0)
-                                                      - 1.0,
-                                                  2.0 * quadrature_points[q_point] (1),
-                                                  2.0 * quadrature_points[q_point] (2)
-                                                      - 1.0);
-
-                                tmp (15) += 2.0 * this->shape_value_component
-                                                  (dof, quadrature_point, 0);
-                                tmp (16) += 2.0 * this->shape_value_component
-                                                  (dof, quadrature_point, 1);
-                                tmp (17) += 2.0 * this->shape_value_component
-                                                  (dof, quadrature_point, 2);
-                              }
-                          }
-
-                        else
-                          if (quadrature_points[q_point] (2) < 0.5)
+                        {
+                          quadrature_point (0) -= 1.0;
+                          n = 1;
+                          
+                          if (quadrature_points[q_point] (1) < 0.5)
                             {
-                              const Point<dim> quadrature_point
-                                               (2.0 * quadrature_points[q_point] (0)
-                                                    - 1.0,
-                                                2.0 * quadrature_points[q_point] (1)
-                                                    - 1.0,
-                                                2.0 * quadrature_points[q_point] (2));
-
-                              tmp (18) += 2.0 * this->shape_value_component
-                                                (dof, quadrature_point, 0);
-                              tmp (19) += 2.0 * this->shape_value_component
-                                                (dof, quadrature_point, 1);
-                              tmp (20) += 2.0 * this->shape_value_component
-                                                (dof, quadrature_point, 2);
+                              m = 0;
+                              
+                              if (quadrature_points[q_point] (2) < 0.5)
+                                l = 0;
+                              
+                              else
+                                {
+                                  quadrature_point (2) -= 1.0;
+                                  l = 1;
+                                }
                             }
-
+                          
                           else
                             {
-                              const Point<dim> quadrature_point
-                                               (2.0 * quadrature_points[q_point] (0)
-                                                    - 1.0,
-                                                2.0 * quadrature_points[q_point] (1)
-                                                    - 1.0,
-                                                2.0 * quadrature_points[q_point] (2)
-                                                    - 1.0);
-
-                              tmp (21) += 2.0 * this->shape_value_component
-                                                (dof, quadrature_point, 0);
-                              tmp (22) += 2.0 * this->shape_value_component
-                                                (dof, quadrature_point, 1);
-                              tmp (23) += 2.0 * this->shape_value_component
-                                                (dof, quadrature_point, 2);
+                              quadrature_point (1) -= 1.0;
+                              m = 1;
+                              
+                              if (quadrature_points[q_point] (2) < 0.5)
+                                l = 0;
+                              
+                              else
+                                {
+                                  quadrature_point (2) -= 1.0;
+                                  l = 1;
+                                }
                             }
-
-                      for (unsigned int i = 0; i < 2; ++i)
-                        for (unsigned int j = 0; j < 2; ++j)
-                          for (unsigned int k = 0; k < 2; ++k)
-                            for (unsigned int l = 0; l <= deg; ++l)
-                              {
-                                tmp (3 * (i + 2 * (j + 2 * k)))
-                                  -= this->restriction[index][2 * (2 * i + j) + k]
-                                     ((4 * i + j + 2) * this->degree + l, dof)
-                                  * this->shape_value_component
-                                    ((4 * i + j + 2) * this->degree + l,
-                                     quadrature_points[q_point], 0);
-                                tmp (3 * (i + 2 * (j + 2 * k)) + 1)
-                                  -= this->restriction[index][2 * (2 * i + j) + k]
-                                     ((4 * i + k) * this->degree + l, dof)
-                                  * this->shape_value_component
-                                    ((4 * i + k) * this->degree + l,
-                                     quadrature_points[q_point], 1);
-                                tmp (3 * (i + 2 * (j + 2 * k)) + 2)
-                                  -= this->restriction[index][2 * (2 * i + j) + k]
-                                     ((2 * (j + 4) + k) * this->degree + l,
-                                      dof)
-                                  * this->shape_value_component
-                                    ((2 * (j + 4) + k) * this->degree + l,
-                                     quadrature_points[q_point], 2);
-
-                                for (unsigned int m = 0; m < deg; ++m)
-                                  {
-                                    tmp (3 * (i + 2 * (j + 2 * k)))
-                                      -= this->restriction[index][2 * (2 * i + j) + k]
-                                         (((2 * j + 5) * deg + m)
-                                          * this->degree + l + n_edge_dofs,
-                                          dof)
-                                      * this->shape_value_component
-                                        (((2 * j + 5) * deg + m)
-                                         * this->degree + l + n_edge_dofs,
-                                         quadrature_points[q_point], 0);
-                                    tmp (3 * (i + 2 * (j + 2 * k)))
-                                      -= this->restriction[index][2 * (2 * i + j) + k]
-                                         ((2 * (i + 4) * this->degree + l)
-                                          * deg + m + n_edge_dofs, dof)
-                                      * this->shape_value_component
-                                        ((2 * (i + 4) * this->degree + l)
-                                         * deg + m + n_edge_dofs,
-                                         quadrature_points[q_point], 0);
-                                    tmp (3 * (i + 2 * (j + 2 * k)) + 1)
-                                      -= this->restriction[index][2 * (2 * i + j) + k]
-                                         ((2 * k * this->degree + l) * deg + m
-                                                                     + n_edge_dofs,
-                                          dof)
-                                      * this->shape_value_component
-                                        ((2 * k * this->degree + l) * deg + m
-                                                                    + n_edge_dofs,
-                                         quadrature_points[q_point], 1);
-                                    tmp (3 * (i + 2 * (j + 2 * k)) + 1)
-                                      -= this->restriction[index][2 * (2 * i + j) + k]
-                                         (((2 * i + 9) * deg + m)
-                                          * this->degree + l + n_edge_dofs,
-                                          dof)
-                                      * this->shape_value_component
-                                        (((2 * i + 9) * deg + m)
-                                         * this->degree + l + n_edge_dofs,
-                                         quadrature_points[q_point], 1);
-                                    tmp (3 * (i + 2 * (j + 2 * k)) + 2)
-                                      -= this->restriction[index][2 * (2 * i + j) + k]
-                                         (((2 * k + 1) * deg + m)
-                                          * this->degree + l + n_edge_dofs,
-                                          dof)
-                                      * this->shape_value_component
-                                        (((2 * k + 1) * deg + m)
-                                         * this->degree + l + n_edge_dofs,
-                                         quadrature_points[q_point], 2);
-                                    tmp (3 * (i + 2 * (j + 2 * k)) + 2)
-                                      -= this->restriction[index][2 * (2 * i + j) + k]
-                                         ((2 * (j + 2) * this->degree + l)
-                                          * deg + m + n_edge_dofs, dof)
-                                      * this->shape_value_component
-                                        ((2 * (j + 2) * this->degree + l)
-                                         * deg + m + n_edge_dofs,
-                                         quadrature_points[q_point], 2);
-                                  }
-                              }
-
-                      tmp *= quadrature.weight (q_point);
-
-                      for (unsigned int i = 0; i <= deg; ++i)
+                        }
+                      
+                      const double weight = 2.0 * quadrature.weight (q_point);
+                      
+                      for (unsigned int i = 0; i < this->degree; ++i)
                         {
                           const double L_i_0
-                            = legendre_polynomials[i].value
-                              (quadrature_points[q_point] (0));
+                            = weight * legendre_polynomials[i].value (quadrature_points[q_point] (0));
                           const double L_i_1
-                            = legendre_polynomials[i].value
-                              (quadrature_points[q_point] (1));
+                            = weight * legendre_polynomials[i].value (quadrature_points[q_point] (1));
                           const double L_i_2
-                            = legendre_polynomials[i].value
-                              (quadrature_points[q_point] (2));
-
+                            = weight * legendre_polynomials[i].value (quadrature_points[q_point] (2));
+                          
                           for (unsigned int j = 0; j < deg; ++j)
                             {
-                              const double l_j_0
-                                = L_i_0 * lobatto_polynomials[j + 2].value
-                                          (quadrature_points[q_point] (1));
-                              const double l_j_1
-                                = L_i_1 * lobatto_polynomials[j + 2].value
-                                          (quadrature_points[q_point] (0));
-                              const double l_j_2
-                                = L_i_2 * lobatto_polynomials[j + 2].value
-                                          (quadrature_points[q_point] (0));
-
+                              const double Le_j_0
+                                = legendre_polynomials[j + 1].value (quadrature_points[q_point] (0));
+                              const double Le_j_1
+                                = legendre_polynomials[j + 1].value (quadrature_points[q_point] (1));
+                              const double lo_j_0
+                                = lobatto_polynomials[j + 2].value (quadrature_points[q_point] (0));
+                              const double lo_j_1
+                                = lobatto_polynomials[j + 2].value (quadrature_points[q_point] (1));
+                              
                               for (unsigned int k = 0; k < deg; ++k)
                                 {
-                                  const double l_k_0
-                                    = l_j_0 * lobatto_polynomials[k + 2].value
-                                              (quadrature_points[q_point] (2));
+                                  const double L_k
+                                    = legendre_polynomials[k + 1].value (quadrature_points[q_point] (2));
                                   const double l_k_1
-                                    = l_j_1 * lobatto_polynomials[k + 2].value
-                                              (quadrature_points[q_point] (2));
+                                    = lobatto_polynomials[k + 2].value (quadrature_points[q_point] (1));
                                   const double l_k_2
-                                    = l_j_2 * lobatto_polynomials[k + 2].value
-                                              (quadrature_points[q_point] (1));
-
-                                  for (unsigned int l = 0; l < 8; ++l)
-                                    {
-                                      system_rhs ((i * deg + j) * deg + k,
-                                                  3 * l)
-                                        += tmp (3 * l) * l_k_0;
-                                      system_rhs ((i * deg + j) * deg + k,
-                                                  3 * l + 1)
-                                        += tmp (3 * l + 1) * l_k_1;
-                                      system_rhs ((i * deg + j) * deg + k,
-                                                  3 * l + 2)
-                                        += tmp (3 * l + 2) * l_k_2;
-                                    }
+                                    = lobatto_polynomials[k + 2].value (quadrature_points[q_point] (2));
+                                  
+                                  this->restriction[index][2 * (2 * l + m) + n]
+                                  ((i * deg + j) * deg + k + n_boundary_dofs, dof)
+                                    += L_i_0 * (this->shape_grad_component (dof, quadrature_point, 0)[2]
+                                                * lo_j_1 * L_k
+                                                + this->shape_grad_component (dof, quadrature_point, 0)[1]
+                                                * Le_j_1 * l_k_2
+                                                + this->shape_value_component (dof, quadrature_point, 0)
+                                                * lo_j_1 * l_k_2);
+                                  this->restriction[index][2 * (2 * l + m) + n]
+                                  ((i + (j + deg) * this->degree) * deg + k + n_boundary_dofs, dof)
+                                    += L_i_1 * (this->shape_grad_component (dof, quadrature_point, 1)[2]
+                                                * lo_j_0 * L_k
+                                                + this->shape_grad_component (dof, quadrature_point, 1)[0]
+                                                * Le_j_0 * l_k_2
+                                                + this->shape_value_component (dof, quadrature_point, 1)
+                                                * lo_j_0 * l_k_2);
+                                  this->restriction[index][2 * (2 * l + m) + n]
+                                  (i + (k + (j + 2 * deg) * deg) * this->degree + n_boundary_dofs, dof)
+                                    += L_i_2 * (this->shape_grad_component (dof, quadrature_point, 2)[1]
+                                                * lo_j_0
+                                                * legendre_polynomials[k + 1].value (quadrature_points[q_point] (1))
+                                                + this->shape_grad_component (dof, quadrature_point, 2)[0]
+                                                * Le_j_0 * l_k_1
+                                                + this->shape_value_component (dof, quadrature_point, 2)
+                                                * lo_j_0 * l_k_1);
                                 }
                             }
                         }
                     }
-
-                  system_matrix_inv.mmult (solution, system_rhs);
-
-                  for (unsigned int i = 0; i < 2; ++i)
-                    for (unsigned int j = 0; j < 2; ++j)
-                      for (unsigned int k = 0; k < 2; ++k)
-                        for (unsigned int l = 0; l <= deg; ++l)
-                          for (unsigned int m = 0; m < deg; ++m)
-                            for (unsigned int n = 0; n < deg; ++n)
-                              {
-                                if (std::abs (solution
-                                              ((l * deg + m) * deg + n,
-                                                   3 * (i + 2 * (j + 2 * k))))
-                                      > 1e-14)
-                                  this->restriction[index][2 * (2 * i + j) + k]
-                                  ((l * deg + m) * deg + n + n_boundary_dofs,
-                                   dof) = solution ((l * deg + m) * deg + n,
-                                                        3 * (i + 2 * (j + 2 * k)));
-
-                                if (std::abs (solution
-                                              ((l * deg + m) * deg + n,
-                                                   3 * (i + 2 * (j + 2 * k)) + 1))
-                                      > 1e-14)
-                                  this->restriction[index][2 * (2 * i + j) + k]
-                                  ((l + (m + deg) * this->degree) * deg + n
-                                                                  + n_boundary_dofs,
-                                   dof) = solution ((l * deg + m) * deg + n,
-                                                        3 * (i + 2 * (j + 2 * k)) + 1);
-
-                                if (std::abs (solution
-                                              ((l * deg + m) * deg + n,
-                                                   3 * (i + 2 * (j + 2 * k)) + 2))
-                                      > 1e-14)
-                                  this->restriction[index][2 * (2 * i + j) + k]
-                                  (l + ((m + 2 * deg) * deg + n) * this->degree
-                                     + n_boundary_dofs, dof)
-                                    = solution ((l * deg + m) * deg + n,
-                                                    3 * (i + 2 * (j + 2 * k)) + 2);
-                              }
                 }
             }
 
@@ -1992,6 +1077,7 @@ FE_Nedelec<dim>::get_dpo_vector (const unsigned int degree, bool dg)
   return dpo;
 }
 
+
 //---------------------------------------------------------------------------
 // Data field initialization
 //---------------------------------------------------------------------------
@@ -2507,8 +1593,9 @@ const
              + GeometryInfo<dim>::lines_per_face)
               = 1.0;
             interpolation_matrix
-            (i + (j + GeometryInfo<dim>::lines_per_face) * source_fe.degree,
-             i + (j + GeometryInfo<dim>::lines_per_face) * this->degree)
+            (i + (j + GeometryInfo<dim>::lines_per_face + source_fe.degree - 1)
+               * source_fe.degree,
+             i + (j + GeometryInfo<dim>::lines_per_face + this->degree-1) * this->degree)
               = 1.0;
           }
       }
@@ -2530,20 +1617,6 @@ FE_Nedelec<1>::get_subface_interpolation_matrix(
 
                    // In this function we compute the
                    // subface interpolation matrix.
-                   // This is done by a projection-
-                   // based interpolation. Therefore
-                   // we first interpolate the
-                   // shape functions of the higher
-                   // order element on the lowest
-                   // order edge shape functions.
-                   // Then the remaining part of
-                   // the interpolated shape
-                   // functions is projected on the
-                   // higher order edge shape
-                   // functions, the face shape
-                   // functions and the interior
-                   // shape functions (if they all
-                   // exist).
 template <int dim>
 void
 FE_Nedelec<dim>::get_subface_interpolation_matrix(
@@ -2551,9 +1624,9 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix(
   const unsigned int subface,
   FullMatrix<double>& interpolation_matrix) const
 {
-                                   // this is only implemented, if the
-                                   // source FE is also a
-                                   // Nedelec element
+                   // this is only implemented, if the
+                   // source FE is also a
+                   // Nedelec element
   typedef FE_Nedelec<dim> FEN;
   typedef FiniteElement<dim> FEL;
 
@@ -2567,12 +1640,12 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix(
           ExcDimensionMismatch (interpolation_matrix.n (),
                                 this->dofs_per_face));
 
-                                   // ok, source is a Nedelec element, so
-                                   // we will be able to do the work
+                   // ok, source is a Nedelec element, so
+                   // we will be able to do the work
   const FE_Nedelec<dim> &source_fe
     = dynamic_cast<const FE_Nedelec<dim>&> (source);
 
-                                   // Make sure, that the element,
+                   // Make sure, that the element,
                    // for which the DoFs should be
                    // constrained is the one with
                    // the higher polynomial degree.
@@ -2580,101 +1653,41 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix(
                    // also if this assertion is not
                    // satisfied. But the matrices
                    // produced in that case might
-                                   // lead to problems in the
+                   // lead to problems in the
                    // hp procedures, which use this
-                                   // method.
+                   // method.
   Assert (this->dofs_per_face <= source_fe.dofs_per_face,
           typename FEL::ExcInterpolationNotImplemented ());
   interpolation_matrix = 0.0;
-                           // Perform projection-based interpolation
-                                   // as usual.
+                   // Set the degrees of freedom
+                   // by interpolation as usual.
   const QGauss<1> edge_quadrature (source_fe.degree);
   const std::vector<Point<1> >&
     edge_quadrature_points = edge_quadrature.get_points ();
+  const std::vector<Polynomials::Polynomial<double> >&
+    legendre_polynomials
+      = Polynomials::Legendre::generate_complete_basis (source_fe.degree - 1);
   const unsigned int& n_edge_quadrature_points = edge_quadrature.size ();
 
   switch (dim)
     {
       case 2:
         {
-          for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
-            for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
-                 ++q_point)
-              {
-                const Point<dim> quadrature_point (0.0,
-                                                   0.5 * (edge_quadrature_points[q_point] (0)
-                                                          + subface));
-
-                interpolation_matrix (0, dof) += 0.5
-                                              * edge_quadrature.weight (q_point)
-                                              * this->shape_value_component
-                                                (dof, quadrature_point, 1);
-              }
-
-          if (source_fe.degree > 1)
+          for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
+               ++q_point)
             {
-              const std::vector<Polynomials::Polynomial<double> >&
-                legendre_polynomials
-                  = Polynomials::Legendre::generate_complete_basis (source_fe.degree - 1);
-              FullMatrix<double> system_matrix_inv (source_fe.degree - 1,
-                                                    source_fe.degree - 1);
-
-              {
-                FullMatrix<double> assembling_matrix (source_fe.degree - 1,
-                                                      n_edge_quadrature_points);
-
-                for (unsigned int q_point = 0;
-                     q_point < n_edge_quadrature_points; ++q_point)
-                  {
-                    const double weight
-                      = std::sqrt (edge_quadrature.weight (q_point));
-
-                    for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
-                      assembling_matrix (i, q_point) = weight
-                                                       * legendre_polynomials[i + 1].value
-                                                         (edge_quadrature_points[q_point] (0));
-                  }
-
-                FullMatrix<double> system_matrix (source_fe.degree - 1, source_fe.degree - 1);
-
-                assembling_matrix.mTmult (system_matrix, assembling_matrix);
-                system_matrix_inv.invert (system_matrix);
-              }
-
-              Vector<double> solution (source_fe.degree - 1);
-              Vector<double> system_rhs (source_fe.degree - 1);
-
+              const Point<dim> quadrature_point (0.0,
+                                                 0.5 * (edge_quadrature_points[q_point] (0)
+                                                        + subface));
+              const double weight = 0.5 * edge_quadrature.weight (q_point);
+                
               for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
                 {
-                  system_rhs = 0.0;
-
-                  for (unsigned int q_point = 0;
-                       q_point < n_edge_quadrature_points; ++q_point)
-                    {
-                      const Point<dim> quadrature_point_0 (0.0,
-                                                           0.5 * (edge_quadrature_points[q_point] (0)
-                                                                  + subface));
-                      const Point<dim> quadrature_point_1 (0.0,
-                                                           edge_quadrature_points[q_point] (0));
-                      const double tmp = edge_quadrature.weight (q_point)
-                                         * (0.5 * this->shape_value_component
-                                                  (dof, quadrature_point_0, 1)
-                                                - interpolation_matrix (0,
-                                                                        dof)
-                                                * source_fe.shape_value_component
-                                                  (0, quadrature_point_1, 1));
-
-                      for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
-                        system_rhs (i) += tmp
-                                       * legendre_polynomials[i + 1].value
-                                         (edge_quadrature_points[q_point] (0));
-                    }
-
-                  system_matrix_inv.vmult (solution, system_rhs);
-
-                  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
-                    if (std::abs (solution (i)) > 1e-14)
-                      interpolation_matrix (i + 1, dof) = solution (i);
+                  const double shape_value = weight * this->shape_value_component (dof, quadrature_point, 1);
+                  
+                  for (unsigned int i = 0; i < source_fe.dofs_per_face; ++i)
+                    interpolation_matrix (i, dof) += shape_value
+                                                  * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0));
                 }
             }
 
@@ -2686,143 +1699,52 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix(
           const double shifts[4][2] = { { 0.0, 0.0 }, { 1.0, 0.0 },
                                         { 0.0, 1.0 }, { 1.0, 1.0 } };
 
-          for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
-            for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
-                 ++q_point)
-              {
-                const double weight = 0.5 * edge_quadrature.weight (q_point);
-
-                for (unsigned int i = 0; i < 2; ++i)
-                  {
-                    Point<dim>
-                      quadrature_point (0.5 * (i + shifts[subface][0]),
-                                        0.5 * (edge_quadrature_points[q_point] (0)
-                                               + shifts[subface][1]),
-                                        0.0);
-
-                    interpolation_matrix (i * source_fe.degree, dof) += weight
-                                                                     * this->shape_value_component
-                                                                       (this->face_to_cell_index (dof, 4),
-                                                                        quadrature_point,
-                                                                        1);
-                    quadrature_point
-                      = Point<dim> (0.5 * (edge_quadrature_points[q_point] (0)
-                                           + shifts[subface][0]),
-                                    0.5 * (i + shifts[subface][1]), 0.0);
-                    interpolation_matrix ((i + 2) * source_fe.degree, dof)
-                      += weight * this->shape_value_component
-                                  (this->face_to_cell_index (dof, 4),
-                                   quadrature_point, 0);
-                  }
-              }
-
-          if (source_fe.degree > 1)
+          for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
+               ++q_point)
             {
-              const std::vector<Polynomials::Polynomial<double> >&
-                legendre_polynomials
-                  = Polynomials::Legendre::generate_complete_basis (source_fe.degree - 1);
-              FullMatrix<double> system_matrix_inv (source_fe.degree - 1,
-                                                    source_fe.degree - 1);
-
-              {
-                FullMatrix<double> assembling_matrix (source_fe.degree - 1,
-                                                      n_edge_quadrature_points);
-
-                for (unsigned int q_point = 0;
-                     q_point < n_edge_quadrature_points; ++q_point)
-                  {
-                    const double weight
-                      = std::sqrt (edge_quadrature.weight (q_point));
-
-                    for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
-                      assembling_matrix (i, q_point) = weight
-                                                       * legendre_polynomials[i + 1].value
-                                                         (edge_quadrature_points[q_point] (0));
-                  }
-
-                FullMatrix<double> system_matrix (source_fe.degree - 1, source_fe.degree - 1);
+              const double weight = 0.5 * edge_quadrature.weight (q_point);
 
-                assembling_matrix.mTmult (system_matrix, assembling_matrix);
-                system_matrix_inv.invert (system_matrix);
-              }
-
-              FullMatrix<double> solution (source_fe.degree - 1,
-                                           GeometryInfo<dim>::lines_per_face);
-              FullMatrix<double> system_rhs (source_fe.degree - 1,
-                                             GeometryInfo<dim>::lines_per_face);
-              Vector<double> tmp (GeometryInfo<dim>::lines_per_face);
-
-              for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
+              for (unsigned int i = 0; i < 2; ++i)
                 {
-                  system_rhs = 0.0;
-
-                  for (unsigned int q_point = 0;
-                       q_point < n_edge_quadrature_points; ++q_point)
+                  Point<dim>
+                    quadrature_point (0.5 * (i + shifts[subface][0]),
+                                      0.5 * (edge_quadrature_points[q_point] (0)
+                                             + shifts[subface][1]),
+                                      0.0);
+                    
+                  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
                     {
-                      const double weight = edge_quadrature.weight (q_point);
-
-                      for (unsigned int i = 0; i < 2; ++i)
-                        {
-                          Point<dim>
-                            quadrature_point_0
-                            (0.5 * (i + shifts[subface][0]),
-                             0.5 * (edge_quadrature_points[q_point] (0)
-                                    + shifts[subface][1]), 0.0);
-                          Point<dim> quadrature_point_1 (i,
-                                                         edge_quadrature_points[q_point] (0),
-                                                         0.0);
-
-                          tmp (i) = weight
-                                    * (0.5 * this->shape_value_component
-                                             (this->face_to_cell_index (dof, 4),
-                                              quadrature_point_0, 1)
-                                           - interpolation_matrix
-                                             (i * source_fe.degree, dof)
-                                           * source_fe.shape_value_component
-                                             (i * source_fe.degree,
-                                              quadrature_point_1, 1));
-                          quadrature_point_0
-                            = Point<dim> (0.5 * (edge_quadrature_points[q_point] (0)
-                                                 + shifts[subface][0]),
-                                          0.5 * (i + shifts[subface][1]),
-                                          0.0);
-                          quadrature_point_1
-                            = Point<dim> (edge_quadrature_points[q_point] (0),
-                                          i, 0.0);
-                          tmp (i + 2) = weight
-                                        * (0.5 * this->shape_value_component
-                                                 (this->face_to_cell_index (dof, 4),
-                                                  quadrature_point_0, 0)
-                                               - interpolation_matrix
-                                                 ((i + 2) * source_fe.degree,
-                                                  dof)
-                                               * source_fe.shape_value_component
-                                                 ((i + 2) * source_fe.degree,
-                                                  quadrature_point_1, 0));
-                        }
-
-                      for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
-                        {
-                          const double L_i
-                            = legendre_polynomials[i + 1].value
-                              (edge_quadrature_points[q_point] (0));
-
-                          for (unsigned int j = 0;
-                               j < GeometryInfo<dim>::lines_per_face; ++j)
-                            system_rhs (i, j) += tmp (j) * L_i;
-                        }
+                      const double shape_value
+                        = weight * this->shape_value_component (this->face_to_cell_index (dof, 4),
+                                                                quadrature_point, 1);
+                      
+                      for (unsigned int j = 0; j < source_fe.degree; ++j)
+                        interpolation_matrix (i * source_fe.degree + j, dof)
+                          += shape_value
+                          * legendre_polynomials[j].value (edge_quadrature_points[q_point] (0));
+                    }
+                  
+                  quadrature_point
+                    = Point<dim> (0.5 * (edge_quadrature_points[q_point] (0)
+                                         + shifts[subface][0]),
+                                  0.5 * (i + shifts[subface][1]), 0.0);
+                  
+                  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
+                    {
+                      const double shape_value
+                        = weight * this->shape_value_component (this->face_to_cell_index (dof, 4),
+                                                                quadrature_point, 0);
+                      
+                      for (unsigned int j = 0; j < source_fe.degree; ++j)
+                        interpolation_matrix ((i + 2) * source_fe.degree + j, dof)
+                          += shape_value
+                          * legendre_polynomials[j].value (edge_quadrature_points[q_point] (0));
                     }
-
-                  system_matrix_inv.mmult (solution, system_rhs);
-
-                  for (unsigned int i = 0;
-                       i < GeometryInfo<dim>::lines_per_face; ++i)
-                    for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
-                      if (std::abs (solution (j, i)) > 1e-14)
-                        interpolation_matrix (i * source_fe.degree + j + 1,
-                                              dof) = solution (j, i);
                 }
+            }
 
+          if (source_fe.degree > 1)
+            {
               const QGauss<2> quadrature (source_fe.degree);
               const std::vector<Point<2> >&
                 quadrature_points = quadrature.get_points ();
@@ -2830,126 +1752,59 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix(
                 lobatto_polynomials
                   = Polynomials::Lobatto::generate_complete_basis
                     (source_fe.degree);
-              const unsigned int n_boundary_dofs
+              const unsigned int n_edge_dofs
                 = GeometryInfo<dim>::lines_per_face * source_fe.degree;
               const unsigned int& n_quadrature_points = quadrature.size ();
 
-              {
-                FullMatrix<double>
-                  assembling_matrix (source_fe.degree * (source_fe.degree - 1),
-                                     n_quadrature_points);
-
-                for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                     ++q_point)
-                  {
-                    const double weight = std::sqrt (quadrature.weight (q_point));
-
-                    for (unsigned int i = 0; i < source_fe.degree; ++i)
-                      {
-                        const double L_i = weight
-                                           * legendre_polynomials[i].value
-                                             (quadrature_points[q_point] (0));
-
-                        for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
-                          assembling_matrix (i * (source_fe.degree - 1) + j,
-                                             q_point)
-                            = L_i * lobatto_polynomials[j + 2].value
-                                    (quadrature_points[q_point] (1));
-                      }
-                  }
-
-                FullMatrix<double> system_matrix (assembling_matrix.m (),
-                                                  assembling_matrix.m ());
-
-                assembling_matrix.mTmult (system_matrix, assembling_matrix);
-                system_matrix_inv.reinit (system_matrix.m (),
-                                          system_matrix.m ());
-                system_matrix_inv.invert (system_matrix);
-              }
-
-              solution.reinit (system_matrix_inv.m (), 2);
-              system_rhs.reinit (system_matrix_inv.m (), 2);
-              tmp.reinit (2);
-
-              for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
+              for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                   ++q_point)
                 {
-                  system_rhs = 0.0;
-
-                  for (unsigned int q_point = 0;
-                       q_point < n_quadrature_points; ++q_point)
+                  const double weight = 0.5 * quadrature.weight (q_point);
+                  const Point<dim>
+                    quadrature_point
+                    (0.5 * (quadrature_points[q_point] (0)
+                            + shifts[subface][0]),
+                     0.5 * (quadrature_points[q_point] (1)
+                            + shifts[subface][1]), 0.0);
+                  
+                  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
                     {
-                      Point<dim>
-                        quadrature_point
-                        (0.5 * (quadrature_points[q_point] (0)
-                                + shifts[subface][0]),
-                         0.5 * (quadrature_points[q_point] (1)
-                                + shifts[subface][1]), 0.0);
-                      tmp (0) = 0.5 * this->shape_value_component
-                                      (this->face_to_cell_index (dof, 4),
-                                       quadrature_point, 0);
-                      tmp (1) = 0.5 * this->shape_value_component
-                                      (this->face_to_cell_index (dof, 4),
-                                       quadrature_point, 1);
-                      quadrature_point
-                        = Point<dim> (quadrature_points[q_point] (0),
-                                      quadrature_points[q_point] (1), 0.0);
-
-                      for (unsigned int i = 0; i < 2; ++i)
-                        for (unsigned int j = 0; j < source_fe.degree; ++j)
-                          {
-                            tmp (0) -= interpolation_matrix
-                                       ((i + 2) * source_fe.degree + j, dof)
-                                    * source_fe.shape_value_component
-                                      ((i + 2) * source_fe.degree + j,
-                                       quadrature_point, 0);
-                            tmp (1) -= interpolation_matrix
-                                       (i * source_fe.degree + j, dof)
-                                    * source_fe.shape_value_component
-                                      (i * source_fe.degree + j,
-                                       quadrature_point, 1);
-                          }
-
-                      tmp *= quadrature.weight (q_point);
-
+                      const double shape_grad_0
+                        = this->shape_grad_component (this->face_to_cell_index (dof, 4),
+                                                      quadrature_point, 0)[1];
+                      const double shape_grad_1
+                        = this->shape_grad_component (this->face_to_cell_index (dof, 4),
+                                                      quadrature_point, 1)[0];
+                      const double shape_value_0
+                        = this->shape_value_component (this->face_to_cell_index (dof, 4),
+                                                       quadrature_point, 0);
+                      const double shape_value_1
+                        = this->shape_value_component (this->face_to_cell_index (dof, 4),
+                                                       quadrature_point, 1);
+                      
                       for (unsigned int i = 0; i < source_fe.degree; ++i)
                         {
-                          const double L_i_0 = legendre_polynomials[i].value
-                                               (quadrature_points[q_point] (0));
-                          const double L_i_1 = legendre_polynomials[i].value
-                                               (quadrature_points[q_point] (1));
-
+                          const double L_i_0
+                            = weight * legendre_polynomials[i].value (quadrature_points[q_point] (0));
+                          const double L_i_1
+                            = weight * legendre_polynomials[i].value (quadrature_points[q_point] (1));
+                          
                           for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
                             {
-                              system_rhs (i * (source_fe.degree - 1) + j, 0)
-                                += tmp (0) * L_i_0
-                                           * lobatto_polynomials[j + 2].value
-                                             (quadrature_points[q_point] (1));
-                              system_rhs (i * (source_fe.degree - 1) + j, 1)
-                                += tmp (1) * L_i_1
-                                           * lobatto_polynomials[j + 2].value
-                                             (quadrature_points[q_point] (0));
+                              interpolation_matrix (i * (source_fe.degree - 1) + j + n_edge_dofs, dof)
+                                += L_i_0 * (shape_grad_0
+                                            * legendre_polynomials[j + 1].value (quadrature_points[q_point] (1))
+                                            + shape_value_0
+                                            * lobatto_polynomials[j + 2].value (quadrature_points[q_point] (1)));
+                              interpolation_matrix (i + (j + source_fe.degree - 1) * source_fe.degree
+                                                      + n_edge_dofs, dof)
+                                += L_i_1 * (shape_grad_1
+                                            * legendre_polynomials[j + 1].value (quadrature_points[q_point] (0))
+                                            + shape_value_1
+                                            * lobatto_polynomials[j + 2].value (quadrature_points[q_point] (0)));
                             }
                         }
                     }
-
-                  system_matrix_inv.mmult (solution, system_rhs);
-
-                  for (unsigned int i = 0; i < source_fe.degree; ++i)
-                    for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
-                      {
-                        if (std::abs (solution (i * (source_fe.degree - 1) + j, 0))
-                              > 1e-14)
-                          interpolation_matrix (i * (source_fe.degree - 1)
-                                                  + j + n_boundary_dofs, dof)
-                            = solution (i * (source_fe.degree - 1) + j, 0);
-
-                        if (std::abs (solution (i * (source_fe.degree - 1) + j, 1))
-                              > 1e-14)
-                          interpolation_matrix (i + (j + source_fe.degree - 1)
-                                                  * source_fe.degree
-                                                  + n_boundary_dofs, dof)
-                            = solution (i * (source_fe.degree - 1) + j, 1);
-                      }
                 }
             }
 
@@ -2975,16 +1830,21 @@ void FE_Nedelec<dim>::interpolate (std::vector<double>&, const std::vector<doubl
                    // its values at the generalized support
                    // points in the finite element space on the
                    // reference cell.
-                   // This is done as usual by projection-based
-                   // interpolation.
+                   // The interpolation on the edge degrees of
+                   // freedom is done by direct calculation.
+                   // For the interpolation on the remaining
+                   // degrees of freedom we use a projection-
+                   // based interpolation scheme. Therefore
+                   // the remaining part of the interpolated
+                   // function is projected on the face shape
+                   // functions and the interior shape
+                   // functions (if they exist).
 template <int dim>
 void
 FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
                               const std::vector<Vector<double> >& values,
                               unsigned int offset) const
 {
-  const unsigned int deg = this->degree-1;
-
   Assert (values.size () == this->generalized_support_points.size (),
           ExcDimensionMismatch (values.size (),
                                 this->generalized_support_points.size ()));
@@ -2996,122 +1856,157 @@ FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
   std::fill (local_dofs.begin (), local_dofs.end (), 0.);
 
   if (offset < dim)
-    switch (dim)
-      {
-        case 2:
-          {
-            const QGauss<1> reference_edge_quadrature (this->degree);
-            const unsigned int& n_edge_points
-              = reference_edge_quadrature.size ();
-
+    {
+      const unsigned int deg = this->degree-1;
+      const std::vector<Polynomials::Polynomial<double> >&
+        legendre_polynomials
+        = Polynomials::Legendre::generate_complete_basis (deg);
+      const QGauss<1> edge_quadrature (this->degree);
+      const std::vector<Point<1> >&
+        edge_quadrature_points = edge_quadrature.get_points ();
+      const unsigned int& n_edge_points = edge_quadrature.size ();
+      
+      switch (dim)
+        {
+          case 2:
+            {
                                   // Let us begin with the
-                                  // interpolation part.
-            for (unsigned int i = 0; i < 2; ++i)
-              {
-                for (unsigned int q_point = 0; q_point < n_edge_points;
-                     ++q_point)
-                  local_dofs[i * this->degree]
-                    += reference_edge_quadrature.weight (q_point)
-                       * values[q_point + i * n_edge_points] (1);
-
+                                  // edge degrees of freedom.
+              for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
+                {
+                  const double weight = edge_quadrature.weight (q_point);
+                  
+                  for (unsigned int i = 0; i < this->degree; ++i)
+                    {
+                      const double L_i
+                        = weight * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0));
+                      
+                      for (unsigned int j = 0; j < 2; ++j)
+                        {
+                          local_dofs[i + j * this->degree]
+                            += L_i * values[q_point + j * n_edge_points] (1);
+                          
+                          if (offset == 0)
+                            local_dofs[i + (j + 2) * this->degree]
+                              += L_i * values[q_point + (j + 2) * n_edge_points] (0);
+                        }
+                    }
+                }
+            
                                   // Add the computed values
                                   // to the resulting vector
                                   // only, if they are not
                                   // too small.
-                if (std::abs (local_dofs[i * this->degree]) < 1e-14)
-                  local_dofs[i * this->degree] = 0.0;
-              }
+              for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
+                for (unsigned int i = 0; i < this->degree; ++i)
+                  if (std::abs (local_dofs[i + e * this->degree]) < 1e-14)
+                    local_dofs[i + e * this->degree] = 0.0;
 
-            if (offset == 0)
-              for (unsigned int i = 0; i < 2; ++i)
+                                  // If the degree is greater
+                                  // than 0, then we have still
+                                  // some interpolations onto
+                                  // the interior shape
+                                  // functions left.
+              if (this->degree > 1)
                 {
-                  for (unsigned int q_point = 0; q_point < n_edge_points;
-                       ++q_point)
-                    local_dofs[(i + 2) * this->degree]
-                      += reference_edge_quadrature.weight (q_point)
-                         * values[q_point + (i + 2) * n_edge_points] (0);
+                  const std::vector<Polynomials::Polynomial<double> >&
+                    lobatto_polynomials
+                      = Polynomials::Lobatto::generate_complete_basis
+                        (this->degree);
+                  std::vector<Polynomials::Polynomial<double> >
+                    lobatto_polynomials_grad (this->degree);
+
+                  for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
+                       ++i)
+                    lobatto_polynomials_grad[i]
+                      = lobatto_polynomials[i + 1].derivative ();
+
+                                  // We set up the system
+                                  // matrix and use it for
+                                  // both, the horizontal
+                                  // and the vertical
+                                  // interior shape
+                                  // functions.
+                  const QGauss<dim> reference_quadrature (this->degree);
+                  const unsigned int& n_interior_points
+                    = reference_quadrature.size ();
 
-                  if (std::abs (local_dofs[(i + 2) * this->degree]) < 1e-14)
-                    local_dofs[(i + 2) * this->degree] = 0.0;
-                }
+                  FullMatrix<double> system_matrix (deg * this->degree,
+                                                    deg * this->degree);
 
-                                  // If the degree is greater
-                                  // than 0, then we have still
-                                  // some higher order edge
-                                  // shape functions to
-                                  // consider.
-                                  // Here the projection part
-                                  // starts. The dof values
-                                  // are obtained by solving
-                                  // a linear system of
-                                  // equations.
-            if (this->degree > 1)
-              {
-                                          // We start with projection
-                                          // on the higher order edge
-                                          // shape function.
-                const std::vector<Polynomials::Polynomial<double> >&
-                  lobatto_polynomials
-                    = Polynomials::Lobatto::generate_complete_basis
-                      (this->degree);
-                const unsigned int
-                  line_coordinate[GeometryInfo<2>::lines_per_cell]
-                    = {1, 1, 0, 0};
-                std::vector<Polynomials::Polynomial<double> >
-                  lobatto_polynomials_grad (this->degree);
-
-                for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
-                     ++i)
-                  lobatto_polynomials_grad[i]
-                    = lobatto_polynomials[i + 1].derivative ();
-
-                                  // Set up the system matrix.
-                                  // This can be used for all
-                                  // edges.
-                FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
-
-                for (unsigned int i = 0; i < system_matrix.m (); ++i)
-                  for (unsigned int j = 0; j < system_matrix.n (); ++j)
-                    for (unsigned int q_point = 0; q_point < n_edge_points;
-                         ++q_point)
-                      system_matrix (i, j)
-                        += boundary_weights (q_point, j)
-                           * lobatto_polynomials_grad[i + 1].value
-                             (this->generalized_face_support_points[q_point]
-                              (1));
-
-                FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
+                  for (unsigned int i = 0; i < this->degree; ++i)
+                    for (unsigned int j = 0; j < deg; ++j)
+                      for (unsigned int k = 0; k < this->degree; ++k)
+                        for (unsigned int l = 0; l < deg; ++l)
+                          for (unsigned int q_point = 0;
+                               q_point < n_interior_points; ++q_point)
+                            system_matrix (i * deg + j, k * deg + l)
+                              += reference_quadrature.weight (q_point)
+                                 * legendre_polynomials[i].value
+                                   (this->generalized_support_points[q_point
+                                                                     + GeometryInfo<dim>::lines_per_cell
+                                                                     * n_edge_points]
+                                    (0))
+                                 * lobatto_polynomials[j + 2].value
+                                   (this->generalized_support_points[q_point
+                                                                     + GeometryInfo<dim>::lines_per_cell
+                                                                     * n_edge_points]
+                                    (1))
+                                 * lobatto_polynomials_grad[k].value
+                                   (this->generalized_support_points[q_point
+                                                                     + GeometryInfo<dim>::lines_per_cell
+                                                                     * n_edge_points]
+                                    (0))
+                                 * lobatto_polynomials[l + 2].value
+                                   (this->generalized_support_points[q_point
+                                                                     + GeometryInfo<dim>::lines_per_cell
+                                                                     * n_edge_points]
+                                    (1));
 
-                system_matrix_inv.invert (system_matrix);
+                  FullMatrix<double> system_matrix_inv (system_matrix.m (),
+                                                        system_matrix.m ());
+                
+                  system_matrix_inv.invert (system_matrix);
 
-                Vector<double> system_rhs (system_matrix.m ());
-                Vector<double> solution (system_rhs.size ());
+                  Vector<double> solution (system_matrix_inv.m ());
+                  Vector<double> system_rhs (system_matrix.m ());
 
-                for (unsigned int line = 0;
-                     line < GeometryInfo<dim>::lines_per_cell; ++line)
-                  if ((line < 2) || (offset == 0))
+                  if (offset == 0)
                     {
-                                  // Set up the right hand side.
-                      system_rhs = 0;
-
-                      for (unsigned int q_point = 0; q_point < n_edge_points;
-                           ++q_point)
+                                  // Set up the right hand side
+                                  // for the horizontal shape
+                                  // functions.
+                      for (unsigned int q_point = 0;
+                           q_point < n_interior_points; ++q_point)
                         {
-                          const double tmp
-                            = values[line * n_edge_points + q_point]
-                              (line_coordinate[line])
-                              - local_dofs[line * this->degree]
-                              * this->shape_value_component
-                                (line * this->degree,
-                                 this->generalized_support_points[line
-                                                                  * n_edge_points
-                                                                  + q_point],
-                                 line_coordinate[line]);
-
-                          for (unsigned int i = 0; i < system_rhs.size ();
-                               ++i)
-                            system_rhs (i) += boundary_weights (q_point, i)
-                                              * tmp;
+                          double tmp
+                            = values[q_point + GeometryInfo<dim>::lines_per_cell
+                                     * n_edge_points] (0);
+
+                          for (unsigned int i = 0; i < 2; ++i)
+                            for (unsigned int j = 0; j < this->degree; ++j)
+                              tmp -= local_dofs[(i + 2) * this->degree + j]
+                                     * this->shape_value_component
+                                       ((i + 2) * this->degree + j,
+                                        this->generalized_support_points[q_point
+                                                                         + GeometryInfo<dim>::lines_per_cell
+                                                                         * n_edge_points],
+                                        0);
+
+                          for (unsigned int i = 0; i < this->degree; ++i)
+                            for (unsigned int j = 0; j < this->degree-1; ++j)
+                              system_rhs (i * (this->degree-1) + j)
+                                += reference_quadrature.weight (q_point) * tmp
+                                   * lobatto_polynomials_grad[i].value
+                                     (this->generalized_support_points[q_point
+                                                                       + GeometryInfo<dim>::lines_per_cell
+                                                                       * n_edge_points]
+                                      (0))
+                                   * lobatto_polynomials[j + 2].value
+                                     (this->generalized_support_points[q_point
+                                                                       + GeometryInfo<dim>::lines_per_cell
+                                                                       * n_edge_points]
+                                      (1));
                         }
 
                       system_matrix_inv.vmult (solution, system_rhs);
@@ -3120,814 +2015,715 @@ FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
                                   // to the resulting vector
                                   // only, if they are not
                                   // too small.
-                      for (unsigned int i = 0; i < solution.size (); ++i)
-                        if (std::abs (solution (i)) > 1e-14)
-                          local_dofs[line * this->degree + i + 1]
-                            = solution (i);
+                      for (unsigned int i = 0; i < this->degree; ++i)
+                        for (unsigned int j = 0; j < deg; ++j)
+                          if (std::abs (solution (i * deg + j)) > 1e-14)
+                            local_dofs[(i + GeometryInfo<dim>::lines_per_cell)
+                                       * deg + j
+                                       + GeometryInfo<dim>::lines_per_cell]
+                              = solution (i * deg + j);
                     }
 
-                                  // Then we go on to the
-                                  // interior shape
-                                  // functions. Again we
-                                  // set up the system
-                                  // matrix and use it
-                                  // for both, the
-                                  // horizontal and the
-                                  // vertical, interior
-                                  // shape functions.
-                const QGauss<dim> reference_quadrature (this->degree);
-                const std::vector<Polynomials::Polynomial<double> >&
-                  legendre_polynomials
-                  = Polynomials::Legendre::generate_complete_basis (this->degree-1);
-                const unsigned int& n_interior_points
-                  = reference_quadrature.size ();
-
-                system_matrix.reinit ((this->degree-1) * this->degree,
-                                      (this->degree-1) * this->degree);
-                system_matrix = 0;
+                                  // Set up the right hand side
+                                  // for the vertical shape
+                                  // functions.
+                  system_rhs = 0;
 
-                for (unsigned int i = 0; i < this->degree; ++i)
-                  for (unsigned int j = 0; j < this->degree-1; ++j)
-                    for (unsigned int k = 0; k < this->degree; ++k)
-                      for (unsigned int l = 0; l < this->degree-1; ++l)
-                        for (unsigned int q_point = 0;
-                             q_point < n_interior_points; ++q_point)
-                          system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
-                            += reference_quadrature.weight (q_point)
-                               * legendre_polynomials[i].value
-                                 (this->generalized_support_points[q_point
-                                                                   + GeometryInfo<dim>::lines_per_cell
-                                                                   * n_edge_points]
-                                  (0))
-                               * lobatto_polynomials[j + 2].value
+                  for (unsigned int q_point = 0; q_point < n_interior_points;
+                       ++q_point)
+                    {
+                      double tmp
+                        = values[q_point + GeometryInfo<dim>::lines_per_cell
+                                 * n_edge_points] (1);
+
+                      for (unsigned int i = 0; i < 2; ++i)
+                        for (unsigned int j = 0; j < this->degree; ++j)
+                          tmp -= local_dofs[i * this->degree + j]
+                                 * this->shape_value_component
+                                   (i * this->degree + j,
+                                    this->generalized_support_points[q_point
+                                                                     + GeometryInfo<dim>::lines_per_cell
+                                                                     * n_edge_points],
+                                    1);
+
+                      for (unsigned i = 0; i < this->degree; ++i)
+                        for (unsigned int j = 0; j < deg; ++j)
+                          system_rhs (i * deg + j)
+                            += reference_quadrature.weight (q_point) * tmp
+                               * lobatto_polynomials_grad[i].value
                                  (this->generalized_support_points[q_point
                                                                    + GeometryInfo<dim>::lines_per_cell
                                                                    * n_edge_points]
                                   (1))
-                               * lobatto_polynomials_grad[k].value
-                                 (this->generalized_support_points[q_point
-                                                                   + GeometryInfo<dim>::lines_per_cell
-                                                                   * n_edge_points]
-                                  (0))
-                               * lobatto_polynomials[l + 2].value
+                               * lobatto_polynomials[j + 2].value
                                  (this->generalized_support_points[q_point
                                                                    + GeometryInfo<dim>::lines_per_cell
                                                                    * n_edge_points]
-                                  (1));
-
-                system_matrix_inv.reinit (system_matrix.m (),
-                                          system_matrix.m ());
-                system_matrix_inv.invert (system_matrix);
-                solution.reinit (system_matrix_inv.m ());
-                system_rhs.reinit (system_matrix.m ());
-
-                if (offset == 0)
-                  {
-                                  // Set up the right hand side
-                                  // for the horizontal shape
-                                  // functions.
-                    system_rhs = 0;
-
-                    for (unsigned int q_point = 0;
-                         q_point < n_interior_points; ++q_point)
-                      {
-                        double tmp
-                          = values[q_point + GeometryInfo<dim>::lines_per_cell
-                                   * n_edge_points] (0);
-
-                        for (unsigned int i = 0; i < 2; ++i)
-                          for (unsigned int j = 0; j < this->degree; ++j)
-                            tmp -= local_dofs[(i + 2) * this->degree + j]
-                                   * this->shape_value_component
-                                     ((i + 2) * this->degree + j,
-                                      this->generalized_support_points[q_point
-                                                                       + GeometryInfo<dim>::lines_per_cell
-                                                                       * n_edge_points],
-                                      0);
-
-                        for (unsigned int i = 0; i < this->degree; ++i)
-                          for (unsigned int j = 0; j < this->degree-1; ++j)
-                            system_rhs (i * (this->degree-1) + j)
-                              += reference_quadrature.weight (q_point) * tmp
-                                 * lobatto_polynomials_grad[i].value
-                                   (this->generalized_support_points[q_point
-                                                                     + GeometryInfo<dim>::lines_per_cell
-                                                                     * n_edge_points]
-                                    (0))
-                                 * lobatto_polynomials[j + 2].value
-                                   (this->generalized_support_points[q_point
-                                                                     + GeometryInfo<dim>::lines_per_cell
-                                                                     * n_edge_points]
-                                    (1));
-                      }
+                                  (0));
+                    }
 
-                    system_matrix_inv.vmult (solution, system_rhs);
+                  system_matrix_inv.vmult (solution, system_rhs);
 
                                   // Add the computed values
                                   // to the resulting vector
                                   // only, if they are not
                                   // too small.
-                    for (unsigned int i = 0; i < this->degree; ++i)
-                      for (unsigned int j = 0; j < this->degree-1; ++j)
-                        if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14)
-                           local_dofs[(i + GeometryInfo<dim>::lines_per_cell)
-                                      * (this->degree-1) + j
-                                      + GeometryInfo<dim>::lines_per_cell]
-                             = solution (i * (this->degree-1) + j);
-                  }
-
-                                  // Set up the right hand side
-                                  // for the vertical shape
-                                  // functions.
-                system_rhs = 0;
-
-                for (unsigned int q_point = 0; q_point < n_interior_points;
-                     ++q_point)
-                  {
-                    double tmp
-                      = values[q_point + GeometryInfo<dim>::lines_per_cell
-                               * n_edge_points] (1);
-
-                    for (unsigned int i = 0; i < 2; ++i)
-                      for (unsigned int j = 0; j < this->degree; ++j)
-                        tmp -= local_dofs[i * this->degree + j]
-                               * this->shape_value_component
-                                 (i * this->degree + j,
-                                  this->generalized_support_points[q_point
-                                                                   + GeometryInfo<dim>::lines_per_cell
-                                                                   * n_edge_points],
-                                  1);
-
-                    for (unsigned i = 0; i < this->degree; ++i)
-                      for (unsigned int j = 0; j < this->degree-1; ++j)
-                        system_rhs (i * (this->degree-1) + j)
-                          += reference_quadrature.weight (q_point) * tmp
-                             * lobatto_polynomials_grad[i].value
-                               (this->generalized_support_points[q_point
-                                                                 + GeometryInfo<dim>::lines_per_cell
-                                                                 * n_edge_points]
-                                (1))
-                             * lobatto_polynomials[j + 2].value
-                               (this->generalized_support_points[q_point
-                                                                 + GeometryInfo<dim>::lines_per_cell
-                                                                 * n_edge_points]
-                                (0));
-                  }
+                  for (unsigned int i = 0; i < this->degree; ++i)
+                    for (unsigned int j = 0; j < deg; ++j)
+                      if (std::abs (solution (i * deg + j)) > 1e-14)
+                        local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
+                                          + deg) * this->degree]
+                          = solution (i * deg + j);
+                }
 
-                system_matrix_inv.vmult (solution, system_rhs);
+              break;
+            }
 
+          case 3:
+            {
+                                  // Let us begin with the
+                                  // edge degrees of freedom.
+              for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
+                {
+                  const double weight = edge_quadrature.weight (q_point);
+                  
+                  for (unsigned int i = 0; i < this->degree; ++i)
+                    {
+                      const double L_i
+                        = weight * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0));
+                    
+                      for (unsigned int j = 0; j < 4; ++j)
+                        local_dofs[i + (j + 8) * this->degree]
+                          += L_i * values[q_point + (j + 8) * n_edge_points] (2);
+                    
+                      if (offset + 1 < dim)
+                        {
+                          for (unsigned int j = 0; j < 2; ++j)
+                            for (unsigned int k = 0; k < 2; ++k)
+                              local_dofs[i + (j + 4 * k) * this->degree]
+                                += L_i * values[q_point + (j + 4 * k) * n_edge_points] (1);
+                          
+                          if (offset == 0)
+                            for (unsigned int j = 0; j < 2; ++j)
+                              for (unsigned int k = 0; k < 2; ++k)
+                                local_dofs[i + (j + 4 * k + 2) * this->degree]
+                                  += L_i * values[q_point + (j + 4 * k + 2) * n_edge_points] (0);
+                        }
+                    }
+                }
+            
                                   // Add the computed values
                                   // to the resulting vector
                                   // only, if they are not
                                   // too small.
+              for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
                 for (unsigned int i = 0; i < this->degree; ++i)
-                  for (unsigned int j = 0; j < this->degree-1; ++j)
-                    if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14)
-                      local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
-                                      + this->degree-1) * this->degree]
-                        = solution (i * (this->degree-1) + j);
-              }
-
-            break;
-          }
+                  if (std::abs (local_dofs[i + e * this->degree]) < 1e-14)
+                    local_dofs[i + e * this->degree] = 0.0;
 
-        case 3:
-          {
-            const QGauss<1>
-              reference_edge_quadrature (this->degree);
-            const unsigned int&
-              n_edge_points = reference_edge_quadrature.size ();
+                                  // If the degree is greater
+                                  // than 0, then we have still
+                                  // some interpolation to the
+                                  // face and interior shape
+                                  // functions left.
+              if (this->degree > 1)
+                {
+                  const std::vector<Polynomials::Polynomial<double> >&
+                    lobatto_polynomials
+                      = Polynomials::Lobatto::generate_complete_basis
+                        (this->degree);
+                  FullMatrix<double> system_matrix (deg * this->degree,
+                                                    deg * this->degree);
+                  std::vector<Polynomials::Polynomial<double> >
+                    lobatto_polynomials_grad (this->degree);
+
+                  for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
+                       ++i)
+                    lobatto_polynomials_grad[i]
+                      = lobatto_polynomials[i + 1].derivative ();
+
+                                  // We set up the system
+                                  // matrix and use it for
+                                  // both, the horizontal
+                                  // and the vertical, shape
+                                  // functions.
+                  const unsigned int
+                    n_face_points = n_edge_points * n_edge_points;
 
-                                  // Let us begin with the
-                                  // interpolation part.
-            for (unsigned int i = 0; i < 4; ++i)
-              {
-                for (unsigned int q_point = 0; q_point < n_edge_points;
-                     ++q_point)
-                  local_dofs[(i + 8) * this->degree]
-                    += reference_edge_quadrature.weight (q_point)
-                       * values[q_point + (i + 8) * n_edge_points] (2);
+                  for (unsigned int i = 0; i < this->degree; ++i)
+                    for (unsigned int j = 0; j < deg; ++j)
+                      for (unsigned int k = 0; k < this->degree; ++k)
+                        for (unsigned int l = 0; l < deg; ++l)
+                          for (unsigned int q_point = 0; q_point < n_face_points;
+                               ++q_point)
+                            system_matrix (i * deg + j, k * deg + l)
+                              += boundary_weights (q_point + n_edge_points,
+                                                   2 * (k * deg + l))
+                                 * legendre_polynomials[i].value
+                                   (this->generalized_face_support_points[q_point
+                                                                          + 4
+                                                                          * n_edge_points]
+                                    (0))
+                                 * lobatto_polynomials[j + 2].value
+                                   (this->generalized_face_support_points[q_point
+                                                                          + 4
+                                                                          * n_edge_points]
+                                    (1));
 
-                                  // Add the computed values
-                                  // to the resulting vector
-                                  // only, if they are not
-                                  // too small.
-                if (std::abs (local_dofs[(i + 8) * this->degree]) < 1e-14)
-                  local_dofs[(i + 8) * this->degree] = 0.0;
-              }
+                  FullMatrix<double> system_matrix_inv (system_matrix.m (),
+                                                        system_matrix.m ());
+                  
+                  system_matrix_inv.invert (system_matrix);
+                  
+                  Vector<double> solution (system_matrix.m ());
+                  Vector<double> system_rhs (system_matrix.m ());
 
-            if (offset + 1 < dim)
-              {
-                for (unsigned int i = 0; i < 2; ++i)
-                  for (unsigned int j = 0; j < 2; ++j)
+                  for (unsigned int face = 0;
+                       face < GeometryInfo<dim>::faces_per_cell; ++face)
                     {
-                      for (unsigned int q_point = 0; q_point < n_edge_points;
-                           ++q_point)
-                        local_dofs[(i + 4 * j) * this->degree]
-                          += reference_edge_quadrature.weight (q_point)
-                             * values[q_point + (i + 4 * j) * n_edge_points]
-                               (1);
+                      switch (face)
+                        {
+                          case 0:
+                            {
+                              if (offset + 1 < dim)
+                                {
+                                      // Set up the right hand side
+                                      // for the horizontal shape
+                                      // functions.
+                                  system_rhs = 0;
 
-                                  // Add the computed values
-                                  // to the resulting vector
-                                  // only, if they are not
-                                  // too small.
-                      if (std::abs (local_dofs[(i + 4 * j) * this->degree])
-                            < 1e-14)
-                        local_dofs[(i + 4 * j) * this->degree] = 0.0;
-                    }
+                                  for (unsigned int q_point = 0;
+                                       q_point < n_face_points; ++q_point)
+                                    {
+                                      double tmp
+                                        = values[q_point
+                                                 + GeometryInfo<dim>::lines_per_cell
+                                                 * n_edge_points] (1);
 
-                if (offset == 0)
-                  for (unsigned int i = 0; i < 2; ++i)
-                    for (unsigned int j = 0; j < 2; ++j)
-                      {
-                        for (unsigned int q_point = 0;
-                             q_point < n_edge_points; ++q_point)
-                          local_dofs[(i + 4 * j + 2) * this->degree]
-                            += reference_edge_quadrature.weight (q_point)
-                               * values[q_point + (i + 4 * j + 2)
-                                        * n_edge_points] (0);
+                                      for (unsigned int i = 0; i < 2; ++i)
+                                        for (unsigned int j = 0; j < this->degree; ++j)
+                                          tmp
+                                            -= local_dofs[4 * i * this->degree
+                                                            + j]
+                                               * this->shape_value_component
+                                                 (4 * i * this->degree + j,
+                                                  this->generalized_support_points[q_point
+                                                                                   + GeometryInfo<dim>::lines_per_cell
+                                                                                   * n_edge_points],
+                                                  1);
 
-                                  // Add the computed values
-                                  // to the resulting vector
-                                  // only, if they are not
-                                  // too small.
-                        if (std::abs (local_dofs[(i + 4 * j + 2)
-                                                 * this->degree]) < 1e-14)
-                          local_dofs[(i + 4 * j + 2) * this->degree] = 0.0;
-                      }
-              }
+                                      for (unsigned int i = 0; i < this->degree; ++i)
+                                        for (unsigned int j = 0; j < deg; ++j)
+                                          system_rhs (i * deg + j)
+                                            += boundary_weights
+                                               (q_point + n_edge_points,
+                                                2 * (i * deg + j)) * tmp;
+                                    }
 
-                                  // If the degree is greater
-                                  // than 0, then we have still
-                                  // some higher order shape
-                                  // functions to consider.
-                                  // Here the projection part
-                                  // starts. The dof values
-                                  // are obtained by solving
-                                  // a linear system of
-                                  // equations.
-            if (this->degree > 1)
-              {
-                                          // We start with projection
-                                          // on the higher order edge
-                                          // shape function.
-                const std::vector<Polynomials::Polynomial<double> >&
-                  lobatto_polynomials
-                    = Polynomials::Lobatto::generate_complete_basis
-                      (this->degree);
-                const unsigned int
-                  line_coordinate[GeometryInfo<3>::lines_per_cell]
-                    = {1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
-                FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
-                FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
-                std::vector<Polynomials::Polynomial<double> >
-                  lobatto_polynomials_grad (this->degree);
-
-                for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
-                     ++i)
-                  lobatto_polynomials_grad[i]
-                    = lobatto_polynomials[i + 1].derivative ();
-
-                Vector<double> system_rhs (system_matrix.m ());
-                Vector<double> solution (system_rhs.size ());
-
-                                  // Set up the system matrix.
-                                  // This can be used for all
-                                  // edges.
-                for (unsigned int i = 0; i < system_matrix.m (); ++i)
-                  for (unsigned int j = 0; j < system_matrix.n (); ++j)
-                    for (unsigned int q_point = 0; q_point < n_edge_points;
-                         ++q_point)
-                      system_matrix (i, j)
-                        += boundary_weights (q_point, j)
-                           * lobatto_polynomials_grad[i + 1].value
-                             (this->generalized_face_support_points[q_point]
-                              (1));
+                                  system_matrix_inv.vmult (solution, system_rhs);
 
-                system_matrix_inv.invert (system_matrix);
+                                      // Add the computed values
+                                      // to the resulting vector
+                                      // only, if they are not
+                                      // too small.
+                                  for (unsigned int i = 0; i < this->degree; ++i)
+                                    for (unsigned int j = 0; j < deg; ++j)
+                                      if (std::abs (solution (i * deg + j))
+                                            > 1e-14)
+                                        local_dofs[(i
+                                                    + GeometryInfo<dim>::lines_per_cell)
+                                                   * deg + j
+                                                   + GeometryInfo<dim>::lines_per_cell]
+                                          = solution (i * deg + j);
+                                }
 
-                for (unsigned int line = 0;
-                     line < GeometryInfo<dim>::lines_per_cell; ++line)
-                  {
-                                  // Set up the right hand side.
-                    system_rhs = 0;
+                                  // Set up the right hand side
+                                  // for the vertical shape
+                                  // functions.
+                              system_rhs = 0;
 
-                    if ((((line == 0) || (line == 1) || (line == 4) ||
-                          (line == 5)) && (offset + 1 < dim)) ||
-                        (((line == 2) || (line == 3) || (line == 6) ||
-                          (line == 7)) && (offset == 0)) || (line > 7))
-                      {
-                        for (unsigned int q_point = 0; q_point < n_edge_points;
-                             ++q_point)
-                          {
-                            double tmp
-                              = values[line * n_edge_points + q_point]
-                                (line_coordinate[line])
-                                - local_dofs[line * this->degree]
-                                * this->shape_value_component
-                                  (line * this->degree,
-                                   this->generalized_support_points[line
-                                                                    * this->degree
-                                                                    + q_point],
-                                   line_coordinate[line]);
-
-                            for (unsigned int i = 0; i < system_rhs.size ();
-                                 ++i)
-                              system_rhs (i)
-                                += boundary_weights (q_point, i) * tmp;
-                          }
+                              for (unsigned int q_point = 0;
+                                   q_point < n_face_points; ++q_point)
+                                {
+                                  double tmp
+                                    = values[q_point
+                                             + GeometryInfo<dim>::lines_per_cell
+                                             * n_edge_points] (2);
+
+                                  for (unsigned int i = 0; i < 2; ++i)
+                                    for (unsigned int j = 0; j < this->degree; ++j)
+                                      tmp -= local_dofs[2 * (i + 4)
+                                                        * this->degree + j]
+                                             * this->shape_value_component
+                                               (2 * (i + 4) * this->degree + j,
+                                                this->generalized_support_points[q_point
+                                                                                 + GeometryInfo<dim>::lines_per_cell
+                                                                                 * n_edge_points],
+                                                2);
+
+                                  for (unsigned i = 0; i < this->degree; ++i)
+                                    for (unsigned int j = 0; j < deg; ++j)
+                                      system_rhs (i * deg + j)
+                                        += boundary_weights
+                                           (q_point + n_edge_points,
+                                            2 * (i * deg + j) + 1)
+                                           * tmp;
+                                }
 
-                        system_matrix_inv.vmult (solution, system_rhs);
+                              system_matrix_inv.vmult (solution, system_rhs);
 
                                   // Add the computed values
                                   // to the resulting vector
                                   // only, if they are not
                                   // too small.
-                        for (unsigned int i = 0; i < solution.size (); ++i)
-                          if (std::abs (solution (i)) > 1e-14)
-                            local_dofs[line * this->degree + i + 1]
-                              = solution (i);
-                      }
-                  }
-
-                                  // Then we go on to the
-                                  // face shape functions.
-                                  // Again we set up the
-                                  // system matrix and
-                                  // use it for both, the
-                                  // horizontal and the
-                                  // vertical, shape
-                                  // functions.
-                const std::vector<Polynomials::Polynomial<double> >&
-                  legendre_polynomials
-                    = Polynomials::Legendre::generate_complete_basis (this->degree-1);
-                const unsigned int
-                  n_face_points = n_edge_points * n_edge_points;
-
-                system_matrix.reinit ((this->degree-1) * this->degree,
-                                      (this->degree-1) * this->degree);
-                system_matrix = 0;
-
-                for (unsigned int i = 0; i < this->degree; ++i)
-                  for (unsigned int j = 0; j < this->degree-1; ++j)
-                    for (unsigned int k = 0; k < this->degree; ++k)
-                      for (unsigned int l = 0; l < this->degree-1; ++l)
-                        for (unsigned int q_point = 0; q_point < n_face_points;
-                             ++q_point)
-                          system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
-                            += boundary_weights (q_point + n_edge_points,
-                                                 2 * (k * (this->degree-1) + l))
-                               * legendre_polynomials[i].value
-                                 (this->generalized_face_support_points[q_point
-                                                                        + 4
-                                                                        * n_edge_points]
-                                  (0))
-                               * lobatto_polynomials[j + 2].value
-                                 (this->generalized_face_support_points[q_point
-                                                                        + 4
-                                                                        * n_edge_points]
-                                  (1));
-
-                system_matrix_inv.reinit (system_matrix.m (),
-                                          system_matrix.n ());
-                system_matrix_inv.invert (system_matrix);
-                solution.reinit (system_matrix.m ());
-                system_rhs.reinit (system_matrix.m ());
-
-                for (unsigned int face = 0;
-                     face < GeometryInfo<dim>::faces_per_cell; ++face)
-                  {
-                    switch (face)
-                      {
-                        case 0:
-                          {
-                            if (offset + 1 < dim)
-                              {
-                                  // Set up the right hand side
-                                  // for the horizontal shape
-                                  // functions.
-                                system_rhs = 0;
-
-                                for (unsigned int q_point = 0;
-                                     q_point < n_face_points; ++q_point)
-                                  {
-                                    double tmp
-                                      = values[q_point
-                                               + GeometryInfo<dim>::lines_per_cell
-                                               * n_edge_points] (1);
-
-                                    for (unsigned int i = 0; i < 2; ++i)
-                                      for (unsigned int j = 0; j < this->degree; ++j)
-                                        tmp
-                                          -= local_dofs[4 * i * this->degree
-                                                        + j]
-                                             * this->shape_value_component
-                                               (4 * i * this->degree + j,
-                                                this->generalized_support_points[q_point
-                                                                                 + GeometryInfo<dim>::lines_per_cell
-                                                                                 * n_edge_points],
-                                                1);
+                              for (unsigned int i = 0; i < this->degree; ++i)
+                                for (unsigned int j = 0; j < deg; ++j)
+                                  if (std::abs (solution (i * deg + j)) > 1e-14)
+                                    local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
+                                                      + deg)
+                                                 * this->degree]
+                                      = solution (i * deg + j);
 
-                                    for (unsigned int i = 0; i < this->degree; ++i)
-                                      for (unsigned int j = 0; j < this->degree-1; ++j)
-                                        system_rhs (i * (this->degree-1) + j)
-                                         += boundary_weights
-                                            (q_point + n_edge_points,
-                                             2 * (i * (this->degree-1) + j)) * tmp;
-                                  }
+                              break;
+                            }
 
-                                system_matrix_inv.vmult (solution, system_rhs);
+                          case 1:
+                            {
+                              if (offset + 1 < dim)
+                                {
+                                      // Set up the right hand side
+                                      // for the horizontal shape
+                                      // functions.
+                                  system_rhs = 0;
 
-                                  // Add the computed values
-                                  // to the resulting vector
-                                  // only, if they are not
-                                  // too small.
-                                for (unsigned int i = 0; i < this->degree; ++i)
-                                  for (unsigned int j = 0; j < this->degree-1; ++j)
-                                    if (std::abs (solution (i * (this->degree-1) + j))
-                                          > 1e-14)
-                                      local_dofs[(i
-                                                  + GeometryInfo<dim>::lines_per_cell)
-                                                 * (this->degree-1) + j
-                                                 + GeometryInfo<dim>::lines_per_cell]
-                                        = solution (i * (this->degree-1) + j);
-                              }
+                                  for (unsigned int q_point = 0;
+                                       q_point < n_face_points; ++q_point)
+                                    {
+                                      double tmp
+                                        = values[q_point
+                                                 + GeometryInfo<dim>::lines_per_cell
+                                                 * n_edge_points
+                                                 + n_face_points] (1);
 
-                                  // Set up the right hand side
-                                  // for the vertical shape
-                                  // functions.
-                            system_rhs = 0;
+                                      for (unsigned int i = 0; i < 2; ++i)
+                                        for (unsigned int j = 0; j < this->degree; ++j)
+                                          tmp -= local_dofs[(4 * i + 1)
+                                                            * this->degree + j]
+                                                 * this->shape_value_component
+                                                   ((4 * i + 1) * this->degree
+                                                    + j,
+                                                    this->generalized_support_points[q_point
+                                                                                     + GeometryInfo<dim>::lines_per_cell
+                                                                                     * n_edge_points
+                                                                                     + n_face_points],
+                                                    1);
 
-                            for (unsigned int q_point = 0;
-                                 q_point < n_face_points; ++q_point)
-                              {
-                                double tmp
-                                  = values[q_point
-                                           + GeometryInfo<dim>::lines_per_cell
-                                           * n_edge_points] (2);
-
-                                for (unsigned int i = 0; i < 2; ++i)
-                                  for (unsigned int j = 0; j < this->degree; ++j)
-                                    tmp -= local_dofs[2 * (i + 4)
-                                                      * this->degree + j]
-                                           * this->shape_value_component
-                                             (2 * (i + 4) * this->degree + j,
-                                              this->generalized_support_points[q_point
-                                                                               + GeometryInfo<dim>::lines_per_cell
-                                                                               * n_edge_points],
-                                              2);
-
-                                for (unsigned i = 0; i < this->degree; ++i)
-                                  for (unsigned int j = 0; j < this->degree-1; ++j)
-                                    system_rhs (i * (this->degree-1) + j)
-                                      += boundary_weights
-                                         (q_point + n_edge_points,
-                                          2 * (i * (this->degree-1) + j) + 1)
-                                         * tmp;
-                              }
+                                      for (unsigned int i = 0; i < this->degree; ++i)
+                                        for (unsigned int j = 0; j < deg; ++j)
+                                          system_rhs (i * deg + j)
+                                            += boundary_weights
+                                               (q_point + n_edge_points,
+                                                2 * (i * deg + j)) * tmp;
+                                    }
 
-                            system_matrix_inv.vmult (solution, system_rhs);
+                                  system_matrix_inv.vmult (solution, system_rhs);
 
-                                  // Add the computed values
-                                  // to the resulting vector
-                                  // only, if they are not
-                                  // too small.
-                            for (unsigned int i = 0; i < this->degree; ++i)
-                              for (unsigned int j = 0; j < this->degree-1; ++j)
-                                if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14)
-                                  local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
-                                                  + this->degree-1)
-                                             * this->degree]
-                                    = solution (i * (this->degree-1) + j);
-
-                            break;
-                          }
+                                      // Add the computed values
+                                      // to the resulting vector
+                                      // only, if they are not
+                                      // too small.
+                                  for (unsigned int i = 0; i < this->degree; ++i)
+                                    for (unsigned int j = 0; j < deg; ++j)
+                                      if (std::abs (solution (i * deg + j))
+                                            > 1e-14)
+                                        local_dofs[(i + GeometryInfo<dim>::lines_per_cell
+                                                      + 2 * this->degree) * deg + j
+                                                      + GeometryInfo<dim>::lines_per_cell]
+                                          = solution (i * deg + j);
+                                }
 
-                        case 1:
-                          {
-                            if (offset + 1 < dim)
-                              {
                                   // Set up the right hand side
-                                  // for the horizontal shape
+                                  // for the vertical shape
                                   // functions.
-                                system_rhs = 0;
+                              system_rhs = 0;
 
-                                for (unsigned int q_point = 0;
-                                     q_point < n_face_points; ++q_point)
-                                  {
-                                    double tmp
-                                      = values[q_point
-                                               + GeometryInfo<dim>::lines_per_cell
-                                               * n_edge_points
-                                               + n_face_points] (1);
+                              for (unsigned int q_point = 0;
+                                   q_point < n_face_points; ++q_point)
+                                {
+                                  double tmp
+                                    = values[q_point
+                                             + GeometryInfo<dim>::lines_per_cell
+                                             * n_edge_points + n_face_points]
+                                      (2);
 
-                                    for (unsigned int i = 0; i < 2; ++i)
-                                      for (unsigned int j = 0; j <= deg; ++j)
-                                        tmp -= local_dofs[(4 * i + 1)
-                                                          * this->degree + j]
-                                               * this->shape_value_component
-                                                 ((4 * i + 1) * this->degree
-                                                  + j,
-                                                  this->generalized_support_points[q_point
-                                                                                   + GeometryInfo<dim>::lines_per_cell
-                                                                                   * n_edge_points
-                                                                                   + n_face_points],
-                                                  1);
+                                  for (unsigned int i = 0; i < 2; ++i)
+                                    for (unsigned int j = 0; j < this->degree; ++j)
+                                      tmp -= local_dofs[(2 * (i + 4) + 1)
+                                                        * this->degree + j]
+                                             * this->shape_value_component
+                                               ((2 * (i + 4) + 1) * this->degree
+                                                + j,
+                                                this->generalized_support_points[q_point
+                                                                                 + GeometryInfo<dim>::lines_per_cell
+                                                                                 * n_edge_points
+                                                                                 + n_face_points],
+                                                2);
 
-                                    for (unsigned int i = 0; i <= deg; ++i)
-                                      for (unsigned int j = 0; j < deg; ++j)
-                                        system_rhs (i * deg + j)
-                                          += boundary_weights
-                                             (q_point + n_edge_points,
-                                              2 * (i * deg + j)) * tmp;
-                                  }
+                                  for (unsigned i = 0; i < this->degree; ++i)
+                                    for (unsigned int j = 0; j < deg; ++j)
+                                      system_rhs (i * deg + j)
+                                        += boundary_weights
+                                           (q_point + n_edge_points,
+                                            2 * (i * deg + j) + 1) * tmp;
+                                }
 
-                                system_matrix_inv.vmult (solution, system_rhs);
+                              system_matrix_inv.vmult (solution, system_rhs);
 
                                   // Add the computed values
                                   // to the resulting vector
                                   // only, if they are not
                                   // too small.
-                                for (unsigned int i = 0; i <= deg; ++i)
-                                  for (unsigned int j = 0; j < deg; ++j)
-                                    if (std::abs (solution (i * deg + j))
-                                          > 1e-14)
-                                      local_dofs[(i + GeometryInfo<dim>::lines_per_cell
-                                                  + 2 * this->degree) * deg + j
-                                                  + GeometryInfo<dim>::lines_per_cell]
-                                        = solution (i * deg + j);
-                              }
-
-                                  // Set up the right hand side
-                                  // for the vertical shape
-                                  // functions.
-                            system_rhs = 0;
-
-                            for (unsigned int q_point = 0;
-                                 q_point < n_face_points; ++q_point)
-                              {
-                                double tmp
-                                  = values[q_point
-                                           + GeometryInfo<dim>::lines_per_cell
-                                           * n_edge_points + n_face_points]
-                                    (2);
-
-                                for (unsigned int i = 0; i < 2; ++i)
-                                  for (unsigned int j = 0; j <= deg; ++j)
-                                    tmp -= local_dofs[(2 * (i + 4) + 1)
-                                                      * this->degree + j]
-                                           * this->shape_value_component
-                                             ((2 * (i + 4) + 1) * this->degree
-                                              + j,
-                                              this->generalized_support_points[q_point
-                                                                               + GeometryInfo<dim>::lines_per_cell
-                                                                               * n_edge_points
-                                                                               + n_face_points],
-                                              2);
-
-                                for (unsigned i = 0; i <= deg; ++i)
-                                  for (unsigned int j = 0; j < deg; ++j)
-                                    system_rhs (i * deg + j)
-                                      += boundary_weights
-                                         (q_point + n_edge_points,
-                                          2 * (i * deg + j) + 1) * tmp;
-                              }
-
-                            system_matrix_inv.vmult (solution, system_rhs);
+                              for (unsigned int i = 0; i < this->degree; ++i)
+                                for (unsigned int j = 0; j < deg; ++j)
+                                  if (std::abs (solution (i * deg + j)) > 1e-14)
+                                    local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
+                                                      + 3 * deg)
+                                                 * this->degree]
+                                      = solution (i * deg + j);
 
-                                  // Add the computed values
-                                  // to the resulting vector
-                                  // only, if they are not
-                                  // too small.
-                            for (unsigned int i = 0; i <= deg; ++i)
-                              for (unsigned int j = 0; j < deg; ++j)
-                                if (std::abs (solution (i * deg + j)) > 1e-14)
-                                  local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
-                                                  + 3 * deg)
-                                             * this->degree]
-                                    = solution (i * deg + j);
-
-                            break;
-                          }
+                              break;
+                            }
 
-                        case 2:
-                          {
-                            if (offset == 0)
-                              {
-                                  // Set up the right hand side
-                                  // for the horizontal shape
-                                  // functions.
-                                system_rhs = 0;
+                          case 2:
+                            {
+                              if (offset == 0)
+                                {
+                                      // Set up the right hand side
+                                      // for the horizontal shape
+                                      // functions.
+                                  system_rhs = 0;
 
-                                for (unsigned int q_point = 0;
-                                     q_point < n_face_points; ++q_point)
-                                  {
-                                    double tmp
-                                      = values[q_point
-                                               + GeometryInfo<dim>::lines_per_cell
-                                               * n_edge_points + 2 * n_face_points]
-                                        (2);
+                                  for (unsigned int q_point = 0;
+                                       q_point < n_face_points; ++q_point)
+                                    {
+                                      double tmp
+                                        = values[q_point
+                                                 + GeometryInfo<dim>::lines_per_cell
+                                                 * n_edge_points + 2 * n_face_points]
+                                          (2);
 
-                                    for (unsigned int i = 0; i < 2; ++i)
-                                      for (unsigned int j = 0; j <= deg; ++j)
-                                        tmp -= local_dofs[(i + 8) * this->degree
-                                                          + j]
-                                            * this->shape_value_component
-                                              ((i + 8) * this->degree + j,
-                                               this->generalized_support_points[q_point
-                                                                                + GeometryInfo<dim>::lines_per_cell
-                                                                                * n_edge_points
-                                                                                + 2
-                                                                                * n_face_points],
-                                               2);
-
-                                    for (unsigned i = 0; i <= deg; ++i)
-                                      for (unsigned int j = 0; j < deg; ++j)
-                                        system_rhs (i * deg + j)
-                                          += boundary_weights
-                                             (q_point + n_edge_points,
-                                              2 * (i * deg + j)) * tmp;
-                                  }
+                                      for (unsigned int i = 0; i < 2; ++i)
+                                        for (unsigned int j = 0; j < this->degree; ++j)
+                                          tmp -= local_dofs[(i + 8) * this->degree
+                                                                    + j]
+                                              * this->shape_value_component
+                                                ((i + 8) * this->degree + j,
+                                                 this->generalized_support_points[q_point
+                                                                                  + GeometryInfo<dim>::lines_per_cell
+                                                                                  * n_edge_points
+                                                                                  + 2
+                                                                                  * n_face_points],
+                                                 2);
+
+                                      for (unsigned i = 0; i < this->degree; ++i)
+                                        for (unsigned int j = 0; j < deg; ++j)
+                                          system_rhs (i * deg + j)
+                                            += boundary_weights
+                                               (q_point + n_edge_points,
+                                                2 * (i * deg + j)) * tmp;
+                                    }
 
-                                system_matrix_inv.vmult (solution, system_rhs);
+                                  system_matrix_inv.vmult (solution, system_rhs);
 
-                                  // Add the computed values
-                                  // to the resulting vector
-                                  // only, if they are not
-                                  // too small.
-                                for (unsigned int i = 0; i <= deg; ++i)
-                                  for (unsigned int j = 0; j < deg; ++j)
-                                    if (std::abs (solution (i * deg + j))
-                                          > 1e-14)
-                                      local_dofs[(i + GeometryInfo<dim>::lines_per_cell
-                                                  + 4 * this->degree) * deg
-                                                 + j
-                                                 + GeometryInfo<dim>::lines_per_cell]
-                                        = solution (i * deg + j);
-                              }
+                                      // Add the computed values
+                                      // to the resulting vector
+                                      // only, if they are not
+                                      // too small.
+                                  for (unsigned int i = 0; i < this->degree; ++i)
+                                    for (unsigned int j = 0; j < deg; ++j)
+                                      if (std::abs (solution (i * deg + j))
+                                            > 1e-14)
+                                        local_dofs[(i + GeometryInfo<dim>::lines_per_cell
+                                                      + 4 * this->degree) * deg
+                                                      + j
+                                                      + GeometryInfo<dim>::lines_per_cell]
+                                          = solution (i * deg + j);
+                                }
 
                                   // Set up the right hand side
                                   // for the vertical shape
                                   // functions.
-                            system_rhs = 0;
+                              system_rhs = 0;
 
-                            for (unsigned int q_point = 0;
-                                 q_point < n_face_points; ++q_point)
-                              {
-                                double tmp
-                                  = values[q_point
-                                           + GeometryInfo<dim>::lines_per_cell
-                                           * n_edge_points
-                                           + 2 * n_face_points] (0);
-
-                                for (unsigned int i = 0; i < 2; ++i)
-                                  for (unsigned int j = 0; j <= deg; ++j)
-                                    tmp -= local_dofs[(4 * i + 2)
-                                                      * this->degree + j]
-                                        * this->shape_value_component
-                                          ((4 * i + 2) * this->degree
-                                                       + j,
-                                           this->generalized_support_points[q_point
-                                                                            + GeometryInfo<dim>::lines_per_cell
-                                                                            * n_edge_points
-                                                                            + 2
-                                                                            * n_face_points],
-                                           0);
-
-                                for (unsigned int i = 0; i <= deg; ++i)
-                                  for (unsigned int j = 0; j < deg; ++j)
-                                    system_rhs (i * deg + j)
-                                      += boundary_weights
-                                         (q_point + n_edge_points,
-                                          2 * (i * deg + j) + 1) * tmp;
-                              }
+                              for (unsigned int q_point = 0;
+                                   q_point < n_face_points; ++q_point)
+                                {
+                                  double tmp
+                                    = values[q_point
+                                             + GeometryInfo<dim>::lines_per_cell
+                                             * n_edge_points
+                                             + 2 * n_face_points] (0);
+
+                                  for (unsigned int i = 0; i < 2; ++i)
+                                    for (unsigned int j = 0; j < this->degree; ++j)
+                                      tmp -= local_dofs[(4 * i + 2)
+                                                        * this->degree + j]
+                                          * this->shape_value_component
+                                            ((4 * i + 2) * this->degree
+                                                         + j,
+                                             this->generalized_support_points[q_point
+                                                                              + GeometryInfo<dim>::lines_per_cell
+                                                                              * n_edge_points
+                                                                              + 2
+                                                                              * n_face_points],
+                                             0);
+
+                                  for (unsigned int i = 0; i < this->degree; ++i)
+                                    for (unsigned int j = 0; j < deg; ++j)
+                                      system_rhs (i * deg + j)
+                                        += boundary_weights
+                                           (q_point + n_edge_points,
+                                            2 * (i * deg + j) + 1) * tmp;
+                                }
 
-                            system_matrix_inv.vmult (solution, system_rhs);
+                              system_matrix_inv.vmult (solution, system_rhs);
 
                                   // Add the computed values
                                   // to the resulting vector
                                   // only, if they are not
                                   // too small.
-                            for (unsigned int i = 0; i <= deg; ++i)
-                              for (unsigned int j = 0; j < deg; ++j)
-                                if (std::abs (solution (i * deg + j)) > 1e-14)
-                                 local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
-                                                 + 5 * deg) * this->degree]
-                                   = solution (i * deg + j);
-
-                            break;
-                          }
+                              for (unsigned int i = 0; i < this->degree; ++i)
+                                for (unsigned int j = 0; j < deg; ++j)
+                                  if (std::abs (solution (i * deg + j)) > 1e-14)
+                                    local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
+                                                       + 5 * deg) * this->degree]
+                                      = solution (i * deg + j);
 
-                        case 3:
-                          {
-                            if (offset == 0)
-                              {
-                                  // Set up the right hand side
-                                  // for the horizontal shape
-                                  // functions.
-                                system_rhs = 0;
+                              break;
+                            }
 
-                                for (unsigned int q_point = 0;
-                                     q_point < n_face_points; ++q_point)
-                                  {
-                                    double tmp
-                                      = values[q_point
-                                               + GeometryInfo<dim>::lines_per_cell
-                                               * n_edge_points + 3 * n_face_points]
-                                        (2);
+                          case 3:
+                            {
+                              if (offset == 0)
+                                {
+                                      // Set up the right hand side
+                                      // for the horizontal shape
+                                      // functions.
+                                  system_rhs = 0;
 
-                                    for (unsigned int i = 0; i < 2; ++i)
-                                      for (unsigned int j = 0; j <= deg; ++j)
-                                        tmp -= local_dofs[(i + 10) * this->degree
-                                                                   + j]
-                                            * this->shape_value_component
-                                              ((i + 10) * this->degree + j,
-                                               this->generalized_support_points[q_point
-                                                                                + GeometryInfo<dim>::lines_per_cell
-                                                                                * n_edge_points
-                                                                                + 3
-                                                                                * n_face_points],
-                                               2);
-
-                                    for (unsigned i = 0; i <= deg; ++i)
-                                      for (unsigned int j = 0; j < deg; ++j)
-                                        system_rhs (i * deg + j)
-                                          += boundary_weights
-                                             (q_point + n_edge_points,
-                                              2 * (i * deg + j)) * tmp;
-                                  }
+                                  for (unsigned int q_point = 0;
+                                       q_point < n_face_points; ++q_point)
+                                    {
+                                      double tmp
+                                        = values[q_point
+                                                 + GeometryInfo<dim>::lines_per_cell
+                                                 * n_edge_points + 3 * n_face_points]
+                                          (2);
 
-                                system_matrix_inv.vmult (solution, system_rhs);
+                                      for (unsigned int i = 0; i < 2; ++i)
+                                        for (unsigned int j = 0; j < this->degree; ++j)
+                                          tmp -= local_dofs[(i + 10) * this->degree
+                                                                     + j]
+                                              * this->shape_value_component
+                                                ((i + 10) * this->degree + j,
+                                                 this->generalized_support_points[q_point
+                                                                                  + GeometryInfo<dim>::lines_per_cell
+                                                                                  * n_edge_points
+                                                                                  + 3
+                                                                                  * n_face_points],
+                                                 2);
+
+                                      for (unsigned i = 0; i < this->degree; ++i)
+                                        for (unsigned int j = 0; j < deg; ++j)
+                                          system_rhs (i * deg + j)
+                                            += boundary_weights
+                                               (q_point + n_edge_points,
+                                                2 * (i * deg + j)) * tmp;
+                                    }
 
-                                  // Add the computed values
-                                  // to the resulting vector
-                                  // only, if they are not
-                                  // too small.
-                                for (unsigned int i = 0; i <= deg; ++i)
-                                  for (unsigned int j = 0; j < deg; ++j)
-                                    if (std::abs (solution (i * deg + j))
-                                          > 1e-14)
-                                      local_dofs[(i + GeometryInfo<dim>::lines_per_cell
-                                                  + 6 * this->degree) * deg + j
-                                                  + GeometryInfo<dim>::lines_per_cell]
-                                        = solution (i * deg + j);
-                              }
+                                  system_matrix_inv.vmult (solution, system_rhs);
+
+                                      // Add the computed values
+                                      // to the resulting vector
+                                      // only, if they are not
+                                      // too small.
+                                  for (unsigned int i = 0; i < this->degree; ++i)
+                                    for (unsigned int j = 0; j < deg; ++j)
+                                      if (std::abs (solution (i * deg + j))
+                                            > 1e-14)
+                                        local_dofs[(i + GeometryInfo<dim>::lines_per_cell
+                                                      + 6 * this->degree) * deg + j
+                                                      + GeometryInfo<dim>::lines_per_cell]
+                                          = solution (i * deg + j);
+                                }
 
                                   // Set up the right hand side
                                   // for the vertical shape
                                   // functions.
-                            system_rhs = 0;
+                              system_rhs = 0;
 
-                            for (unsigned int q_point = 0;
-                                 q_point < n_face_points; ++q_point)
-                              {
-                                double tmp
-                                  = values[q_point
-                                           + GeometryInfo<dim>::lines_per_cell
-                                           * n_edge_points + 3
-                                           * n_face_points] (0);
-
-                                for (unsigned int i = 0; i < 2; ++i)
-                                  for (unsigned int j = 0; j <= deg; ++j)
-                                    tmp -= local_dofs[(4 * i + 3)
-                                                      * this->degree + j]
-                                        * this->shape_value_component
-                                          ((4 * i + 3) * this->degree
-                                                       + j,
-                                           this->generalized_support_points[q_point
-                                                                            + GeometryInfo<dim>::lines_per_cell
-                                                                            * n_edge_points
-                                                                            + 3
-                                                                            * n_face_points],
-                                           0);
-
-                                 for (unsigned int i = 0; i <= deg; ++i)
-                                   for (unsigned int j = 0; j < deg; ++j)
-                                     system_rhs (i * deg + j)
-                                       += boundary_weights
-                                          (q_point + n_edge_points,
-                                           2 * (i * deg + j) + 1) * tmp;
-                              }
+                              for (unsigned int q_point = 0;
+                                   q_point < n_face_points; ++q_point)
+                                {
+                                  double tmp
+                                    = values[q_point
+                                             + GeometryInfo<dim>::lines_per_cell
+                                             * n_edge_points + 3
+                                             * n_face_points] (0);
+
+                                  for (unsigned int i = 0; i < 2; ++i)
+                                    for (unsigned int j = 0; j < this->degree; ++j)
+                                      tmp -= local_dofs[(4 * i + 3)
+                                                        * this->degree + j]
+                                          * this->shape_value_component
+                                            ((4 * i + 3) * this->degree
+                                                         + j,
+                                             this->generalized_support_points[q_point
+                                                                              + GeometryInfo<dim>::lines_per_cell
+                                                                              * n_edge_points
+                                                                              + 3
+                                                                              * n_face_points],
+                                             0);
+
+                                   for (unsigned int i = 0; i < this->degree; ++i)
+                                     for (unsigned int j = 0; j < deg; ++j)
+                                       system_rhs (i * deg + j)
+                                         += boundary_weights
+                                            (q_point + n_edge_points,
+                                             2 * (i * deg + j) + 1) * tmp;
+                                }
 
-                            system_matrix_inv.vmult (solution, system_rhs);
+                              system_matrix_inv.vmult (solution, system_rhs);
 
                                   // Add the computed values
                                   // to the resulting vector
                                   // only, if they are not
                                   // too small.
-                            for (unsigned int i = 0; i <= deg; ++i)
-                              for (unsigned int j = 0; j < deg; ++j)
-                                if (std::abs (solution (i * deg + j)) > 1e-14)
-                                  local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
-                                                  + 7 * deg) * this->degree]
-                                    = solution (i * deg + j);
-
-                            break;
-                          }
+                              for (unsigned int i = 0; i < this->degree; ++i)
+                                for (unsigned int j = 0; j < deg; ++j)
+                                  if (std::abs (solution (i * deg + j)) > 1e-14)
+                                    local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
+                                                      + 7 * deg) * this->degree]
+                                      = solution (i * deg + j);
 
-                        case 4:
-                          {
+                              break;
+                            }
+
+                          case 4:
+                            {
+                              if (offset + 1 < dim)
+                                {
+                                      // Set up the right hand side
+                                      // for the horizontal shape
+                                      // functions.
+                                  if (offset == 0)
+                                    {
+                                      system_rhs = 0;
+
+                                      for (unsigned int q_point = 0;
+                                           q_point < n_face_points; ++q_point)
+                                        {
+                                          double tmp
+                                            = values[q_point
+                                                     + GeometryInfo<dim>::lines_per_cell
+                                                     * n_edge_points + 4
+                                                     * n_face_points] (0);
+
+                                          for (unsigned int i = 0; i < 2; ++i)
+                                            for (unsigned int j = 0; j < this->degree; ++j)
+                                              tmp -= local_dofs[(i + 2)
+                                                                * this->degree
+                                                                + j]
+                                                     * this->shape_value_component
+                                                       ((i + 2) * this->degree
+                                                        + j,
+                                                        this->generalized_support_points[q_point
+                                                                                         + GeometryInfo<dim>::lines_per_cell
+                                                                                         * n_edge_points
+                                                                                         + 4
+                                                                                         * n_face_points],
+                                                        0);
+
+                                          for (unsigned int i = 0; i < this->degree; ++i)
+                                            for (unsigned int j = 0; j < deg; ++j)
+                                              system_rhs (i * deg + j)
+                                                += boundary_weights
+                                                   (q_point + n_edge_points,
+                                                    2 * (i * deg + j)) * tmp;
+                                        }
+
+                                      system_matrix_inv.vmult
+                                        (solution, system_rhs);
+
+                                          // Add the computed values
+                                          // to the resulting vector
+                                          // only, if they are not
+                                          // too small.
+                                      for (unsigned int i = 0; i < this->degree; ++i)
+                                        for (unsigned int j = 0; j < deg; ++j)
+                                          if (std::abs (solution (i * deg + j))
+                                                > 1e-14)
+                                            local_dofs[(i + GeometryInfo<dim>::lines_per_cell
+                                                          + 8 * this->degree) * deg
+                                                          + j
+                                                          + GeometryInfo<dim>::lines_per_cell]
+                                              = solution (i * deg + j);
+                                    }
+
+                                      // Set up the right hand side
+                                      // for the vertical shape
+                                      // functions.
+                                  system_rhs = 0;
+
+                                  for (unsigned int q_point = 0;
+                                       q_point < n_face_points; ++q_point)
+                                    {
+                                      double tmp
+                                        = values[q_point
+                                                 + GeometryInfo<dim>::lines_per_cell
+                                                 * n_edge_points + 4
+                                                 * n_face_points] (1);
+
+                                      for (unsigned int i = 0; i < 2; ++i)
+                                        for (unsigned int j = 0; j < this->degree; ++j)
+                                          tmp -= local_dofs[i * this->degree + j]
+                                                 * this->shape_value_component
+                                                   (i * this->degree + j,
+                                                    this->generalized_support_points[q_point
+                                                                                     + GeometryInfo<dim>::lines_per_cell
+                                                                                     * n_edge_points
+                                                                                     + 4
+                                                                                     * n_face_points],
+                                                    1);
+
+                                      for (unsigned i = 0; i < this->degree; ++i)
+                                        for (unsigned int j = 0; j < deg; ++j)
+                                          system_rhs (i * deg + j)
+                                            += boundary_weights
+                                               (q_point + n_edge_points,
+                                                2 * (i * deg + j) + 1) * tmp;
+                                    }
+
+                                  system_matrix_inv.vmult (solution, system_rhs);
+
+                                      // Add the computed values
+                                      // to the resulting vector
+                                      // only, if they are not
+                                      // too small.
+                                  for (unsigned int i = 0; i < this->degree; ++i)
+                                    for (unsigned int j = 0; j < deg; ++j)
+                                      if (std::abs (solution (i * deg + j))
+                                            > 1e-14)
+                                        local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
+                                                          + 9 * deg)
+                                                     * this->degree]
+                                          = solution (i * deg + j);
+                                }
+
+                              break;
+                            }
+
+                          default:
                             if (offset + 1 < dim)
                               {
                                   // Set up the right hand side
@@ -3943,25 +2739,23 @@ FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
                                         double tmp
                                           = values[q_point
                                                    + GeometryInfo<dim>::lines_per_cell
-                                                   * n_edge_points + 4
-                                                   * n_face_points] (0);
+                                                   * n_edge_points
+                                                   + 5 * n_face_points] (0);
 
                                         for (unsigned int i = 0; i < 2; ++i)
-                                          for (unsigned int j = 0; j <= deg; ++j)
-                                            tmp -= local_dofs[(i + 2)
-                                                              * this->degree
-                                                              + j]
+                                          for (unsigned int j = 0; j < this->degree; ++j)
+                                            tmp -= local_dofs[(i + 6)
+                                                              * this->degree + j]
                                                    * this->shape_value_component
-                                                     ((i + 2) * this->degree
-                                                      + j,
+                                                     ((i + 6) * this->degree + j,
                                                       this->generalized_support_points[q_point
                                                                                        + GeometryInfo<dim>::lines_per_cell
                                                                                        * n_edge_points
-                                                                                       + 4
+                                                                                       + 5
                                                                                        * n_face_points],
                                                       0);
 
-                                        for (unsigned int i = 0; i <= deg; ++i)
+                                        for (unsigned int i = 0; i < this->degree; ++i)
                                           for (unsigned int j = 0; j < deg; ++j)
                                             system_rhs (i * deg + j)
                                               += boundary_weights
@@ -3972,18 +2766,18 @@ FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
                                     system_matrix_inv.vmult
                                       (solution, system_rhs);
 
-                                  // Add the computed values
-                                  // to the resulting vector
-                                  // only, if they are not
-                                  // too small.
-                                    for (unsigned int i = 0; i <= deg; ++i)
+                                      // Add the computed values
+                                      // to the resulting vector
+                                      // only, if they are not
+                                      // too small.
+                                    for (unsigned int i = 0; i < this->degree; ++i)
                                       for (unsigned int j = 0; j < deg; ++j)
                                         if (std::abs (solution (i * deg + j))
                                               > 1e-14)
                                           local_dofs[(i + GeometryInfo<dim>::lines_per_cell
-                                                      + 8 * this->degree) * deg
-                                                      + j
-                                                      + GeometryInfo<dim>::lines_per_cell]
+                                                        + 10 * this->degree)
+                                                     * deg + j
+                                                     + GeometryInfo<dim>::lines_per_cell]
                                             = solution (i * deg + j);
                                   }
 
@@ -3998,22 +2792,23 @@ FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
                                     double tmp
                                       = values[q_point
                                                + GeometryInfo<dim>::lines_per_cell
-                                               * n_edge_points + 4
+                                               * n_edge_points + 5
                                                * n_face_points] (1);
 
                                     for (unsigned int i = 0; i < 2; ++i)
-                                      for (unsigned int j = 0; j <= deg; ++j)
-                                        tmp -= local_dofs[i * this->degree + j]
+                                      for (unsigned int j = 0; j < this->degree; ++j)
+                                        tmp -= local_dofs[(i + 4)
+                                                          * this->degree + j]
                                                * this->shape_value_component
-                                                 (i * this->degree + j,
+                                                 ((i + 4) * this->degree + j,
                                                   this->generalized_support_points[q_point
                                                                                    + GeometryInfo<dim>::lines_per_cell
                                                                                    * n_edge_points
-                                                                                   + 4
+                                                                                   + 5
                                                                                    * n_face_points],
                                                   1);
 
-                                    for (unsigned i = 0; i <= deg; ++i)
+                                    for (unsigned i = 0; i < this->degree; ++i)
                                       for (unsigned int j = 0; j < deg; ++j)
                                         system_rhs (i * deg + j)
                                           += boundary_weights
@@ -4027,508 +2822,397 @@ FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
                                   // to the resulting vector
                                   // only, if they are not
                                   // too small.
-                                for (unsigned int i = 0; i <= deg; ++i)
+                                for (unsigned int i = 0; i < this->degree; ++i)
                                   for (unsigned int j = 0; j < deg; ++j)
                                     if (std::abs (solution (i * deg + j))
                                           > 1e-14)
                                       local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
-                                                      + 9 * deg)
-                                                 * this->degree]
+                                                        + 11 * deg) * this->degree]
                                         = solution (i * deg + j);
                               }
+                        }
+                    }
 
-                            break;
-                          }
-
-                        default:
-                          if (offset + 1 < dim)
-                            {
-                                  // Set up the right hand side
-                                  // for the horizontal shape
+                                  // Finally we project
+                                  // the remaining parts
+                                  // of the function on
+                                  // the interior shape
                                   // functions.
-                              if (offset == 0)
-                                {
-                                  system_rhs = 0;
-
-                                  for (unsigned int q_point = 0;
-                                       q_point < n_face_points; ++q_point)
-                                    {
-                                      double tmp
-                                        = values[q_point
-                                                 + GeometryInfo<dim>::lines_per_cell
-                                                 * n_edge_points
-                                                 + 5 * n_face_points] (0);
-
-                                      for (unsigned int i = 0; i < 2; ++i)
-                                        for (unsigned int j = 0; j <= deg; ++j)
-                                          tmp -= local_dofs[(i + 6)
-                                                            * this->degree + j]
-                                                 * this->shape_value_component
-                                                   ((i + 6) * this->degree + j,
-                                                    this->generalized_support_points[q_point
-                                                                                     + GeometryInfo<dim>::lines_per_cell
-                                                                                     * n_edge_points
-                                                                                     + 5
-                                                                                     * n_face_points],
-                                                    0);
+                  const QGauss<dim> reference_quadrature (this->degree);
+                  const unsigned int&
+                    n_interior_points = reference_quadrature.size ();
 
-                                      for (unsigned int i = 0; i <= deg; ++i)
-                                        for (unsigned int j = 0; j < deg; ++j)
-                                          system_rhs (i * deg + j)
-                                            += boundary_weights
-                                               (q_point + n_edge_points,
-                                                2 * (i * deg + j)) * tmp;
-                                    }
+                                  // We create the
+                                  // system matrix.
+                  system_matrix.reinit (this->degree * deg * deg,
+                                        this->degree * deg * deg);
+                  system_matrix = 0;
 
-                                  system_matrix_inv.vmult
-                                    (solution, system_rhs);
+                  for (unsigned int i = 0; i < this->degree; ++i)
+                    for (unsigned int j = 0; j < deg; ++j)
+                      for (unsigned int k = 0; k < deg; ++k)
+                        for (unsigned int l = 0; l < this->degree; ++l)
+                          for (unsigned int m = 0; m < deg; ++m)
+                            for (unsigned int n = 0; n < deg; ++n)
+                              for (unsigned int q_point = 0;
+                                   q_point < n_interior_points; ++q_point)
+                                system_matrix ((i * deg + j) * deg + k,
+                                               (l * deg + m) * deg + n)
+                                  += reference_quadrature.weight (q_point)
+                                     * legendre_polynomials[i].value
+                                       (this->generalized_support_points[q_point
+                                                                         + GeometryInfo<dim>::lines_per_cell
+                                                                         * n_edge_points
+                                                                         + GeometryInfo<dim>::faces_per_cell
+                                                                         * n_face_points]
+                                        (0)) * lobatto_polynomials[j + 2].value
+                                               (this->generalized_support_points[q_point
+                                                                                 + GeometryInfo<dim>::lines_per_cell
+                                                                                 * n_edge_points
+                                                                                 + GeometryInfo<dim>::faces_per_cell
+                                                                                 * n_face_points]
+                                                (1))
+                                             * lobatto_polynomials[k + 2].value
+                                               (this->generalized_support_points[q_point
+                                                                                 + GeometryInfo<dim>::lines_per_cell
+                                                                                 * n_edge_points
+                                                                                 + GeometryInfo<dim>::faces_per_cell
+                                                                                 * n_face_points]
+                                                (2))
+                                             * lobatto_polynomials_grad[l].value
+                                               (this->generalized_support_points[q_point
+                                                                                 + GeometryInfo<dim>::lines_per_cell
+                                                                                 * n_edge_points
+                                                                                 + GeometryInfo<dim>::faces_per_cell
+                                                                                 * n_face_points]
+                                                (0))
+                                             * lobatto_polynomials[m + 2].value
+                                               (this->generalized_support_points[q_point
+                                                                                 + GeometryInfo<dim>::lines_per_cell
+                                                                                 * n_edge_points
+                                                                                 + GeometryInfo<dim>::faces_per_cell
+                                                                                 * n_face_points]
+                                                (1))
+                                             * lobatto_polynomials[n + 2].value
+                                               (this->generalized_support_points[q_point
+                                                                                 + GeometryInfo<dim>::lines_per_cell
+                                                                                 * n_edge_points
+                                                                                 + GeometryInfo<dim>::faces_per_cell
+                                                                                 * n_face_points]
+                                                (2));
 
-                                  // Add the computed values
-                                  // to the resulting vector
-                                  // only, if they are not
-                                  // too small.
-                                  for (unsigned int i = 0; i <= deg; ++i)
-                                    for (unsigned int j = 0; j < deg; ++j)
-                                      if (std::abs (solution (i * deg + j))
-                                            > 1e-14)
-                                        local_dofs[(i + GeometryInfo<dim>::lines_per_cell
-                                                    + 10 * this->degree)
-                                                   * deg + j
-                                                   + GeometryInfo<dim>::lines_per_cell]
-                                          = solution (i * deg + j);
-                                }
+                  system_matrix_inv.reinit (system_matrix.m (),
+                                            system_matrix.m ());
+                  system_matrix_inv.invert (system_matrix);
+                  system_rhs.reinit (system_matrix_inv.m ());
+                  solution.reinit (system_matrix.m ());
 
-                                  // Set up the right hand side
-                                  // for the vertical shape
-                                  // functions.
-                              system_rhs = 0;
+                  if (offset + 1 < dim)
+                    {
+                      if (offset == 0)
+                        {
+                                  // Set up the right hand side.
+                          system_rhs = 0;
 
-                              for (unsigned int q_point = 0;
-                                   q_point < n_face_points; ++q_point)
+                          for (unsigned int q_point = 0;
+                               q_point < n_interior_points; ++q_point)
+                            {
+                              double tmp
+                                = values[q_point
+                                         + GeometryInfo<dim>::lines_per_cell
+                                         * n_edge_points
+                                         + GeometryInfo<dim>::faces_per_cell
+                                         * n_face_points] (0);
+
+                              for (unsigned int i = 0; i < this->degree; ++i)
                                 {
-                                  double tmp
-                                    = values[q_point
-                                             + GeometryInfo<dim>::lines_per_cell
-                                             * n_edge_points + 5
-                                             * n_face_points] (1);
-
-                                  for (unsigned int i = 0; i < 2; ++i)
-                                    for (unsigned int j = 0; j <= deg; ++j)
-                                      tmp -= local_dofs[(i + 4)
-                                                        * this->degree + j]
+                                  for (unsigned int j = 0; j < 2; ++j)
+                                    for (unsigned int k = 0; k < 2; ++k)
+                                      tmp -= local_dofs[i + (j + 4 * k + 2)
+                                                          * this->degree]
                                              * this->shape_value_component
-                                               ((i + 4) * this->degree + j,
+                                               (i + (j + 4 * k + 2)
+                                                  * this->degree,
                                                 this->generalized_support_points[q_point
                                                                                  + GeometryInfo<dim>::lines_per_cell
                                                                                  * n_edge_points
-                                                                                 + 5
+                                                                                 + GeometryInfo<dim>::faces_per_cell
                                                                                  * n_face_points],
-                                                1);
+                                                0);
 
-                                  for (unsigned i = 0; i <= deg; ++i)
-                                    for (unsigned int j = 0; j < deg; ++j)
-                                      system_rhs (i * deg + j)
-                                        += boundary_weights
-                                           (q_point + n_edge_points,
-                                            2 * (i * deg + j) + 1) * tmp;
+                                  for (unsigned int j = 0; j < deg; ++j)
+                                    for (unsigned int k = 0; k < 4; ++k)
+                                      tmp -= local_dofs[(i + 2 * (k + 2)
+                                                           * this->degree
+                                                           + GeometryInfo<dim>::lines_per_cell)
+                                                        * deg + j
+                                                        + GeometryInfo<dim>::lines_per_cell]
+                                             * this->shape_value_component
+                                               ((i + 2 * (k + 2) * this->degree
+                                                   + GeometryInfo<dim>::lines_per_cell)
+                                                * deg + j
+                                                + GeometryInfo<dim>::lines_per_cell,
+                                                this->generalized_support_points[q_point
+                                                                                 + GeometryInfo<dim>::lines_per_cell
+                                                                                 * n_edge_points
+                                                                                 + GeometryInfo<dim>::faces_per_cell
+                                                                                 * n_face_points],
+                                                0);
                                 }
 
-                              system_matrix_inv.vmult (solution, system_rhs);
-
-                                  // Add the computed values
-                                  // to the resulting vector
-                                  // only, if they are not
-                                  // too small.
-                              for (unsigned int i = 0; i <= deg; ++i)
+                              for (unsigned int i = 0; i < this->degree; ++i)
                                 for (unsigned int j = 0; j < deg; ++j)
-                                  if (std::abs (solution (i * deg + j))
-                                        > 1e-14)
-                                    local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
-                                                    + 11 * deg) * this->degree]
-                                      = solution (i * deg + j);
-                            }
-                      }
-                  }
-
-                                  // Finally we project
-                                  // the remaining parts
-                                  // of the function on
-                                  // the interior shape
-                                  // functions.
-                const QGauss<dim> reference_quadrature (this->degree);
-                const unsigned int&
-                  n_interior_points = reference_quadrature.size ();
-
-                                  // We create the
-                                  // system matrix.
-                system_matrix.reinit (this->degree * deg * deg,
-                                      this->degree * deg * deg);
-                system_matrix = 0;
-
-                for (unsigned int i = 0; i <= deg; ++i)
-                  for (unsigned int j = 0; j < deg; ++j)
-                    for (unsigned int k = 0; k < deg; ++k)
-                      for (unsigned int l = 0; l <= deg; ++l)
-                        for (unsigned int m = 0; m < deg; ++m)
-                          for (unsigned int n = 0; n < deg; ++n)
-                            for (unsigned int q_point = 0;
-                                 q_point < n_interior_points; ++q_point)
-                              system_matrix ((i * deg + j) * deg + k,
-                                             (l * deg + m) * deg + n)
-                                += reference_quadrature.weight (q_point)
-                                   * legendre_polynomials[i].value
-                                     (this->generalized_support_points[q_point
-                                                                       + GeometryInfo<dim>::lines_per_cell
-                                                                       * n_edge_points
-                                                                       + GeometryInfo<dim>::faces_per_cell
-                                                                       * n_face_points]
-                                      (0)) * lobatto_polynomials[j + 2].value
+                                  for (unsigned int k = 0; k < deg; ++k)
+                                    system_rhs ((i * deg + j) * deg + k)
+                                      += reference_quadrature.weight (q_point)
+                                         * tmp
+                                         * lobatto_polynomials_grad[i].value
+                                           (this->generalized_support_points[q_point
+                                                                             + GeometryInfo<dim>::lines_per_cell
+                                                                             * n_edge_points
+                                                                             + GeometryInfo<dim>::faces_per_cell
+                                                                             * n_face_points]
+                                            (0))
+                                           * lobatto_polynomials[j + 2].value
                                              (this->generalized_support_points[q_point
                                                                                + GeometryInfo<dim>::lines_per_cell
                                                                                * n_edge_points
                                                                                + GeometryInfo<dim>::faces_per_cell
                                                                                * n_face_points]
-                                              (1))
+                                            (1))
                                            * lobatto_polynomials[k + 2].value
                                              (this->generalized_support_points[q_point
                                                                                + GeometryInfo<dim>::lines_per_cell
                                                                                * n_edge_points
                                                                                + GeometryInfo<dim>::faces_per_cell
                                                                                * n_face_points]
-                                              (2))
-                                           * lobatto_polynomials_grad[l].value
-                                             (this->generalized_support_points[q_point
-                                                                               + GeometryInfo<dim>::lines_per_cell
-                                                                               * n_edge_points
-                                                                               + GeometryInfo<dim>::faces_per_cell
-                                                                               * n_face_points]
-                                              (0))
-                                           * lobatto_polynomials[m + 2].value
-                                             (this->generalized_support_points[q_point
-                                                                               + GeometryInfo<dim>::lines_per_cell
-                                                                               * n_edge_points
-                                                                               + GeometryInfo<dim>::faces_per_cell
-                                                                               * n_face_points]
-                                              (1))
-                                           * lobatto_polynomials[n + 2].value
-                                             (this->generalized_support_points[q_point
-                                                                               + GeometryInfo<dim>::lines_per_cell
-                                                                               * n_edge_points
-                                                                               + GeometryInfo<dim>::faces_per_cell
-                                                                               * n_face_points]
-                                              (2));
+                                            (2));
+                            }
 
-                system_matrix_inv.reinit (system_matrix.m (),
-                                          system_matrix.m ());
-                system_matrix_inv.invert (system_matrix);
-                system_rhs.reinit (system_matrix_inv.m ());
-                solution.reinit (system_matrix.m ());
+                          system_matrix_inv.vmult (solution, system_rhs);
+
+                                  // Add the computed values
+                                  // to the resulting vector
+                                  // only, if they are not
+                                  // too small.
+                          for (unsigned int i = 0; i < this->degree; ++i)
+                            for (unsigned int j = 0; j < deg; ++j)
+                              for (unsigned int k = 0; k < deg; ++k)
+                                if (std::abs (solution ((i * deg + j) * deg + k))
+                                      > 1e-14)
+                                  local_dofs[((i + 2
+                                                 * GeometryInfo<dim>::faces_per_cell)
+                                              * deg + j
+                                              + GeometryInfo<dim>::lines_per_cell
+                                              + 2
+                                              * GeometryInfo<dim>::faces_per_cell)
+                                             * deg + k
+                                             + GeometryInfo<dim>::lines_per_cell]
+                                    = solution ((i * deg + j) * deg + k);
+                        }
 
-                if (offset + 1 < dim)
-                  {
-                    if (offset == 0)
-                      {
                                   // Set up the right hand side.
-                        system_rhs = 0;
+                      system_rhs = 0;
 
-                        for (unsigned int q_point = 0;
-                             q_point < n_interior_points; ++q_point)
-                          {
-                            double tmp
-                              = values[q_point
-                                       + GeometryInfo<dim>::lines_per_cell
-                                       * n_edge_points
-                                       + GeometryInfo<dim>::faces_per_cell
-                                       * n_face_points] (0);
-
-                            for (unsigned int i = 0; i <= deg; ++i)
+                      for (unsigned int q_point = 0; q_point < n_interior_points;
+                           ++q_point)
+                        {
+                          double tmp
+                            = values[q_point + GeometryInfo<dim>::lines_per_cell
+                                             * n_edge_points
+                                             + GeometryInfo<dim>::faces_per_cell
+                                             * n_face_points] (1);
+
+                          for (unsigned int i = 0; i < this->degree; ++i)
+                            for (unsigned int j = 0; j < 2; ++j)
                               {
-                                for (unsigned int j = 0; j < 2; ++j)
-                                  for (unsigned int k = 0; k < 2; ++k)
-                                    tmp -= local_dofs[i + (j + 4 * k + 2)
+                                for (unsigned int k = 0; k < 2; ++k)
+                                  tmp -= local_dofs[i + (4 * j + k)
                                                       * this->degree]
-                                           * this->shape_value_component
-                                             (i + (j + 4 * k + 2)
-                                              * this->degree,
-                                              this->generalized_support_points[q_point
-                                                                               + GeometryInfo<dim>::lines_per_cell
-                                                                               * n_edge_points
-                                                                               + GeometryInfo<dim>::faces_per_cell
-                                                                               * n_face_points],
-                                              0);
-
-                                for (unsigned int j = 0; j < deg; ++j)
-                                  for (unsigned int k = 0; k < 4; ++k)
-                                    tmp -= local_dofs[(i + 2 * (k + 2)
-                                                       * this->degree
-                                                       + GeometryInfo<dim>::lines_per_cell)
-                                                      * deg + j
-                                                      + GeometryInfo<dim>::lines_per_cell]
-                                           * this->shape_value_component
-                                             ((i + 2 * (k + 2) * this->degree
-                                               + GeometryInfo<dim>::lines_per_cell)
-                                              * deg + j
-                                              + GeometryInfo<dim>::lines_per_cell,
-                                              this->generalized_support_points[q_point
-                                                                               + GeometryInfo<dim>::lines_per_cell
-                                                                               * n_edge_points
-                                                                               + GeometryInfo<dim>::faces_per_cell
-                                                                               * n_face_points],
-                                              0);
-                              }
+                                         * this->shape_value_component
+                                           (i + (4 * j + k) * this->degree,
+                                            this->generalized_support_points[q_point
+                                                                             + GeometryInfo<dim>::lines_per_cell
+                                                                             * n_edge_points
+                                                                             + GeometryInfo<dim>::faces_per_cell
+                                                                             * n_face_points],
+                                            1);
 
-                            for (unsigned int i = 0; i <= deg; ++i)
-                              for (unsigned int j = 0; j < deg; ++j)
                                 for (unsigned int k = 0; k < deg; ++k)
-                                  system_rhs ((i * deg + j) * deg + k)
-                                    += reference_quadrature.weight (q_point)
-                                       * tmp
-                                       * lobatto_polynomials_grad[i].value
-                                         (this->generalized_support_points[q_point
-                                                                           + GeometryInfo<dim>::lines_per_cell
-                                                                           * n_edge_points
-                                                                           + GeometryInfo<dim>::faces_per_cell
-                                                                           * n_face_points]
-                                          (0))
-                                         * lobatto_polynomials[j + 2].value
-                                           (this->generalized_support_points[q_point
+                                  tmp -= local_dofs[(i + 2 * j * this->degree
+                                                       + GeometryInfo<dim>::lines_per_cell)
+                                                    * deg + k
+                                                    + GeometryInfo<dim>::lines_per_cell]
+                                         * this->shape_value_component
+                                           ((i + 2 * j * this->degree
+                                               + GeometryInfo<dim>::lines_per_cell)
+                                            * deg + k
+                                            + GeometryInfo<dim>::lines_per_cell,
+                                            this->generalized_support_points[q_point
                                                                              + GeometryInfo<dim>::lines_per_cell
                                                                              * n_edge_points
                                                                              + GeometryInfo<dim>::faces_per_cell
-                                                                             * n_face_points]
-                                          (1))
-                                         * lobatto_polynomials[k + 2].value
-                                           (this->generalized_support_points[q_point
+                                                                             * n_face_points],
+                                            1)
+                                         + local_dofs[i + ((2 * j + 9) * deg + k
+                                                           + GeometryInfo<dim>::lines_per_cell)
+                                                        * this->degree]
+                                         * this->shape_value_component
+                                           (i + ((2 * j + 9) * deg + k
+                                                 + GeometryInfo<dim>::lines_per_cell)
+                                              * this->degree,
+                                            this->generalized_support_points[q_point
                                                                              + GeometryInfo<dim>::lines_per_cell
                                                                              * n_edge_points
                                                                              + GeometryInfo<dim>::faces_per_cell
-                                                                             * n_face_points]
-                                          (2));
-                          }
-
-                        system_matrix_inv.vmult (solution, system_rhs);
-
-                                  // Add the computed values
-                                  // to the resulting vector
-                                  // only, if they are not
-                                  // too small.
-                        for (unsigned int i = 0; i <= deg; ++i)
-                          for (unsigned int j = 0; j < deg; ++j)
-                            for (unsigned int k = 0; k < deg; ++k)
-                              if (std::abs (solution ((i * deg + j) * deg + k))
-                                    > 1e-14)
-                                local_dofs[((i + 2
-                                             * GeometryInfo<dim>::faces_per_cell)
-                                            * deg + j
-                                            + GeometryInfo<dim>::lines_per_cell
-                                            + 2
-                                            * GeometryInfo<dim>::faces_per_cell)
-                                           * deg + k
-                                           + GeometryInfo<dim>::lines_per_cell]
-                                = solution ((i * deg + j) * deg + k);
-                      }
-
-                                  // Set up the right hand side.
-                    system_rhs = 0;
-
-                    for (unsigned int q_point = 0; q_point < n_interior_points;
-                         ++q_point)
-                      {
-                        double tmp
-                          = values[q_point + GeometryInfo<dim>::lines_per_cell
-                                   * n_edge_points
-                                   + GeometryInfo<dim>::faces_per_cell
-                                   * n_face_points] (1);
-
-                        for (unsigned int i = 0; i <= deg; ++i)
-                          for (unsigned int j = 0; j < 2; ++j)
-                            {
-                              for (unsigned int k = 0; k < 2; ++k)
-                                tmp -= local_dofs[i + (4 * j + k)
-                                                  * this->degree]
-                                       * this->shape_value_component
-                                         (i + (4 * j + k) * this->degree,
-                                          this->generalized_support_points[q_point
-                                                                           + GeometryInfo<dim>::lines_per_cell
-                                                                           * n_edge_points
-                                                                           + GeometryInfo<dim>::faces_per_cell
-                                                                           * n_face_points],
-                                          1);
+                                                                             * n_face_points],
+                                            1);
+                              }
 
+                          for (unsigned int i = 0; i < this->degree; ++i)
+                            for (unsigned int j = 0; j < deg; ++j)
                               for (unsigned int k = 0; k < deg; ++k)
-                                tmp -= local_dofs[(i + 2 * j * this->degree
-                                                   + GeometryInfo<dim>::lines_per_cell)
-                                                  * deg + k
-                                                  + GeometryInfo<dim>::lines_per_cell]
-                                       * this->shape_value_component
-                                         ((i + 2 * j * this->degree
-                                           + GeometryInfo<dim>::lines_per_cell)
-                                          * deg + k
-                                          + GeometryInfo<dim>::lines_per_cell,
-                                          this->generalized_support_points[q_point
-                                                                           + GeometryInfo<dim>::lines_per_cell
-                                                                           * n_edge_points
-                                                                           + GeometryInfo<dim>::faces_per_cell
-                                                                           * n_face_points],
-                                          1)
-                                       + local_dofs[i + ((2 * j + 9) * deg + k
-                                                    + GeometryInfo<dim>::lines_per_cell)
-                                                    * this->degree]
-                                       * this->shape_value_component
-                                         (i + ((2 * j + 9) * deg + k
-                                          + GeometryInfo<dim>::lines_per_cell)
-                                          * this->degree,
-                                          this->generalized_support_points[q_point
-                                                                           + GeometryInfo<dim>::lines_per_cell
-                                                                           * n_edge_points
-                                                                           + GeometryInfo<dim>::faces_per_cell
-                                                                           * n_face_points],
-                                          1);
-                            }
-
-                        for (unsigned int i = 0; i <= deg; ++i)
-                          for (unsigned int j = 0; j < deg; ++j)
-                            for (unsigned int k = 0; k < deg; ++k)
-                              system_rhs ((i * deg + j) * deg + k)
-                                += reference_quadrature.weight (q_point) * tmp
-                                   * lobatto_polynomials_grad[i].value
-                                     (this->generalized_support_points[q_point
-                                                                       + GeometryInfo<dim>::lines_per_cell
-                                                                       * n_edge_points
-                                                                       + GeometryInfo<dim>::faces_per_cell
-                                                                       * n_face_points]
-                                      (1))
-                                   * lobatto_polynomials[j + 2].value
-                                     (this->generalized_support_points[q_point
-                                                                       + GeometryInfo<dim>::lines_per_cell
-                                                                       * n_edge_points
-                                                                       + GeometryInfo<dim>::faces_per_cell
-                                                                       * n_face_points]
-                                      (0))
-                                   * lobatto_polynomials[k + 2].value
-                                     (this->generalized_support_points[q_point
-                                                                       + GeometryInfo<dim>::lines_per_cell
-                                                                       * n_edge_points
-                                                                       + GeometryInfo<dim>::faces_per_cell
-                                                                       * n_face_points]
-                                      (2));
-                      }
+                                system_rhs ((i * deg + j) * deg + k)
+                                  += reference_quadrature.weight (q_point) * tmp
+                                     * lobatto_polynomials_grad[i].value
+                                       (this->generalized_support_points[q_point
+                                                                         + GeometryInfo<dim>::lines_per_cell
+                                                                         * n_edge_points
+                                                                         + GeometryInfo<dim>::faces_per_cell
+                                                                         * n_face_points]
+                                        (1))
+                                     * lobatto_polynomials[j + 2].value
+                                       (this->generalized_support_points[q_point
+                                                                         + GeometryInfo<dim>::lines_per_cell
+                                                                         * n_edge_points
+                                                                         + GeometryInfo<dim>::faces_per_cell
+                                                                         * n_face_points]
+                                        (0))
+                                     * lobatto_polynomials[k + 2].value
+                                       (this->generalized_support_points[q_point
+                                                                         + GeometryInfo<dim>::lines_per_cell
+                                                                         * n_edge_points
+                                                                         + GeometryInfo<dim>::faces_per_cell
+                                                                         * n_face_points]
+                                        (2));
+                        }
 
-                    system_matrix_inv.vmult (solution, system_rhs);
+                      system_matrix_inv.vmult (solution, system_rhs);
 
                                   // Add the computed values
                                   // to the resulting vector
                                   // only, if they are not
                                   // too small.
-                    for (unsigned int i = 0; i <= deg; ++i)
-                      for (unsigned int j = 0; j < deg; ++j)
-                        for (unsigned int k = 0; k < deg; ++k)
-                          if (std::abs (solution ((i * deg + j) * deg + k))
-                                > 1e-14)
-                            local_dofs[((i + this->degree + 2
-                                         * GeometryInfo<dim>::faces_per_cell)
-                                        * deg + j
-                                        + GeometryInfo<dim>::lines_per_cell + 2
-                                        * GeometryInfo<dim>::faces_per_cell)
-                                       * deg + k
-                                       + GeometryInfo<dim>::lines_per_cell]
-                              = solution ((i * deg + j) * deg + k);
-                  }
+                      for (unsigned int i = 0; i < this->degree; ++i)
+                        for (unsigned int j = 0; j < deg; ++j)
+                          for (unsigned int k = 0; k < deg; ++k)
+                            if (std::abs (solution ((i * deg + j) * deg + k))
+                                  > 1e-14)
+                              local_dofs[((i + this->degree + 2
+                                             * GeometryInfo<dim>::faces_per_cell)
+                                          * deg + j
+                                          + GeometryInfo<dim>::lines_per_cell + 2
+                                          * GeometryInfo<dim>::faces_per_cell)
+                                         * deg + k
+                                         + GeometryInfo<dim>::lines_per_cell]
+                                = solution ((i * deg + j) * deg + k);
+                    }
 
                                   // Set up the right hand side.
-                system_rhs = 0;
-
-                for (unsigned int q_point = 0; q_point < n_interior_points;
-                     ++q_point)
-                  {
-                    double tmp
-                      = values[q_point + GeometryInfo<dim>::lines_per_cell
-                               * n_edge_points
-                               + GeometryInfo<dim>::faces_per_cell
-                               * n_face_points] (2);
+                  system_rhs = 0;
 
-                    for (unsigned int i = 0; i <= deg; ++i)
-                      for (unsigned int j = 0; j < 4; ++j)
-                        {
-                          tmp -= local_dofs[i + (j + 8) * this->degree]
-                                 * this->shape_value_component
-                                   (i + (j + 8) * this->degree,
-                                    this->generalized_support_points[q_point
-                                                                     + GeometryInfo<dim>::lines_per_cell
-                                                                     * n_edge_points
-                                                                     + GeometryInfo<dim>::faces_per_cell
-                                                                     * n_face_points],
-                                    2);
+                  for (unsigned int q_point = 0; q_point < n_interior_points;
+                       ++q_point)
+                    {
+                      double tmp
+                        = values[q_point + GeometryInfo<dim>::lines_per_cell
+                                         * n_edge_points
+                                         + GeometryInfo<dim>::faces_per_cell
+                                         * n_face_points] (2);
 
-                          for (unsigned int k = 0; k < deg; ++k)
-                            tmp -= local_dofs[i + ((2 * j + 1) * deg + k
-                                                   + GeometryInfo<dim>::lines_per_cell)
-                                              * this->degree]
+                      for (unsigned int i = 0; i < this->degree; ++i)
+                        for (unsigned int j = 0; j < 4; ++j)
+                          {
+                            tmp -= local_dofs[i + (j + 8) * this->degree]
                                    * this->shape_value_component
-                                     (i + ((2 * j + 1) * deg + k
-                                           + GeometryInfo<dim>::lines_per_cell)
-                                      * this->degree,
+                                     (i + (j + 8) * this->degree,
                                       this->generalized_support_points[q_point
                                                                        + GeometryInfo<dim>::lines_per_cell
                                                                        * n_edge_points
                                                                        + GeometryInfo<dim>::faces_per_cell
                                                                        * n_face_points],
                                       2);
-                        }
 
-                    for (unsigned int i = 0; i <= deg; ++i)
-                      for (unsigned int j = 0; j < deg; ++j)
-                        for (unsigned int k = 0; k < deg; ++k)
-                          system_rhs ((i * deg + j) * deg + k)
-                            += reference_quadrature.weight (q_point) * tmp
-                               * lobatto_polynomials_grad[i].value
-                                 (this->generalized_support_points[q_point
-                                                                   + GeometryInfo<dim>::lines_per_cell
-                                                                   * n_edge_points
-                                                                   + GeometryInfo<dim>::faces_per_cell
-                                                                   * n_face_points]
-                                  (2))
-                               * lobatto_polynomials[j + 2].value
-                                 (this->generalized_support_points[q_point
-                                                                   + GeometryInfo<dim>::lines_per_cell
-                                                                   * n_edge_points
-                                                                   + GeometryInfo<dim>::faces_per_cell
-                                                                   * n_face_points]
-                                  (0))
-                               * lobatto_polynomials[k + 2].value
-                                 (this->generalized_support_points[q_point
-                                                                   + GeometryInfo<dim>::lines_per_cell
-                                                                   * n_edge_points
-                                                                   + GeometryInfo<dim>::faces_per_cell
-                                                                   * n_face_points]
-                                  (1));
-                  }
+                            for (unsigned int k = 0; k < deg; ++k)
+                              tmp -= local_dofs[i + ((2 * j + 1) * deg + k
+                                                     + GeometryInfo<dim>::lines_per_cell)
+                                                  * this->degree]
+                                     * this->shape_value_component
+                                       (i + ((2 * j + 1) * deg + k
+                                             + GeometryInfo<dim>::lines_per_cell)
+                                          * this->degree,
+                                        this->generalized_support_points[q_point
+                                                                         + GeometryInfo<dim>::lines_per_cell
+                                                                         * n_edge_points
+                                                                         + GeometryInfo<dim>::faces_per_cell
+                                                                         * n_face_points],
+                                        2);
+                          }
+
+                      for (unsigned int i = 0; i < this->degree; ++i)
+                        for (unsigned int j = 0; j < deg; ++j)
+                          for (unsigned int k = 0; k < deg; ++k)
+                            system_rhs ((i * deg + j) * deg + k)
+                              += reference_quadrature.weight (q_point) * tmp
+                                 * lobatto_polynomials_grad[i].value
+                                   (this->generalized_support_points[q_point
+                                                                     + GeometryInfo<dim>::lines_per_cell
+                                                                     * n_edge_points
+                                                                     + GeometryInfo<dim>::faces_per_cell
+                                                                     * n_face_points]
+                                    (2))
+                                 * lobatto_polynomials[j + 2].value
+                                   (this->generalized_support_points[q_point
+                                                                     + GeometryInfo<dim>::lines_per_cell
+                                                                     * n_edge_points
+                                                                     + GeometryInfo<dim>::faces_per_cell
+                                                                     * n_face_points]
+                                    (0))
+                                 * lobatto_polynomials[k + 2].value
+                                   (this->generalized_support_points[q_point
+                                                                     + GeometryInfo<dim>::lines_per_cell
+                                                                     * n_edge_points
+                                                                     + GeometryInfo<dim>::faces_per_cell
+                                                                     * n_face_points]
+                                    (1));
+                    }
 
-                system_matrix_inv.vmult (solution, system_rhs);
+                  system_matrix_inv.vmult (solution, system_rhs);
 
                                   // Add the computed values
                                   // to the resulting vector
                                   // only, if they are not
                                   // too small.
-                for (unsigned int i = 0; i <= deg; ++i)
-                  for (unsigned int j = 0; j < deg; ++j)
-                    for (unsigned int k = 0; k < deg; ++k)
-                      if (std::abs (solution ((i * deg + j) * deg + k))
-                            > 1e-14)
-                        local_dofs[i + ((j + 2
-                                         * (deg + GeometryInfo<dim>::faces_per_cell))
-                                        * deg + k
-                                        + GeometryInfo<dim>::lines_per_cell)
-                                   * this->degree]
-                          = solution ((i * deg + j) * deg + k);
-              }
+                  for (unsigned int i = 0; i < this->degree; ++i)
+                    for (unsigned int j = 0; j < deg; ++j)
+                      for (unsigned int k = 0; k < deg; ++k)
+                        if (std::abs (solution ((i * deg + j) * deg + k))
+                              > 1e-14)
+                          local_dofs[i + ((j + 2
+                                             * (deg + GeometryInfo<dim>::faces_per_cell))
+                                          * deg + k
+                                          + GeometryInfo<dim>::lines_per_cell)
+                                       * this->degree]
+                            = solution ((i * deg + j) * deg + k);
+                }
 
-            break;
-          }
+              break;
+            }
 
-        default:
-          Assert (false, ExcNotImplemented ());
-      }
+          default:
+            Assert (false, ExcNotImplemented ());
+        }
+    }
 }
 
 
@@ -4544,7 +3228,6 @@ FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
                               const VectorSlice<const std::vector<std::vector<double> > >& values)
 const
 {
-  const unsigned int deg = this->degree-1;
   Assert (values.size () == this->n_components (),
           ExcDimensionMismatch (values.size (), this->n_components ()));
   Assert (values[0].size () == this->generalized_support_points.size (),
@@ -4554,53 +3237,57 @@ const
           ExcDimensionMismatch (local_dofs.size (), this->dofs_per_cell));
   std::fill (local_dofs.begin (), local_dofs.end (), 0.0);
 
+  const unsigned int deg = this->degree-1;
+  const std::vector<Polynomials::Polynomial<double> >&
+    legendre_polynomials
+      = Polynomials::Legendre::generate_complete_basis (deg);
+  const QGauss<1> edge_quadrature (this->degree);
+  const std::vector<Point<1> >&
+    edge_quadrature_points = edge_quadrature.get_points ();
+  const unsigned int& n_edge_points = edge_quadrature.size ();
+      
   switch (dim)
     {
       case 2:
         {
                                   // Let us begin with the
-                                  // interpolation part.
-          const QGauss<dim - 1> reference_edge_quadrature (this->degree);
-          const unsigned int&
-            n_edge_points = reference_edge_quadrature.size ();
-
-          for (unsigned int i = 0; i < 2; ++i)
-            for (unsigned int j = 0; j < 2; ++j)
-              {
-                for (unsigned int q_point = 0; q_point < n_edge_points;
-                     ++q_point)
-                  local_dofs[(i + 2 * j) * this->degree]
-                    += reference_edge_quadrature.weight (q_point)
-                       * values[1 - j][q_point + (i + 2 * j) * n_edge_points];
-
+                                  // edge degrees of freedom.
+          for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
+            {
+              const double weight = edge_quadrature.weight (q_point);
+              
+              for (unsigned int i = 0; i < this->degree; ++i)
+                {
+                  const double L_i
+                    = weight * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0));
+                  
+                  for (unsigned int j = 0; j < 2; ++j)
+                    for (unsigned int k = 0; k < 2; ++k)
+                      local_dofs[i + (j + 2 * k) * this->degree]
+                        += L_i * values[1 - k][q_point + (j + 2 * k) * n_edge_points];
+                }
+            }
+            
                                   // Add the computed values
                                   // to the resulting vector
                                   // only, if they are not
                                   // too small.
-                if (std::abs (local_dofs[(i + 2 * j) * this->degree]) < 1e-14)
-                  local_dofs[(i + 2 * j) * this->degree] = 0.0;
-              }
+          for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
+            for (unsigned int i = 0; i < this->degree; ++i)
+              if (std::abs (local_dofs[i + e * this->degree]) < 1e-14)
+                local_dofs[i + e * this->degree] = 0.0;
 
                                   // If the degree is greater
                                   // than 0, then we have still
-                                  // some higher order edge
-                                  // shape functions to
-                                  // consider.
-                                  // Here the projection part
-                                  // starts. The dof values
-                                  // are obtained by solving
-                                  // a linear system of
-                                  // equations.
-          if (this->degree-1 > 1)
+                                  // some interpolations onto
+                                  // the interior shape
+                                  // functions left.
+          if (this->degree > 1)
             {
-                                          // We start with projection
-                                          // on the higher order edge
-                                          // shape function.
               const std::vector<Polynomials::Polynomial<double> >&
                 lobatto_polynomials
                   = Polynomials::Lobatto::generate_complete_basis
                     (this->degree);
-              FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
               std::vector<Polynomials::Polynomial<double> >
                 lobatto_polynomials_grad (this->degree);
 
@@ -4609,91 +3296,26 @@ const
                 lobatto_polynomials_grad[i]
                   = lobatto_polynomials[i + 1].derivative ();
 
-                                  // Set up the system matrix.
-                                  // This can be used for all
-                                  // edges.
-              for (unsigned int i = 0; i < system_matrix.m (); ++i)
-                for (unsigned int j = 0; j < system_matrix.n (); ++j)
-                  for (unsigned int q_point = 0; q_point < n_edge_points;
-                       ++q_point)
-                     system_matrix (i, j)
-                       += boundary_weights (q_point, j)
-                          * lobatto_polynomials_grad[i + 1].value
-                            (this->generalized_face_support_points[q_point]
-                             (1));
-
-              FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
-
-              system_matrix_inv.invert (system_matrix);
-
-              const unsigned int
-                line_coordinate[GeometryInfo<2>::lines_per_cell]
-                  = {1, 1, 0, 0};
-              Vector<double> system_rhs (system_matrix.m ());
-              Vector<double> solution (system_rhs.size ());
-
-              for (unsigned int line = 0;
-                   line < GeometryInfo<dim>::lines_per_cell; ++line)
-                {
-                                  // Set up the right hand side.
-                  system_rhs = 0;
-
-                  for (unsigned int q_point = 0; q_point < n_edge_points;
-                       ++q_point)
-                    {
-                      const double tmp
-                        = values[line_coordinate[line]][line * n_edge_points
-                                                        + q_point]
-                          - local_dofs[line * this->degree]
-                          * this->shape_value_component
-                            (line * this->degree,
-                             this->generalized_support_points[line
-                                                              * n_edge_points
-                                                              + q_point],
-                             line_coordinate[line]);
-
-                      for (unsigned int i = 0; i < system_rhs.size (); ++i)
-                        system_rhs (i) += boundary_weights (q_point, i) * tmp;
-                    }
-
-                  system_matrix_inv.vmult (solution, system_rhs);
-
-                                  // Add the computed values
-                                  // to the resulting vector
-                                  // only, if they are not
-                                  // too small.
-                  for (unsigned int i = 0; i < solution.size (); ++i)
-                    if (std::abs (solution (i)) > 1e-14)
-                      local_dofs[line * this->degree + i + 1] = solution (i);
-                }
-
-                                  // Then we go on to the
+                                  // We set up the system
+                                  // matrix and use it for
+                                  // both, the horizontal
+                                  // and the vertical
                                   // interior shape
-                                  // functions. Again we
-                                  // set up the system
-                                  // matrix and use it
-                                  // for both, the
-                                  // horizontal and the
-                                  // vertical, interior
-                                  // shape functions.
+                                  // functions.
               const QGauss<dim> reference_quadrature (this->degree);
-              const unsigned int&
-                n_interior_points = reference_quadrature.size ();
-              const std::vector<Polynomials::Polynomial<double> >&
-                legendre_polynomials
-                = Polynomials::Legendre::generate_complete_basis (this->degree-1);
+              const unsigned int& n_interior_points
+                = reference_quadrature.size ();
 
-              system_matrix.reinit ((this->degree-1) * this->degree,
-                                    (this->degree-1) * this->degree);
-              system_matrix = 0;
+              FullMatrix<double> system_matrix (deg * this->degree,
+                                                deg * this->degree);
 
               for (unsigned int i = 0; i < this->degree; ++i)
-                for (unsigned int j = 0; j < this->degree-1; ++j)
+                for (unsigned int j = 0; j < deg; ++j)
                   for (unsigned int k = 0; k < this->degree; ++k)
-                    for (unsigned int l = 0; l < this->degree-1; ++l)
+                    for (unsigned int l = 0; l < deg; ++l)
                       for (unsigned int q_point = 0;
                            q_point < n_interior_points; ++q_point)
-                        system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
+                        system_matrix (i * deg + j, k * deg + l)
                           += reference_quadrature.weight (q_point)
                              * legendre_polynomials[i].value
                                (this->generalized_support_points[q_point
@@ -4716,14 +3338,13 @@ const
                                                                  * n_edge_points]
                                 (1));
 
-              system_matrix_inv.reinit (system_matrix.m (),
-                                        system_matrix.m ());
+              FullMatrix<double> system_matrix_inv (system_matrix.m (),
+                                                    system_matrix.m ());
+                
               system_matrix_inv.invert (system_matrix);
-                                  // Set up the right hand side
-                                  // for the horizontal shape
-                                  // functions.
-              system_rhs.reinit (system_matrix_inv.m ());
-              system_rhs = 0;
+
+              Vector<double> solution (system_matrix_inv.m ());
+              Vector<double> system_rhs (system_matrix.m ());
 
               for (unsigned int q_point = 0; q_point < n_interior_points;
                    ++q_point)
@@ -4830,61 +3451,50 @@ const
       case 3:
         {
                                   // Let us begin with the
-                                  // interpolation part.
-          const QGauss<1> reference_edge_quadrature (this->degree);
-          const unsigned int&
-            n_edge_points = reference_edge_quadrature.size ();
-
+                                  // edge degrees of freedom.
           for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
             {
-              for (unsigned int i = 0; i < 4; ++i)
-                local_dofs[(i + 8) * this->degree]
-                  += reference_edge_quadrature.weight (q_point)
-                     * values[2][q_point + (i + 8) * n_edge_points];
-
-              for (unsigned int i = 0; i < 2; ++i)
-                for (unsigned int j = 0; j < 2; ++j)
-                  for (unsigned int k = 0; k < 2; ++k)
-                    local_dofs[(i + 2 * (2 * j + k)) * this->degree]
-                      += reference_edge_quadrature.weight (q_point)
-                         * values[1 - k][q_point + (i + 2 * (2 * j + k))
-                                         * n_edge_points];
+              const double weight = edge_quadrature.weight (q_point);
+              
+              for (unsigned int i = 0; i < this->degree; ++i)
+                {
+                  const double L_i
+                    = weight * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0));
+                    
+                  for (unsigned int j = 0; j < 4; ++j)
+                    local_dofs[i + (j + 8) * this->degree]
+                      += L_i * values[2][q_point + (j + 8) * n_edge_points];
+                  
+                  for (unsigned int j = 0; j < 2; ++j)
+                    for (unsigned int k = 0; k < 2; ++k)
+                      for (unsigned int l = 0; l < 2; ++l)
+                        local_dofs[i + (j + 4 * k + 2 * l) * this->degree]
+                          += L_i * values[1 - l][q_point + (j + 4 * k + 2 * l) * n_edge_points];
+                }
             }
-
+            
                                   // Add the computed values
                                   // to the resulting vector
                                   // only, if they are not
                                   // too small.
-          for (unsigned int i = 0; i < 4; ++i)
-            if (std::abs (local_dofs[(i + 8) * this->degree]) < 1e-14)
-              local_dofs[(i + 8) * this->degree] = 0.0;
-
-          for (unsigned int i = 0; i < 2; ++i)
-            for (unsigned int j = 0; j < 2; ++j)
-              for (unsigned int k = 0; k < 2; ++k)
-                if (std::abs (local_dofs[(i + 2 * (2 * j + k)) * this->degree])
-                      < 1e-14)
-                  local_dofs[(i + 2 * (2 * j + k)) * this->degree] = 0.0;
+          for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
+            for (unsigned int i = 0; i < this->degree; ++i)
+              if (std::abs (local_dofs[i + e * this->degree]) < 1e-14)
+                local_dofs[i + e * this->degree] = 0.0;
 
                                   // If the degree is greater
                                   // than 0, then we have still
-                                  // some higher order shape
-                                  // functions to consider.
-                                  // Here the projection part
-                                  // starts. The dof values
-                                  // are obtained by solving
-                                  // a linear system of
-                                  // equations.
+                                  // some interpolation to the
+                                  // face and interior shape
+                                  // functions left.
           if (this->degree > 1)
             {
-                                          // We start with projection
-                                          // on the higher order edge
-                                          // shape function.
               const std::vector<Polynomials::Polynomial<double> >&
                 lobatto_polynomials
                   = Polynomials::Lobatto::generate_complete_basis
                     (this->degree);
-              FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
+              FullMatrix<double> system_matrix (deg * this->degree,
+                                                deg * this->degree);
               std::vector<Polynomials::Polynomial<double> >
                 lobatto_polynomials_grad (this->degree);
 
@@ -4893,90 +3503,23 @@ const
                 lobatto_polynomials_grad[i]
                   = lobatto_polynomials[i + 1].derivative ();
 
-                                  // Set up the system matrix.
-                                  // This can be used for all
-                                  // edges.
-              for (unsigned int i = 0; i < system_matrix.m (); ++i)
-                for (unsigned int j = 0; j < system_matrix.n (); ++j)
-                  for (unsigned int q_point = 0; q_point < n_edge_points;
-                       ++q_point)
-                    system_matrix (i, j)
-                      += boundary_weights (q_point, j)
-                         * lobatto_polynomials_grad[i + 1].value
-                           (this->generalized_face_support_points[q_point]
-                            (1));
-
-              FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
-
-              system_matrix_inv.invert (system_matrix);
-
-              const unsigned int
-                line_coordinate[GeometryInfo<3>::lines_per_cell]
-                  = {1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
-              Vector<double> system_rhs (system_matrix.m ());
-              Vector<double> solution (system_rhs.size ());
-
-              for (unsigned int line = 0;
-                   line < GeometryInfo<dim>::lines_per_cell; ++line)
-                {
-                  // Set up the right hand side.
-                  system_rhs = 0;
-
-                  for (unsigned int q_point = 0; q_point < this->degree; ++q_point)
-                    {
-                      const double tmp
-                        = values[line_coordinate[line]][line * this->degree
-                                                        + q_point]
-                          - local_dofs[line * this->degree]
-                          * this->shape_value_component
-                            (line * this->degree,
-                             this->generalized_support_points[line
-                                                              * this->degree
-                                                              + q_point],
-                             line_coordinate[line]);
-
-                      for (unsigned int i = 0; i < system_rhs.size (); ++i)
-                           system_rhs (i) += boundary_weights (q_point, i)
-                                             * tmp;
-                    }
-
-                  system_matrix_inv.vmult (solution, system_rhs);
-
-                                  // Add the computed values
-                                  // to the resulting vector
-                                  // only, if they are not
-                                  // too small.
-                  for (unsigned int i = 0; i < solution.size (); ++i)
-                    if (std::abs (solution (i)) > 1e-14)
-                      local_dofs[line * this->degree + i + 1] = solution (i);
-                }
-
-                                  // Then we go on to the
-                                  // face shape functions.
-                                  // Again we set up the
-                                  // system matrix and
-                                  // use it for both, the
-                                  // horizontal and the
-                                  // vertical, shape
+                                  // We set up the system
+                                  // matrix and use it for
+                                  // both, the horizontal
+                                  // and the vertical, shape
                                   // functions.
-              const std::vector<Polynomials::Polynomial<double> >&
-                legendre_polynomials
-                  = Polynomials::Legendre::generate_complete_basis (this->degree-1);
-              const unsigned int n_face_points = n_edge_points * n_edge_points;
-
-              system_matrix.reinit ((this->degree-1) * this->degree,
-                                    (this->degree-1) * this->degree);
-              system_matrix = 0;
+              const unsigned int
+                n_face_points = n_edge_points * n_edge_points;
 
               for (unsigned int i = 0; i < this->degree; ++i)
-                for (unsigned int j = 0; j < this->degree-1; ++j)
+                for (unsigned int j = 0; j < deg; ++j)
                   for (unsigned int k = 0; k < this->degree; ++k)
-                    for (unsigned int l = 0; l < this->degree-1; ++l)
+                    for (unsigned int l = 0; l < deg; ++l)
                       for (unsigned int q_point = 0; q_point < n_face_points;
                            ++q_point)
-                        system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
+                        system_matrix (i * deg + j, k * deg + l)
                           += boundary_weights (q_point + n_edge_points,
-                                               2 * (k * (this->degree-1) + l))
+                                               2 * (k * deg + l))
                              * legendre_polynomials[i].value
                                (this->generalized_face_support_points[q_point
                                                                       + 4
@@ -4988,12 +3531,11 @@ const
                                                                       * n_edge_points]
                                 (1));
 
-              system_matrix_inv.reinit (system_matrix.m (),
-                                        system_matrix.m ());
+              FullMatrix<double> system_matrix_inv (system_matrix.m (),
+                                                    system_matrix.m ());
+              
               system_matrix_inv.invert (system_matrix);
-              solution.reinit (system_matrix.m ());
-              system_rhs.reinit (system_matrix.m ());
-
+              
               const unsigned int
                 face_coordinates[GeometryInfo<3>::faces_per_cell][2]
                   = {{1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}};
@@ -5001,6 +3543,8 @@ const
                 edge_indices[GeometryInfo<3>::faces_per_cell][GeometryInfo<3>::lines_per_face]
                   = {{0, 4, 8, 10}, {1, 5, 9, 11}, {8, 9, 2, 6},
                      {10, 11, 3, 7}, {2, 3, 0, 1}, {6, 7, 4, 5}};
+              Vector<double> solution (system_matrix.m ());
+              Vector<double> system_rhs (system_matrix.m ());
 
               for (unsigned int face = 0;
                    face < GeometryInfo<dim>::faces_per_cell; ++face)

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.