]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Simple style changes: indent according to our usual style, break lines
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 14 Aug 2010 02:44:11 +0000 (02:44 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 14 Aug 2010 02:44:11 +0000 (02:44 +0000)
at 80 characters if possible.

git-svn-id: https://svn.dealii.org/trunk@21654 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/vectors.templates.h

index 0260ff79de0f4a1db9894a590b88e0571257d322..6b1504ac559a67c7e06afad873989f2a4f4e01e1 100644 (file)
@@ -2922,816 +2922,1142 @@ namespace internal
 }
 
 
-namespace internals {
-  namespace VectorTools {
-    
-                                        // This function computes the projection of the
-                                        // boundary function on edges for 3D.
+namespace internals
+{
+  namespace VectorTools
+  {
+
+                                    // This function computes the
+                                    // projection of the boundary
+                                    // function on edges for 3D.
     template<int dim, typename cell_iterator>
     void
-    compute_edge_projection (const cell_iterator& cell, const unsigned int face, const
-      unsigned int line, FEValues<dim>& fe_values, const Quadrature<dim>& quadrature,
-      const Function<dim>& boundary_function, const unsigned int first_vector_component,
-      std::vector<double>& dof_values) {
-       fe_values.reinit (cell);
-       
-                                                // Initialize the required objects.
-       std::vector<Tensor<2, dim> > jacobians = fe_values.get_jacobians ();
-       std::vector<Point<dim> > tangentials (fe_values.n_quadrature_points);
-       std::vector<Point<dim> > quadrature_points = fe_values.get_quadrature_points ();
-       std::vector<Vector<double> > values (fe_values.n_quadrature_points,
-         Vector<double> (dim));
-       
-                                                // Get boundary function values at quadrature points.
-       boundary_function.vector_value_list (quadrature_points, values);
-       quadrature_points = quadrature.get_points ();
-       
-       const unsigned int superdegree = cell->get_fe ().degree;
-       const unsigned int degree = superdegree - 1;
-       Point<dim> shifted_reference_point_1;
-       Point<dim> shifted_reference_point_2;
-       unsigned int edge_coordinate_direction[4];
-       
-                                                // Get coordinate directions of the edges of the face.
-       switch (face) {
-          case 0: case 1: {
-                edge_coordinate_direction[0] = 2;
-                edge_coordinate_direction[1] = 2;
-                edge_coordinate_direction[2] = 1;
-                edge_coordinate_direction[3] = 1;
-             break;
-          }
-          
-          case 2: case 3: {
-                edge_coordinate_direction[0] = 0;
-                edge_coordinate_direction[1] = 0;
-                edge_coordinate_direction[2] = 2;
-                edge_coordinate_direction[3] = 2;
-             break;
-          }
-          
-          default: {
-                edge_coordinate_direction[0] = 1;
-                edge_coordinate_direction[1] = 1;
-                edge_coordinate_direction[2] = 0;
-                edge_coordinate_direction[3] = 0;
-          }
-       }
-       
-                                                // The interpolation for the lowest order edge shape
-                                                // functions is just the mean value of the tangential 
-                                                // components of the boundary function on the edge. 
-       for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
-         ++q_point) {
-                                                // Therefore compute the tangential of the edge at the
-                                                // quadrature point.
-          for (unsigned int d = 0; d < dim; ++d) {
-             shifted_reference_point_1 (d) = quadrature_points[q_point] (d);
-             shifted_reference_point_2 (d) = quadrature_points[q_point] (d);
-          }
-          
-          shifted_reference_point_1 (edge_coordinate_direction[line]) += 1e-13;
-          shifted_reference_point_2 (edge_coordinate_direction[line]) -= 1e-13;
-          tangentials[q_point] = 2e13 *
-            (fe_values.get_mapping ().transform_unit_to_real_cell (cell,
-            shifted_reference_point_1)
-            - fe_values.get_mapping ().transform_unit_to_real_cell (cell,
-            shifted_reference_point_2));
-          tangentials[q_point] /= std::sqrt (tangentials[q_point].square ());
-                                        // Compute the mean value.
-          dof_values[line * superdegree] += fe_values.JxW (q_point)
-            * (values[q_point] (0) * tangentials[q_point] (0) + values[q_point] (1)
-            * tangentials[q_point] (1) + values[q_point] (2) * tangentials[q_point] (2))
-            / (jacobians[q_point][0][edge_coordinate_direction[line]]
-            * jacobians[q_point][0][edge_coordinate_direction[line]]
-            + jacobians[q_point][1][edge_coordinate_direction[line]]
-            * jacobians[q_point][1][edge_coordinate_direction[line]]
-            + jacobians[q_point][2][edge_coordinate_direction[line]]
-            * jacobians[q_point][2][edge_coordinate_direction[line]]);
-       }
-       
-                                                // If there are also higher order shape functions we have
-                                                // still some work left.
-       if (degree > 0) {
-          const FEValuesExtractors::Vector vec (first_vector_component);
-          FullMatrix<double> assembling_matrix (degree, fe_values.n_quadrature_points);
-          Tensor<1, dim> shape_value;
-          Tensor<1, dim> tmp;
-          Vector<double> assembling_vector (fe_values.n_quadrature_points);
-          
-                                        // We set up a linear system of equations to get the values
-                                        // for the remaining degrees of freedom associated with
-                                        // the edge.
-          for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
-            ++q_point) {
-                                        // The right hand side of the corresponding problem is the
-                                        // tangential components of the residual of the boundary
-                                        // function and the interpolated part above.
-             tmp = std::sqrt (fe_values.JxW (q_point)
-               / (jacobians[q_point][0][edge_coordinate_direction[line]]
-               * jacobians[q_point][0][edge_coordinate_direction[line]]
-               + jacobians[q_point][1][edge_coordinate_direction[line]]
-               * jacobians[q_point][1][edge_coordinate_direction[line]]
-               + jacobians[q_point][2][edge_coordinate_direction[line]]
-               * jacobians[q_point][2][edge_coordinate_direction[line]]))
-               * tangentials[q_point];
-             shape_value
-               = fe_values[vec].value (cell->get_fe ().face_to_cell_index (line * superdegree, face),
-               q_point);
-                                // In the weak form the right hand side function is multiplicated
-                                // by the higher order shape functions.
-             assembling_vector (q_point) = (values[q_point] (0)
-               - dof_values[line * superdegree] * shape_value[0]) * tmp[0] + (values[q_point] (1)
-               - dof_values[line * superdegree] * shape_value[1]) * tmp[1] + (values[q_point] (2)
-               - dof_values[line * superdegree] * shape_value[2]) * tmp[2];
-             
-             for (unsigned int i = 0; i < degree; ++i)
-                assembling_matrix (i, q_point)
-                  = fe_values[vec].value (cell->get_fe ().face_to_cell_index (i + line * superdegree + 1, face),
-                  q_point) * tmp;
-          }
-          
-          FullMatrix<double> cell_matrix (degree, degree);
-          
-                                        // Create the system matrix by multiplying the assembling
-                                        // matrix with its transposed.
-          assembling_matrix.mTmult (cell_matrix, assembling_matrix);
-          
-          Vector<double> cell_rhs (degree);
-          
-                                        // Create the system right hand side vector by multiplying
-                                        // the assembling matrix with the assembling vector.
-          assembling_matrix.vmult (cell_rhs, assembling_vector);
-          
-          PreconditionJacobi<FullMatrix<double> > precondition;
-          
-                                        // Use Jacobi preconditioner with the PCG method to solve the
-                                        // problem. 
-          precondition.initialize (cell_matrix);
-          
-          SolverControl solver_control (degree, 1e-15, false, false);
-          SolverCG<> cg (solver_control);
-          Vector<double> solution (degree);
-          
-          cg.solve (cell_matrix, solution, cell_rhs, precondition);
-          
-                                        // Store the computed values.
-          for (unsigned int i = 0; i < degree; ++i)
-             dof_values[i + line * superdegree + 1] = solution (i);                                       
-       }
+    compute_edge_projection (const cell_iterator& cell,
+                            const unsigned int face,
+                            const unsigned int line,
+                            FEValues<dim>& fe_values,
+                            const Quadrature<dim>& quadrature,
+                            const Function<dim>& boundary_function,
+                            const unsigned int first_vector_component,
+                            std::vector<double>& dof_values)
+    {
+      fe_values.reinit (cell);
+
+                                      // Initialize the required
+                                      // objects.
+      std::vector<Tensor<2, dim> > jacobians = fe_values.get_jacobians ();
+      std::vector<Point<dim> > tangentials (fe_values.n_quadrature_points);
+      std::vector<Point<dim> > quadrature_points = fe_values.get_quadrature_points ();
+      std::vector<Vector<double> > values (fe_values.n_quadrature_points,
+                                          Vector<double> (dim));
+
+                                      // Get boundary function values
+                                      // at quadrature points.
+      boundary_function.vector_value_list (quadrature_points, values);
+      quadrature_points = quadrature.get_points ();
+
+      const unsigned int superdegree = cell->get_fe ().degree;
+      const unsigned int degree = superdegree - 1;
+      Point<dim> shifted_reference_point_1;
+      Point<dim> shifted_reference_point_2;
+      unsigned int edge_coordinate_direction[4];
+
+                                      // Get coordinate directions of
+                                      // the edges of the face.
+      switch (face)
+       {
+         case 0:
+         case 1:
+         {
+           edge_coordinate_direction[0] = 2;
+           edge_coordinate_direction[1] = 2;
+           edge_coordinate_direction[2] = 1;
+           edge_coordinate_direction[3] = 1;
+           break;
+         }
+
+         case 2:
+         case 3:
+         {
+           edge_coordinate_direction[0] = 0;
+           edge_coordinate_direction[1] = 0;
+           edge_coordinate_direction[2] = 2;
+           edge_coordinate_direction[3] = 2;
+           break;
+         }
+
+         default:
+         {
+           edge_coordinate_direction[0] = 1;
+           edge_coordinate_direction[1] = 1;
+           edge_coordinate_direction[2] = 0;
+           edge_coordinate_direction[3] = 0;
+         }
+       }
+
+                                      // The interpolation for the
+                                      // lowest order edge shape
+                                      // functions is just the mean
+                                      // value of the tangential
+                                      // components of the boundary
+                                      // function on the edge.
+      for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
+          ++q_point)
+       {
+                                          // Therefore compute the
+                                          // tangential of the edge at
+                                          // the quadrature point.
+         for (unsigned int d = 0; d < dim; ++d)
+           {
+             shifted_reference_point_1 (d) = quadrature_points[q_point] (d);
+             shifted_reference_point_2 (d) = quadrature_points[q_point] (d);
+           }
+
+         shifted_reference_point_1 (edge_coordinate_direction[line]) += 1e-13;
+         shifted_reference_point_2 (edge_coordinate_direction[line]) -= 1e-13;
+         tangentials[q_point]
+           = (2e13 *
+              (fe_values.get_mapping ()
+               .transform_unit_to_real_cell (cell,
+                                             shifted_reference_point_1)
+               -
+               fe_values.get_mapping ()
+               .transform_unit_to_real_cell (cell,
+                                             shifted_reference_point_2)));
+         tangentials[q_point]
+           /= std::sqrt (tangentials[q_point].square ());
+
+                                          // Compute the mean value.
+         dof_values[line * superdegree]
+           += (fe_values.JxW (q_point)
+               * (values[q_point] (0) * tangentials[q_point] (0)
+                  + values[q_point] (1) * tangentials[q_point] (1)
+                  + values[q_point] (2) * tangentials[q_point] (2))
+               / (jacobians[q_point][0][edge_coordinate_direction[line]]
+                  * jacobians[q_point][0][edge_coordinate_direction[line]]
+                  +
+                  jacobians[q_point][1][edge_coordinate_direction[line]]
+                  * jacobians[q_point][1][edge_coordinate_direction[line]]
+                  +
+                  jacobians[q_point][2][edge_coordinate_direction[line]]
+                  * jacobians[q_point][2][edge_coordinate_direction[line]]));
+       }
+
+                                      // If there are also higher
+                                      // order shape functions we
+                                      // have still some work left.
+      if (degree > 0)
+       {
+         const FEValuesExtractors::Vector vec (first_vector_component);
+         FullMatrix<double> assembling_matrix (degree, fe_values.n_quadrature_points);
+         Tensor<1, dim> shape_value;
+         Tensor<1, dim> tmp;
+         Vector<double> assembling_vector (fe_values.n_quadrature_points);
+
+                                          // We set up a linear system
+                                          // of equations to get the
+                                          // values for the remaining
+                                          // degrees of freedom
+                                          // associated with the edge.
+         for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
+              ++q_point)
+           {
+                                              // The right hand side of
+                                              // the corresponding
+                                              // problem is the
+                                              // tangential components of
+                                              // the residual of the
+                                              // boundary function and
+                                              // the interpolated part
+                                              // above.
+             tmp =
+               std::sqrt (fe_values.JxW (q_point)
+                          / (jacobians[q_point][0][edge_coordinate_direction[line]]
+                             * jacobians[q_point][0][edge_coordinate_direction[line]]
+                             +
+                             jacobians[q_point][1][edge_coordinate_direction[line]]
+                             * jacobians[q_point][1][edge_coordinate_direction[line]]
+                             +
+                             jacobians[q_point][2][edge_coordinate_direction[line]]
+                             * jacobians[q_point][2][edge_coordinate_direction[line]]))
+               * tangentials[q_point];
+             shape_value
+               = fe_values[vec].value (cell->get_fe ()
+                                       .face_to_cell_index (line * superdegree, face),
+                                       q_point);
+                                              // In the weak form the
+                                              // right hand side function
+                                              // is multiplicated by the
+                                              // higher order shape
+                                              // functions.
+             assembling_vector (q_point)
+               = ((values[q_point] (0)
+                   -
+                   dof_values[line * superdegree] * shape_value[0]) * tmp[0]
+                  +
+                  (values[q_point] (1)
+                   -
+                   dof_values[line * superdegree] * shape_value[1]) * tmp[1]
+                  +
+                  (values[q_point] (2)
+                   -
+                   dof_values[line * superdegree] * shape_value[2]) * tmp[2]);
+
+             for (unsigned int i = 0; i < degree; ++i)
+               assembling_matrix (i, q_point)
+                 = fe_values[vec].value (cell->get_fe ()
+                                         .face_to_cell_index (i + line * superdegree + 1,
+                                                              face),
+                                         q_point) * tmp;
+           }
+
+         FullMatrix<double> cell_matrix (degree, degree);
+
+                                          // Create the system matrix
+                                          // by multiplying the
+                                          // assembling matrix with its
+                                          // transposed.
+         assembling_matrix.mTmult (cell_matrix, assembling_matrix);
+
+         Vector<double> cell_rhs (degree);
+
+                                          // Create the system right
+                                          // hand side vector by
+                                          // multiplying the assembling
+                                          // matrix with the assembling
+                                          // vector.
+         assembling_matrix.vmult (cell_rhs, assembling_vector);
+
+         PreconditionJacobi<FullMatrix<double> > precondition;
+
+                                          // Use Jacobi preconditioner
+                                          // with the PCG method to
+                                          // solve the problem.
+         precondition.initialize (cell_matrix);
+
+         SolverControl solver_control (degree, 1e-15, false, false);
+         SolverCG<> cg (solver_control);
+         Vector<double> solution (degree);
+
+         cg.solve (cell_matrix, solution, cell_rhs, precondition);
+
+                                          // Store the computed values.
+         for (unsigned int i = 0; i < degree; ++i)
+           dof_values[i + line * superdegree + 1] = solution (i);
+       }
     }
-    
-                                        // This function computes the projection of the
-                                        // boundary function on the interior of faces in
-                                        // 3D.
+
+
+
+                                    // This function computes the
+                                    // projection of the boundary
+                                    // function on the interior of
+                                    // faces in 3D.
     template<int dim, typename cell_iterator>
     void
-    compute_face_projection (const cell_iterator& cell, const unsigned int face,
-      FEValues<dim>& fe_values, const Function<dim>& boundary_function, const unsigned
-      int first_vector_component,
-      std::vector<double>& dof_values) {
-       fe_values.reinit (cell);
-       
-                                                // Initialize the required objects.
-       std::vector<Tensor<2, dim> > jacobians = fe_values.get_jacobians ();
-       std::vector<Point<dim> > quadrature_points = fe_values.get_quadrature_points ();
-       std::vector<Vector<double> > values (fe_values.n_quadrature_points,
-         Vector<double> (dim));
-       
-                                                // Get boundary function values at quadrature points.
-       boundary_function.vector_value_list (quadrature_points, values);
-       
-       const FEValuesExtractors::Vector vec (first_vector_component);
-       const unsigned int superdegree = cell->get_fe ().degree;
-       const unsigned int degree = superdegree - 1;
-       double JxW;
-       FullMatrix<double> assembling_matrix (degree * superdegree,
-         dim * fe_values.n_quadrature_points);
-       Vector<double> assembling_vector (assembling_matrix.n ());
-       Vector<double> cell_rhs (assembling_matrix.m ());
-       FullMatrix<double> cell_matrix (assembling_matrix.m (), assembling_matrix.m ());
-       Vector<double> solution (cell_matrix.m ());
-       SolverControl solver_control (cell_matrix.m (), 1e-15, false, false);
-       SolverCG<> cg (solver_control);
-       PreconditionJacobi<FullMatrix<double> > precondition;
-       Tensor<1, dim> tmp;
-       Tensor<1, dim> shape_value;
-       unsigned int global_face_coordinate_directions[2];
-       unsigned int local_face_coordinate_directions[2];
-       
-                                                // Get coordinate directions of the face.
-       switch (face) {
-          case 0: case 1: {
-             global_face_coordinate_directions[0] = 1;
-             global_face_coordinate_directions[1] = 2;
-                local_face_coordinate_directions[0] = 1;
-                local_face_coordinate_directions[1] = 0;
-             break;
-          }
-          
-          case 2: case 3: {
-             global_face_coordinate_directions[0] = 0;
-             global_face_coordinate_directions[1] = 2;
-                local_face_coordinate_directions[0] = 0;
-                local_face_coordinate_directions[1] = 1;
-             break;
-          }
-          
-          default: {
-             global_face_coordinate_directions[0] = 0;
-             global_face_coordinate_directions[1] = 1;
-                local_face_coordinate_directions[0] = 1;
-                local_face_coordinate_directions[1] = 0;
-          }
-       }
-       
-                                        // The projection is divided into two steps. In the first step we 
-                                        // project the boundary function on the horizontal shape functions.
-                                        // Then the bounary function is projected on the vertical shape
-                                        // functions.
-                                        // We begin with the horizontal shape functions and set up a linear
-                                        // system of equations to get the values for degrees of freedom
-                                        // associated with the interior of the face.
-       for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
-         ++q_point) {
-                                        // The right hand side of the corresponding problem is the
-                                        // residual of the boundary function and the already
-                                        // interpolated part on the edges.
-          for (unsigned int d = 0; d < dim; ++d)
-             tmp[d] = values[q_point] (d);
-          
-          for (unsigned int i = 0; i < 2; ++i)
-             for (unsigned int j = 0; j <= degree; ++j)
-                tmp -= dof_values[(i + 2 * local_face_coordinate_directions[0]) * superdegree + j]
-                  * fe_values[vec].value (cell->get_fe ().face_to_cell_index
-                  ((i + 2 * local_face_coordinate_directions[0]) * superdegree + j,
-                  face), q_point);
-          
-          JxW = std::sqrt (fe_values.JxW (q_point)
-            / ((jacobians[q_point][0][global_face_coordinate_directions[0]]
-            * jacobians[q_point][0][global_face_coordinate_directions[0]]
-            + jacobians[q_point][1][global_face_coordinate_directions[0]]
-            * jacobians[q_point][1][global_face_coordinate_directions[0]]
-            + jacobians[q_point][2][global_face_coordinate_directions[0]]
-            * jacobians[q_point][2][global_face_coordinate_directions[0]])
-            * (jacobians[q_point][0][global_face_coordinate_directions[1]]
-            * jacobians[q_point][0][global_face_coordinate_directions[1]]
-            + jacobians[q_point][1][global_face_coordinate_directions[1]]
-            * jacobians[q_point][1][global_face_coordinate_directions[1]]
-            + jacobians[q_point][2][global_face_coordinate_directions[1]]
-            * jacobians[q_point][2][global_face_coordinate_directions[1]])));
-          
-                                // In the weak form the right hand side function is multiplicated
-                                // by the horizontal shape functions defined in the interior of the
-                                // face.          
-          for (unsigned int d = 0; d < dim; ++d)
-             assembling_vector (dim * q_point + d) = JxW * tmp[d];
-          
-          for (unsigned int i = 0; i <= degree; ++i)
-             for (unsigned int j = 0; j < degree; ++j) {
-                shape_value = JxW
-                  * fe_values[vec].value (cell->get_fe ().face_to_cell_index
-                  ((i + GeometryInfo<dim>::lines_per_face) * degree + j
-                    + GeometryInfo<dim>::lines_per_face, face), q_point);
-                
-                for (unsigned int d = 0; d < dim; ++d)
-                   assembling_matrix (i * degree + j, dim * q_point + d) = shape_value[d];
-             }
-       }
-       
-                                        // Create the system matrix by multiplying the assembling
-                                        // matrix with its transposed and the right hand side vector
-                                        // by mutliplying the assembling matrix with the assembling
-                                        // vector. The problem is solved by the PCG method.      
-       assembling_matrix.mTmult (cell_matrix, assembling_matrix);
-       assembling_matrix.vmult (cell_rhs, assembling_vector);
-       precondition.initialize (cell_matrix);
-       cg.solve (cell_matrix, solution, cell_rhs, precondition);
-       
-                                        // Store the computed values.
-       for (unsigned int i = 0; i <= degree; ++i)
-          for (unsigned int j = 0; j < degree; ++j)
-             dof_values[(i + GeometryInfo<dim>::lines_per_face) * degree + j
-               + GeometryInfo<dim>::lines_per_face] = solution (i * degree + j);
-       
-                                                // Now we do the same as above with the vertical shape functions
-                                                // instead of the horizontal ones.
-       for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
-         ++q_point) {
-          for (unsigned int d = 0; d < dim; ++d)
-             tmp[d] = values[q_point] (d);
-          
-          for (unsigned int i = 0; i < 2; ++i)
-             for (unsigned int j = 0; j <= degree; ++j)
-                tmp
-                  -= dof_values[(i + 2 * local_face_coordinate_directions[1]) * superdegree + j]
-                  * fe_values[vec].value (cell->get_fe ().face_to_cell_index
-                  ((i + 2 * local_face_coordinate_directions[1]) * superdegree + j,
-                  face), q_point);
-          
-          JxW = std::sqrt (fe_values.JxW (q_point)
-            / ((jacobians[q_point][0][global_face_coordinate_directions[0]]
-            * jacobians[q_point][0][global_face_coordinate_directions[0]]
-            + jacobians[q_point][1][global_face_coordinate_directions[0]]
-            * jacobians[q_point][1][global_face_coordinate_directions[0]]
-            + jacobians[q_point][2][global_face_coordinate_directions[0]]
-            * jacobians[q_point][2][global_face_coordinate_directions[0]])
-            * (jacobians[q_point][0][global_face_coordinate_directions[1]]
-            * jacobians[q_point][0][global_face_coordinate_directions[1]]
-            + jacobians[q_point][1][global_face_coordinate_directions[1]]
-            * jacobians[q_point][1][global_face_coordinate_directions[1]]
-            + jacobians[q_point][2][global_face_coordinate_directions[1]]
-            * jacobians[q_point][2][global_face_coordinate_directions[1]])));
-          
-          for (unsigned int d = 0; d < dim; ++d)
-             assembling_vector (dim * q_point + d) = JxW * tmp[d];
-          
-          for (unsigned int i = 0; i < degree; ++i)
-             for (unsigned int j = 0; j <= degree; ++j) {
-                shape_value = JxW
-                  * fe_values[vec].value (cell->get_fe ().face_to_cell_index
-                  ((i + degree + GeometryInfo<dim>::lines_per_face) * superdegree + j,
-                  face), q_point);
-                
-                for (unsigned int d = 0; d < dim; ++d)
-                   assembling_matrix (i * superdegree + j, dim * q_point + d) = shape_value[d];
-             }
-       }
-       
-       assembling_matrix.mTmult (cell_matrix, assembling_matrix);
-       assembling_matrix.vmult (cell_rhs, assembling_vector);
-       precondition.initialize (cell_matrix);
-       cg.solve (cell_matrix, solution, cell_rhs, precondition);
-       
-       for (unsigned int i = 0; i < degree; ++i)
-          for (unsigned int j = 0; j <= degree; ++j)
-             dof_values[(i + degree + GeometryInfo<dim>::lines_per_face) * superdegree + j]
-               = solution (i * superdegree + j);
+    compute_face_projection (const cell_iterator& cell,
+                            const unsigned int face,
+                            FEValues<dim>& fe_values,
+                            const Function<dim>& boundary_function,
+                            const unsigned int first_vector_component,
+                            std::vector<double>& dof_values)
+    {
+      fe_values.reinit (cell);
+
+                                      // Initialize the required objects.
+      std::vector<Tensor<2, dim> > jacobians = fe_values.get_jacobians ();
+      std::vector<Point<dim> > quadrature_points = fe_values.get_quadrature_points ();
+      std::vector<Vector<double> > values (fe_values.n_quadrature_points,
+                                          Vector<double> (dim));
+
+                                      // Get boundary function values
+                                      // at quadrature points.
+      boundary_function.vector_value_list (quadrature_points, values);
+
+      const FEValuesExtractors::Vector vec (first_vector_component);
+      const unsigned int superdegree = cell->get_fe ().degree;
+      const unsigned int degree = superdegree - 1;
+      double JxW;
+      FullMatrix<double> assembling_matrix (degree * superdegree,
+                                           dim * fe_values.n_quadrature_points);
+      Vector<double> assembling_vector (assembling_matrix.n ());
+      Vector<double> cell_rhs (assembling_matrix.m ());
+      FullMatrix<double> cell_matrix (assembling_matrix.m (),
+                                     assembling_matrix.m ());
+      Vector<double> solution (cell_matrix.m ());
+      SolverControl solver_control (cell_matrix.m (), 1e-15, false, false);
+      SolverCG<> cg (solver_control);
+      PreconditionJacobi<FullMatrix<double> > precondition;
+      Tensor<1, dim> tmp;
+      Tensor<1, dim> shape_value;
+      unsigned int global_face_coordinate_directions[2];
+      unsigned int local_face_coordinate_directions[2];
+
+                                      // Get coordinate directions of
+                                      // the face.
+      switch (face)
+       {
+         case 0:
+         case 1:
+         {
+           global_face_coordinate_directions[0] = 1;
+           global_face_coordinate_directions[1] = 2;
+           local_face_coordinate_directions[0] = 1;
+           local_face_coordinate_directions[1] = 0;
+           break;
+         }
+
+         case 2:
+         case 3:
+         {
+           global_face_coordinate_directions[0] = 0;
+           global_face_coordinate_directions[1] = 2;
+           local_face_coordinate_directions[0] = 0;
+           local_face_coordinate_directions[1] = 1;
+           break;
+         }
+
+         default:
+         {
+           global_face_coordinate_directions[0] = 0;
+           global_face_coordinate_directions[1] = 1;
+           local_face_coordinate_directions[0] = 1;
+           local_face_coordinate_directions[1] = 0;
+         }
+       }
+
+                                      // The projection is divided
+                                      // into two steps. In the first
+                                      // step we project the boundary
+                                      // function on the horizontal
+                                      // shape functions.  Then the
+                                      // bounary function is
+                                      // projected on the vertical
+                                      // shape functions.  We begin
+                                      // with the horizontal shape
+                                      // functions and set up a
+                                      // linear system of equations
+                                      // to get the values for
+                                      // degrees of freedom
+                                      // associated with the interior
+                                      // of the face.
+      for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
+          ++q_point)
+       {
+                                          // The right hand side of the
+                                          // corresponding problem is
+                                          // the residual of the
+                                          // boundary function and the
+                                          // already interpolated part
+                                          // on the edges.
+         for (unsigned int d = 0; d < dim; ++d)
+           tmp[d] = values[q_point] (d);
+
+         for (unsigned int i = 0; i < 2; ++i)
+           for (unsigned int j = 0; j <= degree; ++j)
+             tmp -= dof_values[(i + 2 * local_face_coordinate_directions[0]) * superdegree + j]
+                    * fe_values[vec].value (cell->get_fe ().face_to_cell_index
+                                            ((i + 2 * local_face_coordinate_directions[0]) * superdegree + j,
+                                             face), q_point);
+
+         JxW = std::sqrt (fe_values.JxW (q_point)
+                          / ((jacobians[q_point][0][global_face_coordinate_directions[0]]
+                              * jacobians[q_point][0][global_face_coordinate_directions[0]]
+                              +
+                              jacobians[q_point][1][global_face_coordinate_directions[0]]
+                              * jacobians[q_point][1][global_face_coordinate_directions[0]]
+                              +
+                              jacobians[q_point][2][global_face_coordinate_directions[0]]
+                              * jacobians[q_point][2][global_face_coordinate_directions[0]])
+                             *
+                             (jacobians[q_point][0][global_face_coordinate_directions[1]]
+                              * jacobians[q_point][0][global_face_coordinate_directions[1]]
+                              +
+                              jacobians[q_point][1][global_face_coordinate_directions[1]]
+                              * jacobians[q_point][1][global_face_coordinate_directions[1]]
+                              +
+                              jacobians[q_point][2][global_face_coordinate_directions[1]]
+                              * jacobians[q_point][2][global_face_coordinate_directions[1]])));
+
+                                          // In the weak form the right
+                                          // hand side function is
+                                          // multiplicated by the
+                                          // horizontal shape functions
+                                          // defined in the interior of
+                                          // the face.
+         for (unsigned int d = 0; d < dim; ++d)
+           assembling_vector (dim * q_point + d) = JxW * tmp[d];
+
+         for (unsigned int i = 0; i <= degree; ++i)
+           for (unsigned int j = 0; j < degree; ++j)
+             {
+               shape_value
+                 = (JxW
+                    * fe_values[vec].value (cell->get_fe ().face_to_cell_index
+                                            ((i + GeometryInfo<dim>::lines_per_face)
+                                             * degree
+                                             + j
+                                             + GeometryInfo<dim>::lines_per_face,
+                                             face),
+                                            q_point));
+
+               for (unsigned int d = 0; d < dim; ++d)
+                 assembling_matrix (i * degree + j,
+                                    dim * q_point + d)
+                   = shape_value[d];
+             }
+       }
+
+                                      // Create the system matrix by
+                                      // multiplying the assembling
+                                      // matrix with its transposed
+                                      // and the right hand side
+                                      // vector by mutliplying the
+                                      // assembling matrix with the
+                                      // assembling vector. The
+                                      // problem is solved by the PCG
+                                      // method.
+      assembling_matrix.mTmult (cell_matrix, assembling_matrix);
+      assembling_matrix.vmult (cell_rhs, assembling_vector);
+      precondition.initialize (cell_matrix);
+      cg.solve (cell_matrix, solution, cell_rhs, precondition);
+
+                                      // Store the computed values.
+      for (unsigned int i = 0; i <= degree; ++i)
+       for (unsigned int j = 0; j < degree; ++j)
+         dof_values[(i + GeometryInfo<dim>::lines_per_face) * degree + j
+                    + GeometryInfo<dim>::lines_per_face]
+           = solution (i * degree + j);
+
+                                      // Now we do the same as above
+                                      // with the vertical shape
+                                      // functions instead of the
+                                      // horizontal ones.
+      for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
+          ++q_point)
+       {
+         for (unsigned int d = 0; d < dim; ++d)
+           tmp[d] = values[q_point] (d);
+
+         for (unsigned int i = 0; i < 2; ++i)
+           for (unsigned int j = 0; j <= degree; ++j)
+             tmp
+               -= dof_values[(i + 2 * local_face_coordinate_directions[1]) * superdegree + j]
+               * fe_values[vec].value (cell->get_fe ().face_to_cell_index
+                                       ((i + 2 * local_face_coordinate_directions[1]) * superdegree + j,
+                                        face), q_point);
+
+         JxW = std::sqrt (fe_values.JxW (q_point)
+                          / ((jacobians[q_point][0][global_face_coordinate_directions[0]]
+                              * jacobians[q_point][0][global_face_coordinate_directions[0]]
+                              +
+                              jacobians[q_point][1][global_face_coordinate_directions[0]]
+                              * jacobians[q_point][1][global_face_coordinate_directions[0]]
+                              +
+                              jacobians[q_point][2][global_face_coordinate_directions[0]]
+                              * jacobians[q_point][2][global_face_coordinate_directions[0]])
+                             *
+                             (jacobians[q_point][0][global_face_coordinate_directions[1]]
+                              * jacobians[q_point][0][global_face_coordinate_directions[1]]
+                              +
+                              jacobians[q_point][1][global_face_coordinate_directions[1]]
+                              * jacobians[q_point][1][global_face_coordinate_directions[1]]
+                              +
+                              jacobians[q_point][2][global_face_coordinate_directions[1]]
+                              * jacobians[q_point][2][global_face_coordinate_directions[1]])));
+
+         for (unsigned int d = 0; d < dim; ++d)
+           assembling_vector (dim * q_point + d) = JxW * tmp[d];
+
+         for (unsigned int i = 0; i < degree; ++i)
+           for (unsigned int j = 0; j <= degree; ++j)
+             {
+               shape_value = JxW
+                             * fe_values[vec].value (cell->get_fe ().face_to_cell_index
+                                                     ((i + degree + GeometryInfo<dim>::lines_per_face) * superdegree + j,
+                                                      face), q_point);
+
+               for (unsigned int d = 0; d < dim; ++d)
+                 assembling_matrix (i * superdegree + j, dim * q_point + d)
+                   = shape_value[d];
+             }
+       }
+
+      assembling_matrix.mTmult (cell_matrix, assembling_matrix);
+      assembling_matrix.vmult (cell_rhs, assembling_vector);
+      precondition.initialize (cell_matrix);
+      cg.solve (cell_matrix, solution, cell_rhs, precondition);
+
+      for (unsigned int i = 0; i < degree; ++i)
+       for (unsigned int j = 0; j <= degree; ++j)
+         dof_values[(i + degree + GeometryInfo<dim>::lines_per_face) * superdegree + j]
+           = solution (i * superdegree + j);
     }
-    
-    
-                                        // This function computes the projection of the
-                                        // boundary function on the faces in 2D.  
+
+
+
+
+                                    // This function computes the
+                                    // projection of the boundary
+                                    // function on the faces in 2D.
     template<int dim, typename cell_iterator>
     void
-    compute_face_projection (const cell_iterator& cell, const unsigned int face,
-      FEValues<dim>& fe_values, const Quadrature<dim>& quadrature, const Function<dim>&
-        boundary_function, const unsigned int first_vector_component, std::vector<double>&
-        dof_values) {
-       fe_values.reinit (cell);
-       
-                                                // Initialize the required objects.
-       std::vector<Point<dim> > tangentials (fe_values.n_quadrature_points);
-       std::vector<Tensor<2, dim> > jacobians = fe_values.get_jacobians ();
-       std::vector<Point<dim> > quadrature_points = fe_values.get_quadrature_points ();
-       std::vector<Vector<double> > values (fe_values.n_quadrature_points,
-         Vector<double> (dim));
-       
-                                                // Get boundary function values at quadrature points.
-       boundary_function.vector_value_list (quadrature_points, values);
-       quadrature_points = quadrature.get_points ();
-       
-       const unsigned int degree = cell->get_fe ().degree - 1;
-       Point<dim> shifted_reference_point_1;
-       Point<dim> shifted_reference_point_2;
-       Tensor<1, dim> tmp;
-       Tensor<1, dim> shape_value;
-       unsigned int face_coordinate_direction;
-       
-                                                // Get coordinate directions of the face.
-       switch (face) {
-          case 0: case 1: {
-                face_coordinate_direction = 1;
-             break;
-          }
-          
-          default:
-                face_coordinate_direction = 0;
-       }
-       
-                                                // The interpolation for the lowest order face shape
-                                                // functions is just the mean value of the tangential 
-                                                // components of the boundary function on the edge. 
-       for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
-         ++q_point) {
-                                                // Therefore compute the tangential of the face at the
-                                                // quadrature point.
-          for (unsigned int d = 0; d < dim; ++d) {
-             shifted_reference_point_1 (d) = quadrature_points[q_point] (d);
-             shifted_reference_point_2 (d) = quadrature_points[q_point] (d);
-          }
-          
-          shifted_reference_point_1 (face_coordinate_direction) += 1e-13;
-          shifted_reference_point_2 (face_coordinate_direction) -= 1e-13;
-          tangentials[q_point] = 2e13
-            * (fe_values.get_mapping ().transform_unit_to_real_cell (cell, shifted_reference_point_1)
-            - fe_values.get_mapping ().transform_unit_to_real_cell (cell, shifted_reference_point_2));
-          tangentials[q_point] /= std::sqrt (tangentials[q_point].square ());
-                                        // Compute the mean value.
-          dof_values[0] += fe_values.JxW (q_point) * (values[q_point] (0)
-            * tangentials[q_point] (0) + values[q_point] (1) * tangentials[q_point] (1))
-            / (jacobians[q_point][0][face_coordinate_direction]
-            * jacobians[q_point][0][face_coordinate_direction]
-            + jacobians[q_point][1][face_coordinate_direction]
-            * jacobians[q_point][1][face_coordinate_direction]);
-       }
-       
-                                                // If there are also higher order shape functions we have
-                                                // still some work left.
-       if (degree > 0) {
-          const FEValuesExtractors::Vector vec (first_vector_component);
-          FullMatrix<double> assembling_matrix (degree, fe_values.n_quadrature_points);
-          Vector<double> assembling_vector (fe_values.n_quadrature_points);
-          
-                                        // We set up a linear system of equations to get the values
-                                        // for the remaining degrees of freedom associated with
-                                        // the face.
-          for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
-            ++q_point) {
-                                        // The right hand side of the corresponding problem is the
-                                        // tangential components of the residual of the boundary
-                                        // function and the interpolated part above.
-             tmp = std::sqrt (fe_values.JxW (q_point)
-               / std::sqrt (jacobians[q_point][0][face_coordinate_direction]
-               * jacobians[q_point][0][face_coordinate_direction]
-               + jacobians[q_point][1][face_coordinate_direction]
-               * jacobians[q_point][1][face_coordinate_direction])) * tangentials[q_point];
-             shape_value
-               = fe_values[vec].value (cell->get_fe ().face_to_cell_index (0, face), q_point);    
-             assembling_vector (q_point) = (values[q_point] (0)
-               - dof_values[0] * shape_value[0]) * tmp[0] + (values[q_point] (1)
-               - dof_values[1] * shape_value[1]) * tmp[1];
-             
-                                // In the weak form the right hand side function is multiplicated
-                                // by the higher order shape functions.
-             for (unsigned int i = 0; i < degree; ++i)
-                assembling_matrix (i, q_point)
-                  = fe_values[vec].value (cell->get_fe ().face_to_cell_index (i + 1, face),
-                  q_point) * tmp;
-          }
-          
-          FullMatrix<double> cell_matrix (degree, degree);
-          
-                                        // Create the system matrix by multiplying the assembling
-                                        // matrix with its transposed.
-          assembling_matrix.mTmult (cell_matrix, assembling_matrix);
-          
-          Vector<double> cell_rhs (degree);
-          
-                                        // Create the system right hand side vector by multiplying
-                                        // the assembling matrix with the assembling vector.
-          assembling_matrix.vmult (cell_rhs, assembling_vector);
-          
-          PreconditionJacobi<FullMatrix<double> > precondition;
-          
-                                        // Use Jacobi preconditioner with the PCG method to solve the
-                                        // problem. 
-          precondition.initialize (cell_matrix);
-          
-          SolverControl solver_control (degree, 1e-15, false, false);
-          SolverCG<> cg (solver_control);
-          Vector<double> solution (degree);
-          
-          cg.solve (cell_matrix, solution, cell_rhs, precondition);
-          
-                                        // Store the computed values.
-          for (unsigned int i = 0; i < degree; ++i)
-             dof_values[i + 1] = solution (i);
-       }
+    compute_face_projection (const cell_iterator& cell,
+                            const unsigned int face,
+                            FEValues<dim>& fe_values,
+                            const Quadrature<dim>& quadrature,
+                            const Function<dim>& boundary_function,
+                            const unsigned int first_vector_component,
+                            std::vector<double>& dof_values)
+    {
+      fe_values.reinit (cell);
+
+                                      // Initialize the required objects.
+      std::vector<Point<dim> > tangentials (fe_values.n_quadrature_points);
+      std::vector<Tensor<2, dim> > jacobians = fe_values.get_jacobians ();
+      std::vector<Point<dim> > quadrature_points = fe_values.get_quadrature_points ();
+      std::vector<Vector<double> > values (fe_values.n_quadrature_points,
+                                          Vector<double> (dim));
+
+                                      // Get boundary function values
+                                      // at quadrature points.
+      boundary_function.vector_value_list (quadrature_points, values);
+      quadrature_points = quadrature.get_points ();
+
+      const unsigned int degree = cell->get_fe ().degree - 1;
+      Point<dim> shifted_reference_point_1;
+      Point<dim> shifted_reference_point_2;
+      Tensor<1, dim> tmp;
+      Tensor<1, dim> shape_value;
+      unsigned int face_coordinate_direction;
+
+                                      // Get coordinate directions of the face.
+      switch (face)
+       {
+         case 0:
+         case 1:
+         {
+           face_coordinate_direction = 1;
+           break;
+         }
+
+         default:
+               face_coordinate_direction = 0;
+       }
+
+                                      // The interpolation for the
+                                      // lowest order face shape
+                                      // functions is just the mean
+                                      // value of the tangential
+                                      // components of the boundary
+                                      // function on the edge.
+      for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
+          ++q_point)
+       {
+                                          // Therefore compute the
+                                          // tangential of the face at
+                                          // the quadrature point.
+         for (unsigned int d = 0; d < dim; ++d)
+           {
+             shifted_reference_point_1 (d) = quadrature_points[q_point] (d);
+             shifted_reference_point_2 (d) = quadrature_points[q_point] (d);
+           }
+
+         shifted_reference_point_1 (face_coordinate_direction) += 1e-13;
+         shifted_reference_point_2 (face_coordinate_direction) -= 1e-13;
+         tangentials[q_point] = 2e13
+                                * (fe_values.get_mapping ()
+                                   .transform_unit_to_real_cell (cell,
+                                                                 shifted_reference_point_1)
+                                   -
+                                   fe_values.get_mapping ()
+                                   .transform_unit_to_real_cell (cell,
+                                                                 shifted_reference_point_2));
+         tangentials[q_point] /= std::sqrt (tangentials[q_point].square ());
+                                          // Compute the mean value.
+         dof_values[0] += fe_values.JxW (q_point)
+                          * (values[q_point] (0)
+                             * tangentials[q_point] (0)
+                             +
+                             values[q_point] (1) * tangentials[q_point] (1))
+                          / (jacobians[q_point][0][face_coordinate_direction]
+                             * jacobians[q_point][0][face_coordinate_direction]
+                             + jacobians[q_point][1][face_coordinate_direction]
+                             * jacobians[q_point][1][face_coordinate_direction]);
+       }
+
+                                      // If there are also higher
+                                      // order shape functions we
+                                      // have still some work left.
+      if (degree > 0)
+       {
+         const FEValuesExtractors::Vector vec (first_vector_component);
+         FullMatrix<double> assembling_matrix (degree, fe_values.n_quadrature_points);
+         Vector<double> assembling_vector (fe_values.n_quadrature_points);
+
+                                          // We set up a linear system
+                                          // of equations to get the
+                                          // values for the remaining
+                                          // degrees of freedom
+                                          // associated with the face.
+         for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
+              ++q_point)
+           {
+                                              // The right hand side of
+                                              // the corresponding
+                                              // problem is the
+                                              // tangential components of
+                                              // the residual of the
+                                              // boundary function and
+                                              // the interpolated part
+                                              // above.
+             tmp = std::sqrt (fe_values.JxW (q_point)
+                              / std::sqrt (jacobians[q_point][0][face_coordinate_direction]
+                                           * jacobians[q_point][0][face_coordinate_direction]
+                                           + jacobians[q_point][1][face_coordinate_direction]
+                                           * jacobians[q_point][1][face_coordinate_direction])) * tangentials[q_point];
+             shape_value
+               = fe_values[vec].value (cell->get_fe ().face_to_cell_index (0, face), q_point);
+             assembling_vector (q_point) = (values[q_point] (0)
+                                            - dof_values[0] * shape_value[0]) * tmp[0] + (values[q_point] (1)
+                                                                                          - dof_values[1] * shape_value[1]) * tmp[1];
+
+                                              // In the weak form the
+                                              // right hand side function
+                                              // is multiplicated by the
+                                              // higher order shape
+                                              // functions.
+             for (unsigned int i = 0; i < degree; ++i)
+               assembling_matrix (i, q_point)
+                 = fe_values[vec].value (cell->get_fe ()
+                                         .face_to_cell_index (i + 1, face),
+                                         q_point) * tmp;
+           }
+
+         FullMatrix<double> cell_matrix (degree, degree);
+
+                                          // Create the system matrix
+                                          // by multiplying the
+                                          // assembling matrix with its
+                                          // transposed.
+         assembling_matrix.mTmult (cell_matrix, assembling_matrix);
+
+         Vector<double> cell_rhs (degree);
+
+                                          // Create the system right
+                                          // hand side vector by
+                                          // multiplying the assembling
+                                          // matrix with the assembling
+                                          // vector.
+         assembling_matrix.vmult (cell_rhs, assembling_vector);
+
+         PreconditionJacobi<FullMatrix<double> > precondition;
+
+                                          // Use Jacobi preconditioner
+                                          // with the PCG method to
+                                          // solve the problem.
+         precondition.initialize (cell_matrix);
+
+         SolverControl solver_control (degree, 1e-15, false, false);
+         SolverCG<> cg (solver_control);
+         Vector<double> solution (degree);
+
+         cg.solve (cell_matrix, solution, cell_rhs, precondition);
+
+                                          // Store the computed values.
+         for (unsigned int i = 0; i < degree; ++i)
+           dof_values[i + 1] = solution (i);
+       }
     }
   }
 }
-    
-    
-                                                // Projection-based interpolation is performed in two (in 2D)
-                                                // respectively three (in 3D) steps. First the tangential
-                                                // component of the function is interpolated on each edge.
-                                                // This gives the values for the degrees of freedom corresponding
-                                                // to the lowest order edge shape functions. Then the interpolated
-                                                // part of the function is subtracted and we project the tangential
-                                                // component of the residual onto the space of the remaining
-                                                // (higher order) edge shape functions. This is done by building
-                                                // a linear system of equations of dimension <tt>degree</tt>. The
-                                                // solution gives us the values for the degrees of freedom
-                                                // corresponding to the remaining edge shape functions. Now we are
-                                                // done for 2D, but in 3D we possibly have also degrees of freedom,
-                                                // which are located in the interior of the faces. Therefore we
-                                                // compute the residual of the function describing the boundary
-                                                // values and the interpolated part, which we have computed in the
-                                                // last two steps. On the faces there are two kinds of shape
-                                                // functions, the horizontal and the vertical ones. Thus we have
-                                                // two solve two linear systems of equations of size
-                                                // <tt>degree * (degree + 1)<tt> to obtain the values for the 
-                                                // corresponding degrees of freedom.
+
+
 
 template <int dim>
-void VectorTools::project_boundary_values_curl_conforming (const DoFHandler<dim>& dof_handler,
-     const unsigned int first_vector_component,
-     const Function<dim>& boundary_function,
-     const unsigned char boundary_component,
-     ConstraintMatrix& constraints,
-     const Mapping<dim>& mapping)
+void
+VectorTools::
+project_boundary_values_curl_conforming (const DoFHandler<dim>& dof_handler,
+                                        const unsigned int first_vector_component,
+                                        const Function<dim>& boundary_function,
+                                        const unsigned char boundary_component,
+                                        ConstraintMatrix& constraints,
+                                        const Mapping<dim>& mapping)
 {
-   std::vector<double> dof_values;
-   std::vector<unsigned int> face_dof_indices;
-   typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active ();
-   unsigned int dofs_per_face;
-   unsigned int superdegree;
-   
-   switch (dim) {
-      case 2: {
-         for (; cell != dof_handler.end (); ++cell)
-            if (cell->at_boundary ())
-               for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
-                  if (cell->face (face)->boundary_indicator () == boundary_component) {
-                     // this is only implemented, if the FE is a Nedelec element
-                     typedef FiniteElement<dim> FEL;
-
-                     AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0),
-                       typename FEL::ExcInterpolationNotImplemented ());
-                     
-                     dofs_per_face = cell->get_fe ().dofs_per_face;
-                     dof_values.resize (dofs_per_face);
-                     
-                     for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
-                        dof_values[dof] = 0.0;
-                     
-                     superdegree = cell->get_fe ().degree;
-                     
-                     QGauss<dim - 1> reference_face_quadrature (2 * superdegree);
-                     Quadrature<dim> face_quadrature
-                       = QProjector<dim>::project_to_face (reference_face_quadrature, face);
-                     FEValues<dim> fe_face_values (mapping, cell->get_fe (), face_quadrature,
-                       update_jacobians | update_JxW_values | update_quadrature_points | update_values);
-                     
-                     // Compute the projection of the boundary function on the edge.
-                     internals::VectorTools::compute_face_projection (cell, face, fe_face_values,
-                       face_quadrature, boundary_function, first_vector_component, dof_values); 
-                     face_dof_indices.resize (dofs_per_face);
-                     cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ());
-                     
-                     // Add the computed constraints to the constraint matrix.
-                     for (unsigned int dof = 0; dof < dofs_per_face; ++dof) {
-                        constraints.add_line (face_dof_indices[dof]);
-                        
-                        if (std::abs (dof_values[dof]) > 1e-14)
-                           constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]);
-                     }
-                  }
-         
-         break;
+                                  // Projection-based interpolation
+                                  // is performed in two (in 2D)
+                                  // respectively three (in 3D)
+                                  // steps. First the tangential
+                                  // component of the function is
+                                  // interpolated on each edge.  This
+                                  // gives the values for the degrees
+                                  // of freedom corresponding to the
+                                  // lowest order edge shape
+                                  // functions. Then the interpolated
+                                  // part of the function is
+                                  // subtracted and we project the
+                                  // tangential component of the
+                                  // residual onto the space of the
+                                  // remaining (higher order) edge
+                                  // shape functions. This is done by
+                                  // building a linear system of
+                                  // equations of dimension
+                                  // <tt>degree</tt>. The solution
+                                  // gives us the values for the
+                                  // degrees of freedom corresponding
+                                  // to the remaining edge shape
+                                  // functions. Now we are done for
+                                  // 2D, but in 3D we possibly have
+                                  // also degrees of freedom, which
+                                  // are located in the interior of
+                                  // the faces. Therefore we compute
+                                  // the residual of the function
+                                  // describing the boundary values
+                                  // and the interpolated part, which
+                                  // we have computed in the last two
+                                  // steps. On the faces there are
+                                  // two kinds of shape functions,
+                                  // the horizontal and the vertical
+                                  // ones. Thus we have two solve two
+                                  // linear systems of equations of
+                                  // size <tt>degree * (degree +
+                                  // 1)<tt> to obtain the values for
+                                  // the corresponding degrees of
+                                  // freedom.
+  std::vector<double> dof_values;
+  std::vector<unsigned int> face_dof_indices;
+  typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active ();
+  unsigned int dofs_per_face;
+  unsigned int superdegree;
+
+  switch (dim)
+    {
+      case 2:
+      {
+       for (; cell != dof_handler.end (); ++cell)
+         if (cell->at_boundary ())
+           for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+             if (cell->face (face)->boundary_indicator () == boundary_component)
+               {
+                                                  // this is only
+                                                  // implemented, if the
+                                                  // FE is a Nedelec
+                                                  // element
+                 typedef FiniteElement<dim> FEL;
+
+                 AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0),
+                              typename FEL::ExcInterpolationNotImplemented ());
+
+                 dofs_per_face = cell->get_fe ().dofs_per_face;
+                 dof_values.resize (dofs_per_face);
+
+                 for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+                   dof_values[dof] = 0.0;
+
+                 superdegree = cell->get_fe ().degree;
+
+                 QGauss<dim - 1> reference_face_quadrature (2 * superdegree);
+                 Quadrature<dim> face_quadrature
+                   = QProjector<dim>::project_to_face (reference_face_quadrature, face);
+                 FEValues<dim> fe_face_values (mapping, cell->get_fe (),
+                                               face_quadrature,
+                                               update_jacobians |
+                                               update_JxW_values |
+                                               update_quadrature_points |
+                                               update_values);
+
+                                                  // Compute the
+                                                  // projection of the
+                                                  // boundary function on
+                                                  // the edge.
+                 internals::VectorTools
+                   ::compute_face_projection (cell, face, fe_face_values,
+                                              face_quadrature, boundary_function,
+                                              first_vector_component, dof_values);
+                 face_dof_indices.resize (dofs_per_face);
+                 cell->face (face)->get_dof_indices (face_dof_indices,
+                                                     cell->active_fe_index ());
+
+                                                  // Add the computed
+                                                  // constraints to the
+                                                  // constraint matrix.
+                 for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+                   {
+                     constraints.add_line (face_dof_indices[dof]);
+
+                     if (std::abs (dof_values[dof]) > 1e-14)
+                       constraints.set_inhomogeneity (face_dof_indices[dof],
+                                                      dof_values[dof]);
+                   }
+               }
+
+       break;
       }
-      
-      case 3: {
-        const unsigned int n_dofs = dof_handler.n_dofs ();
-         std::vector<double> computed_constraints (n_dofs);
-        std::vector<int> projected_dofs (n_dofs);
-        unsigned int degree;
-        unsigned int superdegree;
-        
-        for (unsigned int dof = 0; dof < n_dofs; ++dof)
-           projected_dofs[dof] = -1;
-        
-         for (; cell != dof_handler.end (); ++cell)
-            if (cell->at_boundary ())
-               for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
-                  if (cell->face (face)->boundary_indicator () == boundary_component) {
-                     // this is only implemented, if the FE is a Nedelec element
-                     typedef FiniteElement<dim> FEL;
-
-                     AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0),
-                       typename FEL::ExcInterpolationNotImplemented ());
-                     
-                     superdegree = cell->get_fe ().degree;
-                     degree = superdegree - 1;
-                     
-                     QGauss<dim - 2> reference_edge_quadrature (2 * superdegree);
-                     
-                     dofs_per_face = cell->get_fe ().dofs_per_face;
-                     dof_values.resize (dofs_per_face);
-                     
-                     for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
-                        dof_values[dof] = 0.0;
-                     
-                     face_dof_indices.resize (dofs_per_face);
-                     cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ());
-                     
-                     // First we compute the projection on the edges.
-                     for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) {
-                     // If we have reached this edge through another cell before, we do
-                     // not do here anything unless we have a good reason, i.e. a higher
-                     // polynomial degree.
-                        if (projected_dofs[face_dof_indices[line * superdegree]] < (int) degree) {
-                           Quadrature<dim> edge_quadrature
-                             = QProjector<dim>::project_to_face (QProjector<dim - 1>::project_to_face
-                             (reference_edge_quadrature, line), face);
-                           FEValues<dim> fe_edge_values (mapping, cell->get_fe (), edge_quadrature,
-                             update_JxW_values | update_jacobians | update_quadrature_points | update_values);
-                     // Compute the projection of the boundary function on the edge.
-                           internals::VectorTools::compute_edge_projection (cell, face, line,
-                             fe_edge_values, edge_quadrature, boundary_function, first_vector_component,
-                             dof_values);
-                     // Mark the projected degrees of freedom.
-                           for (unsigned int dof = line * superdegree; dof < (line + 1) * superdegree;
-                             ++dof)
-                              projected_dofs[face_dof_indices[dof]] = degree;
-                        }
-                        
-                     // If we have computed the values in a previous step of the loop,
-                     // we just copy the values in the local vector.
-                        else
-                           for (unsigned int dof = line * superdegree; dof < (line + 1) * superdegree;
-                             ++dof)
-                              dof_values[dof] = computed_constraints[face_dof_indices[dof]];
-                     }
-                     
-                     // If there are higher order shape functions, there is still some
-                     // work left.
-                     if (degree > 0) {
-                        QGauss<dim - 1> reference_face_quadrature (2 * superdegree);
-                        Quadrature<dim> face_quadrature
-                          = QProjector<dim>::project_to_face (reference_face_quadrature, face);
-                        FEValues<dim> fe_face_values (mapping, cell->get_fe (), face_quadrature,
-                          update_JxW_values | update_jacobians | update_quadrature_points | update_values);
-                        
-                     // Compute the projection of the boundary function on the interior
-                     // of the face.
-                        internals::VectorTools::compute_face_projection (cell, face, fe_face_values,
-                          boundary_function, first_vector_component, dof_values);
-                        
-                     // Mark the projected degrees of freedom.
-                        for (unsigned int dof = GeometryInfo<dim>::lines_per_face * superdegree;
-                          dof < dofs_per_face; ++dof)
-                           projected_dofs[face_dof_indices[dof]] = degree;
-                     }
-                     
-                     // Store the computed values in the global vector.
-                     for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
-                        if (std::abs (dof_values[dof]) > 1e-14)
-                           computed_constraints[face_dof_indices[dof]] = dof_values[dof];
-                  }
-         
-                     // Add the computed constraints to the constraint matrix.
-         for (unsigned int dof = 0; dof < n_dofs; ++dof)
-            if (projected_dofs[dof] != -1) {
-               constraints.add_line (dof);
-               constraints.set_inhomogeneity (dof, computed_constraints[dof]);
-            }
+
+      case 3:
+      {
+       const unsigned int n_dofs = dof_handler.n_dofs ();
+       std::vector<double> computed_constraints (n_dofs);
+       std::vector<int> projected_dofs (n_dofs);
+       unsigned int degree;
+       unsigned int superdegree;
+
+       for (unsigned int dof = 0; dof < n_dofs; ++dof)
+         projected_dofs[dof] = -1;
+
+       for (; cell != dof_handler.end (); ++cell)
+         if (cell->at_boundary ())
+           for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+             if (cell->face (face)->boundary_indicator () == boundary_component)
+               {
+                                                  // this is only
+                                                  // implemented, if the
+                                                  // FE is a Nedelec
+                                                  // element
+                 typedef FiniteElement<dim> FEL;
+
+                 AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0),
+                              typename FEL::ExcInterpolationNotImplemented ());
+
+                 superdegree = cell->get_fe ().degree;
+                 degree = superdegree - 1;
+
+                 QGauss<dim - 2> reference_edge_quadrature (2 * superdegree);
+
+                 dofs_per_face = cell->get_fe ().dofs_per_face;
+                 dof_values.resize (dofs_per_face);
+
+                 for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+                   dof_values[dof] = 0.0;
+
+                 face_dof_indices.resize (dofs_per_face);
+                 cell->face (face)->get_dof_indices (face_dof_indices,
+                                                     cell->active_fe_index ());
+
+                                                  // First we compute the
+                                                  // projection on the
+                                                  // edges.
+                 for (unsigned int line = 0;
+                      line < GeometryInfo<3>::lines_per_face; ++line)
+                   {
+                                                      // If we have reached
+                                                      // this edge through
+                                                      // another cell
+                                                      // before, we do not
+                                                      // do here anything
+                                                      // unless we have a
+                                                      // good reason,
+                                                      // i.e. a higher
+                                                      // polynomial degree.
+                     if (projected_dofs[face_dof_indices[line * superdegree]]
+                         <
+                         (int) degree)
+                       {
+                         Quadrature<dim> edge_quadrature
+                           = QProjector<dim>::project_to_face (QProjector<dim - 1>::project_to_face
+                                                               (reference_edge_quadrature, line), face);
+                         FEValues<dim> fe_edge_values (mapping, cell->get_fe (),
+                                                       edge_quadrature,
+                                                       update_JxW_values |
+                                                       update_jacobians |
+                                                       update_quadrature_points |
+                                                       update_values);
+                                                          // Compute the
+                                                          // projection of
+                                                          // the boundary
+                                                          // function on the
+                                                          // edge.
+                         internals::VectorTools
+                           ::compute_edge_projection (cell, face, line,
+                                                      fe_edge_values,
+                                                      edge_quadrature,
+                                                      boundary_function,
+                                                      first_vector_component,
+                                                      dof_values);
+                                                          // Mark the
+                                                          // projected
+                                                          // degrees of
+                                                          // freedom.
+                         for (unsigned int dof = line * superdegree;
+                              dof < (line + 1) * superdegree; ++dof)
+                           projected_dofs[face_dof_indices[dof]] = degree;
+                       }
+
+                                                      // If we have
+                                                      // computed the
+                                                      // values in a
+                                                      // previous step of
+                                                      // the loop, we just
+                                                      // copy the values in
+                                                      // the local vector.
+                     else
+                       for (unsigned int dof = line * superdegree;
+                            dof < (line + 1) * superdegree;
+                            ++dof)
+                         dof_values[dof] = computed_constraints[face_dof_indices[dof]];
+                   }
+
+                                                  // If there are higher
+                                                  // order shape
+                                                  // functions, there is
+                                                  // still some work
+                                                  // left.
+                 if (degree > 0)
+                   {
+                     QGauss<dim - 1> reference_face_quadrature (2 * superdegree);
+                     Quadrature<dim> face_quadrature
+                       = QProjector<dim>::project_to_face (reference_face_quadrature,
+                                                           face);
+                     FEValues<dim> fe_face_values (mapping, cell->get_fe (),
+                                                   face_quadrature,
+                                                   update_JxW_values |
+                                                   update_jacobians |
+                                                   update_quadrature_points |
+                                                   update_values);
+
+                                                      // Compute the
+                                                      // projection of the
+                                                      // boundary function
+                                                      // on the interior of
+                                                      // the face.
+                     internals::VectorTools
+                       ::compute_face_projection (cell, face, fe_face_values,
+                                                  boundary_function,
+                                                  first_vector_component,
+                                                  dof_values);
+
+                                                      // Mark the projected
+                                                      // degrees of
+                                                      // freedom.
+                     for (unsigned int dof = GeometryInfo<dim>::lines_per_face * superdegree;
+                          dof < dofs_per_face; ++dof)
+                       projected_dofs[face_dof_indices[dof]] = degree;
+                   }
+
+                                                  // Store the computed
+                                                  // values in the global
+                                                  // vector.
+                 for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+                   if (std::abs (dof_values[dof]) > 1e-14)
+                     computed_constraints[face_dof_indices[dof]] = dof_values[dof];
+               }
+
+                                        // Add the computed constraints
+                                        // to the constraint matrix.
+       for (unsigned int dof = 0; dof < n_dofs; ++dof)
+         if (projected_dofs[dof] != -1)
+           {
+             constraints.add_line (dof);
+             constraints.set_inhomogeneity (dof, computed_constraints[dof]);
+           }
       }
-   }
+    }
 }
 
 
+
 template <int dim>
-void VectorTools::project_boundary_values_curl_conforming (const hp::DoFHandler<dim>& dof_handler,
-     const unsigned int first_vector_component,
-     const Function<dim>& boundary_function,
-     const unsigned char boundary_component,
-     ConstraintMatrix& constraints,
-     const hp::MappingCollection<dim>& mapping_collection)
+void
+VectorTools::
+project_boundary_values_curl_conforming (const hp::DoFHandler<dim>& dof_handler,
+                                        const unsigned int first_vector_component,
+                                        const Function<dim>& boundary_function,
+                                        const unsigned char boundary_component,
+                                        ConstraintMatrix& constraints,
+                                        const hp::MappingCollection<dim>& mapping_collection)
 {
-   std::vector<double> dof_values;
-   std::vector<unsigned int> face_dof_indices;
-   typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active ();
-   unsigned int dofs_per_face;
-   
-   switch (dim) {
-      case 2: {
-         for (; cell != dof_handler.end (); ++cell)
-            if (cell->at_boundary ())
-               for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
-                  if (cell->face (face)->boundary_indicator () == boundary_component) {
-                     // this is only implemented, if the FE is a Nédélec element
-                     typedef FiniteElement<dim> FEL;
-
-                     AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0),
-                       typename FEL::ExcInterpolationNotImplemented ());
-                     
-                     dofs_per_face = cell->get_fe ().dofs_per_face;
-                     dof_values.resize (dofs_per_face);
-                     
-                     for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
-                        dof_values[dof] = 0.0;
-                     
-                     QGauss<dim - 1> reference_face_quadrature (2 * (cell->get_fe ().degree));
-                     Quadrature<dim> face_quadrature
-                       = QProjector<dim>::project_to_face (reference_face_quadrature, face);
-                     FEValues<dim> fe_face_values (mapping_collection[cell->active_fe_index ()],
-                       cell->get_fe (), face_quadrature, update_jacobians | update_JxW_values
-                       | update_quadrature_points | update_values);
-                     
-                     internals::VectorTools::compute_face_projection (cell, face, fe_face_values,
-                       face_quadrature, boundary_function, first_vector_component, dof_values);
-                     face_dof_indices.resize (dofs_per_face);
-                     cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ());
-                     
-                     for (unsigned int dof = 0; dof < dofs_per_face; ++dof) {
-                        constraints.add_line (face_dof_indices[dof]);
-                        
-                        if (std::abs (dof_values[dof]) > 1e-14)
-                           constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]);
-                     }
-                  }
-         
-         break;
+  std::vector<double> dof_values;
+  std::vector<unsigned int> face_dof_indices;
+  typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active ();
+  unsigned int dofs_per_face;
+
+  switch (dim)
+    {
+      case 2:
+      {
+       for (; cell != dof_handler.end (); ++cell)
+         if (cell->at_boundary ())
+           for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+             if (cell->face (face)->boundary_indicator () == boundary_component)
+               {
+                                                  // this is only
+                                                  // implemented, if the
+                                                  // FE is a Nedelec
+                                                  // element
+                 typedef FiniteElement<dim> FEL;
+
+                 AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0),
+                              typename FEL::ExcInterpolationNotImplemented ());
+
+                 dofs_per_face = cell->get_fe ().dofs_per_face;
+                 dof_values.resize (dofs_per_face);
+
+                 for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+                   dof_values[dof] = 0.0;
+
+                 QGauss<dim - 1> reference_face_quadrature (2 * (cell->get_fe ().degree));
+                 Quadrature<dim> face_quadrature
+                   = QProjector<dim>::project_to_face (reference_face_quadrature, face);
+                 FEValues<dim> fe_face_values (mapping_collection[cell->active_fe_index ()],
+                                               cell->get_fe (), face_quadrature,
+                                               update_jacobians |
+                                               update_JxW_values |
+                                               update_quadrature_points |
+                                               update_values);
+
+                 internals::VectorTools
+                   ::compute_face_projection (cell, face, fe_face_values,
+                                              face_quadrature,
+                                              boundary_function,
+                                              first_vector_component,
+                                              dof_values);
+                 face_dof_indices.resize (dofs_per_face);
+                 cell->face (face)->get_dof_indices (face_dof_indices,
+                                                     cell->active_fe_index ());
+
+                 for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+                   {
+                     constraints.add_line (face_dof_indices[dof]);
+
+                     if (std::abs (dof_values[dof]) > 1e-14)
+                       constraints.set_inhomogeneity (face_dof_indices[dof],
+                                                      dof_values[dof]);
+                   }
+               }
+
+       break;
       }
-      
-      case 3: {
-        const unsigned int n_dofs = dof_handler.n_dofs ();
-         std::vector<double> computed_constraints (n_dofs);
-        std::vector<int> projected_dofs (n_dofs);
-        unsigned int degree;
-         unsigned int superdegree;
-        
-        for (unsigned int dof = 0; dof < n_dofs; ++dof)
-           projected_dofs[dof] = -1;
-        
-         for (; cell != dof_handler.end (); ++cell)
-            if (cell->at_boundary ())
-               for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
-                  if (cell->face (face)->boundary_indicator () == boundary_component) {
-                     // this is only implemented, if the FE is a Nédélec element
-                     typedef FiniteElement<dim> FEL;
-
-                     AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0),
-                       typename FEL::ExcInterpolationNotImplemented ());
-                     
-                     superdegree = cell->get_fe ().degree;
-                     degree = superdegree - 1;
-                     
-                     QGauss<dim - 2> reference_edge_quadrature (2 * superdegree);
-                     
-                     dofs_per_face = cell->get_fe ().dofs_per_face;
-                     dof_values.resize (dofs_per_face);
-                     
-                     for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
-                        dof_values[dof] = 0.0;
-                     
-                     face_dof_indices.resize (dofs_per_face);
-                     cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ());
-                     
-                     for (unsigned int line = 0; line < GeometryInfo<dim>::lines_per_face; ++line) {
-                        if (projected_dofs[face_dof_indices[line * superdegree]] < (int) degree) {
-                           Quadrature<dim> edge_quadrature = QProjector<dim>::project_to_face
-                             (QProjector<dim - 1>::project_to_face (reference_edge_quadrature, line),
-                             face);
-                           FEValues<dim> fe_edge_values (mapping_collection[cell->active_fe_index ()],
-                             cell->get_fe (), edge_quadrature, update_JxW_values | update_jacobians
-                             | update_quadrature_points | update_values);
-                           
-                           internals::VectorTools::compute_edge_projection (cell, face, line,
-                             fe_edge_values, edge_quadrature, boundary_function, first_vector_component,
-                             dof_values);
-                           
-                           for (unsigned int dof = line * superdegree; dof < (line + 1) * superdegree;
-                             ++dof)
-                              projected_dofs[face_dof_indices[dof]] = degree;
-                        }
-                        
-                        else
-                           for (unsigned int dof = line * superdegree; dof < (line + 1) * superdegree;
-                             ++dof)
-                              dof_values[dof] = computed_constraints[face_dof_indices[dof]];
-                     }
-                     
-                     if (degree > 0) {
-                        QGauss<dim - 1> reference_face_quadrature (2 * superdegree);
-                        Quadrature<dim> face_quadrature
-                          = QProjector<dim>::project_to_face (reference_face_quadrature, face);
-                        FEValues<dim> fe_face_values (mapping_collection[cell->active_fe_index ()],
-                          cell->get_fe (), face_quadrature, update_JxW_values | update_jacobians
-                          | update_quadrature_points | update_values);
-                        
-                        internals::VectorTools::compute_face_projection (cell, face, fe_face_values,
-                          boundary_function, first_vector_component, dof_values);
-                        
-                        for (unsigned int dof = GeometryInfo<dim>::lines_per_face * superdegree;
-                          dof < dofs_per_face; ++dof)
-                           projected_dofs[face_dof_indices[dof]] = degree;
-                     }
-                     
-                     for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
-                        if (std::abs (dof_values[dof]) > 1e-14)
-                           computed_constraints[face_dof_indices[dof]] = dof_values[dof];
-                  }
-         
-         for (unsigned int dof = 0; dof < n_dofs; ++dof)
-            if (projected_dofs[dof] != -1) {
-               constraints.add_line (dof);
-               constraints.set_inhomogeneity (dof, computed_constraints[dof]);
-            }
+
+      case 3:
+      {
+       const unsigned int n_dofs = dof_handler.n_dofs ();
+       std::vector<double> computed_constraints (n_dofs);
+       std::vector<int> projected_dofs (n_dofs);
+       unsigned int degree;
+       unsigned int superdegree;
+
+       for (unsigned int dof = 0; dof < n_dofs; ++dof)
+         projected_dofs[dof] = -1;
+
+       for (; cell != dof_handler.end (); ++cell)
+         if (cell->at_boundary ())
+           for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+             if (cell->face (face)->boundary_indicator () == boundary_component)
+               {
+                                                  // this is only
+                                                  // implemented, if the
+                                                  // FE is a Nedelec
+                                                  // element
+                 typedef FiniteElement<dim> FEL;
+
+                 AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0),
+                              typename FEL::ExcInterpolationNotImplemented ());
+
+                 superdegree = cell->get_fe ().degree;
+                 degree = superdegree - 1;
+
+                 QGauss<dim - 2> reference_edge_quadrature (2 * superdegree);
+
+                 dofs_per_face = cell->get_fe ().dofs_per_face;
+                 dof_values.resize (dofs_per_face);
+
+                 for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+                   dof_values[dof] = 0.0;
+
+                 face_dof_indices.resize (dofs_per_face);
+                 cell->face (face)->get_dof_indices (face_dof_indices,
+                                                     cell->active_fe_index ());
+
+                 for (unsigned int line = 0;
+                      line < GeometryInfo<dim>::lines_per_face; ++line)
+                   {
+                     if (projected_dofs[face_dof_indices[line * superdegree]]
+                         <
+                         (int) degree)
+                       {
+                         Quadrature<dim> edge_quadrature = QProjector<dim>::project_to_face
+                                                           (QProjector<dim - 1>::project_to_face (reference_edge_quadrature, line),
+                                                            face);
+                         FEValues<dim> fe_edge_values (mapping_collection[cell->active_fe_index ()],
+                                                       cell->get_fe (),
+                                                       edge_quadrature,
+                                                       update_JxW_values |
+                                                       update_jacobians |
+                                                       update_quadrature_points |
+                                                       update_values);
+
+                         internals::VectorTools
+                           ::compute_edge_projection (cell, face, line,
+                                                      fe_edge_values,
+                                                      edge_quadrature,
+                                                      boundary_function,
+                                                      first_vector_component,
+                                                      dof_values);
+
+                         for (unsigned int dof = line * superdegree;
+                              dof < (line + 1) * superdegree;
+                              ++dof)
+                           projected_dofs[face_dof_indices[dof]] = degree;
+                       }
+
+                     else
+                       for (unsigned int dof = line * superdegree;
+                            dof < (line + 1) * superdegree;
+                            ++dof)
+                         dof_values[dof] = computed_constraints[face_dof_indices[dof]];
+                   }
+
+                 if (degree > 0)
+                   {
+                     QGauss<dim - 1> reference_face_quadrature (2 * superdegree);
+                     Quadrature<dim> face_quadrature
+                       = QProjector<dim>::project_to_face (reference_face_quadrature,
+                                                           face);
+                     FEValues<dim> fe_face_values (mapping_collection[cell->active_fe_index ()],
+                                                   cell->get_fe (),
+                                                   face_quadrature,
+                                                   update_JxW_values |
+                                                   update_jacobians |
+                                                   update_quadrature_points |
+                                                   update_values);
+
+                     internals::VectorTools
+                       ::compute_face_projection (cell, face, fe_face_values,
+                                                  boundary_function,
+                                                  first_vector_component,
+                                                  dof_values);
+
+                     for (unsigned int dof = GeometryInfo<dim>::lines_per_face * superdegree;
+                          dof < dofs_per_face; ++dof)
+                       projected_dofs[face_dof_indices[dof]] = degree;
+                   }
+
+                 for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+                   if (std::abs (dof_values[dof]) > 1e-14)
+                     computed_constraints[face_dof_indices[dof]] = dof_values[dof];
+               }
+
+       for (unsigned int dof = 0; dof < n_dofs; ++dof)
+         if (projected_dofs[dof] != -1)
+           {
+             constraints.add_line (dof);
+             constraints.set_inhomogeneity (dof,
+                                            computed_constraints[dof]);
+           }
       }
-   }
+    }
 }
 
 
+
 template <int dim, template <int, int> class DH, int spacedim>
 void
-VectorTools::compute_no_normal_flux_constraints (const DH<dim,spacedim>         &dof_handler,
-                                                const unsigned int     first_vector_component,
-                                                const std::set<unsigned char> &boundary_ids,
-                                                ConstraintMatrix      &constraints,
-                                                const Mapping<dim, spacedim>    &mapping)
+VectorTools::
+compute_no_normal_flux_constraints (const DH<dim,spacedim>         &dof_handler,
+                                   const unsigned int     first_vector_component,
+                                   const std::set<unsigned char> &boundary_ids,
+                                   ConstraintMatrix      &constraints,
+                                   const Mapping<dim, spacedim>    &mapping)
 {
   Assert (dim > 1,
          ExcMessage ("This function is not useful in 1d because it amounts "

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.